V907 Sco Switched to the Eclipsing Mode Again

V907 Scorpii is a unique triple system in which the inner binary component has been reported to have switched on and off eclipses several times in modern history. In spite of its peculiarity, observational data on this system are surprisingly scarce. Here we make use of the recent Transiting Exoplanet Survey Satellite observations, as well as our own photometric and spectroscopic data, to expand the overall data set and study the V907 Sco system in more detail. Our analysis provides both new and improved values for several of its fundamental parameters: (i) the masses of the stars in the eclipsing binary are 2.74 +/- 0.02 M_0 and 2.56 +/- 0.02 M_0; and (ii) the third component is a solar-type star with mass 1.06 +/-0.11 M_0 (90% C.L.), orbiting the binary on an elongated orbit with an eccentricity of 0.47 +/- 0.02 and a period of 142.01 +/- 0.05 days. The intermittent intervals of time when eclipses of the inner binary are switched on and off are caused by a mutual 26.2 (+/- 2.6) inclination of the inner- and outer-orbit planes, and a favorable inclination of about 71 deg of the total angular momentum of the system. The nodal precession period is Pv = 63.5 +/- 3.3 yr. The inner binary will remain eclipsing for another approx 26 yr, offering an opportunity to significantly improve the parameters of the model. This is especially true during the next decade when the inner-orbit inclination will increase to nearly 90 degrees. Further spectroscopic observations are also desirable, as they can help to improve constraints on the system's orbital architecture and its physical parameters.


INTRODUCTION
Space-borne observations from the Kepler and Transiting Exoplanet Survey Satellite (TESS) missions have assisted with the discover and study a particularly interesting class of compact hierarchical triples (CHTs). These stellar systems are characterized by a smallenough ratio of the outer and inner orbital periods. A plethora of dynamical effects at various timescales can be detected with accurate photometric measurements for CHTs, and if properly interpreted with theoretical modeling, they may bring interesting constraints on many physical parameters of the system. An excellent review of this topic was recently published by Borkovits (2022).
Few archetype CHT systems have been known even before the above-mentioned revolution. In these cases, though, only the secular perturbations related to nodal and pericenter precession were typically detected. The short periods of these effects, sometimes even only a few decades, made them belong to the CHT class (see, e.g., Juryšek et al. 2018, for a review).
V907 Sco is a particularly interesting member of this class of objects. While occasionally observed in the second half of the 20th century, and first cryptically reported to be an eclipsing binary by Strohmeier et al. (1964), a truly comprehensive study of V907 Sco was presented by Lacy et al. (1999). We encourage the interested reader to learn about the sometimes tangled his-tory of the V907 Sco analysis from the accumulated photometric data in this reference (principally Sec. 1). Here we only mention that Lacy et al. (1999) were the first to fully understand the background story of V907 Sco by proposing that it is in fact a triple in which the previously known binary is accompanied by a third star in a sufficiently close orbit. This was necessary to explain the surprisingly short changing periods during which the binary ceased to be eclipsing. Unfortunately Lacy et al. (1999) could not continue to collect accurate photometry of the eclipsing component, since it was not eclipsing during the period of their work. Their important contribution to the observational evidence about V907 Sco thus consisted in taking the first series of the spectroscopic data. While useful in several aspects (e.g., confirming precisely the period of the eclipsing binary, determining the semiamplitudes of the radial velocity curves, and noting the narrowness of the spectral lines), Lacy et al. (1999) were unlucky by having available only a sparse set of spectra spread over an interval of 11 years. Their analysis was therefore ill-suited for such a fastevolving system and led to a few mistakes. First, they attempted to determine masses of the stars in the inner binary but chose an incorrect orbital inclination at the reference epoch of their observations. As a result, their reported masses were slightly incorrect (Sec. 2). More importantly though, the sparsity of their data led them to determine basic parameters of the outer orbit, the period and orbital eccentricity, incorrectly (see again Sec. 2). To their credit, however, we finally note that Lacy et al. (1999) also considered the possibility that V907 Sco belongs to the open cluster M7. Comparing the estimated distance to both M7 and V907 Sco, though the latter arguably rather approximate at that moment, they rightly rejected the association. Their conclusion was recently confirmed with astrometric measurements from Gaia DR3 showing that V907 Sco has a parallax of 2.279 ± 0.033 mas (Gaia Collaboration et al. 2022). This translates to a distance of 439 ± 6 pc, well beyond that of M7.
The rather rough and partially incorrect analysis of Lacy et al. (1999) let these authors to predict that the eclipses of V907 Sco will start again in 2030 ± 5 in their concluding Section 4. However, Borkovits (2022) noted that the recent observations of TESS have already revealed an onset of V907 Sco eclipses in the Sector 13 data (in mid 2019), and confirmed their expected trend towards deeper eclipses in the Sector 39 data (in mid 2021). Clearly, the situation demands a new analysis of V907 Sco, since the recent emergence of its eclipsing state alone proves that some assumptions and/or data treatment in Lacy et al. (1999) have not been cor-rect. This is the primary purpose of this work. However, realizing that past data (especially the pre-TESS photometry) could be improved upon, we also aim to extend the empirical data set on this unusual triple system with new observations. The paper is organized as follows. In Sec. 2 we describe photometric and spectroscopic observations available to us, including archival, previously published, and new data we obtained in the 2022 season. Realizing that several orbital parameters are evolving quickly, we split the analysis into two steps. First, still in Sec. 2, we focus on the bulk of the 2021-2022 data. They conveniently cover a period slightly longer than the outer orbit period, allowing us to properly characterize it, but short enough that the changes of the orbital parameters are minimal. We can thus use plain Keplerian orbits for the data analysis here. This allows us to determine orbital and physical parameters of the inner binary in the V907 Sco system (included the masses and radii of the stars). With this information available, we can thus use the archival photometry and derive inclination changes of the eclipsing binary over more than a century. In Sec. 3 we develop a model accounting for mutual gravitational interactions of all three stars in the system to describe the observed variations of the orbital elements over long-term period (many cycles of the outer orbit). Confrontation with the data allows us to determine some parameters of the third star (including its mass and tilt of its orbital plane with respect to the orbital plane of the inner binary). Finally, in Sec. 4, we outline future prospects to study V907 Sco in even more detail.

OBSERVATIONS AND SHORT-TERM ANALYSIS
Both photometric and spectroscopic observations are available for V907 Sco. The photometry falls into the following three categories: (i) archival data sparsely spanning more than a century, some of which have also been mentioned by Lacy et al. (1999), (ii) the abovementioned TESS observations, and (iii) our new data from the 2022 season. As for the spectroscopy, we conducted an intensive observing campaign in 2022, knowing how quickly the V907 Sco system evolves. Such data have the advantage of minimizing long-term evolutionary effects in spectra. There are also radial velocity measurements published by Lacy et al. (1999) (their Table 3). However, they are spread over a time interval of nearly 11 years, and thus their analysis strongly depends on the ability to consistently include the evolutionary model of the whole system. For that reason, we do not consider them here, but save their analysis for the future work. , each spanning about 27 days and folded on the 3.7762049 day period of the eclipsing binary. The ordinate is differential magnitude, which has been shifted by 0.01 mag for the Sector 13 data for better visualization. The red lines are the respective best-fits obtained using the PHOEBE code (Table 1). In the Sector 13 data, only a very subtle dip is noticed at the zero phase of the primary eclipse, implying the stellar configuration seen by the observer was nearly grazing.  Table 2). The abscissa is the phase of P1 cycle, and the ordinate is the differential magnitude, with a shift for sake of visibility. In what follows, we describe each photometric data set, in order of decreasing precision.

Photometric observations
TESS photometry.-We searched all TESS sectors for field-of-views containing V907 Sco and found that the object was observed twice, in Sectors 13 and 39 (V907 Sco = TIC 261030551; see also Borkovits 2022). Each of them represents about 27 day long period of continuous data, with an approximately day-long data downlink gap in the middle. Target photometry was extracted using the publicly available lightkurve package (e.g., Lightkurve Collaboration et al. 2018). Coincidentally, and of high value to this study, Sector 13 observations began at the onset of the eclipses, as the light curve shows only a very shallow signal/dip at the primary eclipse phase. The exposure cadence was about 30 minutes. The Sector 39 data are superior in their quality and value for our analysis. Not only do they have a more rapid cadence (10 minutes), but they were taken at the moment when the stars in the inner binary of V907 Sco exhibited sufficiently deep eclipses (about 0.06 mag in the case of the primary eclipse) for a thorough analysis. With the expected millimagnitude-level accuracy of the TESS photometry (e.g., Ricker et al. 2015), the Sector 39 data are by far the most precise photometry we possess on the V907 Sco system. Consequently, they provide important constraints on fundamental orbital and physical parameters of the inner binary. Figure 1 shows the TESS light curves from both sectors, phase-folded on the best-fitting mean period P 1 of the binary. We note that our definition of the zero phase is shifted by half an orbital period compared to the definition of Lacy et al. (1999). This is because we prefer to associate the zero phase with the deeper eclipse, when the smaller (less massive) star in the binary eclipses the larger (heavier) star.
We used the well-tested software package PHOEBE (Prša & Zwitter 2005), based on the Wilson & Devinney (1971) code, to obtain fundamental parameters of the inner binary system using (i) TESS Sector 39 photometric data, and (ii) our new spectroscopic data described in Sec. 2.2. Ideally, both datasets should have been taken at the same epoch, but they were not. Luckily, the one year difference between them is short. Figure 6 shows that the binary inclination i 1 changed by only ≃ 1 • from 2021 to 2022, and this makes the semiamplitudes of the radial velocities change by about 0.3%. This is a negligible change in our solution. Our modeling required only minimum assumptions, such as rotational synchronicity of the stars in the binary (F i = 1), while other parameters were left adjusted. This concerns not only the primary ones, such as stellar masses and radii, but also more subtle ones, such as the gravity brightening or albedo. The latter primarily influences the ellipsoidal modulation of the photometric signal outside the eclipses (Fig. 1). The results are summarized in Table 1. The Table 1. Orbital and physical parameters of the inner binary in V907 Sco as derived from the combined lightcurve and radial velocity fitting. As for the orbital inclination values, see  Lacy et al. (1999). This is because these authors assumed an inaccurate value of the inclination i 1 of the binary orbital plane at the effective epoch of their spectroscopic observations. The corresponding radii are 2.53 R ⊙ and 2.68 R ⊙ , consistent with the values expected for the B9V and B8V main sequence stars. The contribution of the third light in the TESS passband is about 6%, and the binary has a small, but well detectable, orbital eccentricity of about 0.008. Both values are found to be consistent with our modeling discussed below (see Sec. 4).
FRAM photometry.-The most recent data, taken between April and September 2022, had the purpose of confirming the trend of the eclipse depth increase (and thus the inclination of the binary orbital plane). They were obtained with the F/(Ph)otometric Robotic Atmospheric Monitor (FRAM) 30 cm telescope at Pierre Auger Observatory, Argentina (e.g., Aab et al. 2021). This is an Orion ODK 300/2040mm system equipped with the CCD camera MII G4-16000; standard reduction routines using dark frames and flat fields were used. The small aperture of the instrument results in the FRAM photometry having less precision than the TESS photometry; a typical uncertainty in an individual frame is ≃ 0.01 magnitude. However, given the fact that the current eclipse depth is already about ten times larger, the FRAM data are still very useful (see the top curve on Fig. 2 and Fig. 10). We obtained observations at four epochs in 2022, clearly demonstrating the expected progression towards deeper eclipses. Importantly, when this trend is mapped to the increasing inclination of the binary orbital plane (Table 2), the time evolution is not linear. Indeed, in Sec. 3.4 we find reasons for this nonlinearity in the dynamical model.
Archival photometry.-In addition to our recent photometric data, we also utilized archival resources. Obviously, this is motivated by the reported intermissions in V907 Sco eclipses (Lacy et al. 1999), and thus we seek the ability to describe them. We used data from: (i) the Archives of Photographic PLates for Astronomical USE (APPLAUSE) project of digitalizing the old photographic plates available online 1 (Groote et al. 2014), and (ii) the Digital Access to a Sky Century at Harvard (DASCH) database of photographic plates from Harvard observatory 2 (Tang et al. 2013). Given the fact that DASCH plates cover an interval of several decades, we could even parse them into several independent intervals in between 1900 and 1920. In spite of their rather poor quality (see the bottom curves on Fig. 2), these data are still very useful as they approximately tell us trends in the inclination variations over a long timespan. Obviously, we downscale the weight of these data by assigning large uncertainty to the derived inclination values (Table 2). Despite our efforts, we unfortunately did not find more photometric datasets from existing survey programs that would provide useful information about V907 Sco. An example is ASAS-sn in 2022, in which observations in the V907 Sco field are subject to a spurious instrumental scatter of about 1.5 mag.  Compilation of the inclination values at different epochs.-Having calibrated the physical parameters of stars in the inner binary from the TESS photometry, we can now fix them and use other photometric data to determine only the inclination of the eclipsing binary orbital plane at their epoch. For this task we used the PHOEBE code as above. Note also that ephemerides for the eclipses are determined accurately enough for our purpose: a fractional uncertainty of ∆P 1 /P 1 ≃ 3 × 10 −6 expands to a ≃ 0.03 phase uncertainty after a centurylong interval (i.e., roughly representing 10 4 cycles of the binary). This is due to the extreme accuracy of the TESS photometry. Nevertheless, we adjusted the epoch of the primary minimum individually in each of the datasets, such that the fit to the model is optimized. As expected, this adjustment turned out to be minimal.
The results are summarized in Table 2 for both the less accurate archival data and the much more precise modern data. From the TESS Sector 13 the value was only roughly estimated, while for the TESS Sector 39 the data were split evenly into two parts (roughly two weeks each), and the inclination was determined for each independently. The two i 1 values are indeed statistically different and match well the expected trend (Sec. 3.4). We are aware of the fact that the TESS data suffer from light contamination from nearby sources and its level may change between different sectors, or even during one sector of data. Due to strong correlation between the inclination and additional stray light this may play some role. However, the difference would be only very small, and the principal results of the analysis when using one Table 3. Parameters of the outer orbit in V907 Sco from available radial velocity data in 2022. The semi-amplitude K corresponds to the motion of the inner binary center of mass with respect to the center of mass of the whole system.

Parameter
Value  single averaged point or two instead would remain almost unchanged. All inclination values are later used in Sec. 3 to constrain parameters of the third component in the system from its gravitational perturbations of the inner binary.

Spectroscopic observations
Our new spectroscopic observations of V907 Sco were obtained using the CHIRON echelle high-resolution spectrograph (R ≃ 28000) mounted on a 1.5 m SMARTS telescope at CTIO observatory located on Cerro Tololo, Chile (e.g., Tokovinin et al. 2013). Altogether, we have 66 spectra taken in between March and October 2022. Typical exposure times were 480 s, and the usable spectral range covered a wavelength range of 450 to 890 nm.  Table 3; the grey horizontal line is at K e2 cos ω2 ≃ −6 km s −1 , a center to which semi-amplitude K is referred to. The abscissa is the Julian date with 2450000 subtracted. Bottom part: residuals of the radial velocities after motion of the inner and outer orbit has been subtracted from the data.
The whole reduction process including rectification of the spectra, identifying the components, and measuring their radial velocities, all of which were accomplished using the reSPEFO 2 package 3 (Harmanec et al. 2020).
The resulting values at each epoch observed are listed in supplementary materials (Table A1). For most of the spectra, we were able to identify a pair of lines that belong to the two stars in the inner binary of V907 Sco (corresponding mainly to Fe and Si; see Fig. 3 for two examples at different phases of the binary motion). Unfortunately, we could not identify spectral lines of the third star in the system (see also Lacy et al. 1999). This non-detection can be explained by a low luminosity (and thus mass) compared to the stars in the inner binary. We estimate a ≃ 5 magnitude difference between the third star and the binary, or a few percent contribution to the total luminosity of the system at maximum. Given the masses derived in Table 1 for the components in the binary, we estimate an upper mass limit on the third star of ≃ 1.05 M ⊙ . This nicely matches the results in Sec. 3.4, where we infer the third component mass from its dynamical effects on the binary.
As mentioned in the previous Section, we used the code PHOEBE to analyse simultaneously the TESS-Sector 39 photometry and our derived radial velocities in order to solve for the orbital and physical parameters of the eclipsing binary (i.e., all signal at the P 1 period). This resulted in the values presented in Table 1 and the velocity curve fits shown in Fig. 4. The data are matched by the model quite well. Next we analysed the residual signal in the measured radial velocities and expected to find the ≃ 99 day period of the outer orbit reported by Lacy et al. (1999). To our surprise, we obtained something else (see Fig. 5). When fitting a Keplerian orbit to the velocity residuals, we find a long period of P 2 ≃ 142 days with a significant eccentricity of e 2 ≃ 0.47. Table 3 presents the resulting parameters. The model fit is rather good, with the final residuals having a dispersion of only 0.64 km s −1 (at the level of the radial velocity measurement uncertainty) and no further signals being detectable (see Supplementary materials, Table A1). We unfortunately lack data near the minimum and maximum of the radial velocity curve for the outer orbit (this happened partly due to bad weather and partly due to overload of the telescope by other projects). This leaves us with parameter uncertainties in Table 3 slightly larger than we would wished. As an example, the minimum mass of the third component, as derived from the mass function, ranges between 0.95 M ⊙ and 1.09 M ⊙ . We only note that this value is satisfactorily compatible both (i) with the maximum third component mass from the absence of its spectroscopic evidence and the third light contribution, and (ii) also with the solution based on photometric data in Sec. 3.4.
Finally, we note that the large value of the eccentricity of the outer orbit is still comfortably within stability limits. Assuming stellar masses in the eclipsing binary derived above, and m 2 ≃ 1 M ⊙ (Sec. 3.4), we find that the system is stable up to e stab 2 ≃ 0.63 according to the criteria in Eggleton & Kiseleva (1995) and Mardling & Aarseth (2001).

LONG-TERM EVOLUTION OF THE SYSTEM
In this section, we analyse the observational results of the V907 Sco system collected above and summarized in Table 2 using a simple, but adequately accurate model. This allows us to describe changes in the system's architecture and its parameters, such as orbital inclinations of both inner and outer orbits, over a timescale of decades. Before discussing the results, we review the basic aspects of our theoretical approach and relevant variables. Additional details can be found in Vokrouhlický & Zasche (2022) or Juryšek et al. (2018), while an interested reader may also want to revisit a classical series of papers by Soderhjelm (1975Soderhjelm ( , 1982Soderhjelm ( , 1984. We find it useful to start with a simple analytic secular model in Sec. 3.1. Effects with long period P 2 are not described by this part, but its strength stems from analyticity, which effectively helps us pre-constrain the model parameters. Next, in Sec. 3.2, we also include the long-period effects using a computationally more demanding numerical model.

Simple analytical model for the secular part
We first introduce notation and important parameters used throughout this section: (i) m 1a and m 1b are masses of the two stars in the eclipsing binary (inner orbit), defining their total and reduced masses simply as M 1 = m 1a + m 1b and µ 1 = m 1a m 1b /M 1 , (ii) m 2 is the mass of the third star in the system (constituting the outer orbit with respect to the barycenter of the inner orbit), and the total mass is M 2 = M 1 + m 2 . We denote also the mean orbital periods P 1 and P 2 of the inner and outer orbits, or alternatively also the corresponding mean motions n 1 = 2π/P 1 and n 2 = 2π/P 2 . Finally, we assume the inner orbit is circular (e 1 = 0), and the outer orbit may have an arbitrary eccentricity e 2 , defining also a convenient parameter η 2 = 1 − e 2 2 . The principal variables of interest are the orbital angular momenta L 1 and L 2 of the inner and outer orbits, which secularly change due to gravitational interactions in the system. When the latter is restricted to the quadrupole coupling, justifiable in the V907 Sco system for which the mass ratio q of the stars in the inner system has been found close to unity, (i) the angular momentum magnitudes L 1 and L 2 are secularly conserved, and (ii) the unit vectors l 1 and l 2 , defining orientation of the respective angular momenta in space, steadily precess about the conserved total angular momentum of the triple L = L l = L 1 + L 2 (here we neglect small contributions of the rotational angular momenta of the three stars, again justifiable for the V907 Sco system, where the narrow spectral lines of the components in the eclipsing system imply their slow rotation). Additionally, as the angular momentum vectors of the inner and outer system roll on the respective cones about L, their mutual configuration is fixed such that their mutual angle J defined by cos J = l 1 · l 2 (also depicting the angle between the orbital planes of the inner and outer orbits), remains secularly constant.
It should be recalled that within the point-mass, threebody problem adopted here, the simple solution outlined above holds when J ≤ 40 • or J ≥ 140 • . In between these values a more complicated, Kozai-Lidov regime occurs (e.g., Soderhjelm 1982;Farago & Laskar 2010). Luckily, this complication is of no concern for V907 Sco. First, we show below that the observations require configuration for which J is restricted to a narrow range of values about ≃ 26.2 • (Sec. 3.4). Second, the long-term stable triples often rearrange such that the Kozai-Lidov regime is effectively eliminated by mutual (tidal) interaction of the stars in the inner orbit (e.g., Soderhjelm 1984); this even allows systems with J = 90 • to exist, the Algol case representing the archetype member of this class (e.g., Baron et al. 2012).
Mathematical representation of the steady precession state of the inner and outer angular momenta reads simply l 1 (t) =l 1 cos α + l ×l 1 sin α + l l ·l 1 (1 − cos α) , (1) and similarly for l 2 (t) (recall l is the unit vector of the conserved total angular momentum of the system).
Here,l 1 is the unit vector of the inner orbit angular momentum at some arbitrary reference time T 0 (i.e., l 1 = l 1 (T 0 )), and α = ν(t − T 0 ), where ν is the fundamental precession frequency. Some authors are satisfied with ν as a free, solved-for parameter, but the analytical theory offers an important step further, namely to interpret ν in terms of more fundamental parameters such as J and masses of stars in the triple system. Remaining still on the grounds of the quadrupole coupling approximation, we have (e.g., Soderhjelm 1975;Breiter & Vokrouhlický 2015) ν is the ratio of the inner and outer orbit angular momenta (i.e., γ = L 1 /L 2 ). Typical triple systems have γ ≪ 1, as the outer orbit carries much larger share of the system's angular momentum. However, V907 Sco belongs to a small group of outlier cases in this respect. This is because (i) it is rather compact with P 2 /P 1 ≃ 37.6, and (ii) the third component is most likely a solar-type star (Sec. 3.4). Adopting the canonical value M 1 ≃ 5.3 M ⊙ , the plausible m 2 ≃ 1 M ⊙ mass implies γ ≃ 0.47, thus a partial share of the inner and outer orbits on the total angular momentum. Large variations of the inner orbit inclination i 1 (witnessed by the intermittent periods of eclipses and no-eclipses) must therefore be accompanied with certain variations of the inclination i 2 of the outer orbit. The analytical solution (1) provides this result readily: orienting the reference system as usual, such that the XY plane coincides with the sky-plane, the inclination i 1 (t) of the eclipsing system is simply cos i 1 (t) = l 1 (t) · e 3 , where e T 3 = (0, 0, 1) (and similarly for the outer orbit, with just replacing l 1 (t) by l 2 (t)).
It is useful at this moment to overview free parameters of the analytical model whose values are to be adjusted (solved for) when confronting with the data (Sec. 3.4). Choosing the reference epoch T 0 , we first havel 1 andl 2 (normal unit vectors to the planes of the inner and outer orbits) to select (note their composition with magnitudes L 1 and L 2 of the angular momenta define the unit vector l of the total angular momentum of the system). Being unit vectors, each of the two would apparently depend on two solved-for parameters. But given the nature of the data we have, photometric and spectroscopic sets, there is an arbitrariness in selection of the reference direction in the sky-plane. As a result, there are only three solved-for parameters needed to characterize the unit vectorsl 1 andl 2 . For instance, we may assumel whereî 1 is the inner orbit inclination at T 0 , and whereî 2 is the outer orbit inclination at T 0 , andΩ 2 is the nodal phase tilt about e 3 of the outer orbit with respect to the inner orbit. Therefore, we have three solved-for model parameters (î 1 ,î 2 ,Ω 2 ) of the geometrical nature. If only photometric data were available, one could complement this set of three parameters with the precession frequency ν to be the fourth solved-for parameter, and use the inner system inclination changes to determine their values (a formulation dating back to Soderhjelm 1975). A partial step forward was introduced by Juryšek et al. (2018), who used the photometrically-derived constraint on ν, and formula Eq.
(2) to set constraints on further physical parameters, such as m 2 and often unknown P 2 for their systems.
Here, though, we are in a slightly more comfortable situation with both photometric and spectroscopic data available, and we can thus afford more general formulation (though a really ambitious approach needs to wait when more accurate data are available in the future). At this moment, our procedure is similar to that in Vokrouhlický & Zasche (2022). We consider the frequency ν as an auxiliary variable only, and use Eq. (2) to connect it with parameters of more interest, such as the stellar masses or mutual inclination J of the inner and outer orbital planes. This does not mean that the new parameters would have a fully correlated solution. This is because (i) J, for instance, depends onl 1 and l 2 on a more fundamental level (recall cos J =l 1 ·l 2 ), and (ii) constraints following from the spectroscopic observations help to decorrelate the stellar masses. So in principle, all quantities on the right hand side of the formula (2) may be adjustable in the process of data fitting. However, for some of them we assume that their pre-constrained values from the analysis in Sec. 2 are precise to such a level that the data describing secular evolution of the system would not result in improvement (and in the same time, their uncertainty would not significantly affect solution of the parameters which are fitted at this stage). For the sake of simplicity, we thus adopt their best-fit values and consider them constant. This is the case of the mean orbital periods P 1 and P 2 of the inner and outer orbits, and the eccentricity e 2 of the outer orbit. In fact, the latter appears in ν and L 2 only via the η 2 factor, and therefore depends on e 2 2 rather than e 2 . We are then left with the stellar masses, of which m 2 is perhaps the most interesting. This is because the analysis of the spectroscopic data in Sec. 2.2 helps to directly constrain both m 1a and m 1b , since i 1 is well known at the epoch of the spectroscopic observations from the conjoined photometric observations of the eclipsing system. A similar determination of m 2 from spectroscopy would need i 2 to be known, but this is available from the long-term evolution model only (and not directly from observations). Given the still limited quality of available data, we thus resolve to provide a preliminary solution for the V907 Sco system at this moment. We fixed m 1a and m 1b values determined in Sec. 2.1, and we adopt m 2 as the fourth solved for parameter.

Numerical model for the long-period part
We used the secular model from Sec. 3.1 to fit the photometric observations (Sec. 3.4 and Fig. 6 in particular) and reached the following conclusions: (i) the model helps us to match the mean trend of the data adequately, and allows us to characterize a large portion of the parameter space which cannot lead to satisfactory results (thus pre-constrain the parameter space), but (ii) formally the best fit is statistically not acceptable (the χ 2 from Eq. (6) is simply too large). While looking into details, we noticed the reason for the situation: on top of the secular terms, there are also long-period perturbations (periodicity P 2 and its multiples) that contribute in a non-negligible way to the long-term evolution of i 1 and i 2 values. Since the post-2019 photometric data are rather accurate, it is important to include these terms into the analysis. At the lowest-orders in the outer eccentricity e 2 they can be expressed analytically, e.g., Soderhjelm (1975) (but see already Brown 1936), and represent a signal with periodicity P 2 /2. For larger e 2 values, which is the case of V907 Sco, an accurate analytic formulation is not available to us. For that reason, we decided to resolve the situation by using a fully numerical model, and profit from the a priori parameter constraints from the secular model described above. Interestingly, we find that at high e 2 regime the leading periodicity of the long-term perturbation of both i 1 and i 2 becomes P 2 (Fig. 6). 4 The set-up of the numerical model uses a point mass configuration and Jacobi coordinates (e.g., Soderhjelm 1982, Eqs. (1)). The near-coplanarity of the inner-and outer-orbits in the V907 Sco system also does not require us to include the gravitational interaction of the stars in the eclipsing system as extended bodies (beyond the point-mass level) that would be otherwise needed for long-term stability if J is large. None of the observable quantities depends on these three-body effects. Because the long period perturbations depend on the mean longitude in orbit of the outer system at the reference epoch T 0 , we must add it to the set of adjustable parameters. Taking the argument of pericenter ω 2 from the analysis of radial velocities (Sec. 2.2), we let the mean anomaly ℓ 2 be the fifth solved-for parameter. There is no a priori constraint on ℓ 2 , and its value is allowed to sample the full interval of 360 • .

Dataset for long-term analysis
Data suitable for long-term analysis of the evolution of V907 Sco consist of results from our analysis of the photometric and spectroscopic observations at individual and separated intervals of time presented in Sec. 2. Each of them provided at a certain epoch t j , usually at the center of the relevant time interval, a value of a specific variable v j together with its uncertainty ∆v j . The triples (t j , v j , ∆v j ), with j = 1, . . . , N data , are the basis of our work in this section. The model described in Secs. 3.1 and 3.2 allows us to compute the valuesv j of these variables at the same epochs t j , and we can thus use the standard χ 2 metric defined by as a measure of the model quality. As expected,v j = v j (p), where p = (î 1 ,î 2 ,Ω 2 , m 2 , ℓ 2 ) constitute the fivedimensional parametric space. Our procedure obviously seeks minimization of χ 2 (p), reaching a minimum value χ 2 min for a certain, best-fit set of parameters p min . Making sure that the fit is statistically acceptable (to that Figure 6. Inclination i1 (label 1) and i2 (label 2) of the inner-(binary) and the outer-component orbits of V907 Sco over the past century and extending towards the next decade (the abscissa in calendar years). Red symbols are data with uncertainties determined in Sec. 2 (Table 2). Right panel provides a zoom on the modern, most accurate data (depicted using the rectangle on the right panels) and shows only the i1 values. The best-fitting model from Figure 7   . Left: distribution of solutions for which χ 2 ≤ χ 2 ⋆ ≃ 13.13 projected onto the plane of mass m2 (abscissa) and precession period Pν (ordinate). Solution density is indicated by greyscale, white for no possible solution, black for the largest number of solutions. The dashed lines delimit the zone of acceptable solutions. The blue star at m2 = 1.062 M⊙ and Pν = 63.5 yr shows location of the best-fitting solution with χ 2 min = 3.89. Right: the minimum χ 2 (abscissa) for a given Pν value (ordinate). The dashed vertical line is the limit for accepted solutions χ 2 = χ 2 ⋆ , the solid vertical line is the best achieved value χ 2 min ≃ 3.89.
end we use the Q-function test, see Press et al. 2007, Chap. 15.1), we seek confidence limits on estimated model parameters p. This is achieved by selecting a certain volume V in the parameter space characterized by χ 2 ≤ χ 2 ⋆ = χ 2 min + ∆χ 2 . For a five-dimensional parameter space, and a 90% confidence level limit, we have to choose ∆χ 2 = 9.24 (see, e.g., Press et al. 2007, Chap. 15.6). Projection of V onto lower-dimensional sections of the parameter space allows us to re-map the statistical analysis onto other variables of interest, which do not constitute the basic set of dimensions by p (for instance, we may characterize this way acceptable values of the precession period P ν = 2π/ν, as shown in Figs. 7 and 9, or the mutual angle J of the orbital planes, as shown in Fig. 9). Even simpler is just projecting V onto the individual axes defined by each of the parameters p (for instance to infer limits in the companion mass m 2 , Figs. 7 and 8).
Having outlined our procedure, we now return to the choice of observables v j . Their sector based on photometric observations consists of data summarized in Table 2, namely the inferred inclination values i 1 of the inner binary orbit from mutual eclipses. This provides 14 datapoints. Next, the spectroscopic observations taken at CTIO in 2022 provide further constrains. First, the analysis of the radial velocity semi-amplitudes of the stars in the eclipsing binary (Fig. 4), and the known inclination i 1 at their epoch, provides values of masses m 1a and m 1b . We consider them fixed. The analysis of the spectroscopic observations further provides radial velocity data for the outer orbit with period P 2 , namely how it affects the systemic velocity of the barycenter of the eclipsing binary (Fig. 5). This provides a constraint on the argument of pericenter ω 2 and the eccentricity e 2 that we also consider fixed in our analysis (we only tested the influence of making e 2 smaller or larger than nominal). Being most interested in constraining the mass m 2 of the third, unseen star, we observe the semi-amplitude of the radial velocity curve in Fig. 5, from which we can derive the projected semimajor axis a 2,proj = a 2 (m 2 /M 2 ) sin i 2 (the semimajor axis of the outer orbit is given by the Kepler's third law n 2 2 a 3 2 = GM 2 ). We obtain, a 2,proj,obs (t s ) = 34.14 R ⊙ with an uncertainty ∆a 2,proj,obs = 2.04 R ⊙ at a reference HJD epoch t s = 59760. This constitutes our additional data to be fit by our model. Finally, we chose the reference epoch T 0 = 59383.9 heliocentric Julian date, namely the mid epoch of the TESS sector 39b observations (Tab. 2).

Results
We ran two sets of simulations: (i) first using only the photometric data, and (ii) second using both photometric data and spectroscopic constraints on the projected semimajor axis of the outer orbit. Our intention was to investigate how good the solution was using the photometric data alone, and how it potentially improves if constraints from spectroscopy are added. Figure 6 shows the best-fitting solution of the inclination i 1 data, including also coupled variations of the outer orbit inclination i 2 . The inclination of the conserved total angular momentum (i.e., the invariable plane) is ≃ 71 • , indicating the solution belongs to the class 2 discussed by Lacy et al. (1999) in Sec. 4 of their analysis (which they rightly favoured). The archival, pre-2000 data are formally inaccurate, nevertheless they help to significantly constrain the solution. The modern data, especially those after 2021, are very accurate and help constrain the solution even more. The important advantage of V907 Sco over the HS Hya case (Vokrouhlický & Zasche 2022) is the fact that the data cover a time interval much longer than the precession period P ν = 2π/ν. This makes the free-parameter values much better constrained: compare Fig. 7 here with Figs. 3 and 4 in Vokrouhlický & Zasche (2022). Figure 8 presents the probability density distribution of the third mass m 2 . We obtain m 2 = 1.06 +0.11 −0.10 M ⊙ , representing a 90% confidence level interval. We briefly tested the sensitivity of the solution to changing the value of the eccentricity e 2 of the outer orbit (which has been fixed in the illustrated solution to its nominal value 0.467). We found that slightly smaller or larger e 2 values would shift the m 2 solution to slightly larger or smaller values. This is to be expected. For instance, smaller e 2 would place the outer star more distant to the eclipsing binary. In order to produce the observed variations of i 1 values, one would thus need larger mass m 2 (and vice versa). Figure 9 presents the probability density distribution of the precession period P ν , reflected in periodicity of variations of i 1 and i 2 in Fig. 6, and the mutual angle J of the inner and outer orbits orbital planes. The 90% confidence values are P ν = 63.5 +3.3 −2.6 yr and J = 26.2 +2.6 −2.2 degrees. The value of J is obviously the sum of the semiamplitudes of i 1 and i 2 variations seen in Fig. 6, while the ratio of their sine values is Figure 9. Probability density distribution of solutions for precession period Pν and mutual inclination J of inner and outer orbits of the V907 Sco system. The available dataset of photometric observations was used. The non-zero values correspond to the projection of five-dimensional parameter-space zone V containing all admissible solutions within 90% confidence limit. The vertical solid line is the median value of the distribution. γ ≃ 0.47 from Eq. (3). The J value is just large enough to produce interesting orbital effects, but smaller than the critical value of ≃ 40 • for the onset of Kozai-Lidov effects (e.g., Soderhjelm 1982;Farago & Laskar 2010). The use of a simple point-mass model in our simulations is therefore well justified. We also note that the argument of pericenter of the outer orbit exhibits a drift of ≃ 1.35 degrees per year (this is not an observationally detected value yet, but it follows from our numerical simulation; see also Eq. (8) in Vokrouhlický & Zasche (2022)). This large value stems from compactness of the system (P 2 /P 1 not too large) and large mass of the binary. This effect needs to be taken into account, apart from the changing inclination values i 1 and i 2 , if we were to interpret the radial velocity measurements reported in Lacy et al. (1999) (Table 3; see also Sec. 4).
Next, we added to the dataset the projected semiamplitude a 2,proj constraint from the spectroscopic data in 2022. We found that this single data point did not improve the solution meaningfully at this moment. This is because (i) the photometric data already constrain the model very well, and (ii) the spectroscopic information is not very accurate yet, principally because the observations missed to detect precise moment of the pericenter passage of the outer orbit (Fig. 5).

DISCUSSION AND CONCLUSIONS
The recent onset of eclipses of the binary component of the V907 Sco system recorded in the TESS data opened a new opportunity to study this interesting triple. This work is only an initial attempt in this direction, but it nevertheless reveals important corrections to prior results. They are as follows: • masses of stars in the eclipsing binary, 2.74 ± 0.02 M ⊙ and 2.56 ± 0.02 M ⊙ , are about 8% larger than previously thought; this is because Lacy et al. (1999) assumed in their work an inaccurate value of the inclination i 1 of the binary orbital plane at the effective epoch of their spectroscopic observations; • our single-season spectroscopic data and a more precise modeling work allowed us to constrain parameters of the orbit of the third star about the center of mass of the binary; unlike previously thought, we find it has a significant eccentricity of 0.47 ± 0.02 and a substantially longer period of 142.01 ± 0.05 days; • secular and long-period perturbations in i 1 produced by the third component in the system allow us to constrain its mass of 1.06 +0.11 −0.10 M ⊙ and determine mutual inclination of the two orbital planes of J = 26.2 +2.6 −2.2 degrees. The sum of all known parameters, both orbital and physical, is given in Tables 1 and 3. In principle, one may use those data to estimate, at least at zero order, the age of stars in the eclipsing binary. However, we find that the result may sensitively depend on their yet unknown metalicity. In any case, the values range in between ≃ 300 and ≃ 500 Myr. This would mean that V907 Sco is a rather young stellar system. More detailed analysis of this issue, including consistency with the estimate of a circularization timescale for the inner binary orbit, is left for future efforts.
The archival and incidental photometry data of V907 Sco proved to be valuable for description of the long-term changes in its architecture. However, their ac- year 2025 year 2030 year 2035 Figure 10. Predicted increase of the primary eclipse depths for the V907 Sco inner binary in the next decade. The effect is due to increasing inclination i1 of the orbital plane (Fig. 6).
The red line is adjusted to the actual measurements in 2022 reported in Sec. 2. The abscissa gives the phase of the P1 ≃ 3.776 day cycle, the ordinate is magnitude.
curacy is obviously low. Modern instrumentation, even of a small aperture, may provide much better results if the star is purposely targeted. Additionally, the situation of V907 Sco is fortuitous since the inclination i 1 will be increasing to a nearly 90 • value for the next decade and the system will be eclipsing until the year ∼ 2049 (Fig. 6). Therefore the upcoming years may offer a chance to significantly improve our understanding of this unique system. As an example, Fig. 10 shows prediction of the increasing depth of primary eclipses during the next decade, a trend reflecting the increase of the inclination i 1 from Fig. 6. This progression holds the promise for i 1 to be determined more precisely in the coming years and for the results of this paper to be improved upon rapidly. The unfortunate miss of the outer orbit radial velocity minimum (Fig. 5) in our work strongly motivates further efforts to observe the V907 Sco system spectroscopically. An intense campaign around the epoch of this minimum (the nearest chance in May 2023) would be especially valuable. These data will help to pin down values of both outer orbit period P 2 and its eccentricity e 2 . Additionally, if the projected semimajor axis of the outer orbit a 2,proj is better constrained from these future data, the analysis developed in this paper may readily provide much better constraints of the mass m 2 .
We also give a brief look at expected, but not detected yet, eclipse time variations (ETVs) in the inner binary of V907 Sco. As the eclipse depth will continually increase in the future (Fig. 10), determination of the mid-eclipse epochs will improve, thus permitting the discovery of variations with respect to the mean period P 1 . If accurate enough, proper interpretation of ETVs may help to further constrain the physical parameters of the system. There are two kinds of ETVs, and they are distinguished by their physical origins. First, the light-time (Roemer) effect is due to the finite speed of light and variations of radial distance (projected along the line of sight) to the eclipsed component in the binary due to its motion about the barycenter of the whole triple system. Second, the dynamical (physical) effect has to do with gravitational perturbation of the third star of the mean motion of stars in the eclipsing binary. Analytical description of the former is rather straightforward, while the same for the latter requires a more involved technique of the perturbation theory to be used. Over the past decade or so, a thorough analysis has been developed by Borkovits and collaborators (e.g., Borkovits 2022, and references therein). Here we specifically used two components of the dynamical ETV signal given in Eqs. (8) and (9) in (Rappaport et al. 2013), and we verified that additional parts are very small. To make sure that this formulation is sufficient, we also determined the dynamical part of the ETVs using direct numerical integration. This is basically the same simulation as shown in Fig. 6, but here we paid special attention to the epochs of eclipses of the two stars in the inner binary.
Results are shown in Fig. 11. We note the dynamical effect dominates by a factor of ∼ 5 over the light-time effect. This is because the system is compact (P 2 /P 1 small enough) and the mass of the third star not too small (m 2 /M 2 large enough). The full amplitude of the dynamical effect, more than 10 minutes, appears large enough to be safely detectable. Especially in the near future when the eclipses will be deep enough, and also the ground-based data would provide adequately precise times of eclipses (which is not the case nowadays, and the number of TESS data is still too limited).
At this moment, we did not attempt to further constrain our solution for the orbital architecture of V907 Sco by using radial velocities published by Lacy et al. (1999). In principle, this is possible by simple propagation of our numerical model back to the epochs of the spectroscopic observations from Table 3 in Lacy et al. (1999), though one would need to assume some uncertainty of the reported radial velocity values. We plan to return to this issue after the system is better constrained by well controlled and accurate photometric and spectroscopic observations in the next year. In particular, V907 Sco will be in the field of view of TESS again in Sector 66 of its operations. This will provide Figure 11. Predicted value of ETVs (∆T ) for the primary component in the V907 Sco. The abscissa is the cycle count of the eclipsing binary (period P1 = 3.776 day; the whole timespan shown is therefore little more than two years), the ordinate is in minutes. The red line is the dynamical (physical) delay due to the third-star perturbations of the mean motion of the binary, the blue line is the light-time (Roemer) effect due to motion of the binary center of mass about the barycenter of the whole system. Both are expressed analytically (e.g., Rappaport et al. 2013). Symbols show the dynamical delay determined numerically for sake of comparison. The perihelion of the outer orbit motion occurs where the gradient of the dynamical effect change is the largest.
very accurate photometric data. Second, we now have a fairly good constraint when the system will pass through the pericenter of the outer orbit. We can thus make efforts to give a high priority to observation during that time period. Such targeted spectroscopic observations will help to improve parameters of the outer orbit, and consequently to constrain the mass m 2 of the third star in the system.
Assuming our solution for m 2 ≃ 1.06 M ⊙ is correct, one may question the role of the radiative flux of the third component in the V907 Sco system in our analysis of the lightcurves during eclipses of the inner binary (the so called third light problem). Given the masses m 1a and m 1b in the range 2.54 to 2.74 M ⊙ , we can estimate the third light contribution at the level 5 %. This is in a good accordance with the TESS lightcurve analysis in Sec. 2, which indicated the third light should be around this value. In the same way, the small ≃ 0.008 eccentricity of the inner binary is close to the forced value by the companion star of approximately solar mass and the orbital parameters listed in Table 3 (see, e.g., Georgakarakos 2009, end references therein).
To conclude with a speculative comment, we note that there is a field star UCAC2 17118369 = 2MASS J17565485-3444584 separated from V907 Sco by ≃ 9 ′′ (possibly the object shown by Lacy et al. 1999, in their Fig. 1). Interestingly, the WDS catalogue records the two stars to be a proper-motion pair (Mason et al. 2001). More importantly though, the Gaia DR3 provides (i) the same distances to both objects within their uncertainty, and (ii) confirms very similar values of their proper motion (though in this case slightly outside their uncertainty overlap, Gaia Collaboration et al. 2022). A detailed look into the putative gravitational bound between the two stars is beyond the scope of the present paper. Here we only mention that (i) their ≃ 3.5 magnitude difference contributes insignificantly, if any, to our TESS photometry analysis in Sec. 2, and (ii) the large distance of the companion cannot have influence on our analysis of the V907 Sco triple dynamics.
We would like to thank an anonymous referee for his/her helpful and critical suggestions improving the level of the manuscript. We thank T. Borkovits for drawing our attention to the TESS 2019 observations of the V907 Sco eclipses onset, P. Harmanec for useful discussions and checks, and the referee, whose suggestions helped to improve the final version of this paper. The work of DV was partially supported by the Czech Science Foundation (grant 21-11058S). BB was supported by National Science Foundation grant AST-1812874. The DASCH project at Harvard is grateful for partial support from NSF grants AST-0407380, AST-0909073, and AST-1313370. Funding for APPLAUSE has been provided by DFG (German Research Foundation, Grant), Leibniz Institute for Astrophysics Potsdam (AIP), Dr. Remeis Sternwarte Bamberg (University Nuernberg/Erlangen), the Hamburger Sternwarte (University of Hamburg) and Tartu Observatory. Plate material also has been made available from Thüringer Landessternwarte Tautenburg, and from the archives of the Vatican Observatory. We would also like to thank the Pierre Auger Collaboration for the use of its facilities. The operation of the robotic telescope FRAM is supported by the grant of the Ministry of Education of the Czech Republic LM2018102. The data calibration and analysis related to the FRAM telescope is supported by the Ministry of Education of the Czech Republic MSMT-CR LTT18004, MSMT/EU funds CZ.02.1.01/0.0/0.0/16 013/0001402 and CZ.02.1.01/0.0/0.0/18 046/0016010. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.