© The Authors 2025

1. Introduction

Triple-star systems represent the simplest multi-star configurations beyond binaries, offering a valuable opportunity to compare stellar model results with observational data. When the third body in a triple system produces eclipses, significant information can be gleaned. These events occur when the inner eclipsing binary occults the third star or vice versa. By analysing the eclipsing binary light curves, extracted measurements of eclipse timing variations, and energy distribution data from archival surveys, it is often possible to determine key stellar parameters, including masses, radii, and effective temperatures, along with orbital parameters. This can be achieved with limited reliance on extensive radial velocity (RV) measurements (see e.g. Carter et al. 2011; Alonso et al. 2015; Rappaport et al. 2023, 2024). Precision of a few percent in mass and radius determinations is often possible for these systems.

While the precision of these parameters may be lower than that achievable for eclipsing binaries through extensive spectroscopic analysis, triple systems are still a valuable test bed for stellar models. Testing and calibrating stellar models with binary systems often demands exceptional observational precision, typically below 1% (see e.g. Valle et al. 2017; Miller et al. 2020; Hełminiak et al. 2021; Anders & Pedersen 2023). The additional constraints provided by a third star in triple systems can aid in refining estimates of stellar evolutionary parameters, even with lower precisions in mass and radius determinations. Studies of triply eclipsing triple systems (e.g. Carter et al. 2011; Feiden et al. 2011; Derekas et al. 2011) exemplify the valuable role a third star can play in refining system solutions.

The census of triply eclipsing triple systems has steadily grown in the last decade. Rappaport et al. (2024) present observations of seven triply eclipsing triple systems. Their analysis leveraged Transiting Exoplanet Survey Satellite data along with archival data from different sources. System parameters were extracted via a photodynamical analysis using LIGHTCURVEFACTORY (Borkovits et al. 2019). No spectroscopically determined temperatures or metallicities were available for the systems; therefore, absolute temperatures and radii were estimated in part using spectral energy distribution fitting. The parameters are, however, astrophysical model dependent because PARSEC isochrones (Bressan et al. 2012) were used as part of input information for the fit to determine the absolute radii and metallicities of the systems. While the mass, radius, and effective temperature determinations were obtained without RV measurements, it has been suggested (Borkovits et al. 2022; Rappaport et al. 2024) that the available excellent eclipse timing variation measurements make their estimation extremely robust.

Investigating these triply eclipsing systems is interesting for testing the reliability of this claim under different scenarios. Stellar models are expected to show the best agreement for low-mass main-sequence (MS) stars that lack a convective core (Valle et al. 2013; Stancliffe et al. 2015). The underlying assumptions within the models have a minimal impact in these cases. However, for more evolved and massive stars, the treatment of convective core overshooting significantly affects their evolution. Consequently, the choice of the astrophysical model becomes increasingly important in determining the system’s parameters for these stars.

Among the Rappaport et al. (2024) systems, two offer a good opportunity to test stellar models accuracy and to check the impact on the fit of the convective core overshooting parameter. A first system, TIC 650024463, is composed of three MS stars with clearly different masses in the range [0.68, 0.88] M_⊙. The system is nearly flat and characterised by an inner binary period of about 7.2 d (semi-major axis of the orbit of 18 R_⊙) and an outer period of about 109 d (semi-major axis of the tertiary orbit of 127 R_⊙). The outer eccentricity is e ≈ 0.33. This system provides a good opportunity to check the reliability of simultaneous stellar model predictions for stars without a convective core.

The second system, TIC 323486857, hosts more massive stars in the range [1.2, 1.6] M_⊙; two are MS stars, and the third, the most massive, is on the red giant branch (RGB). The system has a short outer orbital period of 41 d, a low outer eccentricity e = 0.0066, and is nearly flat. The period of the inner binary is 0.88 d, with a semi-major axis of 5.4 R_⊙, while the semi-major axis of the tertiary orbit is about 81 R_⊙. This system is a good target to investigate how the underlying astrophysical assumptions on the convective core overshooting efficiency influences the system fit when different stellar models are adopted. While a more comprehensive analysis of a larger sample of triply eclipsing triple systems would be ideal, the computational burden imposed by the fit of these system, partly due to the observational uncertainties being significantly greater than those for double-lined eclipsing binary systems, is so high that we were forced to restrict the investigation to these two specific targets. The dimensionality of the stellar model grid required to accurately fit a system is primarily determined by the observational precision of the component masses. For well-characterised double-lined eclipsing binaries with extensive RV measurements, these uncertainties are often of the order of 0.001 M_⊙, allowing the stellar mass to be effectively fixed in the grid computation. This is not the case when observational errors exceed 0.01 M_⊙, as a broader mass range encompassing approximately 3σ around the observational constraint must be covered, necessitating a significantly larger grid.

The structure of the paper is as follows. In Sect. 2 we present the grid of stellar models and the fitting method used in the estimation process. The results for both systems are presented in Sect. 3. Concluding remarks can be found in Sect. 4.

2. Methods

2.1. Stellar model grid

The grid of stellar evolutionary models was calculated for the mass range 0.60 to 1.00 M_⊙ for the first system and 1.10 to 1.70 M_⊙ for the second one, spanning the evolutionary stages from the pre-MS to the RGB tip. The initial metallicity [Fe/H] was varied from −0.4 dex to 0.4 dex with a step of 0.02 dex. We adopted the solar heavy-element mixture by Asplund et al. (2009). For each metallicity, we considered a range of initial helium abundances based on the commonly used linear relation $Y=Y_p+\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}Z$ with the primordial helium abundance Y_p = 0.2471 from Planck Collaboration VI (2020). The helium-to-metal enrichment ratio ΔY/ΔZ was varied from 1.0 to 3.0 with a step of 0.25.

The models were computed with the FRANEC code, in the same configuration as was adopted to compute the Pisa Stellar Evolution Data Base¹ for low-mass stars (Dell’Omodarme et al. 2012). The only difference with respect to those models is that the outer boundary conditions were set by the Vernazza et al. (1981) solar semi-empirical T(τ), which approximate well results obtained using the hydro-calibrated T(τ) (Salaris & Cassisi 2015; Salaris et al. 2018). The models were calculated assuming the solar-scaled mixing-length parameter α_ml = 2.02. The extension of the extra-mixing region beyond the Schwarzschild border was parametrised in terms of the pressure scale height H_p: l_ov = βH_p, with β from 0.00 to 0.28 with a step of 0.01. The code adopts step overshooting assuming an instantaneous mixing in the overshooting treatment. The radiative temperature gradient is adopted in the overshooting region (see Degl’Innocenti et al. 2008, for more details of the overshooting implementation). Atomic diffusion was included adopting the coefficients given by Thoul et al. (1994) for gravitational settling and thermal diffusion. To prevent extreme variations in the surface chemical abundances, the diffusion velocities were multiplied by a suppression parabolic factor that takes a value of 1 for 99% of the mass of the structure and 0 at the base of the atmosphere (Chaboyer et al. 2001).

Raw stellar evolutionary tracks were reduced to a set of tracks with the same number of homologous points according to the evolutionary phase. Details about the reduction procedure are reported in the appendix of Valle et al. (2013). Overall, to achieve the required grid density for all the investigated parameters, more than 10⁵ stellar tracks were computed.

2.2. Fitting technique

The analysis was performed adopting a modified SCEPtER pipeline², a well-tested technique for fitting single and binary systems (e.g. Valle et al. 2017, 2023b). The pipeline estimates the parameters of interest (i.e. the system age, its initial chemical abundances, the convective core overshooting parameter) adopting a grid maximum likelihood approach.

The method we used is explained in detail in Valle et al. (2015); here, we discuss only the modifications needed to fit a triple system. For every j-th point in the fitting grid of pre-computed stellar models, a likelihood estimate is obtained for all stars,

$$\begin{array}{c}\hfill {{\mathcal{L}}^{1,2,3}}_j=\left({\displaystyle \prod _{i=1}^n}\frac{1}{\sqrt{2\pi }{\sigma }_i}\right)\times exp(-\frac{{\chi }^2}{2}),\end{array}$$(1)

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^n}{\left(\frac{o_i-g_i^j}{{\sigma }_i}\right)}^2,\end{array}$$(2)

where o_i are the n system constraints (masses, radii, effective temperatures and metallicities for the three stars), g_i^j are the j-th grid point corresponding values, and σ_i are the respective uncertainties. The adoption of the χ² statistic in Eq. (2) to judge the goodness-of-fit is widely used. However, careful consideration is necessary when evaluating the associated degrees of freedom (dofs). In our case, each star is characterised by four parameters, but these are not independent due to the constraints imposed by stellar evolution tracks. Valle et al. (2021) demonstrate that when parameters are linked by an isochrone (and similarly by a stellar track), one dof is lost for each object. This, however, assumes independent observations for each star, which is not the case for multiple systems. Therefore, accurately determining the actual number of dofs is challenging. In the following, we assume three dofs for each star, and account for a loss of one dof for each estimated parameter, such as initial chemical composition, age, and β.

The joint likelihood of the system is then computed as the product of the three single star likelihood functions. The whole system fit is obtained by imposing that three objects must have a common age (with a tolerance of 1 Myr), identical initial helium abundance, and initial metallicity. This step gives a pool of models with corresponding likelihood values ℒ_p. For these models the temperature (q_T; 1, 2 and q_T; 1, 3) and masses ratio (q_1, 2 and q_1, 3) are computed. The corresponding fitting values ${\stackrel{\sim }{q}}_{T;1,2}$, ${\stackrel{\sim }{q}}_{T;1,3}$, ${\stackrel{\sim }{q}}_{1,2}$, and ${\stackrel{\sim }{q}}_{1,3}$ are accounted for as follows:

$$\begin{array}{c}\hfill \mathcal{L}={\mathcal{L}}_{p}exp(-\frac{{\tau }^2}{2})\end{array}$$(3)

$$\begin{array}{cc}\hfill {\tau }^2& ={\left(\frac{q_{T;1,2}-{\stackrel{\sim }{q}}_{T;1,2}}{{\sigma }_{T;1,2}}\right)}^2+{\left(\frac{q_{T;1,3}-{\stackrel{\sim }{q}}_{T;1,3}}{{\sigma }_{T;1,3}}\right)}^2+{\left(\frac{q_{1,2}-{\stackrel{\sim }{q}}_{1,2}}{{\sigma }_{1,2}}\right)}^2+\hfill \\ \hfill & +{\left(\frac{q_{1,3}-{\stackrel{\sim }{q}}_{1,3}}{{\sigma }_{1,3}}\right)}^2.\hfill \end{array}$$(4)

Since accounting for the radii ratios did not alter the solutions of the studied systems, we excluded it from the final algorithm. Stellar parameters for models with a likelihood ℒ exceeding 95% of the maximum likelihood value³ are averaged and returned as the solution. The error in the fit parameter was evaluated by means of a Monte Carlo procedure, by simulating 5000 artificial triple systems, These systems were obtained by Gaussian perturbations of the observational values, taking the correlation between temperatures, metallicities, and masses into account.

3. Results

3.1. TIC 650024463

Constraints for the TIC 650024463 system, derived from Rappaport et al. (2024), are presented in Table 1. In the absence of a spectroscopic determination, errors in effective temperature and metallicity, [Fe/H], are increased to 100 K and 0.1 dex, respectively. Errors in the masses and radii were obtained by computing the geometric means of the upper and lower error bars given by Rappaport et al. (2024).

As often happens with binary systems (e.g. Valle et al. 2017, 2023b), the posterior distribution of the system parameters exhibits a multimodal density. This is evident in the left panel of Fig. 1, which shows the probability density in the initial helium abundance Y versus age plane. A first solution, S1, peaks at low initial helium abundance Y ≈ 0.26 and an age of $9.{35}_{-1.18}^{+1.42}$ Gyr, while a second solution, S2, has a higher initial helium abundance Y ≈ 0.28 and a slightly younger age of $8.{76}_{-1.07}^{+1.26}$ Gyr. Both solutions provide a good fit of the system, with χ² values of 1.65 and 0.78, respectively, with six dofs. The small values of the computed χ² are primarily attributable to the relatively large uncertainties associated with the system constraints. Specifically, the 1% precision in the stellar masses and radii not only contributes to a better agreement with the proposed system solution but also affords the algorithm greater flexibility to adjust the fitted masses and radii, improving the agreement with the temperatures and metallicities constraints, thereby enhancing the overall goodness-of-fit. Detailed fit results are collected in Table 2. The central panel in Fig. 1 shows that the isochrones resulting from the fit are quite similar. As expected from theoretical investigations (e.g. Valle et al. 2015) the age of the system is mainly set by the more evolved primary star. This is shown in the right panel in Fig. 1 where the kernel density estimators of the age of the stars, from independent fit neglecting the coevality constraints, are shown. The density of the primary star is clearly less dispersed than the other two distributions, leading to an estimated age for the single primary star of $9.{15}_{-1.19}^{+1.40}$ Gyr.

Fig. 1. Fit of the TIC 650024463 system. Left: 2D kernel density estimator of the posterior distribution of the system parameters in the Y vs. age plane. Middle: Best fitting isochrones in the mass vs. radius plane. The points indicate the derived system constraints, while solid red and blue lines correspond to S1 and S2, respectively. Right: Kernel density estimators of the age of the three stars from independent fits.

The posterior density in Fig. 1 demonstrates the well-known degeneracy between age and initial chemical composition, a common challenge in MS stars age calibration. The probability density gradually decreases away from the peaks, indicating that multiple combinations of age and initial chemical composition can yield nearly equivalent fits for this system. Assuming coevality, the averaged age of the system, $9.0_{-1.1}^{+1.4}$ Gyr, is robust because the constraints provided by the three stars with clearly different masses narrow the age range compatible with all of them. The fairly large variability in other parameters, such as the initial helium abundance, only slightly affects the global age. This aligns with the findings of Valle et al. (2024a) in binary systems. As expected the result of the fit agrees with the underlying astrophysical assumption adopted in Rappaport et al. (2024), who obtained the system parameters relying on a PARSEC isochrone with age $10.0_{-0.7}^{+0.9}$ Gyr. It is relevant to note that the 14% precision of our age estimates was obtained with mass constraints at 2% level. Adopting the same methods, we obtained more precise estimates for binary systems, such as TZ For (6%) or CPD-54 810 (4%; Valle et al. 2017, 2023a), but they required masses and radii observational uncertainties well below 1%.

3.2. TIC 323486857

Constraints for the TIC 323486857, derived from Rappaport et al. (2024), are presented in Table 3. As for the previous system, the errors in effective temperature and [Fe/H] were fixed at 100 K and 0.1 dex, respectively.

As anticipated, the solution for this system exhibits some limitations and requires supplementary assumptions. As a first fit attempt we tried to constrain the system to have a common core overshooting efficiency for all the stars. This assumption will be further challenged by performing different fits relaxing this as well as other constraints.

As for the first system, the fit assuming common convective core overshooting efficiency has a bimodal posterior density, but the two solution islands for TIC 323486857 are well separated (left panel in Fig. 2). Two solutions are found: S1 at a low β = 0.02 ± 0.02, and S2 at the edge of the grid, β = 0.28. The first solution is significantly more probable, and the second solution should be viewed with caution. This is because a broader range of allowed overshooting values in the grid would likely have resulted in a different optimal solution. Table 4 summarises the full fit parameters, including the convective core masses M_cc for the primary and secondary stars (the tertiary has already reached the RGB evolutionary phase, where the convective core disappears). Both solutions yield high χ² values of around 8.5 (with five dofs), largely due to a poor fit for the tertiary star. With a precision of about 4% in the stellar masses and radii a better agreement is expected, as discussed for the first system. Focusing solely on the inner binary system improves the goodness-of-fit, with χ_1, 2² values around 2. Overall, we get a system age of $2.{29}_{-0.16}^{+0.17}$ Gyr. The right panel of Fig. 2 highlights the problematic fit for the tertiary. Its fitted position on the isochrone deviates significantly in effective temperature, with the observational constraint being nearly 150 K hotter than the model prediction. This discrepancy could stem from an inaccurate representation of effective temperature in the RGB phase, potentially due to a not perfect treatment of the superadiabatic convection efficiency (Trampedach et al. 2014; Magic et al. 2015). Notably, the PARSEC isochrones used by Rappaport et al. (2024) closely matches the effective temperature of our isochrones in the RGB phase, suggesting that the discrepancy with this constraint of the system is not a peculiarity of our models. Alternatively, systematic errors in the system derived parameters could also explain the discrepancy. However, definitive confirmation requires spectroscopic data to constrain the system’s metallicity and effective temperature.

Fig. 2. Fit of the TIC 323486857 system assuming the same core overshooting efficiency for all stars. Left: 2D kernel density estimator of the posterior distribution of the system parameters in the age vs. β plane. Right: Best fitting isochrone in the radius vs. effective temperature plane. The open dots indicate the derived system constraints, while the solid red and blue lines correspond to S1 and S2, respectively. The triangle identifies the position of the fitted tertiary star for S1.

An interesting difference exists between our fit and the derived radius of the tertiary star under both solutions. This difference can be attributed to the choice of stellar models used for the observational characterisation of the system by Rappaport et al. (2024). Specifically, the PARSEC isochrones utilise a different overshooting scheme compared to our models. In fact the dependence of the core overshooting parameter on the stellar mass is still debated and got relevant attention in the literature, with opposite claims (see Anders & Pedersen 2023, for a review). PARSEC isochrones adopt an increasing core overshooting parameter with stellar mass, starting at 1.0 M_⊙ and reaching a maximum value at 1.5 M_⊙. Moreover, PARSEC parametrisation of the overshooting mechanism is different than our, their maximum values corresponding roughly to β = 0.25 in our implementation (Bressan et al. 2012). This choice in the core overshooting efficiency impacts the predicted radius of the tertiary star, leading to a larger value in the Rappaport et al. (2024) system solution compared to ours, which assumes a constant core overshooting efficiency across all stellar masses.

To further investigate this point, we performed a fit where the algorithm was allowed to choose independent β values for each star, with the only constraint being a common initial chemical composition and age. Figure 3 shows that the solution for the tertiary star peaks at the edge of the grid, while the primary star is well fitted with a low β ≈ 0.02. The poorly evolved secondary star does not provide significant constraints on the overshooting parameter. As in the previous analysis, the β for the tertiary stars should be viewed with caution. This is because a better solution may exist at higher overshooting values beyond the range explored in the grid. The system age estimate of $2.{20}_{-0.14}^{+0.19}$ Gyr closely matches that obtained in the common overshooting scenario. However, the goodness-of-fit for this system remains questionable, χ² = 6.1 (with three dofs) mainly because of the effective temperature tertiary contribution (χ_1, 2² = 1.8). While the high β value for the tertiary star aligns with what is assumed in the PARSEC isochrones, it is surprisingly high for a 1.6 M_⊙ star. For comparison with the literature, adopting the same methodology for binary systems, a maximum estimated overshooting parameter of β ≈ 0.17 was obtained for the more massive TZ For system (Valle et al. 2017). Independently, Claret & Torres (2016) suggested that for a star of mass 1.6 M_⊙ a value of β ≈ 0.1 is expected.

Fig. 3. Kernel density estimators of the individual β parameters for the three stars in the TIC 323486857 system. Solid lines represent the fit that adopts the grid with variable surface [Fe/H] and a 100 K uncertainty in the tertiary effective temperature. Dashed lines were obtained by fixing the track grid surface [Fe/H] to the initial values and adopting 600 K as the uncertainty in the tertiary effective temperature.

Overall, none of the analysed fits provide satisfactory agreement with the constraints from Rappaport et al. (2024). Furthermore, they are potentially affected by several limitations. First, a poor fit of the tertiary effective temperature, as previously discussed, may significantly influence the results. The discrepancy between theoretical models and the observed temperature can distort the fit, forcing the algorithm to stretch the solution towards the grid edge to accommodate this constraint. Second, the assumption of a common metallicity for all three stars, derived from PARSEC isochrones, may not be accurate because the stars are in different evolutionary phases. As extensively discussed in Valle et al. (2023b), the tertiary star’s surface metallicity is expected to be close to its initial value. However, the primary and secondary stars, still on the MS, may experience changes in surface metallicity due to competing mixing processes. The adopted grid includes microscopic diffusion, and its effects are evident in the spread of S1 and S2 best-fit [Fe/H] values in Table 4, exceeding 0.1 dex.

Given the limitations in both the isochrone-derived metallicity and the theoretical predictions of surface [Fe/H] for stars with thin convective envelopes (e.g. Moedas et al. 2022), we performed a supplementary fit to assess the robustness of our results. This fit incorporates a larger 600 K uncertainty in the tertiary effective temperature and fixes the grid surface [Fe/H] for all stellar tracks to the initial value. Therefore, this grid still includes the effect of microscopic diffusion on the evolutionary timescale but neglects its impact on the surface chemical abundances. The results of this fit, adopting independent overshooting for the three stars, are presented in Table 4 under the S3 scenario. The posterior densities of the convective core overshooting parameters are shown in Fig. 3, demonstrating agreement with the results discussed previously. The age of the system $2.{33}_{-0.16}^{+0.18}$ Gyr is in excellent agreement with that determined in the previous scenario. This analysis suggest a limited impact of alternative decisions on the recovery of the main system parameters.

Moreover, unlike previous solutions, scenario S3 provides a value for the tertiary β parameter (0.26) within the grid, supporting its reliability. Overall, scenario S3 exhibits a satisfactory fit with a χ² value of 1.9, partially attributed to the relaxed constraint on the tertiary T_eff, but also to a better fit for the tertiary radius. However, since the models are not nested, a formal significance test comparing the likelihoods of different scenarios is not theoretically justified (see e.g. Härdle & Simar 2012; Feigelson & Babu 2012).

Encouragingly, all our age estimates and that adopted by Rappaport et al. (2024) in their system solution, $2.{37}_{-0.20}^{+0.07}$ Gyr, agree very well. This consistency suggests that the system age is a robust determination, even when using different stellar models in the fitting process, provided a sufficiently large parameter grid is employed. However, this finding may not apply to systems with stars in different evolutionary phases. In this specific case, the age is mainly dictated by the inner binary system, in particular by the primary star, which contributes with more weight than the secondary to the age estimate. An analysis restricted to the inner binary system shows an age of $2.1_{-0.24}^{+0.33}$ Gyr, close to the estimate obtained when considering also the tertiary star constraint. While the contribution of the tertiary star to fixing the system age is not determinant, the comparison of the error range of this fit with that reported above highlights its importance in getting a precise estimate.

4. Conclusions

Leveraging recently released data on triply eclipsing triple systems (Rappaport et al. 2024), we investigated the possibility of adopting them to constrain stellar model parameters. In particular, we focused on evaluating the influence of the astrophysical assumptions employed for the observational characterisation of the system by Rappaport et al. (2024).

We utilised dense grids of pre-computed stellar models to fit the data for the triply eclipsing systems with a modified version of the SCEPtER pipeline adapted for triple systems. For the first system, TIC 650024463, which comprises three low-mass MS stars, we achieved an excellent agreement with observational data at $9.0_{-1.1}^{+1.4}$ Gyr. As commonly observed in binary system fitting, the algorithm identified a multimodal posterior density for the fit parameters, suggesting the presence of two plausible solutions with differing initial helium abundances. A good agreement with the system constraints derived by Rappaport et al. (2024) was expected because the difference between the theoretical stellar models adopted by Rappaport et al. (2024) and the present ones used for the system fit are minor for low-mass MS stars (Valle et al. 2013; Stancliffe et al. 2015).

Fitting the second system, TIC 323486857, proved more challenging. This system comprises stars with masses ranging from about 1.2 M_⊙ to 1.6 M_⊙; the tertiary and most massive star is in the RGB phase. We tested three fitting scenarios: (i) assuming a constant efficiency for convective core overshooting across all stellar masses, (ii) allowing independent overshooting efficiencies for each star, and (iii) assuming a different effective temperature precision for the tertiary star and a different efficiency for the mixing processes regarding the evolution of the surface chemical composition, keeping the independent overshooting efficiencies for each star.

In the first scenario, with a constant overshooting efficiency, the fit for the tertiary star’s parameters, particularly its effective temperature and radius, displayed significant discrepancies with the observed data. Notably, our derived effective temperature aligns with the astrophysical models used in Rappaport et al. (2024) to determine the system parameters. This suggests that the discrepancy is not attributable to the underlying difference in the RGB temperature scale between our models and PARSEC isochrones. As a matter of fact, systematic differences in the effective temperature between models and observations in the RGB phase are reported in literature (see Valle et al. 2024b, and references therein for a recent analysis). Further investigation on this topic will require spectroscopic data for the system. The discrepancy in the tertiary star’s radius is particularly intriguing, as it likely arises from differing assumptions regarding the convective core overshooting efficiency between our models and those employed by Rappaport et al. (2024). Their analysis utilises PARSEC isochrones, which incorporate an increasing convective core overshooting efficiency with stellar mass and imply a larger radius for the RGB star.

When we allowed the overshooting parameter to vary freely, as in second and third scenarios, we achieved a close match to the assumptions of the PARSEC models. A key difference between these scenarios is that the second includes the effect of microscopic diffusion on the surface metallicity, while the third one does not. The third scenario offers the most robust system parameter determination given the fact that the metallicity constraint is derived from an isochrone and not from spectroscopic individual determinations. PARSEC models assume a convective core overshooting parameter of β ≈ 0.25 for the tertiary star, and our analysis of the third scenario yielded β ≈ 0.26. While this value aligns well with the PARSEC model assumptions, it is notably higher than what is typically reported in the literature for a star of 1.6 M_⊙. For instance, Claret & Torres (2016) propose a value of β ≈ 0.1. Furthermore, the third scenario, which mimics the assumptions of PARSEC isochrones as closely as possible, provides a clean fit of the system, significantly improving the tertiary radius fit. This finding suggests that the underlying astrophysical assumptions in the Rappaport et al. (2024) derivation of the system constraints exert a dominant influence on the fitted β parameter. While analysing other triply eclipsing triple systems would be extremely valuable for assessing the generality of this result, the substantial computational cost of fitting these systems limited our investigation. Further independent research on this topic could provide deeper insights into the impact of astrophysical assumptions on the fitting of triply eclipsing triple systems. Encouragingly, despite the relevant differences in the adopted stellar models, our most robust age determination, $2.{33}_{-0.16}^{+0.18}$ Gyr, agrees well with that from the models in Rappaport et al. (2024).

In conclusion, our independent analysis suggests a minimal influence of the astrophysical models used by Rappaport et al. (2024) to obtain stellar system data on the derived parameters for MS stars, particularly those of low mass. However, the suitability of these systems for testing stellar evolution in more massive and evolved stars is still an open question. Encouragingly, the system ages derived here show good agreement with those used in the Rappaport et al. (2024) system parameter determination, indicating a relatively robust age determination. While this conclusion may not hold for systems with stars in different evolutionary stages, it appears that triply eclipsing triple system data, obtained without RV measurements, hold promise to constrain stellar ages, provided a sufficiently extensive parameter grid is employed during the fitting process.

Acknowledgments

We thank the referees for the insightful comments and constructive suggestions. G.V., P.G.P.M. and S.D. acknowledge INFN (Iniziativa specifica TAsP) and support from PRIN MIUR2022 Progetto “CHRONOS” (PI: S. Cassisi) finanziato dall’Unione Europea – Next Generation EU.

References

All Tables

All Figures

Fig. 1. Fit of the TIC 650024463 system. Left: 2D kernel density estimator of the posterior distribution of the system parameters in the Y vs. age plane. Middle: Best fitting isochrones in the mass vs. radius plane. The points indicate the derived system constraints, while solid red and blue lines correspond to S1 and S2, respectively. Right: Kernel density estimators of the age of the three stars from independent fits.

In the text

Fig. 2. Fit of the TIC 323486857 system assuming the same core overshooting efficiency for all stars. Left: 2D kernel density estimator of the posterior distribution of the system parameters in the age vs. β plane. Right: Best fitting isochrone in the radius vs. effective temperature plane. The open dots indicate the derived system constraints, while the solid red and blue lines correspond to S1 and S2, respectively. The triangle identifies the position of the fitted tertiary star for S1.

In the text

Fig. 3. Kernel density estimators of the individual β parameters for the three stars in the TIC 323486857 system. Solid lines represent the fit that adopts the grid with variable surface [Fe/H] and a 100 K uncertainty in the tertiary effective temperature. Dashed lines were obtained by fixing the track grid surface [Fe/H] to the initial values and adopting 600 K as the uncertainty in the tertiary effective temperature.

In the text