Six New Compact Triply Eclipsing Triples Found With TESS

With the advent of long-term, wide field, precision photometry from space with such missions as Kepler (Borucki et al., 2010), K2 (Howell et al., 2014), and TESS (Ricker et al., 2015), it has become relatively easy to discover triply eclipsing triple star systems. These are often found when an extra, isolated pair of eclipses appear in the lightcurve of an ordinary eclipsing binary (EB), or a long exotic-looking extra eclipse appears that cannot be produced in a simple binary (see the recent extensive review of Borkovits 2022). We refer to these as ‘third-body’ events where either the EB occults the third star (hereafter, the ‘tertiary’) in its outer orbit, or vice versa. Typical eclipse periods for the inner EBs are days, while the period for the extra third-body eclipses range from a month to about a year. When one of these triples is found, no further vetting of the object is typically needed before concluding that this is a bound triple system (or possible a higher stellar multiple). By contrast, when a pair of EBs is found in the same photometric aperture, it is not immediately clear whether the two EBs are physically bound or are simply close together in the sky by chance (see the quadruples catalog of Kostov et al. 2022), and further vetting in the form of, e.g., radial velocity measurements (RVs) or eclipse timing variations (ETVs) is required.

Once a compact triple system has been identified via its third body eclipses, several additional characteristics can quickly become apparent about the system. First, if more than one outer eclipse of the same type¹¹1As in ordinary EBs, outer eclipses come in two flavors—primary and secondary eclipses. are seen in succession, then the outer orbital period of the triple is immediately revealed. If both the primary and secondary outer eclipses are seen, then, just as in an ordinary EB, the quantity $e_{\rm out}\cos\omega_{\rm out}$ can be measured, where $e_{\rm out}$ and $\omega_{\rm out}$ are the eccentricity and argument of periastron of the outer orbit, respectively. Finally, the presence of both inner and outer eclipses gives some good indication that the binary orbital plane and the outer orbital plane are at least somewhat aligned (i.e., the systems tend to be ‘flat’) or else the likelihood of detecting both sets of eclipses is relatively lower.

As we and others have shown in a number of previous papers (see, e. g. Carter et al., 2011; Borkovits et al., 2013; Masuda et al., 2015; Orosz, 2015; Alonso et al., 2015; Borkovits et al., 2019a, 2020b, 2022), the outer eclipses have encoded in them a substantial amount of information which, when combined with supplementary data, can ultimately lead to a complete description of the stellar properties (masses, radii, $T_{\rm eff}$, age, and metallicity) as well as the complete orbital configuration of the system. The supplemental material can involve the EB lightcurve itself, the extracted ETV curve, spectral energy distribution (SED) measurements from archival surveys, ground-based photometric surveys, and RVs.

Another feature of compact triple star systems that makes them fascinating objects to study is the relatively short timescales for dynamical interactions. These can occur over a year, a few months, or even weeks. Interesting effects to look for include dynamical as well as light-travel time delays in the ETVs of the EBs, forced apsidal motion in the EB, orbital plane precession if the two orbital planes are misaligned, and even large amplitude von Zeipel-Kozai-Lidov cycles (von Zeipel 1910; Lidov 1962; Kozai 1962) in the case of strongly misaligned orbital planes.

Triple star systems are also interesting in terms of the insight they provide about the formation and subsequent evolution of multistellar systems (see, e.g., Tokovinin 2021; Borkovits 2022; Sect. 7). They are the next simplest entity after binary star systems. In some ways they are analogous to studying He atoms after mastering H atoms. However, while there are more than a million eclipsing binary systems known (see, e.g., Sect. 2; Powell et al. 2021; E. Kruse, 2022 in preparation), the number of triply eclipsing triple systems in the literature, is currently under 20.

Here we present the discovery and detailed analyses of six new compact triply eclipsing triple star systems. In Section 2 we discuss how the discovery of the third body events were made using the TESS data, and present plots of the third body events. Archival spectral energy distributions (SED) are then used in Section 3 to make first estimates of the constituent stellar masses, radii, and $T_{\rm eff}$. We then use archival photometric data from a number of ground-based surveys to determine the outer orbital period of the triples via the third-body eclipses (see Sect. 4). The detailed photodynamical model by which we analyze jointly the photometric lightcurves, eclipse timing variations, and spectral energy distributions is reviewed in Section 5. The system parameters for each of the six triple systems are presented in Section 6 in the form of comprehensive tables, including extracted masses, radii, and effective temperatures, as well as the orbital parameters for both the inner and outer orbits. We summarize our results and discuss a few of the salient findings from our study in Section 7. We also discuss how our compact triple systems may inform us about multistellar formation and evolution.

Our ‘Visual Survey Group’ (VSG; Kristiansen et al. 2022) continues to search for multi-stellar systems in the TESS lightcurves. We estimate that, thus far, we have visually inspected some 10 million lightcurves from TESS. Such visual searches are a complement to more automated ones using machine learning algorithms (see, e.g., Powell et al. 2021; Kostov et al. 2021; Kostov et al. 2022). The lightcurves are displayed with Allan Schmitt’s LcTools and LcViewer software (Schmitt et al., 2019), which allows for an inspection of a typical lightcurve in just $\sim$5 seconds. It is important to note that 9 million of the studied lightcurves were of anonymous stars, while 1 million lightcurves were of preselected eclipsing binaries that were found in the TESS data via machine learning searches (see Powell et al. 2021; E. Kruse, 2022 in preparation).

For our survey work we largely made use of lightcurves from the following sources: Science Processing Operations Center (SPOC, Jenkins et al. 2016); the Difference Imaging Pipeline (Oelkers & Stassun, 2018); the PSF-based Approach to TESS High quality data Of Stellar clusters (PATHOS, Nardiello et al. 2019); the Cluster Difference Imaging Photometric Survey (CDIPS, Bouma et al. 2019); the MIT Quick Look Pipeline (QLP, Huang et al. 2020); the TESS Image CAlibrator Full Frame Images (TICA, Fausnaugh et al. 2020); and the Goddard Space Flight Center (GSFC, see Sect. 2, Powell et al. 2021).

The first signatures that are looked for in terms of identifying triply eclipsing triples are an eclipsing binary lightcurve with an additional strangely shaped extra eclipse or rapid succession of isolated eclipses. One gratifying aspect of finding triply eclipsing triples is that they are in a sense ‘self-vetted’. In particular, there is no way for a single binary, or sets of independent stars or binaries to produce such ‘extra’ eclipsing events. Therefore, additional vetting becomes largely unnecessary in proving that these are indeed triples (or possibly higher order multiples).

While searching through the lightcurves obtained from the first three full years of TESS observations we have found more than $\sim$50 of these triply eclipsing triples. Of these we have determined the outer orbital period for 20 of them. We have previously reported on four of these systems (Borkovits et al. 2020b; Borkovits et al. 2022). He we present the discovery and analysis of six new triply eclipsing triples from among this set: TIC 37743815, TIC 4256558, TIC 54060695, TIC 178010808, TIC 242132789, and TIC 456194776. (See Table LABEL:tbl:mags for the main catalog data of the targets.)

All six targets were measured in full frame images from TESS with either 30-min or 10-min cadence. A portion of the TESS lightcurves for all six sources are shown in Fig.1. The sectors during which these sources were observed with TESS are summarized in Table LABEL:tbl:sectors. TICs 37743815, 42565581, 178010808, and 242132789 were observed during two widely separated sectors each, and a third body event was observed for each source in both sectors, except for TIC 178010808, for which only one third-body event was detected. TIC 456194776 was observed in only one sector. Finally, TIC 54060695 was observed during three sectors, with a third body event detected in each—two secondary and one primary outer eclipses.

Given that the sectors are only approximately a month long and the outer orbital periods range from 42 to 123 days, it is somewhat fortuitous that we managed to observe third body events in nearly all the sectors. However, most of these systems exhibit two eclipses per outer orbit, and we have selected these for presentation in this work precisely because the outer periods are relatively short and therefore easy to detect even in archival data sets. In other words, there are several selection effects at work here.

Once we have discovered a triply eclipsing triple, we would like to develop some initial estimates of the nature of the three stars in the system. To do this we make use of an analysis of the spectral energy distribution (SED). We utilize the VizieR (Ochsenbein et al. 2000; A.-C. Simon & T. Boch: http://vizier.unistra.fr/vizier/sed/) SED service which, in turn, utilizes systematic sky coverage of such surveys as Skymapper (Wolf et al., 2018), Pan-STARRS (Chambers et al., 2016), SDSS (Gunn et al., 1998), 2MASS (Skrutskie et al., 2006), WISE (Cutri et al., 2013), and in some cases Galex (Bianchi et al., 2017). These typically provide $\sim$20 fluxes over the range 0.35 to 21 microns.

Unless we have specific information to the contrary, we assume for our preliminary analysis that the three stars in the system have evolved in a coeval fashion since their birth as a triple system. We further assume that there has been no mass transfer among the three stars, in particular between the binary components. Under these assumptions, there are only four parameters that need to be fitted via a Markov chain Monte Carlo approach (see, e.g., Ford 2005; Rappaport et al. 2021): $M_{\rm Aa}$, $M_{\rm Ab}$, $M_{\rm B}$, and the age of the system, where $Aa$ and $Ab$ refer to the stars in the inner binary, while $B$ is the tertiary star in the outer orbit. We also make use of MIST stellar evolution tracks (Paxton et al. 2011; Paxton et al. 2015; Paxton et al. 2019; Dotter 2016; Choi et al. 2016) for an assumed solar composition²²2Adopted for this preliminary analysis only. See Sect. 5 for a description of the full, and more general, photodynamical analysis., as well as stellar atmosphere models from (Castelli & Kurucz, 2003). If these four parameters can be determined, then the evolution tracks simultaneously determine the stellar radii and effective temperature of all three stars.

In order to fit an SED, one typically requires an accurate distance to the source and the corresponding interstellar extinction, $A_{V}$. Since Gaia (Gaia collaboration, 2021) provides a secure distance with typically better than 5% accuracy, and we can find information on the extinction from a variety of sources (e.g., Bayestar19; Green et al. 2019), in principle we do not need to fit for these parameters. However, we usually add these two quantities to the fitted MCMC parameters, but with priors limited to just $\pm$ 4 times the listed uncertainties on them.

In spite of having some 20 SED data points to work with, and only 4-6 free parameters to fit for, this is most often quite insufficient for a decent solution. The reasons are that (i) many of the SED points are sufficiently close to each other in wavelength so that they are not really independent, and (ii) a typical SED curve contains essentially only 4 or 5 defining characteristics, e.g., flux at any given wavelength on the Rayleigh-Jeans tail, wavelength of the peak in the SED curve, sharpness of the falloff at short wavelengths, etc. Therefore, we find it extremely helpful to add a few supplementary constraints in the MCMC fit. These include estimates of the radius and $T_{\rm eff}$ of the tertiary and the ratio $T_{\rm eff,Aa}/T_{\rm eff,Ab}$. Because all of the tertiary stars in this work are giants³³3The least evolved tertiary is the one in TIC 178010808, a 1.6 M_⊙ star of radius 2.9 R_⊙ whose luminosity exceeds that of the combined EB stars by nearly a factor of 2., they tend to dominate the light from the system. Therefore, we make use of the Gaia DR2 and TESS Input Catalog (TIC v8.2) estimates of the ‘composite’ radius and $T_{\rm eff}$ of the entire system measured as a single object (see Table LABEL:tbl:mags). However, because those estimates are hardly a perfect representation of the tertiary, i.e., there are two other stars in the system, though of considerably lower luminosity, we take uncertainties on $T_{\rm eff}$ to be $\pm 300$ K and on $R_{B}$ to be $\pm 2\,R_{\odot}$. Finally, the temperature ratio of the two EB stars can be estimated approximately from the ratio of their eclipse depths.

Using only this limited information, we fit the SED for all the properties of the stars in our six systems. The results are shown in Fig. 2. In each case, the measured SED points (orange circles) have been corrected for interstellar extinction. The continuous green curves represent the model fits for the stars in the EB, while the cyan curve is for the tertiary. The heavier black curve is the total flux from all three stars. One can get a very good sense from these plots just how much the tertiary dominates the system light. We list on the plots only the nominal best-fitting values for $M$, $R$, and $T_{\rm eff}$ for each of the three stars. We do not list the corresponding uncertainties on these parameters here because we employ a more comprehensive photodynamical analysis in Sect. 5 which utilizes the EB and third-body eclipse contributions to the lightcurve, as well as the SED, to determine the final and more accurate stellar parameters. Nonetheless, these first estimates of the stellar parameters provide a very quick estimate of what kinds of stars we are dealing with. The formal uncertainties on the masses are typically 10%. These parameters can then be utilized as the initial input guesses to the photodynamical analysis.

We utilize the system parameters for the stars found from this preliminary SED analysis to show in Fig. 3 the locations of the stars superposed on the MIST evolution tracks. The tertiaries range from $\sim$3 to 15 R_⊙, with $T_{\rm eff}$ ranging from 6000 K to 4500 K, respectively. By contrast, most of the binary stars are still on the main sequence, and, with only one exception, range in mass from 1.0-1.8 M_⊙. Within four of the systems the two EB stars tend to have similar masses.

As a final note on the SED analysis, we point out that we have not considered pre-MS solutions. These are considered and rejected in the photodynamical analysis.

Once we have discovered a triply eclipsing triple system with TESS, the most important question to answer after determining the basic parameters of the constituent stars is the nature of its outer orbit, in particular the period and eccentricity. For this purpose, in the absence of RV data, we most often make use of the ASAS-SN (Shappee et al. 2014; Kochanek et al. 2017) and ATLAS (Tonry et al. 2018; Smith et al. 2020) archival data sets. The ASAS-SN data sets typically have $\sim$1500-3000 photometric measurements of a given target, while the ATLAS archives often have approximately 1700 photometric measurements. The ATLAS data have the advantage of going somewhat deeper than the ASAS-SN data, but the disadvantage of saturating on brighter stars where ASAS-SN may still perform well. When these two data sets are of comparable quality, we typically add them. Naturally, we also check for KELT (Pepper et al., 2007, 2012), WASP (Pollacco et al., 2006), HAT (Bakos et al., 2002) and MASCARA data (Talens et al., 2017) to see whether they are available for a particular source.

For all the sources, with the exception of TIC 242132789, there were WASP archival data available. We found these to be quite useful in helping to determine the long-term average EB periods. But, it turns out that these data were generally too noisy to aid in the search for the outer third body eclipses, except in the case of TIC 178010808 where the WASP data nicely complemented the ATLAS and ASAS-SN data. For TIC 242132789 there was a set of KELT data in addition to ATLAS and ASAS-SN data. The KELT data marginally detected the spot and ellipsoidal light modulations associated with the outer orbit, but not the outer eclipses.

We do a blind search for the outer eclipses (either outer primary or outer secondary eclipse) using a Box Least Squares transform (Kovács et al., 2002). Before doing the BLS search we remove the lightcurve of the inner EB by Fourier means (as described in Powell et al. 2021). In the process we remove between 5 and 100 orbital harmonics depending on the sharpness of the features in the EB lightcurve. This cleaning process requires knowing the orbital period of the EB very accurately. In turn, we determine the long-term average binary period from the TESS data or from the archival data, whichever yield a more precise result.

We show the results of the above procedure for each of our six triply eclipsing triple stars in Fig. 4. In each panel we show a folded, binned, and averaged lightcurve of the archival data about the period corresponding to the largest and most significant peak in the BLS transform. In all cases, the zero phase for the outer orbit is taken to be the time of one of the third-body eclipses observed in the TESS data. In four of the six sources, both the primary and secondary outer eclipses are clearly detected. In those cases, the value of $e\cos\omega_{\rm out}$ is also accurately determined, in addition to the outer period of the triple. In the fifth source, TIC 178010808, the secondary outer eclipse is only barely detected, if at all. In the sixth source, TIC 242132789, the outer orbital lightcurve has a clear undulating structure superposed on the very clear primary eclipse. Because the slow modulations have the same period as the outer eclipses, we attribute this to starspot(s) on a giant tertiary that is corotating with the binary orbit as well as to ellipsoidal light variations from the giant. Note that the giant in this system has $R\simeq 15\,{\rm R}_{\odot}$. Not only is this star large, but its outer orbital period of 42 days is the shortest among our sample, and one of the shortest period triples known.

For all six of our triply eclipsing triples, we have carried out a photodynamical analysis with the software package Lightcurvefactory (see, e.g. Borkovits et al., 2019a, 2020a, and references therein). As described in earlier work, the code contains (i) a built-in numerical integrator to calculate the three-body perturbed coordinates and velocities of the three stars in the system; (ii) emulators for the TESS light curve, the ETVs extracted therefrom, and radial velocity curve (if available), and (iii) an MCMC-based search routine for the system parameters. The latter utilizes an implementation of the generic Metropolis-Hastings algorithm (see e.g. Ford, 2005). The use of this software package and the consecutive steps of the entire analysis process have been previously explained in detail as the code was applied to a wide range of multistellar systems (Borkovits et al., 2018, 2019a, 2019b, 2020a, 2020b, 2021; Mitnyan et al., 2020). These included tight and wider triple systems (with and without outer eclipses), as well as quadruple systems with either a 2+2 or 2+1+1 configuration. Here we discuss only a few specific points related to the current triples.

In relatively close systems, as we are studying here, perturbations to the EB orbit and the detailed profiles of the third body eclipses carry important information about the system parameters, including constraints on the masses and orbital elements. However, with only one exception in this current work, we have no RV measurements to help constrain the system parameters. Therefore, we adopt a somewhat different strategy. In the analysis we utilize some a priori knowledge of stellar astrophysics and evolution with the use of PARSEC isochrones and evolutionary tracks (Bressan et al., 2012). We make use of tabulated three-dimensional grids in triplets of {age, metallicity, initial stellar mass} of PARSEC isochrones that have stellar temperatures, radii, surface gravities, luminosities, and magnitudes in different passbands of several photometric systems. Then, we allow the three parameters {age, metallicity, initial stellar mass} to vary as adjustable MCMC variables. The stellar temperature, radius, and actual passband magnitude are calculated through trilinear interpolations from the grid points and these values are used to generate synthetic lightcurves and an SED that can be compared to their observational counterparts. This process, which is also built into Lightcurvefactory, is described in detail in Borkovits et al. (2020a). In our prior work, we have termed these solutions ‘MDN’ (model-dependent-no-RV solutions) which is what we implement here.

Regarding the technical details, in the case of the stellar evolution model dependent runs, the freely adjusted (i.e., trial) parameters were as follows:

• (i)

  Stars: Three stellar masses and the ‘extra light’ contamination, $\ell_{4}$, from a possible fourth star (or any other contaminating sources in the TESS aperture). Additionally, the metallicity of the system ([$M/H$]), the (logarithmic) age of the three stars ($\log\tau$), the interstellar reddening $E(B-V)$ toward the given triple, and its distance, were also varied.

• (ii)

  Orbits: Three of six orbital-element related parameters of the inner, and six parameters of the outer orbits, i.e., the components of the eccentricity vectors of the two orbits $(e\sin\omega)_{\mathrm{in,out}}$, $(e\cos\omega)_{\mathrm{in,out}}$, the inclinations relative to the plane of the sky ($i_{\mathrm{in,out}}$), and moreover, three other parameters for the outer orbit, including the period ($P_{\mathrm{out}}$), time of the first (inferior or superior) conjunction of the tertiary star observed in the TESS data ($\mathcal{T}_{\mathrm{out}}^{\mathrm{inf,sup}})$ and, finally, the longitude of the node relative to the inner binary’s node ($\Omega_{\mathrm{out}}$).

Here we add some additional notes about the ‘age’ and the ‘distance’. First, regarding the ‘age’ parameter, our previous experience has led us to believe that, in some cases, it is better to allow the ages of the three stars to be adjusted individually instead of requiring strict coeval evolution. This issue was briefly discussed in Rowden et al. (2020) and Borkovits et al. (2021), and we discuss it below in the case of some individual sources. Regarding the distance of the triple system, one can argue that the accurate trigonometric distances obtained with Gaia (Bailer-Jones et al., 2021) should be used as Gaussian priors to penalize the model solutions. But, given that neither DR2 nor the recently released eDR3 Gaia parallaxes have been corrected for the multistellar nature of the objects, we consider the published parallaxes and corresponding distances to be not necessarily accurate for our systems. Therefore, we decided not to utilize the Gaia distances. Instead, we constrained the distance by minimizing the $\chi^{2}_{\mathrm{SED}}$ value a posteriori, at the end of each trial step (for results see Sect. 6.1).

A couple of other parameters were constrained instead of being adjusted or held constant during our analyses. Specifically, the orbital period of the inner binary ($P_{\mathrm{in}}$) and its orbital phase (through the time of an arbitrary primary eclipse or, more strictly, the time of the inferior conjunction of the secondary star – $\mathcal{T}^{\mathrm{inf}}_{\mathrm{in}}$) are in this category. They were constrained internally through the ETV curves.

Regarding the atmospheric parameters of the stars, we handled them in a similar manner as in our previous photodynamical studies. We utilized a logarithmic limb-darkening law (Klinglesmith & Sobieski, 1970) for which the passband-dependent linear and non-linear coefficients were interpolated in each trial step via the tables from the original version of the Phoebe software (Prša & Zwitter, 2005). We set the gravity darkening exponents for all late type stars to $\beta=0.32$ in accordance with the classic model of Lucy (1967) valid for convective stars and hold them constant. For several of our systems, however, the analysis of the net SED has revealed that the EBs contain hotter stars, having radiative envelopes. For these stars, we set $\beta=1.0$. The choice of this parameter, however, has only minor consequences, as the stars under the present investigation are close to spheroids.

For the photodynamical analysis we utilized a more sophisticated processing for the TESS photometric data by employing the convolution-based differential image analysis tasks of the FITSH package Pál (2012). Furthermore, we note that in preparing the observational data for analysis, to save computational time we dropped out the out-of-eclipse sections of the 30-min cadence TESS lightcurves, retaining only the $\pm 0\aas@@fstack{p}15$ phase-domain regions around the binary eclipses themselves. However, during sections of the data containing the third-body (i.e., ‘outer’) eclipses, we kept the data for an entire binary period both before and after the first and last contacts of the given third-body eclipse.

Moreover, the mid-eclipse times of the inner binaries, used to define the ETV curves, were calculated in a manner that was described in detail in Borkovits et al. (2016). The ETVs were one of the inputs to the photodynamical analysis. We tabulate all the eclipse times in Appendix C as online only tables.

Finally, note also some system specific departures from the standard procedures described above in the case of two of our triples. These are as follows:

(1) In the case of TIC 242132789, the lightcurve of both TESS sectors display non-linear, somewhat erratic variations in the mean out-of-eclipse flux levels. This can partly be attributed to the ellipsoidal light variations (ELV) of the red giant tertiary, which is handled internally by our lightcurve emulator and, therefore, nothing special has to be done to model the ELVs. However, an additional, more erratic contribution to these variations might arise from either time-varying spot activity on the surface of the red giant or any stray light in the aperture (or both). In any case, independent of the origin(s) of this time-varying, irregular contaminating flux, we modelled it by fitting an eighth order polynomial, separately for the two sectors, simultaneously with the triple star lightcurve modelling, during each MCMC step.

(2) The other triple that was handled somewhat differently is TIC 456194776. This system was originally the fourth target of the ground-based photometric follow-up campaign that was described in detail in Sect. 2.3 of Borkovits et al. (2022). In contrast to the other three triply eclipsing triples, of which the detailed analyses were published in Borkovits et al. (2022), unfortunately, we were unable to catch any further third-body eclipses during our observing runs. This is the reason why this system was not included the above-mentioned study, but rather appears in this work. On the other hand, however, we were able to observe 9 regular inner eclipses from the ground between 12 August 2020 and 12 December 2021. Thus, we decided to include these eclipsing lightcurves (in Cousins $R_{C}$-band) and also the extracted eclipse times into our analysis. In addition, near the final stages of our analysis we also acquired RV data for this target. Therefore, we carried out a second kind of analysis for this target which we have termed ‘MDR’ (model-dependent-with-RV) with the inclusion of these RV points. We discuss this latter, MDR solution and compare it with the MDN results in Appendix A.2. This second type of analysis might also be called a ‘spectro-photodynamical solution’.