Eclipse-timing study of new hierarchical triple star candidates in the northern continuous viewing zone of TESS

Aims. We compiled a list of more than 3500 eclipsing binaries located in and near the northern continuous viewing zone (NCVZ) of the TESS space telescope that have su ﬃ cient TESS photometry to search for additional hidden components in these systems. In addition to discovering their hierarchical nature, we also determined their orbital parameters and analyzed their distributions. Methods. We obtained the TESS light curves of all targets in an automated way by applying convolution-aided di ﬀ erential photometry on the TESS full-frame images from all available sectors up to sector 60. Using a new self-developed Python GUI, we visually conﬁrmed all of these light curves, determined the eclipsing periods of the objects, and calculated their eclipse-timing variations (ETVs). The ETV curves were used in order to search for nonlinear variations that could be attributed to a light travel-time e ﬀ ect (LTTE) or dynamical perturbations caused by additional components in these systems. We preselected 351 such candidates and modeled their ETVs with the analytic formulae of pure LTTE or with a combination of LTTE and dynamical perturbations. Results. We were able to ﬁt a model solution for the ETVs of 135 hierarchical triple candidates, 10 systems of which were known from the literature, and the remaining 125 systems are new discoveries. These systems include some more noteworthy ones, such as ﬁve tight triples that are very close to their dynamical stability limit with a period ratio lower than 20, and three newly discovered triply eclipsing triples. We point out that dynamical perturbations occur in GZ Dra, which we found to be a triple, and that the system is one of the most strongly inclined systems known in the literature, with i m = 58 ◦ ± 7 ◦ . We also compared the distributions of some orbital parameters from our solutions with those from a previous Kepler sample. Finally, we veriﬁed the correlations between the available parameters for systems that have Gaia non-single star orbital solutions with those from our ETV solutions.


Introduction
In hierarchical triples both the inner binary and the outer, more distant third star are orbiting around the common center of mass of the system.They are especially interesting systems as their exact formation scenarios and their contributions to various evolutionary channels (e. g. stellar mergers or Type Ia supernovae) are still under debate (Toonen et al. 2020(Toonen et al. , 2022;;Tokovinin 2021).They are also the most likely systems in which we can detect and study dynamical perturbations which can cause variations in their orbital elements (Borkovits 2022).
There are several ways to discover hierarchical triple systems using photometric, spectroscopic and/or astrometric observations.Photometric methods can be applied only for eclipsing binaries (EBs) in which one can detect: i) the presence of an extra light in the light curve (LC) of the inner binary; ii) eclipse timing variations (ETVs) caused by light-travel time effect (LTTE); iii) eclipse depth variations (EDVs) caused by dynamically forced ⋆ Email: mtibor@titan.physx.u-szeged.huprecession of the orbital planes; iv) extra eclipses of the components of the inner binary with the tertiary.Spectroscopy can either reveal directly the triple nature of an object by resolving three sets of atomic lines in the spectra or indirectly by detecting periodic variations in the systemic radial velocity (gamma velocity) of a binary star.Lastly, astrometric observations can detect changes in the position of a binary star in the plane of the sky while orbiting the common center of mass of the triple system.Of all of these options, photometry is the most convenient technique because it can cover larger areas on the sky than spectroscopic studies and with orders of magnitude more objects within them, while one also usually needs larger telescopes and more telescope time for spectroscopy than for photometry in case of the same object.Astrometry is also a bit problematic because currently the Gaia space telescope is our only tool for the purpose of searching for hierarchical triples via this technique, and it just became such a tool very recently as Czavalinga et al. (2023) demonstrated after Gaia's 3rd Data Release (Gaia Collaboration et al. 2023;Babusiaux et al. 2023).That means, currently, photometric techniques are our best option to discover new hierarchical triple systems.Apart from the first listed effect (extra light) which can be unreliable on its own and tells us only the flux ratio of the supposed tertiary relative to the total flux of the system, the other three (ETVs, EDVs, extra eclipses, or their combination) can be excellent indicators.The analysis of ETV (or more traditionally, O-C) curves has a century long history, nevertheless, the first triply eclipsing systems were discovered only very recently, in the last decade (Carter et al. 2011;Derekas et al. 2011), due to the Kepler mission (Borucki et al. 2010) and, EDVs in such objects are observed routinely only since the advent of the era of space photometry.The Kepler space telescope, originally built to find transiting exoplanets, collected a 4 years long, almost uninterrupted photometric data set from a ∼100 square degree area of the sky with unprecedented precision that made it possible to detect the previously mentioned effects.As its successor, TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2015) is conducting a nearly allsky survey of exoplanets and its observations are also a treasure chest for almost all fields of stellar astronomy that should be exploited in many different ways, including searches for compact hierarchical triples.
There are a handful of ETV studies in the literature for large sky surveys such as CoRoT (Hajdu et al. 2017(Hajdu et al. , 2022b)), Kepler (Rappaport et al. 2013;Conroy et al. 2014;Borkovits et al. 2015Borkovits et al. , 2016) ) or OGLE (Hajdu et al. 2019(Hajdu et al. , 2022a) that found hundreds of such compact hierarchical triples (CHTs) collectively.This number can be improved further significantly by the analysis of the whole collection of the available and future TESS data.These systems could be of great help in searching for rare circumstances (e.g., super flat, super circular, super tight, super dynamically interacting, retrograde orbits, von Zeipel-Kozai-Lidov cycles in operation, etc), and fortuitous systems where we can make really precise measurements of the stars, their evolutionary state, and the orbits.
The TESS space telescope during its 2-year-long primary mission performed a unique sky survey collecting a highprecision, almost uninterrupted, ≳ 1 month long photometric data set for the majority of the objects located both in the southern and northern hemispheres.Because of the observational strategy of TESS, for the objects closer to the ecliptic poles, there are longer data sets available, as they can be involved in multiple sectors of observations with as long as a whole year for objects found in the Continuous Viewing Zones (CVZs).In its extended mission, TESS had re-observed almost the same parts of the sky extending the length of these data sets.These temporal samples give us an opportunity to analyze the eclipse times of EBs with relatively short orbital periods (up to around a few tens of days) and search for the signal of additional components with short outer orbital periods (up to around a few hundreds of days) in the ETVs.In this way, we can discover new CHTs and also reanalyze previously known systems in order to obtain more accurate parameters for them.
In this paper, we carry out a comprehensive analysis of the ETV curves of those EBs which were observed by the TESS space telescope in its Northern Continuous Viewing Zone (NCVZ) in between Sectors 14 and 60.As a close predecessor of this current study we refer to the work of Borkovits et al. (2016) who analyzed the four-year-long data set of the prime Keplermission to identify and study CHTs amongst more than 2600 EBs via their ETV curves.In this work we follow a similar (but not identical) method for the identification and ETV analyses of CHT candidates amongst TESS EBs.
The structure of the paper is the following.In Sect. 2 we describe briefly the methods and data sets to identify EBs in the TESS NCVZ sample, then the data reduction and the formation of the ETV curves are discussed.The basics of the mathematical modeling of the ETVs are summarized in Sect.3, while the results of our comprehensive, detailed analysis are discussed and tabulated in Sect. 4 and then, our work is concluded in Sect. 5.

Observational data and data reduction
We used a list of objects showing eclipse-like variations during the Year 2 observations of TESS found in the TESS data using a machine learning approach further described by Powell et al. (2023).From this list, we selected those objects that had at least 8 Sectors of observations available during Year 2, i. e. those that are located in or near the NCVZ of TESS.That yielded 5139 such targets which had long enough and almost uninterrupted data sets that could be useful for our survey.Moreover, we have added one more object (TIC 377105433) to our list from the list of Czavalinga et al. (2023) which, although located in the NCVZ, was too faint to be found by the search algorithm used to construct the GSFC EB Catalog.After filtering out duplicate sources, 3553 objects have made the final set of targets for further analysis.
We obtained all TESS light curves of the objects in our sample up to Sector 60 in an automated way using a convolutionaided differential photometric pipeline based on the FITSH software package (Pál 2012).We removed the remaining trends using the WOTAN (Hippke et al. 2019) Python package or using the Principal Component Analysis method implemented in the lightkurve (Lightkurve Collaboration et al. 2018) Python package in some more complicated cases.In order to determine the eclipsing periods of the systems, first we ran an initial Box Least Squares (BLS; Kovács et al. 2002) search on the full light curves.Then, using the best BLS period, we transformed the light curves to the orbital phase domain and visually checked if the initial period found by the BLS is correct.As a final step, we used the Phase Dispersion Minimization (PDM; Stellingwerf 1978) period search method to obtain more accurate eclipsing periods.
After we determined the periods of the systems, as a next step, we calculated eclipse timing variations in a similar way as (Borkovits et al. 2015) using the following equation: where ∆ is the time difference of the observed and calculated mid-eclipse times (O-C) for the given orbital cycle, T (E) is the mid-eclipse time of the Eth orbital cycle, T 0 is the mid-eclipse time of the zeroth orbital cycle (reference epoch), P s is the reference period and E is the cycle number.Finally, we visually vetted all ETV curves in order to select those which show clear non-linear variations for further analysis.
For this purpose, we developed an interactive program with a graphical user interface (GUI) utilizing the Tkinter module of Python which can be seen in Fig. 1.This program allows the user to quickly load and view any raw LC with several useful interactive features that help to analyze the LC and quickly calculate the ETV curve of an object using the algorithm described above.The already implemented features are: i) interactive detrending of the raw LC; ii) multiple period searching methods, e.g.BLS Periodogram and PDM in order to determine the eclipsing period of the system automatically; iii) reference epoch determination with the calculation of the (folded) orbital phase curve; iv) determination of the times of eclipsing minima or out-ofeclipse maxima with the fitting of higher-order polynomials; v) displaying the ETV curve of the object calculated from the previously determined minima/maxima using the previously found orbital period and reference epoch; vi) save the latest results in a database that allows one to do more complex external analyses on the candidate systems, and also ensure reproducibility.

Theoretical basics of the ETV analysis
The ETVs were modeled in a very similar manner as described in Borkovits et al. (2015Borkovits et al. ( , 2016)).The full analytic model takes the form: where the polynomial term contains corrections for the calculated eclipse times (in the case of i = 0, 1), constant-rate and, linear period variations (in the case of i = 2, 3).The constant and linear (i = 0, 1) terms were used for all the sources in our analyses.In contrast to this, the i = 2, 3 polynomial terms were fitted only in a few cases, where the presence of a quadratic or even more complicated additional variations in the ETV curves were evident.Moreover, the (∆ LTTE,dyn,apse ) terms represent the contributions of LTTE, P 2 -timescale dynamical perturbations and apsidal motion, respectively, while the square bracket denotes the difference of each terms measured in the Eth and zeroth (i.e. reference) cycle.Here, for consistency, we give only very brief descriptions and the formulae of these latter terms as they were used, while their detailed explanation can be found in Borkovits et al. (2015Borkovits et al. ( , 2016)).

Light-travel time effect (LTTE)
The classical LTTE (or Rømer-delay) contribution comes from the periodically varying distance of the inner binary during its orbital revolution around the common center of mass of the triple system that causes a periodic variation in the light travel time towards the observer detected as periodic eclipse timing variation.Based on Irwin (1952), it was taken into account in the form of where i 2 , e 2 , v 2 and ω 2 refer to the inclination, eccentricity, true anomaly and the argument of the periastron of the outer orbit1 (i.e. the relative orbit, on which the tertiary star orbits around the center of mass of the inner binary), respectively, while a AB stands for the semi-major axis of the (physical) orbit of the center of mass of the inner pair around the center of mass of the whole triple system.Finally, c denotes the speed of light.
It can be shown easily that, as far as the orbital elements of the outer orbit remain constant, this effect can be described with a pure sine in the eccentric anomaly domain, of which the amplitude is By analogy with the mass function used for single-lined spectroscopic binaries (SB1), it is usual to introduce the mass function as where, in the last form, the numerical constant is correct when time is expressed in days and masses in solar mass.(Note, here and hereafter, m A,B,C denote the individual masses of the three constituent stars, while m ABC stands for the total mass of the triple system.) In regard to the wide, or outer orbit, an LTTE solution carries exactly the same information that can be extracted from the radial velocity measurements of a single-line spectroscopic binary (SB1).Therefore, besides the orbital period of the outer orbit (P 2 ) one can determine three of the six orbital elements of the outer orbit (e 2 , ω 2 , τ 2 ) and, moreover, the projected semi-major axis of the inner pair's orbit around the common centre of mass of the triple (a AB sin i 2 ) from which, in the absence of any external information, one can determine the mass of the third companion as a unique function of two additional parameters, namely, the total mass of the inner binary m AB and the orbital inclination (i 2 ).Such external information might have come from the observations of radial velocities of the third star (producing a quasi SB2-like system) and/or from detection of third-body eclipses (a good example for both situations is e. g., HD 181068, Borkovits et al. 2013) or from the detection of dynamically driven ETVs in tight triple systems.Note, when such external information was not available for the total mass of the inner pair, in the case of W UMa-type, i. e., overcontact binaries, we estimated the binary's mass with the use of the empirical mass-period relations of Gazeas & Ste ¸pień (2008), while for other systems we simply set m AB = 2 M ⊙ .

P 2 -timescale dynamical perturbations
The contribution of dynamical perturbations comes from the possible gravitational interaction between the outer component and the inner binary which can result in the change of the orbital parameters of the inner binary and hence period variation.While LTTE terms were fitted by default during our analysis, the P 2timescale dynamical contributions were taken into account only when the tightness of the system under investigation made it necessary.The most elaborate presentation of the theory of the P 2 − (medium-) timescale perturbations and their manifestations in the ETV of a tight EB was presented and discussed in Borkovits et al. (2015).For the sake of brevity, here we repeat only the lowest order (i.e.most significant) quadrupole terms (modifying slightly Eqs.5-11 of Borkovits et al. 2015), while higher order (and more lengthy) formulae can be found in the appendices of Borkovits et al. (2015).According to this theory, the quadrupolelevel dynamical contribution to an ETV curve in a CHT is the Article number, page 3 of 28 A&A proofs: manuscript no.TESS_NCVZ_clean Fig. 1.Current look of our self-developed, unique light and ETV curve analyzer GUI.On the top left corner, the raw lightcurve can be seen (blue) alongside the fitted trends (orange).On the top right corner, the detrended lightcurve is plotted after the normalization by the previously fitted trends.On the lower left corner, the previously detrended and normalized lightcurve binned and folded with the specified bin numbers and ephemeris is plotted alongside the polinomial templates fitted on the two types of eclipses.Finally, on the lower right corner, the ETVs calculated from the primary (red) and secondary (blue) mid-eclipse times using the indicated ephemeris on the top are plotted.Below each panel there are several useful features (see the text for details) which enable a fast interactive analysis of any set of lightcurves. following: where stands for the characteristic amplitude of the dynamical effects on a timescale equal to the outer period (P 2 ).Furthermore, I = cos i m denotes the cosine of the mutual or relative inclination (i m ), i. e., the angle between the inner and outer orbital planes, while u 2 = v 2 + ω 2 stands for the true longitude of the third star measured from the ascending node of the wide orbit with respect to the tangential plane of the sky.Moreover, the node-like quantities n 1,2 represent the angular distances of the intersection of the two orbital planes (i.e. the dynamical node) from the intersections of the respective orbital planes with the tangential plane of the sky (i.e. the observable nodes), while α = n 2 − n 1 and β = n 2 + n 1 (See Fig. 1 of Borkovits et al. 2015, for an easy visualization).Furthermore, the trigonometric terms, depending on the parameters of the outer orbit (and the current position of the third star) are as follows: and, finally, ∆ * 1 (sin i m cot i 1 ) stands for the usually negligibly small precession terms which are not shown here, however, they are taken into account during the analytic fitting process.(Note, that the upper signs refer to the primary eclipses, while the lower ones for the secondary events.) As one can see, in the case of a circular inner orbit (which occurs frequently for close binary stars and, hence, for most of the EBs), all the K functions disappear and, moreover, f 1 reduces to unity, resulting in a much simpler equation for ∆ dyn .
The ∆ dyn term was included into the ETV analysis only when the tightness (i.e. the ratio P 2 /P 1 ) of the triple system made it necessary.In several cases, typically for the tightest, eccentric and/or inclined triples, the shapes of the ETV curves have shown already at a first inspection that the dynamical effects are important/dominant.For some other cases, however, the need to include dynamical effects was far from evident.For this reason, in the case of all the LTTE-only fitted ETVs the software automatically calculated the nominal amplitude ratio of A dyn /A LTTE under the assumption that m AB = 2 (discussed above) and i 2 = 90 • , the latter of which yields the minimum value of the amplitude ratio.If this ratio was found to be higher than 0.25, we repeated the fitting process, switching on the dynamical contribution.
The dynamical ETV terms, in theory, give much more orbital, geometric and other dynamical parameters, than does a pure LTTE-dominated ETV solution.This is so because the dynamical terms depend at the same time not only upon the orientations of both orbits relative to the observer, but also upon the orientations of the two orbits relative to each other (as the gravitational perturbations are driven by the relative positions of the three stars with respect to each other) and, hence, the dynamical terms give the full, three-dimensional, spatial configuration of the triple system, including such a key-parameter as the mutual (or, relative) inclination (i m ) of the inner and outer orbits.Moreover, in the amplitude of the dynamical ETV the outer mass ratio (q 2 = m C /m AB ) also occurs, and moreover, the outer inclination i 2 can also be obtained.Thus, with the combination of the LTTE and dynamical terms, m C and m AB can also be determined directly.Unfortunately, however, the parameter space of a combined LTTE and dynamical ETV solution is so complex and multiply degenerate, that in the absence of any other, external information, one cannot expect accurate mass determination from such a combination of contributions to the ETV.A detailed discussion of the parameters which can be mined out from a dynamical LTTE solution can be found in Borkovits et al. (2015), while the degeneracies of such combined ETV solutions are discussed in Sect.6.3 of Rappaport et al. (2013) and Sect.5.2 of Borkovits et al. (2015).

Apsidal motion
Apsidal motion (i.e. the rotation of the orbital ellipse within the orbital plane) may have several different origins.The most frequent and well known variations are the relativistic, classical tidal (arising from the non-spherical mass distribution of the tidally deformed bodies) and, the dynamical (third-body perturbed) apsidal motions (see, e.g., Borkovits et al. 2019a, for a short review).Irrespective on the origin of this phenomena, its effect on the ETV curves can be described mathematically as follows: (The form of this equation in widespread use is given as a trigonometric series in ω 1 .The slight inclination angle dependence is also taken into account in Gimenez & Garcia-Pelayo 1983.)In general, for tight triples and, especially, for such triples with less compact inner pairs, the dynamically forced apsidal motion substantially dominates over the other two contributions.We will show in Sect.4, that for all of the currently investigated triples with eccentric inner binaries, this is the actual scenario.Hence, instead of fitting a (linear) apsidal advance rate (∆ω 1 ), hidden in Eq. 15, as an additional free parameter, we calculated the current values of ω 1 (and also of e 1 ) inherently for each time during the fitting process, with the use of the perturbation equations, and in the manner, as described in Borkovits et al. (2015).

Methods of the ETV analysis
The software package which realizes the mathematical model discussed above is described in detail in Borkovits et al. (2015), and it has been applied previously in several studies analyzing ETV curves of several Kepler (Borkovits et al. 2015(Borkovits et al. , 2016)), K2 (Borkovits et al. 2019b), CoRoT (Hajdu et al. 2017), OGLE (Hajdu et al. 2019) and TESS (Borkovits et al. 2021;Pribulla et al. 2023) EBs.For those systems where there was no need to include the medium-period dynamical ETV contribution (Sect.3.1.2)we applied a simple Levenberg-Marquardt-type (LM) differential least squares parameter optimization to find the most probable solution.(The pros and cons of the use of this method for the current problem, and some further details are discussed in Sect.5.2 of Borkovits et al. 2015.)When the inclusion of the dynamical contribution was necessary, we followed a different parameter optimization technique.Because of the large number of parameters to be fitted, as well as their highly non-linear correlations, this is not an ideal situation for the LM method.Therefore, we switched to Markov-Chain Monte Carlo (MCMC) approach.This is a novelty of the software package relative to its formerly used versions.In such a manner we were able to explore the whole (or, at least, physically realistic) part of the parameters and obtain more realistic uncertainties.

Results and discussion
As was mentioned in Sect. 2 in total we identified 351 EBs of which the ETVs exhibited promising non-linear variations.In order to identify already known triple systems or candidates in Article number, page 5 of 28 Fig. 2. Locations of the analyzed hierarchical triple star candidates together with other triples, detected with space telescopes Kepler and TESS in the P 1 − P 2 plane.Black and red dots and green squares represent the current NCVZ candidates.For the black and red systems simple LTTE solutions were satisfactory, while for the green ones the dynamical effects were also taken into account.(The red symbols, however, denote the most uncertain solutions which were omitted from the statistical investigations.)Gray dots and squares denote the same classes of those hierarchical triple star candidates in the original Kepler-field, which were analyzed in Borkovits et al. (2016).Gray triangles represent recently discovered K2 and TESS systems with accurately known photodynamical solutions.(These are mostly triply eclipsing triple systems.)The vertical red line at the left shows the lower limit of the period of overcontact binaries.The horizontal and sloped blue lines are boundaries that roughly separate detectable ETVs from the undetectable ones.The detection limits again, were set to 50 sec.These amplitudes were calculated following the same assumptions as in Fig. 8 of Borkovits et al. (2016) The arrows indicate the direction of the increasing of the respective amplitudes.The shaded regions from left to right represent (i) the W UMa desert, i.e. the (almost) empty domain where a tight third companion of a short-period EB would be certainly detectable through its LTTE, even in the absence of a measurable dynamical delay; (ii) the purely dynamical region, i. e., where the dynamical effect should be detectable, while the LTTE would not and; (iii) the dynamically unstable region in the sense of the Mardling & Aarseth (2001) formula.Note, that while the border of this latter shaded area was also calculated with e 2 = 0.35, we also give the limit for e 2 = 0.1, as the thinner red line within this (light blue) region.
our sample, we thoroughly checked the literature of these 351 EBs.We found five formerly known EBs, where previous ETV studies (based on former, ground-based observed eclipse times) indicate the presence of an additional, more distant component.These are TIC 219109908 (EF Dra, Yang 2012), TIC 219738202 (BX Dra, Park et al. 2013), TIC 233532554 (RR Dra, Wang & Zhu 2021;Şenavcı et al. 2022), TIC 298734307 (HL Dra, Shi et al. 2021) and TIC 288734990 (NSVS 01286630, Wolf et al. 2016;Zhang et al. 2018).In the first four cases, however, the previously found third-body periods are in the range of decades and, hence, the currently found, short-term non-linear ETV behaviors cannot be accounted for by those putative companions, while in the last system although the outer orbital period is on the same order as in our solution, the TESS ETV data and, hence, our solution are not in accord with the former results.For 9 other objects in our candidate sample we found some external archival eclipse times, but there were no actual earlier ETV studies of these objects.
In conclusion, out of the promising 351 EBs, we identified 135 hierarchical triple system candidates, of which the CHT nature looks almost certain, and another 77 systems which are at least moderately likely CHT candidates.The remaining 58 EBs are less likely cases.In Fig. 2 we plot the locations of these 77 likely CHT candidates in the P 1 vs P 2 plane where, for a comparison, we also plot some previously known CHTs from some other surveys.Moreover, we plot the ETV curves together with Article number, page 6 of 28 the best analytic LTTE or, combined LTTE+dynamical fits for all the 135 systems in Appendix A.
Here, in what follows, first we present our findings both for the pure LTTE and the combined LTTE+dynamical solution systems in general, and then we discuss some specific subgroups of the newly discovered CHT candidates.

Triple system candidates with pure LTTE solutions
From the ETVs of the narrowed sample of the 135 good CHT candidate systems, we identified 94 EBs in which case we were able to find pure LTTE solutions.According to the reliability of the solutions we divided these systems into three groups.Into the first group we ranked those 19 EBs where we judge our LTTE solutions to be certain or, at least highly likely.These are typically such close likely triple systems, where the period of the third component (i.e. the outer period) does not exceed half of the nominal ∼ 1278-day length of the complete TESS NCVZ data set, the amplitude of the ETV curve(s) is substantially higher than the scatter of the individual ETV points and, finally, all the sensitive sections of the light-time orbits are well covered.Note, while the first two criteria look trivial, we had to introduce the third one because of the large intermediate gaps in the observations between Sectors 26 and 40, and Sectors 41 and 47.An example of those CHT candidates where this last criterion was applied is TIC 159398028.This system, despite its outer period of P 2 = 671 ± 3 d, was ranked only in the second (less certain) group because none of the lower ETV extrema -and only one of the two upper extrema -was covered by the observations.The outer P 2 periods of the third companions are in the range of 109 and 628 days.All these compact triples have very close inner EBs (there are only three inner pairs with P 1 > 1 days).This is naturally not an unexpected fact, since for longer inner periods such compact systems would have dynamically dominated ETVs.We tabulate the main parameters of the LTTE solutions of this group in Table 1.The uncertainties of each parameters are given in parentheses and refer to the last digit(s).As the pure LTTE solutions were calculated with non-linear LM optimization, these are formal uncertainties obtained via matrix-inversion and the resultant correlation matrices after reducing the χ 2 -s to unity.In the last three columns we give the minimum third-body mass (m C ) min , the nominal dynamical to LTTE amplitude ratio, A dyn A LTTE , calculated with the use of this minimum mass and, finally, the assumed mass of the binary, m AB , which was used for both the calculation of the minimum of the outer mass, and the amplitude ratio.When the total mass of the inner binary is known from the literature, then that value is given together with the current references, but in most cases we give only some realistic estimate (which we denote with : after the numerical value).
We classified 17 additional systems with less certain but likely LTTE solutions.These systems are CHTs with a high certainty, but do not fulfill all the criteria which were required for the first group.Their outer periods are in between 507 and 1011 days, i.e. more than one full orbit is covered at least partly, in all cases.All but four systems have inner periods less than 1 days.These are listed in Table 2.
Finally, in Table 3 we tabulate the parameters of the LTTE 'solutions' of the remaining 58 systems.This is a very inhomogeneous group.In the case of some EBs the LTTE solution looks quite reasonable, however, when the period is certainly longer than the duration of the dataset, the parameters we obtained are likely not very robust.In other cases, however, even the interpretation of the non-linear ETVs as LTTE signals seem suspect.
Strictly speaking, from the fact, that we were able to obtain any LTTE solutions do not follow that there is any third component in the given systems.Therefore, this third set should be considered only with great caution, and we do not recommend using the numerical results of this subgroup for any statistical studies.Instead, the group of these systems should be considered only as EBs which display clearly non-linear ETVs and hence, are worthy of further attention.

Triple system candidates with combined LTTE and dynamical solutions
For 41 EBs we obtained combined LTTE and dynamical ETV solutions.The inferred outer periods in this group range from 68 days to ∼ 4000 days.Similar to the pure LTTE systems, again, we consider the solutions of those 27 triples certain where the outer period does not exceed half the length of the given data set (see Table 4).(Naturally, the same caution that was discussed in case of the LTTE systems, must be taken into account here as well.)The nominally less certain solutions of the remaining 14 systems are tabulated in Table 5.The structure of these Tables 4, 5 departs slightly from that of the three tables of LTTE triples (Tables 1-3).This is in accord with the fact that from a combined solution one can determine individual masses (at least m C and m AB ), and, hence the semi-major axis of the outer orbit (a 2 ) can also be readily calculated.Moreover, in these tables instead of a theoretically calculated, nominal dynamical to LTTE contribution amplitude ratio, we give the true, 'measured' ratio of A dyn /A LTTE , which can be determined directly from the analytic ETV solution curves.Because of the large number of the parameters that can be obtained from such a dynamical solution, we also introduce Table 6 where we give additional parameters for the combined solution systems.These are mainly the apsidal motion and/or orientation and nodal precession parameters of these triples.In this regard, however, some caution should be taken as follows.
First, as was discussed in Sect.5.2 of Borkovits et al. (2016) in the context of the mutual inclination (i m ) parameter, practically one cannot distinguish between, e.g. an i m = 0 • and i m = 10 • ETV solution, at least not on the timescale of TESS observations.From an observational point of view, however, there may be strong differences between the two solutions, as the second one may result in rapid inner inclination (i 1 ) and, hence eclipse depth variations and, therefore, they should be clearly distinguishable.Hence, in all cases, when our solution resulted in such a mutual inclination value, which was not in accord with the rate of the observed EDVs (including non-EDV cases), we did not accept that solution, and forced a coplanar (i m = 0 • ) solution.
Second, the tabulated inner inclinations (i 1 ) should also be considered with great caution.The dynamical ETV solution, especially for coplanar cases, is almost insensitive to this parameter.Hence, it is better to fix its value a priori.Since for the vast majority of the combined LTTE+dyn systems no light curve solution is available, generally i 1 is an unknown parameter.Fortunately, however, most of the dynamically dominated systems have longer inner periods (see below) and, hence, one should not get large errors by assuming inclinations of i 1 ∼ 87  TIC No.  Consequently, the pure LTTE triple star candidates should be offset towards the shorter inner period, i.e., closest systems.

Systems with special interests
In what follows, we discuss some small subgroups of the identified triple star candidates which are worthy of special interest for different reasons.

Tightest triples
We found 5 systems with P 2 /P 1 ≲ 20.These are TICs 236774836 (with P 2 /P 1 ∼ 9.35), 219885468 (14.7), 389966039 (16.3), 235934882 (18.6), 233684019 (19.9).(In Fig. 2 these systems are represented with the five green rectangles which are the closest to the border of the dynamically unstable domain at the lower right part of the plot.)Despite the weakly hierarchical nature of these triples, our analytic, hierarchical three-body models fit and describe these ETVs fairly well.However, according to our prior experiences with Kepler systems, for such tight triples, our analytic hierarchical triple star model becomes less accurate on a somewhat longer timescale of a decade (which, in general, is about the timescale of the "apse-node" type or, long-period perturbations of these triples), suggesting that the analytic model of this latter class should be improved).In this regard, TIC 219885468 serves as a good case study of that effect, because recently Borkovits & Mitnyan (2023) had carried out a more complex photodynamical analysis of this triple star system.
During such an analysis the motion of the triple star system is integrated numerically and, hence, we can compare our analytical results with the outputs of the more exact numerical integration.Regarding the masses, Borkovits & Mitnyan (2023)  = 83 • ± 8 • ).However, this difference is close to 180 • , and can be understood from that point of view, that the analytic perturbation formulae depend upon only 2ω out and, hence, they are insensitive to a 180 • difference in ω 2 .On the other hand, the LTTE-term depends on ω 2 , which fact, in general, could have resolved this ambiguity, but in the current situation the LTTE contribution is only about 2 − 3% of the dynamical one in the ETV and, thus, our analytical solution still suffers from this ambiguity.In conclusion, we can say that our approximate, analytical ETV solutions look quite acceptable even for such tight triples.On the other hand, however, turning to the some decade-long evolution (i.e. apsidal motion) of the current (and also, other, similarly tight) systems, we again refer to Borkovits & Mitnyan (2023) who pointed out, that there are some discrepancies between the analytically predicted and numerical integrated apsidal motion periods (for the current system P analytical apse,1 = 20 yr versus P photdyn apse,1 = 13.6 yr).Hence, future follow-up observations and then numerical integration and modeling of the triple star motion of these tight systems would be especially important.

Triple star candidates with eclipse depth variations (EDVs)
We identified eight triple star candidates where the inner EBs exhibit EDVs.These systems are: TICs 233684019, 233729038,236774836,259004910,288611133,357686232,417701432,441768471.Amongst them, the most dramatic variation can be found in TIC 236774836, which is the tightest triple system candidate in our entire sample (with P 2 /P 1 ∼ 9.35).While during the Year 2 observations the eclipse depths int his source remain constant, a dramatic decrease in the eclipse depths takes place in the Year 4 data and, the eclipses disappear completely near the beginning of Year 5 (Sector 56).The most likely explanation for such a rapid EDV is a non-coplanar configuration of the given triple system which forces precession of the orbital planes and, hence, variation of the inclination of the EB.As the timescale of such an effect is proportional to P2 2 /P 1 , one can expect the detection of such a rapid effect in the tightest and most compact triple systems.Finally, note that TIC 230002837 also exhibits slight, suspicious EDVs, but we were unable to find any reliable (or, at least, uncertain) third-body solution for its ETVs.

Eccentric EBs (eEBs) with detectable apsidal motion
We found eccentric EB orbits in 27 candidate systems (see Table 6).All of them belong to the systems with combined LTTE+dynamical ETV solution.In an eEB one can expect apsidal motion.While the eccentricity of an EB's orbit can be discerned, in general, simply by a visual inspection of a possibly unequal time interval between the odd-even and even-odd eclipse events and/or the different durations of the even and odd eclipses, the apsidal motion most readily reveals itself via the divergence or convergence of the primary and secondary ETV curves (which is the manifestation of the anticorrelated behaviour of the primary and secondary ETV curves on a time interval, that is substantially shorter than the period of the apsidal precession).Because the duration of the nominally ∼ 1300 d TESS observations is very short relative to the decades to millennia long tidal and relativistic apsidal motion periods, one cannot expect the detection of apsidal motion caused primarily by these effects in our ETV curves.Thus, where we clearly detect apsidal motion in our short-term ETVs, we can be certain that this effect is due to rapid, dynamically forced apsidal motion.
We identified clear apsidal motion in 10 systems.These are TICs 199616648, 219771659, 219885468, 233684019, 233729038, 233738966, 236774836, 256514937, 357686232, 417701432.Note, there are some other systems, where the ETVs suggest apsidal motion but the detection is less certain since the P 2 -period dynamically forced ETVs may hide the convergence or divergence of the two ETVs.These 7 systems are TICs 236769201, 243337122, 259271740, 272705788, 288611133, 376976908, 441768471.In Table 6 we give the theoretical dynamical apsidal motion periods (P apse ) which are calculated by our software using the quadrupole order perturbation theory of the hierarchical three-body problem (see Borkovits et al. 2015, Appendix C, for a detailed description).As one can see, the shortest theoretical apsidal motion period is P apse = 9.6 yr.It was detected in TIC 236774836, i. e., in that notable system where the eclipses completely disappeared for 2022 (see in Sect.4.3.2above).
Note also, that our software package calculates the theoretical, dynamical apsidal motion for all non-circular, dynamically dominated ETVs and, hence, P apse is tabulated in Table 6 for all eEBs, irrespective of whether the apsidal motion manifests itself during the short interval of TESS observations.Hence, one can realize that there are four systems, for which the given P apse value is negative.Here, the negative sign denotes retrograde apsidal motion.While the tidal and relativistic apsidal motions are always prograde (i.e., the orbital ellipse revolves in the same direction as that of the orbiting bodies), dynamically forced apsidal motion can be retrograde either for dynamical or geometric reasons (see, e. g.Borkovits et al. 2022).

Triply eclipsing triples
We found third-body eclipses in four systems.These are TICs 229785001, 233684019, 288430645, 441738417.Their third-body periods range from 145 to 472 days.The periods deduced from the ETV curves are in agreement with the locations of the extra eclipses in all the four systems, confirming that the same objects are the origin of the extra eclipses and the timing variations.Note, a detailed photodynamical analysis of TIC 229785001 was published earlier by Rappaport et al. (2023).The results of our analytic ETV solution are in good accord with their results.We plan complex photodynamical investigations of the remaining three systems in the near future.

Interpretation of the peculiar ETV curve of GZ Draconis
While GZ Draconis does not appear to be an extremely interesting, peculiar triple system; nonetheless we feel it worthwhile to discuss the cautionary tale of this quite average EB.This relatively bright (V=9.4)object first was known as a visual double star from the 1980's (see, e. g.Muller 1984;Heintz 1987) and then, the EB nature of the brighter component was discovered by the Hipparcos satellite (Kazarovets et al. 1999;Fabricius & Makarov 2000).The system was observed continuously in 2-min cadence mode during Years 2, 4 and 5 of TESS observations.From our point of view, only the ETV curves calculated from the TESS-observed primary and secondary eclipses have importance.It was clear for us already after the Year 2 observations that the timing diagrams show clear, third-body induced ETVs.By the end of the Year 2 observations, and using some former, less accurate times of primary eclipses, calculated from the publicly available SuperWASP observations, we were able to find reasonably good LTTE solution with a period of P 2 = 476 ± 1 d.We realized, however, that there was an unexpected dip on the primary ETV curve around BJD ∼ 2 458800 with an amplitude of ∼ 0 d .0003and duration of ∼ 40 days, but did not pay careful attention to this feature, assuming that it was caused by some light curve distortions (for example, spottedness).Note, this dip cannot be seen on the much less accurate secondary ETV curve.
Then the new NCVZ monitoring of TESS began during the second year of the first extended mission, and the new eclipse times confirmed the former ETV solution sectors by sectors.The uninteresting characteristics of such a straightforward ETV curve changed immediately when the Sector 53 observations became available, and we calculated the new set of eclipse times.The same dip in the ETV curve as in the Year 2 data occurred again and, at the very same phase of the third-body orbit as two outer orbital cycles earlier (Fig. 3).By this time we realized that even though the P2/P1 ratio of ∼211 is large (i.e., not a very tight triple) the dynamical effect might still actually produce a nonnegligible contribution to the ETV curve and explain the unexpected "extra" dips. 3A further, more in-depth consideration of the shape of the ETV curve makes this possibility even more exciting.Taking a look at on the other (bottom) extrema of the ETV curves, one realizes that during these phases both in the Year 2 (in between cca.JDs 2 458 930 and 2 459 000) and the Year 5 (JDs 2 459 880 -2 459 950) the curvatures of the ETVs practically disappear, which suggests the presence of a second bump, which counteracts the LTTE.Therefore, we may assume two bumps within one outer period (for the dynamical curve) which, in the case of a circular inner orbit (as is the case in GZ Dra), may occur only if the two orbital planes are substantially inclined to each other (see Eq. 46 of Borkovits et al. 2003).Thus, the shape of the ETV suggests that there is an exceptional chance in the case of GZ Draconis for determining the mutual inclination purely from the ETV for such a non-tight triple star system.Then, our search for an inclined, combined LTTE and dynamical ETV solution became quite fruitful, and we realized that this triple star must be one of the currently most inclined triple star systems known with i m = 58 • ± 7 • .Finally, in this regard, we make two additional notes.First, we refer to the fact that that this mutual inclination is well within the high mutual inclination domain of the original von Zeipel-Lidov-Kozai phenomena (see, e.g.Naoz 2016, and further references therein), where large amplitude eccentricity cycles may occur.Despite this fact, one cannot expect such large amplitude eccentricity variations in the currently circular inner pair, because of the tidal interactions between the two components, which cancels the ZLK cycles.
Second, such a large mutual inclination implies very large amplitude orbital plane precession cycles, which would lead to substantial EDVs and also the disappearance (and later, reappearance) of the eclipses.While this may be true, as one can see in Table 6, the theoretical period of such variations in GZ Draconis (with the measured parameters) is about P node ∼ 3 800 yr, which explains the absence of EDV during the few-year-long TESS observations.et al. (2016) published the largest collection of CHTs with known orbital periods using a sample of 222 objects found in the Kepler field.We have about a third of that sample size with 77 objects that are likely or very likely to be triples, but since we applied similar detection and analysis techniques, it may be worthwhile to perform a similar statistical analysis on our newly found systems in the TESS NCVZ and make a comparison.

Borkovits
In Fig. 4, we show the distribution of the outer orbital periods in our sample.It can be seen that most systems have outer orbital periods of 100-800 days.In our photometric data sets there are smaller and larger gaps due to the observing strategy of TESS: the former are caused by the fact that the telescope changes its observational field of view in every ∼27 days and it is not guaranteed that the object will continue to be observed in Table 4. Orbital Elements from combined dynamical and LTTE solutions with outer period shorter than half of the length of the data sets.
In Fig. 5., we plot the distribution of the outer eccentricities of systems in our sample.Except six objects, all of them have an eccentricity value smaller than 0.55 and the distribution looks more or less flat with some small fluctuations.It does not have similar characteristics to those of the distribution from the Kepler sample.Nevertheless, if we compare the cumulative eccentricity distribution (Fig. 6.) with those from earlier studies, they are very similar to each other4 and completely different from simple linear or thermal distributions from theory.It resembles more the distribution that of found by Duchêne & Kraus (2013) for populations of field EBs.Fig. 7. shows the relation between the outer orbital period and the outer eccentricity for both studies.The black line denotes a fit for a linear relation to the TESS NCVZ sample with  a correlation of 0.05 (i.e., not very correlated) which is significantly lower than the value of 0.34 found by Borkovits et al. (2016) for the Kepler sample.We emphasize again, however, that our TESS sample is much smaller than the Kepler sample of Borkovits et al. (2016) and, this is more expressively true for the longer outer period systems, where one would expect higher outer eccentricities, as they cannot be detected certainly in the current stage of the TESS sample.
In Fig. 8., we plot the distribution of the tertiary component masses (M C ) for both our and the Kepler sample.For systems with measurable dynamical perturbations, exact masses can be determined, nevertheless for systems with LTTE only minimum masses can be estimated from the mass function using the constraints on the inner binary mass explained in Sect.3.1.1.Similarly to the Kepler sample, the majority of the systems with measured masses have tertiaries with less than 1.8 M ⊙ .
Fig. 9. shows the correlation between the exact masses of the binary and the tertiary components for those systems in which measurable dynamical perturbations occur for both our TESS NCVZ and the Kepler sample.The orange and blue lines represent the cases where all three components have the same mass, Article number, page 13 of 28  (3) 84(60) 58683.0(1)  3(1) 87.7(6) 89.23(4) 3(1)  230012179 8.9235427(4) 25(2) 0.031(5) 91.0(3) 58670.490(9) 5103  0   89  0  −3880  233073872 d 14.142051(9) 41(1) 0.039(2) 139(3) 58665.5(1and where the mass of the tertiary equals the mass of the inner binary, respectively.The majority of the systems are located under or close to the orange line, while only 4 systems are above the blue line indicating that the tertiary has a higher mass than the inner binary components.It may also be worthwhile pointing out that in the Kepler sample, the majority of the systems have an inner binary mass of less than 2 M ⊙ , while in our TESS NCVZ sample, there are only a handful of systems with such inner bi-nary masses and the majority of the systems are above this value.Hence these systems may serve as an extension to the regime of higher binary masses.On Fig. 10., we plot the distribution of the mutual inclination for systems in which the measurable dynamical perturbations allowed us to determine this quantity for both samples.The vast majority of the systems in our TESS NCVZ sample have a (nearly) coplanar configuration and only a few systems can Article number, page 14 of 28     be found with higher mutual inclination values.This is different from the Kepler results which showed a small peak centered at i mut = 40 • .However, in both cases the results are somewhat hampered by small number statistics.
Finally, in Fig. 11., we show the distribution of apsidal motion timescales for both samples.

Comparison with Gaia NSS candidates
As it was recently demonstrated in Czavalinga et al. (2023), it is also possible to identify new hierarchical triple candidates using Gaia data comparing the orbital periods from the astrometric or spectroscopic Gaia NSS orbital solutions with the eclipsing periods from lightcurves and finding a difference of at least 5 times.For this purpose, we also crossmatched our full sample of 5139 objects with the Gaia DR3 Non-Single Star catalog (Gaia Collaboration 2022).That yielded 32 and 116 matches with Gaia astrometric and spectroscopic orbital solutions, respectively.After the comparison of the orbital periods of the Gaia orbital solutions with the binary periods from the TESS light curves, we excluded 8 and 115 objects as the two types of periods were identical or exactly integer or half integer multiples of each other (i.e Gaia picked up the EB period or its possible aliases because of limited Gaia data).That left a total of 25 triple candidates in our sample with available Gaia orbital solution for the outer orbit of a possible tertiary component.As a final step, we analyzed the TESS light and ETV curves of these systems one-by-one and excluded 10 further systems for various other reasons (e.g. the source was not an eclipsing binary, or the light curve and the Gaia NSS solution actually belonged to another close star without a Gaia NSS solution).After all this, in total we found 15 systems that have a Gaia NSS orbital solution, and, moreover, the ETVs also confirm the validity of these solutions for all of them.For these objects, we tabulated the orbital parameters from the two types of orbital solutions in Table 7. Five of these objects (TICs 229762991, 232606864, 237234024, 320324245, 377105433) were already found and published recently by Czavalinga et al. (2023) and our ETV solutions for these triples are in a good agreement with theirs, while the other 10 systems are new discoveries that also demonstrate the potential of the Gaia NSS catalog in the field of discovering new hierarchical triples.In Fig. 12. we plotted the correlations of different orbital parameters determined from the two types of orbital solutions (Gaia NSS vs. ETVs) for the 22 systems found by Czavalinga et al. (2023) and the 10 systems newly found in our NCVZ analysis in a similar way as Czavalinga et al. (2023).From these figures, it can be seen that the 10 new systems found by us can further confirm the conclusions made by Czavalinga et al. (2023), i. e. Gaia does a good job picking up the correct outer orbital period, but the ETV results are often more accurate, especially in case of the outer eccentricity, projected semi-major axis and argument of periastron.

Summary
In the hope of discovering new compact hierarchical triples, we have conducted a survey of eclipsing binaries based on TESS data.We have selected 3553 objects from a list of EBs found in the TESS data using a machine learning approach which are in or near the NCVZ of TESS and had a sufficient amount of photometric data.We have obtained the light curves of every source from all available TESS sectors up to Sector 60 in an automated way.After that, using our newly developed Python GUI, we have visually vetted all of the light curves of our targets along with determining their eclipsing periods and calculating their ETV curves.In the end of this process, we have selected 351 objects for further detailed analysis that seemed to have non-linear ETVs which could be caused by previously unknown third bodies in these systems.Finally, we tried to fit these ETVs with analytical formulae of LTTE or LTTE + dynamical perturbations.
In total we were able to derive models for 94 objects with pure LTTE and 41 objects with LTTE + dynamical solutions which we divided into three and two subgroups based on the quality of the models, respectively.The first subgroup contains the most reliable, highly certain solutions for 19 and 27 objects with the two types of ETV solutions, respectively.The second subgroup has 17 and 14 objects with likely valid solutions, but with probably higher uncertainties in their orbital parameters.For the pure LTTE systems only, the third subgroup consists of 58 objects having very uncertain solutions that are still likely triple systems, but their parameters should be used with considerable caution.That makes a total of 135 hierarchical triple candidates in our NCVZ sample of which 10 objects were already known from previous studies, while the rest of the 125 systems are new discoveries.Among these systems there are a number that are individually interesting.First, there are five very tight systems with a period ratio of less than 20 which are the closest to the dynamical stability limit.It is worth noting that, although our analytic models describe these systems well on a shorter timescale of the available TESS data set, they are inaccurate on decade-long timescales and numerical integration of the orbits should be applied to model these systems correctly.There are also four triply eclipsing triples in our sample, out of which one (TIC 29785001) was already known and for that one our ETV solution is consistent with the previous one of Rappaport et al. (2023).The other three systems are new discoveries and we plan a separate, more detailed comprehensive analysis in the near future for them.
We have outlined the case of GZ Dra in which we have recognized and successfully modelled the previously unknown effect of dynamical perturbations.This analysis showed that GZ Dra is actually one of the most inclined triple systems known with i m = 58 • ± 7 • .Although this inclination is well within the range of the original von Zeipel-Lidov-Kozai effect, we cannot expect large amplitude eccentricity cycles because of the presence of strong tidal interactions between the components of the highly circularized inner pair.Nevertheless, one would expect large amplitude orbital plane precession in the system which would cause EDVs in its light curve, but the timescale of this variation is on the order of a thousand years and that is the reason we cannot see its effect in the few years long TESS data set.
Because we have also determined all the same parameters for every triple candidates in our sample, we are able to make a comparison with the distributions of Borkovits et al. (2016) found for the full Kepler sample.For this purpose, we used only those 77 systems that have a certain or likely valid ETV solutions (i.e. the first two subgroups).In this sample of ours, there is a slightly elevated number of systems with an outer orbital period of 100-800 days that is most likely caused by the observing strategy of TESS.In the case of outer eccentricities we have found that their cumulative distribution is similar to those found by earlier surveys (Kepler, OGLE, Gaia), but very different from theoretical (flat and thermal) distributions.We have not found any correlation between the outer periods and eccentricities for systems in our sample.Regarding the masses of the outer components in our systems, we have found that the majority of them are below 1.8 M ⊙ which is similar to the Kepler sample.Among those 41 systems in which dynamical perturbations occur, and exact masses can be derived for the components, we have found that only four have a more massive tertiary component than their inner binaries.In the Kepler sample the majority of the systems have an inner binary mass lower than 2 M ⊙ , while in our sample it is the opposite complementing the regime of triples with higher inner binary masses.Similarly to what our other studies of TESS binaries have shown (Rappaport et al. 2023), the vast majority of these systems also have negligible mutual inclination (i.e. coplanar configuration) and only a few systems have distinctly non-zero mutual inclination.The latter also implies that we cannot see any peak at around 40 • as found by Borkovits et al. (2016) for the Kepler sample.
Finally, we have crossmatched our list of targets with the Gaia NSS catalog and found 15 systems that have either astrometric or spectroscopic orbital solutions with substantially different orbital period than their eclipsing period.For all of these objects, our ETV solutions can confirm that they are triple systems, nevertheless 5 of them were already found and published by Czavalinga et al. (2023).We have also checked the same correlations as Czavalinga et al. (2023) for the orbital parameters of these systems derived from the two different methods and we can draw the same conclusions as they did.The Gaia NSS catalog has great potential in discovering new hierarchical triples, the most reliable parameter in these solutions is the outer orbital Article number, page 16 of 28 Table 7. Outer orbital parameters from the Gaia NSS and our TESS ETV solutions for those 15 triple candidates that have both.For each object, the upper values belong to the former, while the lower values correspond to the latter solution, respectively.period, while the other parameters can have larger uncertainties than found via the ETV fits.
In this study, we have demonstrated that TESS data -in addition to many other applications -is an excellent tool for finding and analyzing new hierarchical triple stars.Since a lot of systems in our sample have uncertain solutions and more observations would be needed in order to improve their quality, we hope that TESS will continue monitoring these objects and their reanalysis will become possible in the not too distant future.

Fig. 3 .
Fig. 3.Primary ETV curve of GZ Draconis, formed from the 2-min TESS measurements (red dots), together with the combined LTTE+dyn solution (black curve).We plot separately both the LTTE (light blue) and quadrupole-order dynamical (light red) contributions.The unusual bumps near to the upper extrema of the observed curves are clearly visible and, indicate a significant, non-coplanar dynamical contribution.See text for further details.

Fig. 6 .
Fig.6.Cumulative distribution of outer orbital eccentricities for our TESS NCVZ (red) systems along with samples from earlier studies as Kepler (blue), Gaia (magenta), OGLE (orange) and theories as flat (green) or thermal (black) distributions.
Fig. 7. Correlation of outer orbital periods and eccentricities for the Kepler (blue) and TESS NCVZ (red) samples.Note, the correlation value pertains to the NCVZ set only.

Fig. 8 .
Fig. 8. Distribution of tertiary masses for the Kepler and TESS NCVZ samples directly measured from combined dynamical+LTTE ETV solutions (magenta and green, respectively) and lower limits estimated from pure LTTE solutions (blue and red, respectively).

Fig. 9 .
Fig. 9. Correlation of tertiary and inner binary masses for the Kepler (blue) and TESS NCVZ (red) samples with combined dynamical + LTTE ETV solutions.

Fig. 10 .
Fig. 10.Distribution of mutual orbital inclinations for the Kepler (blue) and TESS NCVZ (red) samples with combined dynamical + LTTE ETV solutions.

Fig. 12 .
Fig. 12. Correlations between the orbital periods (upper left panel), projected semi-major axes (upper right panel), eccentricities (lower left panel) and argument of periastrons (lower right panel) coming from the Gaia NSS solutions and fitting LTTE models to the TESS ETVs of the NCVZ sample.Different colors show different types of Gaia NSS orbital solutions.The yellow dashed lines represent a ratio of unity, while for the bottom right panel the green dashed line shows the ±180 • difference between the argument of periastrons coming from the two kinds of solutions.In order to avoid breaking the dashed lines in the bottom right panel, we employ a ±360 • shift for some of the ω out,NSS values which is equivalent because of the periodicity of the argument of periastron.The faded dots are additional previously known systems fromCzavalinga et al. (2023)

Table 1 .
• − 89 • , as we have done.Comparing these triple systems with the pure LTTE triples, one can see, that in general they contain less close inner binaries.(The mean and the median of the inner periods are P 1 = 9 d .341and P med 1 = 7 d .539,respectively.)This property may come about Article number, page 7 of 28 A&A proofs: manuscript no.TESS_NCVZ_clean Triple system candidates with certain or very likely LTTE solutions.

Table 2 .
Triple system candidates with less certain but likely LTTE solutions.

Table 3 .
Triple system candidates with very uncertain, questionable LTTE solutions.

Table 6 .
Apsidal motion and/or orientation parameters from AME and dynamical fits.
: Allowing non-coplanar solution resulted in i m = 43 • ± 10 • , but we do not trust in it; b : Adjusted mutual inclination resulted in i m = 23 • ± 7 • which would lead to ∆i 1 ∼ 27 • during TESS observations and hence, certainly overestimated; c : Third-body eclipses; 2+1+1 configuration; photodynamical solution from Rappaport et al. (2023); d : Adjusted mutual inclination resulted in i m = 22 • ± 7 • which would lead to ∆i 1 ∼ 4 • during TESS observations and hence, certainly overestimated; e : Grazing third-body eclipses + EDV f : Eclipse depth variation; g : Adjusted mutual inclination resulted in i m = 11 • ± 2 • which would lead to ∆i 1 ∼ 6 • during TESS observations and hence, certainly overestimated; h : Disappearing eclipses; i : Adjusted mutual inclination resulted in i m = 13 • ± 7 • which would lead to ∆i 1 ∼ 9 • during TESS observations and hence, certainly overestimated; j : No secondary ETV curve; k : Third-body eclipses a