Observational Detection of Higher Order Secular Perturbations in Tight Hierarchical Triple Stars

In this work, we search for observational evidence of higher-order secular perturbations in three eclipsing binaries. These are slightly eccentric binaries, and they form the inner pairs of tight, compact, hierarchical triple star systems. We analyze simultaneously the high precision satellite ($Kepler$ and $TESS$) light curves, eclipse timing variations, combined spectral energy distributions (through catalog passband magnitudes) and, where available, radial velocities of KICs 9714358, 5771589 and TIC 219885468. Besides the determination of robust astrophysical and dynamical properties of the three systems, we find evidence that the observed unusual eclipse timing variations of KIC 9714358 are a direct consequence of the octupole-order secular eccentricity perturbations forced by an unusual, resonant behaviour between the lines of the apsides of the inner and outer orbital ellipses. We also show that, despite its evident cyclic eclipse depth variations, KIC~5771589 is an almost perfectly coplanar system (to within $0.3^\circ$), and we explain the rapid eclipse depth variations with the grazing nature of the eclipses. Finally, we find that the inner pair of TIC~219885468 consists of two twin stars and, hence, in this triple there are no octupole order three-body perturbations. Moreover, we show that this triple is also coplanar on the same level as the former one, but due to its deep eclipses, it does not exhibit eclipse depth variations. We intend to follow this work up with further analyses and a quantitative comparison of the theoretical and the observed perturbations.


Introduction
The new window which was opened up by Kepler spacecraft [1] to the exotic world of the tight, compact hierarchical triple stellar systems (TCHTs) has now become more widely opened and transparent with the ongoing operation of the TESS space telescope [2] mission.While Kepler had an inevitable primacy in the discovery and principal characterizations of this rare new class of stellar systems, the multiply repeated revisitations of TESS, over five years, have allowed us to extend our knowledge about the properties of such systems.In a wider perspective, our studies of the stellar three-body problem have enhanced our understanding not only quantitatively, but even qualitatively, for the reasons which we discuss briefly in the forthcoming text.
Similar to our previous works [e.g. 3] we consider hierarchical triple or multiple (stellar) systems to be 'tight', if the orbital period ratio of the outer orbit to the inner orbit remains below ∼ 100 (i.e., P out /P in ≲ 100), while we call a triple or multiple star system 'compact', when the outer period does not exceed, let's say, 1000 days. 1 As is widespreadly accepted [see, e. g. [4][5][6], the only possible long-term stable configuration amongst three masses, being of the same order of mass, is the hierarchical configuration, i.e., when one of the three bodies is continuously located much farther from the other two than these latter two from each other.In such a system the motion, in general, can be well approximated with two binaries, i.e., two Keplerian two-body motions.The inner binary is formed by the two closer objects, while the outer binary consists of the third, most distant component, and the center of mass of the inner binary.Hence, tightness chiefly characterizes the strength of the third-body perturbation of such systems.This is so since the magnitude of the gravitational perturbations relative to the pure Keplerian two-body motions primarily depend on the ratio of the semi-major axes or, more strictly, of the instantaneous separations of the outer to the inner binaries.The former can readily be converted to such a direct observable as the period ratio.
On the other hand, compactness is directly connected to the physical dimensions of the system as a whole.This characteristic size of the given triple star, naturally, has very crucial astrophysical implications, not only in regard to the early and the late evolutionary phases of such systems [see, e.g.7-9, and references therein], but how tidal and other non third-body perturbations relate to the gravitational perturbations of the third mass.From our point of view, compactness primarily has some practical significance.This significance has two, somewhat counteracting aspects.First, the largest amplitude, and most interesting, orbital perturbations in tight triple star systems are the so-called long-period or apse-node perturbations -see below -which are of primary interest in this paper.The characteristic periods of these perturbations scale with the ratio of P2 out /P in , i.e., these are essentially driven by the product of the measures of the tightness and the compactness.Hence, the more compact a tight triple is, the shorter the observational window or the length of the necessary dataset for the direct detection and investigation of such effects.On the other hand, however, such triples can be serendipitously discovered and characterized most easily through the medium-period third-body perturbations (see again, below) that drive eclipse timing variations (ETVs) of such eclipsing binaries (EBs) which form the inner binary of a tight hierarchical triple system.The amplitude of such medium-timescale dynamically driven ETVs, however, scales as P 2  in /P out and, in such a way, for the most frequently occurring and shortest period EBs (some hours to few days), the dynamically driven ETVs may remain under the detection limit even in the tightest triples.This latter fact explains why TCHTs were almost completely unknown before the fouryear-long observations of the prime Kepler mission.The vast majority of the EBs known before Kepler had periods less than 3-4 days, in which case there was no chance of detecting the gravitational third-body perturbations via eclipse timing variations with the accuracies of the highly inhomogeneous ground-based observations. 2The quasi-continuous, almost four-year-long observations of the Kepler space telescope has led to the discovery of more than a thousand new EBs with eclipsing periods longer than 5 days.Amongst them, Borkovits et al. [11,12] have identified 46 tight hierarchical triple star candidates through the manifestations of the medium-term (P out -period) third-body perturbations in their ETVs.(Note, they found 16 further tight triple star candidates amongst the almost 2000 additional, shorter period Kepler EBs, but this was mainly a consequence of the unprecedentedly precise observations of the space telescope and/or other effects, such as, e.g., third-body eclipses.)From this sample of the newly discovered 62 tight triple systems, 39 can be ranked as TCHT, i.e., a triple star system which is not only tight, but compact, as well.While these former studies depend primarily on the P out -timescale, the so-called medium period perturbations, the four-year-long ETVs of most of the eccentric TCHTs also exhibited signals of rapid, dynamically driven apsidal motion.Hence, a simplified, quadrupole approximation of the dynamically forced apsidal motion (and also of the nodal regression) was included into the analytic ETV model.
The TESS spacecraft, Kepler's successor, first reobserved the prime Kepler-field almost exactly a decade after Kepler began its survey.Since then, the vast majority of the original Kepler targets and, hence, the Kepler discovered TCHTs have been revisited typically 4-7 times, each for a four-week-long observing session.Thus, currently the lengths of the available datasets are nearly one and a half decades.Naturally, 14 years is tiny in the context of most interesting astronomical timescales, but in the case of TCHTs, this interval is long enough for the robust detection and quantitative analysis of the 'long-period' (or, secular) third-body perturbations.In this paper we have selected a few such TCHTs from the original Kepler sample, which are exceptional even amongst the TCHTs.In particular, their special properties make possible the clear detection and investigation of even the higher order secular perturbations with the exclusive use of the high-precision satellite data.In what follows, first we give a brief review in Sect. 2 of the analytic theory of the doubly averaged or, secular stellar three-body problem, focusing on the observable quadrupole and octupole perturbation terms in the low mutual inclination domain.Then, in Sect. 3 we introduce the three TCHTs that were selected for the current analysis.The details of the observational materials, the data preparation, and the complex photodynamical analysis are described in Sects.4 and 5, while the results are discussed in Sect.6.

Dynamics of TCHTs
The hierarchical three-body problem is sometimes called the 'stellar three-body problem', however, the fields where the hierarchical treatment can be applied are much wider than the domain of triple and multiple stars.For example, the hierarchical three-body treatment can be used for investigations of the perturbations of the Earth-Moon system due to the Sun [13][14][15], Jupiter's effects on the motion of the main belt asteroids and/or, comets [16,17], and the perturbations on a binary comprised of the Earth and an artificial satellite due to the Moon [18].
According to our knowledge, Slavenas [19,20] was the first one who investigated the hierarchical three-body problem in the context of a stellar triple (λ Tauri, the very first and, for a long-time, the only TCHT).Due to the fact, however, that there were almost no known tight, or even compact, triple stellar systems, the hierarchical or, stellar three-body problem, apart from a very few exceptions [see, e. g. [21][22][23][24][25][26], has remained beyond of the scope of the scientific community.The situation did change completely, and quickly, when the Kozai mechanism (nowadays called the 'Von Zeipel-Lidov-Kozai' -ZLK -mechanism) or cycles were essentially 'rediscovered' in regard to the formation and evolution of stellar triples during the late 1990's.At that time the first observational evidence was found in triple stellar and planetary systems for the effectiveness of secular third-body perturbations [see, e. g. [27][28][29] and, practically, at the same time, in an epochal paper Kiseleva et al. [30] proposed a combined mechanism of third-body perturbations with tidal interactions to explain the formation of the close binary systems.In our view, this last work transformed the ZLK mechanism into its current, widespread and acknowledged position.
The mechanism, which is currently known as ZLK theory, is practically nothing more than the secular, doubly averaged theory of the hierarchical or, stellar three-body problem.This theory describes the long-term dynamical evolution of such a three body system, where the evolution is driven, in its pure form, only by the mutual gravitational interactions of the three bodies.
If we slightly modify the original classification of Brown [14], in a hierarchical triple system the perturbations are effective over three different timescales, and they can be classified according to these characteristic timescales as follows: • short-period perturbations, for which the typical period is in the order of P in , and the amplitude is related to (P in /P out ) 2 , • medium-period perturbations, with a characteristic period of P out , and amplitude of P in /P out and, • long-period perturbations3 having period about P 2 out /P in , and an amplitude that may reach unity.
It was Harrington [21,22] who, for the first time, described the equations of the stellar three-body problem in Hamiltonian formalism for any arbitrary values of outer eccentricities and mutual inclinations and, gave third-order (sometimes called: 'octupole') solutions for the long-period or, apse-node time-scale variations of the orbital elements, with the application of the double-averaging method of von Zeipel [31,32,33,34].
While the current discoveries of TCHTs (including the currently analyzed systems, as well) are mainly based on the medium-period perturbations of the ETVs of their inner EBs, the scope of this paper is connected to the long-period perturbations, and thus, in what follows, we concentrate only on the latter.The original ZLK theorem is restricted to the lowest order perturbative terms of the doubly averaged Hamiltonian (which is quadratic in the small parameter α = a in /a out , hence, it is called the 'quadrupole approximation').At that level, the problem has only one degree of freedom, and thus, it is integrable, and the solution can be given with the use of elliptic functions of Weierstrass P. Note, although the original analytic solution of Kozai [17] is strictly valid only if the outer orbit coincides with the invariable plane of the triple (i.e., all the angular momentum of the system is stored in the outer orbit), as was shown, e.g. by Harrington [21] and Söderhjelm [24], the same solution may remain valid, with small modifications, even in the case where the orbital angular momentum of the inner orbit is non-negligible, but small enough.On the other hand, however, as was shown by Ford et al. [35] via numerical integrations and, later analytically by Naoz et al. [36] (amongst others), in the case of an eccentric outer orbit, the octupole order terms may substantially alter the long-term behaviour of a hierarchical triple system.
As was excellently reviewed by Naoz [37], the vast majority of the studies of the doubly averaged stellar three-body problem (being either analytical, numerical or, both) are directed toward the investigation of the very long-timescale evolution of such systems.Specifically their interests are mainly focused on configurations where large amplitude eccentricity and (mutual) inclination cycles (including prograde-retrograde flip-flops) may occur.Much less effort has been directed toward the short-term (human timescale), directly observable effects of the stellar three-body problem.In this latter context, besides the above mentioned few exceptions, from the point of view of the present work, the paper of Krymolowski & Mazeh [38] is extremely relevant.These authors compared the outcomes of the quadrupole and octupole order perturbation theories of the doubly averaged problem to each other (and also to the results of direct numerical integrations) for different mutual inclinations (from the coplanar case up to i mut = 60 • ).They found that the octupole perturbations may be very significant even in the coplanar (or, nearly coplanar) case.This is so because, in an exactly coplanar configuration, the long-period quadrupole perturbations, for example, in the inner eccentricity completely disappear.Despite this fact, numerical integrations show that there exist long-period perturbations in the inner eccentricity in the coplanar configurations, and this effect is chiefly determined by the octupole order perturbations.
As those TCHTs that we investigate in this paper are actually nearly coplanar systems, in what follows, first we discuss the octupole effects quantitatively.Following the formulae of Harrington [21], but using partly different notations, the doubly averaged Hamiltonian up to the octupole order takes the form: where e in,out and ω dyn in,out denote the eccentricities and dynamical arguments of periastron of the inner and outer orbits, respectively.These latter quantities (sometimes denoted as the Delaunay variables g in,out ) -should not be confused with the observable arguments of periastrons (the latter of which will hereafter be denoted as ω in,out ) -give the angular distances of the pericenter points of the orbits from the ascending nodes of the intersections of the orbital planes.Furthermore, I = cos i mut stands for the cosine of the mutual (or, relative) inclination of the two orbital planes.Finally, where we have introduced the parameters A G and A oct G which are more closely related to observable quantities, and can be expressed as Additional, newly introduced quantities in the equations above are the gravitational constant G, the individual masses of the inner, close binary stars, as m Aa , m Ab , the total mass of the inner pair, m A = m Aa + m Ab , the individual mass of the third component, m B , and the total mass of the triple system m AB = m A + m B .Moreover, the mass ratios are denoted by q in,out for the inner and outer binaries, respectively.In regard to the amplitude-like quantity of A oct G , as one can see, it disappears for q in = 1, i e., when the inner pair consists of two equal mass stars. 4It is common also to introduce the quantity of which gives the strength of the octupole order terms relative to the quadrupole ones.
Finally, note, that we use the usual Delaunay action variables as while the conjugate angular variables are the mean anomalies (l in,out ), dynamical arguments of pericenters (ω dyn in -or, traditionally, g in,out ), and dynamical longitudes of the nodes (Ω dyn in,out -or, h in,out ) of the inner and outer orbits, respectively.One should also keep in mind that, as was noted by Naoz et al. [36], although Eq. (1) formally does not contain the nodes (Ω dyn in,out ), which would suggest the constancy of the conjugate variables H in,out , this arises only from an incorrect application of the elimination of the nodes.For a correct treatment, one should take it into account that, in the original Hamiltonian, both the inner and outer dynamical nodes are present, but only through the sines and cosines of their differences, which is a constant ∆Ω dyn = Ω dyn in − Ω out in = π and, therefore, sin ∆Ω dyn = 0, and cos ∆Ω dyn = 1.If one substitutes these values into the original Hamiltonian, it leads to the usual form of Eq. (1).Hence, strictly speaking, from the simplified (Eq. 1) form of the Hamiltonian, the constancy of H in,out does not follow.This fact, however, does not affect the calculation of the perturbation equations of the other elements through the usual, formal way, and hence, in what follows, we do not take it into account.
The perturbation equations that are interesting for us, take the following form: In what follows, we consider only the prograde, coplanar scenario, i.e., when i mut = 0 • .This restriction can be justified by the fact that, as will be shown in Sect.6, the mutual inclination in none of the three considered systems exceed i mut = 0.5 • .In the coplanar case, when I = 1, Eqs. (12)(13)(14)(15)(16)(17)(18) become much simpler.Before we give these simpler expressions, we note that, for coplanar orbits, the dynamical arguments of periastron, and the dynamical nodes lose their meanings.The quantities that continue to remain meaningful in this situation are the dynamical longitudes of periastron (ϖ dyn in,out = ω dyn in,out + Ω dyn in,out ), and also the angle between the directions of the two periastron points, i.e., ω where in the last expression we took into account the fact that ∆Ω dyn = π.Furthermore, we notice that in the coplanar case ω in,out = ϖ dyn in,out .Thus, as can be easily seen, for the prograde, coplanar case, the perturbation equations up to the octupole order have the following forms: and, finally, e in e out As can be easily seen, in this scenario there are no quadrupole level long-period perturbations in the inner eccentricity, and the observable dynamical apsidal motion rates of both orbits remain constant (again, at the quadrupole level).In contrast to this, the octupole order perturbations of these elements depend on trigonometric functions of the angle between the two periastron directions.One should keep in mind, that this angle, in general, varies at a lower rate than the individual apsidal lines themselves and, hence, one expect a longer period cyclic, octupole variation, than the usual (quadrupole) dynamically forced apsidal motion period.
A detailed analytic and/or numerical analysis of these perturbation equations are beyond the scope of the current paper.Our aim remains only to identify the consequences of these orbital perturbations directly from the high-precision observations of selected TCHTs, and then, to determine accurate dynamical (and astrophysical) parameters of these systems, and finally, to compare our findings at least qualitatively with the predictions of the theory, discussed briefly above.A more detailed study including the quantitative comparison of the theoretical and the observed perturbations will be published later.

Selected systems
We selected two TCHTs from the well-observed primary Kepler-sample and, moreover, a third system, which serves as a counterexample in the sense that, although this system is very similar to the previous two in several aspects, the octupole order perturbations are almost nulled out due to a near unit mass ratio of the inner binary.Unfortunately, we did not find such an illustrative TCHT amongst the Kepler systems and, hence, we took one from the TESS northern continuous viewing zone (NCVZ) sample.The main catalog data for the three systems are collected in Table 1, where, in addition to the J2000.0coordinates and the catalog magnitudes from near ultraviolet to infrared passbands (i.e., form GALEX NUV to Wise W4 magnitudes) we also tabulate the effective temperatures (T eff ), photogeometric distances, metallicities ([M/H]), interstellar reddenings (E(B − V)) and proper motions (µ α,δ ) from different catalogs.Finally, we also give the specific RUWE (renormailzed unit weight error) parameters, which was introduced in Gaia DR2 as an indicator of the quality of the astrometric solutions.This parameter is relevant in the context of the current study because it was found that a value greater than ≳ 1.4 could indicate the multiplicity of the source [see, e. g . 39].In what follows, we discuss briefly the available basic information on these three TCHTs.
KIC 9714359 = KOI-6073 (P in = 6.48 d; P out = 103.8d; P out /P in = 16.0,ϵ ≈ 0.016) is the compact, tight triple star system which chiefly inspired this study.This target was identified first as an EB in the first release of the Kepler EB catalog [48].The hierarchical triple star nature of this source was reported first by Rappaport et al. [49] analyzing the ETVs of the first 13 quarters of data of Kepler EBs, who also pointed out that the ETV is dominated almost exclusively by the dynamical, third-body perturbation effects over the classic light-travel time effect (LTTE).Then the system was included in the sample of 26 eccentric EBs which were comprehensively reanalyzed by Borkovits et al. [11] with the use of their most complex analytic ETV model.Despite the extreme tightness of the triple, they found a quite satisfactory analytical, coplanar ETV solution with a slightly eccentric inner orbit, and moderately eccentric outer orbit (e in = 0.015 ± 0.001 and e out = 0.30 ± 0.01, respectively), and with a dynamically forced apsidal motion period of P apse,in = 30.6yr.The Notes.(a) TESS Input Catalog (TIC v8.2) [40].(b) Gaia EDR3 [41].(c) PanSTARRS, [42].
(g) WISE point source catalog [46].(h) Photogeometric distances from Bailer-Jones et al. [47].Note also, that for the SED analysis in Sect. 5 the uncertainties of the passband magnitudes were set to σ mag = max(σ catalog , 0.030) to avoid the strong overdominance of the extremely accurate Gaia magnitudes over the other measurements.
model also indicated inner and outer mass ratios of q in = 0.45 ± 0.01 and q out = 0.21 ± 0.01, respectively.Further studies, however, resulted in slightly discrepant values for some of the parameters above.First, Kjurkchieva et al. [50] found a substantially (∼ 5 − 6 times) larger inner eccentricity of e in = 0.085 from the Kepler light curves, and, moreover, based on astrophysical reasons, they estimated a lower inner mass ratio (q in = 0.33).Then, Windemuth et al. [51] also reanalyzed the Kepler light curves, adding SED information and theoretical PARSEC [52] evolutionary tracks to them in an automated way.While their eccentricity of e in = 0.016 ± 0.002 was really in good agreement with that of the dynamically inferred value of Borkovits et al. [11], they obtained an even much more extreme inner mass ratio (q in = 0.234) and, correspondingly, they found that the system was in a pre-MS state with an age of log τ ≈ 7.44 ± 0.04.
While these discrepancies, noted above, might be interesting in and of themselves, what makes this system especially important for us is the unusual behavior of the ETV curves during the TESS observations.KIC 9714358 was reobserved by TESS in Sectors 14, 15, 40, 54 and 55 in the lower cadence full frame image (FFI) mode.As one can readily see in the upper left panel of Fig. 1 the new ETV points, determined from TESS-observed eclipses, do not follow the former trend, and the long-term behaviour of the red (primary) and blue (secondary) ETV curves can no longer be described with two antiphased cosine curves, as one would expect for an EB with low eccentricity and (almost constant) apsidal motion rate (the latter of which is a clear prediction of the quadrupole-order hierarchical three-body problem).Hence, our immediate suspicion was that what we are seeing is nothing more than a direct, observational manifestation of the octupole-order perturbation effects.We will return to this question in Sect.6.
Finally, note that this is the only system in our sample for which radial velocity (RV) data are available from the literature.Data release 16 of the APOGEE-2 project [53] tabulates 16 individual RV observations, which we used for our analysis (see Fig. 3).
KIC 5771589 = KOI-6625 (P in = 10.79 d; P out = 113.1 d; P out /P in = 10.5, ϵ ≈ 0.004).This system was identified as an EB, and then as a triple star candidate in the same sequence as for the previous target.It was also included in the detailed analytical ETV studies of Borkovits et al. [11,12].This is the second tightest known triple in the whole studied Kepler sample, and the only one for which the ETV curves, determined purely from the ∼ 4 yr-long original Kepler observations, cover more than half of an apsidal cycle. 5The analytical ETV studies suggested a large outer mass ratio (q out = 0.98 ± 0.06), which was strongly supported by the very shallow eclipses (∼ 0.2 − 0.4%).
The most remarkable feature of the Kepler light curve is, however, the varying depths of the eclipses (see Fig. 4), implying orbital plane precession forced by the not-exactly coplanar tertiary.According to our knowledge, there are only two EBs in the whole primary Kepler sample which exhibit a clear reversal of the trend of the eclipse depth variations (EDVs) during the four-year-long observations.The other such target is the triply eclipsing triple star system KIC 6964043 [3], which displayed an opposite trend in its EDV compared to KIC 5771589.In the former case the eclipse depths increased first and then decreased, while in the currently considered KIC 5771589 the eclipses became shallower and shallower during the first half of the Kepler observations, and then they started to increase again.These two situations are not symmetric with respect to each other.This is because an EB can reach the maximum of its eclipse depths in two different ways during a nodal precession cycle.The EB may either pass through the exact edge-on view (i.e.i = 90 • -which, as was pointed out by Borkovits et al. 3, really did happen in KIC 6964043) without reaching an extremum in its observable inclination or, it may reach its nearest position to an edge-on view (i.e.| cos i| has a minimum, but not zero value or, of course, the inclination itself has an extremum).In contrast to this, when eclipse depths reach a minimum value, one can be certain that the visible inclination has reached its extremum farthest from i = 90 •6 .In such a way, the current light curve behaviour offers more strict constraints on the orbital configuration.
KIC 5771589 was reobserved with the TESS space telescope in Sectors 14, 26, 40, 41, 53 and 54.In contrast to KIC 9714358, for this target 2-min cadence observations are also available for all sectors.Unfortunately, despite the relative brightness of the system, due to the very shallow eclipse depths, the TESS data have much lower quality.Despite the low S/N ratios, however, very shallow eclipses can still be identified at the expected times, but in the Year 2 data they are almost lost in the noise.By contrast, the eclipses can be seen clearly at the end of the Year 4 observations, indicating that the eclipse depths must have increased again in between the Year 2 and Year 4 TESS observations (see the lower right panel of Fig. 4).
TIC 219885468 has P in = 7.54 d; P out = 110.8d; P out /P in = 14.7, and ϵ ≈ 0.0004.This system was not observed with the Kepler spacecraft, but instead, it is located near to the northern ecliptic pole, in the NCVZ of the TESS space telescope.But, we did not choose it primarily for its relatively long data train, but rather for didactic reasons.In particular, we wanted to illustrate that, despite the clear similarity to the above-mentioned KIC 97154358 (in both their inner and outer periods), due to the twin-nature of the inner binary star, and hence, in the absence of octupole order perturbations, the short-term orbital evolution of the two triples differ significantly.
TIC 219885468, a previously unknown EB, was identified as a likely triple star candidate during our ongoing study of TESS EBs in the NCVZ which surveys the ETVs of these objects to identify triple star candidates through their eclipse timing variations (Mitnyan et al., to be submitted soon).The ETVs of this EB display typical timing variations that are dominated by third-body perturbations, including medium-term effects as well as dynamically forced rapid apsidal motion (see bottom row of Fig. 1).The primary and secondary eclipses look very similar in depth, which indicates that the inner mass ratio should be near unity-and this does null out the octupole perturbation terms.Moreover, no eclipse depth variations are observed during the ∼ 1300 days of TESS observations, which suggests a strong coplanarity (see, Fig. 5).

Observational data, and its preparation for the analysis
Two of our three targets were observed nearly continuously in long cadence (LC; time resolution of 29.4 min) mode during the prime Kepler mission.Furthermore, for KIC 5771589 one month of short cadence (SC; time resolution of 58.9 sec) data in Quarter 4.1 is also available.For our analysis we used those preprocessed LC datasets which can be publicly downloaded from the Villanova website 7 of the third revision of the Kepler Eclipsing Binary Catalog [56,57].Quarters 4, 12 and 13 data of KIC 9714358, however, are not available at this site, hence, these were downloaded directly from the MAST database 8 .
Regarding the TESS observations, we downloaded the FFIs of all three targets for all sectors, when they were observed.(As is well known, Year 2 FFIs were obtained with 30-min cadence time, Year 4 ones with 10 minutes, while Year 5 FFIs -available only for TIC 219885468 among our targets -have 200-sec cadence time.)We then processed these FFIs using a convolution-based differential photometric pipeline implemented in the FITSH package [58].Finally, we detrended these raw light curves, which were affected by differing amounts of scattered light and other effects, with smoothing polynomials of various degrees to each light curve section.
After the detrending processes, as a first step, we determined mid-eclipse times for each individual eclipse, in the very same manner as was described in Borkovits et al.  [12]. 9Then, preparing the dataset for the light curve analysis section of our complex photodynamical treatment, first, we binned the Year 4-5 data into 1800 sec bins for a homogeneous sampling (not only with respect to Year 2 TESS data, but also to the Kepler LC data).Then, in order to reduce computational costs, and also to give higher statistical weight to the eclipses (which carry most of the relevant astrophysical and dynamical information amongst the light curve points), we dropped out the larger parts of the outof-eclipse light curves and, kept only small light-curve sections around the primary and secondary eclipses.
In the case of KIC 9714358 we downloaded the 16 heliocentric RV data points, published in APOGEE DR 16 [53], from the Vizier site 10 and also used them.
Finally, note that, similar to our previous works, we included a composite spectral energy distribution (SED) analysis into our joint photodynamical runs (see below) for which we mainly used those catalog magnitudes which are tabulated in Table 1.

Photodynamical modeling
We have carried out complex photodynamical analyses of our three targets with the software package LIGHTCURVEFACTORY [see, e.g.59,60, and references therein].This code contains (i) a built-in numerical integrator to calculate the gravitationally (three-body effects), tidally, and (optionally) relativistically perturbed Jacobian coordinates and velocities of the three stars in the system; (ii) emulators for multiband light curves, ETV curves, and RV curves, (iii) built-in, tabulated PARSEC grids [52] for interpolating fundamental stellar parameters (e.g.radii and effective temperatures) as well as a composite SED model for the three stars, and (iv) an MCMC-based search routine for fitting the parameters.The latter applies our own implementation of the generic Metropolis-Hastings algorithm [see e.g.61].The use of this software package and the consecutive steps of the entire analysis process have been previously explained in detail in several former papers [see, e. g. 59,60,62,63].In this regard we note that, despite the fact that LIGHTCURVEFACTORY has a built-in subroutine for correcting the emulated light curves for finite cadence times [see 59, for the description of the process], due to the relatively long cadence times of 1800 sec, we did not applied such corrections.This decision can be justified a posteriori by the fact, that applying cadence time corrections in the case of KIC 9714358, resulted in parameters that remained within the ±1σ domains of the uncorrected results (tabulated in Table 2) -even   in the case of the most sensitive parameters (the inclinations).However, since the analysis runs with the cadence time corrections result in computation times about four-five times longer (on average), we declined to use this feature routinely.
In the absence of an RV curve for each of the three stars in any of our triple systems (or even for two of the stars), we combined a composite SED analysis and the use of precalculated PARSEC grids as proxies to determine the masses of the constituent stars.Naturally, such a photodynamical model solution is no longer astrophysically modelindependent.Nonetheless, in a previous work [64] we have shown that for compact, dynamically interacting triple stars (as is the case for the currently investigated systems), such a solution may result in stellar masses and radii within 5-10% of the results of an astrophysical model-independent solution, the latter of which is based on dynamical masses from RV data.The geometric and dynamical parameters (i.e. the orbital elements and also the mass ratios), as well as dimensionless quantities such as the fractional radii and temperature ratios of the stars remain practically unchanged between the two kinds of solutions.
In the case of these model-dependent runs, the adjusted parameters were as follows: (i) Stars: Three stellar mass related parameters: the mass of the most massive component (being either m Aa , the primary of the inner pair or, m B , the tertiary star), and the inner and outer mass ratios (q in,out ).Moreover, the metallicity of the system ([M/H]), the (logarithmic) age of the three coeval stars (log τ), and the interstellar reddening E(B − V) for the given triple were also varied.Additionally, the 'extra light' contamination, ℓ 4 parameter was also freely adjusted for two of the three systems.(ii) Orbits: Three of six orbital-element related parameters of the inner, and six parameters of the outer orbits, i.e., the eccentricity vector components of both orbits (e sin ω) in,out , (e cos ω) in,out , the inclinations relative to the plane of the sky (i in,out ), and moreover, three other parameters for the outer orbit, including the period (P out ), the longitude of the node relative to the inner binary's node (Ω out ), and the periastron passage times (τ out ) were adjusted.
A couple of other parameters were constrained instead of being adjusted or held constant during our analyses, as follows: (i) Stars: The radii and temperatures of the three stars were calculated with the use of three-linear interpolations from the precomputed 3D (metallicity; logarithmic age; stellar mass) PARSEC grids.Additionally, the distance of the system (which is necessary for the SED fitting) was calculated a posteriori at the end of each trial step, by minimizing the value of χ 2 SED .(ii) Atmospheric parameters of the stars: we handled them in a similar manner as in our previous photodynamical studies.We utilized a logarithmic limb-darkening law [65] for which the passband-dependent linear and non-linear coefficients were interpolated at each trial step via the tables from the original version of the Phoebe software [66].
We set the gravity darkening exponents for all late type stars to β = 0.32 in accordance with the classic model of Lucy [67] valid for convective stars and held them constant.Note, however, that the choice of this parameter has only minor consequences, since the stars in the present study are close to spheroids.(iii) Orbits: The orbital period of the inner binary (P in ) and its orbital phase (through the time of an arbitrary primary eclipse or, more strictly, the time of the inferior conjunction of the secondary star -T inf in ) were constrained internally through the ETV curves.Finally, in the case of KIC 9714358, the systemic radial velocity (V γ ) was constrained in a similar way, as was done with the distances, described above, i.e. this parameter was calculated a posteriori with a minimization of the value of χ 2 RV .The median values of the orbital and physical parameters, as well as some derived quantities, of the three triple systems, computed from the MCMC posteriors and their 1σ uncertainties are tabulated in Tables 2 -4.Furthermore, the observed vs. model lightcurves are plotted in Figs. 2, 4, 5, while the observed vs. model ETV curves are shown in Fig. 1.In the case of KIC 9714358 we also plot the APOGEE-2 RV data vs. the best photodynamical model, both in the time and the phase domains in Fig. 3.
While the majority of the tabulated parameters are self-explanatory, here we add some notes on the apsidal and nodal motion related parameters.We give the theoretically estimated apsidal motion periods both in the observational (P apse ) and the dynamical (P dyn apse ) frame of references.For their calculation we do not use exclusively the theoretical (quadrupole) model of the point-mass third-body perturbations, but the tidal and general relativistic effects are also taken into account.In this regard, we also tabulate the individual contributions of the third-body, tidal and relativistic effects to the apsidal advance rate (∆ω 3b,tide,GR ).Moreover, besides the apsidal motion parameters, we also give the theoretical, quadrupole, nodal regression period (P dyn node ), i.e., the time needed for a 360 • regression of the parameters Ω dyn in,out .Further details of these derivations are discussed in Section 6.2 of Kostov et al. [68].

Discussion
In what follows we discuss our results for the three investigated systems individually.Though we are primarily interested in the inferred dynamical behaviour of our triples, first we consider briefly the astrophysical implications of our results and then discuss the dynamics of the systems in detail.

KIC 9714358
In accord with the previous results of Windemuth et al. [51], we found this triple to be very young.Our inferred age of τ ≈ 46 Myr implies that the two less massive components of the triple must still be in the pre-Main Sequence contraction stage.This conclusion is imposed by the fact that the secondary is too large, while at the same time is too cool relative to the primary.To be more specific, the surface brightness (and, hence,   indirectly, the temperature ratio) and the ratio of the stellar radii are strongly constrained through the shape, durations and depth ratios of the primary and secondary eclipses, and our LIGHTCURVEFACTORY code was unable to find any reliable coeval PARSEC isochrone MS solutions for the two stars.Moreover, in the present situation the inner mass ratio (q in = 0.406 ± 0.006) is also very robustly constrained dynamically both from the ETV curve (where, however, the inner mass ratio occurs only at the octupole levels of both the medium-and long-term perturbations) and also from the RV curve of the primary component of the EB.Strictly speaking, as is well-known, in a single-lined spectroscopic binary (SB1) the amplitude of the RV curve constrains directly only the spectroscopic mass function, which can easily be written in the following form: This coupled with the known (inner) inclination of i in = 87.61• ± 0.06 • , and the primary's mass (at least to within a few percent uncertainty, as in the current situation), of m Aa = 0.83 ± 0.03 M ⊙ , uniquely gives the (inner) mass ratio q in = 0.406 ± 0.006.Thus, this very robust mass ratio clearly excludes any post-MS solutions and, hence, insofar as we assume that the triple (or, at least the two inner stars) have formed together at the same time, and without any substantial interactions (in particular, mass exchange) during their prior evolution, we can conclude that this system is at the very beginning of its life.
Comparing our results with former survey results, our SED solution clearly supports a hotter primary component with T Aa = 5 300 ± 100 K, than the one listed in either the TIC v8.2 catalog (4 764 ± 109 K, see in Table 1) or inferred from the APOGEE-2 spectra (4 807 ± 95 K, Jönsson et al. 53).In this context some further caution is needed because our inferred (photometric) distance was found to be d = 666 ± 11 pc, which substantially exceeds the Gaia EDR3 value of d EDR3 = 577 ± 7 pc, [47].One should keep in mind, however, that Gaia EDR3 parallaxes and, hence, the calculated distances have not yet been corrected for multiplicity.Thus, naturally, due to this discrepancy, the astrophysical implications of our results should be considered only with some caution.
On the other hand, these discrepancies do not influence the validity of our dynamical results.Besides the robust inner and outer mass ratios (the latter being q out = 0.265 ± 0.003), and the eccentricities of the two orbits (e in = 0.0286 ± 0.0003 and e out = 0.252 ± 0.002), here we mention the complete absence of any eclipse depth variations not only during the four years of the Kepler data, but also during the extended 13-year-long interval of the Kepler and TESS observations.This latter point very strongly supports our findings about the nearly exact coplanarity of the inner and outer orbital planes, which was found to be i mut = 0.07 • ± 0.03 • .(It will be shown in the next subsection that even a few tenths of a degree departure from exact coplanarity may lead to a robust detection of EDVs for such a tight system.) The most interesting dynamical feature of the current system is, however, the unusual behaviour of the lines of the apsides of both orbits, which, in our interpretation, is a clear manifestation of the octupole order perturbation effects and, it has clearly observable consequences.
The theoretical observable apsidal motion period of the inner EB of KIC 9714358, according to the quadrupole-level perturbation theory, should be P theo apse,in = 26.4 ± 0.2 yr (see the middle of Table 2).As one can readily see in the upper left panel of Fig. 1, the ∼5 000day-long observational dataset (which is longer than half the duration of the predicted full apsidal cycle) contradicts the quadrupole theory.In the left panel of Fig. 6 we plot the variations of the observable and dynamical arguments of periastron according to our best-fitted solution, extending the numerical integration of LIGHTCURVEFACTORY to the end of the current century.As one can see, the most characteristic effect is that the lines of the apsides of the inner and outer orbits (red and blue curves) revolve with the same speed, resulting in similar inner and outer apsidal motion periods of P meas apse ≈ 78 yr (which is quite close to the theoretical quadrupole outer apsidal motion period, P theo apse,out = 76.8 ± 0.6 yr).However, they do so in such a manner that, while the outer orbital ellipse rotates with a constant rate, the major axis (i.e. the apsidal line) of the inner orbit oscillates or, librates, around the direction of the outer apsidal line with a period of P osc apse ∼ P meas apse /2.5 and, with a half-amplitude of ∼ 15 • .(In other words, the difference of ω in − ω out oscillates around 0 • ).Naturally, the same is true for the dynamical arguments of periastrons with the difference here being that ω dyn in − ω dyn out oscillates around 180 • , as is shown nicely in the right panel of Fig. 6.
In the right panel of Fig. 6, besides the difference of ω dyn in − ω dyn out , we also plot the variations of the inner eccentricity (e in ) at the same times.As one can see, the inner eccentricity oscillates with the very same period and, with a 0.25 phase offset relative to the apsidal oscillations.This is in nice accord with the octupole order long-timescale perturbation equation for the eccentricity (see above, in Eq. 19).And, one can conclude that the currently observable uneven variations in the ETV curves of KIC 9714358, i.e. the momentarily rapidly diverging primary and secondary ETV curves, are the direct manifestations of the current growth of the eccentricity of the binary orbit due to the octupole-order terms.
At this point, however, one should note an important caveat.As was mentioned above, our analysis has shown that this triple system is quite young, and its components (or, at least, the less massive ones) are most probably in pre-MS stages.Such components may potentially exhibit strong spot activity and may rotate non-synchronously (and perhaps even in a non-aligned way).Both effects can affect the ETV curves either indirectly, through the spot-induced light curve distortions which may mimic rapid quasi-periodic and anticorrelated variations in the primary and secondary ETV curves [69][70][71] or, directly influencing the apsidal motion rate of the inner binary [72,73].In the current situation, however, we are convinced that the peculiar behaviour of the ETVs cannot be explained with these latter effects.First of all, in all the known cases, the time-scales of the quasi-periodic spot-induced variations range from a few days to a few months [see, again, the plots in 70,71] and, moreover, their amplitudes do not exceed 5-10 minutes.Moreover, as these ETVs are virtual and caused only by light curve distortions which affect the determinations of the mid-eclipse times, these distortions would appear in the residuals of our model light curves and, hence, may become directly detectable. 11On the other hand, regarding the possibility of the non-synchronous rotation, this may actually affect the apsidal motion rate, and, hence, the apsidal motion period, but would not change the eccentricity of the inner pair.It is the latter which is the main source of that bump in the ETVs, which now looks evident in the TESS measured eclipse times of KIC 9714358.
Here we note, that KIC 9714358 is scheduled to be reobserved by TESS in Sectors 74, 75 (January -February, 2024) and 81 (July/August, 2024), when the offset between the primary and secondary eclipses is expected to be near the next maximum (see upper right panel of Fig. 1).Thus, these findings and the corresponding predictions of our model will be verifiable very soon.

KIC 5771589
As mentioned above in Sect.3, the two noteworthy observables of this triple system are the very rapid apsidal motion and the very quickly varying eclipse depths.Despite our best efforts, we were unable to model the variations of the eclipse depths perfectly.In general, as one can see in Fig. 4, LIGHTCURVEFACTORY returns the main properties of the EDVs in the case of both the Kepler and the TESS data.The problem, however, shows up in the second half of the Kepler data where our models fail to reproduce the correct depth ratios of the secondary vs primary eclipses.In other words, while the Kepler observations suggest that the ratio of the depths of the secondary and primary eclipses remains near constant during the four years of the observations, our model fits tend to equalize the depths of the two kinds of eclipses toward the end of the Kepler observations.
To understand the origin of this problem, one should keep in mind that for an eccentric (e ̸ = 0) EB viewed not exactly edge-on (i ̸ = 90 • ), the ratio of the eclipse depths, in addition to being a function of the ratio of the surface brightnesses, also depends strongly upon the position of the observable argument of periastron.This is so because, for the eclipse that occurs closer to periastron, the projected separation of the two stellar disks will be smaller than for the other kind of eclipse (closer to the apastron).Thus, a larger portion of the surface of the eclipsed star will be occulted.Naturally, the more eccentric the orbit, and the more grazing the eclipse, the more significant this effect may be.And, naturally, one can reach a situation where the eclipse closer to apastron disappears entirely, as was observed recently e.g., in the case of KIC 5731312 [3].Note, we illustrate graphically these effects for a hypothetical eccentric EB in Fig. 7.
In the case of the current system, both the position of the argument of periastron and the inclination have varied during the four years of Kepler observations.As we mentioned above, but it can also be seen directly from the ETV curve (middle left panel of Fig. 1), the inner orbital ellipse made a bit more than half revolution during the Kepler observations.According to our photodynamical model, the observable argument of periastron attained a value of (ω = 270 • ) around the very beginning of the Year 2010 (circa 2 455 200, see also the left panel of Fig. 8).This means that during that time, the primary eclipses occurred during the periastron passages, while the secondaries occurred at the apastron points, thereby causing the largest differences in the sizes of the eclipsed areas (bottom row of Fig. 7).Then, toward the end of the Kepler era, the situation reversed.During the last days of the perfectly operating three gyroscopes on the Kepler spacecraft, the observable argument of periastron reached ω = 90 • , i.e., the last Kepler-observed primary eclipses occurred around the apastron points, while the secondaries were near periastron (second row of Fig. 7).This is the reason why the secondary eclipse depths tend to be much closer to the primary ones by the end of the Kepler data in our photodynamical models.We note, however, that the same argument can be made purely from the observations, without the use of our photodynamical model results and, hence, one cannot claim this discrepancy as a failure of the model.This is so, because the ETV curves show in themselves that the  , respectively.The left panels represent the synthetic light curves, which in all four cases are phased to the inferior conjunction of the less massive secondary, while the right panels show the corresponding geometry, as seen from the Earth, and also as it would be seen from the pole of the orbit (small, inserted figures).(In the case of the insert plots, the direction of the Earthly observer is along the negative Y axis.Note also that the secondary revolves counter-clockwise in all the orbit plots.)As one can see on the light curve plots, the depth ratio of the primary and secondary eclipses depend strongly on the orientation of the ellipse.When the ellipse seen along the minor axis (ω = 0 • , 180 • ) the eclipsed disk areas are similar during the inferior and superior conjunctions and, hence, the ratio of the eclipse depths -similar to the circular orbit case -primarily depend only upon the ratio of the surface brightnesses.On the other hand, when the ellipse is seen along its major axis (ω = 90 • and 270 • ), the eclipsed surface area is much larger during the periastron-eclipse than during the apastron-eclipse.In the current illustration, due to the relatively large eccentricity of e = 0.3, the difference in the eclipsed surface areas is so large that the eclipse depths reverse.system was seen from the direction of the major axis at the very beginning and the end of the Kepler observations.(This can be deduced from the fact, that the primary and secondary curves intersect each other only when the orbital ellipse is seen from the direction of the major axis, i.e., when one of the eclipses occurs at periastron, while the other at apastron.)Then, the fact, that the secondary (blue) ETV curve is located below the primary (red) one during the entire Kepler observations reveals that during this interval, the secondary eclipses were "in a hurry" in contrast to the primary ones (or, in other words, the secondary eclipses occurred before photometric phase of 0.5).This means that the periastron passage happened in between the primary and the secondary eclipses (i.e., ω had values in between 270 • [or, −90 • ] and 90 • , reaching ω = 0 • at the mid-time -see the top row of Fig. 7).But, even if one would assume retrograde apsidal motion, which would mean that the role of the apastrons and periastrons were exchanged, this cannot explain why the observed ratios of the secondary to primary eclipse depths would have immunity to the apsidal revolution.
In our understanding, the only possibility for explaining the near constant depth ratios during half an apsidal cycle is that the true observable inclination of the system should be closer to i = 90 • than our findings of i in = 86.05• ± 0.09 • would suggest (since the more edge-on an eccentric EB is viewed, the smaller the apsidal-phase-dependent difference of the eclipsed disk areas and, for orbits viewed perfectly edge-on, this effect entirely disappears).This fact, however, would require more contaminating light in the system.The main possibilities to fulfill this requirement are as follows: (i) the system has a fourth, more distant, undetected component; (ii) the tertiary is more luminous than would be inferred from the PARSEC stellar isochrones that we utilized or, (iii) the outer mass ratio (q out ) could be larger and, hence, the tertiary would be more massive (and, consequently brighter) relative to the binary members than our solution suggests.Regarding this last possibility, we made efforts to find acceptable solution with higher q out but all of our trials failed.So, in what follows, we discuss our findings in the imperfect case, where the (non-)variation of the eclipse depth ratios is imperfectly modelled.Despite this fact, however, we believe that our findings, to be discussed below, are basically correct.
According to our joint photodynamical solution, in KIC 5771589 the most massive star is the third, distant component, with m B = 1.02 ± 0.01 M ⊙ , though, the primary of the inner binary is only slightly less massive (m Aa = 0.94 ± 0.01 M ⊙ ), while its binary companion has a lower mass of m Ab = 0.73 ± 0.01 M ⊙ .Here we note that, despite the lack of RV data, our solutions yield masses with surprisingly small uncertainties of only 1 − 2%.Due to the problematic nature of the light curve solution discussed above, however, we are not convinced that these small uncertainties are realistic.The quantity which is certainly more robust is the outer mass ratio, being q out = 0.612 ± 0.003.This value is substantially lower than what was inferred from the previous analytic ETV fits of Borkovits et al. [11,12].One should keep in mind, however, that those former fits were based only on the Kepler data (long before the launch of the TESS space telescope) and, because of the extremely tight nature of the current system, as the authors noted, the analytical ETV models they used were somewhat inadequate.
Our solution suggests a metal deficient ([M/H] ≈ −0.50), old (log τ = 9.83 ± 0.02), and slightly evolved system.According to the accepted solution, the massive tertiary is clearly on its way toward the giant branch, having log g B = 3.7 ± 0.1.One should keep in mind, however, that in the absence of third-body eclipses, the direct effects of the tertiary on the joint, complex solutions (apart from the SED fitting) manifest themselves only in the outer mass ratio (via timings) and the contaminated light (via the eclipse depths).Hence, any further statements about the astrophysical condition of the tertiary depend strongly upon the pre-assumption that the whole system formed and evolved in a coeval manner.Our solution gives a distance of d = 709 +33 −55 pc, which differs substantially from the Gaia EDR3-based value of d EDR3 = 870 ± 100 pc.This fact might again suggest that the system could have an additional, fourth stellar component (not considered in the photodynamical solution), and the high value of the RUWE parameter (11.34) indeed indicates that the Gaia astrometric solution is probably significantly influenced by the multiplicity of the system.
Turning to the readily observable long-term or, secular perturbations in this triple, the rapid, dynamically forced apsidal motion is evident.The doubly averaged secular theory predicts apsidal motion periods of the inner and outer orbits to be (P apse = 11.22 ± 0.04 and 51.7 ± 0.2 yrs, respectively).From the fact that the ∼ 4-yr-long Kepler dataset covers a bit more than half of an apsidal cycle, one can infer directly that the quadrupole model, again, gives an incorrect result for this system.This can also be seen in the left panel of Fig. 8 which shows that the true apsidal motion period is about P meas apse ≈ 6.3 yr.On the other hand, however, note that the outer apsidal motion period according to our numerical integration is P meas apse ≈ 58.1 yr which is, again, in much better agreement with the theoretical quadrupole value.
The significance of the octupole order perturbations in the variations of the inner eccentricity is, again, nicely demonstrated in the right panel of Fig. 8.As one can see, the cyclic variations of e in in between ∼ 0.01 and ∼ 0.02 have exactly the same period as for the quantity ω dyn in − ω dyn out , which, in the current system varies between 0 • and 360 • with a period of ∼ 8 yrs.
The most interesting feature of KIC 5771589, however, is its continuous eclipse depth variations.What may be surprising, therefore, at first sight is the fact that according to our results the inner and outer orbital planes are almost coplanar.The non-zero mutual inclination, which is the origin of the nodal precession and, hence, the eclipse depth variations is very small, being only i mut = 0.29 • ± 0.03 • .The question naturally arises as to how is it possible that such a small mutual inclination can cause such a readily observable eclipse depth variation?The answer lies in the extreme grazing nature of the eclipses, as is shown in the left panel of Fig. 9.
Finally, in Fig. 10 we plot the predicted eclipse depth variations for the present century.As one can see, for the expected near-future observations of TESS (in Sectors 74, 80, 81) the inner inclination will be larger and, hence, one can predict somewhat deeper and better detectable eclipses.The numerical integrations give a nodal regression period of P meas node ≈ 9.2 yr which is in perfect accord with the theoretical value tabulated in Table 3.

TIC 219885468
Our third triple is a good counterexample in regard to both (i) the octupole perturbations driving apsidal motion and eccentricity variations, and (ii) the remarkable eclipse depth variations in an almost flat system.
Both the tightness and the compactness of this system is in between our previous two example triples.Moreover, similar to the other two triples, the inner orbit has a small eccentricity (e in = 0.0423 ± 0.0006), while the outer orbit is moderately eccentric (e out = 0.3903 ± 0.0007).The outer mass ratio, again, makes this triple quite similar to KIC 9714358 (q out = 0.319 ± 0.002 vs. q out = 0.265 ± 0.003 for the TIC and the KIC systems, respectively), while the mutual inclination is almost identical to that of the EDV system Variations of the observable inclinations of the inner and outer orbital planes (i in -red; i out -blue, respectively) of KIC 5771589 (left) and TIC 219885468 (right).The green area represents the domain where the inner binary produces at least one eclipse during a revolution, while the cyan area covers that inclination domain where both eclipses are present.Furthermore, the nearly horizontal.black lines denote the lower and upper limits of that domain of the outer inclination (i out ), where third-body eclipses might occur.As one can see, in the case of KIC 5771589, the inner inclination (i in ) shows small amplitude oscillations near the border of the EB domain, implying that the eclipses are grazing.Moreover, the blue curve (i out ) remains continuously below the third-body eclipse inclination limit, hence, no third-body eclipses can occur.In contrast to this, in the case of TIC 219885468 the small amplitude oscillations of i in are located quite deeply in the (cyan/green) EB domain, indicating that the eclipses are quite deep and, hence, the eclipse depths are insensitive to such small inclination variations.Note also that this is an almost triply eclipsing system, as the outer inclination is just below the triply eclipsing limit.KIC 5771589 (i mut = 0.25 • ± 0.14 • vs. i mut = 0.29 • ± 0.03 • , respectively).The similarities, however, end at this point.From a dynamical point of view, the main difference between the present and the two KIC systems is that the inner EB of TIC 219885468 consists of two twin stars (q in = 0.988 ± 0.017) and, hence, the octupole (and any other even order) perturbation terms can be neglected.This is demonstrated nicely in the two panels of Fig. 11, where it can be readily seen that both the inner and outer ellipses rotate with constant rates, as is expected for coplanar systems in the quadrupole theory.Moreover, the theoretically calculated quadrupole apsidal advance rates (P theo apse,in = 19.49± 0.09 yr and P theo apse,out = 51.2 ± 0.1 yr) agree quite well with the apsidal motion periods 'measured' via numerical integration (P meas apse,in ≈ 13.6 yr and P meas apse,out ≈ 54.5 yr).The agreement is even better in the nodal regression period being P meas node ≈ 13.9 yr vs. P theo node = 14.09 ± 0.07 yr.Turning to this latter effect, despite the similar amplitude of the precession cone, in the case of TIC 219885468 one cannot detect any EDVs.The two panels of Fig. 9 illustrate nicely the substantial differences amongst the eclipse geometries of the two EBs.While in the case of KIC 5771578 there are grazing eclipses, in TIC 219885468 the eclipses are much deeper.Hence, the areas of the eclipsed surfaces are less sensitive to such small inclination variations.
Regarding the astrophysical parameters of the three stars, the twin EB members are found to be late F/early G-type, slightly evolved MS stars (log τ = 9.49 ± 0.

Conclusions
We investigated three such EBs which form the inner pairs of TCHT stellar systems.Two systems (KICs 9714358 and 5771589) were observed with the Kepler space telescope in its 4-yr-long prime mission, while the third target (TIC 219885468) is located in the NCVZ of TESS.Besides the evident medium-period third-body perturbations, the ETVs of all three EBs show rapid dynamically driven apsidal motion, and one of the three targets (KIC 5771589) also exhibits cyclic EDVs.We carried out complex photodynamical analyses for the three systems analyzing jointly their Kepler and TESS light curves (naturally, for TIC 219885468 only the TESS light curve was used), ETV curves, composite SEDs and, in the case of KIC 9714358, APOGEE-2 RV data as well.We determined reliable, robust astrophysical and orbital parameters for all three systems and their constituent stars.
In our study, however, we focused mainly on the dynamical properties, and higher order third-body perturbations of the systems.We found clear evidence that in the case of the almost perfectly coplanar KIC 9714358 the inner and outer orbital ellipses are oriented in the same directions and while, on average, they rotate evenly with the same apsidal revolution rate (with a period of P meas apse ≈ 78 yr), the axis of the inner orbit librates around this average direction with an oscillation period of P osc apse,in ≈ 30 yr and with an amplitude of ∼ 15 • .The variations of the inner eccentricity of KIC 9714358 show the same period as the apsidal oscillations, but with a phase shift of ∼ 0.25 with respect to the libration.We connect this behaviour to the octupole perturbations; however, we leave the quantitative investigation of this issue for a future study.We expect, however, that the near-future TESS observations will provide us with additional observational evidence for this behaviour.
Regarding KIC 5771589, despite the evident cyclic EDVs of this system (detected both in the Kepler and TESS eras), we found, somewhat surprisingly, that the triple was almost coplanar, with a mutual inclination of i mut = 0.29 • ± 0.03 • .We explain this unusual finding with the grazing nature of the eclipses, as the depths of grazing eclipses are much more sensitive to even very small variations in the observable inclination.This triple is the second tightest Kepler-triple, and also this EB has the second shortest apsidal motion period.Despite this fact, we did not find observational evidence for the operation of the octupole-order perturbations, however, our numerical integration did reveal that they are present.The weakness of the octupole effects in the present system could be in accord with the fact that the relative strength of the octupole order effects relative to the quadrupole ones is much lower than in KIC 9714358.
Finally, TIC 219885468 was selected as a counterexample to the two KIC sources.This triple, in both compactness and tightness, is very similar to the other two systems, but the nearly equal primary and secondary eclipse depths suggested that the two stars of the inner EB should be similar in mass and, hence, one cannot expect octupole order third-body perturbations in this system.While our analysis justified this pre-assumption, finding q in = 0.988 ± 0.017 revealed another aspect in which TIC 219885468 serves as counterexample.The mutual inclination of the system was found to be almost identical to that of the former system (i mut = 0.25 • ± 0.14 • vs. 0.29 • ± 0.03 • ).In the current system, however, the eclipses are deep and, hence, they are insensitive to such small variations in inclination.
As our primary aim was to detect observationally the higher order secular perturbations, we did not investigate the past and future evolution of these three triple systems as this question is beyond the focus of the present paper.Currently all three systems are stable in the sense of the semi-empirical dynamical stability limit of Mardling & Aarseth [6].This fact, however, does not imply that the systems also remain stable in the distant future.Such investigations may also be the part of a future study.

Figure 1 .
Figure 1.Observed and modelled ETVs of KICs 9714358 (upper), 5771578 (middle) and TIC 219885468 (bottom).The left-hand panels cover the time interval approximately from the beginning of the Kepler observations to the end of the second northern ecliptic hemisphere scan of TESS, while the right-hand panels show longer, 100 yr-long intervals, which offer direct predictions for future observations.Larger red and blue symbols represent primary and secondary ETV points determined from the observations, while the smaller points connected with straight lines stand for ETV points calculated according to the best-fitting complex photodynamical models.Finally, the black, almost horizontal lines represent the LTTE contributions to the curves.Note, in the middle right panel we downshifted the originally overlapping secondary ETV curve of KIC 5771589 by 0.2 days, for a better view of the small effect of variations in the inner eccentricity forced by the octupole perturbations.See text for further details.

Figure 2 .
Figure 2. Some satellite-observed primary and secondary eclipses of KIC 9714358 (blue dots) with the best-fitting photodynamical model (red curves).Left panel: Consecutive primary and secondary eclipses from the beginning, mid-time, and near to the end of the prime Kepler mission.Right panel: Characteristic primary and secondary eclipses from the 2019 and 2022 visits of TESS.Each of the sections is 0.4 day long.

Figure 3 .
Figure 3. APOGEE-2 RVs of the primary component of KIC 9714358 (red dots) together with the best photodynamical model (blue curves).Upper panel: A characteristic section of the observed RV points and the corresponding model fit in the time domain.Both the inner and the outer orbits are well visible.Lower left panel: Phase folded RV data according to the inner period, after the removal of the signals of the outer orbit.Lower right panel: Phase folded RV data according to the outer period, after the removal of the effects of the inner orbital motion.The varying shape of the folded solution, caused primarily by the rapid apsidal motion of the outer orbit is readily visible.Note, the thin horizontal black lines at V = −66.91km s −1 denote the systemic radial velocity of the triple star.

Figure 4 .
Figure 4. Upper panel: The ∼ 4-year-long Kepler light curve of KIC 5771589 (blue dots) with the bestfitted photodynamical model (red).Alternating gray-white stripes denote the consecutive observing quarters.Note, Q5, Q9 and Q13 data are missing, as in these quarters the stellar light had fallen onto a bad CCD camera.Lower left panel: Consecutive primary and secondary eclipses from the beginning, mid-time, and near to the end of the prime Kepler mission.Lower right panel: Characteristic primary and secondary eclipses from the 2019, 2021 and 2022 visits of TESS.Each of the sections is 1 day long.

Figure 5 .
Figure 5.Primary and secondary eclipses from the beginning of Year 2 and Year 4, as well as from the end of Year 5 TESS-observations of TIC 219885468.Each of the sections is 1 day long.

Figure 6 .
Figure 6.The evolution of the orbits in KIC 9714358 from its discovery to the end of the current century.Left panel: The variations of the observable and the dynamical arguments of periastron of the inner and outer orbits (ω in,out and ω dyn in,out , respectively).Right panel: The cyclic variation of the inner eccentricity, and the very similar variations of the difference of the dynamical arguments of periastrons of the two orbits.The gray shaded area represents the interval of the Kepler observations.See text for further details.

Figure 7 .
Figure 7. Illustration of the effect of the orientation of the orbit on the properties of the eclipses.For illustration purposes we generated light curves for a hypothetical binary with parameters: m Aa = 1.00 M ⊙ , m Ab = 0.80 M ⊙ ; R Aa = 1.03 R ⊙ , m Ab = 0.799 R ⊙ ; T Aa = 5878 K, T Ab = 5274 K; P = 5.00 d; e = 0.30; i = 85 • .The only difference amongst the four rows is the orientation of the apsidal line of the ellipse, which are ω = 0 • , 90 • , 180 • and 270• , respectively.The left panels represent the synthetic light curves, which in all four cases are phased to the inferior conjunction of the less massive secondary, while the right panels show the corresponding geometry, as seen from the Earth, and also as it would be seen from the pole of the orbit (small, inserted figures).(In the case of the insert plots, the direction of the Earthly observer is along the negative Y axis.Note also that the secondary revolves counter-clockwise in all the orbit plots.)As one can see on the light curve plots, the depth ratio of the primary and secondary eclipses depend strongly on the orientation of the ellipse.When the ellipse seen along the minor axis (ω = 0 • , 180 • ) the eclipsed disk areas are similar during the inferior and superior conjunctions and, hence, the ratio of the eclipse depths -similar to the circular orbit case -primarily depend only upon the ratio of the surface brightnesses.On the other hand, when the ellipse is seen along its major axis (ω = 90 • and 270 • ), the eclipsed surface area is much larger during the periastron-eclipse than during the apastron-eclipse.In the current illustration, due to the relatively large eccentricity of e = 0.3, the difference in the eclipsed surface areas is so large that the eclipse depths reverse.

Figure 8 .
Figure 8.The same plots as in Fig. 6, but for KIC 5771589.See text for further details.

Figure 9 .
Figure 9. Variations of the observable inclinations of the inner and outer orbital planes (i in -red; i out -blue, respectively) of KIC 5771589 (left) and TIC 219885468 (right).The green area represents the domain where the inner binary produces at least one eclipse during a revolution, while the cyan area covers that inclination domain where both eclipses are present.Furthermore, the nearly horizontal.black lines denote the lower and upper limits of that domain of the outer inclination (i out ), where third-body eclipses might occur.As one can see, in the case of KIC 5771589, the inner inclination (i in ) shows small amplitude oscillations near the border of the EB domain, implying that the eclipses are grazing.Moreover, the blue curve (i out ) remains continuously below the third-body eclipse inclination limit, hence, no third-body eclipses can occur.In contrast to this, in the case of TIC 219885468 the small amplitude oscillations of i in are located quite deeply in the (cyan/green) EB domain, indicating that the eclipses are quite deep and, hence, the eclipse depths are insensitive to such small inclination variations.Note also that this is an almost triply eclipsing system, as the outer inclination is just below the triply eclipsing limit.

Figure 10 .
Figure 10.Simulated light curve of KIC 5771589 in Kepler's photometric band for all of the 21st century.The gray, shaded area, as before, represents the interval of the Kepler observations, while the other vertical lines stand for the past and scheduled future observations of TESS.

Figure 11 .
Figure 11.The same plots as in Fig. 6, but for TIC 219885468.The filled brown area represents the interval in between the beginning of the very first (S14) and the end of the momentarily last (S60) NCVZ observations of TESS.See text for further details.
06) with masses of m Aa = 1.19 ± 0.03 M ⊙ and m Ab = 1.18 ± 0.04 M ⊙ and, effective temperatures of T Aa = 6 410 ± 40 K and T Ab = 6 400 ± 40 K.The less massive third component is a K dwarf with m B = 0.75 ± 0.02 M ⊙ and T B = 4 885 ± 50 K.The inferred distance of the triple star agrees with the Gaia EDR3 distance within two sigma (d = 1 150 ± 20 pc vs d EDR3 = 1 111 ± 13 pc).

Table 1 .
Main properties of the three systems from different catalogs

Table 2 .
Orbital and astrophysical parameters of KIC 97143586 from the joint photodynamical lightcurve, ETV, RV, SED and PARSEC isochrone solution.Meanings of most of the parameters are explained in the text, with the exceptions of i inv and Ω inv , which quantities give the position of the invariable plane with respect to the tangent plane of the sky.The osculating orbital elements are given for epoch t 0 = 2 454 953.0.

Table 3 .
Orbital and astrophysical parameters of KIC 5771589 from the joint photodynamical lightcurve, ETV, SED and PARSEC isochrone solution.The osculating orbital elements are given for epoch t 0 = 2 454 953.0.

Table 4 .
Orbital and astrophysical parameters of TIC 219885468 from the joint photodynamical lightcurve, ETV, SED and PARSEC isochrone solution.The osculating orbital elements are given for epoch t 0 = 2 458 683.0.