The Compact Triply Eclipsing Triple Star TIC 209409435 Discovered with TESS

We report the discovery in $TESS$ Sectors 3 and 4 of a compact triply eclipsing triple star system. TIC 209409435 is a previously unknown eclipsing binary with a period of 5.717 days, and the presence of a third star in an outer eccentric orbit of 121.872 day period was found from two sets of third-body eclipses and from eclipse timing variations. The latter exhibit signatures of strong 3rd-body perturbations. After the discovery, we obtained follow-up ground-based photometric observations of several binary eclipses as well as another of the third-body eclipses. We carried out comprehensive analyses, including the simultaneous photodynamical modelling of $TESS$ and ground-based lightcurves (including both archival WASP data, and our own follow-up measurements), as well as eclipse timing variation curves. Also, we have included in the simultaneous fits multiple star spectral energy distribution data and theoretical PARSEC stellar isochrones. We find that the inner binary consists of near twin stars of mass 0.90 $M_\odot$ and radius 0.88 $R_\odot$. The third star is just 9% more massive and 18% larger in radius. The inner binary has a rather small eccentricity while the outer orbit has $e = 0.40$. The inner binary and outer orbit have inclination angles within 0.1$^\circ$ and 0.2$^\circ$ of 90$^\circ$, respectively. The mutual inclination angle is $\lesssim 1/4^\circ$. All of these results were obtained without radial velocity observations.

factor is important when trying to determine system parameters on the timescale of satellite observations (months to years) for CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010), K2 (Howell et al. 2014), and TESS (Ricker et al. 2015), not to mention other human timescales involving theses and careers in astronomy.
In this and previous work on hierarchical stellar systems we have focused on what can be learned from photometric measurements, as opposed to radial velocity studies. There are four principle photometric effects that contribute to a successful search for, and subsequent study of, these compact hierarchical triple systems: (1) binary eclipse timing variations ('ETVs') due to the classical Roemer delays (or light travel time effects); (2) ETVs due to what is referred to as "dynamical delays", i.e., where the presence of the third body physically changes the orbital period of the inner binary; (3) the shapes of the binary eclipses; and (4) when very lucky, eclipses of the third body by the binary and vice versa.
The Roemer delay has an amplitude of ARoem G 1/3 c(2π) 2/3 MC M 2/3 ABC P 2/3 2 (1 − e 2 2 cos ω2) 1/2 sin i2 (1) (see, e. g. Rappaport et al. 2013), where MC and MABC are the masses of the third body and the whole triple system, respectively, while P2, i2, e2, and ω2 are the orbital period, inclination angle, eccentricity, and argument of periastron of the outer orbit of the triple star system, respectively. By contrast, the amplitude of the dynamical delay is given approximately (1 − e 2 2 ) −3/2 e2 (2) (see Rappaport et al. 2013), for the case of near co-planarity of the orbital planes, where P1 is the orbital period of the inner binary. We illustrate the magnitudes of these delays in Fig. 1. To keep the plot simple, we fix all the masses at 1 M , the eccentricity e2 = 0.3, the orbit inclination angle i2 = 60 • , and ω2 = 45 • . The plot shows contours of constant delay amplitudes for both the Roemer and physical delays at 50, 100, 200, 500, and 1000 s. In addition, we show contours of constant probability for outer, or third-body, eclipses. The latter are based on assumed circular orbits and stars of 1 R .
In all, with the CoRoT, Kepler, K2, and now TESS missions, there are more than 200 compact triple (or quadruple) systems where ETV measurements have revealed the hierarchical nature of these systems (see, Borkovits et al. 2016;Hajdu et al. 2017 and references therein). In addition, in a very small subset of these cases, the orbital planes of the outer third star (or binary) are fortuitously aligned well enough with our line of sight so that there are so-called "third body" eclipsing events. These third body eclipses greatly enhance our ability to diagnose the system parameters. To our knowledge, there are only 17 such triply eclipsing systems 1 References: (1) Carter et al. (2011); (2) Masuda et al. (2015) with known outer periods (see Table 1). We show these in Fig. 1 as heavy filled circles. The outer orbits range from 34 days to 1100 days, while the binary periods cover an interval of 0.25 to 32 days. We also show in Fig. 1, as fainter points, the remaining ∼200 triple systems measured with Kepler that have known or estimated outer orbital periods (see Borkovits et al. 2016).
In the triply eclipsing systems, when eclipse timing variations are combined with a photodynamical analysis (see, e.g., Borkovits et al. 2019a), including spectral energy distributions ('SEDs'), stellar isochrone models, and the Gaia distance, many of the system parameters can be determined. This includes all three masses, periods, eccentricities, and inclination angles of the various orbital planes.
In this paper, we report on TIC 209409435, the first of the compact triply eclipsing triples found with TESS, plus a full orbital solution for the system. The system consists of a 5.7-day binary in orbit with a third star in a 122-day outer orbit. The photometric properties of the composite system are summarized in Table 2. In Section 2 we describes all the available observational data, as well as their preparation for the analysis. Then, Section 3 provides a full explanation of the steps of the joint physical and dynamical modeling of the light-and ETV curves, SED, parallax and stellar isochrones. In Section 4 we discuss the results from astrophysical and dynamical points of views. Finally, in Sect. 5 we summarize our findings and draw conclusions from our work. The results are shown in the plane of binary orbital period vs. the period of the outer third body (i.e., of the triple star system). The colored horizontal lines are contours of typical amplitudes for the Roemer delay, while the diagonal lines are contours of typical physical delays (Borkovits et al. 2015). The horizontal dashed lines are crude estimates of the eclipse probability for a third body orbiting the binary. The 16 filled black circles are known triply eclipsing triple systems (see Table 1), and the red square is for TIC 209409435, the system reported in this work. The fainter dots are some 200 other triple systems found with Kepler that have measured or estimated outer orbital periods (see Borkovits et al. 2016). The TESS lightcurve from Sector 3, in which it was first discovered, and Sector 4. The full frame images have 30-minute cadence. The four arrows mark the times of the anomalous third-body eclipses. The first two of these correspond to the inner binary eclipsing the third star, while the second set is the opposite situation where the 3rd star eclipses the inner binary.
Notes. (a) TIC-8 catalog (Stassun et al. 2018). (b) AAVSO Photometric All Sky Survey (APASS) DR9, (Henden et al. 2015), http://vizier.u-strasbg.fr/viz-bin/VizieR?-source= II/336/apass9. (c) Gaia DR2 (Gaia collaboration 2018). (d) 2MASS catalog (Skrutskie et al. 2006). (e) WISE point source catalog (Cutri et al. 2013 In addition to conventional computer searches through the TESS data for periodic signals, e.g., due to transiting planets and eclipsing binaries, a group of us has been visually surveying all the the full-frame image ('FFI') stars down to about TESS magnitude 13.5. The lightcurves that we use for surveying the data are from the MIT Quicklook pipeline ('QLP'; Huang et al. 2019). For each star, five lightcurves are extracted from apertures with different sizes ranging from 1.75 − 8 pixels 2 . The best aperture is chosen for stars in a particular magnitude bin based on the photometric precision after detrending. Fainter stars have relatively smaller photometric apertures. The 1.75-pixel apertures are almost never used because our method is difference-image-based, and we automatically deblend the source flux based on the target's T magnitude. This requires that we have an aperture that includes most of the light from the source, and the 1.75-pixel aperture is not really sufficient for that. TIC 209409435 was observed by the TESS spacecraft (Ricker et al. 2015) during Year 1 in Sectors 3 and 4. No 2 minute cadence data were available, thus, this triple-star system was observed serendipitously as part of the FFI coverage of those sectors. For the FFIs, the observational cadence is 30 minutes.  . Expanded view of the two sets of anomalous thirdbody eclipses. The dark blue circles represent those data points that were used for the photodynamical model (see Sect. 3), while the other out-of-eclipse data (not used in the modelling) are plotted as pale blue circles. The red curve is the photodynamical model solution corrected for the 30-minute integration time; the residuals to the model are also shown below the lightcurves. The top panel shows the first two events where the inner binary passes in front of the third star and produces two eclipses in addition to the regular binary eclipses spaced every 5.717 ÷ 2 days. The broader of the anomalous eclipses occurs in between two regular eclipses where the transverse motion of the binary stars across the sky is at a relative minimum. The bottom panel shows the two anomalous eclipses when the third star passes in front of the inner binary. By chance circumstances, the pattern of the two anomalous events in the top and bottom panels, with respect to the regular eclipses, is nearly identical.
While surveying these images using the LcTools software system (Schmitt et al. 2019), two of us (RG and TJ) independently discovered the triple nature of this source (on 19 and 24 November, 2019). The QLP TESS lightcurve is shown in Fig. 2. Aside from an unremarkable 5.7-day binary with two nearly equal eclipses per orbit, it is apparent that there are four anomalous eclipses (marked with arrows in the figure). The two sets of closely spaced anomalous eclipses are separated by ∼32 days. The anomalous eclipses are shown in more detail in Fig. 3. The red curve is a model fit which will be discussed later in the paper. Because of the temporal pattern of the anomalous eclipses, we immediately suspected that these were due to a third body orbiting the binary.
We then computed an eclipse timing variations ('ETV') . Eclipse timing variations of TIC 209409435 during the TESS discovery observations. The red and blue points are for the primary and secondary eclipses, respectively. The highly significant non-linear behavior is a clear indication that the third star is perturbing the binary. The smooth curves are fits to an ad hoc third-order polynomial just to guide the eye. The pale and dark blue points are the original and the 1-hr averaged WASP photometric data for this target, respectively. Only the dark blue points were used in the fit. The solid red curve is the photodynamical model back-projected more than a decade from the recent TESS observations. curve for the 17 available eclipses (see Fig. 2). The results are shown in Fig. 4, and tabulated in Table 3. (Note, due to a lack of data points on the ingress of the very first eclipse, it was dropped from the ETV curve used for the complex photodynamical analysis in Sect. 3, and not listed in the table.) The red and blue points represent the primary and secondary eclipses, respectively, while the smooth curves are fits to a cubic function just to guide the eye. The offset between the two curves by ∼20 minutes indicates that the binary orbit has a small eccentricity. The non-linear behavior of the ETVs, over just the 50-day interval of the TESS observations, confirms that there are dynamical interactions with the third body.

WASP Observations
Fortuitously for our study, TIC 209409435 was observed within the field of the WASP-South project Collier Cameron et al 2006) during three seasons between June 2006 and January 2012. (There was a huge gap in the observations between January 2008 and August 2011.) The WASP instruments each consist of an array of 8 cameras with Canon 200-mm f/1.8 lenses and 2k×2k e2V CCD detectors providing images with a field-of-view of 7.8 • × 7.8 • at an image scale of 13.7 arcsec/pixel. Images are obtained through a broad-band filter covering 400-700 nm. Fluxes are measured in an aperture with a radius of 48 arcsec for the WASP data. The data are processed with the SYSRem algorithm (Tamuz et al. 2005) to remove instrumental effects.
Almost two dozen eclipses can be identified in the WASP archival observations of TIC 209409435. Although many of these eclipses were observed only partially and, therefore, we were able to determine mid-eclipse times for only a portion of them (see in Table 3), these data were essential for our analysis. In particular, they helped to constrain: (i) the outer orbital period -not only through the ETV of the regular eclipses, but also by narrowing the possible locations of the outer eclipses; (ii) the 'flatness' of the triple system through the presence, or the lack of eclipse depth variations which might be forced by nodal precession in case of non-coplanarity of the inner and outer orbital planes; and also (iii) the dynamically forced apsidal motion, and in such a way, the dynamical properties of this strongly gravitationally interacting triple. Therefore, we used the entire WASP lightcurve for the complex photodynamical analysis (see Sect. 3). We plot two sections of the WASP data in Fig. 5.

Follow-up Ground-Based Observations
Following the discovery of this triply eclipsing triple star system we organized a photometric follow-up observational campaign with the participation of three amateur astronomers operating their own private observatories. We describe the observatories and the observations in the following subsections.

Ground-Based Observatories
PEST Observatory: The observatory is owned and operated by Thiam-Guan (TG) Tan. PEST is equipped with a 12-inch The early February event was similar to those that were observed by TESS, i. e., two separate outer eclipses (around times 8882.0 and 8883.5) which occurred between two regular eclipses, indicating that the two members of the inner binary eclipsed the outer component separately. By contrast the next outer eclipse event, between 8914 and 8916 displayed a long duration, large amplitude, irregular shaped dip, resembling the one that was observed by Kepler in the triple system EPIC 249432662 (Borkovits et al. 2019a). In this event the outer star eclipsed the binary members during their conjunction and, therefore, due to the almost coplanar configuration, the secondary of the inner binary had a small velocity relative to the outer star, resulting in the long-lasting eclipse, while the inner pair's primary, moving into the opposite direction, produced the sharp dip at the bottom of the main event.
Meade LX200 SCT f/10 telescope, an SBIG ST-8XME CCD camera, a BVRI filter wheel, a focal reducer yielding f/5, and an Optec TCF-Si focuser controlled by the observatory computer. PEST has a 31 × 21 field of view and a 1.2 per pixel scale. PEST is located in a suburb of the city of Perth, Western Australia. Observations of TIC 209409435 were carried out in Cousins RC-band, with an exposure time of 60 sec. A custom pipeline based on C-Munipack was used to calibrate the images and extract the differential time-series photometry using a photometric aperture of 6.2 radius.

Raemor Vista Observatory and Junk Bond Observatory:
The observations were conducted by T.G. Kaye and the image sets were calibrated and measured by B. Gary and T.G. Kaye. Raemor Vista has a 0.4-m Dreamscope with an Apogee CG 16 M CCD camera. 30-s exposures were used, and the images were spectrally unfiltered and binned into 2×2 pixels. The Junk Bond Observatory has a 0.8-m Ritchey Chrétien telescope with an SBIG STL6303E CCD camera. 30-s exposures were used, and the images were also unfiltered and binned into 2 × 2 pixels.
ROAD Observatory (Remote Observatory Atacama Desert): The observatory is located in San Pedro de Atacama, Chile and operated by F.-J. Hambsch remotely from Belgium. It uses a 16-inch telescope with a FLI ML16803 CCD camera having a 4096 × 4096 Kodak KAF-16803 image sensor. The exposure were 60 s with a Johnson V filter. Aperture photometry using the freely available software LESVEPHOTOMETRY was use together with reference and check stars. Dark and sky flat calibration frames were applied to the images.

Measured Ground-Based Eclipses
Observations were taken on 13 nights between 30 November 2019 and 7 March 2020. During the first six weeks of this campaign we were able to observe three regular eclipses of the inner binary (see Table 3), and an almost complete ingress portion of an additional secondary inner eclipse. These observations were essential for the preliminary photodynamical modelling (see Sect. 3) which made it possible to constrain the outer period of the triple and to predict the occurrence time of the forthcoming extra (third-body) eclipses with an accuracy of a few hours. Thanks to these predictions we were able to successfully observe and detect the upcoming third-body eclipse events which occurred relatively close to their predicted times (see Fig. 6). As it turned out, the first predicted third-body eclipse on 1 February 2020 was late in arriving by 6 hours, and this fortuitously gave us the longitude coverage on the Earth to detect it.

Overview of the Triple System
Based on the lightcurve information discussed in the previous subsections, we have carried out a detailed photodynamical evaluation of the system parameters. This is described in detail in Section 3, which follows next. For now we simply present a to-scale drawing of what the TIC 209409435 triple-star system looks like-see Fig.7. The red and blue lines in Fig. 7 represent the inner binary with its 5.717-d period, while the green curve is the orbit of the outer third body with a period of 121.87 days. The observer is viewing the system from a distant position along the −Y axis. The   outer orbit is sufficiently eccentric that, at least according to the drawing, it is easy to understand why the intervals between 3rd body events are not nearly equally spaced. We next describe how the photometry data alone, without any radial velocity measurements, are sufficient to deduce the constituent stellar masses and system geometry with remarkable fidelity.

JOINT PHYSICAL AND DYNAMICAL MODELLING OF ALL THE AVAILABLE OBSERVATIONAL DATA
We carried out complex photodynamical modelling of TIC 209409435 with the use of the software package Lightcurvefactory (see Borkovits et al. 2019aBorkovits et al. , 2020, and further references therein). In order to obtain a comprehensive, as well as physically and dynamically consistent model for this triple we simultaneously analysed: (i) four sets of photometric lightcurves (the two sectors of TESS measurements, the historical WASP observations, and T. G. Tan's Cousins RC-band and F.-J. Hambsch's Johnson V -band observations obtained during our photometric follow-up campaign); (ii) both the primary and secondary ETV curves deduced from these lightcurves (see Table 3); and also (iii) the archived stellar energy distribution (SED) for this system, in the form of different photometric passband magnitudes (see in Table 2).
The advantages of such a simultaneous, comprehensive analysis, as well as the consecutive steps of the complete procedure were explained in detail in a series of papers (Rappaport et al. 2017;Borkovits et al. 2018Borkovits et al. , 2019aBorkovits et al. ,b, 2020 and we believe that it is not necessary to repeat them here. Therefore, here we discuss only those details that are specific to the current study.
Before carrying out the photodynamical analysis we further prepared the four sets of lightcurves as described in the following three paragraphs.
For the photodynamical analysis of the TESS data we processed the original TESS full-frame images using a convolution-based differential photometric pipeline, based on the various tasks of the Fitsh package (Pál 2012). Then, from the raw lightcurve obtained from this process we removed what are likely instrumental effects using the software package Wōtan (Hippke et al. 2019). We also restricted the intervals of time within the lightcurves to be modelled to the orbital phase domain ±0.04 orbital cycles around each eclipse. The major exceptions to this were the regions around the third-body eclipses. For these, we also utilized the flat, out-of-eclipse lightcurves between the two consecutive regular eclipses preceding and following the two pairs of thirdbody eclipses (see in Fig. 3). Regarding the WASP measurements, we utilized data from the entire time domain because even the out-of-eclipse regions for the inner binary also carried potentially significant information about where the third body eclipses might have occurred. However, in order to reduce computational costs we formed 1-hr averages from these out-of-eclipse data points, and these binned data were then used for the analysis.
Similarly, in the case of the ground-based photometric follow-up observations, instead of the original data, we used their 15-minute averages.
In the first stage of our study we carried out a joint, simultaneous photodynamical analysis of the prepared TESS and WASP lightcurves, as well as the first few groundbased regular eclipse observations of T. G. Tan, together with the ETV curves derived from them. From this analysis, we were able to determine not only the outer period, together with the other orbital elements of the outer orbit, but also well constrained relative (i.e., dimensionless) stellar parameters (i.e., fractional radii and ratios of temperatures and masses). With these parameters in hand we were thereby able to predict the occurrence times of the forthcoming extra eclipses with sufficient accuracy to guide us in making further ground-based follow-up observations.
In the absence of radial velocity measurements, however, we were unable to determine model-independent, dynamical masses for each component. Therefore, in the second stage of the analysis, in order to obtain physical quantities within the framework of a self-consistent model, similar to the method followed in Borkovits et al. (2020), we added the observed composite SED of the triple system to the fit, and made efforts to find consistent, coeval PARSEC evolutionary tracks (Bressan et al. 2012) for the three stars.
For this latter purpose we used machine readable PAR-SEC isochrone tables generated via the web based tool CMD 3.3 3 . These tables contain theoretically computed fundamental stellar parameters and absolute passband magnitudes in several different photometric systems, for a large three-dimensional grid of ages, metallicities and initial stellar masses (see also Chen et al. 2019, and further references therein). The preselection of the most appropriate subsets of the tabulated isochrones for the current analysis, as well as the interpolation of fundamental stellar parameters and absolute passband magnitudes from the grid of values for an actual trial run is described in Section 3 of Borkovits et al. (2020).
The PARSEC tracks during each trial step were used as follows. The actual trial values of stellar mass, age, and metallicity determined the position of each star on the set of PARSEC tracks. Then, using a trilinear interpolation based on the closest grid points of the pre-calculated tables, the code interpolated the radius and temperature of each star as well as their absolute passband magnitudes for the SED fitting. These stellar radii and temperatures were used as input parameters to the light curve modelling section of the code. Furthermore, the interpolated absolute passband magnitudes transformed into model passband magnitudes with the use of the extinction parameter and the system's distance. Then, their sum was compared to the observed magnitudes in each passband. Note, however, in the last steps, the distance (d) was not a free parameter, but was constrained a posteriori in each trial step by minimizing the value of χ 2 SED . We intentionally did not use the Gaia DR2 parallax for constraining the distance. The reason is that the Gaia DR2 data might contain systematics for compact triple systems (see, e.g. the discussion of this problem in Borkovits et al. 2020). However, as we will discuss above, our results are in perfect agreement with the Gaia DR2 catalog data.
Therefore, in most of the runs we adjusted the following parameters: -Three parameters related to the orbital elements of the inner binary as follows: the eccentricity and the argument of periastron via e1 cos ω1 and e2 sin ω2, and the inclination, i1.
-Six parameters related to the orbital elements of the wide orbit of the third component: P2, e2 sin ω2, e2 cos ω2, i2, the time of the superior conjunction of the outer star, 3 http://stev.oapd.inaf.it/cgi-bin/cmd T sup 2 , and the position angle of the node of the wide orbit, Ω2. 4 -Three mass-related parameters: the mass of the primary of the inner binary, mA, and the mass ratios of the two orbits q1,2.
-Four lightcurve-dependent parameters as the passbanddependent 'third light(s)' T ESS , WASP, R C and V .
-Three parameters for the PARSEC isochrone and SED fitting: the age, log τ , and the metallicity, [M/H], of the system and the extinction coefficient E(B − V ). Furthermore, the following parameters were internally constrained: -The instantaneous orbital period, P1, of the inner binary and the inferior conjunction time, T inf 1 , of the secondary component of the inner binary, i.e., the mid-primaryeclipse-time at the zero epoch, were constrained via the use of the ETV curves in the manner explained in Appendix A of Borkovits et al. (2019a); -The stellar radii, RA,B,C, and effective temperatures, TA,B,C, were calculated from interpolation at each trial step from the appropriate grid elements of the PARSEC isochrone tables; -and, the distance of the system was constrained a posteriori by minimizing the value of χ 2 SED .
Regarding the other lightcurve-related parameters, we utilized a logarithmic limb-darkening law. The limbdarkening coefficients were interpolated from passbanddependent tables in the Phoebe software (Prša & Zwitter 2005). In turn, the Phoebe tables were derived from the Castelli & Kurucz (2004) stellar atmospheric models. We have found that due to the nearly spherical shapes of the stars in the inner binary, accurate settings of the gravity darkening coefficients have almost no effect on the lightcurve solution. Therefore, we simply adopted a fixed value of g = 0.32 which, according to the venerable model of Lucy (1967), is appropriate for stars with a convective envelope. Furthermore, we neglected the reradiation/illumination effect in order to save computation time. By contrast, the Doppler-boosting effect (Loeb & Gaudi 2003;van Kerkwijk et al. 2010), which is also found to be negligible for this system, requires only very minor additional computational time, and was therefore included in the model calculations.
The orbital and astrophysical parameters derived from the photodynamical analysis are tabulated in Table 4 and will be discussed in the subsequent Section 4. The corresponding model lightcurves are presented in Figs. 3, 5, and 6, while the model ETV curves plotted against the observed ETVs are shown in Fig. 8. Finally, the results of the stellar SED are plotted in Fig. 9.

DISCUSSION OF THE RESULTS
According to our results, TIC 209409435 is comprised of three very similar solar type stars. The two components of Table 4. Orbital and astrophysical parameters of TIC 209409435 from the joint photodynamical lightcurve, ETV, SED and PARSEC isochrone solution. Besides the usual observational system of reference related angular orbital elements (ω, i, Ω), their counterparts in the system's invariable plane related dynamical frame of reference are also given (ω dyn , i dyn , Ω dyn ). Moreover, im denotes the mutual inclination of the two orbital planes, while i inv and Ω inv give the position of the invariable plane with respect to the tangential plane of the sky (i. e., in the observational frame of reference Notes. a: Instantaneous, osculating orbital elements, calculated for epoch t 0 = 2458382.0000 (BJD); b: Interpolated from the PARSEC isochrones.  the inner binary pair are perfect twins, having a mass ratio of q1 = 1.002 ± 0.003, while the outer star has a larger mass only ∼9% higher than the inner components. These mass ratios are very robust and model independent due to the photodynamical part of our analysis. By contrast, in the absence of radial velocity data, the actual physical masses of the stars (and, correspondingly, other fundamental parameters such as the physical dimensions of the components), were obtained in an astrophysical model-dependent way, with the use of theoretical stellar evolutionary tracks in the form of PARSEC isochrones. Our combined analysis, however, has led to a self-consistent and robust result whose high fidelity is confirmed by the comparison of the observed net stellar SED to the model SED. This yields a system distance of d = 990±20 pc, which is about to 2-σ from the Gaia DR2 distance of dGaiaDR2 = 949±15 pc. Note, however, that the triple nature of this source might have resulted in systematic discrepancies in the Gaia DR2 parallax (see, e. g. Benedict et al. 2018).
From this same analysis we can state with high confidence that TIC 209490435 is a system having three MS stars with masses mA,B = 0.90 ± 0.03 M and mC = 0.98 ± 0.04 M , with effective temperatures TA,B = 5800 ± 100 K and TC = 6070 ± 90 K. The system is ∼ 5 +1 −2 Gyr-old, and slightly metal-deficient with [M/H] = −0.16 ± 0.1.
We have found that the triple is very flat, where the mutual inclination of the inner and outer orbits is im = 0. • 24 ± 0. • 08. The flatness of the system together with the almost circular inner, and moderately eccentric outer orbit, should make the configuration of this triple stable over the nuclear lifetime of the stars. We have taken a step toward verifying this with a 10 Myr-long numerical integration (carried out with the same integrator that was used for the photodynamical modelling). From an observational perspective, in the near complete absence of orbital plane precession, the inner pair will remain an eclipsing binary, and the outer star will also continue to produce the extra third-body eclipses over any human time scales. During the 10 Myr-long inte-gration, the inner and outer orbital inclinations oscillated between 89. • 64 i1 90. • 01, and 89. • 79 i2 89. • 86 with a period of Pprec ≈ 17 yr. 5 The apsidal precession of the low-eccentricity inner binary orbit has a very similar period of Papse1 ≈ 15 yr, while the outer orbit has a century-long Papse2 ≈ 103 yr apsidal period. The short-term variations in some of the orbital elements of the inner and outer orbits, obtained from our numerical integration, are plotted in Fig. 10.
Turning back to the most prominent characteristic of our triple, i.e. its flatness, other similarly very flat (im < 5 • ) and compact (P2 < 200 d) triple systems with accurately known parameters that were reported previously are HD 181068 (Borkovits et al. 2013), HD 144548 (Alonso et al. 2015), ξ Tau (Nemravová et al. 2016), EPIC 249432662 (Borkovits et al. 2019a), HIP 41431 (Borkovits et al. 2019b) and TIC 220397947 (Borkovits et al. 2020). 6 All of these systems show remarkable similarities, as one can see in Table 5. 5 One should keep in mind, however, that in the absence of any information on the axial rotation of the stars in the binary, we assumed synchronized and aligned rotation for the integration. Inclined spin axes might result in some orbital plane precession, and even drastic variations in the configuration of the triple system on a longer time scale (see, e.g. Correia et al. 2016). 6 Note, Borkovits et al. (2016) lists 8 additional compact triples with P 2 < 200 days in the original Kepler -field with mutual inclinations most probably less than 5 • . But, in the absence of accurate photodynamical solutions, we do not list these systems. Similarly, we omit λ Tau, the triple with the shortest outer period known (Ebbighausen & Struve 1956;Fekel & Tomkin 1981). Though, in the absence of any eclipse depth variations over the last 110 years one can suppose with a great certainty that this triple is also extremely flat (see, e. g. Söderhjelm 1975;Kiseleva et al. 1998;Berdyugin et al. 2018). Finally we note that we made a concerted effort to collect accurate photodynamical solutions from the literature for similarly flat systems with somewhat longer outer periods of up to P 2 ≤ 1000 days, but still relatively compact. Unfortunately, we were unsuccessful in finding any. Mean orbital phase of the wide orbit Figure 10. Variations of the instantaneous (osculating) anomalistic periods, eccentricities and (observable) arguments of periastron of the inner and outer binaries (grey curves in all panels), obtained via numerical integration, together with the one-inner/outer-orbit averages (black). Furthermore, the red and blue lines connect the instantaneous orbital elements sampled at those integration points that are closest to the primary and secondary mid-eclipse points of the inner binary, respectively. (Note, for the outer orbit the red and blue lines coincide almost perfectly, therefore, we plot only the red one.) For example, all systems have high mass-ratio inner binaries (q1 ≥ 0.88), and in five of the seven the inner pair is formed by near twin stars (i.e. q1 > 0.95). Furthermore, in all but one of the triples the outer third component is the more massive. Moreover, all inner orbits are almost circular, while the outer orbits, apart from HD 181068, have moderate eccentricities. 7 Such flat triple systems were most likely formed by disk fragmentation and accretion-driven migration. However, the quantitative details of these effects have not been fully explored (see, e. g. Tokovinin & Moe 2020), and it is also unclear what the role is of the third star in this process. Therefore, it is especially important to improve the sample of compact flat triple systems with accurately known system parameters.

SUMMARY AND CONCLUSIONS
In this work we report on the discovery and analysis of a new triply eclipsing triple star system, TIC 209409435. The discovery observations were made with TESS during Sectors 3 and 4. We also made use of archival WASP data as well as follow-up ground-based observations by amateur astronomers. We carried out a photodynamical analysis of the observational data which included the lightcurves, ETVs extracted from the lightcurves, the SEDs, and PARSEC evolution tracks. From this analysis we were able to obtain a comprehensive set of system masses and orbital parameters, all without being able to make RV observations.
In spite of the hundreds of thousands of eclipsing binaries that have been discovered over the years, including many whose light curves have been very well studied with CoRoT, Kepler, K2, and TESS, there are only 17 known EBs that 7 In this latter triple the outer component is a red giant, therefore, the circularization of the outer orbit can probably be explained with the more effective tidal damping effects in its present red giant state.
have eclipsing third bodies orbiting them, including this new discovery (see Table 1). TIC 209409435, among the 17 triply eclipsing triples, has the 6th shortest outer orbital period, the 4th highest ratio of P1/P2, and the 6th highest ratio of P 2 1 /P2, the latter being a measure of the dynamical time delays that are produced.
Thus TIC 209409435 is an impressive interactive system with pronounced ETVs and exotic-looking third-body eclipses twice every 122 days. The presence of the outer body eclipses helped enable us to make very robust determinations of the orbital parameters and the masses (see Table 4). The stellar masses and radii are good to about 3% accuracy on average. The triple system is extremely flat and aligned with the observers with i1 = 89.98 • ± 0.08 • , i2 = 89.79 • ± 0.09 • , and mutual inclination im = 0.24 • ± 0.08 • . The system is 5 Gyr old and is manifestly highly dynamically stable. The photometric distance of 990 ± 20 pc matches the Gaia distance of 949 ± 15 pc within 2-σ.
The photodynamical analysis has led to an orbital ephemeris that should be able to accurately predict when future third-body events can be observed from the ground, including with small telescopes. Moreover, TIC 209409435 is scheduled to be observed again in TESS Sectors 30 and 31 (i.e., nominally, between 22 Sept and 19 Nov 2020) during the extended Year 3 TESS mission. During this interval, extra third-body eclipsing events are expected on 1-3 October and 3-4 November. Though these dates are somewhat inauspicious in the sense that they are close to the mid-times of the two sectors (i.e. to the data downloading time), some portions of these events will hopefully be observed, allowing us to further refine the photodynamical model.
Overall we have found that TIC 209409435 has been a very gratifying system to have worked on.

ACKNOWLEDGMENTS
T. B. and T. M. acknowledge the financial support of the Hungarian National Research, Development and Innovation Office -NKFIH Grant KH-130372.