© The Authors 2024

1. Introduction

Before the epoch of the Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2015) missions there were very few well-studied compact hierarchical triple star systems. Here we define ‘compact triple’ as a triple star system with an outer period of ≲1000 days, down to the currently shortest known outer period of 33.0 d (λ Tau; Ebbighausen & Struve 1956). Thanks to the precision space-based photometry of these missions, as well as Gaia astrometry and radial velocity measurements (Gaia Collaboration 2016; Babusiaux et al. 2023), and the Optical Gravitational Lensing Experiment (OGLE) ground-based survey (Udalski et al. 1997), we now know of some thousand compact triple systems (see, e.g., Rappaport et al. 2013; Borkovits et al. 2015, 2016; Czavalinga et al. 2023a; Hajdu et al. 2019). Among the best studied of these are about 50 that have years-long eclipse timing variation (ETV) curves and/or third body eclipses where the tertiary star eclipses the inner binary or vice versa (Borkovits et al. 2015; Rappaport et al. 2022, 2023).

In general, the triple systems with third body eclipses end up providing the most detailed information about their respective triple systems, including a complete understanding of the stellar parameters (masses, radii, T_eff, and ages) as well as the orbital parameters, including the inclination angles of the two orbital planes, and the mutual inclination angle. This detailed information can often be obtained without the need for extensive radial velocity (RV) measurements. The reasons for this success stem from the fact that triply eclipsing triples tend to be quite compact, with orbital periods typically in the range of 40–300 days. In turn, the compactness of the triples can lead to interesting and important dynamical interactions, including so-called ‘dynamical’ delays, driven apsidal motion in the inner binary, possible precession of the orbital planes, as well as some new and fascinating anomalous dynamical interactions, (see, e.g., Rappaport et al. 2013; Borkovits et al. 2015; Borkovits & Mitnyan 2023).

The items that allow for robust system solutions involve the following five ingredients: (1) precision space-based photometry of the inner binary eclipses and the third-body eclipses; (2) eclipse timing variation curves (ETV) which are distinct from measuring the shapes and depths of the eclipses; (3) archival ground-based photometry that has the potential to determine the outer orbital period via the third-body eclipses; (4) archival spectral energy distributions (SED); and (5) a comprehensive complex photodynamical code which puts all this information together in a global analysis. Of course, supplemental RV data are always useful for further reducing the uncertainties and are much appreciated.

Knowing the system parameters is important in trying to understand the complex formation scenarios of triple and multiple star systems (see, e.g., Tokovinin & Moe 2020), as well as for investigating expected and unexpected dynamical interactions in the present-day systems. The latter category includes dynamical interactions on timescales intermediate between the binary and outer period (see, e.g., Appendix A1 of Rappaport et al. 2022), and between the outer period and the apse-node long term period (Pribulla et al. 2023).

In this work we report on seven new triply eclipsing triples. In Sect. 2 we discuss how these systems are found, and the basic properties of the seven systems presented in this work. Segments of the TESS light curves exhibiting third body eclipses are shown in Sect. 3 along with model fits. The folded outer orbital light curves obtained from archival ground-based data are presented in Sec. 4. In Sect. 5 we show ETV curves extracted from the TESS light curves, as well as RV curves for a few sources where we were able to acquire them. The photodynamical analysis code is described briefly in Sect. 6. The individual systems are discussed in Sect. 7. Finally, in Sect. 8 we summarize our work and draw some broader conclusions.

2. The seven triply eclipsing triples

2.1. Discovery of the triples

Most of the triply eclipsing triples reported by us over the past several years (see, e.g., Borkovits et al. 2020b,a; Rappaport et al. 2022, 2023) were found by our ‘Visual Survey Group’ (VSG; Kristiansen et al. 2022). This involved visual inspection of millions of TESS light curves. The vast majority of the light curves were produced from TESS full-frame images (FFI) of anonymous stars down to TESS magnitude T ≲ 15. However, most of the discovered triply eclipsing triples were actually found among a much smaller number (∼1 million) of light curves of preselected eclipsing binaries (EBs). The latter were identified in the TESS data via machine learning (ML) searches¹ (see Powell et al. 2021) of more than 100 million TESS FFI light curves from Sectors 1–48. We made our own light curves from the raw FFIs on the NASA Center for Climate Simulation (NCCS) Discover supercomputer using the ELEANOR (Feinstein et al. 2019) code. The search of these ML-selected binaries was therefore the more productive of the two approaches, but such searches are less likely to discover new classes of objects because we are specifically selecting binaries to look at.

The visual surveying of the light curves is done with Allan Schmitt’s LcTools and LcViewer software (Schmitt et al. 2019), which facilitates the inspection of a typical light curve in just a matter of a few seconds. Our searches for triply eclipsing triples involve looking for an eclipsing binary light curve with an additional extra eclipse that is typically strangely shaped and of longer duration than the regular EB eclipses, or a rapid succession of isolated eclipses. Once such ‘extra’ eclipses, or ‘third body’ events are found, there is usually very little doubt that they are due to a tertiary star eclipsing the EB or vice versa.

From searches through light curves obtained from the first three full years of TESS observations, we have found more than 70 of these triply eclipsing triples. Of these 70 triply eclipsing systems, we have been able to determine the outer orbital period for about half of them via a combination of TESS and archival ground-based photometry, and, occasionally, from Gaia-determined orbits. We have previously reported five of these 30 systems in Borkovits et al. (2020b,a), and Mitnyan et al. (2020) six more in Rappaport et al. (2022); and an additional nine systems in Rappaport et al. (2023). Here we present a detailed analysis and study of seven of the remaining triply eclipsing triples from among this set. In particular, we have chosen those triples which are not scheduled for additional TESS observations, at least not until the end of the ongoing second extended mission. Hence, it is not surprising that six of our seven targets are located in the southern hemisphere.

2.2. Summary of available data

Fortunately, there exist excellent archival data for all of these sources, in addition to the TESS data. In particular, we made use of archival data from (1) Gaia Data Release 3 (Gaia Collaboration 2023); (2) the Mikulski Archive for Space Telescopes (MAST)², in particular the TESS Input Catalog v8.2 (Paegert et al. 2021); (3) the All-Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017); (4) Asteroid Terrestrial-impact Last Alert System (ATLAS³; Tonry et al. 2018; Smith et al. 2020); and (5) the online VizieR SED viewer (Ochsenbein et al. 2000)⁴. Last but not least, (6) for one target, TIC 321978218, we also used SuperWASP data (Pollacco et al. 2006). This relatively bright system was observed with SuperWASP cameras during five consecutive seasons between 2010 and 2014, and these observations partly cover 11 third-body eclipses, which we have included in our analysis (see below). We note that SWASP observations are also available for another target, TIC 176713425, but in this latter case, we did not find these data to be useful for our analysis and, hence, neglected them.

The Gaia and MAST data were used primarily to construct a table of the basic properties of the seven systems reported in this work (Table 1 with, e.g., coordinates, magnitudes, distances, proper motions, and so forth). Amongst the Gaia data we also found the outer orbits for two of our systems. We used the ASAS-SN and ATLAS data primarily to obtain an independent and accurate determination of the outer orbital period for 6 of the 7 systems via the third-body eclipses over a baseline of ∼10 yr. The SED data were obtained from VizieR; these points were used in conjunction with other data to make some of the initial fits to the system parameters.

The sectors in which the sources were observed by TESS, and the sectors during which third-body events were actually seen, are summarized in Table 2. In Table 3 we point out the two sources whose outer orbits were found with Gaia⁵ and references to another source that had previous reports in the literature, and we list the outer orbital parameters of period, eccentricity, and argument of periastron we find from this work.

In addition to the above archival data, we were able to acquire new RV data for two of the sources with the CHIRON spectrometer mounted on the 1.5-m SMARTS telescope at CTIO in Chile (Tokovinin et al. 2013). Spectra with a resolution of 25 000 were taken in the service mode and processed by the standard pipeline (Paredes et al. 2021). We note that a set of archival RV data were also collected for a third target – TIC 185615681.

Finally, in Table 4 we summarize the information we have available for all of the seven sources. This includes eclipse timing variation (ETV) curves that exhibit non-linear behavior, RV data, and spectroscopic outer orbits from Gaia NSS catalog (in two cases).

3. Light curve and model Fits

In Figs. 1–7 we show short segments of the TESS light curves that exhibit third body outer eclipses where the tertiary star eclipses the inner eclipsing binary or vice versa. The blue points are the TESS photometric measurements. For the light curve modeling we used 30-min cadence data even in those sectors where shorter cadence data were available. Hence, all these blue points are at a cadence of 30-min. The regular binary eclipses are either self-evident (or called out in the respective figure captions), while the ‘extra’ eclipses are, for the most part, all the dips in flux that cannot be ascribed to the regular EB eclipses. In some cases the third-body eclipses are quite irregular and anomalous looking, while in other cases, the eclipses look almost ‘normal’ but occur in rapid succession and/or are clearly out of place. Moreover, in the eleven panels of Fig. 8 we plot those third-body eclipses of TIC 321978218 which were observed, at least in part, during the SuperWASP survey, and were used in our analysis.

Fig. 1. Light curves (blue points) and model fits (smooth red curves) for the third-body eclipses of TIC 133771812. The third body eclipses are marked with vertical arrows. The sector numbers are indicated in the lower left corner of each panel. Calculation of the model curves will be discussed in Sect. 6.

Fig. 2. Light curve (blue points) and model fit (smooth red curves) for the only third-body eclipse of TIC 176713425 observed with TESS. The third body eclipses are marked with vertical arrows. The other notation is the same as in Fig. 1.

Fig. 3. Light curves (blue points) and model fits (smooth red curves) for the third-body eclipses of TIC 185615681 (marked with vertical arrows). The other notation is the same as in Fig. 1.

Fig. 4. Light curves (blue points) and model fits (smooth red curves) for the third-body eclipses of TIC 287756035 (marked by horizontal arrows). The other notation is the same as in Fig. 1.

Fig. 5. Light curves (blue points) and model fits (smooth red curves) for the TESS-observed third-body eclipses of TIC 321978218 (marked by horizontal arrows). The other notation is the same as in Fig. 1.

Fig. 6. Light curves (blue points) and two kinds of model fits for TIC 323486857. The light-gray fit represents the pure three-body model light curve, while the red fit is the net fit of the triple star model and a four-frequency Fourier model of the light curve distortions (see text for details). The region of the third-body eclipses is indicated with the horizontal arrows. The other notation is the same as in Fig. 1.

Fig. 7. Light curves (blue points) and model fits (smooth red curves) for TIC 650024463. The vertical and horizontal arrows mark the times of third body eclipses. The other notation is the same as in Fig. 1.

Fig. 8. Eleven partially observed third-body eclipses of TIC 321978218 from the five seasons of the SWASP project. Gray dots represent the individual observations, while larger blue circles denote their 10-min. averages which were used for the joint photodynamical analysis. Red curves stand for the best-fit model.

The smooth red curves are models taken from the full photodynamical analysis. These will be discussed in Sect. 6.

These third-body eclipses contain crucial information about both the orbits (in particular, the mutual inclination angle and outer orbital period) as well as the properties of the stars themselves (e.g., the relative sizes and effective temperatures).

4. Archival outer orbit folds

For each of the seven triply eclipsing triple systems, we looked up the archival ASAS-SN (Shappee et al. 2014; Kochanek et al. 2017) and ATLAS (Tonry et al. 2018; Smith et al. 2020) data sets that were available. The goal here is to find additional third body eclipses in the archival record and thereby possibly determine the outer orbital period of these systems. In general the ASAS-SN data typically yield between 3000 and 6000 photometric points spanning the past decade (over the whole sky). The ATLAS data have only a somewhat smaller number of points than ASAS-SN in the northern hemisphere, but often ≲1000 points in the south. The ATLAS photometry is superior for the fainter targets (e.g., ≳15th magnitude).

We median normalize the V and g bands of ASAS-SN to each other, and likewise the o and c bands of ATLAS. Finally, we median normalize the ASAS-SN and ATLAS data to each other to form a data set typically consisting of between 4000 and 8000 archival points. Next, we remove the EB orbital light curve by subtracting out typically 100 harmonics of the EB period from the raw data (see, e.g., Powell et al. 2021). We then use a Box Least Squares (BLS) algorithm⁶ (Kovács et al. 2002) to do a blind search for the third-body eclipses of the outer orbit covering periods from 20 to 1000 days. These third body events usually occur over an interval of about a day. The outer periods in this sample of seven sources range from 41 days to a year. Therefore, the ‘duty cycle’ for the third body events ranges from ∼0.02 to ∼0.003. For an archival sample of, for example, 5000 points, we can then expect somewhere between 15 and 100 points to land within the outer orbital phase range of the third body events, and thereby measure an average drop in flux at that phase. Thus, for systems with longer outer periods, the task of uniquely identifying the outer orbital period from the archival data becomes ever more challenging.

The results of our search for the outer orbital periods are shown in Fig. 9. For six of the seven sources (all except TIC 133771812) there is a clear BLS signal at the outer period. For TIC 133771812, with an outer period of 244 days, and third-body eclipse depths of only ∼6%, we were not able to detect the outer eclipses at a significant level. On the other hand, we were able to robustly detect the outer eclipses in TIC 287756035 with a year-long outer period. The only difficulty with that one is that due to the annual observing schedule of the ground-based observations, we cannot tell whether the outer period is 368 days or half that. However, other information that is available allows us to break this degeneracy and select the 368 day period.

Fig. 9. Archival outer orbit folds. The ASAS-SN and ATLAS archival data were used to find the outer period of these triple systems independently of what we were able to learn from ETV or RV curves, or third-body eclipses observed with TESS. Typically, there are 3000–6000 archival photometric points spanning a decade. In the case of TIC 287756035, there is about 1/4 of the outer orbital phase which is missing due to the fact that the outer period nearly exactly matches that of the Earth year, and the consequential observing seasons. The light red curves are fits to a non-physical function consisting of a modified hyperbolic secant (e.g., Eq. (1) of Rappaport et al. 2016). The fit covers only one of the two cycles for clarity in viewing the eclipses, and is used only to find the phases of the eclipses.

We note that four of the seven triples studied in this work (TICs 185615681, 321978218, 323486857, and 650024463) exhibit both types of outer eclipses (i.e., secondary vs. primary) in the TESS data – thanks to our nearly edge-on view of their outer orbits. Additionally, in the cases of TICs 176713425 and 287756035, although TESS has observed only one type of outer eclipse, our photodynamical analyses (see later in Sect. 6) predict the existence of the other type of outer eclipse, as well. In contrast to these, only TIC 321978218 exhibits detectable third body primary and secondary eclipses in our archival folded light curves (Fig. 9).

In general, these outer eclipses we see in the archival folds are detected at the ∼8–20σ confidence level. Therefore, if the secondary outer eclipses have several times less area under their normalized light curve, then they can become very difficult to detect. In fact, we do not detect the secondary eclipse in five of the six above-mentioned sources. In four of the six cases, it is obvious that the shallow secondary eclipses would be undetectable in the archival data. In the specific case of TIC 650024463, where the secondary outer eclipses have ∼2/3 the depths of the primary eclipses, one might wonder whether these should have been detected in the archival fold for that source. However, when we take into account that the primary outer eclipses have a duration that is ∼2–3 times longer than the secondary outer eclipses, then the expected detection significance of the secondary eclipses will drop by about a factor of 3–4. This is a sufficient decrease so we cannot expect to detect it robustly–and we do not. Finally, in the case of TIC 287756035, the above mentioned seasonal gaps in the ground-based observations would explain the absence of the archival fold detection of the other type of third-body eclipse⁷.

5. ETV and RV curves

The eclipse timing variation curves form a crucial piece of input information for the comprehensive photodynamical analysis we discuss in Sect. 6. These are based on eclipse midtimes for both the primary and secondary eclipses of the inner EB of the system. We extract these mid-eclipse times in the manner discussed previously in Borkovits et al. (2015, 2016). The non-linear behavior we see in the ETV curves comes predominantly from three basic effects, as follows. First, there is the classical light-travel-time effect (LTTE; Roemer 1677) due to the changing distance to the EB as it is pulled around in its outer orbit by the tertiary star. Second is the ‘dynamical’ or ‘physical’ delays which, in nearly coplanar orbits, are caused largely by the lengthening of the EB orbit due to the presence of the tertiary star (see, e.g., Rappaport et al. 2013; Borkovits et al. 2015). This effect depends on the instantaneous separation between the EB and the tertiary star, and therefore varies with the phase of the outer orbit if that orbit is eccentric. (3) Finally, there is the well-known apsidal motion exhibited by eccentric binaries. This is a long timescale effect which has three different origins: (i) non-spherical mass distributions of the tidally distorted binary components; (ii) effects of general relativity and; (iii) perturbations forced by the tertiary star. The timescale of this latter effect, dynamically driven apsidal motion, is on the order of $P_{\text{out}}^2/P_{\text{in}}$ (Söderhjelm 1975), where P_out and P_in are the outer and inner binary periods, respectively. The two dominant drivers of apsidal motion in the systems reported here are due to the mutual tidal deformations of the two EB stars and forced apsidal motion by the third body.

We see all three of these effects, namely, LTTE, dynamical delays, and apsidal motion, in the ETV curves of the seven triple systems reported in this work. The curves themselves are shown in Figs. 10–15, while the mid-eclipse times used for the derivation of these ETV curves are tabulated in seven separate tables in Appendix A.

Fig. 10. TESS primary and secondary ETV curves (red and blue circles, respectively), with the best-fit photodynamical solution (small red and blue points) for TIC 133771812 (see Sect. 6). The horizontally centered black curve represents the pure LTTE contribution. Vertical lines mark the times of the observed outer eclipses (dashed green-the binary occulting the tertiary star).

Fig. 11. ETV curves formed from the TESS and ground-based follow up observations with the best-fit photodynamical solution for TIC 176713425. The horizontally centered black curve represents the pure LTTE contribution. Dashed vertical green line marks the position of the only (two-dipped) third-body eclipse event which was observed with TESS.

Fig. 12. TESS ETV curves with the best-fitted photodynamical solution for TIC 185615681 (TZ Pyx). The horizontally centered black curve represents the pure LTTE contribution. Vertical dashed green lines represent the times of the observed outer eclipses when the binary occults the tertiary star, and vice versa for the dot-dashed red lines.

Fig. 13. ETV curves for TIC 321978218 with the best-fit photodynamical solution superposed. Left panel: ETV curve derived from TESS data only. We note that only those sections of the ETV curves, where TESS data are available, are shown. Vertical black lines denote the different sector boundaries; dashed green lines represent the times of the observed outer eclipses when the binary occults the tertiary star, and vice versa for the dot-dashed red lines. Right panel: overall ETV curve after adding average ‘seasonal’ ETV points (see text for explanation) calculated from SWASP observations. The longer-term ETV curve also nicely shows an additional variation. We model it with an fourth stellar component, making the system a 2+1+1 quadruple. See text for details.

Fig. 14. TESS ETV curves with the best-fit photodynamical solution for TIC 323486857. We note that only those sections of the ETV curves, where TESS ETV data are available, are shown. Vertical lines have the same meaning as in Fig. 13.

Fig. 15. TESS ETV curves with the best-fitted photodynamical solution for TIC 650024463. The horizontally centered black curve represents the pure LTTE contribution. Vertical lines have the same meaning as in Fig. 13.

TIC 133771812. The ETV curve for TIC 133771812 with a 12-day EB and 244-day outer orbit (see Fig. 10) shows the sparse TESS sampling–essentially two sectors of data every two years. The smaller red and blue dots, connected with straight lines are the model fits (described in detail in Sect. 6) showing the periodic dynamical delays on the timescale of the outer orbit, and the convergence between the primary and secondary eclipse times, due to the apsidal motion, which in this case is being driven by the tertiary star. The variations in the ETV curve with the outer period of 244 d are ∼90% dynamical delays and ∼10% LTTE effects. As another way of illustrating the ratio of the LTTE contribution to the overall ETV curve, we superpose the pure LTTE contribution on the ETV curve for TIC 133771812. We then also show the LTTE contribution to the ETV curves for each of the other sources we investigate.

TIC 176713425. Figure 11 shows the ETV curve for TIC 176713425. Here there are 2+1 sectors of data separated by 3 years and, moreover, six ground-based eclipse times obtained within the framework of our photometric follow up campaigns. Fortunately we know the outer orbital period rather well from the archival data (see Fig. 9), namely 52.58 days. Thus, the Sector 16, 17, and 57 data are sufficient to determine satisfactorily the outer orbital parameters. The variations in the ETV curve with the outer orbital period of 52.58 d are ∼90% dynamical delays and ∼10% LTTE effects.

TIC 185615681. The ETV curve for TIC 185615681 shown in Fig. 12 has sparse points taken at roughly two-year intervals. Fortunately, the outer orbital period is also well measured both from the five sets of TESS-observed third-body eclipses and the archival data to be 55.9 d. The ETV variations with the period of the outer orbit are quite evident, as is the rapid apsidal motion forced by the third body. The former variations in the ETV curve are ∼65% dynamical delays and ∼35% LTTE effects. We also note that this is the only triple in our sample where ETV data prior to the TESS mission are also available in the literature (see Sect. 7.2.3).

This is the only system in our sample where the inner EB is a known double-lined spectroscopic binary (SB2) and, hence, RV data for the EB members are available (see Fig. 16) from the literature (Duerbeck & Rucinski 2007) and, we used these data for our analysis.

Fig. 16. RV curves and model fits for TIC 185615681. The thin, horizontal line at V = −5.95 km s−1 represents the systemic radial velocity of the triple star system. If TIC 185615681 were a single binary, the RVs of the two components would intersect each other exactly on that line. The offset of consecutive intersections of the RV curves from this line is due to the revolution of the inner pair around the center of mass of the whole triple system. The RV points were taken from Duerbeck & Rucinski (2007).

TIC 287756035. This is the only system in our sample where the ETV curve, derived from the TESS observations, is of no value. This is the case because the scatter of the measured mid-eclipse times exceeds the full amplitude of the model-ETV by more than an order of magnitude. On the other hand, however, we have obtained five RV points for the tertiary star which provides us with valuable information about the outer orbit in this system (see Fig. 17 and Table B.1 as well).

Fig. 17. RV curve and model fit for the third component of TIC 287756035. The thin, horizontal line at V = −3.9 km s−1 represents the systemic radial velocity of the triple star system.

TIC 321978218. The left panel of Fig. 13 displays the ETV curve of TIC 321978218 as derived exclusively from the TESS observations. As one can see, the 57.5-day period outer orbit is completely covered both during the Year 3 (Sectors 28 and 29) and Year 5 (Sectors 68 and 69) TESS observations. This fact makes the determination of the outer orbital period and the other orbital elements quite robust without the use of the extra information carried by the third-body eclipses. Our analytic ETV fit also revealed the fact that the ETV is dominated by the LTTE effect (∼73%) over the dynamical delays (∼27%). Besides these data, we also obtained 17 RV points of the dominant tertiary star with the CHIRON instrument. These RV data also nicely cover the outer orbit (see Fig. 18 and Table B.2). A combination of the LTTE-dominated ETV curve, which describes the motion of the inner EB around the centre of mass of the triple system, and the RV curve, which reflects the tertiary’s motion, makes the outer pair a quasi-SB2 binary. In turn, one can determine both the outer mass ratio (q_out) and the total (projected) mass of the inner EB (M_Asin³i_out = M_Aasin³i_out + M_Absin³i_out) and the individual (projected) mass of the tertiary (M_Bsin³i_out) component purely from the amplitudes of the ETV and RV curves. However, the actual system turned out to be not so simple. The reason is that, after receiving the Sector 68 observations and then deriving ETV points from these data, we realized immediately that the ETV curve exhibits a further, longer timescale, non-linear variation. Sector 69 data confirmed these findings.

Fig. 18. RV curve and model fit for the third component of TIC 321978218. The thin, horizontal line at V = 1.6 km s−1 represents the systemic radial velocity of the quadruple star system.

Then, we were able to extend our ETV study thanks to the archival SWASP observations. As was mentioned above, five seasons of SWASP observations are available for TIC 321978218. The lesser quality of these observations do not allow us to determine individual times of eclipses with the accuracies we achieve with TESS which would be necessary to study the low-amplitude (≲10⁻³ d) ETV variations. Hence, instead of using individual eclipse times, we formed seasonal folded, binned light curves from each of the five seasons of SWASP observations (using an EB period of P_in = 0.5698127 d). With the use of these seasonal average light curves, we determined five times of seasonal ‘eclipse points’ (both for the primary and secondary eclipses). Each of these results in a representative ETV value which we locate near to the mid-time of the observations used to make the average light curve. Moreover, since these seasonal folded EB light curves essentially average out the 57.5-day outer period ETV variations, we chose a location for the mean ETV point that is not only near the SWASP seasonal midtime, but one that is also near the time when the outer 57.5-day outer ETV value would be near zero. (In practice, this means that the tertiary’s orbital phase, measured from the middle of first TESS-observed third-body eclipse, was around $0\stackrel{\text{p}}{.}17$ or $0\stackrel{\text{p}}{.}73$.) These nominal seasonal ETV points extend the length of the available ETV data by almost 5000 days and confirm that TIC 321978218 exhibits additional period variation(s) (see right panel of Fig. 13).

At this time we have insufficient information to say anything certain about the nature of these variations either qualitatively or quantitatively. In this study, however, we decided to model this extra ETV variation under the assumption that TIC 321978218 is a 2+1+1 type quadruple star system. We emphasize, however, that, though this assumption was necessary for a consistent modeling of the available datasets, there is no guarantee that the parameters of the outermost orbit and the fourth star mean much unless the outermost orbit is ultimately confirmed.

TIC 323486857. The very noisy ETV curve of TIC 323486857 is shown in Fig. 14. The large uncertainties of the points, which make the ETV curve nearly unusable, are due to the highly distorted light curve. Fortunately, in addition to the Gaia spectroscopic solution for the outer orbit, nine individual third body eclipses were detected during the six TESS sectors and, hence, the outer period, which is the shortest among our sample, is well characterized. As can be seen in Fig. 14, the ETV curve is mainly dominated by the LTTE, which is expected for such a nearly flat, doubly circular triple star system.

TIC 650024463. Figure 15 shows an ETV curve for TIC 650024463 that is fairly well populated with data points from 9 sectors of data spanning more than four years of TESS observations. In addition to the well mapped out variations with the outer orbital period of 108 days, the better part of 2/3 of a complete apsidal motion cycle is observed. The ETV variations with the outer period of 108 d are ∼90% dynamical delays and ∼10% LTTE effects.

6. Photodynamical models

All seven of the triple systems in this work have been subjected to a detailed photodynamical analysis⁸ in order to extract all of the system parameters. For the photodynamical analysis we utilize the software package LIGHTCURVEFACTORY (see, e.g., Borkovits et al. 2019a, 2020a, and references therein). This code has been developed over the past decade and has been described in a number of previous papers. The code contains three basic features: (i) emulators for multi-passband light curve(s), the corresponding ETVs, and radial velocity curve (this feature is present whether or not we have actual RV data); (ii) a built-in numerical integrator (a seventh-order Runge-Kutta-Nyström algorithm) to calculate the perturbed three-, or multiple-body orbits, specifically, the coordinates and velocities of the three or more stars in the system; and (iii) a Markov chain Monte Carlo (MCMC)-based search routine for determining the best-fit system parameters, and the statistical uncertainties, as well. The latter feature uses our own implementation of the generic Metropolis-Hastings algorithm (see, e.g., Ford 2005).

LIGHTCURVEFACTORY was developed to analyze multiple star systems, including binaries, triples, and quadruple stars of both the 2+2 and 2+1+1 architecture. The workings of the code, the steps involved in the analysis procedure, and its application to a wide range of multi-stellar systems have been explained in detail (Borkovits et al. 2018, 2019a,b, 2020a,b, 2021; Mitnyan et al. 2020). Specifically LIGHTCURVEFACTORY has been used successfully to analyse compact and wider triple systems, both with and without outer eclipses.

Essentially all the details of how this code was used to analyze the compact triply eclipsing triple systems that were found with TESS, were described in Rappaport et al. (2022). Therefore, here we will provide only a high-level overview of the inputs to the code and the parameters that are either fitted or constrained by the MCMC fit. Altogether, in the case of a hierarchical triple configuration, there are 25–27 system parameters that result directly from the analysis. Specifically, these are all the stellar parameters (9 = 3 × mass, radius, T_eff), all 12 of the elements of the inner and outer orbits, to which, in the case of radial velocity data, one can add the systemic radial velocity as a thirteenth orbital parameter, as well as the 4 system parameters: distance to the source and the interstellar extinction, as well as the system metallicity and age. Finally, one may optionally adjust the amount of any passband-dependent contaminated (extra) light ℓ₄ if it is necessary.

The LIGHTCURVEFACTORY code obviously requires a substantial amount of ‘input’ information in order to allow for the extraction of 26 independent system parameters. This input information can be divided into two basic categories. First, there are the ‘data’. These include: the (i) EB eclipse profiles, (ii) third-body eclipse profile(s), (iii) times of the EB eclipses which are distinct from the shape and depths of eclipses, (iv) archival SED values, and (v) radial velocities (available only for three of our systems). Second, we utilize PARSEC model stellar evolution tracks and isochrones as well as model stellar atmospheres (Bressan et al. 2012). The evolution tracks enable us to find the stellar radius and T_eff for a given stellar mass, age and metallicity, while the isochrones allow us to compute stellar magnitudes in different passbands in order to fit the SED curve. The available input information for each of the seven triples is summarized in Table 4.

The EB and third-body eclipse profiles and the ETV curves used in the photodynamical analysis were taken from the TESS full-frame images (‘FFI’). The photometry on the FFIs was done with Andras Pal’s FITSH package (Pál 2012). In order to save on computational time we binned the 10-min cadence data to 30-min cadence⁹, and dropped the out-of-eclipse regions of these light curves, keeping only the $\pm 0\stackrel{\text{p}}{.}15$ phase-domain regions around the EB eclipses themselves. However, whenever a data segment contains a third-body (i.e., ‘outer’) eclipse, we keep the data for an entire binary period before and after the first and last contacts of that particular third-body eclipse.

As noted in Table 3, we have RV data for only three of the seven sources. This may raise the question of how we are able to derive absolute stellar masses, temperatures, and radii. Here we provide a qualitative answer to this question. There are in fact several pieces of information that involve combinations of the masses. First, the ETV curve contains information about the light travel time effect and that, in turn, is equivalent to an SB1 RV measurement of the outer orbit. We have useable ETV data for six of the seven sources (see Table 3). And, for that one case where we have no ETV curve, there are RV data (for the tertiary). We note that the same information is also encoded into the light curves, as the light curves are fitted in the time-domain. Finally in this regard, the ETV curve (and the light curve, as well) also contain signatures of the dynamical delays which, in turn, carry information about primarily q_out and, in a less certain way (through the higher order perturbations), about q_in (see, e.g., Borkovits et al. 2015).

Additionally, the geometry and timing of the outer eclipses carry significant further information about the ratio of q_out/q_in, as was elaborated on, for example, in Appendix A of Borkovits et al. (2013). It is even intuitive that, while the inner eclipses provide information about R_Aa, Ab/a_in, the outer eclipses provide information about R_{Aa, Ab, B}/a_out and, hence, their combination leads to the above mentioned ratio of mass ratios.

We did not use spectroscopically determined temperatures and metallicities for most of the systems (as they are unavailable) but, instead, we use SED fits to determine or constrain absolute temperatures. This information is also combined with the results of the simultaneous light curve fits, the latter of which also give combinations of the ratios of the surface brightnesses and, hence, indirectly the absolute temperature of the three stars. Primarily these may be obtained from combinations of the eclipse depths of the inner and outer eclipses. Here one should keep in mind that the ratio of the surface brightnesses of the inner EB components can be obtained not only from the mutual eclipses of the inner components, but also from those events where the inner EB members eclipse (or are eclipsed by) the tertiary star separately. Hence from a light curve containing both inner and outer eclipses, the information which can be mined out is not simply the sum of the parameters that can be determined from two independent eclipsing light curves, but much more. The use of SED information to obtain temperatures is also employed in several other recent papers, for example in Miller et al. (2020) and Stassun & Torres (2016).

Finally, we combine this SED information with theoretical, coeval stellar isochrones which provide information not only on the radii and T_eff of the stars, but also the masses for a given age (but of course, such masses are no longer independent of the astrophysical assumptions and, hence, they may be somewhat inferior to those dynamical masses which can be directly inferred from high quality RV data). We discuss this question in detail in Borkovits et al. (2020a, 2022a). All the stars are typically tied together via a “coeval” assumption. And, of course, knowledge of the masses sets the size scales of the system, which then provides for absolute determinations of semi-major axes and stellar radii.

7. System parameters

7.1. Tables of fitted parameters

With the use of the 25–27 directly fitted system parameters, one can also determine a number of other astrophysically and/or dynamically relevant additional parameters for each star in the triple. Thus, we can include all the basic stellar parameters in our tabulated results, as well as the inner and outer orbital parameters discussed above, including the relative orientations of the two orbits. We also give several relative values (e.g., R/a, T_B/T_Aa), some global system parameters such as distance, [M/H], and E(B − V), and a number of parameters that are derived from the fitted parameters. Each of the results tables contains a few dozen different parameters which are, for the most part, written as symbols only. These symbols were defined in Rappaport et al. (2023), but for the sake of completeness we repeat these definitions here in Table 5. Then, the system parameters that are derived from the photodynamical analyses are listed in Tables 6 through 9.

Before discussing the results individually for each system, we make one other important point. The orbital parameters (including also the inner and outer periods) tabulated below are instantaneous, osculating orbital elements, which are valid strictly only for the specific epoch t₀ which is given in the very first row of each table. Therefore, the given periods and conjunction times are not applicable for observational predictions (i.e., for determining times of future inner or outer eclipses). For these latter purposes we provide Table 10 which contains, among other things, the inner and outer eclipsing periods of the sources, and should be used for planning occasional future observations.

7.2. Results for the individual systems

7.2.1. TIC 133771812

This is the only target in our sample for which we were unable to obtain a reasonable third-body orbit fold from the archival data. Thus, in order to determine the outer period we had to restrict ourselves to just the three sets of third-body eclipses observed with TESS in years 1, 3 and 5, and to the high quality ETV data obtained from these observations. Our analytic ETV solutions preferred an outer period of P_out ∼ 244 d, but allowed also for half that period, namely, P_out ∼ 122 d. We produced preliminary photodynamical solutions for both outer periods, and we found that the true period should actually be about P_out ∼ 244 days. In other words, we were unable to find both dynamically and astrophysically consistent and reliable solutions with the shorter outer period. We have tabulated the median posteriors of the adjusted parameters and several further derived quantities of the complex photodynamical analysis in the first three columns of Table 6.

According to our results, all three components of TIC 133771812 are somewhat more massive than our Sun. The inner EB consists of the most and least massive (M_Aa = 1.58 ± 0.04 M_⊙ and M_Ab = 1.20 ± 0.03 M_⊙, respectively) and, correspondingly, the hottest and coolest (T_Aa = 7700 ± 160 K and T_Ab = 6410 ± 100 K) stars of the triple. The third component lies in between the two EB members both in its mass and temperature, but is closer to the secondary of the inner EB with M_B = 1.25 ± 0.03 M_⊙ and T_B = 6550 ± 110 K. The inferred masses, in the absence of RV data, have relatively large uncertainties of ∼3 − 4%. Relative quantities, such as the mass ratios, however, can be deduced with better accuracies. These are found to be q_in = 0.763 ± 0.006 and q_out = 0.449 ± 0.010.

Amongst all the triply eclipsing triple stars which were analyzed photodynamically in our present and former papers, TIC 133771812 has by far the longest inner EB period with $P_{\mathrm{in}}=12\stackrel{\text{d}}{.}33$. Thus, despite the fact that the inferred outer period ($P_{\mathrm{out}}=243\stackrel{\text{d}}{.}9$) is the second longest in the current set, this target is still the second tightest (P_out/P_in ≃ 19.8) among these seven triply eclipsing triples. Due the relative large size of the inner EB and, hence, the small fractional radii (r_Aa = R_Aa/a_in = 0.0491 ± 0.0005 and r_Ab = 0.0368 ± 0.0008), tidal effects are certainly negligible. Hence, not surprisingly, this inner EB has the largest eccentricity with e_in = 0.0337 ± 0.0004 among our present sample. The combination of inner and outer (e_out = 0.217 ± 0.006) eccentricities, combined with the tightness of the system, still leaves it safely within the dynamically stable region according to the Mardling & Aarseth (2001) stability criterion.

7.2.2. TIC 176713425

This is the only northern (ecliptic) hemisphere target amongst the currently investigated triples. Its original FITSH light curves were partially blended with the nearby overcontact binary ATO J353.9138+42.3630 (TIC 176713436). Hence, to eliminate the signal of the latter object from the light curve of our target, we applied a principal component analysis (PCA) method implemented in the lightkurve (Lightkurve Collaboration 2018) Python package and its dependencies: astropy (Astropy Collaboration 2018), astroquery (Ginsburg et al. 2019) and TESScut (Brasseur et al. 2019).

We achieved this by creating 60 × 60 pixel size cutouts around the target and choosing an aperture such that the contaminating star also has some contribution in the image region outside the aperture. After that, we created a design matrix containing regressors from all pixels surrounding the aperture and calculated the principal component vectors using lightkurve’s built-in PCA method. These principal components are a combination of signals from scattered light, spacecraft motion and the nearby contaminating star. As a final step, we used the built-in RegressionCorrector method of lightkurve that uses linear algebra to find a combination of these principal component vectors that will transform the input light curve closest to zero. As the mean flux level of the light curve is different from zero, an additional offset term has to be used by appending a constant term to the design matrix. In this way, all of the previously mentioned signals were modeled and eliminated from the light curve of the target.

Due to the northerly position of the target, we were able to obtain photometric follow up observations in Hungarian observatories, and thereby observed eight additional regular eclipses (four primary and four secondary minima)¹⁰ with the two identical 80 cm RC telescopes of Baja Astronomical Observatory and Gothard Astrophysical Observatory, Szombathely (for short descriptions of the instruments used see Borkovits et al. 2022a). Despite the fact that TESS observed only one set of third-body events, the parameters of the outer orbit can be determined robustly not only due to the archival observations which provide the outer orbital period (P_out), but also because of the very high quality ETV data (Fig. 11). In this regard, we note that the ground-based follow up eclipse times were found to be very helpful in excluding some possible alternative outer periods.

Turning to the astrophysical implications of the photodynamical results (see last columns of Table 6), we found that the system consists of nearly identical triplet stars. The masses of the EB primary and the tertiary agree to within their 1-σ statistical uncertainties (M_Aa = 1.17 ± 0.03 M_⊙ vs. M_B = 1.16 ± 0.03 M_⊙), and the secondary of the inner EB has a lower mass by only 6 − 7% (M_Ab = 1.09 ± 0.03 M_⊙) (referring to the more precise relative quantities, we find: q_in = 0.933 ± 0.005 and q_out = 0.515 ± 0.004, respectively). What makes this system more interesting from an astrophysical and, especially, evolutionary point of view is that the two more massive components (stars Aa and B) are very likely at the end of their TAMS or, just at the beginning of their evolution from the MS toward the giant branch. This means that, although these stars are more massive, the slightly less massive secondary actually looks a bit hotter, though, the small difference is well within the 1-σ posterior uncertainties and, thus, looks statistically insignificant ($T_{\mathrm{Aa}}={6033}_{-192}^{+78}$ K and $T_{\mathrm{B}}={6033}_{-160}^{+97}$ K vs. $T_{\mathrm{Ab}}={6051}_{-153}^{+77}$ K). At the same time, despite the quite similar masses of the inner EB members, the radii of the EB components differ substantially, being R_Aa = 1.555 ± 0.015 R_⊙ vs. R_Ab = 1.137 ± 0.045 R_⊙ (or, again citing the precise fractional radii, we find r_Aa = R_Aa/a_in = 0.2060 ± 0.0018 vs. r_Ab = 0.1525 ± 0.0025).

These results emphasize that a more robust and precise determination of the stellar masses with the addition of RV observations would be very important, since the inner EB of this triple may be used to strongly constrain stellar evolutionary models. Moreover, another aspect that should be considered is that, because of the nearly equal surface brightnesses of the primary and secondary components of the inner circular EB, one can assume that due to the total eclipses, the secondary-to-primary eclipse depth ratio is dominated by the limb darkening law of the binary components. Therefore, in order to obtain a better characterization of limb darkening law(s) – a recently intensively investigated area, having crucial importance in the precision analysis of transiting exoplanets, or even of eclipsing binaries (see e.g., Csizmadia et al. 2013; Maxted 2023) – high-precision photometric observations with good time resolution of future regular (binary) eclipses would also be very useful.

Turning to the dynamical properties of the triple, the period ratio P_out/P_in = 27.9 indicates that TIC 176713425 is a moderately tight system, where one can expect the effects of significant third-body perturbations. This can be seen nicely in the ETV curve, where the dynamical delays strongly dominate over the pure, geometric LTTE. On the other hand, however, due to the small physical size of the inner binary (and, hence, the significant fractional radii of stars Aa, Ab) the tidal effects are non-negligible and the tidally forced nominal¹¹ apsidal advance rate in the inner EB is of the same order as the dynamically forced one ($\mathrm{\Delta }{\omega }_{\mathrm{in}}^{\mathrm{tide}}={387}^{″}\pm {15}^{″}{\mathrm{cycle}}^{-1}$ vs. $\mathrm{\Delta }{\omega }_{\mathrm{in}}^{3\mathrm{b}}={1230}^{″}\pm {15}^{″}{\mathrm{cycle}}^{-1}$).

7.2.3. TIC 185615681 = TZ Pyx

This is the only system in our current sample of triples where the inner EB has been known already for a long time. The variability of the system was discovered by Strohmeier (1966). Duerbeck & Rucinski (2007) obtained and analyzed 8 RV points for each of the two stellar components of the EB. They found that the inner EB consists of two very similar stars (quantitatively, they found a spectroscopic mass ratio of q^spec = 0.996 ± 0.020) that are hot and likely of late A spectral type. These authors, however, did not take note of the slight eccentricity of the binary’s orbit. It was Otero (2007) who reported for the first time not only the binary eccentricity, but also evidence for the apsidal motion in the EB. This effect was then studied first by Zasche (2012) who, with the use of 16 badly scattered ground-based eclipse times (spanning an interval of ∼7000 days), determined an apsidal motion period of P_apse = 157 ± 37 yr. But, as they note: “The apsidal motion fit is still not very convincing and would have even shorter period [...], but only further observations would confirm this hypothesis.” From this apsidal motion period they inferred that the (mean) apsidal motion or, internal structure constant, k₂, is in accord with the theoretically expected value. Naturally, they were not aware of the presence of the third stellar component in the system, whose perturbations, as our analysis clearly reveals, completely dominates over the tidal effects. Hence, even if their apsidal motion period had been correct, it would not be conclusive in regard to the internal structure constants of the stars. Most recently Kim et al. (2018) repeated the ETV analysis, using the available 20 data points, and got a much shorter apsidal motion period of P_apse = 22.89 ± 0.30 yr, which is much closer to our findings.

Turning to our analysis of the available former ground-based observational material, we used only the RV data of Duerbeck & Rucinski (2007). Neither ground-based light curves, nor ETV data were used for our photodynamical runs. While neglecting the ground-based light curve data is well justifiable given the availability of the far superior TESS light curves, the ETV points, in theory, might be usable to obtain a more precise detection of the apsidal motion period. Unfortunately, however, the historical ETV data have such bad scatter that, from our perspective, these are unusable to mine out any additional information.

The twin components of TZ Pyx (M_Aa = 2.013 ± 0.034 M_⊙; M_Aa = 2.008 ± 0.033 M_⊙) are the hottest, fully radiative stars in our sample. Moreover, this is the only triple in our sample, where the third, outer component is substantially less massive (M_B = 0.888 ± 0.014 M_⊙) than either of the binary members. (Referring again to the more precise relative quantities, we find: q_in = 0.998 ± 0.002 and q_out = 0.221 ± 0.003, respectively.) In such a way, the contribution of the third component to the total system flux is negligibly small, only ∼1%. On the other hand, the depths of the primary third-body eclipses may reach ∼6 − 8%, where, for such a relatively bright (V = 10.7) variable star, such extra dips would be serendipitously detectable even with sub-metre size ground-based telescopes. Despite this fact, no indications of the triple nature of this system have been previously reported. This fact emphasizes that one should not be surprised to occasionally find a triply eclipsing tertiary hidden even amongst the best observed and brightest EBs.

Regarding the orbital properties, we found a small, but significant, eccentricity for the inner orbit, being e_in = 0.01067 ± 0.00009. This value agrees well with the results of Kim et al. (2018) who obtained e_in = 0.01068 ± 0.00008¹². The outer orbit was also found to be slightly eccentric with e_out = 0.0982 ± 0.0008. This is the second lowest outer eccentricity in our current sample.

Despite the moderately tight nature of the triple system (P_out/P_in = 24.2) the full-amplitudes of the primary and secondary ETVs, on the P_out-timescale, remain under 1.5 min (∼0.001 d). This effect is hardly detectable with incidental, highly inhomogeneous eclipse times measured with ground-based photometry. The ETVs, determined from the precise TESS photometry, however, clearly show the presence of typical, third-body induced dynamical ETVs (see Fig. 12). Our solution reveals that the ETVs on this timescale are nearly equally contributed by the third-body perturbations and the geometric LTTE.

On a longer timescale, however, the main effect seen in the ETV curves is due to apsidal motion. In turn, the third-body-driven apsidal motion dominates over the classic tidal component, but the latter contribution is non-negligible ($\mathrm{\Delta }{\omega }_{\mathrm{in}}^{3\mathrm{b}}={738}^{″}\pm 5^{″}{\mathrm{cycle}}^{-1}$ vs. $\mathrm{\Delta }{\omega }_{\mathrm{in}}^{\mathrm{tide}}={144}^{″}\pm 2^{″}{\mathrm{cycle}}^{-1}$). In this regard, we note that even though three inner eccentric EBs of the current sample of seven systems exhibit well observed apsidal motion, TZ Pyx is the only one where the classic tidal effects are important. Hence, this is the only system in our sample where the proper settings of the internal structure constants have real significance. Therefore, we set the first two apsidal motion constants to k₂ = 0.005 and k₃ = 0.001 for both members of the inner EB, in accord with the recent tables of Claret (2023). Taking into account all three different contributions of the apsidal motion (though the relativistic contribution was found to be much smaller than the other two, with $\mathrm{\Delta }{\omega }_{\mathrm{in}}^{\mathrm{GR}}=2\stackrel{″}{.}83\pm 0\stackrel{″}{.}04{\mathrm{cycle}}^{-1}$) our solution gives a theoretical apsidal motion period of $P_{\mathrm{apse}}^{\mathrm{theo}}=18.2\pm 0.2$ yr.

7.2.4. TIC 287756035

This system has, by far, the longest of the outer periods among this current set of seven triply eclipsing triples at 368 days. In fact, this outer period is so close to that of an Earth year, that we were able to see only ∼3/4 of the phases in the outer orbital fold of the archival data (see left panel, middle row of Fig. 9). The system is quite flat with the inner and outer orbital inclinations near 89° and a mutual inclination angle of i_mut = 1.4°.

The outer orbit has a quite typical eccentricity for these systems of e_out = 0.23 ± 0.02, while the observable argument of the outer periastron was found to be ${\omega }_{\mathrm{out}}=53\stackrel{°}{.}7\pm 0\stackrel{°}{.}7$. One should keep in mind, however, that both TESS and the formerly mentioned ground-based surveys observed only one type of third-body eclipse, namely, those where the tertiary star occulted the inner binary. Therefore we are not able to constrain the quantity e_out cos ω_out through the offset of the other type of third-body eclipse from an orbital phase of $0\stackrel{\text{p}}{.}5$. In addition to this fact, due the lack of any usable ETV curve (which would offer additional constraints on e_out and ω_out), our primary source for information about e_out and ω_out is only the poorly covered RV curve of the tertiary (see Fig. 17). Thus, the obtained e_out and ω_out values are less robust than in the case of the other systems studied in this work.

One can, however, make, an indirect, inverse check on the reliability of the obtained outer eccentricity and argument of periastron values. Our photodynamical solution (of which e_out cos ω_out was a direct, adjustable parameter) resulted in a median posterior of e_out cos ω_out = 0.14 ± 0.01. Taking into account that the phase of the superior conjunction relative to the inferior conjunction can be approximated as ${\varphi }_{\sup }\simeq \frac{1}{2}-\frac{2}{\pi }ecos\omega$ (see, e.g., Sterne 1939), one can calculate that the other type of third-body eclipse should occur around phase $0\stackrel{\text{p}}{.}411\pm 0\stackrel{\text{p}}{.}007$. Taking a look at the archival folded light curve of TIC 287756035 (Fig. 9), one can see that this phase value is very close to the left edge of the seasonal gap of the observations, but the target was definitely observable during this orbital phase. Despite this, no dip can be noticed in the fold around that phase. One should keep in mind, however, that this phase is very close to the beginning of the seasonal gaps, and one may expect that the target was observable only for a short time just after sunset and hence, only observations reduced in number and quality are available during these intervals. Another possibility might be that, due to the weaker constraints on e_out cos ω_out, the obtained posteriors might have larger statistical uncertainties than we found. Thus, further RV observations would be very important to clarify this question.

Apart from this issue discussed above, what makes this system very interesting from an astrophysical point of view is that the two more massive stellar components, with very similar masses of M_B = 1.133 ± 0.04 M_⊙ and M_Ab = 1.116 ± 0.04 M_⊙, are evolved to very different levels¹³. The slightly more massive tertiary is quite evolved at R_B = 6.7 ± 0.3 R_⊙, and dominates the light of the system. On the other hand, the secondary component of the inner EB is just at the very beginning of its evolution toward the giant branch with R_B = 2.36 ± 0.06 R_⊙. For such near solar mass stars, they must be fairly old at 5.7 Gyr in order to be ascending the giant branch. Interestingly, this shows what will become of our Sun within a short (astronomical) time. We also note that even the third, less massive star with M_Aa = 0.93 ± 0.02 M_⊙ and R_Aa = 0.98 ± 0.02 R_⊙ is currently located on the TAMS just before its evolution toward the giant branch. Thus, at the current epoch, this star is the hotter one and, hence, it is occulted by its more massive companion in the inner binary during the primary eclipses. Thus, we call star Aa the primary of the inner EB, even though it has the lower mass of the pair. It is somewhat surprising that Gaia has not yet measured the spectroscopic orbit of the tertiary given that it dominates the system light, is reasonably bright at G = 13.5, and has relatively narrow lines with T_eff ≃ 4830 K. Only the tertiary in TIC 323486857 is somewhat more evolved (R ≃ 7.5 R_⊙) among our current set of triples.

The system is the opposite of ‘tight’ with P_out/P_in ≃ 368/2.1 ≃ 175, and therefore, given the flatness of the system and small-to-modest eccentricities for the inner and outer orbits, respectively, no large dynamical effects are expected.

Finally, we note that star Ab, at R = 2.36 R_⊙, is only modestly underfilling its Roche lobe which is currently about 3.45 R_⊙. Thus, since Ab is already well evolved off the main sequence it will overflow its Roche lobe relatively soon. By contrast, the tertiary star, at R = 6.75 R_⊙, underfills its Roche of 104 R_⊙ by a large margin. However, since the tertiary is more evolved than star Ab, and therefore has a big evolutionary lead over Ab, it is an interesting question as to which star will first overflow its Roche lobe. It appears that the tertiary will actually grow to 104 R_⊙ within 100 Myr and overflow its Roche lobe, while star Ab will require about 40% more time to overflow its Roche lobe in the EB.

7.2.5. TIC 321978218

TIC 321978218 has the largest archival observational material among our current sample of triples. Besides the above mentioned four seasons of SWASP observations, this target was also spectroscopically surveyed by the Radial Velocity Experiment (RAVE), and detailed spectroscopic parameters (T_eff, log g, different abundances, etc.) were published in their data release 5 (Kunder et al. 2017). Our own spectroscopic follow up has also resulted in 17 high-quality RV points. The system was found to be single lined in the RAVE analysis, and we also found that only the lines of the tertiary component are observable. These are in nice accord with the findings of the photodynamical analysis, which shows that the system’s light is largely dominated by the tertiary component. Hence, we can conclude that the quantitative spectroscopic results of RAVE basically describe the tertiary component. Consequently, the effective temperature of T_eff = 5974 ± 150 K given in RAVE DR5 yields the effective temperature of the dominant tertiary. This result becomes especially interesting when we combine the preliminary analytic RV and ETV solutions to find an estimate for the mass of this tertiary star. A first analytic solution for the high-quality RV curve (see Fig. 18) yields a spectroscopic mass function of f_sp(M_A) = 0.506 ± 0.003 M_⊙. On the other hand, the preliminary analytic LTTE ETV solution provides a mass function of f_ETV(M_B) = 0.290 ± 0.030 M_⊙. Writing these out explicitly

$$\begin{array}{c}\hfill f_{\mathrm{ETV}}(M_{\mathrm{B}})=M_{\mathrm{B}}{sin}^3{i}_{\mathrm{out}}{\left(\frac{q_{\mathrm{out}}}{1+q_{\mathrm{out}}}\right)}^2\end{array}$$(1)

and,

$$\begin{array}{c}\hfill f_{\mathrm{sp}}(M_{\mathrm{A}})=M_{\mathrm{A}}{sin}^3{i}_{\mathrm{out}}{\left(\frac{q_{\mathrm{out}}^{-1}}{1+q_{\mathrm{out}}^{-1}}\right)}^2=M_{\mathrm{A}}{sin}^3{i}_{\mathrm{out}}{\left(\frac{1}{1+q_{\mathrm{out}}}\right)}^2,\end{array}$$(2)

one can readily see that

$$\begin{array}{c}\hfill M_{\mathrm{B}}{sin}^3{i}_{\mathrm{out}}=f_{\mathrm{ETV}}^{1/3}{(f_{\mathrm{ETV}}^{1/3}+f_{\mathrm{sp}}^{1/3})}^2=1.41\pm 0.10M_{\odot }.\end{array}$$(3)

Thus, one may immediately realize that a star with an effective temperature of T_eff ≃ 5900 K and stellar mass of m_B ≈ 1.4 M_⊙ cannot be a main sequence object (we would expect a range of T_eff ≈ 6500 − 7000 K for this mass). This suspicion became even stronger when we made efforts to model the light curve without taking into account any astrophysical constraints, that is, by directly adjusting the fractional radii of the stars and their temperature ratios. In such a way we found a quite good fit, but the resultant log g’s for the inner EB stars were found to be ∼4.4 together with temperatures of less than 5000 K. In other words, the fractional radii of the EB members¹⁴ (this latter quantity is mainly dictated by the phase-durations of the regular eclipses) were found to be too large.

After an initial consideration, these findings call to mind three plausible explanations; however, two of them can be readily refuted. First, for a main-sequence tertiary with a mass of M_B ≈ 1.8 M_⊙ and corresponding effective temperature of T_eff ≈ 8000 K, a consistent model can be obtained with three coeval MS stellar components, where the total mass of the inner EB is well constrained by the spectroscopic mass function. In this case, however, (i) the ETV fit becomes considerably weaker, (ii) the model system SED clearly lies far above the measured catalogued SED points and, (iii) the observed spectrum is clearly inconsistent with such a hot source. Second, one can assume that the EB members are, again, main sequence stars, but they have radii that are inflated by 10–15% due to tidal interactions with each other (see, e.g., Spada et al. 2013). The problem in this case, however, is that if one assumes the tertiary is a quite evolved¹⁵ ∼1.4 M_⊙ star, the EB components with the proper masses would be far too bright for the observed eclipse depths. This means that for relatively brighter EB components, the deeper third-body eclipses, which occur when the EB stars transit in front of the tertiary would be too shallow, while the other third-body eclipses, when the EB stars are occulted by the third tertiary star, would be too deep. Moreover, the regular EB eclipses would also be deeper than are observed.

After the above arguments, we found only one really plausible explanation for the system properties. According to our interpretation TIC 321978218 is a very young, pre-main sequence system. This assumption would explain both the low temperature of the third star relative to its mass and the large fractional radii of the less massive EB members. An additional problem, however, arises even in this pre-MS scenario. We were unable to find strictly coeval stars which would satisfy all the observational constraints. Instead we allowed a different age for the tertiary component than that of the EB members, the latter two of which were considered to be strictly coeval. As one can see in Table 8, we found satisfactory solutions with a slightly younger tertiary t_B ≈ 11 Myr vs. t_{Aa, Ab, C} ≈ 21 Myr. Such a small difference of ≈10 Myr might be comparable with even the life time of a circumbinary disk, and also may imply the operation of a sequential disk instability mechanism as a possible explanation of multiple stellar system formation (see, e.g., Tokovinin 2021, and further references therein). It may also amount to only some minor inconsistencies in the pre-main sequence stellar evolution tracks.

Having established that TIC 321978218 is a young system, we return to the interesting system structure initially introduced in Sect. 5. The securely understood part of the orbital architecture is that of a flat edge-on triple with a 0.57-day inner EB and a tertiary in a 57.5 day outer period. The mutual inclination angle is only ∼2°. The tertiary contains some 90% of the system bolometric luminosity.

In addition, as we found from the ETV curves in Fig. 18 there is a fairly clear additional long-term structure superposed on the ETV of the triple system itself. If we assume that this additional perturbation to the light curve is due to another body orbiting the triple farther out, we find evidence for a quadruple system with a 2+1+1 hierarchy. If this is the correct interpretation, the outer orbital period is ≈4400 d, the outer eccentricity is ∼0.37, and the mutual inclination with the inner triple is ∼10°.

The only astrophysical parameter associated with the fourth star that is actually constrained by the available observational material is its (very low) mass of M_C = 0.16 M_⊙. Because of its negligibly small contribution to the system’s total flux, as well as the absence of any other effects on the light curve, all other parameters such as stellar radius and effective temperature remain basically undetermined. Due the consistency of the model, however, LIGHTCURVEFACTORY calculates all the astrophysical parameters of this fourth star in the same manner as for the other three components, namely, via the use of PARSEC isochrones. During all the runs, for technical reasons, the age of this fourth star was considered to be equal to the age of the inner EB system (21 Myr). The resultant effective temperature (T_C = 2790 ± 100 K) follows from its mass and age. The highly bloated radius (R_C = 0.49 ± 0.02 R_⊙) can be understood in the long time it takes for a star this low in mass to contract onto the main sequence. However, none of these parameters for the fourth star should be taken too seriously until the existence of the fourth star is confirmed by further observations.

7.2.6. TIC 323486857

This is another target where the TESS light curve was strongly contaminated by nearby, and likely variable, stars. Thus, we had to apply our PCA analysis again in order to separate the signal of the triply eclipsing triple system from other stars. Due to the strongly variable nature of the contaminating light, this was not a simple task and, during the process, we probably lost a significant fraction of the flux coming from our target system. Moreover, even after the PCA process, we obtained a substantially distorted light curve. We assume that this remaining signal is real, and comes from our system, revealing strong chromospheric activity in the late type stars. We made some further efforts to remove, or at least reduce, the effects of these extra distortions to the geometric triple star signal (the latter of which we are investigating), but these efforts met with only partial success. Finally, we carried out our usual complex photodynamical analysis, but simultaneous to the triple-star light-curve fitting, we also fit the extra distortions together with harmonic functions with four frequencies, as we have done previously in some other cases (see Borkovits et al. 2018, for a detailed description of this process). In such a way we were able to obtain a quite reliable solution, but the parameter uncertainties are somewhat larger than in the case of the other triple systems.

TIC 323486857 is comprised of three comparable-mass stars in highly circular inner and outer orbits that are nearly coplanar. The mutual inclination angle between these two orbits is only $2\stackrel{°}{.}5\pm 1°$. The tertiary star has a mass of 1.6 M_⊙ compared to 1.4 M_⊙ for the primary star of the inner binary. Just this small difference in mass, however, is sufficient for the tertiary star to be substantially evolved (R_B = 7.5 R_⊙), while the primary EB star has not evolved much.

This system has three related properties that make it rather interesting. It has a combination of (i) a very short outer orbital period (41 d), (ii) a quite low outer eccentricity (0.0066 ± 0.0040), and (iii) a fairly large tertiary fractional radius r_B = R_B/a_out = 0.093 ± 0.007. In Fig. 19 we plot the essence of this information in the plane of e_out–r_B for all 33 triply eclipsing triples that we have studied recently (see Sect. 8 for more information), including TIC 323486857. With the exception of the two sources at the bottom left of the plot (marked with a faint green ellipse) there is a strong inverse correlation between the outer eccentricity and the tertiary’s fractional radius. We believe this is caused by dissipation in the tides raised in the evolving giants by the companion EB. The most circularized systems (with the two aforementioned exceptions) have outer orbital separations of only ∼7–13 times the radius of the tertiary, while the tidally induced decay of the eccentricity goes as a high power of R_B/a_out (Zahn 1977). This is completely analogous to the circularization process that operates in ordinary binary stars. The other two systems with low eccentricity, but relatively small values of r_B, seem otherwise perfectly normal systems with well measured parameters–we presume these were born circular.

Fig. 19. Parameter eout plotted against the fractional radius of the tertiary component rB = RB/aout for 33 triply eclipsing triple systems that we have studied recently (see Sect. 8 for details). The red line sketches out the inverse correlation between these two parameters, with a functional form of eout ≃ exp(−42rB). Exceptions are the two systems marked with a faint green ellipse with very low eccentricities, but not large values of rB. Numbers next to individual sources are the TIC numbers. We omit error bars on the systems with eout ≳ 0.1 for increased clarity.

7.2.7. TIC 650024463 = Gaia DR3 4618836249918572544

The stellar identification of this system is somewhat problematic. This is because, according to Gaia DR3, there are three potential candidate stars within 10″ of each other (see Fig. 20). These are Gaia DR3 4618836249918572928 = TIC 302965294 with G = 13.28 mag; 4618836249918572544 = TIC 650024463 with G = 14.84 mag at 2″ from the former target; and 4618836249918735488 = TIC 302965293 with G = 15.13 mag, a bit farther at 9″. The three targets have completely different DR3 parallaxes and proper motions, hence, there are no physical relations amongst them, at least not with any certainty. On the other hand, it is also clear that the images of the three stars cannot be disentangled due to the large pixel size of TESS and, naturally, the blended light curve of the three stars is totally dominated by the brightest target star.

Fig. 20. Positional information for TIC 650024463 and its two neighbors. Upper panel: A 4 × 4 pixels TESS image near the position of TIC 650024463 (marked by the red arrow) from Sector 1, highlighting nearby field stars (orange dots). Lower panel: photocenter measurements for the primary eclipses detected in Sector 1 compared to the relative pixel position of TIC 302965293. The small red dots represent the measurements for each eclipse (four in Sector 1); the large red circle represents the average photocenter; and the grey contours represent the corresponding 1-, 2- and 3-σ confidence intervals. As seen from the panel, the measurements clearly rule out TIC 302965293 as a potential source of the detected eclipses.

Thus, our first assumption was that the triply eclipsing triple stellar system is hosted by the brightest stellar image, namely that it belongs to TIC 302965294. The very first light curve fitting trials, however, revealed that a consistent EB and triply eclipsing model can be obtained only with the assumption of 80 − 85% contaminated light. However, for the brightest star, TIC 302965294, it is clearly not possible that this much contaminated light comes from one of the fainter neighboring stars. On the other hand, if the triple resides in one of the fainter stars (i.e., TICs 650024463 or 302965293), then the extra light is nicely accounted for by the brightest star which does not host the triple. These two fainter stars have G_RP band fluxes of ≈0.26 and ≈0.19, respectively, relative to the brightest star.

Because TIC 302965293 is 9″ from the brightest star TIC 302965294, we carried out a photo-center analysis using the eclipses observed in the source (see, e.g., Kostov et al. 2024). The top panel of Fig. 20 shows the layout of the neighboring stars at the TESS pixel level. The bottom panel indicates the relative location of TIC 302965293 (G = 15.13), while the red points are the photocenter locations of the primary EB eclipses in comparison. The grey contours are the confidence levels for the eclipse locations. Thus we see that TIC 302965293 is robustly ruled out as the source of the eclipses. This photo-center analysis is not quite sufficient to distinguish TIC 302965294 from TIC 650024463 because they are only 2″ apart. However, TIC 302965294 has already been ruled out via the shallow depth of the eclipses (see above).

Thus, even though we now know that TIC 650024463 is the correct host of the triple system, due its extreme proximity to the brighter star TIC 302965294, there are no available cataloged passband magnitudes for this target – with the exception of the three Gaia DR3 magnitudes. Hence, in this case, we departed from our usual SED + PARSEC fitting procedure and simply fixed the Gaia DR3 parallax (and, hence, distance). We also set the interstellar reddening E(B − V) parameter to the value given in TIC v8.2, and kept it fixed.

Apart from the light contamination issue, the complex photodynamical analysis resulted in quite robust results for the system parameters. In particular, the orbital dynamics and 3D configuration of the entire stellar system were well determined and, moreover, the relative (dimensionless) stellar quantities (i.e., the mass, temperature and radius ratios of the stars) are also determined with high accuracies. This success is due to the fact that TESS observed the target during nine sectors (the most amongst the currently investigated seven systems). These observations recorded not only six third-body eclipses (Fig. 7), but the derived ETV curves (see in Fig. 15) show marked, large amplitude, dynamically dominated P_out-period cycles, and cover more than half of a full apsidal revolution cycle (the latter of which is also clearly dominated by third-body effects).

These ETV properties are consequences of the fact that this system is by far the tightest triple in our sample with $P_{\mathrm{in}}=7\stackrel{\text{d}}{.}20$, $P_{\mathrm{out}}=108\stackrel{\text{d}}{.}73$ and, hence, P_out/P_in ≈ 15.1. From a dynamical perspective, this triple star system is very similar to the three triples (KIC 9714358 – $P_{\mathrm{in}}=6\stackrel{\text{d}}{.}47$, $P_{\mathrm{out}}=104\stackrel{\text{d}}{.}08$; KIC 5771589 – $10\stackrel{\text{d}}{.}68$, $113\stackrel{\text{d}}{.}87$; TIC 219885468 – $7\stackrel{\text{d}}{.}51$, $111\stackrel{\text{d}}{.}51$) which were recently analyzed in detail by Borkovits & Mitnyan (2023). Besides the very similar inner and outer periods, further similarities are that all of the four systems are (i) quite flat (well within $i_{\mathrm{mut}}=0\stackrel{°}{.}5$), and have (ii) small, but significantly non-zero inner eccentricities (0.0036 ≲ e_in ≲ 0.042), and (iii) moderate outer (0.16 ≲ e_out ≲ 0.39) eccentricities.

As it was pointed out by Borkovits & Mitnyan (2023), such tight, eccentric and flat systems are ideal for observational detections of higher order third-body perturbations on timescales as short as a few years. The different apsidal motion amplitudes of the ETV curves in the first and the second half of the TESS observations of TIC 650024463 (Fig. 15) suggest the additional presence of observable effects of higher, octupole order third-body perturbations in the current system.

To further study this effect, we numerically integrated the motion of the triple system for the forthcoming few centuries. From these integrations we show and discuss some representative results for the current century. In the upper panel of Fig. 21 we plot the numerically generated ETV curve for the entire 21st century. The $P_{\mathrm{apse}}^{\mathrm{obs}}=10.9$ yr-period, dynamically forced apsidal motion is readily visible from the anticorrelated nature of the primary (red) and secondary (blue) ETV curves. These variations correspond nicely to the period of revolution of the inner orbital ellipse in the observational frame of reference, that is, to the variation of the observational argument of periastron (ω_in) plotted with red dots in the lower left panel of Fig. 21. This period of 10.9 yr, however, is substantially shorter than the theoretical value of $P_{\mathrm{apse}}^{\mathrm{theo}}=15.45\pm 0.08$ yr (see in Table 9) which was calculated from the simple, quadrupole-order theory. This discrepancy indicates that the quadrupole-order perturbation theory may well be insufficient in the present situation.

Fig. 21. Numerically integrated models for TIC 650024463 for the entire 21st century. Upper panel: model ETV curve: red and blue lines represent the primary and secondary curves, respectively, while red circles and blue squares stand for the TESS-observed ETV points. Besides the rapidly oscillating vertical ‘spikes’, which come from the Pout period of 108 d, namely, the so-called medium term third-body perturbations, the 10.9 year third-body forced apsidal motion is also clearly visible. The amplitudes of consecutive half-cycles vary substantially which is the consequence of the octupole-order eccentricity perturbations. Lower left panel: variations of the observable and dynamical arguments of periastron (i.e., the apsidal motion of the triple) until year 2100. Lower right panel: cyclic variations of the inner eccentricity (ein – red) and the corresponding, similar-period variations of the differences of the dynamical inner and outer arguments of periastrons (${\omega }_{\mathrm{in}}^{\mathrm{dyn}}-{\omega }_{\mathrm{out}}^{\mathrm{dyn}}$ – blue) during the present century. See text for further details.

Besides this weak, indirect indication for the presence of the higher-order perturbation effects, however, the very same ETV curve carries more clear and well observable evidence about the presence of higher order effects. As one can see in the upper panel of Fig. 21, the amplitudes of the consecutive half-apsidal motion cycles in the ETV curve are really highly variable during the entire 21st century. These amplitude variations do not correspond strictly to the apsidal motion cycles, (and hence, the amplitude varies cycle-by-cycle), but instead, correspond to the ∼12.9 yr-period (inner) eccentricity cycles. The lower right panel of Fig. 21 illustrates that the inner eccentricity (again, apart from the P_out-period mid-timescale variations) oscillates between (e_in)_min ∼ 0.015 and (e_in)_max ∼ 0.028 with the period mentioned above. As is known from the theory of the secular hierarchical three-body problem, for such an almost perfectly flat ($i_{\mathrm{m}}=0\stackrel{°}{.}41\pm 0\stackrel{°}{.}27$) triple, the lowest order (quadrupole) perturbations in the inner eccentricity vanish (see, e.g., Naoz 2016) and, consequently, these e-cycles cannot be explained with the usual, quadrupole-order perturbations. On the other hand, these oscillations are in phase with the variations of the differences of the inner and outer dynamical arguments of periastron ($\mathrm{\Delta }\equiv {\omega }_{\mathrm{in}}^{\mathrm{dyn}}-{\omega }_{\mathrm{out}}^{\mathrm{dyn}}$) and reach their extrema (i.e., minima and maxima of e_in’s) when Δ = 0° or 180° and, hence, the octupole-order dominant perturbation term in e_in, which is proportional to sin Δ (see, e.g., Borkovits & Mitnyan 2023, Eq. (12)), disappears. Thus, one can readily conclude that these e-cycles, and the concomitant variations in the shape of the apsidal ETV curves are direct, observable manifestations of the higher, octupole order third-body perturbation effects.

8. Summary, discussion, and conclusions

In this work we have presented seven new compact triply eclipsing triples with a full solution to their stellar and orbital parameters. All of these systems were found by searching through millions of TESS photometric light curves looking for third body eclipses in binary systems (see e.g., Kristiansen et al. 2022; Rappaport et al. 2022).

We utilized the TESS photometric light curves, the ETV points derived from the light curves, archival SED data, archival photometry from ASAS-SN and ATLAS, and in some cases, ground-based follow up RV observations as well as eclipse photometry. These were combined in a complex photodynamical analysis where we solve for all the system parameters, as well as the distance to the source. Typical uncertainties on the masses and radii are in the range of a couple of per cent to not much more than ∼5%. Uncertainties on the angles associated with the orbital planes (e.g., i_out and i_mut) range from a fraction of a degree to about a degree. See Tables 6 through 9.

In Fig. 22 we present a set of four correlation plots among some of the physically interesting parameters associated with our collection of 33 recently analyzed triply eclipsing triples. To the seven sources presented in this work, we have added the 15 triples from our previous closely related papers (Rappaport et al. 2022, 2023), as well as eleven triples studied in Borkovits et al. (2019a, 2020b, 2022a,b), Czavalinga et al. (2023b), and Mitnyan et al. (2020) using largely the same selection criteria and methods of analysis.

Fig. 22. Statistical plots for properties of 33 triply eclipsing triples uniformly analyzed (see text for references). Top-row panels: eout vs. Pout and iout vs. imut. In the upper right panel the vertical lines denote the transition to the ZLK cycles, and to retrograde orbits, respectively. Bottom-row panels: qin vs. qout and RB vs. MB. The sloped dashed lines in the latter plot are for RB = 1 MB and = 5 MB (both expressed in solar units), as rough guides of unevolved and quite evolved stars, respectively.

The top left panel in Fig. 22 shows how e_out varies with P_out. In general, there is little correlation, except for the fact that the most circular outer orbits tend to have the relatively shorter outer periods, where the circularization may be brought about by tidal damping with an evolved tertiary. See also Fig. 19 and the associated discussion.

The top right panel of Fig. 22 shows the inclination angle of the outer orbit vs. the mutual inclination angle of the inner binary with respect to the outer orbit. The fact that most of the values of i_out lie near 90° is a selection effect since these triples were, in fact, discovered from the presence of third-body eclipses. To a lesser extent, the same selection also holds for the small values of i_mut, otherwise third body eclipses would be somewhat more difficult to detect. Two of the systems have large enough i_mut (20° and 40°) to undergo substantial precession of their orbital planes. Finally, there is one triple (TIC 276162169) in a nearly flat system, but where the outer orbit is retrograde with respect to the inner EB. These are rare systems.

The q_in vs. q_out plot (lower left panel in Fig. 22) shows an interesting feature, namely that the ratio of q_out ≡ M_B/M_EB never exceeds unity (and not just by definition) and has a median value of 0.62. This result is not obviously caused by any selection effects. For purposes of this plot only we define q_in as the ratio of the lower mass to the higher mass star in the EB.

This result for ${\overline{q}}_{\mathrm{out}}$ having a mean of ∼0.62 and never exceeding unity, favors a formation scenario for tight triples involving sequential disk fragmentation and subsequent accretion proposed by Tokovinin & Moe (2020; see also Tokovinin 2021; Offner et al. 2023, for recent reviews). In this scenario, the accretion rate onto the protostellar cores has an inverse relation with their mass-ratio. Hence, the less massive core of the progenitor EB secondary accretes faster while it has lower mass. Similarly, the tertiary, with the youngest core, will accrete faster while it has lower mass than the sum of the elder, inner cores. Consequently, under the assumption that there is enough matter, the triple tends toward a ‘double twin’ scenario, which means that both q_in and q_out would be about unity. Naturally, the accretion may stop earlier, resulting in q_out < 1, but it does not allow for q_out > 1.

The bottom right panel of Fig. 22 shows the relation between the tertiary radius and its mass. About a dozen of the stars lie fairly near the main sequence, roughly represented by the R = M line in the figure. However, the other ∼20 of the tertiary stars have evolved distinctly away from the main sequence. That is due both to an observational selection effect involving the detectability of third-body eclipses, and the fact that the EB stars cannot grow as large as the tertiaries without overflowing their Roche lobes.

A couple of the statistical plots rely on the following assumption made in our photodynamical analyses. Here we have generally adopted the assumption that the inner EB stars are coeval with the tertiary stars and, specifically, that no mass has been transferred among the stars or lost from the system. In turn, this includes the assumption that the tertiary star itself was always a single star, and not the merger product of another binary star.

Another item worthy of a brief discussion is the photometric distance we derive vs. the parallactic distance found by Gaia. Our distances are plotted in Fig. 23 where they are compared to Gaia’s distances. These are the same 33 triply eclipsing systems discussed above (of which we have photometric distances for 31), and augmented by 9 other triple systems analyzed in the same manner, but which are not triply eclipsing. These latter sources are taken from Borkovits et al. (2020a, 2022b), Gaulme et al. (2022), and Borkovits & Mitnyan (2023).

Fig. 23. Comparison of Gaia distances (Bailer-Jones et al. 2021) to 40 triple systems with distances found from our photodynamical fits to the system parameters. The systems marked in red are 6 from the current work with fitted distances. The blue curve is the line where the Gaia and our photometric distances would match. The formal correlation coefficient between the two sets of data is 0.993, and the fitted slope is 1.028 ± 0.004. There is one system, TIC 52041148, at a Gaia distance of 5931 ± 400 pc and photometric distance of 1357 ± 30, that is off the plot.

The overall general agreement between the Gaia distances and our photometric distances is quite striking, and with generally comparable error bars. The formal correlation coefficient between the two sets of distances is 0.993 and the fitted slope between them is 1.028 ± 0.004. Aside from this general agreement, there are a fair number of points where the two distances differ by more than a few statistical error bars. In these cases, it is not obvious from an inspection of this figure which distance is the more accurate one. This may be a case where one or both sets have underestimated uncertainties. In the case of Gaia, the outer orbit, which can be of the order of a 1/2–1 year, may be slightly distorting the parallax measurement¹⁶. For our photometric distances, there may be issues with inhomogeneous SED data, some of which may be taken during eclipses, the PARSEC isochrones are based on isolated and non-rotating stars, and there are uncertainties in the wavelength-dependent interstellar extinction that we use.

Now that the 33 triply eclipsing triples discussed above (see Fig. 22) have their basic parameters determined, and all but three have outer periods shorter than a year, this would make for an interesting follow-up ground-based eclipse timing project. Most of these systems are accessible with modest-size amateur telescopes since the majority of the 33 objects have G magnitudes of ≲13. The ETV data from TESS itself was typically quite fundamental in the determination of the systems parameters via the photodynamical analyses. Therefore, future timing observations of the regular EB eclipses in these systems could significantly improve the parameter determinations. The ETV delays are typically in the minute range, so readily within the reach of amateur observations. Since there are a large number of these triply eclipsing triples, there should be a system to observe nearly every night.

Acknowledgments

This project has received funding from the HUN-REN Hungarian Research Network. T. B., T. M., I. B. B. and A. P. acknowledge the financial support of the Hungarian National Research, Development and Innovation Office – NKFIH Grants K-147131 and K-138962. This paper makes extensive use of data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. Distances and other astrometric properties for all targets, and outer orbits for two of our sources, were taken from the archives of 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. Some of the SED fluxes and magnitudes were obtained with the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Additionally, some of the SED fluxes and magnitudes were obtained with the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. We used the Simbad service operated by the Centre des Données Stellaires (Strasbourg, France) and the ESO Science Archive Facility services (data obtained under request number 396301). This research has also made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in Ochsenbein et al. (2000). Our searches for outer orbital periods from photometry utilized archival data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile. Our searches for outer orbital periods also made use of archival data from ASAS-SN data for this work. The ASAS-SN Sky Patrol (All-Sky Automated Survey for Supernovae) is an automated program to search for supernovae and other astronomical transients, headed by astronomers from Ohio State University (US). It has 20 robotic telescopes in both the northern and southern hemispheres and surveys the entire sky approximately once every day.

References

Appendix A: Eclipse times of the inner EBs of the seven triples

In this appendix we tabulate the times of the individual primary and secondary eclipses of the inner EBs of the triples considered in this study. These naturally include mostly eclipses from TESS, plus a few that were observed from the ground (Tables A.1–A.7).

Appendix B: Radial velocities of TICs 287756035 and 321978218

In this appendix we list the RV data obtained for two tertiary components out of our seven triples (Tables B.1-B.2). All of these data were obtained with the CHIRON spectrometer on the SMARTS telescope.

All Tables

All Figures

	Fig. 1. Light curves (blue points) and model fits (smooth red curves) for the third-body eclipses of TIC 133771812. The third body eclipses are marked with vertical arrows. The sector numbers are indicated in the lower left corner of each panel. Calculation of the model curves will be discussed in Sect. 6.
In the text

	Fig. 2. Light curve (blue points) and model fit (smooth red curves) for the only third-body eclipse of TIC 176713425 observed with TESS. The third body eclipses are marked with vertical arrows. The other notation is the same as in Fig. 1.
In the text

	Fig. 3. Light curves (blue points) and model fits (smooth red curves) for the third-body eclipses of TIC 185615681 (marked with vertical arrows). The other notation is the same as in Fig. 1.
In the text

	Fig. 4. Light curves (blue points) and model fits (smooth red curves) for the third-body eclipses of TIC 287756035 (marked by horizontal arrows). The other notation is the same as in Fig. 1.
In the text

	Fig. 5. Light curves (blue points) and model fits (smooth red curves) for the TESS-observed third-body eclipses of TIC 321978218 (marked by horizontal arrows). The other notation is the same as in Fig. 1.
In the text

Fig. 6. Light curves (blue points) and two kinds of model fits for TIC 323486857. The light-gray fit represents the pure three-body model light curve, while the red fit is the net fit of the triple star model and a four-frequency Fourier model of the light curve distortions (see text for details). The region of the third-body eclipses is indicated with the horizontal arrows. The other notation is the same as in Fig. 1.

In the text

	Fig. 7. Light curves (blue points) and model fits (smooth red curves) for TIC 650024463. The vertical and horizontal arrows mark the times of third body eclipses. The other notation is the same as in Fig. 1.
In the text

	Fig. 8. Eleven partially observed third-body eclipses of TIC 321978218 from the five seasons of the SWASP project. Gray dots represent the individual observations, while larger blue circles denote their 10-min. averages which were used for the joint photodynamical analysis. Red curves stand for the best-fit model.
In the text

Fig. 9. Archival outer orbit folds. The ASAS-SN and ATLAS archival data were used to find the outer period of these triple systems independently of what we were able to learn from ETV or RV curves, or third-body eclipses observed with TESS. Typically, there are 3000–6000 archival photometric points spanning a decade. In the case of TIC 287756035, there is about 1/4 of the outer orbital phase which is missing due to the fact that the outer period nearly exactly matches that of the Earth year, and the consequential observing seasons. The light red curves are fits to a non-physical function consisting of a modified hyperbolic secant (e.g., Eq. (1) of Rappaport et al. 2016). The fit covers only one of the two cycles for clarity in viewing the eclipses, and is used only to find the phases of the eclipses.

In the text

	Fig. 10. TESS primary and secondary ETV curves (red and blue circles, respectively), with the best-fit photodynamical solution (small red and blue points) for TIC 133771812 (see Sect. 6). The horizontally centered black curve represents the pure LTTE contribution. Vertical lines mark the times of the observed outer eclipses (dashed green-the binary occulting the tertiary star).
In the text

	Fig. 11. ETV curves formed from the TESS and ground-based follow up observations with the best-fit photodynamical solution for TIC 176713425. The horizontally centered black curve represents the pure LTTE contribution. Dashed vertical green line marks the position of the only (two-dipped) third-body eclipse event which was observed with TESS.
In the text

	Fig. 12. TESS ETV curves with the best-fitted photodynamical solution for TIC 185615681 (TZ Pyx). The horizontally centered black curve represents the pure LTTE contribution. Vertical dashed green lines represent the times of the observed outer eclipses when the binary occults the tertiary star, and vice versa for the dot-dashed red lines.
In the text

Fig. 13. ETV curves for TIC 321978218 with the best-fit photodynamical solution superposed. Left panel: ETV curve derived from TESS data only. We note that only those sections of the ETV curves, where TESS data are available, are shown. Vertical black lines denote the different sector boundaries; dashed green lines represent the times of the observed outer eclipses when the binary occults the tertiary star, and vice versa for the dot-dashed red lines. Right panel: overall ETV curve after adding average ‘seasonal’ ETV points (see text for explanation) calculated from SWASP observations. The longer-term ETV curve also nicely shows an additional variation. We model it with an fourth stellar component, making the system a 2+1+1 quadruple. See text for details.

In the text

	Fig. 14. TESS ETV curves with the best-fit photodynamical solution for TIC 323486857. We note that only those sections of the ETV curves, where TESS ETV data are available, are shown. Vertical lines have the same meaning as in Fig. 13.
In the text

	Fig. 15. TESS ETV curves with the best-fitted photodynamical solution for TIC 650024463. The horizontally centered black curve represents the pure LTTE contribution. Vertical lines have the same meaning as in Fig. 13.
In the text

Fig. 16. RV curves and model fits for TIC 185615681. The thin, horizontal line at V = −5.95 km s−1 represents the systemic radial velocity of the triple star system. If TIC 185615681 were a single binary, the RVs of the two components would intersect each other exactly on that line. The offset of consecutive intersections of the RV curves from this line is due to the revolution of the inner pair around the center of mass of the whole triple system. The RV points were taken from Duerbeck & Rucinski (2007).

In the text

	Fig. 17. RV curve and model fit for the third component of TIC 287756035. The thin, horizontal line at V = −3.9 km s−1 represents the systemic radial velocity of the triple star system.
In the text

	Fig. 18. RV curve and model fit for the third component of TIC 321978218. The thin, horizontal line at V = 1.6 km s−1 represents the systemic radial velocity of the quadruple star system.
In the text

Fig. 19. Parameter eout plotted against the fractional radius of the tertiary component rB = RB/aout for 33 triply eclipsing triple systems that we have studied recently (see Sect. 8 for details). The red line sketches out the inverse correlation between these two parameters, with a functional form of eout ≃ exp(−42rB). Exceptions are the two systems marked with a faint green ellipse with very low eccentricities, but not large values of rB. Numbers next to individual sources are the TIC numbers. We omit error bars on the systems with eout ≳ 0.1 for increased clarity.

In the text

Fig. 20. Positional information for TIC 650024463 and its two neighbors. Upper panel: A 4 × 4 pixels TESS image near the position of TIC 650024463 (marked by the red arrow) from Sector 1, highlighting nearby field stars (orange dots). Lower panel: photocenter measurements for the primary eclipses detected in Sector 1 compared to the relative pixel position of TIC 302965293. The small red dots represent the measurements for each eclipse (four in Sector 1); the large red circle represents the average photocenter; and the grey contours represent the corresponding 1-, 2- and 3-σ confidence intervals. As seen from the panel, the measurements clearly rule out TIC 302965293 as a potential source of the detected eclipses.

In the text

Fig. 21. Numerically integrated models for TIC 650024463 for the entire 21st century. Upper panel: model ETV curve: red and blue lines represent the primary and secondary curves, respectively, while red circles and blue squares stand for the TESS-observed ETV points. Besides the rapidly oscillating vertical ‘spikes’, which come from the Pout period of 108 d, namely, the so-called medium term third-body perturbations, the 10.9 year third-body forced apsidal motion is also clearly visible. The amplitudes of consecutive half-cycles vary substantially which is the consequence of the octupole-order eccentricity perturbations. Lower left panel: variations of the observable and dynamical arguments of periastron (i.e., the apsidal motion of the triple) until year 2100. Lower right panel: cyclic variations of the inner eccentricity (ein – red) and the corresponding, similar-period variations of the differences of the dynamical inner and outer arguments of periastrons (${\omega }_{\mathrm{in}}^{\mathrm{dyn}}-{\omega }_{\mathrm{out}}^{\mathrm{dyn}}$ – blue) during the present century. See text for further details.

In the text

Fig. 22. Statistical plots for properties of 33 triply eclipsing triples uniformly analyzed (see text for references). Top-row panels: eout vs. Pout and iout vs. imut. In the upper right panel the vertical lines denote the transition to the ZLK cycles, and to retrograde orbits, respectively. Bottom-row panels: qin vs. qout and RB vs. MB. The sloped dashed lines in the latter plot are for RB = 1 MB and = 5 MB (both expressed in solar units), as rough guides of unevolved and quite evolved stars, respectively.

In the text

Fig. 23. Comparison of Gaia distances (Bailer-Jones et al. 2021) to 40 triple systems with distances found from our photodynamical fits to the system parameters. The systems marked in red are 6 from the current work with fitted distances. The blue curve is the line where the Gaia and our photometric distances would match. The formal correlation coefficient between the two sets of data is 0.993, and the fitted slope is 1.028 ± 0.004. There is one system, TIC 52041148, at a Gaia distance of 5931 ± 400 pc and photometric distance of 1357 ± 30, that is off the plot.

In the text