Four new compact triply eclipsing triples found with Gaia and TESS

ABSTRACT

This paper presents a comprehensive analysis of four triply eclipsing triple star systems, namely TIC 88206187, TIC 14839347, TIC 298714297, and TIC 66893949. The four systems with third-body eclipses were found in the TESS light curves from among a sample of ∼400 matches between known eclipsing binaries and the Gaia DR3 Non-Single Star solution data base. We combined photometric light curves, eclipse timing variations, archival spectral energy distributions, and theoretical evolution tracks in a robust photodynamical analysis to determine the orbital and system parameters. The triples have outer periods of 52.9 d, 85.5 d, 117 d, and 471 d, respectively. All dozen stars have masses ≲ 2.6 M_⊙. The systems are quite flat with mutual inclination angles between the inner and outer orbital planes that are all ≲ 4°. The outer mass ratios (q ≡ M₃/M_bin) range from 0.39 to 0.76, consistent with our earlier collection of compact triply eclipsing triples. TIC 88206187 exhibits a fractional radius of the outer tertiary component (r_B ≡ R_B/a_out) exceeding 0.1 (only the third such system known), and we consider its future evolution. Finally, we compare our photodynamical analysis results and the orbital parameters given in the Gaia DR3 NSS solutions, indicating decent agreement, but with the photodynamical results being more accurate.

1 INTRODUCTION

Triply eclipsing triple star systems are outstanding objects in which to study the properties and dynamics of multistellar systems. These are triple systems in which the inner binary eclipses the tertiary star, or vice versa, during the course of their outer orbit. The outer orbital periods range from a month to a year or two, and therefore dynamical interactions can usually be observed within the duration of a PhD study. Furthermore, even in the absence of radial velocity (RV) observations, readily available spaced-based photometric observations, such as with Kepler (Borucki et al. 2010) or TESS (Ricker et al. 2015) can yield precise eclipse timing measurements, as well as the photometry of the third-body eclipses. These can be combined with the archival spectral energy distribution (SED), and processed in a comprehensive spectro-photodynamical analysis to yield a full set of accurately determined stellar and orbital parameters for the system. In terms of the spatial configuration of the system, this analysis can provide such quantities as the ratio of orbital periods, the orbital eccentricities, and the mutual inclination of the two orbits. In turn, these quantities play a crucial role in the long-term evolution of the whole system. These are therefore excellent laboratories for testing theories regarding the formation, evolution, and final evolutionary states of multistar systems.

The first such system was discovered in Kepler data by Carter et al. (2011), and the second one followed very shortly thereafter (Derekas et al. 2011). Within a few years, more than a dozen publications had reported the discovery of 15 more similar objects, among them the first detections using other observing facilities than Kepler, including CoRoT (Hajdu et al. 2017) and OGLE (Hajdu et al. 2022). Since then, TESS observations have accelerated the search for these objects and have led the way in such studies. Borkovits et al. (2020b) reported the first identification and analysis of a triply eclipsing triple based on TESS data, and in less than 3 yr since then, 20 additional similar objects have been identified and analysed in detail (Mitnyan et al. 2020; Borkovits et al. 2022; Powell et al. 2022; Rappaport et al. 2022, 2023). That means, thanks to Kepler and then nowadays TESS, the discovery of triply eclipsing triples has become more or less routine, and the numbers of known triply eclipsing triple systems have been steadily growing.

In our previous work (Czavalinga et al. 2023), we utilized the Gaia Data Release 3 (DR3; Babusiaux et al. 2023; Gaia Collaboration et al. 2023) Non-Single Star (NSS; Gaia Collaboration 2022; Pourbaix et al. 2022) solutions to search for tertiary stars in eclipsing binary (EB) systems that had previously been reported in the literature. Gaia DR3 provides a unique opportunity to search for close hierarchical triple stars across the entire sky based on long-term astrometric and spectroscopic observations. If we compare the eclipsing period with the period obtained from the Gaia NSS solutions, and the period ratio is higher than 5, we can infer that the object is likely a triple-star system candidate with the inner period being that of the EB and the outer period of the tertiary star the one determined by Gaia’s NSS solutions (see Czavalinga et al. 2023, for details). Using this method, we identified 403 potential compact hierarchical triple-star system candidates among ∼1 million known EBs, including four newly identified triply eclipsing triples. Here, we provide a detailed analysis of the latter four systems. In Section 2, we describe the observational data and the methods utilized for the discovery and validation of the corresponding objects. In Section 3, we construct a photodynamical model for each system, and then in Section 4, we give a detailed discussion of the resulting models. In Section 5, we compare the orbital parameters from our photodynamical models with those from the Gaia NSS solutions, and finally we summarize our findings in Section 6.

2 DISCOVERY AND OBSERVATIONS

2.1 TESS observations

To validate the triple nature of our candidate systems we started by constructing a TESS photometric light curve (LC) for each one. Our aim was to search for eclipse timing variations (ETVs) and possible third-body eclipses in the LC. After downloading all the available Full-Frame Images from the MAST portal,¹ we applied the same convolution-aided image subtraction photometry pipeline based on FITSH (Pál 2012) along with the same detrending process based on the Wōtan (Hippke et al. 2019) and lightkurve (Lightkurve Collaboration et al. 2018) Python packages that we used in our previous paper (see details in Czavalinga et al. 2023).

We identified six LCs with potential extra third-body eclipses. Two of them had already been discovered by other teams, with one system being HD 181068, the second triply eclipsing triple found in Kepler observations (Derekas et al. 2011), and the other one being TIC 229785001 recently analysed in detail by Rappaport et al. (2023). The remaining four systems were completely unknown in the literature. These are as follows:

TIC 14839347 was observed by TESS in Sectors 14, 15, 41, and 55. We found a well-defined third-body eclipse in Sector 41, and then, after making use of the outer period from the Gaia NSS solution, we additionally found two very shallow secondary eclipses in Sectors 14 and 55 (see upper left panel of Fig. 1).

TIC 66893949 was observed in Sectors 15, 41, and 55. In Sector 15, we identified an extra dip in the LC, which we interpreted as third-body eclipse. No eclipses occurred in Sector 41, but we observed two extra eclipsing events with different depths in Sector 55 (upper right panel of Fig. 1). We note that recently Rowan et al. (2023) also independently discovered this object as a triply eclipsing triple, but they have not analysed the system in any detail.

TIC 88206187 was observed by TESS in Sectors 19 and 59. We discovered a very deep and long, flat bottomed eclipse in the middle of Sector 19. The source also produced another, similarly long in duration but much shallower, third-body eclipse in Sector 59 (lower left panel of Fig. 1).

Finally, TIC 298714297 was observed in Sectors 15, 55, and 56 by TESS. We found evidence for one third-body eclipse in Sector 55, accompanied by what are likely prominent, stellar-activity induced variations in the LC (lower right panel of Fig. 1).

For the ETV calculation, we applied the same process as in our previous paper (Czavalinga et al. 2023), and that is described in detail in Borkovits et al. (2015). In short, we fitted template polynomials to both the primary and secondary eclipses in the phase-folded, binned LC of each system. The mid-eclipse times were thereby determined for each orbital cycle. We list the determined individual mid-eclipse times in Tables A1–A4 separately for all the four systems.

2.2 ASAS-SN, ATLAS observations

We looked up all four triply eclipsing triples in the ASAS-SN (Shappee et al. 2014; Kochanek et al. 2017) and ATLAS (Tonry et al. 2018; Smith et al. 2020) archives. All four systems have good data in the ASAS-SN archives, and two have useful data in the ATLAS archives (TIC 298714297 and TIC 66893949 are too bright for ATLAS). We readily see the EB LC in the archival data for all four systems.

After we determine an accurate EB period for data that span of order a decade, we remove the EB LC by fitting for, and then subtracting out, between 50 and 100 orbital harmonics. This allows for a more sensitive BLS search (Box Least Squares; Kovács, Zucker & Mazeh 2002) for the eclipses of the outer orbit.

For TIC 66893949 and TIC 298714297 (with outer periods of 470.6 d and 118.6 d, respectively), there are insufficient archival data, and too few outer eclipses observed, to make a statistically significant detection of the outer eclipses. However, for TIC 14839347 and TIC 88206187, the outer periods of 85.47 d and 52.84 d, respectively, were robustly detected in a combination of the ASAS-SN and ATLAS data. The results are shown in Figs 2 and 3, respectively. In the case of TIC 88206187, Fig. 3 shows a shallow but clear secondary outer eclipse. The phasing of the two outer eclipses yields e_outcos ω_out = 0.003 ± 0.005, with an implication that the outer orbit is likely fairly circular.

3 PHOTODYNAMICAL ANALYSIS

Similar to other TESS-discovered triply eclipsing triple stars, we carried out a joint, simultaneous LC, ETV curve and SED analysis with the software package Lightcurvefactory. The software package itself, and the consecutive steps of the complex analysis, were described in several papers, e.g. in Borkovits et al. (2020a, 2022) and Rappaport et al. (2023), and hence, we will not repeat the details here, but rather restrict ourselves to some particular notes about the current, specific systems.

Our analyses are mainly based on the TESS LCs, which were processed in the manner described in Section 2.1. In contrast to the vast majority of our former analyses of triple and quadruple systems, in the case of three of the four systems reported here, we did not restrict our analyses only to the narrow regions of the inner and outer eclipses, dropping out the majority of the out-of-eclipse portions of the LCs, but rather we kept the complete TESS time series. The reason is that, apart from the somewhat wider 4.8-d inner EB in TIC 66893949, the other three, much more compact inner EBs, display significant ellipsoidal light variations in the out-of-eclipse LC sections, which we found to be worthwhile to retain for the analysis. Moreover, we note the special case of TIC 298714297, in which the LC exhibits further remarkable rotational variations with periods similar to the eclipsing period of the inner EB. We modelled these variations simultaneously with the full, eclipsing LC solution (but separately for Sectors 15 and 55–56), with the use of additional fits of harmonic functions in the manner described, e.g. in Borkovits et al. (2018, 2021).

As usual, we fit the ETV curves of the four EBs simultaneously with their LCs. As our literature search did not turn up any additional previously recorded eclipse times, our ETV curves were restricted to only those times that we determined from the TESS observations themselves.

Due to the lack of publicly available RV timeseries for any of the four currently investigated systems, with the exception of TIC 14839347, we used a combination of (i) observed composite SED values (tabulated from the available catalogue passband magnitudes, and listed in Table 1) and (ii) pre-computed, tabulated PARSEC isochrones (Bressan et al. 2012) to find stellar masses and corresponding effective temperatures, as well. The use of PARSEC isochrone-based SED analysis for this purpose was explained in Borkovits et al. (2020a), while a comparison of the accuracy and efficiency of such an astrophysical model-dependent analysis with the classical, astrophysical model-independent analysis (the latter of which is based on RV data), was carried out and discussed in Borkovits et al. (2022).

In the case of one of the four systems, TIC 14839347, our preliminary analysis run, however, revealed a situation where the secondary is fairly low in mass, but large enough to fill or nearly fill its Roche lobe. We therefore suspect that the secondary star in TIC 14839347 may have already transferred a portion of its envelope to the (current) primary star. As a consequence, for this triple we cannot use this proxy method because the precomputed PARSEC evolution grids are valid only for single stars, i.e. for such binary components, which have not previously exchanged mass. Therefore, in the case of this particular triple, we followed an iterative method as follows. First, we carried out Lightcurvefactory fitting, excluding the SED analysis section. In such a way, we obtained strong constraints for the relative (or, fractional) radii of the three stars, as well as their temperature ratios. Then, using these relative quantities, as well as the Gaia EDR3 derived distance (Bailer-Jones et al. 2021), we made an independent SED analysis without the use of any astrophysical pre-assumptions, with a slightly modified version of the method described in Rappaport et al. (2022), Section 3. In this way, we obtained likely mass and temperature ranges for the three constituent stars. Then we used the SED-determined values for the mass of the primary of the inner binary (m_Aa) and the effective temperature of the tertiary (T_B) with their statistical uncertainties as Gaussian priors, and reiterated the whole simultaneous LC and ETV curve-fitting procedure with Lightcurvefactory. With this iterative procedure, we were able to find the system geometry and dynamics robustly, while also inferring the physical properties of the constituent stars at a reasonable level of accuracy.

In regard to TIC 14839347, we also had to introduce some further modification to the light-curve fitting part of the complex analysis. We found that, in this case, the reflection/irradiation effect gives a significant contribution to the LC of the inner EB. Hence, we ‘switched on’ this effect for the modelling. First we set the bolometric albedos according to the theoretically expected values of A = 1.0 and 0.5 for the primary (Aa) and secondary stars (Ab), (i.e. radiative and convective envelopes), respectively. We found, however, that the fit becomes distinctly better when one sets A_Ab ≈ 0.8. For this reason, finally we decided to adjust, as an exception to our usual approach, the bolometric albedos and gravity darkening exponents (β) for these two stars.

Finally, we note the fact that none of the four investigated triples is a ‘tight’ system (the outer-to-inner period ratios, P_out/P_in are somewhat large, ranging from 45 for TIC 88206187, to 109 TIC 298714297) and, hence, one cannot expect significant third-body perturbations which would produce significant departures from simple Keplerian motions. On the other hand, considering the fact that the third-body eclipses are extremely sensitive to the current system geometry, similar to our former analyses (where much tighter triples were analysed) we integrated the motion of the three stellar components numerically, instead of approximating the motion with the superposition of two Keplerian orbits. The numerical integrator built into Lightcurvefactory is a seventh order Runge–Kutta–Nystrom integrator. This integrator provides the gravitationally (pure Newtonian two- and three-body terms), tidally (within the framework of the equilibrium tide model) and, optionally, the relativistically perturbed Cartesian Jacobian coordinates (and velocities) for each observational instant. A detailed description of the integrator itself, and its implementation in the photodynamical software package can be found in Borkovits, Forgács-Dajka & Regály (2004) and Borkovits et al. (2019a, b).

We initialized and ran several MCMC chains for each of the four triple systems. The best-fitting parameters (median values of the posteriors), their 1−σ uncertainties, as well as several derived parameters, are tabulated in Tables 2 and 3. Here we note, that while most of the given derived parameters (e.g. the semimajor axes, RV amplitudes, bolometric lumonisities, and so on) do not require any further explanations, we discuss briefly the calculations and significances of the different apsidal and nodal motion parameters, tabulated in between the orbital elements and the stellar parameters.

We give the theoretical apsidal motion periods both in the observational and the dynamical frames of reference (P_apse and |$P_\mathrm{apse}^\mathrm{dyn}$|⁠). While the formulae for their calculations, from the initial osculating orbital elements at each accepted MCMC trial step, are discussed in detail in section 6.3 of Kostov et al. (2021), here we emphasize only the fact that the observational and dynamical apsidal motions and, hence, their time-scales are substantially different. From an observational point of view, what is significant is the revolution of the argument of periastron in the observational frame of reference.² This revolution manifests itself, e.g. in the periodic, quasi-sinusoidal and anticorrelated shifts of the primary and secondary eclipses, in the case of eccentric EBs. In contrast to this, the revolution of the argument of pericentre in the dynamical frame³ cannot be observed directly, however, this dynamical argument of periastron and, hence, the dynamical-frame apsidal motion are significant for dynamical studies, as this parameter is one that occurs in the perturbation equations. The observable-frame apsidal motion, in general, is a non-linear combination of the dynamical-frame apsidal motion and the (dynamical) nodal regression, and hence, their periods (i.e. P_apse and |$P_\mathrm{apse}^\mathrm{dyn}$|⁠) may be substantially different (see e.g. Borkovits et al. 2015 for further details). Note also, that besides the apsidal motion and nodal regression (⁠|$P_\mathrm{node}^\mathrm{dyn}$|⁠) periods, we give separately the three components of the apsidal advance rates (classical tidal, Δω_tide; general relativistic, Δω_GR; and dynamical third-body, Δω_3b, respectively) for one orbital revolution of the inner and outer orbits, respectively. These contributions are calculated in the dynamical frame of reference.

Finally, note that the third-body eclipse sections of the model LCs from the best-fitting complex photodynamical solutions are plotted in Fig. 1, while the model ETV curves of the best-fitting solutions are shown in Figs 4–7. We briefly discuss our results for each system separately in Section 4.

4 DISCUSSION

4.1 TIC 14839347

As was mentioned above, the analysis of this triple system had to be carried out iteratively, and with extra care. We adopted |$m_\mathrm{Aa}=2.50\pm 0.13\, \mathrm{M}_\odot$| and T_B = 5000 ± 200 K for the mass of the primary and effective temperature of the tertiary, from the separate SED analysis. The inner mass ratio was found to be as low as q_in = 0.25 ± 0.02, which, together with a near Roche lobe-filling secondary, make this inner binary appear to be a typical Algol-type system with a reversed mass ratio. The distant third star was found to be a little less massive than the current primary of the inner binary (⁠|$m_\mathrm{B}=2.36_{-0.09}^{+0.17}\, \mathrm{M}_\odot$|⁠).

This finding again emphasizes that, most likely, intensive mass exchange has occurred between the two inner binary stars in the near past of this triple system, as the tertiary component is evidently an evolved, red giant star (⁠|$R_\mathrm{B}=8.2\pm 0.3\, \mathrm{R}_\odot$|⁠), in contrast to the similarly massive, hot (T_Aa = 8975 ± 380 K) and much less evolved (⁠|$R_\mathrm{Aa}=2.87\pm 0.03\, \mathrm{R}_\odot$|⁠) primary component.

Here we emphasize again that, in the case of this particular system, we took into account the reflection/irradiation effect and, moreover, both the two gravity-darkening exponents and the bolometric albedos of the inner pair of stars were freely adjusted MCMC parameters. In the case of the hot and, hence, radiative primary star (component Aa) we obtained an albedo of A_Aa = 0.97 ± 0.05, which fits nicely with the theoretically expected value of |$A_\mathrm{Aa}^\mathrm{theo}=1$|⁠.⁴ On the other hand, the bolometric albedo obtained for the convective secondary (A_Ab = 0.76 ± 0.05) significantly differs from the theoretically expected value of |$A_\mathrm{Ab}^\mathrm{theo}=0.5$|⁠. According to our knowledge, such a high albedo for a convective star is quite unusual, however, a deeper investigation of this question is beyond the scope of this paper.

Turning to the orbital properties of this triple system, we find that the inner orbit is circular, as is expected in the case of a nearly semidetached system.⁵ By contrast, the outer orbit displays some small, but significant, eccentricity (e_out = 0.04 ± 0.01). What makes this result, however, a bit less robust is that the secondary outer eclipses (in Sectors 14 and 55) are located exactly at phase 0.5. This is reflected by the fact that the outer argument of pericentre is ω_out = 270° ± 3° and, hence, e_outcos ω_out ≈ 0. In turn, this would suggest that the outer orbit is seen exactly from the direction of its major axis. However, one should keep in mind that, in this particular case, the outer eccentricity is determined only via the weakly determined parameter e_outsin ω_out, which primarily manifests itself in the durations of the outer eclipses and has strong correlations with the outer inclination, the fractional radius of the tertiary, the node of the outer orbit relative to the inner one, and even some other parameters.⁶

As both the inner and outer orbits are fairly circular and almost coplanar (⁠|$i_\mathrm{m}=3{_{.}^{\circ}} 5_{-1}^{+3}$|⁠), one cannot expect large perturbations and, hence, there are only small departures from Keplerian motions for the inner and outer orbits. This is true despite the fact that the nominal apsidal motion period of the inner orbit is shorter than a year (P_apse = 0.80 ± 0.02 yr). As one can see from the different kinds of apsidal advance rates Δω, in the case of the inner orbit the classic tidal (Δω_tide) contribution is the dominant one – larger by one order of magnitude than that of the third-body-driven apsidal motion. We emphasize again, however, that as the inner orbit is practically circular, this very short apsidal motion is only virtual, stemming from the above mentioned slight imperfections in the use of osculating orbital elements in this particular case. On the other hand, for the small but clearly non-zero mutual inclination, one can expect a very low (few degrees) amplitude inclination variation with a period of |$P_\mathrm{node}^\mathrm{dyn}=51\pm 1$| yr, which might manifest itself in very small eclipse depth variations in the future.

Here, we note also that the SED portion of the iterative solution for this system allowed for an age determination of the system of 520 ± 115 Myr,⁷ and an independent fit for the interstellar extinction of A_V = 2.2 ± 0.16.

4.2 TIC 66893949

This system has the longest inner (P_in = 4.805 d) and outer (P_out = 471.0 d) periods. Because the TESS observations cover only 16–17  per cent of the full outer orbit (see Fig. 5) it is quite surprising (and, of course, fortuitous) that both primary and secondary third-body eclipses were observed.

According to our photodynamical results, the system consists of three sun-like, main sequence stars (⁠|$m_\mathrm{Aa}=1.51\pm 0.04\, \mathrm{M}_\odot$|⁠, |$m_\mathrm{Ab}=1.01\pm 0.02\, \mathrm{M}_\odot$|⁠, and |$m_\mathrm{B}=0.97\pm 0.04\, \mathrm{M}_\odot$|⁠). Note, amongst our four triple systems, this is the only one where the tertiary star was found to be the least massive.

In connection with the somewhat larger inner orbital separation, the inner binary pair has a small, but significant, eccentricity (e_in = 0.005 ± 0.001). The outer orbit also has the largest eccentricity in our sample with e_out = 0.402 ± 0.004. Interestingly, despite the somewhat wide configuration, this triple was found to be quite flat (i_mut = 0|${_{.}^{\circ}}$|7 ± 0|${_{.}^{\circ}}$|5).

The dynamical time-scales in this systems exceed several centuries for both the inner and outer subsystems. Moreover, as one can see from the separate apsidal advance rates, both the classic tidal and the relativistic contributions to the apsidal motions are negligible relative to the dynamical (forced by third-body perturbations) apsidal motions of the inner and outer orbits.

Finally, we note that our analysis resulted in a slightly larger distance than that of Bailer-Jones et al. (2021), i.e. d_phot = 625 ± 17 pc, versus d_EDR3 = 587 ± 9 pc, a 2−σ discrepancy.

4.3 TIC 88206187

Despite the fact that this triple was observed only in two TESS sectors, we were able to obtain the most robust solution for this system among those studied in this work. The reason is that, fortunately, TESS observed both primary and secondary third-body eclipses and, moreover, both the inner and outer eclipses are total (i.e. having flat-bottom mid-eclipse sections), which strongly constrain the surface brightness ratios of the constituent stars. Moreover, the combination of the TESS observations with the archive ground-base data (see Section 2.2) we were able to determine the outer period with considerable accuracy (P_out = 52.84 ± 0.01 d), which was found to be in reasonable agreement with that of the somewhat less accurate Gaia SB1 solution (P_out = 53.03 ± 0.06 d; Table 4). This makes this triple the most compact in our sample.

Already, a very first inspection of the 2-d-long primary third-body eclipse around BJD 2 458 829 (at the middle of Sector 19, see lower left panel of Fig. 1) reveals that the tertiary is likely the largest component and, most probably, is a red giant star.⁸ The unequal primary and secondary third-body eclipse depths, with the flat bottom of the much deeper (primary) third-body eclipse, indicate that the members of the inner EB are much smaller, but substantially hotter, than the tertiary star. From this, one can directly infer that the tertiary component should be the most massive amongst the three stars, being the most evolved (assuming, of course, coeval evolution for the three objects). Finally, the fact that, in the case of the primary third-body eclipse (at Sector 19) the tertiary component passes its inferior conjunction point, reveals the outer orbital phase at that (or any other) moment uniquely, despite the absence of any informative ETV data (see Fig. 6).

The photodynamical solution then confirmed all these preliminary assessments, as it was found that the tertiary star was actually a moderately massive (⁠|$m_\mathrm{B}=2.6\pm 0.1\, \mathrm{M}_\odot$|⁠) red giant star with the basic parameters of |$R_\mathrm{B}=11.7\pm 0.3\, \mathrm{R}_\odot$| and T_B = 4950 ± 70 K. The members of the inner binary pair, however, were found to be A and F-type, but still MS stars, with parameters of |$m_\mathrm{Aa,Ab}=2.05\pm 0.05\, \mathrm{M}_\odot$|⁠, |$1.35\pm 0.03\, \mathrm{M}_\odot$|⁠; |$R_\mathrm{Aa,Ab}=2.25\pm 0.05\, \mathrm{R}_\odot$|⁠, |$1.35\pm 0.03\, \mathrm{R}_\odot$|⁠; |$T_\mathrm{Aa,Ab}=8\, 360^{+270}_{-160}$| K, |$6\, 590^{+170}_{-70}$| K for the primary and secondary binary components, respectively.

The future evolution of this triple seems to be quite interesting. Due to the compactness of both the inner and outer subsystems, we may expect that both the currently more massive red giant tertiary and the TAMS-aged primary of the inner binary will fill out their respective Roche lobes. The current size of the tertiary and its Roche lobe radius are 11.7 R_⊙ and 38.3 R_⊙, respectively. The same quantities for the primary EB star, Aa, are 2.25 R_⊙ and 2.94 R_⊙, respectively. Thus, star Aa needs to increase its radius only by |$\simeq 30~{{\ \rm per\ cent}}$| in order to overflow its Roche lobe, whereas the tertiary would need to triple its current radius to overflow its Roche lobe. None the less, according to MIST stellar evolution tracks (Choi et al. 2016; Dotter 2016) it will take star Aa some 260 Myr to reach a state where it overflows its Roche lobe with respect to star Ab. By contrast, the more massive and currently more evolved tertiary star will evolve much more rapidly and will reach the tip of the red giant branch (RGB) in a few Myr. At that time, its radius will not quite be sufficiently large to fill its Roche lobe. However, in less than 180 Myr after that, it will ascend the asymptotic giant branch (AGB) and overfill its Roche lobe for certain.

The detailed evolution of this system, after mass transfer from the giant tertiary to the inner binary commences, is beyond the scope of this paper. However, we can speculate that since the tertiary will be a convective giant at that time, and the outer mass ratio (m_B/(m_Aa + m_Ab)) is 0.76, the mass transfer could be dynamically stable. But, formally, for a completely convective donor star, this ratio should be ≲ 2/3.

Apart from its future evolution with mass transfer, this triple may harbour other interesting effects in its present configuration as well. Recently, Gao et al. (2023) called attention to the fact that such compact triple stars, where the tertiaries are red giants, may produce remarkable tidal effects whereby tidal dissipation could lead to observable orbital shrinkage within some decades-long time-scales. In this regard, we note that amongst all the known analysed compact triple star systems, TIC 88206187 is only the third case in which the fractional radius of the distant tertiary exceeds 0.1 (r_B = R_B/a_out = 0.109 ± 0.002),⁹ and thus, tertiary tides might currently be effective, at least marginally. In this regard, also note that, despite the relatively more significant apsidal advance rate of the outer orbit due to tidal forcing (Δω_tide = 9″ ± 1″/cycle), it is clearly negligible relative to the third-body effects (Δω_tide = 11400″ ± 100″/cycle) and does not alter the time-scale of the apsidal precession (P_apse = 185 ± 1 yr).

Finally, we note that the photometric distance inferred from our analysis is d_phot = 2580 ± 60 pc, which is in quite a good agreement with the Gaia-parallax inferred distance of d_EDR3 = 2490 ± 90 pc (Bailer-Jones et al. 2021).

4.4 TIC 298714297

In contrast to the other systems in this study, TIC 298714297 is found to be a triplet of three low-mass cool red (K and M-type) dwarfs. Besides the three dips of a complex third-body eclipse, the other remarkable feature of the LC is the likely rotational distortions, caused by stellar spots, which are most pronounced in the Sector 55, 56 LCs. Moreover, the system shows three sudden, short ∼|$3-4~{{\ \rm per\ cent}}$| brightenings, which we attribute to stellar flares.

All these features suggest strong chromospheric/photospheric stellar activity on at least one of the stellar components. In the case of the rotational modulation, since the observed period is equal, or very close to, the period of the inner EB, it seems quite likely that its origin is component Aa (which is the most massive and brightest amongst the three stars). In the case of the three stellar flare events, their origin, strictly speaking, is less certain. But, considering the other signals of strong magnetic activity are likely from star Aa, we may tentatively assume that the flare events are also hosted by this star. It is noteworthy, however, that the two larger flares were observed during Sector 15, when the rotational modulation (i.e. due to the spottedness) was much less pronounced.

Regarding the astrophysical parameters, this triple was found to be not only the least massive (⁠|$m_\mathrm{Aa}=0.83\pm 0.03\, \mathrm{M}_\odot$|⁠, |$m_\mathrm{Ab}=0.50\pm 0.02\, \mathrm{M}_\odot$|⁠, and |$m_\mathrm{B}=0.68\pm 0.02\, \mathrm{M}_\odot$|⁠) of our set, but also the most aged (τ ∼ 10 Gyr). Note, however, that some caution is necessary, as our model solution resulted in a reddening of E(B−V) = 0.3 ± 0.08, which looks quite unrealistic for such a close object (d_phot = 138 ± 4 pc versus d_EDR3 = 111.9 ± 0.7 pc; Bailer-Jones et al. 2021).

The inner orbit was found to be practically circular (⁠|$e_\mathrm{in}=0.002_{-0.001}^{+0.007}$|⁠), which is naturally expected for such an old and close binary, while the outer orbit was found to be moderately eccentric (e_out = 0.24 ± 0.03). In this regard, we note that the secondary’s ETV curve in the inner EB is offset a little from the primary’s one (Fig. 7). However, in our interpretation, in this case such an offset is not due a slight orbital eccentricity, but rather is caused by the LC variations due to the starspots on at least one of the stars. Such shifts in the ETVs were found by Tran et al. (2013).

Finally, note also that besides TIC 14839347, this is the other system in our sample where the apsidal motion of the inner orbit is dominated by the tidal forces instead of third-body perturbations (Δω_tide = 105″ ± 10″/cycle versus Δω_3b = 63″ ± 2″/cycle). But, again, due to the very low inner eccentricity one cannot expect strong, non-Keplerian variations during the apsidal motion cycle of P_apse = 28 ± 2 yr.

5 COMPARING ORBITAL PARAMETERS WITH GAIA

In general, we find substantial agreement between the Gaia and photodynamical results for the four systems studied in this work. As far as we can judge, however, our results are generally the more accurate ones. Having said this, however, there are a number of caveats to discuss. First, we note that a comparison of the orbital elements from the Gaia NSS solution and those from the Lightcurvefactory photodynamical models is limited due to the differences in their methodologies. The Gaia NSS solutions adopt a Keplerian two-body orbit model, whereas the calculated orbital elements from Lightcurvefactory represent instantaneous osculating elements. Consequently, the orbital elements obtained from these two approaches are not exactly comparable. However, for the four sources studied in this work, the perturbations to the motions along the outer orbits are negligible on the time-scale of the Gaia observation interval, and therefore they do not result in significant discrepancies on this account. For example, the apsidal motion time-scales of the outer orbits in all our four triples range from 185 to 2820 yr. Hence, the orientations of the astronomical orbits measured by Gaia will vary by only a negligible amount during the full operation of the space telescope. Moreover, due to the flatness of the four systems, no significant variations in the outer inclinations will occur. One can therefore expect that the orbits will remain (almost) pure Keplerian over the entire life cycle of Gaia.

Another issue is that in the Gaia DR3 catalogue, the NSS solutions are limited to spatially unresolved stellar systems, primarily due to the specific data processing constraints employed by the Gaia team. In the case of astrometric solutions (i.e. for TIC 298714297 and TIC 66893949) Gaia measures only the orbit of the photocentre of the three stars, and not of any of the individual constituent members within the triple system (Marcussen & Albrecht 2023; equation C2 of Rappaport et al. 2022). This fact, however, may produce some departures from pure Keplerian astrometric orbits, as during the relatively deep primary eclipses of the four inner pairs, one may expect some quasi-discontinuous shifts in the positions of the photocenter.¹⁰ In order to illustrate these effects, we plot the theoretically expected photocentre orbits (together with the above mentioned potential discontinuities in the small, inserted plots) for all four investigated systems in Fig. 8. Keeping in mind all of the above discussion, we conclude that the only direct comparisons that can be made for these systems with the photodynamical results are the period, the eccentricity, and the inclination of the outer orbit. In the cases involving the Gaia spectroscopic orbits (i.e. TIC 88206187 and TIC 14839347), Gaia primarily observes the spectral lines of the more luminous tertiary component. Thus, we can also compare the Gaia RV semi-amplitude (K1) of the outer tertriary with the value from the photodynamical fit.

A third issue is that generally the argument of periastron of the outer orbit from the Gaia measurements and those of the photodyamical results typically differ by 180°. This is due to the fact that Gaia is generally tracking the brighter tertiary, while the photodynamical analysis utilizes the ETV points from the inner binary.

Finally, with these caveats in mind, we compare the Gaia and photodynamical analysis in Table 4. We find generally good agreement for the outer orbital period and its inclination angle. However, in the case of the outer eccentricity, the differences are mostly larger, but in the case of TIC 66893949, for example, the two values are in good agreement.

6 SUMMARY AND CONCLUSIONS

The most common method for discovering triply eclipsing triple star systems involves searching for evidence of a third-body eclipse in the extended time series of EB systems. These utilize the extensive time series data sets obtained from the TESS or Kepler photometric sky survey telescopes. The searches are usually done by eye or with the assistance of machine learning tools (see e.g. Rappaport et al. 2022, 2023). It goes without saying that the longer the duration of a data set the greater the likelihood of identifying eclipses involving third bodies within a given system.

In our previous paper (Czavalinga et al. 2023), we reported the discovery of four triply eclipsing triple systems using an alternative method that utilizes the Gaia DR3 data set. In that study, we cross-matched a catalogue of approximately 1 million binaries with Gaia stars that had longer-period NSS orbital solutions. Of the ∼400 triple systems found in that study, four of them turned out to also be triply eclipsing, as discovered through the use of TESS data. Though the fractional yield was small (1  per cent), this still demonstrated the effectiveness of alternative methods in the search for triply eclipsing triples.

In this work, we studied these four triply eclipsing systems in detail using precision ETV curves and photometry from TESS data, archival photometry from ASAS-SN and ATLAS, and archival SED curves. These input data sets undergo a comprehensive analysis via a spectro-photodynamics code (LightcurveFactory) to yield both the stellar and orbital parameters of the systems. In general, the four new compact triples are not dissimilar to ones our group has been studying over the past number of years.

TIC 66893949 consists of stars with masses similar to our Sun. This system has the longest outer period among the triples investigated in our sample, and it has the most eccentric outer orbit e_out ≃ 0.4. Even though it has a wide configuration, the mutual inclination angle is fairly low at i_mut = 0.6 ± 0.5 degrees.

TIC 88206187 has the shortest outer period in our sample P_out = 52.92 ± 0.04 d, and contains a tertiary red giant star with radius of 11.7 ± 0.3 R_⊙. This represents only the third system discovered thus far in which the fractional radius (i.e. R/a) of the distant tertiary component exceeds 0.1, suggesting the presence of ongoing tidal effects. Our photodynamical model benefitted from the ATLAS and ASAS-SN LCs, which show both the outer primary and secondary eclipses, and therefore provide a precise period and a value for e_outcos ω_out. From the MIST stellar evolution tracks, we find that it is likely that the outer tertriary will be the first star in the system to fill its Roche lobe when it becomes an AGB star and transfer mass to the inner binary pair.

TIC 298714297 is one of the closest triply eclipsing triple systems ever found with a distance of 138 ± 4 pc. It contains low-mass stars that exhibit pronounced spot activity in one of the components of the inner binary, with several flares identified. The most probable host of these flares is star Aa with the largest luminosity in the system.

We encountered some difficulties with modelling TIC 14839347. The low mass secondary component of the inner pair was found to be strongly non-spherical, filling or nearly filling its Roche lobe, which hints at a previous episode of mass transfer. Consequently, the use of the PARSEC single-star evolutionary tracks within the Lightcurvefactory analysis could not be used. However, we were able to use a separate SED fitting code (Rappaport et al. 2022) that does not assume any physical relation between the radius and T_eff of either binary star, and which can estimate the stellar properties of all three stars. That code, in turn, used estimates of a number of dimensionless ratios such as of the radii and T_eff, which we could get from Lightcurvefactory. Since neither code could proceed without results from the other, we ran them iteratively until a satisfactory solution was achieved. From this analysis we find that the inner binary has a mass ratio of 0.25 ± 0.04 and, as mentioned above, the secondary is nearly filling its Roche lobe. Therefore the inner binary appears as an Algol system with an inverted mass ratio. The outer stellar component is an evolved red giant star, while the primary star of the inner binary appears to be less evolved than the tertiary though these two stars are similar in mass. This feature may also be explained by the fact that a mass exchange event has occurred between the inner binary stars in the near past.

The use of Gaia DR3 NSS solutions have demonstrated the value of exploring triple systems using alternative methodologies. The upcoming release of Gaia DR4 will provide even more extensive data, making it extremely worthwhile to search for signs of triples within these data sets and subsequently perform photodynamical analyses on them. Such a greatly enhanced statistical sample of compact triples should yield a more comprehensive understanding of these important, intriguing, and fun systems.

ACKNOWLEDGEMENTS

AP acknowledges the financial support of the Hungarian National Research, Development and Innovation Office – NKFIH grant K-138962.

This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts.

This work has made use of data from the European Space Agency (ESA) mission Gaia,¹¹ processed by the Gaia Data Processing and Analysis Consortium (DPAC).¹² Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

We have made use of the All-Sky Automated Survey for Supernovae archival photometric data. See Shappee et al. (2014) and Kochanek et al. (2017) for details of the ASAS-SN survey. We also acknowledge use of the photometric archival data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. See Tonry et al. (2018) and Heinze et al. (2018) for specifics of the ATLAS survey.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

We used the Simbad service operated by the Centre des Données Stellaires (Strasbourg, France) and the ESO Science Archive Facility services (data obtained under request number 396301).

This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in (Ochsenbein, Bauer & Marcout 2000).

DATA AVAILABILITY

The TESS data underlying this article were accessed from MAST (Barbara A. Mikulski Archive for Space Telescopes) Portal (https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html), including the data products found in the bulk download website (http://archive.stsci.edu/tess/bulk_downloads/bulk_downloads_ffi-tp-lc-dv.html). Part of the data were derived from sources in public domain as given in the respective footnotes. The derived data generated in this research and the code used for the photodynamical analysis will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: MID-ECLIPSE TIMES OF THE FOUR ECLIPSING BINARIES