11institutetext: HUN–REN–SZTE Stellar Astrophysics Research Group, H-6500
Baja, Szegedi út, Kt. 766, Hungary 22institutetext: Department of Physics,
Kavli Institute for Astrophysics and Space Research, M.I.T., Cambridge, MA
02139, USA
22email<EMAIL_ADDRESS>33institutetext: Baja Astronomical Observatory of
University of Szeged, H-6500 Baja, Szegedi út, Kt. 766, Hungary.
33email<EMAIL_ADDRESS>44institutetext: Konkoly Observatory,
Research Centre for Astronomy and Earth Sciences, H-1121 Budapest, Konkoly
Thege Miklós út 15-17, Hungary 55institutetext: ELTE Eötvös Loránd University,
Gothard Astrophysical Observatory, Szent Imre h. u. 112, 9700 Szombathely,
Hungary 66institutetext: HUN–REN–ELTE Exoplanet Research Group, H-9700
Szombathely, Szent Imre h. u. 112, Hungary 77institutetext: Amateur
Astronomer, Glendale, AZ 85308 88institutetext: Centre For Computation
Astrophysics, NYC, USA 99institutetext: Princeton University, Princeton, USA
1010institutetext: Amateur Astronomer, 12812 SE 69th Place Bellevue, WA 98006,
USA 1111institutetext: Cerro Tololo Inter-American Observatory $|$ NSF’s
NOIRLab, Casilla 603, La Serena, Chile 1212institutetext: NASA Goddard Space
Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA
1313institutetext: SETI Institute, 189 Bernardo Avenue, Suite 200, Mountain
View, CA 94043, USA 1414institutetext: Citizen Scientist, c/o Zooniverse,
Dept, of Physics, University of Oxford, Denys Wilkinson Building, Keble Road,
Oxford, OX1 3RH, UK 1515institutetext: Brorfelde Observatory, Observator
Gyldenkernes Vej 7, DK-4340 Tølløse, Denmark 1616institutetext: National Space
Institute, Technical University of Denmark, Elektrovej 327, DK-2800 Lyngby,
Denmark 1717institutetext: Amateur Astronomer, 7507 52nd Place NE Marysville,
WA 98270, USA 1818institutetext: Astronomical Institute, Slovak Academy of
Sciences, 05960 Tatranská Lomnica, Slovakia 1919institutetext: Astrophysics
Group, Keele University, Staffordshire, ST5 5BG, UK

# Seven new triply eclipsing triple star systems

S. A. Rappaport 1122 T. Borkovits 1133445566 T. Mitnyan 11 R. Gagliano 77 N.
Eisner,, Flatiron Research FellowHenry Norris Russel Fellow8899 T. Jacobs 1010
A. Tokovinin 1111 B. Powell 1212 V. Kostov 12121313 M. Omohundro 1414 M. H.
Kristiansen 15151616 R. Jayaraman 22 I. Terentev 1414 H.M. Schwengeler 1414 D.
LaCourse 1717 Z. Garai 55661818 T. Pribulla 1818 P.F.L. Maxted 1919 I. B. Bíró
1133 I. Csányi 33 A. Pál 44 A. Vanderburg 22

(Received January .., 2024; accepted … , …)

###### Abstract

Aims. We have identified nearly a hundred close triply eclipsing hierarchical
triple star systems from data taken with the space telescope TESS. These
systems are noteworthy in that we can potentially determine their dynamical
and astrophysical parameters with a high precision. In the present paper, we
report the comprehensive study of seven new compact triply eclipsing triple
star systems taken from this larger sample: TICs 133771812, 176713425,
185615681, 287756035, 321978218, 323486857, and 650024463.

Methods. Most of the data for this study come from TESS observations, but two
of them have Gaia measurements of their outer orbits, and we obtained
supplemental radial velocity (RV) measurements for three of the systems. The
eclipse timing variation curves extracted from the TESS data, the photometric
light curves, the RV points, and the spectral energy distribution (SED) are
combined in a complex photodynamical analysis to yield the stellar and orbital
parameters of all seven systems.

Results. Four of the systems are quite compact with outer periods in the range
of 41-56 days. All of the systems are substantially flat, with mutual
inclination angles of $\lesssim 2^{\circ}$. Including the systems reported in
this work, we have now studied in considerable detail some 30 triply eclipsing
triples with TESS, and are accumulating a meaningful census of these systems.

###### Key Words.:

Stars:) binaries (including multiple): close – (Stars:) binaries: eclipsing –
(Stars:) binaries: general – Stars: fundamental parameters – Stars:
individual: TICs 133771812, 176713425, 185615681, 287756035, 321978218,
323486857, 650024463

##  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 $\lesssim 1000$ days, down to the currently shortest
known outer period of 33.0 d ($\lambda$ 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 et al.,
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_{\rm
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, 2022a; 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\lesssim 15$. However, most of the
discovered triply eclipsing triples were actually found among a much smaller
number ($\sim$1 million) of light curves of preselected eclipsing binaries
(EBs). The latter were identified in the TESS data via machine learning (ML)
searches111We use a neural network binary classifier trained to find eclipses
in light curves as an initial screening tool. (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), Mitnyan et al. (2020) and Borkovits et al. (2022a); 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)222https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html, 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
(ATLAS333https://fallingstar-data.com/forcedphot/; Tonry et al. 2018; Smith et
al. 2020); and (5) the online VizieR SED viewer (Ochsenbein et al.
2000)444http://vizier.cds.unistra.fr/vizier/sed/. 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 LABEL:tbl:mags7
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 $\sim$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
LABEL:tbl:sectors. In Table LABEL:tbl:noticed we point out the two sources
whose outer orbits were found with Gaia555https://vizier.cds.unistra.fr/viz-
bin/VizieR-3?-source=I/357&-out.max=50&-out.form=HTML 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 archive RV
data were also collected for a third target – TIC 185615681.

Finally, in Table LABEL:tbl:input 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).

Table 1: Main properties of the seven triple systems from different catalogs

Parameter | 133771812 | 176713425 | 185615681 | 287756035 | 321978218 | 323486857 | 650024463
---|---|---|---|---|---|---|---
RA J2000 | 08:05:16.560 | 23:35:35.42 | 08:41:08.26 | 13:14:51.86 | 00:04:12.35 | 11:45:48.28 | 02:41:20.17
Dec J2000 |  $-$41:49:57.97 | 42:22:17.57 |  $-$32:12:03.06 |  $-$46:44:29.21 |  $-$50:55:00.91 |  $-$65:02:48.99 |  $-$82:25:16.73
$T^{a}$ | $13.30\pm 0.01$ | $13.75\pm 0.01$ | $10.52\pm 0.01$ | $12.86\pm 0.01$ | $11.35\pm 0.01$ | $11.61\pm 0.02$ | $14.4\pm 0.6$
Gb | $13.76\pm 0.00$ | $14.25\pm 0.00$ | $10.65\pm 0.00$ | $13.47\pm 0.01$ | $11.73\pm 0.00$ | $12.31\pm 0.00$ | $14.85\pm 0.00$
GBPb | $14.12\pm 0.00$ | $14.62\pm 0.00$ | $10.75\pm 0.00$ | $14.00\pm 0.01$ | $12.00\pm 0.00$ | $12.93\pm 0.01$ | $15.20\pm 0.02$
GRPb | $13.19\pm 0.00$ | $13.68\pm 0.00$ | $10.46\pm 0.00$ | $12.79\pm 0.00$ | $11.28\pm 0.00$ | $11.54\pm 0.01$ | $14.17\pm 0.01$
Ba | $14.62\pm 0.01$ | $15.40\pm 0.46$ | $10.97\pm 0.16$ | $14.77\pm 0.08$ | $12.42\pm 0.06$ | $13.68\pm 0.02$ | …
Vc | $14.00\pm 0.12$ | $14.60\pm 0.11$ | $10.68\pm 0.03$ | $13.75\pm 0.13$ | $11.92\pm 0.05$ | $12.68\pm 0.07$ | …
Jd | $12.49\pm 0.03$ | $13.06\pm 0.02$ | $10.56\pm 0.03$ | $11.92\pm 0.02$ | $10.77\pm 0.03$ | $10.53\pm 0.02$ | …
Hd | $12.23\pm 0.02$ | $12.73\pm 0.03$ | $10.50\pm 0.02$ | $11.35\pm 0.02$ | $10.50\pm 0.02$ | $9.97\pm 0.02$ | …
Kd | $12.17\pm 0.02$ | $12.66\pm 0.03$ | $10.48\pm 0.02$ | $11.28\pm 0.02$ | $10.44\pm 0.03$ | $9.83\pm 0.02$ | …
W1e | $12.09\pm 0.02$ | $12.65\pm 0.02$ | $10.16\pm 0.02$ | $11.22\pm 0.02$ | $10.40\pm 0.02$ | $9.70\pm 0.02$ | …
W2e | $12.09\pm 0.02$ | $12.67\pm 0.02$ | $10.21\pm 0.02$ | $11.27\pm 0.02$ | $10.44\pm 0.02$ | $9.75\pm 0.02$ | …
W3e | $12.08\pm 0.28$ | $12.02\pm 0.22$ | $10.17\pm 0.07$ | $11.29\pm 0.11$ | $10.41\pm 0.08$ | $9.57\pm 0.04$ | …
$T_{\rm eff}$ [K]b | $5342\pm 300$ | $5355\pm 115$ | $8325\pm 950$ | $4973\pm 230$ | $6506\pm 800$ | $4518\pm 150$ | $4362\pm 125$
$T_{\rm eff}$ [K]a | $8357\pm 123$ | $5937\pm 125$ | $8696\pm 155$ | $5090\pm 122$ | $5949\pm 130$ | $7197\pm 130$ | …
Radius [$R_{\odot}$]b | $3.41\pm 0.41$ | $3.48\pm 0.11$ | … | $5.07\pm 0.47$ | $1.97\pm 0.50$ | $11.28\pm 1.2$ | $14.3\pm 0.8$
Radius [$R_{\odot}$]a | $2.82\pm NA$ | $3.08\pm NA$ | $3.11\pm 0.12$ | $5.20\pm NA$ | $2.35\pm 0.13$ | $9.47\pm NA$ | …
Distance [pc]f | $1803\pm 46$ | $2144\pm 86$ | $833\pm 13$ | $2139\pm 83$ | $570\pm 7$ | $1818\pm 38$ | $1126\pm 26$
$E(B-V)$a | $0.59\pm 0.01$ | $0.15\pm 0.01$ | $0.13\pm 0.01$ | $0.12\pm 0.00$ | $0.01\pm 0.01$ | $0.75\pm 0.02$ | $0.09\pm 0.00$
$\mu_{\alpha}$ [mas/yr]b | $-2.56\pm 0.01$ | $-3.82\pm 0.01$ | $-3.51\pm 0.01$ | $-15.07\pm 0.02$ | $24.38\pm 0.01$ | $10.61\pm 0.01$ | $-2.43\pm 0.02$
$\mu_{\delta}$ [mas/yr]b | $+5.63\pm 0.01$ | $-6.43\pm 0.01$ | $-3.58\pm 0.02$ | $4.73\pm 0.01$ | $-0.83\pm 0.02$ | $1.18\pm 0.01$ | $-5.01\pm 0.02$
RUWEb,g | 0.975 | 1.00 | 1.11 | 1.15 | 1.79 | 0.957 | 1.07
astr_ex_noise [mas]b,h | 0 | 0 | 0.13 | 0.066 | 0.23 | 0.043 | 0.065
astr_ex_noise_sigb,h | 0 | 0 | 29 | 3.3 | 117 | 2.1 | 1.2
$P_{\rm binary}$i [d] | 12.370 | 1.8987 | 2.3185 | 2.0814 | 0.5698 | 0.8840 | 7.1971
$P_{\rm triple}$i [d] | 242.3 | 52.55 | 55.84 | 367.6 | 57.51 | 41.36 | 107.95

Notes. General: “NA” and ellipses in this table indicate that the value is not
available. (a) TESS Input Catalog (TIC v8.2) (Paegert et al., 2021). (b) Gaia
EDR3 (Gaia collaboration, 2021); the uncertainty in $T_{\rm eff}$ and $R$
listed here is 1.5 times the geometric mean of the upper and lower error bars
cited in DR2. Magnitude uncertainties listed as 0.00 are $\lesssim 0.005$. (c)
AAVSO Photometric All Sky Survey (APASS) DR9, (Henden et al., 2015),
http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=II/336/apass9. (d) 2MASS
catalog (Skrutskie et al., 2006). (e) WISE point source catalog (Cutri et al.,
2013). (f) Bailer-Jones et al. (2021), geometric distances. (g) The Gaia
renormalized unit weight error (RUWE) is the square root of the normalized
$\chi^{2}$ of the astrometric fit to the along-scan observations. Values in
excess of about unity are sometimes taken to be a sign of stellar
multiplicity. (h) Abbreviations for “astrometric_excess_noise” and
“astrometric_excess_noise_significance” (Lindegren et al. 2021;
https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html);
these are a measure of “the disagreement, expressed as an angle, between the
observations of a source and the best-fitting standard astrometric model.”
Values of astr_ex_noise_sig $\gtrsim 2$ are considered significant. (i) Binary
and outer orbital periods from this work; given here for reference purposes.

Table 2: TESS observation sectors for the triplesa

Object | Sectors Observed | Third Body Events
---|---|---
TIC 133771812 | S7,S8,S34,S35,S61,S62 | S7,S34,S61
TIC 176713425 | S16,S17,S57 | S16
TIC 185615681 | S8,S34,S35,S61,S62 | S8,S34,S35,S61,S62
TIC 287756035 | S11,S37,S38,S64 | S37,S64
TIC 321978218 | S2,S28,S29,S68,S69 | S28,S29,S68,S69
TIC 323486857 | S10,S11,S37,S38,S64,S65 |  S10b,S11,S37,S38b,S64b,S65
TIC 650024463 | S1,S12,S13,S27,S28,S39,S66,S67,S68 | S12,S27,S28,S39,S66,S67

Notes. (a) None of these sources will be observed in further TESS sectors
until the end of Cycle 6 observations; (b) both primary and secondary third-
body eclipses observed in these sectors.

Table 3: Other detections of the triplesa

Object |  This Workb |  Gaia DR3 Orbitsb,c | Other
---|---|---|---
|  ($P_{\rm out}$, $e_{\rm out}$, $\omega_{\rm out}$) |  ($P_{\rm out}$, $e_{\rm out}$, $\omega_{\rm out}$) |
TIC 133771812 |  243.87 d; 0.217; 195∘ | … | …
TIC 176713425 |  52.57 d; 0.412; 95.2∘ | … | …
TIC 185615681 |  55.86 d; 0.098; 309∘ | … | Rowan et al. (2023)
TIC 287756035 |  368 d; 0.235; $53.7^{\circ}$ | … | …
TIC 321978218 |  57.51 d; 0.258; $123\degr$ |  57.52 d; 0.28; 116∘ | …
TIC 323486857 |  41.36 d; 0.007; $241^{\circ}$ |  41.35 d; 0.045; 180∘ | …
TIC 650024463 |  107.88 d; 0.323; 351∘ | … | …

Notes. (a) Three of the seven triply eclipsing systems were spotted in prior
broad surveys, but no quantitative analysis of the system parameters
(especially of the constituent stars) was undertaken. Two of them have outer
orbits reported by Gaia, but no third body eclipses were reported by Gaia. One
of the systems, TIC 185615681, was included in the Rowan et al. (2023) list of
mutlistellar systems tabulated from ASAS-SN data, and studied for the apsidal
motion of its inner binary by Zasche (2012) and Kim et al. (2018), as well as
for RVs by Duerbeck & Rucinski (2007). (b) Where available, we show for each
source, the period, eccentricity, and argument of periastron of the outer
orbit, separated by semicolons. (c) The Gaia orbital solutions are both
spectroscopic; Babusiaux et al. (2023); Gaia Collaboration (2023).

Table 4: Input information for the system analysis

Object | TESS 3rd-Body | TESS EB | Archival Outer |  SED Pointsc |  ETV Curved |  RV Datae |  Gaia Orbitf
---|---|---|---|---|---|---|---
|  Eclipse(s)a |  Light curvea |  Eclipsesb |  |  |  |
TIC 133771812 | $\checkmark$ | $\checkmark$ |  | $\checkmark$ | $\checkmark$ |  |
TIC 176713425 | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ |  |
TIC 185615681 | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ |
TIC 287756035 | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ |  | $\checkmark$ |
TIC 321978218 | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$
TIC 323486857 | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ |  | $\checkmark$
TIC 650024463 | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ |  |

Notes. (a) See Figs. 1–7. (b) Fig. 9. (c) Obtained from
http://vizier.cds.unistra.fr/vizier/sed/. (d) See Figs. 10–15. (e) See Figs.
16–18. (f) Babusiaux et al. (2023); Gaia Collaboration (2023)

##  3 Light curve and model Fits

In Figures 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.

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).

Figure 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.  Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

Figure 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.

##  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
$\lesssim 1000$ points in the south. The ATLAS photometry is superior for the
fainter targets (e.g., $\gtrsim 15$th 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)
algorithm666https://exoplanetarchive.ipac.caltech.edu/cgi-bin/Pgram/nph-pgram
(Kovács et al., 2002) to do a blind search for the third-body eclipses of the
outer orbit covering periods from 20-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 $\sim$0.02 to $\sim$0.003. For an
archival sample of, for example, 5000 points, we can then expect somewhere
between 15-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 $\sim$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.

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 $\sim$8-20 $\sigma$ 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 $\sim$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 $\sim$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 eclipse777We return to this question in Sect.
7.2.4..

Figure 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., Eqn. 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.

##  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_{\rm out}^{2}/P_{\rm in}$ (Söderhjelm, 1975), where $P_{\rm out}$
and $P_{\rm 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.

Figure 10: TESS primary and secondary ETV curves (red and blue circles,
respectively), with the best-fit photodynamical solution 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). Figure 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. Figure 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.

Figure 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.  Figure 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. Figure 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.

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 Figures 10 – 15, while the mid-eclipse
times used for the derivation of these ETV curves are tabulated in seven
separate tables in Appendix A.

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 $\sim 90$%
dynamical delays and $\sim 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
$\sim 90$% dynamical delays and $\sim 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
$\sim$$65$% dynamical delays and $\sim$$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.

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 18 as well).

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 ($\sim$$73\%$) over the dynamical delays
($\sim$$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 19). 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_{\rm out}$) and the total (projected)
mass of the inner EB
($M_{\mathrm{A}}\sin^{3}i_{\mathrm{out}}=M_{\mathrm{Aa}}\sin^{3}i_{\mathrm{out}}+M_{\mathrm{Ab}}\sin^{3}i_{\mathrm{out}}$)
and the individual (projected) mass of the tertiary
($M_{\mathrm{B}}\sin^{3}i_{\mathrm{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.

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
($\lesssim 10^{-3}\,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_{\rm
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\aas@@fstack{p}17$ or $0\aas@@fstack{p}73$.) These nominal seasonal ETV
points extend the length of the available ETV data by almost 5 000 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 $\sim$$90$% dynamical delays and $\sim$$10$% LTTE effects.

Figure 16: RV curves and model fits for TIC 185615681. The thin, horizontal
line at $V=-5.95\,\mathrm{kms}^{-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
different offsets 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).
Figure 17: RV curve and model fit for the third component of TIC 287756035.
The thin, horizontal line at $V=-3.9\,\mathrm{kms}^{-1}$ represents the
systemic radial velocity of the triple star system. Figure 18: RV curve and
model fit for the third component of TIC 321978218. The thin, horizontal line
at $V=1.6\,\mathrm{kms}^{-1}$ represents the systemic radial velocity of the
quadruple star system.

##  6 Photodynamical models

All seven of the triple systems in this work have been subjected to a detailed
photodynamical analysis888To our knowledge, the first photodynamical analysis
of a triply eclipsing triple system was carried out by Carter et al. (2011)
for the system KOI-126. 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, 2019b, 2020a, 2020b, 2021; Mitnyan
et al., 2020). Specifically Lightcurvefactory has been used successfully to
analyse compact and wider triple systems, both with and without outer
eclipses.

Table 5: Definitions of triple system parameters in Tables 6– 9

Parametera | Definition
---|---
$t_{0}$ | Epoch time for osculating elements
$P$ | Orbital period
$a$ | Orbital semimajor axis
$e$ | Orbital eccentricity
$\omega$ | Argument of periastron (of secondary)
$i$ | Orbital inclination angle
$\mathcal{T}_{0}^{\mathrm{inf/sup}}$ |  Time of conjunction of the secondaryb
$\tau$ | Time of periastron passage
$\Omega$ | Longitude of the node relative to
| the node of the inner orbit
$i_{\rm mut}$ |  Mutual inclination anglec
$q$ | Mass ratio (secondary/primary)
$K_{\mathrm{pri}}$ | “K” velocity amplitude of primary
$K_{\mathrm{sec}}$ | “K” velocity amplitude of secondary
$R/a$ | Stellar radius divided by semimajor axis
$T_{\rm eff}/T_{\rm eff,Aa}$ | Temperature relative to EB primary
fractional flux | Stellar contribution in the given band
$M$ | Stellar mass
$R$ | Stellar radius
$T_{\mathrm{eff}}$ | Stellar effective temperature
$L_{\mathrm{bol}}$ | Stellar bolometric luminosity
$M_{\mathrm{bol}}$ | Stellar absolute bolometric magnitude
$M_{V}$ | Stellar absolute visual magnitude
$\log g$ | log surface gravity (cgs units)
$[M/H]$ | log metallicity abundance to H, by mass
$E(B-V)$ | Color excess in B-V bands
extra light, $\ell_{4}$ | Contaminating flux in the given band
$(M_{V})_{\mathrm{tot}}$ | System absolute visual magnitude
distance | Distance to the source

Notes. (a) The units for the parameters are given in Tables 6– 9. (b) The
superscript of “inf/sup” indicates inferior vs. superior conjunctions. (By
default we give inferior conjunctions. Superior conjunctions are indicated by
asteriks.) (c) More explicitly, this is the angle between the orbital planes
of the inner binary and the outer triple orbit.

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
$\times$ mass, radius, $T_{\rm 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 $\ell_{4}$ 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_{\rm 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
LABEL:tbl:input.

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 cadence999While Lightcurvefactory is able to handle directly
finite exposure times, we found that finite exposure (more, strictly speaking
cadence-) time corrections were unnecessary to apply even for the 30-min
cadence light curves, due to the relatively long durations of both the inner
and outer eclipses compared to the cadence times., and dropped the out-of-
eclipse regions of these light curves, keeping only the $\pm
0\aas@@fstack{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 LABEL:tbl:noticed, 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_{\rm out}$ and, in a less
certain way (through the higher order perturbations), about $q_{\rm 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_{\rm out}/q_{\rm 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_{\rm Aa,Ab}/a_{\rm in}$, the outer eclipses provide information about
$R_{\rm Aa,Ab,B}/a_{\rm 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 at 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_{\rm 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_{0}$ 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.

Table 6: Orbital and astrophysical parameters of TICs 133771812 and 176713425
from the joint photodynamical light curve, ETV, SED and PARSEC isochrone
solution.

| TIC 133771812 | TIC 176713425
---|---|---
Orbital elements
| Subsystem
| Aa–Ab | A–B | Aa–Ab | A–B
$t_{0}$ [BJD - 2400000] | $58491.0$ | $58738.5$
$P$ [days] | $12.33390_{-0.00051}^{+0.00045}$ | $243.8900_{-0.0059}^{+0.0048}$ | $1.895107_{-0.000090}^{+0.000101}$ | $52.943_{-0.017}^{+0.019}$
$a$ [R⊙] | $31.66_{-0.30}^{+0.20}$ | $262.0_{-2.2}^{+1.6}$ | $8.456_{-0.069}^{+0.082}$ | $89.42_{-0.74}^{+0.77}$
$e$ | $0.03368_{-0.00044}^{+0.00043}$ | $0.2174_{-0.0058}^{+0.0058}$ | $0.00592_{-0.00056}^{+0.00058}$ | $0.4118_{-0.0070}^{+0.0072}$
$\omega$ [deg] | $230.39_{-0.61}^{+0.55}$ | $194.91_{-0.81}^{+0.78}$ | $105.1_{-2.2}^{+2.5}$ | $95.22_{-0.58}^{+0.87}$
$i$ [deg] | $89.98_{-0.13}^{+0.11}$ | $89.770_{-0.020}^{+0.019}$ | $88.56_{-0.73}^{+0.66}$ | $88.714_{-0.089}^{+0.084}$
$\mathcal{T}_{0}^{\mathrm{inf/sup}}$ [BJD - 2400000] | $58501.1638_{-0.0005}^{+0.0004}$ | ${58509.928_{-0.010}^{+0.009}}^{*}$ | $58739.8183_{-0.0001}^{+0.0001}$ | ${58740.0972_{-0.0041}^{+0.0040}}^{*}$
$\tau$ [BJD - 2400000] | $58499.8871_{-0.0209}^{+0.0189}$ | $58320.200_{-0.554}^{+0.583}$ | $58738.951_{-0.012}^{+0.013}$ | $58687.448_{-0.042}^{+0.054}$
$\Omega$ [deg] | $0.0$ | $-0.22_{-0.11}^{+0.13}$ | $0.0$ | $0.62_{-0.77}^{+0.96}$
$i_{\mathrm{mut}}$ [deg] | $0.31_{-0.14}^{+0.14}$ | $1.00_{-0.52}^{+0.83}$
$\varpi^{\mathrm{dyn}}$ [deg] | $50.39_{-0.61}^{+0.55}$ | $14.91_{-0.81}^{+0.78}$ | $285.1_{-2.2}^{+2.5}$ | $275.22_{-0.58}^{+0.87}$
$i^{\mathrm{dyn}}$ [deg] | $0.25_{-0.12}^{+0.11}$ | $0.06_{-0.03}^{+0.03}$ | $0.83_{-0.43}^{+0.69}$ | $0.17_{-0.09}^{+0.14}$
$\Omega^{\mathrm{dyn}}$ [deg] | $228_{-14}^{+15}$ | $48_{-14}^{+15}$ | $72_{-29}^{+108}$ | $252_{-29}^{+108}$
$i_{\mathrm{inv}}$ [deg] | $89.81_{-0.04}^{+0.03}$ | $88.68_{-0.11}^{+0.14}$
$\Omega_{\mathrm{inv}}$ [deg] | $-0.18_{-0.09}^{+0.10}$ | $0.52_{-0.64}^{+0.80}$
mass ratio $[q=M_{\mathrm{sec}}/M_{\mathrm{pri}}]$ | $0.7627_{-0.0057}^{+0.0054}$ | $0.4488_{-0.0092}^{+0.0099}$ | $0.9334_{-0.0048}^{+0.0052}$ | $0.5150_{-0.0041}^{+0.0042}$
$K_{\mathrm{pri}}$ [km s-1] | $56.23_{-0.57}^{+0.50}$ | $17.26_{-0.31}^{+0.32}$ | $109.01_{-0.83}^{+0.87}$ | $31.92_{-0.35}^{+0.24}$
$K_{\mathrm{sec}}$ [km s-1] | $73.76_{-0.61}^{+0.45}$ | $38.44_{-0.34}^{+0.33}$ | $116.67_{-0.98}^{+1.31}$ | $61.86_{-0.53}^{+0.70}$
Apsidal and nodal motion related parameters
$P_{\mathrm{apse}}$ [year] | $52.81_{-0.81}^{+0.78}$ | $225.0_{-1.6}^{+1.8}$ | $7.09_{-0.12}^{+0.13}$ | $59.71_{-0.85}^{+0.82}$
$P_{\mathrm{apse}}^{\mathrm{dyn}}$ [year] | $23.64_{-0.32}^{+0.31}$ | $35.96_{-0.34}^{+0.35}$ | $4.15_{-0.05}^{+0.05}$ | $8.57_{-0.10}^{+0.09}$
$P_{\mathrm{node}}^{\mathrm{dyn}}$ [year] | $42.81_{-0.51}^{+0.51}$ | $10.01_{-0.11}^{+0.11}$
$\Delta\omega_{\mathrm{3b}}$ [arcsec/cycle] | $1850_{-24}^{+25}$ | $24064_{-233}^{+231}$ | $1231_{-13}^{+14}$ | $21915_{-239}^{+255}$
$\Delta\omega_{\mathrm{GR}}$ [arcsec/cycle] | $0.731_{-0.014}^{+0.009}$ | $0.134_{-0.002}^{+0.002}$ | $2.207_{-0.036}^{+0.043}$ | $0.381_{-0.007}^{+0.007}$
$\Delta\omega_{\mathrm{tide}}$ [arcsec/cycle] | $0.179_{-0.009}^{+0.009}$ | $0.00051_{-0.00003}^{+0.00003}$ | $387_{-17}^{+15}$ | $0.638_{-0.028}^{+0.027}$
Stellar parameters
| Aa | Ab | B | Aa | Ab | B
Relative quantities
fractional radius [$R/a$] | $0.0491_{-0.0005}^{+0.0004}$ | $0.0368_{-0.0008}^{+0.0008}$ | $0.0047_{-0.0001}^{+0.0002}$ | $0.2060_{-0.0018}^{+0.0018}$ | $0.1525_{-0.0025}^{+0.0024}$ | $0.0189_{-0.0008}^{+0.0011}$
temperature relative to $(T_{\mathrm{eff}})_{\mathrm{Aa}}$ | $1$ | $0.8345_{-0.0058}^{+0.0049}$ | $0.8486_{-0.0095}^{+0.0139}$ | $1$ | $1.0039_{-0.0022}^{+0.0041}$ | $1.0040_{-0.0095}^{+0.0052}$
fractional flux [in TESS-band] | $0.3476_{-0.0121}^{+0.0141}$ | $0.1141_{-0.0015}^{+0.0016}$ | $0.1341_{-0.0120}^{+0.0137}$ | $0.3657_{-0.0107}^{+0.0110}$ | $0.2037_{-0.0024}^{+0.0025}$ | $0.3552_{-0.0256}^{+0.0260}$
fractional flux [in Cousins $R_{C}$-band] | $-$ | $-$ | $-$ | $0.3808_{-0.0122}^{+0.0123}$ | $0.2134_{-0.0052}^{+0.0063}$ | $0.3770_{-0.0321}^{+0.0171}$
Physical quantities
$M$ [M⊙] | $1.588_{-0.044}^{+0.028}$ | $1.209_{-0.033}^{+0.026}$ | $1.254_{-0.035}^{+0.034}$ | $1.166_{-0.029}^{+0.037}$ | $1.091_{-0.026}^{+0.028}$ | $1.164_{-0.031}^{+0.025}$
$R$ [R⊙] | $1.554_{-0.017}^{+0.017}$ | $1.165_{-0.034}^{+0.031}$ | $1.225_{-0.041}^{+0.046}$ | $1.743_{-0.027}^{+0.024}$ | $1.290_{-0.027}^{+0.027}$ | $1.692_{-0.072}^{+0.105}$
$T_{\mathrm{eff}}$ [K] | $7706_{-144}^{+168}$ | $6426_{-89}^{+112}$ | $6561_{-114}^{+103}$ | $6033_{-192}^{+78}$ | $6051_{-153}^{+77}$ | $6033_{-160}^{+97}$
$L_{\mathrm{bol}}$ [L⊙] | $7.66_{-0.61}^{+0.67}$ | $2.08_{-0.22}^{+0.23}$ | $2.49_{-0.29}^{+0.35}$ | $3.63_{-0.55}^{+0.20}$ | $2.01_{-0.23}^{+0.13}$ | $3.43_{-0.43}^{+0.42}$
$M_{\mathrm{bol}}$ | $2.56_{-0.09}^{+0.09}$ | $3.97_{-0.11}^{+0.12}$ | $3.78_{-0.14}^{+0.13}$ | $3.37_{-0.06}^{+0.18}$ | $4.01_{-0.07}^{+0.13}$ | $3.43_{-0.12}^{+0.15}$
$M_{V}$ | $2.52_{-0.09}^{+0.09}$ | $3.97_{-0.12}^{+0.12}$ | $3.77_{-0.15}^{+0.14}$ | $3.39_{-0.06}^{+0.20}$ | $4.03_{-0.07}^{+0.15}$ | $3.44_{-0.12}^{+0.16}$
$\log g$ [dex] | $4.253_{-0.008}^{+0.009}$ | $4.388_{-0.015}^{+0.016}$ | $4.358_{-0.021}^{+0.019}$ | $4.022_{-0.007}^{+0.007}$ | $4.254_{-0.012}^{+0.012}$ | $4.043_{-0.047}^{+0.040}$
Global system parameters
$\log$(age) [dex] | $8.760_{-0.053}^{+0.142}$ | $9.720_{-0.021}^{+0.097}$
$[M/H]$ [dex] | $-0.009_{-0.094}^{+0.049}$ | $0.076_{-0.126}^{+0.080}$
$E(B-V)$ [mag] | $0.425_{-0.023}^{+0.024}$ | $0.158_{-0.045}^{+0.022}$
extra light $\ell_{4}$ [in TESS-band] | $0.407_{-0.026}^{+0.016}$ | $0.075_{-0.022}^{+0.025}$
extra light $\ell_{4}$ [in Cousins $R_{C}$-band] | $-$ | $0.026_{-0.015}^{+0.031}$
$(M_{V})_{\mathrm{tot}}$ | $2.03_{-0.10}^{+0.08}$ | $2.39_{-0.06}^{+0.18}$
distance [pc] | $1428_{-27}^{+32}$ | $2097_{-51}^{+51}$

Table 7: Orbital and astrophysical parameters of TICs 185615681 and 287756035
from the joint photodynamical light curve, ETV, SED and PARSEC isochrone
solution.

| TIC 185615681 | TIC 287756035
---|---|---
Orbital elements
| Subsystem
| Aa–Ab | A–B | Aa–Ab | A–B
$t_{0}$ [BJD - 2400000] | $58516.0$ | 58596.5
$P$ [days] | $2.3180401_{-0.0000066}^{+0.0000064}$ | $56.06661_{-0.00066}^{+0.00071}$ | $2.082184_{-0.000082}^{+0.000102}$ | $367.923_{-0.023}^{+0.022}$
$a$ [R⊙] | $11.718_{-0.057}^{+0.073}$ | $104.78_{-0.52}^{+0.57}$ | $8.725_{-0.112}^{+0.049}$ | $318.1_{-4.5}^{+2.1}$
$e$ | $0.010669_{-0.000071}^{+0.000087}$ | $0.09817_{-0.00084}^{+0.00087}$ | $0.0036_{-0.0012}^{+0.0012}$ | $0.235_{-0.025}^{+0.019}$
$\omega$ [deg] | $207.87_{-0.41}^{+0.58}$ | $308.74_{-0.51}^{+0.54}$ | $175_{-31}^{+27}$ | $53.70_{-0.66}^{+0.63}$
$i$ [deg] | $88.06_{-0.12}^{+0.11}$ | $89.826_{-0.097}^{+0.071}$ | $89.27_{-0.83}^{+0.59}$ | $89.095_{-0.028}^{+0.049}$
$\mathcal{T}_{0}^{\mathrm{inf/sup}}$ [BJD - 2400000] | $58518.17078_{-0.00004}^{+0.00004}$ | ${58517.269_{-0.015}^{+0.015}}^{*}$ | $58599.3048_{-0.0008}^{+0.0009}$ | $59311.672_{-0.021}^{+0.039}$
$\tau$ [BJD - 2400000] | $58517.777_{-0.003}^{+0.004}$ | $58496.431_{-0.079}^{+0.082}$ | $58598.76_{-0.18}^{+0.16}$ | $59071.8_{-1.1}^{+1.4}$
$\Omega$ [deg] | $0.0$ | $2.17_{-0.14}^{+0.13}$ | $0.0$ | $1.22_{-0.54}^{+0.64}$
$i_{\mathrm{mut}}$ [deg] | $2.79_{-0.16}^{+0.17}$ | $1.39_{-0.49}^{+0.60}$
$\varpi^{\mathrm{dyn}}$ [deg] | $27.83_{-0.41}^{+0.57}$ | $128.74_{-0.51}^{+0.54}$ | $355_{-31}^{+27}$ | $233.71_{-0.66}^{+0.63}$
$i^{\mathrm{dyn}}$ [deg] | $1.97_{-0.12}^{+0.12}$ | $0.82_{-0.05}^{+0.05}$ | $1.27_{-0.44}^{+0.55}$ | $0.12_{-0.04}^{+0.05}$
$\Omega^{\mathrm{dyn}}$ [deg] | $50.8_{-1.8}^{+1.9}$ | $230.8_{-1.8}^{+1.9}$ | $98_{-37}^{+27}$ | $278_{-37}^{+27}$
$i_{\mathrm{inv}}$ [deg] | $89.30_{-0.08}^{+0.08}$ | $89.12_{-0.06}^{+0.04}$
$\Omega_{\mathrm{inv}}$ [deg] | $1.53_{-0.10}^{+0.10}$ | $1.12_{-0.50}^{+0.59}$
mass ratio $[q=M_{\mathrm{sec}}/M_{\mathrm{pri}}]$ | $0.998_{-0.002}^{+0.002}$ | $0.221_{-0.003}^{+0.002}$ | $1.193_{-0.018}^{+0.019}$ | $0.552_{-0.005}^{+0.004}$
$K_{\mathrm{pri}}$ [km s-1] | $127.76_{-0.63}^{+0.77}$ | $17.21_{-0.16}^{+0.15}$ | $115.44_{-2.27}^{+1.30}$ | $16.01_{-0.27}^{+0.19}$
$K_{\mathrm{sec}}$ [km s-1] | $127.99_{-0.65}^{+0.80}$ | $77.81_{-0.38}^{+0.57}$ | $96.58_{-0.53}^{+0.49}$ | $28.95_{-0.27}^{+0.24}$
$V_{\gamma}$ [km s-1] | $-5.97_{-0.16}^{+0.15}$ | $-3.90_{-0.20}^{+0.27}$
Apsidal and nodal motion related parameters
$P_{\mathrm{apse}}$ [year] | $18.17_{-0.12}^{+0.15}$ | $64.38_{-0.11}^{+0.09}$ | $6.56_{-0.69}^{+0.68}$ | $6429_{-117}^{+142}$
$P_{\mathrm{apse}}^{\mathrm{dyn}}$ [year] | $9.28_{-0.05}^{+0.07}$ | $14.67_{-0.05}^{+0.07}$ | $6.49_{-0.68}^{+0.66}$ | $515_{-9}^{+10}$
$P_{\mathrm{node}}^{\mathrm{dyn}}$ [year] | $18.99_{-0.13}^{+0.09}$ | $560_{-11}^{+10}$
$\Delta\omega_{\mathrm{3b}}$ [arcsec/cycle] | $738.6_{-5.9}^{+4.3}$ | $13564_{-65}^{+47}$ | $25.22_{-0.49}^{+0.45}$ | $2535_{-50}^{+45}$
$\Delta\omega_{\mathrm{GR}}$ [arcsec/cycle] | $2.832_{-0.028}^{+0.035}$ | $0.391_{-0.004}^{+0.004}$ | $1.946_{-0.050}^{+0.022}$ | $0.087_{-0.002}^{+0.002}$
$\Delta\omega_{\mathrm{tide}}$ [arcsec/cycle] | $144.4_{-1.6}^{+1.7}$ | $0.227_{-0.003}^{+0.003}$ | $1112_{-106}^{+132}$ | $0.131_{-0.012}^{+0.014}$
Stellar parameters
| Aa | Ab | B | Aa | Ab | B
Relative quantities
fractional radius [$R/a$] | $0.1749_{-0.0005}^{+0.0005}$ | $0.1744_{-0.0005}^{+0.0005}$ | $0.00749_{-0.00008}^{+0.00008}$ | $0.1129_{-0.0017}^{+0.0018}$ | $0.2707_{-0.0054}^{+0.0061}$ | $0.0212_{-0.0008}^{+0.0006}$
temperature relative to $(T_{\mathrm{eff}})_{\mathrm{Aa}}$ | $1$ | $0.9986_{-0.0014}^{+0.0014}$ | $0.5633_{-0.0057}^{+0.0042}$ | $1$ | $0.9126_{-0.0117}^{+0.0095}$ | $0.8043_{-0.0053}^{+0.0115}$
fractional flux [in TESS-band] | $0.4840_{-0.0035}^{+0.0039}$ | $0.4806_{-0.0036}^{+0.0039}$ | $0.0126_{-0.0007}^{+0.0007}$ | $0.0372_{-0.0015}^{+0.0016}$ | $0.1537_{-0.0065}^{+0.0065}$ | $0.7708_{-0.0174}^{+0.0170}$
Physical quantities
$M$ [M⊙] | $2.009_{-0.030}^{+0.038}$ | $2.005_{-0.029}^{+0.037}$ | $0.888_{-0.013}^{+0.014}$ | $0.935_{-0.023}^{+0.013}$ | $1.116_{-0.053}^{+0.027}$ | $1.133_{-0.056}^{+0.029}$
$R$ [R⊙] | $2.051_{-0.013}^{+0.012}$ | $2.045_{-0.012}^{+0.012}$ | $0.785_{-0.011}^{+0.011}$ | $0.986_{-0.020}^{+0.017}$ | $2.361_{-0.055}^{+0.060}$ | $6.753_{-0.318}^{+0.205}$
$T_{\mathrm{eff}}$ [K] | $9122_{-110}^{+140}$ | $9108_{-109}^{+137}$ | $5136_{-62}^{+64}$ | $6025_{-166}^{+37}$ | $5501_{-207}^{+54}$ | $4833_{-50}^{+53}$
$L_{\mathrm{bol}}$ [L⊙] | $23.11_{-2.45}^{+2.30}$ | $22.92_{-2.49}^{+2.29}$ | $0.449_{-0.045}^{+0.050}$ | $1.14_{-0.11}^{+0.06}$ | $4.59_{-0.80}^{+0.26}$ | $22.59_{-2.87}^{+1.04}$
$M_{\mathrm{bol}}$ | $1.36_{-0.10}^{+0.12}$ | $1.37_{-0.10}^{+0.12}$ | $5.64_{-0.12}^{+0.11}$ | $4.63_{-0.06}^{+0.11}$ | $3.12_{-0.06}^{+0.21}$ | $1.39_{-0.05}^{+0.15}$
$M_{V}$ | $1.39_{-0.06}^{+0.08}$ | $1.40_{-0.06}^{+0.08}$ | $5.79_{-0.15}^{+0.16}$ | $4.67_{-0.06}^{+0.12}$ | $3.25_{-0.07}^{+0.26}$ | $1.72_{-0.05}^{+0.17}$
$\log g$ [dex] | $4.116_{-0.003}^{+0.003}$ | $4.118_{-0.003}^{+0.003}$ | $4.595_{-0.006}^{+0.006}$ | $4.421_{-0.012}^{+0.011}$ | $3.738_{-0.021}^{+0.020}$ | $2.835_{-0.025}^{+0.029}$
Global system parameters
$\log$(age) [dex] | $8.768_{-0.059}^{+0.030}$ | $9.756_{-0.041}^{+0.081}$
$[M/H]$ [dex] | $-0.021_{-0.112}^{+0.175}$ | $-0.262_{-0.122}^{+0.128}$
$E(B-V)$ [mag] | $0.093_{-0.027}^{+0.022}$ | $0.110_{-0.049}^{+0.032}$
extra light $\ell_{4}$ [in TESS-band] | $0.023_{-0.007}^{+0.006}$ | $0.038_{-0.016}^{+0.017}$
$(M_{V})_{\mathrm{tot}}$ | $0.63_{-0.06}^{+0.08}$ | $1.42_{-0.04}^{+0.20}$
distance [pc] | $942_{-13}^{+12}$ | $2561_{-102}^{+77}$

Table 8: Orbital and astrophysical parameters of TIC 321978218 from the joint
photodynamical light curve, ETV, SED and PARSEC isochrone solution.

TIC 321978218
---
Orbital elements
| Subsystem
| Aa–Ab | A–B | AB–C
$t_{0}$ [BJD - 2400000] | $58354.0$
$P$ [days] | $0.570267_{-0.000025}^{+0.000041}$ | $57.53549_{-0.00020}^{+0.00020}$ | $4395_{-84}^{+52}$
$a$ [R⊙] | $3.473_{-0.011}^{+0.009}$ | $92.46_{-0.37}^{+0.36}$ | $1691_{-19}^{+18}$
$e$ | $0.00100_{-0.00027}^{+0.00043}$ | $0.2583_{-0.0016}^{+0.0016}$ | $0.374_{-0.033}^{+0.049}$
$\omega$ [deg] | $140_{-53}^{+22}$ | $122.56_{-0.27}^{+0.28}$ | $68.9_{-4.3}^{+3.1}$
$i$ [deg] | $87.96_{-1.02}^{+0.37}$ | $90.000_{-0.069}^{+0.089}$ | $95.5_{-6.6}^{+9.2}$
$\mathcal{T}_{0}^{\mathrm{inf/sup}}$ [BJD - 2400000] | $58354.2542_{-0.0003}^{+0.0002}$ | ${59072.6568_{-0.0031}^{+0.0033}}^{*}$ | …
$\tau$ [BJD - 2400000] | $58354.061_{-0.077}^{+0.040}$ | $59018.153_{-0.071}^{+0.037}$ | $56799_{-24}^{+13}$
$\Omega$ [deg] | $0.0$ | $0.06_{-0.57}^{+0.79}$ | $-2.2_{-7.2}^{+7.3}$
$(i_{\mathrm{mut}})_{\mathrm{A-AB,A-C}}$ [deg] | $2.17_{-0.40}^{+1.08}$ | $11.3_{-5.4}^{+7.0}$
$(i_{\mathrm{mut}})_{\mathrm{AB-C}}$ [deg] | … | $9.9_{-4.9}^{+6.5}$
$\varpi^{\mathrm{dyn}}$ [deg] | $320_{-53}^{+22}$ | $302.53_{-0.30}^{+0.29}$ | $249.2_{-4.4}^{+3.1}$
$i^{\mathrm{dyn}}$ [deg] | $5.7_{-2.4}^{+3.3}$ | $4.2_{-2.1}^{+2.8}$ | $5.7_{-2.9}^{+3.7}$
$\Omega^{\mathrm{dyn}}$ [deg] | $350_{-33}^{+54}$ | $338_{-62}^{+68}$ | $159_{-57}^{+68}$
$i_{\mathrm{inv}}$ [deg] | $92.3_{-2.8}^{+4.0}$
$\Omega_{\mathrm{inv}}$ [deg] | $-0.9_{-3.2}^{+3.3}$
mass ratio $[q=M_{\mathrm{sec}}/M_{\mathrm{pri}}]$ | $0.782_{-0.015}^{+0.018}$ | $0.853_{-0.007}^{+0.008}$ | $0.052_{-0.005}^{+0.006}$
$K_{\mathrm{pri}}$ [km s-1] | $135.10_{-1.21}^{+1.54}$ | $38.76_{-0.31}^{+0.34}$ | $1.02_{-0.10}^{+0.11}$
$K_{\mathrm{sec}}$ [km s-1] | $172.85_{-2.15}^{+1.79}$ | $45.42_{-0.10}^{+0.09}$ | $19.80_{-0.53}^{+0.49}$
$V_{\gamma}$ [km s-1] | $1.707_{-0.045}^{+0.044}$
Apsidal and nodal motion related parameters
$P_{\mathrm{apse}}$ [year] | $0.88_{-0.07}^{+0.05}$ | $526.9_{-2.6}^{+2.4}$ | $1004_{-84}^{+1138}$
$P_{\mathrm{apse}}^{\mathrm{dyn}}$ [year] | $0.86_{-0.07}^{+0.05}$ | $35.81_{-0.15}^{+0.14}$ | $38.43_{-0.18}^{+0.19}$
$P_{\mathrm{node}}^{\mathrm{dyn}}$ [year] | $38.28_{-0.22}^{+0.23}$
$\Delta\omega_{\mathrm{3b}}$ [arcsec/cycle] | $101.25_{-0.54}^{+0.49}$ | $5700_{-22}^{+25}$ | $405544_{-7670}^{+6085}$
$\Delta\omega_{\mathrm{GR}}$ [arcsec/cycle] | $4.111_{-0.026}^{+0.022}$ | $0.307_{-0.003}^{+0.002}$ | $0.0192_{-0.0006}^{+0.0008}$
$\Delta\omega_{\mathrm{tide}}$ [arcsec/cycle] | $2251_{-119}^{+197}$ | $0.450_{-0.025}^{+0.038}$ | $0.0016_{-0.0001}^{+0.0001}$
Stellar parameters
| Aa | Ab | B | C
Relative quantities
fractional radius [$R/a$] | $0.2888_{-0.0035}^{+0.0041}$ | $0.2370_{-0.0027}^{+0.0055}$ | $0.0211_{-0.0003}^{+0.0003}$ | $0.00029_{-0.00001}^{+0.00001}$
temperature relative to $(T_{\mathrm{eff}})_{\mathrm{Aa}}$ | $1$ | $0.8522_{-0.0073}^{+0.0075}$ | $1.2755_{-0.0093}^{+0.0095}$ | $0.5596_{-0.0156}^{+0.0148}$
fractional flux [in TESS-band] | $0.0764_{-0.0019}^{+0.0019}$ | $0.0252_{-0.0011}^{+0.0012}$ | $0.6998_{-0.0199}^{+0.0169}$ | $0.0009_{-0.0002}^{+0.0003}$
fractional flux [in SWASP-band] | $0.0672_{-0.0032}^{+0.0030}$ | $0.0169_{-0.0013}^{+0.0013}$ | $0.8795_{-0.0403}^{+0.0254}$ | $0.0004_{-0.0001}^{+0.0001}$
Physical quantities
$M$ [M⊙] | $0.970_{-0.018}^{+0.013}$ | $0.757_{-0.006}^{+0.009}$ | $1.473_{-0.023}^{+0.025}$ | $0.165_{-0.016}^{+0.018}$
$R$ [R⊙] | $1.003_{-0.013}^{+0.014}$ | $0.824_{-0.009}^{+0.017}$ | $1.949_{-0.031}^{+0.027}$ | $0.492_{-0.025}^{+0.024}$
$T_{\mathrm{eff}}$ [K] | $4969_{-48}^{+83}$ | $4242_{-53}^{+63}$ | $6326_{-59}^{+147}$ | $2791_{-95}^{+78}$
$L_{\mathrm{bol}}$ [L⊙] | $0.553_{-0.025}^{+0.035}$ | $0.198_{-0.014}^{+0.017}$ | $5.46_{-0.24}^{+0.45}$ | $0.013_{-0.002}^{+0.003}$
$M_{\mathrm{bol}}$ | $5.41_{-0.07}^{+0.05}$ | $6.53_{-0.09}^{+0.08}$ | $2.93_{-0.09}^{+0.05}$ | $9.47_{-0.20}^{+0.23}$
$M_{V}$ | $5.71_{-0.11}^{+0.07}$ | $7.36_{-0.15}^{+0.15}$ | $2.92_{-0.09}^{+0.05}$ | $13.08_{-0.53}^{+0.65}$
$\log g$ [dex] | $4.421_{-0.013}^{+0.009}$ | $4.485_{-0.019}^{+0.011}$ | $4.026_{-0.011}^{+0.011}$ | $4.269_{-0.016}^{+0.018}$
Global system parameters
$\log$(age) [dex] | $7.326_{-0.049}^{+0.030}$ | $7.046_{-0.033}^{+0.024}$ | $7.326_{-0.049}^{+0.030}$
$[M/H]$ [dex] | $-0.003_{-0.102}^{+0.118}$
$E(B-V)$ [mag] | $0.015_{-0.011}^{+0.032}$
extra light $\ell_{4}$ [in TESS-band] | $0.197_{-0.017}^{+0.021}$
extra light $\ell_{4}$ [in SWASP-band] | $0.035_{-0.028}^{+0.045}$
$(M_{V})_{\mathrm{tot}}$ | $2.82_{-0.09}^{+0.05}$
distance [pc] | $636_{-9}^{+10}$

Table 9: Orbital and astrophysical parameters of TICs 323486857 and 650024463
from the joint photodynamical light curve, ETV, SED and PARSEC isochrone
solution.

| TIC 323486857 | TIC 650024463
---|---|---
Orbital elements
| Subsystem | Subsystem
| Aa–Ab | A–B | Aa–Ab | A–B
$t_{0}$ [BJD - 2400000] | 58569.0 | $58320.0$
$P$ [days] | $0.88502_{-0.00036}^{+0.00015}$ | $41.4268_{-0.0021}^{+0.0016}$ | $7.197817_{-0.000028}^{+0.000026}$ | $108.7251_{-0.0013}^{+0.0013}$
$a$ [R⊙] | $5.362_{-0.053}^{+0.035}$ | $81.51_{-0.60}^{+0.61}$ | $18.21_{-0.12}^{+0.07}$ | $127.09_{-0.87}^{+0.47}$
$e$ | $0.0021_{-0.0009}^{+0.0012}$ | $0.0066_{-0.0037}^{+0.0042}$ | $0.02570_{-0.00026}^{+0.00024}$ | $0.32317_{-0.00011}^{+0.00011}$
$\omega$ [deg] | $151_{-66}^{+77}$ | $241_{-139}^{+48}$ | $28.83_{-0.49}^{+0.53}$ | $351.28_{-0.10}^{+0.09}$
$i$ [deg] | $84.9_{-0.4}^{+2.3}$ | $84.59_{-0.45}^{+0.28}$ | $90.233_{-0.055}^{+0.059}$ | $90.099_{-0.090}^{+0.069}$
$\mathcal{T}_{0}^{\mathrm{inf/sup}}$ [BJD - 2400000] | $58569.7694_{-0.0003}^{+0.0003}$ | $58595.020_{-0.028}^{+0.044}$ | $58330.8185_{-0.0002}^{+0.0002}$ | ${58642.2864_{-0.0086}^{+0.0085}}^{*}$
$\tau$ [BJD - 2400000] | $58569.47_{-0.16}^{+0.20}$ | $58577_{-20}^{+17}$ | $58325.9422_{-0.0096}^{+0.0100}$ | $58623.745_{-0.024}^{+0.025}$
$\Omega$ [deg] | $0.0$ | $1.9_{-2.5}^{+1.2}$ | $0.0$ | $0.33_{-0.26}^{+0.30}$
$i_{\mathrm{mut}}$ [deg] | $2.3_{-0.7}^{+1.2}$ | $0.37_{-0.22}^{+0.31}$
$\varpi^{\mathrm{dyn}}$ [deg] | $339_{-70}^{+73}$ | $76_{-142}^{+47}$ | $208.83_{-0.49}^{+0.53}$ | $171.28_{-0.10}^{+0.09}$
$i^{\mathrm{dyn}}$ [deg] | $2.05_{-0.60}^{+1.07}$ | $0.28_{-0.08}^{+0.14}$ | $0.30_{-0.18}^{+0.24}$ | $0.07_{-0.04}^{+0.06}$
$\Omega^{\mathrm{dyn}}$ [deg] | $98_{-37}^{+27}$ | $285_{-20}^{+90}$ | $112_{-12}^{+25}$ | $292_{-12}^{+25}$
$i_{\mathrm{inv}}$ [deg] | $84.68_{-0.47}^{+0.41}$ | $90.125_{-0.068}^{+0.057}$
$\Omega_{\mathrm{inv}}$ [deg] | $1.7_{-2.2}^{+1.0}$ | $0.27_{-0.21}^{+0.24}$
mass ratio $[q=M_{\mathrm{sec}}/M_{\mathrm{pri}}]$ | $0.877_{-0.021}^{+0.028}$ | $0.606_{-0.014}^{+0.008}$ | $0.780_{-0.003}^{+0.004}$ | $0.490_{-0.004}^{+0.004}$
$K_{\mathrm{pri}}$ [km s-1] | $142.41_{-2.07}^{+3.21}$ | $37.44_{-0.56}^{+0.29}$ | $56.12_{-0.39}^{+0.26}$ | $20.55_{-0.16}^{+0.14}$
$K_{\mathrm{sec}}$ [km s-1] | $162.74_{-2.52}^{+1.82}$ | $61.80_{-0.75}^{+0.59}$ | $71.93_{-0.50}^{+0.40}$ | $41.94_{-0.28}^{+0.23}$
Apsidal and nodal motion related parameters
$P_{\mathrm{apse}}$ [year] | $0.58_{-0.05}^{+0.15}$ | $139.0_{-1.0}^{+1.5}$ | $15.45_{-0.08}^{+0.08}$ | $62.93_{-0.09}^{+0.08}$
$P_{\mathrm{apse}}^{\mathrm{dyn}}$ [year] | $0.56_{-0.04}^{+0.14}$ | $14.80_{-0.07}^{+0.14}$ | $6.87_{-0.03}^{+0.03}$ | $10.34_{-0.03}^{+0.03}$
$P_{\mathrm{node}}^{\mathrm{dyn}}$ [year] | $16.56_{-0.19}^{+0.10}$ | $12.38_{-0.05}^{+0.04}$
$\Delta\omega_{\mathrm{3b}}$ [arcsec/cycle] | $356.5_{-4.2}^{+2.2}$ | $9917_{-94}^{+50}$ | $3716_{-17}^{+15}$ | $37300_{-109}^{+102}$
$\Delta\omega_{\mathrm{GR}}$ [arcsec/cycle] | $4.068_{-0.077}^{+0.053}$ | $0.429_{-0.006}^{+0.007}$ | $0.710_{-0.009}^{+0.006}$ | $0.169_{-0.002}^{+0.001}$
$\Delta\omega_{\mathrm{tide}}$ [arcsec/cycle] | $5243_{-1121}^{+481}$ | $12.6_{-3.2}^{+5.0}$ | $0.378_{-0.024}^{+0.026}$ | $0.0014_{-0.0001}^{+0.0001}$
Stellar parameters
| Aa | Ab | B | Aa | Ab | B
Relative quantities
fractional radius [$R/a$] | $0.3479_{-0.0155}^{+0.0092}$ | $0.2514_{-0.0098}^{+0.0181}$ | $0.0925_{-0.0049}^{+0.0083}$ | $0.0531_{-0.0008}^{+0.0009}$ | $0.0369_{-0.0003}^{+0.0002}$ | $0.00593_{-0.00006}^{+0.00006}$
temperature relative to $(T_{\mathrm{eff}})_{\mathrm{Aa}}$ | $1$ | $0.9816_{-0.0066}^{+0.0109}$ | $0.7741_{-0.0128}^{+0.0145}$ | $1$ | $0.8104_{-0.0038}^{+0.0041}$ | $0.8998_{-0.0067}^{+0.0063}$
fractional flux [in TESS-band] | $0.0972_{-0.0080}^{+0.0071}$ | $0.0496_{-0.0035}^{+0.0050}$ | $0.6820_{-0.0609}^{+0.0569}$ | $0.1013_{-0.0028}^{+0.0029}$ | $0.0207_{-0.0007}^{+0.0008}$ | $0.0411_{-0.0020}^{+0.0021}$
Physical quantities
$M$ [M⊙] | $1.408_{-0.043}^{+0.021}$ | $1.223_{-0.028}^{+0.052}$ | $1.588_{-0.038}^{+0.050}$ | $0.877_{-0.018}^{+0.012}$ | $0.684_{-0.014}^{+0.008}$ | $0.765_{-0.017}^{+0.010}$
$R$ [R⊙] | $1.873_{-0.105}^{+0.046}$ | $1.344_{-0.057}^{+0.110}$ | $7.541_{-0.455}^{+0.750}$ | $0.967_{-0.015}^{+0.015}$ | $0.673_{-0.009}^{+0.005}$ | $0.754_{-0.011}^{+0.008}$
$T_{\mathrm{eff}}$ [K] | $6550_{-44}^{+61}$ | $6434_{-60}^{+86}$ | $5068_{-71}^{+76}$ | $5682_{-21}^{+28}$ | $4608_{-31}^{+31}$ | $5114_{-39}^{+43}$
$L_{\mathrm{bol}}$ [L⊙] | $5.80_{-0.69}^{+0.34}$ | $2.78_{-0.32}^{+0.64}$ | $33.51_{-2.44}^{+5.16}$ | $0.877_{-0.020}^{+0.021}$ | $0.182_{-0.005}^{+0.006}$ | $0.349_{-0.015}^{+0.015}$
$M_{\mathrm{bol}}$ | $2.86_{-0.06}^{+0.14}$ | $3.66_{-0.23}^{+0.13}$ | $0.96_{-0.16}^{+0.08}$ | $4.91_{-0.03}^{+0.03}$ | $6.62_{-0.03}^{+0.03}$ | $5.91_{-0.04}^{+0.05}$
$M_{V}$ | $2.84_{-0.06}^{+0.14}$ | $3.65_{-0.23}^{+0.13}$ | $1.22_{-0.15}^{+0.10}$ | $4.99_{-0.03}^{+0.03}$ | $7.14_{-0.05}^{+0.05}$ | $6.15_{-0.06}^{+0.06}$
$\log g$ [dex] | $4.043_{-0.019}^{+0.036}$ | $4.267_{-0.050}^{+0.028}$ | $2.882_{-0.067}^{+0.049}$ | $4.409_{-0.015}^{+0.014}$ | $4.617_{-0.003}^{+0.004}$ | $4.566_{-0.006}^{+0.006}$
Global system parameters
$\log$(age) [dex] | $9.374_{-0.038}^{+0.013}$ | $10.001_{-0.033}^{+0.039}$
$[M/H]$ [dex] | $0.015_{-0.058}^{+0.018}$ | $-0.107_{-0.055}^{+0.036}$
$E(B-V)$ [mag] | $0.340_{-0.020}^{+0.020}$ | $0.088$
extra light $\ell_{4}$ [in TESS-band] | $0.172_{-0.063}^{+0.064}$ | $0.837_{-0.003}^{+0.003}$
$(M_{V})_{\mathrm{tot}}$ | $0.93_{-0.15}^{+0.07}$ | $4.56_{-0.01}^{+0.01}$
distance [pc] | $1434_{-74}^{+123}$ | $1126$

Table 10: Derived ephemerides for the seven triple systems to be used for
planning future observations.

TIC ID | 133771812 | 176713425 | 185615681 | 287756035 | 321978218 | 323486857 | 650024463
---|---|---|---|---|---|---|---
| Inner binary
$P$ | 12.36949 | 1.898728 | 2.3185738 | 2.0814225 | 0.5698104 | 0.883992 | 7.1971
| 12.36844 |  | 2.3185448 |  |  |  | 7.1969
$\mathcal{T}_{0}$ | 60 035.027 | 60 317.664 | 60 055.382 | 60 127.070 | 60 434.640 | 60 150.345 | 60 252.278
| 60 020.719 |  | 60 056.539 |  |  |  | 60 248.733
$\mathcal{A}_{\mathrm{ETV}}$ | 0.030 | 0.006 | 0.0007 | 0.0025 | 0.0013 | 0.0008 | 0.028
$D$ | 0.322 | 0.228 | 0.269 | 0.283 | 0.085 | 0.168 | 0.307
| Wide binary (third body eclipses)
$P$ | 242.3 | 52.55 | 55.84 | 367.64 | 57.52 | 41.356 | 107.95
$\mathcal{T}_{0}^{\mathrm{inf}}$ | 58 600.1: | (58762.2:) | 58 547.5 | 59 311.55 | 59 095.85 | 58 595.02 | 59 579.5
$D^{\mathrm{inf}}$ | 6.5 | 6.3: | 1.45 | 4.21 | 1.01 | 1.29 | 2.70
$\mathcal{T}_{0}^{\mathrm{sup}}$ | 58 510.4 | 58 740.1 | (58 572.9) | 59 470.28: | 59 072.45 | (58 615.71) | 59 503.8
$D^{\mathrm{sup}}$ | 6.3 | 0.89 | 1.90 | 3.90: | 0.68 | 1.27 | 3.18

Notes. (a) For the inner binaries: $P$, $\mathcal{T}_{0}$,
$\mathcal{A}_{\mathrm{ETV}}$, $D$ are the period, reference time of a primary
minimum, half-amplitude of the ETV curve, and the full duration of an eclipse,
respectively. $\mathcal{T}_{0}$ is given in BJD – 2 400 000, while the other
quantities are in days. For all those triples where the inner eccentricities
are very small and, hence, the shifts of the secondary eclipses relative to
phase 0.5 are negligible (quantitatively, they are much smaller than the full
durations of the individual eclipses), the same reference times and periods
can be used to predict the times of the secondary eclipses. In the case of the
three eccentric EBs we give a separate period and reference time for the
secondary eclipses, listing them below the primary ephemerides. (b) For the
outer orbits we give separate reference times for the third body eclipses
around the inferior and superior conjunctions of the tertiary component. The
eclipse durations, $D$, of the third-body eclipses do not give the extent of
any specific third body events. Rather $D$ represents the time difference
corresponding to the very first and last moments around a given third-body
conjunction when the first/last contact of a third-body event may occur).
Double dots (:) call attention to the less certain superior/inferior
conjunction times at those types of third-body events (i.e., primary vs.
secondary outer eclipses) because they were not observed by TESS. Conjunction
data, in parentheses, indicate that only very shallow third-body eclipses may
occur which can hardly be observed with ground-based instruments.

###  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_{\mathrm{out}}\sim 244$ d, but allowed also
for half that period, namely, $P_{\mathrm{out}}\sim 122$ d. We produced
preliminary photodynamical solutions for both outer periods, and we found that
the true period should actually be about $P_{\mathrm{out}}\sim 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_{\mathrm{Aa}}=1.58\pm 0.04\,\mathrm{M}_{\sun}$ and
$M_{\mathrm{Ab}}=1.20\pm 0.03\,\mathrm{M}_{\sun}$, respectively) and,
correspondingly, the hottest and coolest ($T_{\mathrm{Aa}}=7700\pm 160$ K and
$T_{\mathrm{Ab}}=6410\pm 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_{\mathrm{B}}=1.25\pm
0.03\,\mathrm{M}_{\sun}$ and $T_{\mathrm{B}}=6550\pm 110$ K. The inferred
masses, in the absence of RV data, have relatively large uncertainties of
$\sim$$3-4$%. Relative quantities, such as the mass ratios, however, can be
deduced with better accuracies. These are found to be
$q_{\mathrm{in}}=0.763\pm 0.006$ and $q_{\mathrm{out}}=0.449\pm 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\aas@@fstack{d}33$. Thus,
despite the fact that the inferred outer period
($P_{\mathrm{out}}=243\aas@@fstack{d}9$) is the second longest in the current
set, this target is still the second tightest
($P_{\mathrm{out}}/P_{\mathrm{in}}\simeq 19.8$) among these seven triply
eclipsing triples. Due the relative large size of the inner EB and, hence, the
small fractional radii
($r_{\mathrm{Aa}}=R_{\mathrm{Aa}}/a_{\mathrm{in}}=0.0491\pm 0.0005$ and
$r_{\mathrm{Ab}}=0.0368\pm 0.0008$), tidal effects are certainly negligible.
Hence, not surprisingly, this inner EB has the largest eccentricity with
$e_{\mathrm{in}}=0.0337\pm 0.0004$ among our present sample. The combination
of inner and outer ($e_{\mathrm{out}}=0.217\pm 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$\times$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)101010We note, however, that the last four eclipse times were obtained
only after the photodynamical analysis was finalized. Despite this, they are
plotted in Fig. 11 and also tabulated in Table 12. These additional
measurements are in perfect accord with the photodynamical model. 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_{\mathrm{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-$\sigma$ statistical uncertainties ($M_{\mathrm{Aa}}=1.17\pm
0.03\,\mathrm{M}_{\sun}$ vs. $M_{\mathrm{B}}=1.16\pm
0.03\,\mathrm{M}_{\sun}$), and the secondary of the inner EB has a lower mass
by only $6-7\%$ ($M_{\mathrm{Ab}}=1.09\pm 0.03\,\mathrm{M}_{\sun}$).
(Referring to the more precise relative quantities, we find:
$q_{\mathrm{in}}=0.933\pm 0.005$ and $q_{\mathrm{out}}=0.515\pm 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-$\sigma$ 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_{\mathrm{Aa}}=1.555\pm 0.015\,\mathrm{R}_{\sun}$ vs.
$R_{\mathrm{Ab}}=1.137\pm 0.045\,\mathrm{R}_{\sun}$ (or, again citing the
precise fractional radii, we find
$r_{\mathrm{Aa}}=R_{\mathrm{Aa}}/a_{\mathrm{in}}=0.2060\pm 0.0018$ vs.
$r_{\mathrm{Ab}}=0.1525\pm 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_{\mathrm{out}}/P_{\mathrm{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 nominal111111We say ’nominal’ because, due
to the nearly circular shape of the inner orbit, one cannot really expect
observable apsidal motion. apsidal advance rate in the inner EB is of the same
order as the dynamically forced one
($\Delta\omega_{\mathrm{in}}^{\mathrm{tide}}=387^{\prime\prime}\pm
15^{\prime\prime}\,\mathrm{cycle}^{-1}$ vs.
$\Delta\omega_{\mathrm{in}}^{\mathrm{3b}}=1230^{\prime\prime}\pm
15^{\prime\prime}\,\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^{\mathrm{spec}}=0.996\pm 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 $\sim$7000 days), determined an apsidal
motion period of $P_{\mathrm{apse}}=157\pm 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_{2}$, 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 then available 20 data points, and
got a much shorter apsidal motion period of $P_{\mathrm{apse}}=22.89\pm 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_{\mathrm{Aa}}=2.013\pm
0.034\,\mathrm{M}_{\sun}$; $M_{\mathrm{Aa}}=2.008\pm
0.033\,\mathrm{M}_{\sun}$) 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_{\mathrm{B}}=0.888\pm
0.014\,\mathrm{M}_{\sun}$) than either of the binary members. (Referring again
to the more precise relative quantities, we find: $q_{\mathrm{in}}=0.998\pm
0.002$ and $q_{\mathrm{out}}=0.221\pm 0.003$, respectively.) In such a way,
the contribution of the third component to the total system flux is negligibly
small, only $\sim$$1\%$. On the other hand, the depths of the primary third-
body eclipses may reach $\sim$$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_{\mathrm{in}}=0.01067\pm 0.00009$.
This value agrees well with the results of Kim et al. (2018) who obtained
$e_{\mathrm{in}}=0.01068\pm 0.00008$.121212One should keep in mind, however,
that the two values cannot be perfectly compared, as our orbital elements are
instantaneous osculating elements, while the ones listed in the Kim et al.
(2018) paper are some kind of average elements, deduced from the observations.
The outer orbit was also found to be slightly eccentric with
$e_{\mathrm{out}}=0.0982\pm 0.0008$. This is the second lowest outer
eccentricity in our current sample.

Despite the moderately tight nature of the triple system
($P_{\mathrm{out}}/P_{\mathrm{in}}=24.2$) the full-amplitudes of the primary
and secondary ETVs, on the $P_{\mathrm{out}}$-timescale, remain under $1.5$
min ($\sim 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 ($\Delta\omega_{\mathrm{in}}^{\mathrm{3b}}=738^{\prime\prime}\pm
5^{\prime\prime}\,\mathrm{cycle}^{-1}$ vs.
$\Delta\omega_{\mathrm{in}}^{\mathrm{tide}}=144^{\prime\prime}\pm
2^{\prime\prime}\,\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_{2}=0.005$ and $k_{3}=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
$\Delta\omega_{\mathrm{in}}^{\mathrm{GR}}=2\aas@@fstack{\prime\prime}83\pm
0\aas@@fstack{\prime\prime}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 $\sim$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_{\rm
mut}=1.4^{\circ}$.

The outer orbit has a quite typical eccentricity for these systems of
$e_{\mathrm{out}}=0.23\pm 0.02$, while the observable argument of the outer
periastron was found to be $\omega_{\mathrm{out}}=53\aas@@fstack{\circ}7\pm
0\aas@@fstack{\circ}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_{\mathrm{out}}\cos\omega_{\mathrm{out}}$ through the offset of the other
type of third-body eclipse from an orbital phase of $0\aas@@fstack{p}5$. In
addition to this fact, due the lack of any usable ETV curve (which would offer
additional constraints on $e_{\mathrm{out}}$ and $\omega_{\mathrm{out}}$), our
primary source for information about $e_{\mathrm{out}}$ and
$\omega_{\mathrm{out}}$ is only the poorly covered RV curve of the tertiary
(see Fig. 17). Thus, the obtained $e_{\mathrm{out}}$ and
$\omega_{\mathrm{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_{\mathrm{out}}\cos\omega_{\mathrm{out}}$
was a direct, adjustable parameter) resulted in a median posterior of
$e_{\mathrm{out}}\cos\omega_{\mathrm{out}}=0.14\pm 0.01$. Taking into account
that the phase of the superior conjunction relative to the inferior
conjunction can be approximated as
$\phi_{\mathrm{sup}}\simeq\frac{1}{2}-\frac{2}{\pi}e\cos\omega$ (see, e.g.,
Sterne, 1939), one can calculate that the other type of third-body eclipse
should occur around phase $0\aas@@fstack{p}411\pm 0\aas@@fstack{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_{\rm out}\cos\omega_{\rm 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_{\mathrm{B}}=1.133\pm
0.04\,\mathrm{M}_{\sun}$ and $M_{\mathrm{Ab}}=1.116\pm
0.04\,\mathrm{M}_{\sun}$, are evolved to very different levels131313$M_{B}$ is
only 0.017 M⊙ greater than $M_{\rm Ab}$, but is within 1-$\sigma$ of being
equal. However, $M_{B}$ could be as much as 0.14 M⊙ more massive at the
1.5-$\sigma$ level. Also, using the more accurate mass ratios, we find that
$M_{\mathrm{Ab}}/M_{\mathrm{B}}=0.986\pm 0.009$, which is less than unity at
the 1.5-$\sigma$ level.. The slightly more massive tertiary is quite evolved
at $R_{\mathrm{B}}=6.7\pm 0.3\,\mathrm{R}_{\sun}$, 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_{\mathrm{B}}=2.36\pm 0.06\,\mathrm{R}_{\sun}$. 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_{\mathrm{Aa}}=0.93\pm 0.02\,\mathrm{M}_{\sun}$ and $R_{\mathrm{Aa}}=0.98\pm
0.02\,\mathrm{R}_{\sun}$ 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_{\rm eff}\simeq 4830$ K. Only
the tertiary in TIC 323486857 is somewhat more evolved ($R\simeq 7.5$ R⊙)
among our current set of triples.

The system is the opposite of ‘tight’ with $P_{\rm out}/P_{\rm in}\simeq
368/2.1\simeq 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 modest 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_{\mathrm{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_{\mathrm{eff}}=5974\pm 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_{\mathrm{sp}}(M_{\mathrm{A}})=0.506\pm
0.003\,\mathrm{M}_{\sun}$. On the other hand, the preliminary analytic LTTE
ETV solution provides a mass function of
$f_{\mathrm{ETV}}(M_{\mathrm{B}})=0.290\pm 0.030\,\mathrm{M}_{\sun}$. Writing
these out explicitly

$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}$
(1)

and,

$f_{\mathrm{sp}}(M_{\mathrm{A}})=M_{\mathrm{A}}\sin^{3}i_{\mathrm{out}}\left(\frac{q^{-1}_{\mathrm{out}}}{1+q^{-1}_{\mathrm{out}}}\right)^{2}=M_{\mathrm{A}}\sin^{3}i_{\mathrm{out}}\left(\frac{1}{1+q_{\mathrm{out}}}\right)^{2},$
(2)

one can readily see that

$M_{\mathrm{B}}\sin^{3}i_{\mathrm{out}}=f^{1/3}_{\mathrm{ETV}}\left(f^{1/3}_{\mathrm{ETV}}+f^{1/3}_{\mathrm{sp}}\right)^{2}=1.41\pm
0.10\,\mathrm{M}_{\sun}.$ (3)

Thus, one may immediately realize that a star with an effective temperature of
$T_{\mathrm{eff}}\simeq 5900$ K and stellar mass of $m_{\mathrm{B}}\approx
1.4\,\mathrm{M}_{\sun}$ cannot be a main sequence object (we would expect a
range of $T_{\rm eff}\approx 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 $\sim$$4.4$ together temperatures of less than
$5000$ K. In other words, the fractional radii of the EB members141414For the
connection between the fractional radii and the surface gravity see, e.g.,
Southworth 2004; Hajdu et al. 2017. (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}\approx 1.8\,\mathrm{M}_{\sun}$ and
corresponding effective temperature of $T_{\mathrm{eff}}\approx 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 cumulative 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
evolved151515This qualifier is necessary for $T_{\mathrm{eff}}\approx 6000$ K
in the case of such a massive star. $\sim$$1.4\,\mathrm{M}_{\sun}$ 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_{\mathrm{B}}\approx 11$ Myr vs. $t_{\mathrm{Aa,Ab,C}}\approx 21$ Myr. Such
a small difference of $\approx 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 $\sim$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 $\approx 4400$ d, the outer
eccentricity is $\sim$0.37, and the mutual inclination with the inner triple
is $\sim$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_{\mathrm{C}}=0.16\,\mathrm{M}_{\odot}$. 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_{\mathrm{C}}=2790\pm 100$ K) follows from its mass and age. The highly
bloated radius ($R_{\mathrm{C}}=0.49\pm 0.02\,\mathrm{R}_{\sun}$) 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\aas@@fstack{\circ}5\pm 1\degr$. 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_{\rm
B}=7.5\,\mathrm{R}_{\odot}$), while the primary EB star has not evolved much.

Figure 19: Parameter $e_{\rm out}$ plotted against the fractional radius of
the tertiary component $r_{\mathrm{B}}=R_{\rm B}/a_{\mathrm{out}}$ 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 $e_{\rm out}\simeq\exp(-42r_{\rm B})$.
Exceptions are the two systems marked with a faint green ellipse with very low
eccentricities, but not large values of $r_{\rm B}$. Numbers next to
individual sources are the TIC numbers. We omit error bars on the systems with
$e_{\rm out}\gtrsim 0.1$ for increased clarity.

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\pm 0.0040$), and (iii) a fairly large
tertiary fractional radius $r_{\mathrm{B}}=R_{\rm B}/a_{\rm out}=0.093\pm
0.007$. In Fig. 19 we plot the essence of this information in the plane of
$e_{\rm out}$–$r_{\rm 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 $\sim$7-13 times the radius
of the tertiary, while the tidally induced decay of the eccentricity goes as a
high power of $R_{\rm B}/a_{\rm 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_{\rm B}$, seem otherwise perfectly normal systems with well
measured parameters–we presume these were born circular.

####  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.

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_{\rm RP}$ band
fluxes of $\approx 0.26$ and $\approx 0.19$, respectively, relative to the
brightest star.

Figure 20: Positional information for TIC 650024463 and its two neighbors.
Upper panel: A $4\times 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-$\sigma$
confidence intervals. As seen from the panel, the measurements clearly rule
out TIC 302965293 as a potential source of the detected eclipses.

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_{\mathrm{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).

Figure 21: Numerically integrated models for TIC 650024463 for the entire
$21^{\mathrm{st}}$ 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 $P_{\mathrm{out}}$ 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 ($e_{\mathrm{in}}$ – 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.

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\aas@@fstack{d}20$,
$P_{\mathrm{out}}=108\aas@@fstack{d}73$ and, hence,
$P_{\mathrm{out}}/P_{\mathrm{in}}\approx 15.1$. From a dynamical perspective,
this triple star system is very similar to the three triples (KIC 9714358 –
$P_{\mathrm{in}}=6\aas@@fstack{d}47$, $P_{\mathrm{out}}=104\aas@@fstack{d}08$;
KIC 5771589 – $10\aas@@fstack{d}68$, $113\aas@@fstack{d}87$; TIC 219885468 –
$7\aas@@fstack{d}51$, $111\aas@@fstack{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\aas@@fstack{\circ}5$), and have (ii)
small, but significantly non-zero inner eccentricities ($0.0036\lesssim
e_{\mathrm{in}}\lesssim 0.042$), and (iii) moderate outer ($0.16\lesssim
e_{\mathrm{out}}\lesssim 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 $21^{st}$ 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 ($\omega_{\mathrm{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.

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 $21^{st}$ 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 $\sim$$12.9$ yr-period (inner)
eccentricity cycles. The lower right panel of Fig. 21 illustrates that the
inner eccentricity (again, apart from the $P_{\mathrm{out}}$-period mid-
timescale variations) oscillates between $(e_{\mathrm{in}})_{\mathrm{min}}\sim
0.015$ and $(e_{\mathrm{in}})_{\mathrm{max}}\sim 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\aas@@fstack{\circ}41\pm 0\aas@@fstack{\circ}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
($\Delta\equiv\omega^{\mathrm{dyn}}_{\mathrm{in}}-\omega^{\mathrm{dyn}}_{\mathrm{out}}$)
and reach their extrema (i.e., minima and maxima of $e_{\mathrm{in}}$’s) when
$\Delta=0\degr$ or $180\degr$ and, hence, the octupole-order dominant
perturbation term in $e_{\mathrm{in}}$, which is proportional to $\sin\Delta$
(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.

Figure 22: Statistical plots for properties of 33 triply eclipsing triples
uniformly analyzed (see text for references). Top-row panels: $e_{\rm out}$
vs. $P_{\rm out}$ and $i_{\rm out}$ vs. $i_{\rm mut}$. In the upper right
panel the vertical lines denote the transition to the ZLK cycles, and to
retrograde orbits, respectively. Bottom-row panels: $q_{\rm in}$ vs. $q_{\rm
out}$ and $R_{\rm B}$ vs. $M_{\rm B}$. The sloped dashed lines in the latter
plot are for $R_{\rm B}$ =1 MB and = 5 MB (both expressed in solar units), as
rough guides of unevolved and quite evolved stars, respectively.

##  8 Summary, discussion, and conclusions