© The Authors 2022

1. Introduction

Among the 5000 stars that are visible with the naked eye, about 2000 are known to be multiple-star systems. Naked-eye stars account for a small fraction of the stars in the Milky Way, but are reasonably representative of the incidence of binarity, which is estimated to be between 50 and almost 100% (e.g., Eggleton 2006). Some systems are close enough to be in contact, while others are far apart enough to evolve almost independently. Binary systems are known to have orbital periods as short as 0.2 days or as long as thousands of years. Studying multiple-star systems is thus important in order to understand the evolution of our Galaxy.

Stars in multiple systems are also precious benchmarks for calibrating asteroseismology when it is possible to measure the mass of the oscillating star independently from its oscillation properties. This is the case with double-lined spectroscopic binary (SB2) where the components eclipse each other, or optically resolved binary systems. Hitherto, all of the solar-like oscillators belonging to eclipsing binary stars (EBs) are red giants (RGs), and all have been detected by the Kepler mission (Hekker et al. 2010; Gaulme et al. 2013, 2014; Beck et al. 2014, 2015; Kuszlewicz et al. 2019; Benbakoura et al. 2021). So far, 14 wide SB2 EBs including an oscillating RG have been fully characterized with the help of ground-based radial-velocity support (Frandsen et al. 2013; Rawls et al. 2016; Gaulme et al. 2016; Brogaard et al. 2018; Themeßl et al. 2018; Benbakoura et al. 2021). Three more bona fide oscillating RGs in EBs were reported by Gaulme & Guzik (2019), and are currently under study. Beyond EBs, it is also possible to determine the masses of stars belonging to visual multiple-systems, where individual components are spatially resolved, which allows their projected orbits to be retrieved, provided radial velocities (RVs) are available as well. So far, very few of these systems have included solar-like pulsators (Marcadon et al. 2018; Metcalfe et al. 2020). The drawback of such systems is that the orbits are long, and getting accurate orbital parameters can take decades and rely on heterogeneous datasets. Finally, hierarchical triple systems are promising types of benchmarks; they are composed of a close binary with a relatively distant companion (e.g., Tokovinin 1997; Ford et al. 2000; Borkovits et al. 2003). Depending on their orbital configurations, the presence of eclipses, and significant eclipse timing variations (ETVs), it is possible to determine the masses of their components (Derekas et al. 2011; Carter et al. 2011; Borkovits et al. 2016).

KIC 7955301 is a hierarchical triple system composed of a young RG and a pair of small main-sequence (MS) stars, which was observed for nearly four consecutive years by the original NASA Kepler mission (Borucki et al. 2011) at long (29.4244-min) cadence (Table 1). It was discovered by Gaulme et al. (2013) and Rappaport et al. (2013) who noticed its outstanding eclipse timing, depth, and duration variations. The pair of MS stars shows partial eclipses, but the RG does not eclipse the inner binary. The eclipse timings show a periodicity of about 208.6 days, while the pair of MS stars orbits in 15.3 days. The precession of the orbital plane of the MS stars is quite visible as well, and has a period longer than the time series. This system was part of the cohort of systems with ETVs studied by Borkovits et al. (2016). From eclipse timing only, they estimated the mass of the RG to be 1.5 ± 0.5 M_⊙ and the total mass of the MS pair to be 2.2 ± 0.8 M_⊙. From asteroseismic scaling relations and a quick look at the mixed dipole (l = 1) modes, Gaulme et al. (2013) inferred the RG mass to be 1.2 ± 0.2 M_⊙, the radius to be 5.9 ± 0.2 R_⊙, and the core rotation to be about 30 days. By considering that asteroseismic scaling relations tend to slightly overestimate stellar masses for RGs (Gaulme et al. 2016), the RG is likely a descendent of an F-type star that had a convective envelope during the MS, and the RG is expected to rotate slowly.

With this paper, we firstly aim to test our ability to measure accurate masses in hierarchical triple systems. We know it is possible to determine the mass of a star that belongs to an EB down to 1% (e.g., Maxted et al. 2020). However, no study has specifically been developed to measure the mass of an oscillating star in a hierarchical triple that is a single-lined spectroscopic triple system (ST1). We note that the famous RG in the hierarchical triple system HD 181068 (Derekas et al. 2011), observed by Kepler, is triply eclipsing and does not show any oscillations, likely because of the mode suppression observed in short period multiple systems (Gaulme et al. 2020). Given the clear ETVs and the clear oscillations of the RG, KIC 7955301 appears to be an ideal case for testing our methods on hierarchical triple systems with an oscillating component. Subsequently, the approaches described here should be extended to all known triple systems with ETVs that include an oscillating component (Gaulme et al. 2013).

The second objective is a rather unique opportunity to study an RG star in greater detail, thanks to the multiple approaches used simultaneously to give a clear picture of the system. Our study makes use of the Kepler photometric data and high-resolution optical spectra obtained with the 3.5 m telescope at Apache Point Observatory, and the 1.93-m telescope at the Haute-Provence Observatory (Sect. 2). The spectra are used to determine both the RVs and the atmospheric parameters (Sect. 3). The Kepler light curves are used to analyze both the eclipses for the asteroseismic analysis (Sect. 4) and the dynamical modeling (Sect. 5). Finally, the deduced parameters, such as the oscillation frequencies, atmospheric parameters, and the orbital parameters, are used to optimize a stellar evolution model (Sect. 6).

2. Observations

2.1. The Kepler light curve

We worked with the Kepler public light curves that are available on the Mikulski Archive for Space Telescopes (MAST)¹. Two types of time series are available: the Simple Aperture Photometry (SAP) and the Pre-search Data Conditioning Simple Aperture Photometry (PDC-SAP) light curves. The latter consist of time series that were corrected for discontinuities, systematic errors, and excess flux due to aperture crowding (Twicken et al. 2010). They do not meet our requirements for monitoring possible rotational modulation or eclipse depth, which are often altered during the process (e.g., García et al. 2014; Gaulme et al. 2014). We thus made use of the SAP data to preserve any possible long-term signal. This choice entailed our own detrending and stitching operation on the light curves while ensuring that the rotational modulation was preserved after each interruption of the time series. The methods employed to clean the time series are detailed in Gaulme et al. (2016).

For asteroseismology, light curves with eclipses are an issue because their signal contaminates the Fourier domain, in which we performed the asteroseismic analysis. We tried two different approaches to remove the eclipses from the time series before computing its power density spectrum. The first approach was better in principle: it consisted of modeling each eclipse with a typical eclipse function, for example that from Mandel & Agol (2002) for fitting exoplanetary transits, and then subtracting the model from the light curve. However, this approach, as noted in Gaulme et al. (2016), is less efficient than simply clipping out the eclipses, followed by gap filling with a second-order polynomial. Indeed, an eclipse model is always a little imperfect and when subtracted about 100 times (the number of eclipses during the Kepler run), they still damage the Fourier transform, unlike eclipse removal and filling. Figure 1 (top panel) shows the original time series along with the time series with no eclipses that was used for analyzing the RG oscillations. The bottom panel highlights the eclipses variations in the form of a folded light curve where a shift was introduced between consecutive orbits. The eclipse depth varies from 0.5 to almost 4%, timing varies by about 4 h, and the duration varies from approximately 4.3 to 6.7 h.

Fig. 1. Kepler light curve of KIC 7955301. Top panel: relative flux as a function of time, where t = 0 is the first day of data. The blue curve is the stitched light curve, whereas the orange curve has the eclipses removed and is used for asteroseismic analysis. Bottom panel: stitched light curve folded over the inner binary orbital period (15.32 days), where consecutive eclipses are vertically offset by 0.003 to ease the visibility of the transit depth and timing variations.

2.2. High-resolution optical spectra

From 2012 to 2019, we were granted time on the Astrophysical Research Consortium échelle spectrograph (ARCES) of the 3.5 m ARC telescope at Apache Point observatory (APO), which covers the whole visible domain at an average resolution of 31 000 (Wang et al. 2003). This allowed us to monitor KIC 7955301 23 times, among observations of other RGs in multiple systems. Even though ARCES was not designed for precise RV measurements, it has successfully been used for this purpose in earlier works (e.g., Rawls et al. 2016; Gaulme et al. 2016; Benbakoura et al. 2021). The measurement error reported in these papers is about 0.5 km s⁻¹ for an RG spectrum with a signal-to-noise ratio (S/N) between 10 and 20. In practice, our spectra have a S/N ranging from 10 to 25. The ARCES optical spectra were processed and analyzed in the same way as in Gaulme et al. (2016) and Benbakoura et al. (2021), and we refer the reader to these papers for details.

In addition to the APO spectrograph, we were granted observing time on the SOPHIE échelle spectrograph at the 1.93 m telescope at the Haute-Provence Observatory (OHP). We obtained spectra of KIC 7955301 on June 8 and 9, 2018, and on October 9, 2018. The SOPHIE data processing pipeline directly provides the RVs from the spectra. We refer the reader to Santerne et al. (2011a,b) for details on the data reduction and RV measurements based on SOPHIE spectra.

3. Spectroscopic analysis

3.1. Radial velocities

From the reduced one-dimensional spectra obtained at APO, we computed the RVs with the broadening-function (BF, Rucinski 2002) technique. The fundamental hypothesis of this method is that the observed spectrum is the theoretical spectrum convolved with a broadening function that accounts for stellar rotation and instrumental effects. The BF technique deconvolves the observed spectrum by the theoretical one to extract the BF. For SB2 systems, the BF shows two peaks, one per component.

We employed the theoretical spectra generated by the PHOENIX BT-Settl code (Allard et al. 2003), which were computed with the solar abundances derived by Asplund et al. (2009). As in Gaulme et al. (2016) and Benbakoura et al. (2021), we used templates of MS stars (T_eff = 6200 K) to maximize the chance of detecting the signal from the companion star. Such a choice actually also helped the deconvolution of the RG spectrum because the large number of absorption lines in RG spectra tends to increase the noise in the BF profiles when the original spectra do not have an excellent S/N, which was our case. We computed the BF by making use of the wavelength range 4500−5800 Å. Our final RV data were obtained after correcting the BF profiles from the barycentric corrections, which account for the rotation and revolution of the Earth with respect to the target. We computed the barycentric corrections with the PyAstronomy² routine helcorr, which is based on the Piskunov & Valenti (2002) algorithm.

All the RVs we produced in this work are compiled in Table A.1 and displayed in Fig. 2 together with the series of BF profiles. For illustration purposes we overplotted data points with a fit performed with a simple Keplerian orbit. This allowed us to estimate the standard deviation of the measurements to be about 0.70 km s⁻¹. We note that employing a Keplerian orbit is a very simplistic approach as KIC 7955301 is a triple system. In Sect. 5 we show that the argument of periastron of the outer orbit changed by about 20° during the 4 years of Kepler observations, meaning that a simple Keplerian orbit is not sufficient to accurately model the RG orbit.

Fig. 2. Broadening function of the 23 spectra taken with ARCES. Spectra were sorted by increasing orbital phase. The position along the y axis corresponds with the orbital phase. Barycentric velocity corrections were included in the BF computation. The RV data are indicated with the red plus symbols and their best-fit Keplerian-orbit model by a dashed line.

3.2. Disentangling the spectra

A proper estimate of the atmospheric parameters of the RG and possibly of the inner binary components requires us to disentangle the spectra. For a triple system with such a large flux contrast between the evolved primary component and the inner dwarf binary, it is not an easy task. The spectral lines of the three components are diluted and it is challenging to disentangle the individual contributions to the continuum without knowing the light ratio. It is straightforward to get the ratio from a light curve when all components eclipse each other. Unfortunately, in our case, the RG is not eclipsing with the inner pair. The remaining method for measuring the light ratio and distentangling the spectra consists in scaling the absorption lines from synthetic spectra to the series of observed spectra (Pavlovski et al. 2009, 2018, and references therein).

According to the dynamical model presented in the following sections, the contribution of the MS components to the total luminosity is close to 10%, which implies that we cannot visualize their signatures in either the BF profile or in the spectra. Nevertheless, very faint components have been revealed by spectral disentangling even at a level of only 1−2% of the total light, such as the M dwarf in the EB V530 Ori (Torres et al. 2014), the Roche-lobe filling giant in the inner EB of the Algol triple system (Kolbas et al. 2015), and the MS companions of the RGs in EBs (e.g., Gaulme et al. 2016; Hełminiak et al. 2017; Brogaard et al. 2018; Themeßl et al. 2018; Benbakoura et al. 2021).

In spectral disentangling as originally formulated by Simon & Sturm (1994), the spectra of the individual components are simultaneously reconstructed with an optimization of the orbital elements of the multiple-star system. The problem consists in solving the matrix equation A ⋅ x = b, where the vector x represents the unknown individual spectra of the components – which we aim at extracting – and b represents all of the observed spectra. The design matrix A is constructed from Doppler shifts for a given set of the orbital elements and exposure dates. A could also contain light dilution factors, if known. Since the set of linear equations are overdetermined, the solution may be calculated with the linear-algebra technique known as singular value decomposition. In our particular case, the RVs for the RG component are measured by the BF (Sect. 3.1, and Table A.1). From the dynamical model of the system described in Sect. 5, we computed the expected values of the RVs of the components of the inner binary.

The spectral disentangling was performed with the code CRES (Ilijic 2004), which is based on the Simon & Sturm (1994) method for spectral disentangling in the visible domain. In CRES, the RVs, light dilution factors of each component, and the observed spectra are inputs. The observed spectra should be given in units of wavelength. The observed spectra do not need to be in an uniform scale, and the spectra from different spectrographs do not need to be resampled first. Any portion of the observed spectra could be masked out to exclude undesired wavelength regions.

Altogether, we uses the 23 high-resolution échelle spectra obtained at APO (Sect. 2.2). Spectra have S/Ns from ∼5 to 60, with an average S/N ≃ 27, as measured in several short line-free windows near 5500−5600 Å. The S/Ns were then used to assign weights to the observed spectra. We note that the code actually allows an assignment of the weights to individual pixels, another feature advantageous in the wavelength-domain disentangling, but we did not use this option, assuming it was sufficient to take the weight per spectrum, and not per pixel. The selected spectral segments had various lengths, from about 40 to 90 Å, always taking into account overlapping regions, which serve to double-check the quality of our disentangling.

It is known that without any substantial change in the fractional light of the components in the course of an orbital cycle, there is no unique solution to reconstruct the individual spectra of the components (Pavlovski & Hensberge 2005). One possibility is to optimize the light dilution factors in the calculations of the components’ spectra. However, with the dominant fractional light contribution of the RG component, which is not eclipsing with the inner pair of the MS stars, this would be meaningless. The best option is to perform spectral disentangling in a pure separation mode and then determine fractional light contribution for each component from its spectral characteristics.

The reconstruction of the individual spectra of all three components was performed in the wavelength range 5100−5650 Å. In addition, the spectral range of the disentangled spectrum for the RG component was extended to the range 4700−6825 Å. In the latter case, our intention was to cover the Balmer lines Hα and Hβ, as well as the sodium doublet Na I at λλ 5890 and 5896 Å. Some portions of the reconstructed spectra are shown in Fig. 3. It is evident from the bottom panel of Fig. 3 showing reconstructed spectra for the stars in the inner (eclipsing) pair that the signatures of the stellar spectra are successfully revealed, although the noise level is large. We can estimate the gain in the S/N for the spectra of the disentangled components using simple calculations, since spectral disentangling, in principle, works as co-addition of the observed spectra. By assuming random noise with an average S/N ≃ 27, the 23 spectra should lead to a total S/N of about 130. Of this, the RG spectrum benefits most since it contributes about 91% of the total light, as is determined later in this section. The other two components share the remaining ∼9% – given their contribution is similar (Sect. 5, Table 4) –, meaning that the S/N of their reconstructed spectra is about 5−6. Thus, their spectral signatures are revealed only thanks to the gain of S/N caused by the spectral disentangling. The successful isolation of the spectra for the stars in the inner system, which is based on the RVs calculated from the predicted orbit, is an encouraging confirmation of the correctness of the dynamical model.

Fig. 3. Sample region of the mean spectra of the individual components of the system obtained by disentangling the APO spectra. Top panel: spectrum of the RG component (blue line) and its best fit obtained with the GSSP code (red) from 5340 to 5450 Å. Bottom panel: spectra of the inner binary system Ba and Bb (blue lines) together with their best-fit synthetic spectra (red lines). The optimal values representing the best fits are listed in Table 2.

The deepest absorption lines in the disentangled spectrum of the RG component are due to Mg Ib triplet lines at 5167, 5172, and 5183 Å, and the Na I D doublet at 5890, and 5896 Å. The depths of these sets of spectral lines constrain the fractional light contribution of the RG to the total light of the system. A limit defined by the physical solution is at 89%, which is corroborated with the optimal fitting (Table 2).

It is well established by observational evidence that lithium is depleted in the majority of giants (Brown et al. 1989). There are rare exceptions known as Li-rich giants with abundances A(Li)> 1.5 (a scale where the number of hydrogen atoms log N(H)=12.0). Examination of the disentangled spectrum of the RG component of KIC 7955301 shows almost no trace of the Li I resonance doublet at 6708 Å. An estimate of an upper limit of the lithium abundance was made with a comparison to the synthetic line profiles. The calculations were performed assuming local thermodynamic equilibrium with model atmospheres calculated with Atlas9 (Castelli & Kurucz 2003), while line profiles for Li I 6708 Å were calculated with the spectrum synthesis code UCLSYN (Smith 1992). The atomic data are from Yan et al. (1998). We were able to establish only an upper limit in lithium abundance A(Li)<  − 1.0 dex. This is a common lithium abundance found in modern massive spectroscopic surveys (Luck & Heiter 2007; Buder et al. 2018; Charbonnel et al. 2020).

3.3. Stellar atmospheric parameters

For the analysis of the disentangled spectra of all three components of the KIC 7955301 system, we employed the Grid Search in Stellar Parameters³ (GSSPTkachenko 2015) software package, specifically its GSSP_SINGLE module. It is a grid search-based spectrum analysis algorithm that (on the fly) generates synthetic spectra in an arbitrary wavelength range based on a pre-computed grid of model atmospheres, and performs a comparison between the observed spectrum and each synthetic spectrum from the grid in the χ² statistical framework. The reported 1-σ uncertainties are computed from the χ²-statistics and take into account possible correlations between the free parameters. The GSSP algorithm allows for the simultaneous optimization of the effective temperature, the surface gravity, the micro- and macro-turbulent velocities, the projected rotational velocity, and the metallicity of the star. Optionally, one can also optimize for the degree of the light dilution in the spectrum due to the star being a member of a binary and/or higher-order multiple system.

In this study, we employed the grid of LLMODELS model atmospheres (Shulyak et al. 2004) and opted for the inclusion of the light dilution factor into the optimization to account for the a priori unknown dilution of the disentangled spectrum of each of the three stellar components under analysis. The light dilution factor was assumed to be wavelength independent, which is a fair assumption for the RG component given it is the dominant contributor by far (>90%) to the composite spectrum of the system, as well as for the two MS components, given the limited wavelength interval of some 500 Å that could be used for the disentangling. We additionally note that the macro-turbulent velocity parameter was ignored in the analysis of all three stellar components of the system, and we also fixed log g for both MS components to the values inferred from the light curve solution, we fixed the micro-turbulent velocity to 2 km s⁻¹, and we assumed solar chemical composition ([M/H] = 0.0 dex) for both of them. The choice to fix so many parameters in the spectrum analysis of the MS components was dictated by their small cumulative contribution (<8−10%) to the total light of the system, and hence their very noisy disentangled spectra. The consequence of fixing the macro-turbulent velocity parameter for the RG component was that the v sini value reported was in fact representative of the combined spectral line broadening due to the effects of rotation and macroturbulence. Finally, we note that we explored three different options in the analysis of the RG spectrum, namely fixing the surface gravity to the two values inferred from the light curve solution, reported in Table 2, and treating log g as a free parameter.

4. Asteroseismic analysis

4.1. Global oscillation properties

The stellar granulation and accurate values of ν_max and mode amplitude H_max were estimated by fitting the power spectrum (Fig. 4) as commonly performed in asteroseismology (Kallinger et al. 2014), and already used in Gaulme et al. (2016). Following Kallinger et al. (2014), the power density spectrum is fitted by

(1)

Fig. 4. Oscillation spectrum of KIC 7955301. Top panel: power spectral density of the time series after eclipse removal and filling. The dashed red line represents the stellar background noise: it is the sum of the dashed blue lines (correlated noise Harvey profiles) and the dashed green line (white noise). The solid red line represents the fit of the power spectrum, including the stellar background and the Gaussian envelope of the oscillations. Bottom panel: échelle diagram associated with the frequency spacing Δν = 10.48 μHz. The colored lines indicate the oscillation universal pattern, where blue is l = 0, green l = 1, and red l = 2. Modes of degree l = 3 are visible half way in between the l = 0 and l = 1 ridges. The horizontal dashed line indicates the location of νmax.

where N is the function describing the noise, η is a damping factor originating from the data sampling, B is the sum of three “Harvey” functions (super Lorentzian functions centered on 0), and G is the Gaussian function that accounts for the oscillation excess power:

(2)

The terms ν_max and H_max are the central frequency and height of the Gaussian function. The best-fit values are: ν_max = 124.89 ± 0.33 μHz, H = 354.2 ± 10.5 ppm μHz⁻¹, and σ = 12.5 ± 0.4 μHz.

A first estimate of Δν was performed with the envelope of the autocorrelation function (EACF) developed by Mosser & Appourchaux (2009) from the whitened power spectral density (power spectrum divided by background function). From the EACF, Δν = 10.46 ± 0.05. Then, we used the universal pattern of RGs introduced by Mosser et al. (2011) to correct it. The principle of this method is to compare the measured oscillation frequencies to a theoretical law, the so-called universal pattern, predicting the variations in these frequencies as a function of Δν and the radial order. That way, the value of Δν was revised to 10.49 ± 0.02 μHz.

We then computed a proxy of the stellar masses and radii using the asteroseismic scaling relations that were originally proposed by Kjeldsen & Bedding (1995) for solar-like MS oscillators, and then successfully applied to RGs (e.g., Mosser et al. 2013). We employed the asteroseismic scaling relations as proposed by Mosser et al. (2013) for RGs; in other words, where ν_max, ⊙ = 3104 μHz, Δν_⊙ = 138.8 μHz, T_eff, ⊙ = 5777 K, and where the observed Δν is converted into an asymptotic one, such as Δν_as = 1.038 Δν. By considering the temperature T_{eff, A, sce.A} = 4720 ± 105 K that was found from the disentangled spectrum by fixing the surface gravity at the value obtained from the most complete model log g = 3.0 (scenario A, Table 2), the stellar parameters are R_A = 5.91 ± 0.07 R_⊙, and M_A = 1.27 ± 0.05 M_⊙.

Finally, we checked the amplitude of the oscillations to see whether tidal interactions in the system had altered their properties, as was observed in close systems by Gaulme et al. (2014, 2020), and Benbakoura et al. (2021). From the background fitting, the height of the Gaussian function is 354 ± 11 ppm μHz⁻¹, which is a typical value for an RG with ν_max = 125 μHz, according to the sample of 4500 RGs analyzed by Gaulme et al. (2020, Fig. 7) with the exact same codes. We note that the height of the oscillation envelope is a little underestimated because the RG only contributes approximately 91% of the photometric flux, whereas most of the 4500 stars displayed in Gaulme et al. (2020, Fig. 7) are not in multiple systems and do not suffer from this dilution factor.

An alternative metric on the amplitude of the oscillation modes consists in measuring the amplitude of the largest l = 0 mode, as was performed by Benbakoura et al. (2021). The fitting of the l = 0 peaks with Lorentzian functions leads to a maximum l = 0 amplitude of 7.06 ± 0.08 ppm, which is in agreement with the known RGs in EBs that do not show any alteration of their oscillation properties (Fig. 5). We conclude that the oscillations of KIC 7955301 have regular amplitudes and do not show any global suppression, as was observed for RGs in short period binary systems. This result is consistent with the long-period binary systems observed so far (Benbakoura et al. 2021).

Fig. 5. Oscillation amplitude. Top panel: height of the Gaussian envelope used to model the contribution of the oscillations to the stellar background. The red dot highlights the position of KIC 7955301 with respect to a sample of 4500 RGs analyzed by Gaulme et al. (2020) indicated in blue. The sample splits into three categories: the oscillators with regular amplitude, those with depleted l = 1 modes, and the active RGs that show a global mode suppression (l = 0, 1, 2, 3). KIC 7955301 falls in the regular oscillators. Bottom panel: Amplitude of the largest l = 0 mode of KIC 7955301 compared with those of the 35 RGs in EBs that were studied by Benbakoura et al. (2021). The panel displays the relative difference between expected and measured oscillation amplitudes (%) as a function of orbital period Porb (days). The red line indicates 0, i.e., stars whose oscillations are not detected. The two dashed blue lines represent the region in which relative mode amplitude lies for systems with orbital periods longer than 180 days within two sigma. The size of each symbol represents the amplitude of stellar variability (large means variable, small means not variable), and the gray scale indicates the pulsation mode amplitude (white means no modes; black means a large amplitude).

4.2. Mixed modes

4.2.1. Formalism

As indicated in the Introduction, the RG component of KIC 7955301 shows a rich spectrum of dipolar mixed pressure and gravity modes. Here we briefly summarize some basic properties of the mixed modes (see Mosser et al. 2015, and references therein). The frequency of the mixed modes is an implicit expression given by

(3)

where q is the coupling factor between the p and g modes, θ_p is the phase of the mixed modes with respect to the asymptotic p-mode frequencies separated by the large separation Δν, while θ_g is the phase of the mixed modes with respect to the asymptotic g-mode periods separated by the period spacing ΔΠ₁. In other words,

(4)

where ν is the mixed mode frequency and ϵ_p is a phase offset, and

(5)

where P is the mixed mode period and ϵ_g is a gravity offset. Solving the implicit relation given by Eq. (3) provides the mixed mode frequencies.

The period separation ΔP between two consecutive mixed mode periods is given by

(6)

where ζ is related to the phases given above by

(7)

ζ is also useful for expressing the frequency splitting of these dipole modes δν_rot, since we have

(8)

where δν_core is the rotation mainly sampled by the g modes, while δν_env is the rotation mainly sampled by the p modes. We stress that this expression is not symmetrical in terms of the azimuthal order m since it depends on the mixed mixed frequency

(9)

as discussed in Mosser et al. (2012).

Equations (3) and (8) provide a full description of the mode frequency and its associated rotational splitting. To ensure a robust analysis of the mixed modes, four independent studies were led by co-authors Appouchaux, Gehan, Mosser, and Vrard. All followed the approach developed by Mosser et al. (2015), where the first step consists of stretching the oscillation spectrum to transform the mixed-mode pattern into a comb-like pattern based on the gravity period spacing ΔΠ₁. Despite a common approach, each fitter performed a different analysis with independent routines. The following subsections provide details for each analysis.

4.2.2. Analysis adapted from Appourchaux (2020)

The stretched power spectrum was interpolated onto a regular grid for producing the échelle diagram of Fig. 6 (top panel). The fitted asymptotic frequencies of the multiplet of the dipole mixed modes were computed using the expressions of Mosser et al. (2015), where ζ was replaced by the analytical formulation given above. We note again that the mode splitting was not symmetrical due to ζ depending on frequency. First guesses for ΔΠ₁, rotation, and q were provided by the results of Mosser et al. (2014, 2015, 2017), respectively. These three parameters were then visually adjusted to obtain a figure similar to that of Fig. 6. The asymptotic frequencies obtained by this procedure were very close to the actual peaks, within 0.1 μHz, easing the fitting of the peaks by a maximum likelihood estimation (MLE). Due to the dense forest of dipole mixed modes, the fitting window for these modes was 0.5 μHz, while for l = 0, 2 and 3, it was 3 μHz. Fitted frequencies are displayed in Table B.1. The échelle diagram of the power spectrum with the fitted frequencies is shown in Fig. 7.

Fig. 6. Stretched period échelle diagram for dipole gravity-dominated mixed modes. The top panel shows the work led by co-author Appourchaux with ΔΠ1 = 76.3 s. Fitted ridges with m = {− 1, 1} are represented as red and orange crosses, respectively, while m = 0 is represented as a continuous white line. The black stripes show the location of the l = 0, 2 mode power, which has been removed for clarity. The bottom panel summarizes the findings of co-author Gehan, where peaks with a height-to-background ratio equal to or above ten are represented. The symbol size varies with the measured power spectral density. Ridges with m = {− 1, 0, 1} are represented in green, light blue, and red, respectively. Red and green crosses represent the fit of the ridges with m = ±1. The m = 0 ridge is identified by considering that it is located at the mid-point between the m = ±1 ridges.

Fig. 7. Peak-bagging results, including mixed-mode analysis by co-author Appourchaux, displayed in the form of an échelle diagram of the power spectrum for Δν = 10.44 μHz. Fitted frequencies are shown: l = 0 by white diamonds, l = 2 by orange diamonds, l = 3 by red diamonds, l = 1, m = −1 by white crosses, l = 1, m = 0 by white circles, and l = 1, m = +1 by white pluses.

In a second step, we determined all the parameters of the asymptotic expression of the dipole mode frequencies, namely ΔΠ₁, q, ϵ_g, δ₀₁, ν_core, and ν_env by carrying out an unweighted least squares fit of the fitted frequencies. The optimization was performed by using the same algorithm as in Appourchaux (2020). Data were trimmed to exclude three frequencies out of the 62 dipole mode frequencies (see frequencies in Table B.1). The root mean square value of the optimized difference was 0.00 ± 0.03 μHz. The error bars were then extracted by performing a Monte Carlo simulation using the error bars on the dipole mode frequencies extracted from the MLE. The optimization was then repeated 100 times. The values obtained are reported in Table 3.

The inclination angle of the star was extracted using a single mixed mode (n = −94) for which the triplet was clearly fitted. We used the ratio between heights of modes with different azimuthal orders (Gizon & Solanki 2003). The ratio between the m = 0 and the m = ±1 was 0.08, providing an angle of degrees.

4.2.3. Analysis following Gehan et al. (2018)

Before the mixed-mode analysis, the frequency at maximum amplitude ν_max was derived by using the FRA pipeline. It computes a smoothed power spectrum, then fits a Gaussian envelope accounting for oscillations along with a background contribution to measure ν_max. It also performs a likelihood ratio test to validate the measured ν_max and check whether oscillations are present. The FRA pipeline was used in Huber et al. (2022) and is explained in detail in Gehan et al. (2022). The large spacing Δν was measured by using the EACF.

The power spectrum was then stretched as in Mosser et al. (2015) and Gehan et al. (2018). Pressure-dominated mixed modes were located using the RG universal oscillation pattern (Mosser et al. 2011) and excluded, in order to keep only gravity-dominated modes that are mostly sensitive to the core rotation rate. The stretched period échelle diagram reveals that dipole gravity-dominated mixed modes line up along two ridges, associated with the azimuthal orders m = ±1 (Fig. 6, bottom). Fitting these ridges (Fig. 6) allowed us to derive the mean core rotational splitting, which was found to be δν_rot, core = 380 ± 6 nHz (Gehan et al. 2018, for details). Once the azimuthal order of mixed modes was identified, we measured the inclination angle of the rotation axis from the ratio between heights of modes with different azimuthal orders (Gizon & Solanki 2003). Only the ridges with azimuthal orders m = ±1 were visible, but we knew where to look for the missing m = 0 ridge in the background, at the midpoint between the m = +1 and m = −1 ridges (Fig. 6, bottom). We obtained an inclination of degrees. We refer the reader to Gehan et al. (2021) for the details of the inclination measurement procedure. We note that the maximal possible value according to the uncertainties is i = 90°, so that the star is possibly seen equator-on. This is because the m = 0 ridge is lost in the background, and thus there is a possibility that the peaks with m = 0 are completely absent and that the observed signal is only due to background noise.

4.2.4. Analysis following Mosser et al. (2018)

Stretching the power spectrum first requires a precise fitting of the radial and quadrupole modes l = 0 and 2, thanks to the RG universal oscillation pattern, in order to precisely locate the expected pure pressure dipole modes. Firstly, the modes with the largest height-to-background ratios (HBR) were used to determine the parameters of the fit. Secondly, we picked modes with lower HBRs, still larger than five, provided that they stuck to the asymptotic pattern within a frequency range wide six times the spectral resolution. We applied the method of Mosser et al. (2018) for analyzing the rotational splittings to overcome the difficulties arising from the overlap of rotational multiplets larger than period spacings. The global seismic parameters derived from the fit of the mixed-mode pattern, following the updated approach described by Pinçon et al. (2019), are reported in Table 3. The frequency analysis is displayed in Fig. 8.

Fig. 8. Peak-bagging results, including mixed-mode analysis by co-author Mosser, in the form of a line-by-line échelle diagram, for radial pressure orders eight to 13. Radial and quadrupole modes are highlighted in red and green, respectively. The expected locations of the dipole mixed modes are labeled with their mixed radial orders, and are indicated with a color depending on the azimuthal order: dark blue (m = −1), light blue (m = 0), or purple (m = +1). ℓ = 3 modes, not identified in this plot, are located near the abscissa 0.21. The gray dashed lines correspond to height-to-background ratios of seven and ten.

4.2.5. Analysis adapted from Vrard et al. (2016)

The steps that follow the power-spectrum stretching differ from the analysis led by co-author Mosser. We started by computing the Fourier transform of the stretched spectrum to obtain a first ΔΠ₁ value following the work of Vrard et al. (2016). After that, the coupling parameter q was determined following the method described in Mosser et al. (2017), which involved searching for the maximum of the Fourier Transform of the stretched oscillation spectrum. Finally, an iterative process was performed between ΔΠ₁ and q: the ΔΠ₁ was re-measured with the measured q value, then the inverse was performed until the ΔΠ₁ and q values converged. The final measurement of ΔΠ₁ and q was obtained when those values converged. The final fitted parameters led to ΔΠ₁ = 76.09 ± 0.72 for the gravity mode period spacing and q = 0.17 ± 0.05 for the coupling parameter.

4.2.6. The big picture from the mixed-mode analysis

The four independent approaches used to analyze the mixed dipole modes agree (Table 3). Firstly, the measured mixed-mode period spacing ΔΠ₁ ≈ 76 s associated with the large frequency p-mode spacing Δν ≈ 10.5 μHz definitely classifies the star as a red giant branch (RGB) star (Mosser et al. 2014). Secondly, the fact that the m = 0 modes are not detected indicates a large inclination angle i of the RG rotation axis with respect to the line of sight. The two authors who carefully analyzed the ratio of the amplitudes of the mixed m = 1 to m = 0 dipole modes report a very consistent inclination angle of ≈75°, with a rather large uncertainty ranging from 60° to 90°. Thirdly, according to Mosser et al. (2012, Eqs. (22) and (23)), the frequency spacing measured on the mixed modes provides a proxy of the mean core rotation period with the relation P_core ≈ 1/2δν_core = 15.1 ± 0.2 days by assuming an uncertainty of 5 nHz on δν_core. The average splitting of the envelope determined by co-author Appourchaux is about ten times smaller, leading to a mean rotation period of days. If we add the fact that no spot-related variability was measured from the Kepler light curve and that the oscillation amplitude meets our expectations for that type of RG, we can safely conclude that the RG is a regular slow rotator that does not exhibit any sign of tidally increased spin.

5. Dynamical model

5.1. Method

The combined photodynamical analysis of the Kepler photometry, the ETVs of the inner binary, and the RV curve of the outer giant component offer a robust dynamical determination of the mass of the RG (M_A) and of the total mass of the inner binary (M_bin). The amplitude of the RV curve gives the spectroscopic mass function:

(10)

where q_out = M_bin/M_A and i_out are the mass ratio and the inclination of the outer (or wide) binary formed by the RG and the center of mass of the inner eclipsing pair. In addition, Borkovits et al. (2016) showed that the ETV of the EB is strongly dominated by the gravitational three-body perturbations against the pure light-travel time effect. Detailed analytical investigations of third-body perturbed ETVs (Borkovits et al. 2011, 2016) have shown that such curves tightly constrain almost all of the orbital elements of both the inner and outer orbits, together with the mutual inclination of the two orbital planes. Their amplitude is related to the outer mass ratio as

(11)

Finally, the variation in the eclipse depth and the duration of the inner EB offer a very precise determination of the varying inclination of both the inner orbit of the eclipsing pair and the outer orbit.

The combined photodynamical analysis was carried out with the software package LIGHTCURVEFACTORY (see, e.g. Borkovits et al. 2019, and further references therein). This code contains a built-in numerical integrator to calculate the three-body perturbed coordinates and velocities of the three bodies, a multiband light curve, ETV and RV curve emulators, and a Markov chain Monte Carlo (MCMC)-based parameter search routine for the inverse problem. The latter employs an implementation of the generic Metropolis-Hastings algorithm (see e.g. Ford 2005).

5.2. Best-fit model

We simultaneously analyzed the eclipse photometry from the Kepler light curve, the primary and secondary ETVs deduced from the same light curve, and the RVs of the outer RG component. For the Kepler light curve, we only used the phase sections of each eclipse to reduce computational costs. We now describe the set of parameters that were adjusted with the MCMC, once preliminary tests were over.

The set includes nine of the twelve orbital parameters: regarding the inner orbit of the eclipsing pair, we fit e₁cosω₁, e₁sinω₁, and i₁, where e₁ is the osculating eccentricity, ω₁ the argument of periastron, and i₁ the orbital inclination. About the wide outer orbit, besides e₂cosω₂, e₂sinω₂, and i₂, we fit the anomalistic period P₂, the periastron passage time τ₂ and the longitude of the node Ω₂. Regarding two of the three “missing” parameters describing the inner orbit, namely, the period P₁ and time of a primary eclipse or, in other words, the time of an inferior conjunction of the secondary component of the inner pair , they were constrained internally at each trial step via the ETV curves in the manner described in Borkovits et al. (2019). Finally, for the third “missing” orbital parameter, the longitude of the node of the inner orbit Ω₁ was set to zero. This could be done because, for the complete analysis, only the difference of the nodes, that is, ΔΩ = Ω₂ − Ω₁, has relevance. Therefore, in such a manner, by setting Ω₂, one sets ΔΩ. We note that, because the motion of the three stars are not purely Keplerian, none of the orbital parameters are constant. Therefore, all the trial values either adjusted, constrained, or fixed in the initializing phase of any trial steps refer only to a given epoch t₀.

We also fitted three parameters connected to the masses of the three stars: the RG mass M_A, the spectroscopic mass function f(M_B), and the mass ratio of the EB q₁ = M_Bb/M_Ba.

Finally, the fit includes four parameters almost exclusively connected to the light curve: the first is the duration of the primary eclipse around the epoch t₀ (Δt_pri), which is related to the sum of the fractional radii of the members of the inner EB, that is, the radius divided by the semimajor axis. The second is the ratio of the stellar radii of the inner pair R_Bb/R_Ba. The third is the RG radius R_A. The fourth is the ratio of the effective temperatures of the eclipsing pair T_Bb/T_Ba. We note that in case of a single-band photometric observations, the EB light curve only carries information about the temperature ratio – most strictly speaking, the passband-dependent surface brightness ratio – of the two components, rather than their true temperatures. Therefore, one of the temperatures should be taken from some external sources. In a single binary, it can be obtained from spectroscopic analysis or from stellar energy distribution (SED) data (see, e.g., Miller et al. 2020). In the present situation, most of the total flux comes from the third RG component, which makes the classical approach irrelevant. Therefore, we assumed that the two components of the inner binary are still on the MS, by being less massive than the RG, so that at each trial step the software sets the effective temperature of the primary component of the eclipsing pair T_Ba accordingly to the mass – radius and mass – effective temperature relations of Tout et al. (1996) for MS stars. About the RG component, we assumed that its contribution to the light curve is only a constant extra flux⁴. Therefore, we fixed its temperature to T_A = 4800 K for all runs. In such a way, adjusting the RG radius R_A directly sets its bolometric luminosity, and indirectly its total flux in the Kepler photometric band. Hence, the adjustment of R_A was simply used to set the amount of the extra flux to the EB light curve, or in other terms, it served as a substitute for the usual third light parameter.

Regarding other smaller effects that are mainly light curve related, we used a logarithmic limb-darkening law. The corresponding parameters were calculated internally at each trial step with the use of the publicly available passband-dependent tables of the PHOEBE software (Prša & Zwitter 2005)⁵. For the well detached nature of both the inner and the outer binaries and, therefore, for the practically spherical shape of the stars, the gravity darkening coefficients have no detectable effect on the light curve solution. Therefore, we did not consider the results of the recent, more sophisticated study of Claret & Bloemen (2011), but simply adopted the fixed value of g = 0.32, which is based on the seminal model of Lucy (1967). Finally, we neglected both the reradiation and illumination, and the Doppler-boosting effects (Loeb & Gaudi 2003; van Kerkwijk et al. 2010). While the first effect did not really have any influence on the light curve of our system, this is certainly not true for the second one. In this regard, we note that, as one can see in the left panel of Fig. 9, Doppler-boosting would result in an ∼1000 ppm amplitude variation in the out-of-eclipse flux level with the period of the outer orbit, which, in theory, should be detected with Kepler. Despite this fact, it cannot be seen in the light curves, because this low amplitude variation, having a period longer than two quarters, was filtered out during the data processing, which justifies its negligence in the analysis.

Fig. 9. Long-cadence Kepler light curve of KIC 7955301 together with the joint photodynamical model light curve solution. Upper left panel: complete Q0–Q17 light curve (blue) with the best-fit solution (red). This solution was obtained by neglecting the Doppler-boosting effect. The gray curve shows the same model after “switching on” the Doppler-boosting. While this effect, in theory, should have been observable with the accuracy of Kepler photometry, the reasons for neglecting it are discussed in the text. Upper right panel: one-month-long section of the Kepler observations overplotted with the photodynamical solution. The pale blue circles represent each individual observation but, for the analysis, only dark blue data (located within the phase environment of each eclipses) were considered. Lower panels show the residual data.

The median values of the orbital and physical parameters of the triple system derived from the MCMC posteriors and their 1-σ uncertainties are tabulated in Table 4. The observed versus model photometric, RV, and ETV curves are plotted in Figs. 9–11.

Fig. 10. Radial-velocity data with their best fit models. Left panel: radial-velocity data (red) as a function of time expressed in modified Julian dates, with their best fit model (blue line). Right: same folded over the RG orbital period and displayed as a function of orbital phase. The dispersion of the best-fit curve in blue is caused by the varying shape of the RV curve as a function of time.

Fig. 11. Eclipse timing variations of KIC 7955301 on different timescales: observations vs model and predictions. Left: primary and secondary ETV curves derived from Kepler observations (red circles and blue boxes, respectively), together with the photodynamical model solution (smaller, pale red and blue symbols, connected with straight lines, respectively). Right: eclipse-timing-variation-like curves derived from the times of the inferior and superior conjunctions of the secondary of the inner pair (pale red and blue lines), overplotted with forecasted primary and secondary ETV points (larger red and blue symbols) during the eclipsing phases of the inner binary for the first half of the present century. The dynamically forced short-period apsidal motion of the inner binary is nicely visible. The gray-shaded area stands for the interval of the Kepler measurements.

6. Stellar evolution model

We modeled the system in four independent manners by using parameters that were deduced from the observations. The first two manners modeled the RG component only, by considering the seismic information and atmospheric parameters. The first method (PARAM) made use of the asteroseismic global parameters ν_max, Δν, and ΔΠ₁ to optimize the model from a grid of stellar evolution models. The second method made use of the individual frequencies of the radial modes and ΔΠ₁. The third method determined the age of the system by optimizing isochrones and evolutionary tracks on the classical parameters, without inputs from asteroseismology. The fourth method was a complete model of the triple system, which included both the RG and the inner binary, but that did not make use of the asteroseismic constraints.

6.1. Red-giant model based on the global seismic parameters

We used the code PARAM (Rodrigues et al. 2017) to infer the radius, mass, and age of the RG component by using a combination of seismic and non-seismic constraints. On the one hand, the average large frequency separation was computed using the radial-mode frequencies of the models in the grid, not added as an a posteriori correction to the scaling relation between Δν and the square root of the stellar mean density. On the other hand, ν_max in the model grid was computed using a simple scaling relation (Kjeldsen & Bedding 1995), by considering ν_max, ⊙ = 3090 μHz. The period spacing could also be included in the list of seismic constraints, in which case it was compared with the asymptotic value of ΔΠ₁ coming out of the models (see Rodrigues et al. 2017). The grid of stellar evolution models used in PARAM is the same as the reference grid adopted in Khan et al. (2019) and Miglio et al. (2021) (G2).

We used the individual radial-mode frequencies in Table B.1 to compute an average Δν that could be directly compared with that from the model grid, and obtained a value (Δν = 10.48) compatible with that obtained using the EACF (see Sect. 4). With the observational constraints Δν, ν_max, T_eff, and [M/H], we obtained R_A = 5.84 ± 0.09 R_⊙, M_A = 1.27 ± 0.05 M_⊙, and an age of 4.9 ± 0.9 Gyr. In Fig. 12 we show the posterior probability density function (PDF) of radius, mass, and age, compared with the same results from the asteroseismic modeling based on individual frequencies that is exposed in the next section (Sect. 6.2).

Fig. 12. Mass, radius, and age posterior probability distribution functions derived from seismic constraints. Empty histograms show the distributions obtained using PARAM, and the dotted lines show their associated ±1σ credible regions. Solid histograms and shaded regions indicate the same quantities obtained using the individual mode frequencies.

6.2. Red-giant model based on the oscillation frequencies

We generated models using the MESA code (Paxton et al. 2011, 2013, 2015; Paxton et al. 2018, 2019) with diffusion and settling of heavy elements, with no core overshoot, no mass loss, using an Eddington-gray atmosphere, and with chemical abundances evaluated relative to the solar values of Grevesse & Sauval (1998).

The nonradial mixed modes in RGs can be described as combinations of pure p- and g-modes – or “π” and “γ-modes”, in the sense of Aizenman et al. (1977). We evaluated these pure π- and γ-mode frequencies using the modified wave operator construction described in Ong & Basu (2020), which splits the standard wave operator into two complementary p- and g-mode wave operators, each supporting only one type of wave propagation. We solved for these mode frequencies numerically with the pulsation code GYRE (Townsend & Teitler 2013). The observed p-dominated mixed modes of degree l ≥ 2 in RGs are very well approximated as pure π-modes, and as such we matched the model π-modes against the observed modes of degree l ≥ 2. Likewise, dipole modes far from the p-dominated mixed modes have period spacings close to those of the pure γ-modes; accordingly, rather than evaluating ΔΠ₁ using its asymptotic value, we computed it directly from the pure dipole γ-modes.

The mode frequencies returned from stellar models are known to be systematically offset from those that would be obtained from stars of identical interior structure, owing to systematic errors in the modeling of stellar surfaces – this is referred to as the asteroseismic “surface term”. As such, the model frequencies cannot be used directly to construct likelihood functions against the observed mode frequencies in the usual fashion. Instead, we used a nonparametric description – “ϵ-matching”, via the algorithm of Roxburgh (2016) – to characterize the discrepancy between the model and observed mode frequencies in a surface-independent fashion. For this purpose we used only the model l = 0, 2, 3 frequencies, since those modes that we observed were either radial p-modes, or could be approximated as being entirely p-dominated; here we used the mode frequencies reported by co-author Appourchaux. For the dipole modes, which are more strongly mixed, we used the period spacing only, also as reported by co-author Appourchaux. Since the g-mode cavity is well-localized into the stellar interior, no surface correction was needed.

To supplement these seismic constraints, we used the spectroscopic constraints from scenario C ([M/H] and T_eff). These were used to compute log-likelihood functions as

(12)

where f was initially set to one when performing the optimization. Optimization was performed using the differential-evolution algorithm as implemented in the yabox package (Mier 2017), using the initial mass M, initial helium and metal mass fractions Y_i and Z_i, and the mixing-length efficiency parameter α_MLT as the independent variables for this optimization. After the optimization was complete, f was set to the minimum value of along the optimization trajectory, in order to ensure that the reported results were not entirely dominated by the seismic quantities, which were significantly more precisely constrained than the spectroscopic quantities.

We may consider the optimization trajectory to yield a series of nonuniformly distributed samples of the likelihood function over the input parameter space. As such, if we were to report estimates of stellar properties by naively taking likelihood-weighted averages of the desired quantities over all the sampled points, or find posterior distributions by constructing likelihood-weighted histograms, we would be over-representing parts of the parameter space where the density of samples is high to begin with. Effectively, the sampling function of the optimization trajectory provides an implicit prior distribution, which we estimated using a Gaussian kernel density estimator and divided out to yield a nonuniformly sampled posterior distribution function assuming a uniform prior. Using this uniform-prior posterior distribution, we report the posterior median and ±1σ quantiles of the following of fundamental stellar properties:

(13)

(14)

(15)

(16)

(17)

(18)

(19)

We compare these results with those obtained with the PARAM code, by plotting the posterior distributions of the mass, radius, and age in Fig. 12. The mass and radius appear to be in agreement with the values obtained from only Δν, ν_max, and ΔΠ₁ with the PARAM routine, while the ages are somewhat different, although still within the 1σ error. We discuss this difference in Sect. 7.2.

6.3. Overall system model based on classical parameters

We have compared the classical parameters of the three components of the system with the isochrones and tracks from the YaPSI database (Spada et al. 2017); the results are summarized in Fig. 13. In the two panels on the left of the figure, the measured masses and radii are compared with several isochrones of appropriate metallicity, generated by interpolation from the YaPSI database. The best fit was obtained for an age of 4.5 ± 0.2 Gyr for the primary component (top left panel); the components of the inner binary are marginally compatible with this isochrone (bottom left). In the right panel of Fig. 13, we show the best-fitting isochrone, together with evolutionary tracks constructed with the same code and input physics used for the YaPSI models (see Spada et al. 2017 for details). For the primary star, both the isochrone and the track are in good agreement with the observed radius and luminosity (we note that the isochrone and track for component A essentially coincide on the RGB, as a result of the fast-paced evolution typical of this evolutionary phase). The agreement is less good for the two components of the inner binary. This could be ascribed to nonstandard effects not taken into account in our models.

Fig. 13. Visual fit of observed classical parameters of the system with theoretical tracks and isochrones from the YaPSI database (Spada et al. 2017). Left: comparison in the mass radius plane with isochrones of different ages for the primary (top panel) and the two components of the inner binary (bottom). Right: comparison in the radius-luminosity plane; a 4.5-Gyr isochrone and evolutionary tracks of appropriate mass and metallicity are shown.

6.4. Model based on isochrones and SED

We carried out a second spectro-photodynamical study with the software package LIGHTCURVEFACTORY that includes the light curve, ETVs and RVs, as earlier, but also the SED, a stellar evolution code (PARSEC, Bressan et al. 2012), and the Gaia-based accurate trigonometric distance. The cumulative SED, which consists of pass-band magnitudes, drives the temperature of sources, or at least of the dominant one. That way, effective temperatures together with stellar masses and radii determine the locations of the stars on co-eval PARSEC evolutionary tracks. Besides the light curve, RVs, and ETVs, the publicly available multipassband magnitudes (see Table 1) were simultaneously fitted against the theoretical passband magnitudes interpolated from theoretical PARSEC isochrone tables. We refer the reader to Borkovits et al. (2020) for a detailed description of the method.

With this new model, the free parameters partially differ from Sect. 5. Firstly, three new free parameters are added: the logarithm of the age, the metallicity [M/H] of the system, and the interstellar extinction E(B − V). Furthermore, the distance of the system is constrained a posteriori at each trial step, by minimizing the value of . Secondly, the stellar radii and temperatures that used to be free are no longer fitted. However, the code computes them internally by interpolating the appropriate PARSEC tables according to the trial values of the [mass, age, metallicity] triplets at each trial run. The results are tabulated in the last columns of Table 4.

The comparison of the results from the purely dynamical analysis, which is independent from stellar evolution models, and the analysis including dynamical analysis, SED, and isochrones, produce output parameters that lie within their 1-σ uncertainties. This means that the inclusion of stellar evolution models in the complex analysis did not lead to any undesired bias in the values of the physical properties of the system. It also means that our model-independent light curve solution has led to dimensionless quantities, such as fractional radii and temperature ratios, which are fully consistent with stellar evolutionary models.

The mass of the RG component is in excellent agreement with the asteroseismic analysis. The inner binary components appear to be on the MS with masses , M_Bb = 0.88 ± 0.01 M_⊙, and radii of R_Ba = 0.90 ± 0.01 R_⊙, R_Bb = 0.81 ± 0.02 R_⊙ in both solutions.

In contrast, the radius of the RG component R_A and the temperatures T_Ba, T_Bb of the inner stars differ by about 5 − 10% between the two models. This is expected because these quantities are not functions of the light curve properties and of the dynamics of the system. Only the ratio of the effective temperatures of the inner binary T_Bb/T_Ba is model independent; this is in good agreement (within 1%) between the two models. In the purely dynamical model (Sect. 5), the RG temperature was fixed to the available catalog value T_A = 4800 K. Here, the combined dynamical, SED, and isochrone model confirms this value by finding K, which is within the uncertainty of the value we measured from the disentangled optical spectra (Table 2). In the dynamical model, the temperature T_Ba of the primary star of the EB was internally set with the use of the zero age MS mass-radius and mass-luminosity relations of Tout et al. (1996), and the temperature ratio T_Bb/T_Ba was an adjustable parameter. In the combined model they were constrained from the appropriate PARSEC isochrone. This latter model has resulted in higher temperatures by ∼200 K in accordance with the fact that KIC 7955301 is quite an old system, having an age of 4.5 ± 0.2 Gyr according to the present model. Finally, the combined model predicts the RG radius to be , which is in agreement with the asteroseismic value.

The stellar metallicity [M/H] = and interstellar extinction E(B − V)=0.09 ± 0.02 obtained from the combined model are consistent with catalog values (see Table 1). Finally, the SED fitting provides an estimate of the photometric distance of pc, which is in perfect agreement with trigonometric distance of d = 1375 ± 35 pc derived from the Gaia EDR3 catalog by Bailer-Jones et al. (2021). However, we keep in mind that Gaia EDR3 data are not corrected for binarity and might suffer from systematic errors for such a triple system.

Regarding the dynamical evolution of the system, the small outer-to-inner period ratio P₂/P₁ ≈ 13.8 puts the system among the category of compact hierarchical triple systems. That being said, given the small mutual inclination of and almost circular inner orbit (e₁ = 0.0276 ± 0.0001), the orbital configuration of the system is stable on a dynamical timescale. Consequently, the dynamical evolution of the system will be driven by the stellar evolution of the evolved RG component. Therefore, here we restrict our discussion for some short term effects, and their observational consequences.

In Fig. 14 we plot the variations of the observable arguments of periastrons (ω_1, 2) and inclinations (i_1, 2) of the inner and outer orbits during the present century. The very fast, dynamically forced apsidal motion can nicely be seen on the ETV curves determined from the 4-year-long Kepler observations (see Fig. 11). An observationally more dramatic effect caused by precession of the orbital plane of the inner binary is illustrated in Fig. 15. As the orbital plane precesses with a half amplitude of around the invariable plane of the triple (having an inclination of ) with a period of P_node ≈ 19.2 yr, the inner binary exhibits eclipses (with varying eclipse depths) during a ≈7.3 yr-long interval, while in the remaining ≈11.9 yr, the system no longer exhibits eclipses. Fortunately, Kepler observations caught the first half of the latest eclipsing session and it led to the discovery of this very exciting system. This eclipsing session ended at the very beginning of 2016, and KIC 7955301 will not exhibit any eclipses till the beginning of 2028. In accordance with our findings, the TESS spacecraft did not observe any remarkable light variations of this target when the original Kepler field was revisited three times from 2019 to 2022 (Fig. 16).

Fig. 14. Evolution of some of the orbital elements on a century-long timescale, obtained via numerical integrations. Left: variations of the observable arguments of periastrons (ω1, 2) of the inner and outer orbits (black and red, respectively). The dynamically forced apsidal motions of both orbits are well visible. The gray and orange shaded regions mark the intervals of the Kepler and the ground-based spectroscopic follow up observations, respectively. Right: inclinations of the inner and outer binaries (black and red curves, respectively). The green- and lightblue-shaded horizontal areas denote the inclination domains of the inner and outer orbits where eclipses can occur. The vertical gray area represents the domain of the Kepler observations, while the thin, black, vertical line denotes the interval of the 27-day-long TESS observations (Sector 14). As one can see, during these measurements the inner inclination (i1) was far below the green line, and therefore the inner pair did not present any eclipses. It should be noted that the dark blue line simply stands for i1, 2 = 90° to guide the eye.

Fig. 15. Photodynamical model light curve of KIC 7955301 for the present century. The ∼19.2 yr-period precession cycle is clearly visible. Within a cycle the inner pair exhibits regular eclipses during a ∼7.3 yr-long interval. The last eclipsing session had finished by January 2, 2016 with a short (∼80 ppm) amplitude fading, while the next session is expected to begin on January 18, 2028. The gray vertical area represents the time of Kepler observations, and TESS Sector 14 observations are also denoted with a thin vertical line.

Fig. 16. TESS observations of KIC 7955301 from 2019 to 2022. Top panel: time series processed with the software package FITSH (Pál 2012). Bottom panel: power spectral density of the times series for Sector 14 only. The vertical dot-dashed red line indicates νmax and the vertical dashed blue line indicates Kepler’s Nyquist frequency.

7. Concluding discussion

7.1. Dynamical versus seismic masses

With this paper, we illustrate the first detailed analysis of a hierarchical triple system that includes an oscillating RG by combining asteroseismology, spectroscopy, and dynamical measurements. This system is an important benchmark for calibrating the measurement of RG masses with asteroseismology. Thanks to the acquisition of 23 high-resolution spectra with ARCES at APO, we were able to first classify the system as an SB1 and monitor the RV shift of the RG component along its orbit. Then, with the dynamical model solution, we could disentangle and co-add the optical spectra to derive the atmospheric parameters of the three components, even though the parameters of the inner binary have a low S/N. By fixing the RG surface gravity to that found by the complete LIGHTCURVEFACTORY model, we were able to determine the RG temperature and metallicity at T_A = 4720 ± 105 K and [M/H] = −0.01 ± 0.12 (Table 2).

We computed the asteroseismic mass of the RG in three different manners. The first and most basic one, which consisted in applying the asteroseismic scaling laws tuned for RG stars by Mosser et al. (2013) from the values of ν_max, Δν, and T_A, led to M_A, SL = 1.27 ± 0.05 M_⊙. The second made use of the (PARAM) approach that fits the global asteroseismic parameters, including ΔΠ₁, and by including the metallicity on a grid of stellar-evolution models. The mass appears to be exactly the same as with the asteroseismic scaling relation: M_A, PARAM = 1.27 ± 0.05 M_⊙. Finally, the most sophisticated approach based on fitting individual oscillation frequencies led to a mass of .

The RG mass that we derived from the purely dynamical modeling – independent from both asteroseismology and stellar evolution models – is . The uncertainty is thus about 2%, which is better than what can be obtained from asteroseismology or stellar evolution models based on atmospheric parameters, luminosity, and parallaxes. This type of system, a solar-like oscillator in an SB1 hierarchical triple system, should definitely be considered for helping to calibrate the asteroseismic masses. About ten of them are available in Gaulme et al. (2013) and Gaulme & Guzik (2019), even though none of them exhibit such clear ETVs, which means that a 2% accuracy on mass may not be achievable with those.

It is interesting to notice that in this system, the asteroseismic mass matches the dynamical one, contrary to what has been observed for most RGs in EBs (Gaulme et al. 2016; Brogaard et al. 2018; Benbakoura et al. 2021), where asteroseismology appeared to overestimate masses by about 15% on average. This puzzling fact may just be a case where the sample of RG oscillators with accurate independent mass estimates is small and the dispersion of the observed overestimation is large, but it could also highlight issues with RGs in EBs where the companion is an MS dwarf star. For an RG in an EB, the mass of the RG is driven by the low S/N Doppler shift that is measured by tracking the absorption lines of the companion star. In most cases, the flux coming from the companion is less than 5% of the total flux. We also note that a good agreement between dynamical and seismic masses was found with the double RG KIC 9246715 (Rawls et al. 2016), where there were no issues of noisy RVs.

In contrast, KIC 7955301 is a special target with respect to the bulk of RGs in EBs because it is the least evolved, and thus the smallest RG (≈5.85 R_⊙) for which we have an independent mass measurement. The radii of all of the RGs in EBs listed in Benbakoura et al. (2021), for example, range from 8 to 14 R_⊙. In addition, a significant fraction of them are red clump stars. So far, no significant overestimation of stellar masses and radii by asteroseismology has been observed for MS and subgiant stars. We may have an intermediate case here, for which asteroseismic measurements are unbiased. To clarify the question, our recommendation is to take new RV measurements of the RGs in EBs listed in Benbakoura et al. (2021) and Gaulme & Guzik (2019) with larger S/Ns, that is, longer exposure times, and with spectrometers that are dedicated to high-precision RV measurements.

7.2. Age of the system

One of the goals of the future ESA PLATO mission (Rauer et al. 2014) is to provide ages of solar-like MS stars hosting exoplanets with an accuracy of 10%. Many papers dealing with stellar physics report ages arising from stellar evolution codes with a similar precision, especially when asteroseismic measurements are available. With this paper, we could have run a single stellar evolution code and provided an age with a relatively low uncertainty. We decided to let four different modelers estimate the age in an independent fashion. We obtained 4.9 ± 0.9 Gyr with PARAM, Gyr with MESA, 4.5 ± 0.2 Gyr with YaPSI, and 4.5 ± 0.2 Gyr with LIGHTCURVEFACTORY, that is, with uncertainties of 18.4, 22.7, 4.4, and 4.4% respectively. Both the YaPSI and LIGHTCURVEFACTORY estimates have a better precision, likely because the three components of the system are included into the optimization process. Ages determined for the RG only (component A), with the help of asteroseismic and classical parameters, have uncertainties of about 20%, despite RGs usually being considered good cases for age estimates thanks to their fast evolution along the RG branch.

Our age estimate thus ranges from 3.3 to 5.8 Gyr, which is quite broad given the unusually vast amount of information that we have. By looking more closely, the age is consistent between the three “standard” approaches that use the classical parameters (mass, metallicity, and effective temperature), and the global asteroseismic parameters (PARAM only). On the contrary, the model based on reproducing the individual oscillation frequencies with MESA stands apart. Even though a deeper understanding of the differences between the age estimates goes beyond the scope of this paper, we can already draw some conclusions. A careful inspection of the posterior PDF of the mass, radius, and age (Fig. 12) shows a good agreement between the mass and radius obtained with PARAM and MESA, whereas the age PDF appears to be bimodal with MESA. In fact, the secondary peak of the age distribution matches the PARAM estimate (at ≈4.5 Gyr). A key aspect of the MESA-model based optimization is that the helium abundance Y and metallicity Z are two free independent parameters of the model. The resulting Y value of 0.29 is quite large, considering the slightly subsolar value found for the metallicity (Z ≈ 0.012) and assuming a linear Y = Y(Z) relation, and that the Sun falls on that relation (see Fig. A.4, Miglio et al. 2021). A difference of about 0.02 in Y with respect to what we have in the grid used in PARAM would mean a difference of ∼15% in age, which, considering also the slight difference in mass, could explain the different age estimate.

7.3. Formation of the system

Physical properties and the orbital configuration carry information about the formation process of the systems. We refer the reader to Borkovits et al. (2022) and references therein, where a connection between dynamical and orbital parameters and formation scenarios is discussed. An outer eccentricity of 0.27 for a 209-day period outer orbit is far from being exceptional. For example, Borkovits et al. (2016) reported 16 triple star candidates with a smaller outer period than that of KIC 7955301 and six of them have larger outer eccentricities. Since then, further highly eccentric systems amongst substantially more compact triply eclipsing triples were found in TESS data (Borkovits et al. 2022).

With an outer semimajor axis less than 10 AU, KIC 7955302 is considered to be a compact triple system. The formation of this type of system involves either a fragmentation of the accretion disk or a disk-mediated capture of the outer component. In both cases the outer star forms by being less massive than the inner binary. Then, thanks to its wider and eccentric orbit, it sweeps out a larger orbit, and thus accretes most of the infalling material from the accretion disk, which eventually leads q_out toward larger values. Although it is expected that the total mass of the inner pair remains larger than the mass of the wide component by the end of the complete formation process, the third component may actually become the most massive star of the system (Tokovinin 2021), which is the case of KIC 7955301.

7.4. Evolution of the system

In the short term, we have shown that the inner binary components eclipse each other over a period of about 7.3 years, and stop eclipsing for the next 11.9 years. The RG will never eclipse the inner binary as seen from the Earth. This is why we were able to study the eclipses of the inner binary from the Kepler data, and why we are not able to do so with TESS. The eclipses should be visible again during the ESA PLATO mission, whose nominal operations should take place between 2027 and 2031. The PLATO observations will certainly help fine-tune the model of the system.

In the long term and in general, stars are likely to merge in close multiple-star systems when the most massive component reaches the tip of the RG branch. This is especially true for low-mass stars that ignite helium only at the tip of the RGB, after reaching very large radii (∼200 R_⊙). RGs may swallow planets or stars, converting their orbital momentum into spin and friction loss. There is indirect observational evidence for it. For example, about 3% of the solar-mass stars on the RGB and 11% of the helium-burning stars observed by Kepler display unusually fast rotation and magnetic fields (e.g., Tayar et al. 2015; Ceillier et al. 2017; Tayar & Pinsonneault 2018; Gaulme et al. 2020). This means that a fraction of the RGs in that mass range gain angular momentum between the RGB and the red clump (RC). Until alternative explanations are provided, a likely possibility is that these stars have engulfed a stellar or substellar companion. Recently, Price-Whelan et al. (2020) reported a notable dearth of close companions around the RC – where companions may have been engulfed when the primary star ascended the upper RGB – in the color-magnitude diagram produced by the APO Galactic Evolution Experiment (APOGEE, Majewski et al. 2017), which is part of the Sloan Digital Sky Survey IV (e.g., Blanton et al. 2017). The orbital configuration of KIC 7955301 represents a borderline case for this type of scenario. The current distance at periastron between the giant and the inner binary is about ≈142 R_⊙ according to the data from Table 4.

We can formulate some broad predictions for the long-term evolution of the system, assuming that the dynamical interaction among the components can be neglected until they come into contact with each other because of post-MS evolution. Figure 17 shows the evolution of the radii of the three components according to evolutionary tracks generated with the MIST web interpolator (Choi et al. 2016), setting the masses to the central values given in Table 4 (1.29, 0.94, and 0.88 M_⊙ for components A, Ba, and Bb, respectively), and the metallicity to +0.06. Significant mass loss takes place during the RGB and asymptotic giant branch (AGB) phases, and as a result, all three stars reach a final mass of about 0.5 M_⊙ by the end of their respective AGB phases. Stellar rotation is not included in these models. At an age of approximately 5 Gyr, component A (already an RG at the present time) reaches the RGB tip, and shortly after (relative to the 10⁹ yr timescale) the AGB tip and thermal pulses phase. The peak radius reached at the tip of the RGB (≈160 R_⊙) is mostly equal to the distance at periastron, while the peak radius on some of the AGB pulses is significantly larger than the semimajor axis of the outer orbit, and thus merging can be expected to occur at one of those times. Should merging of component A with the inner binary be avoided because of some process not considered in this crude analysis, at around 15 Gyr, component Ba will evolve onto the RGB and merge with Bb. In light of this discussion, it is comforting that none of our age estimates is larger than 5 Gyr.

Fig. 17. Top panel: radius evolution of the three components according to the MIST tracks (with masses and metallicity set equal to their nominal values from Table 4). The semimajor axis of the inner binary is indicated by the black dotted line, and the semimajor axis of the outer orbit with dashed lines. Bottom panel: zoom on the RG component around its RGB and AGB phases.

Acknowledgments

This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This work is also based on observations obtained with the Apache Point Observatory 3.5-m telescope, which is owned and operated by the Astrophysical Research Consortium. Part of our spectroscopic observations were done at the Observatoire de Haute-Provence. We acknowledge the technical team for their onsite support as well as the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU for their financial support. P.G. and F.S. were supported by the German space agency (Deutsches Zentrum für Luft- und Raumfahrt) under PLATO data grant 50OO1501. P.G. and J.J. acknowledge NASA grant NNX17AF74G for partial support. The research leading to these results has (partially) received funding from the KU Leuven Research Council (grant C16/18/005: PARADISE) and from the Belgian federal Science Policy Office (BELSPO) through PRODEX grant PLATO. M.B. acknowledges support from the “Deutscher Akademischer Austauschdienst” (DAAD) and the Université Paris Diderot.

References

Appendix A: Radial velocity data

Table A.1 contains the RVs obtained at Apache Point and Haute-Provence observatories.

Appendix B: Oscillation frequencies

Table B.1 displays the oscillation frequencies of the RG component, as performed by co-authors Appourchaux and Mosser.

Appendix C: Eclipse properties

Table C.1 presents the eclipses properties obtained with the LIGHTCURVEFACTORY.

All Tables

All Figures

Fig. 1. Kepler light curve of KIC 7955301. Top panel: relative flux as a function of time, where t = 0 is the first day of data. The blue curve is the stitched light curve, whereas the orange curve has the eclipses removed and is used for asteroseismic analysis. Bottom panel: stitched light curve folded over the inner binary orbital period (15.32 days), where consecutive eclipses are vertically offset by 0.003 to ease the visibility of the transit depth and timing variations.

In the text

	Fig. 2. Broadening function of the 23 spectra taken with ARCES. Spectra were sorted by increasing orbital phase. The position along the y axis corresponds with the orbital phase. Barycentric velocity corrections were included in the BF computation. The RV data are indicated with the red plus symbols and their best-fit Keplerian-orbit model by a dashed line.
In the text

Fig. 3. Sample region of the mean spectra of the individual components of the system obtained by disentangling the APO spectra. Top panel: spectrum of the RG component (blue line) and its best fit obtained with the GSSP code (red) from 5340 to 5450 Å. Bottom panel: spectra of the inner binary system Ba and Bb (blue lines) together with their best-fit synthetic spectra (red lines). The optimal values representing the best fits are listed in Table 2.

In the text

Fig. 4. Oscillation spectrum of KIC 7955301. Top panel: power spectral density of the time series after eclipse removal and filling. The dashed red line represents the stellar background noise: it is the sum of the dashed blue lines (correlated noise Harvey profiles) and the dashed green line (white noise). The solid red line represents the fit of the power spectrum, including the stellar background and the Gaussian envelope of the oscillations. Bottom panel: échelle diagram associated with the frequency spacing Δν = 10.48 μHz. The colored lines indicate the oscillation universal pattern, where blue is l = 0, green l = 1, and red l = 2. Modes of degree l = 3 are visible half way in between the l = 0 and l = 1 ridges. The horizontal dashed line indicates the location of νmax.

In the text

Fig. 5. Oscillation amplitude. Top panel: height of the Gaussian envelope used to model the contribution of the oscillations to the stellar background. The red dot highlights the position of KIC 7955301 with respect to a sample of 4500 RGs analyzed by Gaulme et al. (2020) indicated in blue. The sample splits into three categories: the oscillators with regular amplitude, those with depleted l = 1 modes, and the active RGs that show a global mode suppression (l = 0, 1, 2, 3). KIC 7955301 falls in the regular oscillators. Bottom panel: Amplitude of the largest l = 0 mode of KIC 7955301 compared with those of the 35 RGs in EBs that were studied by Benbakoura et al. (2021). The panel displays the relative difference between expected and measured oscillation amplitudes (%) as a function of orbital period Porb (days). The red line indicates 0, i.e., stars whose oscillations are not detected. The two dashed blue lines represent the region in which relative mode amplitude lies for systems with orbital periods longer than 180 days within two sigma. The size of each symbol represents the amplitude of stellar variability (large means variable, small means not variable), and the gray scale indicates the pulsation mode amplitude (white means no modes; black means a large amplitude).

In the text

Fig. 6. Stretched period échelle diagram for dipole gravity-dominated mixed modes. The top panel shows the work led by co-author Appourchaux with ΔΠ1 = 76.3 s. Fitted ridges with m = {− 1, 1} are represented as red and orange crosses, respectively, while m = 0 is represented as a continuous white line. The black stripes show the location of the l = 0, 2 mode power, which has been removed for clarity. The bottom panel summarizes the findings of co-author Gehan, where peaks with a height-to-background ratio equal to or above ten are represented. The symbol size varies with the measured power spectral density. Ridges with m = {− 1, 0, 1} are represented in green, light blue, and red, respectively. Red and green crosses represent the fit of the ridges with m = ±1. The m = 0 ridge is identified by considering that it is located at the mid-point between the m = ±1 ridges.

In the text

	Fig. 7. Peak-bagging results, including mixed-mode analysis by co-author Appourchaux, displayed in the form of an échelle diagram of the power spectrum for Δν = 10.44 μHz. Fitted frequencies are shown: l = 0 by white diamonds, l = 2 by orange diamonds, l = 3 by red diamonds, l = 1, m = −1 by white crosses, l = 1, m = 0 by white circles, and l = 1, m = +1 by white pluses.
In the text

Fig. 8. Peak-bagging results, including mixed-mode analysis by co-author Mosser, in the form of a line-by-line échelle diagram, for radial pressure orders eight to 13. Radial and quadrupole modes are highlighted in red and green, respectively. The expected locations of the dipole mixed modes are labeled with their mixed radial orders, and are indicated with a color depending on the azimuthal order: dark blue (m = −1), light blue (m = 0), or purple (m = +1). ℓ = 3 modes, not identified in this plot, are located near the abscissa 0.21. The gray dashed lines correspond to height-to-background ratios of seven and ten.

In the text

Fig. 9. Long-cadence Kepler light curve of KIC 7955301 together with the joint photodynamical model light curve solution. Upper left panel: complete Q0–Q17 light curve (blue) with the best-fit solution (red). This solution was obtained by neglecting the Doppler-boosting effect. The gray curve shows the same model after “switching on” the Doppler-boosting. While this effect, in theory, should have been observable with the accuracy of Kepler photometry, the reasons for neglecting it are discussed in the text. Upper right panel: one-month-long section of the Kepler observations overplotted with the photodynamical solution. The pale blue circles represent each individual observation but, for the analysis, only dark blue data (located within the phase environment of each eclipses) were considered. Lower panels show the residual data.

In the text

Fig. 10. Radial-velocity data with their best fit models. Left panel: radial-velocity data (red) as a function of time expressed in modified Julian dates, with their best fit model (blue line). Right: same folded over the RG orbital period and displayed as a function of orbital phase. The dispersion of the best-fit curve in blue is caused by the varying shape of the RV curve as a function of time.

In the text

Fig. 11. Eclipse timing variations of KIC 7955301 on different timescales: observations vs model and predictions. Left: primary and secondary ETV curves derived from Kepler observations (red circles and blue boxes, respectively), together with the photodynamical model solution (smaller, pale red and blue symbols, connected with straight lines, respectively). Right: eclipse-timing-variation-like curves derived from the times of the inferior and superior conjunctions of the secondary of the inner pair (pale red and blue lines), overplotted with forecasted primary and secondary ETV points (larger red and blue symbols) during the eclipsing phases of the inner binary for the first half of the present century. The dynamically forced short-period apsidal motion of the inner binary is nicely visible. The gray-shaded area stands for the interval of the Kepler measurements.

In the text

	Fig. 12. Mass, radius, and age posterior probability distribution functions derived from seismic constraints. Empty histograms show the distributions obtained using PARAM, and the dotted lines show their associated ±1σ credible regions. Solid histograms and shaded regions indicate the same quantities obtained using the individual mode frequencies.
In the text

Fig. 13. Visual fit of observed classical parameters of the system with theoretical tracks and isochrones from the YaPSI database (Spada et al. 2017). Left: comparison in the mass radius plane with isochrones of different ages for the primary (top panel) and the two components of the inner binary (bottom). Right: comparison in the radius-luminosity plane; a 4.5-Gyr isochrone and evolutionary tracks of appropriate mass and metallicity are shown.

In the text

Fig. 14. Evolution of some of the orbital elements on a century-long timescale, obtained via numerical integrations. Left: variations of the observable arguments of periastrons (ω1, 2) of the inner and outer orbits (black and red, respectively). The dynamically forced apsidal motions of both orbits are well visible. The gray and orange shaded regions mark the intervals of the Kepler and the ground-based spectroscopic follow up observations, respectively. Right: inclinations of the inner and outer binaries (black and red curves, respectively). The green- and lightblue-shaded horizontal areas denote the inclination domains of the inner and outer orbits where eclipses can occur. The vertical gray area represents the domain of the Kepler observations, while the thin, black, vertical line denotes the interval of the 27-day-long TESS observations (Sector 14). As one can see, during these measurements the inner inclination (i1) was far below the green line, and therefore the inner pair did not present any eclipses. It should be noted that the dark blue line simply stands for i1, 2 = 90° to guide the eye.

In the text

Fig. 15. Photodynamical model light curve of KIC 7955301 for the present century. The ∼19.2 yr-period precession cycle is clearly visible. Within a cycle the inner pair exhibits regular eclipses during a ∼7.3 yr-long interval. The last eclipsing session had finished by January 2, 2016 with a short (∼80 ppm) amplitude fading, while the next session is expected to begin on January 18, 2028. The gray vertical area represents the time of Kepler observations, and TESS Sector 14 observations are also denoted with a thin vertical line.

In the text

	Fig. 16. TESS observations of KIC 7955301 from 2019 to 2022. Top panel: time series processed with the software package FITSH (Pál 2012). Bottom panel: power spectral density of the times series for Sector 14 only. The vertical dot-dashed red line indicates νmax and the vertical dashed blue line indicates Kepler’s Nyquist frequency.
In the text

	Fig. 17. Top panel: radius evolution of the three components according to the MIST tracks (with masses and metallicity set equal to their nominal values from Table 4). The semimajor axis of the inner binary is indicated by the black dotted line, and the semimajor axis of the outer orbit with dashed lines. Bottom panel: zoom on the RG component around its RGB and AGB phases.
In the text