New compact hierarchical triple system candidates identified using Gaia DR3 ⋆

Aims. We introduce a novel way to identify new compact hierarchical triple stars by exploiting the huge potential of Gaia DR3 and also its future data releases. We aim to increase the current number of compact hierarchical triple systems significantly. Methods. We used several eclipsing binary catalogs from different sky surveys that list a total of more than 1 million targets to search for Gaia DR3 non-single-star orbital solutions with periods substantially longer than the eclipsing periods of the binaries. Those solutions in most cases are likely to belong to outer orbits of tertiary stars in those systems. We also attempted to validate some of our best-suited candidates using TESS eclipse timing variations. Results. We find 403 objects with suitable Gaia orbital solutions of which 27 are already known triple systems, leaving 376 newly identified hierarchical triple system candidates in our sample. We find the cumulative probability distribution of the outer orbit eccen-tricities to be very similar to those found in earlier studies based on observations of the Kepler and OGLE missions. We find measurable nonlinear eclipse timing variations or third-body eclipses in the TESS data for 192 objects which we also consider to be confirmed candidates. Of these, we construct analytical light-travel time effect models for the eclipse timing variations of 22 objects with well-sampled TESS observations. We compare the outer orbital parameters from our solutions with those from the Gaia solutions and find that the most reliable orbital parameter is the orbital period, while the values of the other parameters should be used with caution.


Introduction
Hierarchical triple stellar systems consist of three stars that are organized as an inner, binary pair with a more distant, outer (tertiary) component. These systems are exceptional in that their orbital and astrophysical parameters can often be determined precisely, especially for those in which the inner binary either eclipses or is eclipsed by the tertiary during their orbital revolution. These highly accurate parameters can be used to refine stellar evolutionary models and to study the formation and evolution of such systems from which all fields of stellar astrophysics can benefit. These systems can also be of help in our understanding of some peculiar objects with compact components (e.g., Type Ia supernovae, Kushnir et al. 2013). This is because the formation of their progenitor close binaries requires some effective orbit-shrinking mechanism during their lifetime that in many cases could be explained by dynamical interactions with a third ⋆ Full Table A.1 is only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/670/A75 stellar component in the system. The formation and evolutionary scenarios of triple systems have recently been reviewed by Tokovinin (2021), and Toonen et al. (2020Toonen et al. ( , 2022, respectively. These so-called dynamical interactions are more common and easier to detect within the scale of human lifetimes in compact hierarchical triple (CHT) systems where the outer orbital period is less than ∼1000 days. There is still a relatively small number of systems known in this category; the latest version of the Multiple Star Catalog (Tokovinin 2018) lists 201 such systems, but there are also a few quadruple or higher-order hierarchical systems included in this number as well. Nevertheless, nowadays, thanks to the various ground-based and more importantly space-based photometric all-sky surveys, the numbers of compact hierarchical systems are rapidly growing. For example, Hajdu et al. (2019Hajdu et al. ( , 2022 identified 177 and 16 such CHTs, respectively, using data obtained by the OGLE (Optical Gravitational Lensing Experiment) survey. A recent comprehensive review of the different methods and results of the field of discovering compact hierarchical systems can be found in Borkovits (2022).
Most of the CHT stellar systems were discovered by either photometric or spectroscopic studies and detection of the astrometric wobble of any of the components of such systems was not feasible as part of broad searches. The third data release (DR3; Gaia Collaboration 2023; Babusiaux et al. 2023) of the Gaia space telescope (Gaia Collaboration 2016) is unique in the sense that it contains a catalog of objects for which the time-domain, ultra-precise astrometric observations allowed the Gaia team to constrain a simple Keplerian two-body orbital solution for 135 760 objects using their variable astrometric positions on the plane of the sky. Moreover, this data release also contains the first spectroscopic results of Gaia in which the time-domain radial-velocity measurements also resulted in 186 905 objects with a similar kind of two-body orbital solution. Finally, there are some objects for which eclipses could be detected in the Gaia photometric measurements and these eclipses made it possible to construct such a two-body orbital solution for 86 918 objects. The previously mentioned solutions are of single type, but there are also systems with two types of observational data for which a combined orbital solution was found, namely "astrometric+spectroscopic" (33 467 objects) and "eclipsing+spectroscopic" (155 objects). Detailed information about the modeling steps can be found in the corresponding section of the Gaia DR3 online documentation (Pourbaix et al. 2022). All of the orbital solutions and associated parameters for these binary or multiple star candidates are publicly available and are listed in the Vizier tables of the Gaia DR3 non-single-stars (NSS) catalog (Gaia Collaboration 2022).
Our goal is to identify possible NSS models in which the orbital solution may belong to an outer orbit of a potential hierarchical triple (or multiple) system. These are detected through either the motion of the unresolved inner binary or the outer tertiary star, or the light centroid of both. For this purpose, our main idea was to collect as many eclipsing binaries (EBs) as possible with previously determined orbital periods from different catalogs found in the literature and then to search for Gaia orbital solution periods that are at least a factor of five longer than the corresponding EB period. In Sect. 2, we describe the data sets we used for our search along with the steps of their analysis. In Sect. 3, we summarize and discuss the main results of our search for triple candidates, while in Sect. 4, we present a possible validation method of the newly found candidates using eclipse timing variations (ETVs) calculated from TESS data. Finally, we highlight our main conclusions in Sect. 5.

Catalog data and analysis: Cross-matching and identification of candidates
We downloaded the full data sets of all main catalogs containing lists of EBs, namely APASS (AAVSO Photometric All Sky Survey; Munari et al. 2014), ASAS-SN (All-Sky Automated Survey for Supernovae; Shappee et al. 2014;Rowan et al. 2022), GCVS (General Catalog of Variable Stars; Samus' et al. 2017), Kepler (Kepler Eclipsing Binary Stars; Kirk et al. 2016), TESS (TESS Eclipsing Binary Stars; Prša et al. 2022), and VSX (The International Variable Star Index; Watson et al. 2006). The total number of objects exceeded 1 million EBs (see Table 1). We cross-matched these lists with the Gaia EDR3 source catalog 1 by coordinates via TOPCAT (Tool for OPerations on Catalogues  And Tables; Taylor 2005). We used the default 5 ′′ search radius for this purpose, which yielded more sources than the original lists, implying a number of duplicates. We did not immediately filter these out as it was easier and more efficient to remove them at the end of the analysis (see the text later). After cross-matching with the Gaia EDR3 catalog, we downloaded all the Vizier tables in the NSS catalog containing the different types of NSS solutions. We then performed a search for each Gaia source ID that was in our lists and collected all available orbital model parameters for each source. In order to make a more complete survey, we also searched for stars with a relatively small angular separation 2 from the unresolved EBs in our list to see if we could find a bound tertiary component on a wider orbit around them that is resolved by Gaia. However, we could not find any system that met such criteria, mainly because of the primary selection of the astrometric processing steps applied by the Gaia team, which makes it unlikely to find any such systems (Halbwachs et al. 2023;Riello et al. 2021).
As a next step, we kept all objects for which the Gaia orbital period solution was at least five times longer than the orbital period found in the corresponding EB catalog. We chose this condition as a very loose threshold because currently the tightest known triple system has an outer-to-inner period ratio of 5.4 (Xia et al. 2019) and because of dynamical stability issues. It is very unlikely that we will find an object below our chosen period ratio limit (see Borkovits et al. 2022). Finally, in order to filter out duplicate sources, we joined the lists of objects we kept from the different catalogs and filtered out those source IDs that appeared in this complete list multiple times (i.e., those that are listed in multiple catalogs under different names) to have only a single occurrence of each distinct object. We also removed those systems for which multiple different Gaia source IDs were found and kept only those for which the brightness was consistent with the brightness of the EB we searched for. We also noticed that some EBs had been listed with incorrect coordinates, as we found Gaia source IDs for them based on those incorrect coordinates that were different from their real source IDs. We found their real source IDs based on their common names instead of their coordinates listed in the corresponding catalogs. We discarded these systems as well, because neither of them actually had a Gaia NSS solution that fulfilled our search criteria.

Summary of the identified candidates
As a result of our study, we have identified 403 hierarchical triple candidates among formerly known EBs from the literature based  on the periods of their Gaia NSS orbital solutions. The first 35 candidates and the parameters of their corresponding Gaia NSS orbital solutions are listed in Table A.1, while the full list is available at the CDS. There are 100 pure spectroscopic, 267 pure astrometric, and 31 combined Gaia NSS solutions in our sample. There are also three objects with two separate astrometric and spectroscopic solutions along with two other objects with separate astrometric and combined solutions. For objects with two different types of orbital solutions, we chose to use the pure astrometric solution in every case as a convention. In Fig. 1, we illustrate the distribution of the inner and outer orbital periods (the periods of the corresponding NSS solutions) for all of our candidates. Broadly speaking, the outer periods range from somewhat greater than 1000 days down to just a few days, though most of them are concentrated between 300 and 1000 days. The sharp falloff in systems above ∼1000 days is an artefact of the Gaia observational window. The inner binary periods range from a fraction of a day to 10 days, simply representing the observational selection effects of the underlying EB catalogs that we used. Most of the longer outer periods are based on either a Gaia astrometric or spectroscopic solution, while essentially all of the outer orbits with periods of less than ∼100 days are pure Gaia spectroscopic solutions. We also indicate our chosen limiting period ratio of P out /P in = 5 with a dashed line. The closer a system is to this line, the tighter it is. We note that there is one point below this limit for which we provide an explanation in Sect. 4.1. Prior to this study, there were some 400 known triple-star systems with outer periods of ≲1000 days, and only about 30 known with periods of ≲200 days. These systems are discussed and tabulated in Tokovinin (2018), Borkovits et al. (2016), Hajdu et al. (2019), Hajdu et al. (2022), and Fig. 1 of Borkovits et al. (2022). Therefore, this study has the potential to substantially increase the number of known triple-star systems with measured outer orbits.
The upper panel of Fig. 2 shows the differential eccentricity probability distribution for the 403 triple candidates that we identified. The distribution is crudely described by the function dN/de ≃ 8.34 e 0.618 exp [−5.15e 1.618 where e is the outer eccentricity (see below for the origin of this functional form). This curve is superposed (after renormalization with the number of objects and the bin width) on the eccentricity distribution in Fig. 2, and provides a good fit considering the limited statistics. In the lower panel of Fig. 2 we show the cumulative probability distribution for the outer eccentricities that we found. For comparison, we also show the eccentricity distributions for the Kepler sample of triple systems (Borkovits et al. 2016) and the outer orbits for 205 secure triple systems identified in the OGLE sample (Hajdu et al. 2019). In the same plot, we compare these three empirical distributions with those of a flat eccentricity distribution (dN/de = constant) and a so-called "thermal distribution" (dN/de ∝ e). The empirical distributions appear to be fairly similar to each other but distinctly different from a differential distribution that is constant (linear, as in a cumulative distribution). Moreover, we note that the empirical distributions are dramatically distinct from the thermal distribution presented by Jeans (1919) for binary stellar systems that have reached a state of energy equipartition through a large series of dynamical encounters. In fact, the empirical distributions fairly closely resemble the empirical eccentricity distribution for ordinary binary stellar systems (see The cumulative eccentricity distribution for the Gaia-based outer eccentricity distribution can be reasonably well represented by the following analytic expression: The derivative of this expression yields the above analytic approximation to the differential eccentricity distribution. There is one conceivable observational selection effect for this data set that concerns the eccentricity distribution. That is, large eccentricities generally tend to make the triple star system more unstable. From Eq. (16) of Rappaport et al. (2013), which is based on expressions for dynamical stability by Mikkola (2008) and Mardling & Aarseth (2001), and we see that for period ratios of 1000, 300, 100, 30, and 10, the limiting eccentricities are e = 0.94; 0.88; 0.78; 0.58, and 0.28. Therefore, care should be taken when interpreting the relative dearth of high-eccentricity systems.
In the sections below, we discuss what we have done to try to validate the triple-star nature of some of these systems using ETV data. Additionally, as it turns out, some 27 of these systems are previously identified triples or multiples (see Table 2). We also discuss the "truth probability" for the SB1 solutions evaluated from the empirical formula of Bashi et al. (2022). Thus, through these we acquire reasonable confidence that a substantial fraction (e.g., ≳50%, and probably as high as 75%) of the new candidates are actually valid triple stars with a small fraction of quadruples. If this is found to be the case, then we will have substantially increased the number of known compact triple-star systems. If even half of our new candidates are proven to be valid triples, then this study will have increased the number of known compact triples by approximately 1.5 times both below 1000 days and below 200 days. These compact triple systems, with RV orbits, third-body eclipses, and/or measurable dynamical interactions, are the most likely multi-stellar systems to provide a means to comprehensive measure all of their stellar and orbital parameters. Here we hope to demonstrate that finding Gaia outer orbits associated with known EBs is a novel way to find new compact triples.

Systems already known as candidates
We checked the available publications in the literature for all single objects individually using SIMBAD (Wenger et al. 2000) and NASA ADS 3 . Most of these objects had not been involved in extensive studies, apart from their identification as EBs. In Table 2, we list all 27 systems among our candidates for which we find any information about their possible triple-star nature in the literature. Comparing the outer orbital periods from the Gaia solutions and the literature references, there are 14 objects (KIC 8330092, KIC 6525196, KIC 8043961, KIC 9665086, KIC 8719897, KIC 10991989, HD 181068, KIC 7955301, TIC 229785001, TIC 437230910, TIC 99013269, SAO 167450, HD 81919 and HD 190583) for which reasonable agreement can be found. We consider these systems as validated candidates that are confirmed by the Gaia mission as an independent source. The only exception is SAO 167450, which is a visual member of AA Ceti in a bound EB+EB (2+2) hierarchical quadruple system, for which the Gaia spectroscopic orbital solution found belongs to one of the inner EBs and not to the outer orbit of the two EBs around their common center of mass.
For seven systems, although clear signs indicate that they have a multiple nature, their outer orbital period was unknown from previous studies. For V0857 Her, the presence of an additional component was suggested based on the rotational broadening function extracted from its spectra. An obvious sinusoidal variation was detected in the ETVs of TIC 229082842 and TIC 59042291, but the corresponding data sets were insufficient to derive their outer orbital periods. The signals of two EBs were detected in the TESS light curve (LC) of TIC 305635022 that originate at the same pixels, making this source a candidate 2+2 system. The remaining three objects, TIC 24704862, TIC 301956407, and TIC 410654298, are all subsystems of known visual double or multiple stars. In our opinion, the orbital solutions found by Gaia could also be considered as independent confirmation of the multiple nature of these systems.
There are significant differences between the two listed periods of the remaining six systems (Gaia vs. the literature). Two of them, although cataloged as EBs, are in reality single stars with transiting exoplanets (TIC 438071843 and HATS-26). For these systems, it is possible that the Gaia orbital solutions belong to further stellar or exoplanet components on wider orbits. Finally, for the remaining four candidates, the ratios of the two periods   (Gaia vs. the literature) are close to integers or half integers, which could indicate that at least one of the solutions over or underfitted the available sparse data for the same orbits; however, it is also possible that the listed orbital periods belong to different components in these systems.

Validation of some systems using TESS ETVs
In order to validate our list of new triple candidates with independent data and another method, we also collected TESS full-frame images (FFIs) for all objects from every available sector up to and including Sector 55. The LCs were extracted from the FFIs in an automated manner by applying a convolution-aided image-subtraction photometry pipeline (see e.g., Mitnyan et al. 2020) using FITSH (Pál 2012). The remaining nonastrophysical trends were mostly removed from the LCs by the WOTAN package (Hippke et al. 2019). In the case of LCs that were affected with significant stray light or an extra (probably false) signal from a nearby object, we performed a principal component analysis (PCA) that is more powerful in removing signals that originate outside the aperture used for the photometry. For this purpose, we used the built-in methods of the lightkurve (Lightkurve Collaboration 2018) Python package and its corresponding dependencies: astropy (Astropy Collaboration 2018), astroquery (Ginsburg et al. 2019), and TESScut . After that, we calculated the ETVs for each object from all their available eclipses found in the TESS LCs in the same way as described by Borkovits et al. (2015). During this process, we found that some of the EBs are cataloged with incorrect eclipsing periods, and so we list those systems in Table A.2 along with their revised EB periods based on TESS LCs. Next, we calculated a simple analytic light-travel-time effect (LTTE) model for all the available ETVs using the orbital parameters of the corresponding Gaia NSS solutions with the following equation: star in its outer orbit, ν 2 is the true anomaly of the orbit of the tertiary star, and c is the velocity of light. We note that we also estimated the expected contributions of the so-called dynamical delays with Eq. (12) of Borkovits et al. (2016) which were substantially smaller than the amplitude of the LTTE; a simple LTTE model should therefore be able to fit the ETVs well.
After that, we compared these LTTE models pre-calculated with the parameters from the Gaia solutions to the ETVs from TESS LCs. As our epoch and period determination are not always fully accurate, we also fitted a simple linear term for every model. As a first look, we noticed that none of the ETVs can be completely modeled with the parameters from the Gaia orbital solutions without slight changes. Nevertheless, for the majority of the systems, the outer orbital period of the Gaia solution seemed to be in good agreement with the ETVs.
As a next step, we tried to apply modifications to the parameters of the Gaia NSS solution parameters in order to match the model LTTE curves to the observed TESS ETVs. We tried to fit for the orbital parameters using a simple least-squares method, but because of the very limited temporal coverage of the ETV points available for most objects, the majority of these fits resulted in significantly different but equally acceptable solutions, which heavily depended on the initial model parameters. We therefore tried a different approach to fit the ETVs manually. For this purpose, we developed a Python Graphical User Interface (GUI) utilizing the PyQt5 4 module, which allows the user to change the model parameters individually and plot the resulting models interactively in real time in order to find, by visual inspection, the best candidates with sufficient orbital coverage to perform a fit. We then tweaked the values of the parameters from the Gaia solutions manually via a "chi-by-eye" in order to approximately match the ETV points, and then optimized the parameters further with the Nelder-Mead Simplex method using the scipy 5 Python package. Finally, we used these optimized parameters as initial parameters for a Markov-chain Monte Carlo (MCMC) sampling using the emcee (Foreman-Mackey et al. 2013) Python package. 4 https://pypi.org/project/PyQt5/ 5 https://docs.scipy.org/doc/scipy/index.html

Validated candidates with sufficient orbital coverage or third-body eclipses
The observing strategy of TESS strongly limits the available data for objects, as most of them are observed only once or twice over time-windows of ∼27 days within a TESS Year, but objects closer to the ecliptic poles can have more overlapping sectors with better data coverage with a quasi-continuous window of up to 1 yr. Moreover, the same hemispheres are only observed every other year, and therefore an apparent 1 year gap will appear even in the data set of the most well-observed objects. This means that our candidates also suffer from these observational effects and as a result we are not able to confirm the triple nature of all of them, even if they show significant nonlinear ETVs. Nevertheless, we found 22 candidates for which there is sufficient data coverage for their ETVs to be reliably modeled with a simple LTTE solution. We plot the ETVs and the corresponding simple LTTE solutions using the Gaia and our own adjusted parameters for two of them in Fig. 3 as examples, while the rest of them can be found in Figs. A.1-A.3. The corresponding model parameters can be found in Table A.3. As one can see in Fig. 3 (and the other such figures in the Appendix), most noticeably, the LTTE models calculated with the parameters of the Gaia orbital solutions can have either the same argument of periastron (ω NSS = ω ETV ) or be 180 degrees away (ω NSS = ω ETV ± 180 • ) from the LTTE models coming from our modeling process. This is not an error, and the reason for this is that the ETVs calculated from the TESS eclipses are always tracking the motion of the EB around the center of mass of the triple system, while for Gaia, this case is not as simple. For the astrometric solutions, Gaia uses the variations in the position of the photocenter of the object, while for spectroscopic solutions, it uses the variations in one or two detected peaks (for SB1 and SB2 solutions, respectively) in the cross-correlation function (CCF). This means that, without the ETV solutions (or in the case of spectroscopic solutions, the CCFs), it is not straightforward to say which component of the triple system is tracked by Gaia.
In order to compare the orbital parameters of the Gaia NSS and our own ETV solutions, we plot their correlations in Fig. 4 be almost perfectly aligned with the line that represents a ratio of unity. This means Gaia found the outer orbital periods of these triple systems almost perfectly. In the top right panel, the projected semi-major axes of the different solutions are compared and it can be seen that for this parameter, the relation is far from perfect. Only a handful of objects can be found near the line of unity, and almost all are more or less above this line. This could arise from factors related to the nature of the astrometric solutions: For example, for the majority of triple system parameters, the photocenter revolves around a smaller orbit, and therefore the semi-major axis of the outer orbit will be underestimated (see Eq. (C2) of Rappaport et al. 2022). Also, the ETV solutions can overestimate the semi-major axis if the orbits are not completely covered. However, for spectroscopic solutions, we cannot really be certain without reviewing the individual spectra 6 and the corresponding CCFs. In particular, the Gaia semi-major axis will BJD -2400000  In conclusion, of the 22 objects with the best-sampled TESS ETVs, one can state that, for triple candidates, the most reliable parameter from the Gaia orbital solutions is the outer orbital period and the outer eccentricity (although with higher uncertainties), while the other parameters should be used with caution, and only after identifying which component of the system is actually tracked by Gaia.
We also found four additional systems (TIC 66893949, TIC 298714297, TIC 88206187 and TIC 14839347) that show outer-orbit third-body eclipses which inherently confirms the triple nature of these objects. First, we checked the FFIs of the corresponding objects and the extra eclipse events originate from the same pixels as the regular inner binary eclipses. Moreover, we also found evidence in archival ASAS-SN, ATLAS (Tonry et al. 2018;Smith et al. 2020), and WASP (Butters et al. 2010) photometric data that these objects are indeed triply eclipsing triple systems. A comprehensive photodynamical analysis of these newly discovered systems will be presented in a dedicated paper in the near future. Moreover, we also identified four additional systems (TIC 457671987, TIC 224053059, TIC 376606423, and TIC 410654298) that are most likely previously unknown EB+EB quadruples according to their TESS LCs. Nevertheless, we do not consider them as validated systems because although the signal of both EBs for these objects comes from the same pixels, that does not necessary confirm that they are certainly bound in the same system. In particular, the two EBs could simply be seen in the same direction without any direct physical connection between them.
We would like to highlight the special case of TIC 410654298 here, which is the only object we kept despite the ratio of the Gaia and the EB period being slightly less than 5. This object is cataloged with a ∼20 day EB period and there are two separate types of Gaia NSS solutions (astrometric and combined astrometric+spectroscopic) available for it with a similar period of ∼170 day. These resulted in an initial period ratio of above 5. Nevertheless, after analyzing its TESS LC, we found that its actual EB period is ∼36.4 days, which yields a period ratio of below our chosen limit. However, as there are two different types of Gaia NSS solutions available for this system, we would assume that the orbit found by Gaia is valid. A plausible explanation for this object could be that it is a 2+2 system with one eclipsing and one noneclipsing component, and Gaia found the inner orbit of the latter, and therefore the 170 day period would not correspond to the outer orbit of the two binaries around the center of mass of the quadruple system.

Candidates with significant ETVs, but insufficient orbital coverage
There are 192 objects that show significant nonlinear ETVs, but only a very short part of their outer orbit is covered by TESS observations, and therefore they cannot be fitted reliably with a simple LTTE model. Although the LTTE models calculated with the parameters from the Gaia orbital solutions do not perfectly match the observed ETVs, they are at least in accordance regarding the outer orbital periods for most of these systems. These systems also demonstrate that the most trustworthy parameter from the Gaia NSS solutions is the orbital period. We plotted two systems in Fig. 5 as an example to illustrate how the duration of available observational data compares to the typical outer orbital solutions. As nonlinear ETVs are apparent in these systems, we also consider these systems to be confirmed hierarchical triple systems; nevertheless, we could not validate their outer orbital parameters from the Gaia NSS solutions.

Candidates without significant ETVs or external data
Finally, the remaining 158 systems (403 total minus 245 that are validated in one way or another) either do not have sufficient or any TESS observations available, or their ETV curves do not show any nonlinear behavior. In the absence of supplemental data, we are not able to confirm the triple nature of these systems, and we therefore continue to consider them as candidates that should be studied in the future when a sufficient amount of high-quality photometric (or spectroscopic) observational data becomes available.

Objects with short outer periods
The shortest outer period CHT currently known is λ Tau with a 33 day outer orbital period (Ebbighausen & Struve 1956). As can be seen in our Fig. 1, there are some 30 candidates with periods under 35 days. Empirically, such systems should be very rare, and so it is somewhat surprising that a significant number of our systems could belong to this category. However, it is hard to confirm the true nature of these systems using only TESS LCs and ETVs calculated from them. The amplitude of the LTTE depends on the semi-major axis, which is clearly small for such systems, while the dynamical delays are proportional to the outer eccentricity, which is also expected to be small according to formation scenarios. Therefore, the amplitude of the ETVs of such systems are expected to be very small, within the order of our uncertainties. A further difficulty is that there is a very limited number (or lack) of TESS observations of these systems. High-resolution spectroscopy with good temporal coverage could be of help in revealing the true configurations of these objects, but such data are not available for now as Gaia DR3 contains only an average of all spectra for each object. Nevertheless, Bashi et al. (2022), based on the spectroscopic observations of known EBs by LAMOST and GALAH, determined an empirical formula that can estimate the probability that a Gaia spectroscopic SB1 solution is valid. We plot the period distribution of all 92 of our SB1-type candidates along with their empirical "truth probabilities" in Fig. 6. For 64 and 36 systems this probability is at least 50 and 75%, respectively. In general, it can be seen that, the longer the outer period, the more likely the solution is to be valid. This means that among the candidates with short outer periods, we expect there to be several invalid solutions. Also, there is a possibility that if the Gaia orbital solution is valid, the object is actually a quadruple system with a 2+2 configuration with an eclipsing and a noneclipsing binary component, and Gaia actually found the spectroscopic orbital solution of the noneclipsing binary component (see the case of SAO 167450 in Sect. 3.2 and Table 2). However, for now, we cannot confirm or reject these candidates with certainty, because of insufficient supplemental observational data.

Conclusions
In this paper, we performed a search for possible hierarchical triple systems among more than 1 million cataloged EBs. We did this by looking for Gaia astrometric and/or spectroscopic NSS solutions with periods that are at least five times longer than the corresponding EB period. We identified 403 such triple candidates, of which 376 are newly proposed, while 27 are previously known candidates confirmed by Gaia observations. We collected all available TESS observations up to Sector 55 for all of our candidates and calculated their ETV curves. In the process, we discovered four newly found triply eclipsing triple systems simply from the extra eclipsing events found in their TESS LCs, and four additional 2+2 quadruple candidates. For 22 objects well covered by TESS observations, we were able to successfully fit the ETVs with a simple LTTE model using the parameters of the Gaia orbital solutions after slight modifications. These model solutions show that the most reliable parameter from the Gaia solutions is the orbital period. However, for the projected semi-major axis, the eccentricity, and the argument of periastron, we had to make larger changes in order to achieve acceptable fits to the ETVs, and we therefore think that their values taken from the Gaia solutions should be used with caution. There are also 192 objects among our candidates that show significant nonlinear ETVs. However, these could not be fitted reliably with a simple LTTE solution because of insufficient temporal coverage of the TESS data. Nevertheless, we consider these systems to be confirmed candidates as well, because of their nonlinear ETVs, even though they have no LTTE orbital solutions based on the ETVs. This makes a total of 218 reasonably confirmed CHTs, which is more than half the number of formerly known listed systems, that is, the 201 systems in the latest version of the Multiple Star Catalog (Tokovinin 2018) together with the 193 CHTs identified by Hajdu et al. (2019Hajdu et al. ( , 2022. The remaining 158 newly discovered systems remain in the status of triple candidates as there are no available supplemental TESS observations for them or they do not show detectable nonlinear ETVs. These systems should be revisited in the future after sufficient observational data become available for them in order to confirm their triple nature, along with those 192 systems for which the currently available ETV data do not allow us to determine their orbital parameters independently from the Gaia orbital solutions. We summarize these numbers in Table 3. Finally, we would like to emphasize the efficiency of our newly applied method in finding previously unknown hierarchical triple or multiple systems. Out of the 1 million EBs that we investigated, we identified 403 triple or multiple star candidates (0.04%). This is eight times higher compared to the efficiency of the Visual Survey Group (Kristiansen et al. 2022), who visually inspected 10 6 TESS EB LCs to find 50 triply eclipsing triples (∼0.005%). This means our new method is about an order of magnitude more productive in finding CHTs, and will be useful to identify further such objects when future Gaia data releases and more cataloged EBs from ongoing and future large surveys become available. A75, page 9 of 16 A&A 670, A75 (2023)