Eight new, triply eclipsing triple star systems

Introduction

It was only 15 years ago, when the discoveries of the first two triply eclipsing triple star systems were reported [22, 25]. These novel findings were based on investigations of the photometric data from the Kepler-spacecraft [20]. Later, further photometric observations with the Kepler space telescope led to the discovery of more than a dozen additional triply eclipsing systems.

While a few subsequent discoveries of triply eclipsing triples were reported sporadically with the use of other space-borne and ground-based instruments (see the list of Borkovits 2022), the discovery of triply eclipsing triples became routine only with the advent of the TESS space telescope [56].

to its monitoring of the entire sky over the past eight years¹, by now TESS has found more than a hundred triply eclipsing triple star candidates.

Triply eclipsing triple star systems, in general, consist of at least one close² and tight³ triple stellar subsystem, where, moreover, the innermost pair is generally an eclipsing binary (EB). It is not necessary for a triply eclipsing system to be close and/or tight, but the closeness of the third component substantially increases the chance of mutual extra eclipsing events. Similarly it is also not required that the inner pair should be an EB (though, in the absence of regular, inner eclipses, the name ‘triply eclipsing’ would not be exactly correct); however, the discovery of extra eclipsing events, and the investigations of such triples are generally much easier when inner, regular eclipses are also present.

In several former papers, we investigated in a homogeneous manner TESS-discovered triply eclipsing stellar systems using a spectro-photodynamical analysis. Apart from the very first two TESS-discovered triply eclipsing triples (TIC 209409435 - Borkovits et al. 2020b; TIC 278825952 - Mitnyan et al. [42]), we have subsequently tended to group systems into sets of between 3 and 10 triples, and then investigated these subsets in detail (Borkovits et al. [16] (2022a); Rappaport et al. [52, 53, 54]; Czavalinga et al. [24]). In some cases, however, when a particular system was found to be exceptionally interesting and/or a record holder in some way, we dedicated individual papers to such objects (Powell et al. [49]; Kostov et al. [33]; Borkovits et al. [17]; Rappaport et al. [55]). Here we also note the paper of Powell et al. [50] which was dedicated to only a particular triply eclipsing triple stellar system, but where the method of analysis followed a somewhat different approach.

This current paper also joins the former category of the study of a modest-sized group of eight triply eclipsing triple star systems analysed with our complex photodynamical code. For these particular systems, such an analysis has become possible only this year (2026), thanks to the latest TESS observations, as well as to some new additional, ground-based photometric and/or spectroscopic follow-up observations.

In Section 2 we describe the collection of systems that we have selected for this detailed study. We explain briefly how these sources were discovered. In Section 3 we analyse the archival data for these triples. In Section 4, the light curves exhibiting third body eclipses are introduced along with model fits, and briefly discussed. The eclipse timing variation curves (ETVs) are introduced and discussed in general in Section 5. The photodynamical code, with which the system parameters are extracted, is briefly reviewed in Section 6. In Section 7 we summarize the system parameters in a set of comprehensive tables, and the results for each individual system are presented. Finally, some conclusions and statistical assessments are made in Sect. 8.

The eight triply eclipsing triples

Catalogued Information

For this collection of eight formerly unanalyzed triply eclipsing triple systems, we started by gathering their main catalogued properties (see Table A.1). Five of the eight triples are located in the southern ecliptic hemisphere, including one system (TIC 220456364) which is located in the Southern Continuous Viewing Zone (SCVZ), while two systems are located in the northern hemisphere. The eighth triple, TIC 306690194, is also in the northern hemisphere, and is located near to the ecliptic; hence, it was observed in some of the ecliptic sectors of the TESS satellite (see in Table A.2).

Discovery

Before the advent of the eight-year-long monitoring of nearly the entire sky with TESS, the triply eclipsing triple nature of all eight of the currently analysed systems was unknown. Moreover, in the case of the majority of the current systems, the eclipsing binary (EB) nature of the inner pair was either not known or it was discovered only in the last few years due to such ground-based surveys as, e.g. the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. [58]). Regarding the current eight systems, seven of them were identified by our ‘Visual Survey Group’ (VSG; Kristiansen et al. [35]) in the manner described in detail, e.g., in Rappaport et al. [54]. Note, however, that one of these seven (TIC 343763382 = ASASSN-V J181441.19-755530.2) was also independently discovered and reported by Zasche et al. [66], who also gave the correct outer period. Moreover, the triply eclipsing nature of TIC 232874003 was found, and communicated personally by Petr Zasche.

Follow-up photometric observations

Since the majority of the currently investigated eight systems are located far in the southern hemisphere, they are mostly unavailable to European and North-American observers. Therefore, we organized ground-based photometric follow-up campaigns only for the three targets which are located in the northern hemisphere (TICs 306690194, 358795247 and 379084660). We obtained usable ground-based measurements for all of them. The ground-based photometric observations of TIC 3587956247 were carried out with two identical 80-cm telescopes in Gothard Astrophysical Observatory, Szombathely (GAO80) and Baja Astronomical Observatory of Szeged University, Baja (BAO80), both in Hungary. The GAO80 observations were carried out in Sloan $z'$-band, while the BAO80 ones a Sloan $r'$-band filter was used. One regular eclipse observation of TIC 379084660 was also carried out with the BAO80 instrument, using a Sloan $r'$-band filter. The details of both these instruments, and the way in which the observations were made, were described briefly in Borkovits et al. (2022a). These observations were carried out over 13+1 nights (for TICs 358795247 and 379084660, respectively) and resulted in 9+1 further regular eclipse observations as well as, in the case of TIC 358795247, some parts of one additional third-body eclipse event.

Finally, only after completing the spectro-photodynamical analyses, we received 8+16 nights of additional observations of TICs 306690194 and 358795247, respectively. These latter data sets were obtained from the small private observatory of R.U. in Jílové u Prahy (Czech Republic), as well as from some remote telescopes in southern Switzerland, in Carona, operated by the same observer. These observations were carried out with three different telescopes, namely, an unfiltered 34-mm refractor and 150-mm and 200-mm reflectors in Cousins $R_{\mathrm{c}}$-band. When we plotted these photometric light curves against the formerly accepted models, we found that they agreed with the emulated model light curves both in timings and shapes. Therefore, we decided not to rerun our fitting processes, including these new data, but to use them simply to verify the former results. The eclipse times determined from all the ground-based observations (including those that were derived from these late ground-based observations) are listed together with the eclipse times extracted from the TESS-light curve in Tables of H.1-H.8.)

Follow-up spectroscopic observations

Three of the eight targets were bright enough for ground-based spectroscopic follow-up observations with moderate-sized telescopes to obtain radial velocities (RVs). In the case of the brighter northern targets, i.e., TIC 306690194 and TIC 358795247, we obtained high-resolution spectra using the Tillinghast Reflector Echelle Spectrograph (TRES; Fürész [28]; Szentgyorgyi & Fürész [60]) on the 1.5-m Tillinghast reflector at the Fred L. Whipple Observatory (Mount Hopkins, Arizona, USA), which is operated by the Smithsonian Astrophysical Observatory. This instrument has a resolving power of about 44,000, and observes the wavelength range from 390 to 910 nm. For TIC 306690194 we obtained six spectra, from which we were able to extract RVs for the two members of the inner double-lined spectroscopic EB (stars Aa and Ab). The velocities were measured using TODCOR, a two-dimensional cross-correlation technique (Zucker & Mazeh [68]). In the case of TIC 358795247, which the ETV data indicate is a quadruple system (see Sect. 5 below), we were able to identify the spectral lines of all four components, and measure their velocities in most of our 29 spectra. The measurements were made with QUADCOR, which is an extension of TODCOR to four dimensions (Torres et al. [64]). The velocities for these two systems are shown in Appendix E and tabulated in Appendix G.

The third system for which we carried out spectroscopic observations and, hence, obtained RV data, was TIC 343763382. In the case of this far-southern system we obtained six RV points for all the three stellar components (that is, stars Aa, Ab and B). The data were obtained with the CHIRON fiber-fed optical spectrometer (Tokovinin et al. [62]) which is mounted on an 1.5-m telescope at Cerro Tololo Inter-American Observatory (CTIO) in Chile and operated by the SMARTS consortium.⁴ Spectra with a resolution of 28,000 were processed by the standard pipeline (Paredes et al. [47]) and the RVs were measured by cross-correlation with a binary template and fitting three dips. These RV data are also shown and tabulated in Appendices E and G, respectively.

Archival photometric data

As is the usual case in our analysis of triple systems, just after the discovery of any triply eclipsing triple, we gather all the available archival data for the given target and, carry out an analysis of these data. The two main purposes of such an analysis are (i) to strengthen or refute the triple nature of the object, and (ii) to obtain the outer orbital period, if possible. The latter is necessary for a further, more detailed analysis of the system. In some cases, however, at the time of discovery, the available archival data are insufficient to yield a clear detection of the outer eclipses or for a determination of the outer period. In such cases we iterate our analysis every couple of years as more and more archival data have been acquired. Here, we present the results of our current archival data analysis for these eight systems.

The outer period is a crucial parameter to determine in triple-star systems. This can be measured via RV measurements, ETVs of the inner binary, and the timing of the outer eclipses. Some difficulties are that we often do not have good (or any) RV coverage, and both the ETVs and outer eclipses observed with TESS (our main source of photometry for triples) are typically sparsely sampled over the eight years of TESS observations. This can lead to multiple solutions for the outer orbital period.

We therefore make use of archival data which has thousands of samples spread over years to more than a decade. If the third-body eclipses last for $\sim0.3$-$1$ day, and the outer periods are $\sim30$-$300$ days, then the eclipsing intervals occupy of order $1/2\%$ (± a factor of $\sim5$) of the outer period. Therefore, with some 3000 typical photometric samples in an archival data set, we might expect about 15 low points due to outer eclipses. If these can be detected (i.e., separated from low noise points) via a box least squares (BLS; Kovács et al. [34]), then the outer orbital period can be determined independently from the TESS data.

In this work, we have made use of archival photometry from the ASAS-SN Shappee et al. [58]; Kochanek et al. [32], the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. [63]), and the Zwicky Transient Facility (ZTF; Masci et al. [39]; Bellm et al. [2]) to search for the outer periods of the eight triples studied in this work. This procedure was successful for six of the eight sources in our study, and the results are summarized in Fig. C.1.

If the area under the EB light curve dips in flux is substantially greater than under the outer eclipse dips (as is typically the case), then the EB eclipses must be removed in order to have a realistic chance of detecting the outer eclipses in a BLS transform. To this end, we first remove any obvious high outliers or those below what could reasonably be expected from inner or outer eclipses. The photometric measurements from different filter bands within a given archival data set are median normalized to each other. If we combine more than one archival data set, the magnitudes from each set are, in turn, median normalized to each other. The final data preparation step is to force-fit a harmonic series, starting at the base frequency of the EB orbit, and subtract it from the data. We typically use 75 orbital harmonics, and the resultant linear fit then requires inverting a $151 \times151$ matrix, essentially an instantaneous calculation on a modern laptop.

The archival data set, prepared in the manner described above, is then subjected to a BLS transform⁵. If the highest BLS peak (or a harmonic of comparable amplitude) turns out to represent the ultimately correct outer orbital period, then we report this as an independent and robust detection of the outer period. The pointwise folds of the outer orbital periods for the six sources in this study for which this criterion was met are shown in Fig. C.1. In all cases, phase 0 of the folds is taken to be the approximate mid-time of one of the TESS third-body eclipses. Collections of low flux points at common phases are highlighted by a red curve that represents an arbitrary eclipse shape, width, and depth. These are not fits to the data. In some cases, the secondary outer eclipse is not centered between the two primary outer eclipses, and their locations, if there are no obvious low points to guide us, are informed by the outer orbital parameters ($e_{\rm out}$ and $\omega_{\rm out}$) given in Tables B.2 – B.5.

In four of the six panels of Fig. C.1, the deepest outer eclipse is in agreement with the phasing given by the TESS outer eclipse in that it appears at phase 0. For TIC 343763382, both outer eclipses are of about the same average depth, so either one could have been phased to zero. In the case of TIC 306690194, the deepest outer eclipse actually comes at a phase very near 1.5. The smaller red dip at phase zero corresponds to the outer eclipse times detected by TESS (and too shallow to be picked up by the archival data). The times of the deeper outer eclipses were not observed by TESS, and thus not yet detected. However, the model light curve for this source does indeed show that such outer eclipses will eventually be seen with TESS.

Light curve and model fits

The TESS space telescope has observed a total of 39 triply eclipsing events (i.e., outer eclipses) in these eight triple star systems. From among these events, for illustrative purposes, we show 23 of them in Figure 1, D.1 and D.2 (three for each system, except for TIC 306690194 which has only two). We note that each outer eclipsing ‘event’ really should be referred to as ‘event’. This is so, because, as one can see, there are third-body eclipsing events which are multiply dipped. In general, all outer eclipses are unique, even in the same triple star system. This results from the fact that the shapes and the depths, as well as the number of the dips in flux, even in the same triple system, depend very strongly on the relative positions of the three stars. Moreover, we intentionally selected for display third-body eclipsing events from both the superior and inferior conjunctions of the third-body. Third-body eclipsing events at the superior conjunction of the third star (i.e., when the third star is eclipsed by the members of the inner EB either together or separately) or at the inferior conjunction (when the third star eclipses member[s] of the EB) may differ not only due to the varying relative stellar positions and the different surface brightness ratios, but also due to not atypical significant outer eccentricities. The sole exception is TIC 306690194 where TESS observed only two third-body eclipsing events, and both of them occurred at the inferior conjunction of the third star. (In this regard see also our corresponding discussion above, in Sect. 3.)

We show in Figure 1, D.1 and D.2 the best-fit spectro-photodynamical model light curves (as red curves), superimposed on the TESS flux data points (blue dots) with 1800-sec cadence. The spectro-photodynamical models and results will be discussed below in Sects. 6, 7. For now, however, we simply make some general comments about these plots. As one can see, there are some panels where the accepted model eclipse depths depart from the observed ones. Some examples can be seen, e.g. in the first two panels of Figure 1, where some of the third-body (i.e., outer) eclipses have lower depths in the model than in the observed eclipses. But, further examples can also be seen in some panels of Figs. D.1 and D.2. For example, in the case of the last panel of the second line of Fig. D.1 (third-body eclipse event during the Sector 104 superior conjunction of TIC 232874003) not only the third-body eclipses, but even the regular two-body eclipses, have smaller model depths than the observed ones. But, oppositely, in the other two panels of the same row (that is, other third-body eclipses in the very same triple, but during other TESS sectors), the regular eclipses have the correct depths, while the modeled third-body eclipses show a bit greater depths than the observed ones. In our opinion, the main reason for these discrepancies comes from the contaminating fluxes and stray-light, which vary slightly sector-by-sector during the TESS observations, that do not allow a more perfect modeling. (In our experience, some tenth of a degree in the mutual inclination of the inner and outer orbits, and/or in the locations of the relative ascending node – that is, $\Omega_{\mathrm{out}}-\Omega_{\mathrm{in}}$ – may result in such discrepancies. Despite this, we are convinced that these slightly varying contaminating fluxes and stray light affect our solution certainly within its 1-$\sigma$ uncertainty.) Naturally, these statements about the varying contaminating flux do not explain the in-sector discrepancies in the light curves of TIC 358795247 (upper row in Fig. D.2), where the departures between the observed and model light curves may arise from some unmodelled chromospheric activity that causes real, and rapidly varying, irregular light curve distortions.

ETV curves

After the detection of a third-body eclipsing event in any of TESS light curves, in addition to gathering the archival observations, the other immediate task is the determination of the mid-minima times of all the regular eclipses and the production and analysis of the ETV curve. In fortuitous cases, the ETV curve itself can reveal the third-body period. This can happen because a cyclic behaviour with the outer period ($P_{\mathrm{out}}$) must be present in the ETV of any compact and/or tight triple system. The more continuous the data train, the easier it is to uniquely determine the outer period via the ETV curve itself.

The $P_{\mathrm{out}}$-period cycles seen in an ETV curve, in general, have two basic origins. One of them is a purely geometrical effect. This is the well-known Römer-delay, which is also referred to as the light-travel-time effect [57]. The LTTE results from the finite speed of light ($c$), and occurs in every triple system (as long as it is not observed pole-on), where the distance of the inner EB from the observer varies cyclically as the EB revolves (with period $P_{\mathrm{out}}$) around the common centre of mass of the triple system. It is easy to show that the amplitude of this clearly sinusoidal effect is

or,

where, in this latter form, the numerical constant is valid only when the masses are expressed in solar masses, the period is in days, and the amplitude itself is also in days. (Here, and hereafter, $M_{\mathrm{B}}$ indicates the mass of the tertiary component [i.e., star B], while $M_{\mathrm{AB}}$ stands for the total mass of the system [i.e., $M_{\mathrm{AB}} = M_{\mathrm{A}} + M_{\mathrm{B}}$, where $M_{\mathrm{A}} = M_{\mathrm{Aa}} + M_{\mathrm{Ab}}$], and the orbital elements have their usual meanings. See also these same notations in Table B.1.)

The other effect, which has the same period of $P_{\mathrm{out}}$, arises from perturbations by the third star that result in departures from a pure Keplerian two-body orbital motion of the inner binary. Hence, we call these ‘dynamical effects’ (DEs). Such perturbations operate on different timescales, but one of them, the so-called ‘medium-timescale perturbations’ has a characteristic timescales of $P_{\mathrm{out}}$. Currently, the most-detailed analytic investigations of such perturbations for the ETV curves can be found in Borkovits et al. (2015). When the inner and outer orbits are nearly coplanar, then the characteristic, theoretical amplitude of this approximately periodic effect is

(Rappaport et al. 2013).

In most cases, however, a clear detection of the outer orbital period from the ETV results requires several observing sessions with TESS of a given target. This is so because, in general, one TESS observing session is shorter (in many cases, much shorter) than the outer period itself. In this circumstance, the only things we can detect from a single or few sectors of TESS-observed ETVs, are that either the curves exhibit some curvature, or ETV sections which were taken a few years apart have different slopes. Both of these phenomena are evidence of rapid period variations. Moreover, this is the chief reason why we need photometric ground-based follow-up monitoring of the regular, i.e., inner EB eclipses.

In Figure 2 we illustrate both LTTE and DE dominated ETVs for the current eight systems. The larger red circles and blue boxes represent those mid-eclipse times that were calculated from TESS as well as ground-based primary and secondary eclipse observations. Moreover, the connected, smaller red and blue dots are numerically generated mid-eclipse times, and represent the best spectro-photodynamical model (see Sect. 6). While these red and blue curves come mainly from the sum of the LTTE and DE contributions to the ETV, we also show separately the pure LTTE contribution as the continuous black curve. As one can see, in most cases, the total amplitude of the ETV is much higher than the amplitude of the pure LTTE-term (in black). This means that most of the investigated ETV curves are mainly dominated by dynamical terms. Moreover, the fact that the average of the LTTE contribution, with the sole exception of the case of TIC 358795247, is horizontal means that the average distance of the considered triple system is constant, or varies only linearly (at least on the timescale of the observations). Opposite to this scenario, in the case of TIC 358795247 (i.e. right panel, third row) one can see that the pure LTTE contribution (in black) is practically equal to the total ETV. That is, here the DEs are negligible, and what is seen is actually the pure classical LTTE. Moreover, we can see a longer-term non-linearity in the ETV, as well. Before considering the spectroscopic RV data, all that we could conclude is that some additional, longer term periodicity is present in this system. Later, as was mentioned already in Sect. 2 and will also be discussed below in Sect. 7.2.6, we spectroscopically detected the lines of a fourth star, whose RVs vary linearly and, in accord with the ETV. Therefore, we strongly suspect that the longer term variation is the consequence of a longer period LTTE. Consequently, we model TIC 358795247 as a 2+1+1-type quadruple system.

We also show, as green and brown vertical lines, the times of those third-body eclipsing events that were observed by TESS (or, in the case of the last dashed horizontal line in TIC 358795247, by the BAO80 instrument). One can see that the green lines in every case are located close to the lower extrema of the black LTTE curves, while the brown ones are near to the upper extrema of the LTTE. This is so, because third-body eclipsing events may occur either during the superior or the inferior conjunctions of the third star, that is, when the distance of the EB from the observer is close to its minimum or maximum. (Exact equality would hold only for a circular outer orbit.)

Regarding the ETV curves and the spectro-photodynamical solutions for the other systems (see the different panels of Figure 2), one can readily see that none of the other systems is clearly LTTE dominated. In the cases of TICs 232874003 and 343763382 the LTTE contribution are as large as the DE terms, but the remaining five systems are clearly dynamically dominated in the sense that the LTTE terms are almost negligible relative to the DE contributions.

Finally, we note that the ETV curves obtained from the primary (red) and secondary (blue) eclipses do not exactly overlap each other, for example, in the case of TIC 77177148 (first row, left panel of Figure 2, or in some other ETVs shown in the different panels. When the two ETV curves do not overlap, this directly indicates that the difference between the mid-eclipse times of the primary and the consecutive secondary eclipses differ from half of the current eclipsing period. Since the vertical separation between a primary and the consecutive secondary ETV points for small inner eccentricities is primarily related to $\Delta\propto e_{\mathrm{in}}\cos\omega_{\mathrm{in}}$, one can readily see that this effect is a direct consequence of the non-zero eccentricities of the inner EB orbits. Moreover, when this vertical separation varies, or even reverses sign (e.g., the blue curve can be either above or below the red curve), this is a direct consequence of apsidal motion (AM; i.e., the variation of $\omega_{\mathrm{in}}$)⁶ which, in the current systems, are mostly also dynamically dominated.

The mid-eclipse times used for the derivation of these ETV curves are tabulated in Appendix H. These mid-eclipse times were calculated with the use of template eclipses in the very same manner as was described already in several former papers, for example in Borkovits et al. [6].

Photodynamical models

We carried out detailed (spectro-)photodynamical analysis of the current eight triply eclipsing triples. This analysis was done with the software package Lightcurvefactory that we have used previously for numerous other multi-stellar systems, as were listed in the Introduction. Therefore, our current analyses are practically the same as we have done with more than 40 previous triple-star systems that we analyzed. This tends to result in a homogeneously analyzed set of parameters for all these systems. The details of the entire analysis process, including the input data sets, the input parameters (which might be either adjusted, fixed or constrained), as well as the different adjusting strategies have all be extensively explained in several former papers (like, e.g. Borkovits et al. [7, 12] (2019a, 2020a); Mitnyan et al. [42]) and, hence, we did not repeat them in detail.

Here we note only the basics of Lightcurvefactory. In general, the code is able to handle the motions of the stars using a seventh order Runge-Kutte-Nyström integrator⁷, as well as the coupled spin rotations of all the bodies. It also emulates the light curves for rotationally and tidally distorted stars with any kind of eclipses that may occur in the following hierarchical stellar (or, of course, planetary) systems: binaries, triples, both kinds of hierarchical quadruples (i.e., both 2+2 and (2+1)+1 type hierarchies) and, moreover, some kinds of higher-order hierarchies up to six stars. The code is able to simultaneously analyze photometric light curves (with a maximum of five different passbands), RV data, ETV curves, as well as composite stellar energy distribution (SED) data (i.e., observed passband magnitudes). Moreover, the code contains different sets of pre-built PARSEC isochrone tables (Bressan et al. [21]), from which the physical radii, effective temperatures, as well as different internal passband magnitudes of any given stars can be tri-linearly iterated and, hence, constrained. The latter assumes that a given star has evolved as a single object and that there has been no previous mass transfer. Finally, to solve the inverse problem, the software package contains a Markov Chain Monte Carlo (MCMC)-based search routine, which is used not only to determine the best-fit system parameters, but also to obtain the statistical uncertainties. This feature uses our own implementation of the generic Metropolis-Hastings algorithm (see, e.g., Ford [27]).

Regarding the eight current triple systems, we numerically integrated the motions of the stars for all of them. In seven cases we used a 2+1-type hierarchical triple star configuration. The exception is TIC 358795247, where both the ETV curve and the RV data prefer the presence of a fourth, more distant, component, for which we applied a (2+1)+1-type hierarchical quadruple configuration. We used the TESS light curves for all the systems (see below). In the case of TIC 358795247, we also incorporated two additional light curves from ground-based follow up photometry. The GAO80 observations were carried out in Sloan $z'$-band, while the BAO80 photometry was taken in Sloan $r'$-band. For TIC 379084660, due to the very different eclipse depths (see below), we also included a second light curve in passband Sloan $r'$, despite the fact that this latter light curve contained only one ground-based regular eclipse observation. For all eight systems we used both the primary and the secondary ETV curves in the simultaneous, complex analysis.

Regarding fitting the ETVs as a quasi-independent data set, one may object that the timings of the regular eclipse minima are already inherent in the light curves. However, we have found empirically that a fit to the system parameters is far more sensitive when the ETVs are fit independently from, but simultaneously with, the light curve itself. Each eclipse yields its own precise mid-minimum time by fitting effectively for only a single parameter – namely, the time. When mixed in with the 30 or so other parameters to be adjusted or constrained in the fit to the light curves, the time of minimum becomes unnecessarily entangled with the other parameters and the results are not as precise or accurate.

We had RV data for only three of the systems and, therefore, the photodynamical analysis for just these three targets included RV information. These are TICs 306690194, 343763382 and 358795247, where RV curves for two, three and four stars are available, respectively. All these multiple RV curves were included in our analysis. Here we mention that when both stars’ RV curves are known in a binary, that determines the two stellar masses (if the inclination angle is known). Similarly, in the case of the latter two systems, the presence of the RVs of all the stars strongly constrains all the stellar masses. In the case of TIC 306690194, in theory, only the minimum mass of the tertiary is strictly constrained. For the other five targets, however, where no RV curves are available, we had to use some other methods to estimate the stellar masses and, hence, the physical dimensions of the stars.

Fortunately, we have found in our previous studies that the composite SED together with the PARSEC isochrones can be used very effectively as proxies for the masses. Our experience and results in this latter regard are described in detail in Borkovits et al. (2020a) and Borkovits et al. (2022a). Note that, for the sake of a fully homogeneous investigation, this latter method, i.e., the use of built-in tabulated PARSEC stellar isochrones, as well as composite system SED data, to constrain the fundamental stellar parameters, was also used in the case of the three systems for which RVs were available. These three targets were thereby used as ‘control’ systems, since the two methods (using RVs vs. isochrones for constraining masses and other stellar parameters) can result in different stellar masses. The way this works is that if isochrone-deduced masses turn out to be inconsistent with the RV-inferred masses, this will result in unusually high values of $\chi^2$. We do not find this to be the case. Further discussion in this regard can be found in Appendix F.

In conclusion, in most cases, where no RV data were available, we adjusted and constrained 16 and 15 variables, as follows. The adjusted variables include the primary’s (i.e., star Aa) mass and the other two mass ratios; three and six orbital elements or other parameters of the inner and outer orbits, respectively; extra flux in the TESS-band; and age, metallicity and interstellar extinction of the system. The constrained parameters are the stellar radii and effective temperatures from the built-in PARSEC evolutionary tracks of the three stars, as well as the period and inferior conjunction time of the inner orbit, the (photometric) distance of the triple star, and the passband-dependent limb-darkening coefficients of the stars. Where RV data were available, the systematic RV of the entire system ($V_\gamma$) is an additional constrained parameter. We suggest that the interested reader consult with the two schematic flow-charts given as fig:5 of Borkovits et al. (2020a) and Figure 3 of Borkovits et al. [16] for further information about this whole process.

Turning now to the observational material that we used our analysis, we summarize these as follows:

1. We binned all the available TESS light curves to 1800-sec cadence times. Similar to our previous triples’ studies, we did this mainly because the cadence times of the exclusively available full-frame-image (FFI)-type TESS photo-

metric data were different for the various sectors. In Years 1 and 2 (Sectors 1–26) the cadence time of the FFIs was 1800-sec. Then, in Years 3 and 4 (Sectors 27–55) the FFIs had a cadence of 600-sec, while from Year 5, i.e., after Sector 56, 200-sec cadence FFI data are available. By binning all the data to 1800 s, we effectively obtained a homogeneous cadence for all the data used in our light curve analysis. (Though, naturally, we always determined the mid-eclipse times from the original TESS dataset and, hence, in most cases, we used finer cadenced data.) For six of the eight systems we made use of the FITSH package of Pál [48] for the photometry of the FFI data. This pipeline, however, was only available until the Sector 77 data. Photometry of the newer data was carried out with the freely available Lightkurve [36] software. Moreover, in the cases of TICs 220456234 and 379084660 all the FFI photometry was carried out with this latter software.

(ii) In order to save computational time and also to give higher weight to the most valuable sections of the light curves, we used only data from the vicinity of those light curve sections where eclipses were observed. Therefore, we excluded the out-of-eclipse light curve sections, and retained only those parts where the photometric phase was within the ±0.10-phase domains of the inner primary and secondary eclipses. Moreover, we kept 5-10 day-long light curve sections around each of the third-body eclipsing events (including those light-curve sections near to the (outer) inferior and/or superior conjunctions where no eclipses were observed due to the lower and/or varying outer inclinations)⁸

(iii) In the cases of TIC 358795247 (TIC 379084660), as was mentioned above, we also used two (one) additional light curves that were made from the ground-based follow-up observations. Since the original exposure times were different from night-to-night, and even for the same instrument under different weather conditions, we binned these data to 900-sec, and used these as two further datasets. Also, as mentioned earlier, the GAO80 ground-based photometry was carried out in the Sloan $z'$-band, while the BAO80 data were in the Sloan $r'$-band. These ground-based observations were strongly down-weighted, and were mostly used as a check for times in between the sparse TESS-observed data sections. The latter was accomplished by penalizing (in a $\chi^2$ sense) only really bad fits to the ground-based data.

(iv) As mentioned earlier, especially in Sect. 5, we computed mid-eclipse times from the original TESS and ground-based datasets. These times of minima, again, are listed in Appendix § 17. The use and necessity of ETVs in our spectro-photodynamical analysis were explained above.

(v) For three of the sources, we also included RV data in our complex analysis. In these cases the systemic velocity $V_\gamma$ was also constrained with an a-posteriori minimisation of the $\chi^2_{\mathrm{RV}}$ value at the end of each trial step, and the result obtained for $V_\gamma$ is also tabulated.

(vi) Finally, in all cases, a simultaneous SED fit was also included. The observed catalog magnitudes used for this fitting are listed in Table A.1.

System parameters

Tables of fitted parameters

We tabulate the results of the above described (spectro-) photodynamical runs in the tables of Appendix B. The meanings of most of the parameters are explained in Table B.1, while the results on the eight investigated triple systems are tabulated in the following Tables B.2–B.5. We emphasize again that almost all of the orbital parameters, except for the first two lines, give the instantaneous osculating orbital elements (or are computed from them). These are valid strictly only at the epoch $t_0$ (given in the first line of each table). The osculating parameters are computed at any given time from the instantaneous (Jacobian) coordinates and velocities of the given object⁹, and give the elements of that two-body Keplerian orbit on which the object would move in the future—if all the perturbations were to disappear at that moment. Naturally, the true motion of each object lies on the envelope curve of the instantaneous osculating orbits. Therefore, the instantaneous osculating elements are not identical with the sets of orbital elements that appear, e.g., in the theoretical perturbation equations. Importantly, they cannot be observed directly. For example, the period, $P$, tabulated in the third line of each results table is the osculating anomalous period that belongs to the inner and the outer orbits at instant $t_0$. This period should not be confused with the period that can be inferred, e.g., from an ETV or, RV curve, or simply, by observing the time between the consecutive eclipses.

In order to offer a natural comparison with other observational studies, and also for purpose of predicting future eclipses in these systems, we give the average observing periods, denoted by $P_{\mathrm{ETV}}^{\mathrm{avg}}$ (second line of each table). We calculate these quantities in the following manner. First, we emulate a model light curve of the system using the parameters for the best-fit (spectro-)photodynamical model corresponding to the time interval 2010-01-01 to 2030-01-01. Then, we calculate model mid-eclipse times from these emulated light curves for the 20-year-long interval in a manner similar to what we used to compute mid-eclipse times from the observational data. Naturally, Lightcurvefactory directly generates mid-eclipse times by itself; but instead of those, we use the indirect step of determining mid-eclipse times from the emulated light curves for a double-check on the results. Finally, the dense ETV curves, calculated in such a manner for this 20-year-long interval, are analysed analytically with the same code used and described, e.g., in Borkovits et al. (2015, 2016, 2025a). In this way, we obtain the inner and outer observable average periods together with their uncertainties, and we tabulate these values in the line for $P_{\mathrm{ETV}}^{\mathrm{avg}}$. We also note that these periods have clear geometrical meanings. The value obtained for the inner orbit gives the period that should be used for an ETV curve with zero average slope (i.e., horizontal). The value valid for the outer orbit, however, is the average period of the cycles caused by the third-body on the ETV curve. We stress that this is not necessarily equal to the time interval between two ‘primary’ (or ‘secondary’) outer eclipses, as the time between consecutive, similar types of outer eclipses also depends on the locations of the three stars in their orbits, and can vary by even more than one inner orbital cycle.

Turning back the instantaneous osculating orbital elements, and also to the other quantities given in Tables B.2–B.5, we list there the median of all the posteriors obtained in all the individual runs of the MCMC chains. Moreover, we also give the $1-\sigma$ upper and lower statistical uncertainties. We do not necessarily give the adjusted (or constrained) quantities directly but, in some cases, we list only parameters that are computed from the adjusted (or constrained) parameters. For example, though in each MCMC run the adjusted parameters are $e\cos\omega$ and $e\sin\omega$ for each orbit, in the Tables we directly list $e$ and $\omega$. We also note that several derived quantities are listed as well. For example, the code directly calculates the RV amplitudes ($K_{\mathrm{pri/sec}}$) from other model parameters. Therefore, we stress that the given $K$-s are always derived and not fitted, independent of whether RV data are available or not for a given system.

Besides the tabulated results, in what follows we discuss our results for all eight individual triply eclipsing triple systems. First, however, we make some general remarks. In Borkovits et al. [16] we considered ultra-compact hierarchical triples (UCHTs), that is, such systems where the outer period is less than (or equal to) 100 days. In the current work, we did not apply such a restriction, and only two of the eight investigated systems can be considered UCHTs. (Or, if we allow borderline cases, we may add a third system, TIC 77177148 as a marginal UCHT, given that its outer period is $P_{\mathrm{out}}=101\mathord{.}4$ days.) For the current systems, however, all but one of them have another interesting characterization for triples, i.e., they are ‘tight’, i.e., $P_{\mathrm{out}}/P_{\mathrm{in}}$ is relatively small (see the last line of Table A.1 for the tightness parameter). Apart from TIC 358795247 (which is, however, a quite rare, four-lined spectroscopic quadruple), the tightness parameter of the systems is below 31. Interestingly, in the case of TIC 385933895 (which is one of the two or three UCHTs in the current study) this ratio is smaller than 10, making this triple one of the tightest systems known. Therefore, it is far from surprising that the currently investigated systems are dynamically quite active, where the third-body perturbations already play important roles on monthly timescales. These can be seen nicely in the ETV curves, that is, in Figure 2. (We will return to this question below in connection with TIC 220456364, in Subsect. 7.2.2.)

What we can also say globally for the investigated systems is that all eight triples contain only main-sequence stars. Apart from TIC 358795247, all the inner periods have $P_{\mathrm{in}}>3\mathord{.}3$ day. Or, in other words, there are only two inner binaries from the current eight where the fractional radii of the stars exceed $r=R/a>0.1$. And, even the maximum of the fractional radii (which is $r=0.141$ in the case of star Aa of the above mentioned TIC 358795247) is quite small and, therefore, all the stars can be considered as nearly spherical (though we did not apply such an approximation during our analyses). Due to the small $r=R/a$ values, there are no noticeable proximity effects on the light curves—this can nicely be seen in the light curve plots, i.e., in Figure 1 and D.1–D.2). Therefore, any tidal effects can safely be neglected in all eight systems. This is also the reason why $\Delta\omega_{\mathrm{b}}\gg\Delta\omega_{\mathrm{tide}}$ for all but one of the inner binaries (see the appropriate lines of Tables B.2–B.5) or, in other words, the AM in all cases where it is present, is clearly dynamically dominated.

Finally, note that the light curve sections of the (spectro-)photodynamical model fits can be seen in the above mentioned Figure 1 and D.1–D.2); the model ETVs are plotted against the observed points in Figure 2; while the model RVs are shown against the observed data both in the time and phase domains in Figs. E.1–E.3.

Results for the individual systems

TIC 77177148

TESS observed a two-dipped, $\approx20\%$-deep third-body eclipse event back in November, 2018 (in TESS Sector 5, see the left panel of Fig. 1). While this system was one of the very first triply eclipsing triples discovered in the TESS data, unfortunately we were unable to determine the outer period or to find any working photodynamical models. This was mainly due to the fact that the regular eclipses are quite shallow, and do not reach even $\approx2\%$ deep. This is a consequence of the fact that the vast majority ($\approx88\%$) of the system’s total light comes from the third, most massive stellar component (see in Table B.2), which strongly dilutes the depths of the regular eclipses. (Moreover, note that our solution also resulted in a high-contamination flux ratio of $\ell_{\mathrm{x}}=0.46\pm0.02$ that is an additional source of the light dilution. The source of this contamination is likely TIC 77177146 that is only $\approx28^{\prime\prime}$ away from our target, and which is, in the TESS-band, fainter only by $\Delta T\approx0.9$ mag. According to Gaia DR3, and the derived distance of Bailer-Jones et al. 2021, however, the distance between the two stars is about 28 pc or even larger. Therefore, this object is likely only a visual companion, and does not belong to the triple system.)

In the end, we find that this triple is formed by three low-mass, cool stars. The total mass of the inner binary is less than the mass of our Sun, being $M_{\mathrm{Aa}}=0.49\pm0.04\,M_{\odot}$ and $M_{\mathrm{Ab}}=0.46\pm0.03\,M_{\odot}$ respectively, while the tertiary component is more massive with $M_{\mathrm{B}}=0.86\pm0.05\,M_{\odot}$. Considering the more certain mass ratios, these are $q_{\mathrm{in}}=0.94\pm0.01$ and $q_{\mathrm{out}}=0.90\pm0.02$, respectively. We found the inner orbit of the two M dwarfs to be only slightly eccentric with $e_{\mathrm{in}}=0.0029\pm0.0003$, which is reflected most strongly in the dynamically induced AM. The latter AM effect manifests itself in the low-amplitude, anti-correlated variations of the primary (red) and secondary (blue) ETV curves on the $\approx8$-year timescale of the TESS observations, as can be seen in the first row, left panel of Figure 2. On the other hand, the outer orbit has a moderate eccentricity of $e_{\mathrm{out}}=0.244\pm0.001$. However, in spite of this substantial outer eccentricity, the mutual inclination angle was found to be quite small at $2\mathord{.}2\pm0\mathord{.}4$, which, in turn, implies a practically flat triple system¹⁰.

TIC 220456364