1. Introduction
The new window that 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 inevitable primacy in the discovery and principal characterizations of this rare new class of stellar systems, the multiple 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 also 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.,
), while we call a triple or multiple star system ‘compact’ when the outer period does not exceed, let us say, 1000 days
1. As is widely 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 the case because the magnitude of the gravitational perturbations relative to the pure Keplerian two-body motions primarily depends on the ratio of the semi-major axes or, more strictly, on 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,
8,
9] and references therein) but in regard to 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 practical significance. This significance has two somewhat counteracting aspects. The first aspect is connected to the largest amplitude, and most interesting, orbital perturbations in tight triple star systems, 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
, 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. Considering the second aspect, 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) that form the inner binary of a tight hierarchical triple system. The amplitude of such medium-timescale dynamically driven ETVs, however, scales as
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 four-year-long observations of the prime
Kepler mission. The vast majority of the EBs known before
Kepler had periods of 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
2. The quasi-continuous, almost four-year-long observations of the
Kepler space telescope have 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 (
-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 that is not only tight but compact as well. While these former studies depend primarily on the
-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 of the nodal regression) was included in 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
Section 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
Section 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
Section 4 and
Section 5, while the results are discussed in
Section 6.
2. 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 comprising the Earth and an artificial satellite due to the Moon [
18].
To our knowledge, Slavenas [
19,
20] was the first to investigate 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 the scope of the scientific community. The situation changed 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 1990s. 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 , and the amplitude is related to ,
medium-period perturbations, with a characteristic period of and an amplitude of ,
long-period perturbations3 having a period of about
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 [
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
; hence, it is called the ‘quadrupole approximation’). At that level, the problem has only one degree of freedom, and thus it is integrable; the solution can be given with the use of elliptic functions of Weierstrass
. 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 behavior 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 to the results of direct numerical integrations) for different mutual inclinations (from the coplanar case up to
). 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 long-period perturbations exist 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
and
denote the eccentricities and dynamical arguments of periastron of the inner and outer orbits, respectively. These latter quantities (sometimes denoted as the Delaunay variables
), which should not be confused with the observable arguments of periastrons (the latter of which will hereafter be denoted as
), give the angular distances of the pericenter points of the orbits from the ascending nodes of the intersections of the orbital planes. Furthermore,
stands for the cosine of the mutual (or, relative) inclination of the two orbital planes. Finally,
where we have introduced the parameters
and
, 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 , ; the total mass of the inner pair, ; the individual mass of the third component, ; and the total mass of the triple system, . Moreover, the mass ratios are denoted by for the inner and outer binaries, respectively.
In regard to the amplitude-like quantity of
, as one can see, it disappears for
, i.e., when the inner pair consists of two equal mass stars
4. It 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 (
), the dynamical arguments of pericenters (
—or, traditionally,
), and the dynamical longitudes of the nodes (
—or,
) of the inner and outer orbits, respectively. One should also keep in mind that, as was noted by Naoz et al. [
36], although Equation (
1) formally does not contain the nodes (
), which would suggest the constancy of the conjugate variables
, 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
and, therefore,
, and
. If one substitutes these values into the original Hamiltonian, it leads to the usual form of Equation (
1). Hence, strictly speaking, from the simplified (Equation (
1)) form of the Hamiltonian, the constancy of
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 of interest to us take the following form:
and, regarding the perturbations in the observable arguments of periastron:
In what follows, we consider only the prograde, coplanar scenario, i.e., when
. This restriction can be justified by the fact that, as will be shown in
Section 6, the mutual inclination in none of the three considered systems exceeds
. In the coplanar case, when
, Equations (
12)–(
17) 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 (
), as well as the angle between the directions of the two periastron points, i.e.,
, where in the last expression we took into account the fact that
. Furthermore, we notice that in the coplanar case
. 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,
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 the 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; hence, one should 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 is 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; to determine the 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.
3. 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; 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.0 coordinates 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 (
), photogeometric distances, metallicities ([
]), interstellar reddenings (
), and proper motions (
) from different catalogs. Finally, we also give the specific
(renormailzed unit weight error) parameters, which were 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 briefly discuss the available basic information on these three TCHTs.
KIC 9714359 = KOI-6073 (
d;
d;
,
) is the compact, tight triple star system that 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], who analyzed 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 that 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, a moderately eccentric outer orbit (
and
, respectively), and a dynamically forced apsidal motion period of
yr. The model also indicated inner and outer mass ratios of
and
, 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
from the
Kepler light curves, and, moreover, based on astrophysical reasons, they estimated a lower inner mass ratio (
). 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
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 (
), and, correspondingly, they found that the system was in a pre-MS state with an age of
.
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
Figure 1, the new ETV points, determined from
TESS-observed eclipses, do not follow the former trend, and the long-term behavior 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 an (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
Section 6.
Finally, note that this is the only system in our sample for which radial velocity (RV) data are available from the literature. The data release 16 of the APOGEE-2 project [
53] tabulates 16 individual RV observations, which we used for our analysis (see
Figure 2).
KIC 5771589 = KOI-6625 (
d;
d;
,
). 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
5. The analytical ETV studies suggested a large outer mass ratio (
), 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
Figure 3), 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 that 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.,
, 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.,
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
6. In such a way, the current light curve behavior 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 between the Year 2 and Year 4
TESS observations (see the lower right panel of
Figure 3).
TIC 219885468 has d; d; ; and . 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 differs 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
Figure 1). The primary and secondary eclipses look very similar in depth, which indicates that the inner mass ratio should be near unity—and this nulls 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
Figure 4).
5. Photodynamical Modeling
We have carried out complex photodynamical analyses of our three targets with our own software package
Lightcurvefactory, developed by one of us (T. Borkovits, 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 s, 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
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 model-independent. 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 , the primary of the inner pair, or , the tertiary star), and the inner and outer mass ratios (). Moreover, the metallicity of the system ([]), the (logarithmic) age of the three coeval stars (), and the interstellar reddening for the given triple were also varied. Additionally, the ‘extra light’ contamination parameter was also freely adjusted for two of the three systems.
- (ii)
Orbits: Three of six orbital-element related parameters of the inner orbits, and six parameters of the outer orbits, i.e., the eccentricity vector components of both orbits , ; the inclinations relative to the plane of the sky (); and three other parameters for the outer orbit, including the period (), the longitude of the node relative to the inner binary’s node (), and the periastron passage times (), 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, and 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 .
- (ii)
The atmospheric parameters of the stars: we handled them in a similar manner as in our previous photodynamical studies. We utilized a logarithmic limb-darkening law [
65] for which the passband-dependent linear and non-linear coefficients were interpolated at each trial step via the tables from the 0.xx versions of the
Phoebe software, developed and maintained by A. Prša [
66]. We set the gravity darkening exponents for all late type stars to
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 () 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—) were constrained internally through the ETV curves. Finally, in the case of KIC 9714358, the systemic radial velocity () 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 .
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
uncertainties, are tabulated in
Table 2,
Table 3 and
Table 4. Furthermore, the observed vs. model lightcurves are plotted in
Figure 3,
Figure 4 and
Figure 5, while the observed vs. model ETV curves are shown in
Figure 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
Figure 2.
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 (
) and the dynamical (
) frame of references. For their calculation, we do not exclusively use 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 (
). Moreover, besides the apsidal motion parameters, we also give the theoretical, quadrupole, nodal regression period (
), i.e., the time needed for a
regression of the parameters
. Further details of these derivations are discussed in Section 6.2 of Kostov et al. [
68].
6. Discussion
In what follows, we discuss our results for the three investigated systems individually. Though we are primarily interested in the inferred dynamical behavior of our triples, first we briefly consider the astrophysical implications of our results and then discuss the dynamics of the systems in detail.
6.1. 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
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 (
) 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 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 directly contains only the spectroscopic mass function, which can easily be written in the following form:
This, coupled with the known (inner) inclination of and the primary’s mass (at least to within a few percent uncertainty, as in the current situation) of , uniquely gives the (inner) mass ratio . Thus, this very robust mass ratio clearly excludes any post-MS solutions; 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 more clearly supports a hotter primary component with
K than the one listed in either the TIC v8.2 catalog (
K, see in
Table 1) or inferred from the APOGEE-2 spectra (
K, Jönsson et al. [
53]). In this context, some further caution is needed because our inferred (photometric) distance was found to be
pc, which substantially exceeds the Gaia EDR3 value of
pc [
47]. One should keep in mind, however, that Gaia EDR3 parallaxes; 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 ) and the eccentricities of the two orbits ( and ), 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 . (It will be shown in the next subsection that even a few tenths of a degree departure from exact coplanarity may lead to the robust detection of EDVs for such a tight system).
The most interesting dynamical feature of the current system is, however, the unusual behavior 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
yr (see the middle of
Table 2). As one can readily see in the upper left panel of
Figure 1, the ∼5000-day-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
Figure 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
yr (which is quite close to the theoretical quadrupole outer apsidal motion period,
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
and with a half-amplitude of ∼15
. (In other words, the difference of
–
oscillates around
). Naturally, the same is true for the dynamical arguments of periastrons, with the difference here being that
–
oscillates around
, as is shown nicely in the right panel of
Figure 6.
In the right panel of
Figure 6, besides the difference of
–
, we also plot the variations of the inner eccentricity (
) at the same times. As one can see, the inner eccentricity oscillates with the very same period and with a
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 Equation (
18)). 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 behavior 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 min. 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
11. On 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 will not change the eccentricity of the inner pair. It is the latter that 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
Figure 1). Thus, these findings and the corresponding predictions of our model will be verifiable very soon.
6.2. KIC 5771589
As mentioned above in
Section 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
Figure 3,
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 (
) EB viewed not exactly edge-on (
), 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 graphically illustrate these effects for a hypothetical eccentric EB in
Figure 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, which can also be seen directly from the ETV curve (middle left panel of
Figure 1), the inner orbital ellipse made slightly more than half a revolution during the
Kepler observations. According to our photodynamical model, the observable argument of periastron attained a value of (
) around the very beginning of the Year 2010 (circa 2,455,200; see also the left panel of
Figure 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
Figure 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
, i.e., the last
Kepler-observed primary eclipses occurred around the apastron points, while the secondaries were near periastron (second row of
Figure 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; hence, one cannot claim this discrepancy as a failure of the model. This is so because the ETV curves show in themselves that the 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 and 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 between the primary and the secondary eclipses (i.e.,
had values between
[or,
] and
, reaching
at the mid-time—see the top row of
Figure 7). But, even if one assumed retrograde apsidal motion, which would mean that the roles of the apastrons and periastrons were exchanged, this would not 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 than our findings of 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 () could be larger; 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 an acceptable solution with higher , 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 modeled. 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
, though the primary of the inner binary is only slightly less massive (
), while its binary companion has a lower mass of
. 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 that is certainly more robust is the outer mass ratio, being
. 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 (), old (), and slightly evolved system. According to the accepted solution, the massive tertiary is clearly on its way toward the giant branch, having . 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 pc, which differs substantially from the Gaia EDR3-based value of 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 (
and
yrs, respectively). From the fact that the ∼4-yr-long
Kepler dataset covers slightly 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
Figure 8, which shows that the true apsidal motion period is about
yr. On the other hand, however, note that the outer apsidal motion period according to our numerical integration is
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
Figure 8. As one can see, the cyclic variations of
between ∼0.01 and ∼0.02 have exactly the same period as for the quantity
–
, which, in the current system, varies between
and
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 is the origin of the nodal precession; hence, the eclipse depth variations are very small, being only
. 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
Figure 9.
Finally, in
Figure 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; hence, one can predict somewhat deeper and better detectable eclipses. The numerical integrations give a nodal regression period of
yr, which is in perfect accord with the theoretical value tabulated in
Table 3.
6.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 between our previous two example triples. Moreover, similar to the other two triples, the inner orbit has small eccentricity (
), while the outer orbit is moderately eccentric (
). The outer mass ratio, again, makes this triple quite similar to KIC 9714358 (
vs.
for the TIC and the KIC systems, respectively), while the mutual inclination is almost identical to that of the EDV system KIC 5771589 (
vs.
, 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 (
); hence, the octupole (and any other even order) perturbation terms can be neglected. This is demonstrated nicely in the two panels of
Figure 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 (
yr and
yr) agree quite well with the apsidal motion periods ‘measured’ via numerical integration (
yr and
yr).
The agreement is even better in the nodal regression period, being
yr vs.
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
Figure 9 nicely illustrate 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 () with masses of and and effective temperatures of K and K. The less massive third component is a K dwarf with and K. The inferred distance of the triple star agrees with the Gaia EDR3 distance within two sigma ( pc vs. pc).
7. Conclusions
We investigated three such EBs that 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, jointly analyzing 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 yr), the axis of the inner orbit librates around this average direction with an oscillation period of 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 behavior 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 behavior.
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 . 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 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 revealed 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; hence, one cannot expect octupole-order third-body perturbations in this system. While our analysis justified this pre-assumption, finding 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 ( vs. ). In the current system, however, the eclipses are deep; hence, they are insensitive to such small variations in inclination.
As our primary aim was to observationally detect 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 part of a future study.