# Grid-based estimates of stellar ages in binary systems

###### Key Words.:

Binaries: eclipsing – methods: statistical – stars: evolution – stars: low-mass – stars: fundamental parameters^{1}

## Abstract

Context:

Aims: We investigate the performance of grid-based techniques in estimating the age of stars in detached eclipsing binary systems. We evaluate the precision of the estimates due to the uncertainty in the observational constraints, , and the systematic bias caused by the uncertainty in convective core overshooting, element diffusion, mixing-length value, and initial helium content.

Methods: We adopted the SCEPtER grid, which includes stars with mass in the range [0.8; 1.6] and evolutionary stages from the zero-age main sequence to the central hydrogen depletion. Age estimates have been obtained by a generalisation of the maximum likelihood technique described in our previous work.

Results: We showed that the typical random error in age estimates – due only to the uncertainty affecting the observational constraints – is about , which is nearly independent of the masses of the two stars. However, such an error strongly depends on the evolutionary phase and becomes larger and asymmetric for stars near the ZAMS where it ranges from about to . The systematic bias due to the including convective core overshooting – for mild and strong overshooting scenarios – is about 50% and 120% of the error due to observational uncertainties. A variation of in the helium-to-metal enrichment ratio accounts for about of the random error. The neglect of microscopic diffusion accounts for a bias of about 60% of the error due to observational uncertainties. We also introduced a statistical test of the expected difference in the recovered age of two coeval stars in a binary system. We find that random fluctuations within the current observational uncertainties can lead genuine coeval binary components to appear to be non-coeval with a difference in age as high as 60%.

Conclusions:

## 1 Introduction

The possibility of obtaining direct measurements of stellar mass and radius makes eclipsing binary systems the objects of fundamental importance (see e.g. Andersen1991; Torres2010). Nowadays, for stars in double-lined systems, these quantities can be determined with precisions of a few percentage points or better (see e.g. Clausen2008; Pavlovski2014). Thus, detached eclipsing binaries offer a unique test bed for evolutionary models of single stars, giving the opportunity of obtaining accurate estimates of their ages.

However, this age estimation is usually performed by fitting the observed data with isochrones or stellar evolutionary tracks (see among many others Southworth2011; Torres2012; Vos2012; Pavlovski2014). Recently, different techniques have been established for single stars, based upon a maximum likelihood estimation over a grid of pre-computed models (e.g. Gai2011; Basu2012; Gennaro2012; Mathur2012; PradaMoroni2012; scepter1; eta). Such methods are very fast and flexible since they can simultaneously take all the available observables into account, reaching high precision in the estimates and self-consistently evaluating their errors.

Despite these appealing features, very few efforts can be found in the literature that address the performance of these techniques for binary systems and that investigate the possible sources of their internal bias. Two studies that address, among other topics, the problem of binary ages are by Gennaro2012 and Schneider2014. The former is a Bayesian investigation of the age reconstruction of pre-main sequence stars, which also considers binary systems. The authors adopted nine synthetic test systems, subjected each of them to 100 Monte Carlo perturbations, and evaluated the fraction of fake non-coeval estimates at level. The results indicated a fraction of erroneous non-coevality from 0% to 95% with a mean of 13%. Given the restricted sample of explored systems, Gennaro2012 recognised the importance of a full analysis of the expected differences in age estimates in different ranges of mass, metallicity, and evolutionary phases. The last was a study by Schneider2014, who tested a Bayesian technique for age estimation and, st other tasks, compared the inferred stellar ages for the primary and secondary stars of eclipsing binaries in the mass range between 4.5 and 28 . In a sample of 18 binary systems, they found agreement at 1 level in the estimated ages of 17 of them. They also evaluated the mean difference in age of 14 systems (systems with mass ratio higher than 0.97 were excluded) to be Myr at the 95% confidence level.

The aim of this paper is to theoretically analyse the accuracy of grid-based age estimates for detached eclipsing binary systems. We also study the bias on these estimates due to the uncertainties in some physical mechanisms adopted in the models (convective core overshooting and element microscopic diffusion efficiencies, mixing-length value), and in the stellar chemical composition (mainly the initial helium content).

Finally, we introduce a statistical test of the expected difference in the reconstructed ages of two coeval binary components simply due to the observational uncertainties. Such a test is of interest because the recovered non-coevality of stars in binary systems – i.e. the inability to fit both components with a single isochrone – is often used to claim some deficiency in current generation of stellar models. This sometimes leads to introducing and/or calibrating some physical processes in evolutionary codes, such as convective core overshooting, or to varying the external convection efficiency through the mixing-length parameter (see e.g. Andersen1991; Pols1997; Ribas2000a; Torres2006; Morales2009; Clausen2009; Clausen2010; Torres2010; Torres2014, and references therein).

We restrict our analysis to central hydrogen-burning stars with mass in the range from 0.8 to 1.6 . A similar investigation has already been performed (eta) for the age determination of isolated single stars for which the main observables adopted in the reconstruction were different from the present ones. In that case neither the mass nor the radius were available, while in the present study we do not rely on the asteroseismic observables. As discussed in the following, the performances and the systematic biases of grid-based recovery techniques sensitively depend on the specific set of observables actually adopted in the recovery procedure.

The structure of the paper is the following. In Sect. 2 we discuss the method and the grids used in the estimation process. The main results are presented in Sects. 3 and 4. In Sect. 5 we present a statistical test of the expected differences between the single-star age estimates. Some concluding remarks can be found in Sect. 6. The Appendix A describes the sampling strategy followed to build the synthetic dataset.

## 2 Grid-based recovery technique

The basic technique, derived from
Basu2012 and suitable for isolated star estimation, is
described in scepter1; eta, hereafter V14 and V15.
Age estimates are based upon a modified SCEPtER (Stellar CharactEristics Pisa
Estimation gRid) scheme, adapted for binary
stars. The code and the standard grid developed for this
work are available in the R package
SCEPtERbinary^{2}

The adopted implementation assumes that and are detached binary system stars for which the following vectors of observed quantities are available: . We let be the nominal uncertainty in the observed quantities. For each point on the estimation grid of stellar models, we define . Let be the single-star likelihood functions defined as

(1) |

where

(2) |

The single-star likelihood functions are independently evaluated for the two stars for each grid point within of all the variables from ; let be the two maximum values obtained in this step. The single-star ages are estimated by averaging the corresponding quantity of all the models with likelihood greater than . Informative priors can be inserted as a multiplicative factor in Eq. (1), as a weight attached to the grid points.

In the following we assume, implicitly or explicitly, that the stars in the binary system are coeval. This assumption justifies the adoption of a single age as representative of the whole system. A first possible way to estimate the age of the binary system is simply to take the mean of the ages of the single components, computed independently. This method has a computational complexity of , where are the sizes of the grid-point samples.

A second approach explicitly assumes in the likelihood computation that the stars are coeval and computes the joint likelihood function only for the couples of models in the two boxes with ages within 10 Myr. The joint likelihood is computed as the product of the single-star likelihood functions. Let be the maximum value obtained in this step. The joint-star estimated age is obtained by averaging the corresponding quantity of all the couples of models with likelihood greater than . This technique has a computational complexity of . Readers interested in the technical solutions implemented are referred to the code provided in the above-mentioned SCEPtERbinary R package. This estimation method is assumed as our standard in the following.

The described technique is similar to those previously employed in the literature in the framework of isochrone fitting for binary age estimations (see, amongst many, Pols1997; Lastennet2002). The main differences are that the quoted approaches adopt a time-consuming numerical functional minimisation as the best model estimates and use profiles to infer the errors on the estimated parameters, often fixing the value of the masses to their observed values. However these methods can lead to possibly severe error underestimation (see e.g. Basu2010; Quirion2010). The technique adopted in this and similar grid-based approaches (e.g. Gai2011; Basu2012) computes as best estimate a local mean of the best-matching models, as described above. The error on the estimates are then computed by a Monte Carlo simulation. This approach allows a confidence interval to be provided for the age estimate and observational parameter correlations to be considered that are not addressed in profiles (see e.g. Lastennet2002; Basu2010).

Although in the current paper we are not interested in obtaining a statistical confidence interval for the age of observed systems, the adopted method can be easily used to this purpose. To address this point,
we reside on a generation of a synthetic
sample of
binary systems, starting from the observed values and following, for each star, a multivariate normal
distribution with vector of mean and
covariance matrix .
For the Monte Carlo simulations, a value of can be adopted since it provides a fair balance between
computation time and accuracy of the results^{3}

### 2.1 Standard stellar model grid

The standard estimation grid of stellar models is obtained using FRANEC
stellar evolution code (scilla2008; tognelli2011), in the same
configuration as adopted to compute the Pisa Stellar
Evolution Data Base^{4}

The grid consists of 141 680 points (110 points for 1 288 evolutionary
tracks),
corresponding to evolutionary stages from the ZAMS to the central hydrogen
depletion^{5}

## 3 Binary system age estimates: internal accuracy

To evaluate the internal accuracy of grid-based age estimates of binary systems, we analysed the ideal case in which stellar models are in perfect agreement with stars. This is the most favourable case where the error in the estimated ages arises only from the uncertainties affecting the observational constraints used in the recovery procedure.

To achieve such an ideal situation, the first age-estimation test was performed on a synthetic dataset obtained by sampling artificial binary systems from the same standard estimation grid of stellar models used in the recovery procedure. To simulate the effect of observational uncertainties, we added a Gaussian noise in all the observed quantities. We assumed as standard deviations 100 K in , 0.1 dex in [Fe/H], 1% in mass, and 0.5% in radius. Since the observationally inferred values of a given physical quantity for the two binary components are correlated, correlation coefficients should be used in the covariance matrix whenever a realistic noise has to be simulated. In our standard scenario, we assumed a correlation of 0.95 between the primary and secondary effective temperatures, 0.95 between the metallicities, 0.8 between the masses, and no correlation between the radii. Motivations of these choices and an analysis of the impact on the results of different assumptions of correlations amongst stellar quantities, radii included, are presented in Sect. 3.1.

Figure 1 shows the Monte Carlo relative error
in estimated age as a
function of the mass and the relative age of the primary star. In this and in
the
following analogous figures, a positive
value of the age relative error corresponds to an overestimated age.
The relative age is defined as the ratio between
the age of the star and the age of
the same star at central hydrogen exhaustion (the age is conventionally set to
0 at the ZAMS position).
It shows the position of the median
relative error and the position of the envelope of age-relative error as a function of
mass and relative age. The envelope is obtained as in V15,
computing the
16th and 84th quantiles of the
relative errors over a moving window in mass and relative age^{6}

Table 1 reports the position of the envelope as a function of the mass of the two stars, of their relative age, and of the mass ratio of the system^{7}

Several resulting features may depend in principle on the technique employed to obtain the system age and on the adopted perturbation strategy, so we devoted Sections 3.1 and 3.2 to discussing the influence of these possible bias sources. A detailed discussion of the sampling strategy adopted to build the synthetic dataset and an analysis of its effect on the obtained results can be found in Appendix A.

Mass () | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

0.8 | 0.9 | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 | 1.6 | |||

primary star | |||||||||||

-7.1 | -7.0 | -7.1 | -7.2 | -6.7 | -6.6 | -6.8 | -7.1 | -5.2 | |||

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |||

6.5 | 7.0 | 6.6 | 7.6 | 7.5 | 7.4 | 7.0 | 7.0 | 7.4 | |||

Secondary star | |||||||||||

-7.2 | -7.0 | -7.1 | -7.1 | -7.1 | -6.7 | -6.3 | -5.3 | -3.4 | |||

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.2 | |||

7.0 | 7.5 | 7.7 | 7.7 | 7.6 | 7.4 | 6.5 | 6.0 | 5.3 | |||

Relative age | |||||||||||

0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | |

Primary star | |||||||||||

-23.5 | -30.8 | -22.9 | -15.5 | -11.3 | -8.9 | -7.5 | -6.3 | -5.8 | -5.0 | -4.3 | |

5.0 | 0.3 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

92.7 | 46.2 | 24.6 | 16.2 | 11.7 | 9.1 | 6.9 | 6.4 | 6.1 | 5.0 | 4.1 | |

Secondary star | |||||||||||

-23.3 | -17.0 | -9.8 | -7.7 | -6.8 | -6.3 | -5.8 | -5.2 | -4.8 | -4.3 | -3.9 | |

0.2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

47.3 | 19.3 | 10.3 | 7.9 | 6.5 | 6.3 | 5.8 | 5.3 | 4.9 | 4.3 | 3.9 | |

Mass ratio | |||||||||||

0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95 | 1.00 | |

-3.6 | -4.7 | -5.9 | -6.6 | -7.2 | -7.6 | -7.9 | -7.9 | -7.0 | -6.4 | -6.0 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

5.0 | 5.8 | 6.3 | 7.6 | 8.3 | 8.6 | 8.5 | 8.4 | 7.5 | 6.7 | 6.3 | |

Mass ratio (alternative sampling) | |||||||||||

-5.2 | -5.9 | -6.7 | -7.5 | -8.4 | -9.1 | -9.5 | -9.6 | -8.8 | -8.6 | -8.5 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

4.9 | 6.2 | 7.1 | 8.3 | 9.0 | 9.6 | 10.1 | 10.4 | 9.5 | 9.2 | 8.9 | |

Mass ratio (rejection step) | |||||||||||

-3.6 | -4.7 | -5.9 | -6.6 | -7.2 | -7.6 | -7.9 | -8.2 | -8.0 | -7.5 | -7.0 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

5.0 | 5.8 | 6.3 | 7.6 | 8.3 | 8.6 | 8.5 | 8.6 | 8.6 | 8.2 | 7.9 |

### 3.1 Effect of correlations

To simulate the observational uncertainties (see Sect. 2), we added a Gaussian perturbation to the four observables (i.e. , , , [Fe/H]) of the stars in the previously generated synthetic dataset of artificial binary systems. In the standard scenario they are sampled from a multivariate normal distribution assuming correlation of between the two effective temperatures, between the two metallicities, between the two masses, and zero between the two radii. These correlations arise from the fact that the analysis of data for binary stars usually determines some function of the observables at higher accuracy than the observables themselves.

The strongest correlation is expected between the effective temperatures of the two stars, since the difference in or the ratio of the effective temperatures in an eclipsing system is usually determined with relative accuracy that is about three times greater than that of the individual effective temperatures themselves (see e.g. Claret2003; Southworth2007; Southworth2013; Torres2014). A direct computation, by means of error propagation under normality assumption, on the data presented in the cited literature showed that the correlation of the two effective temperatures is reproduced well when assuming .

Regarding [Fe/H], we adopted a high correlation coefficient () since usually only a mean value for the entire system is presented in the analyses. For mass determination, the mass ratio is usually determined with relative accuracy of about 1.6 higher than the mean one of the masses themselves (see e.g. Popper1986; Lacy2008; Torres2009; Brogaard2011; Vos2012; Sandquist2013). Assuming normality, this corresponds to a correlation of between the two star masses. However, even significantly lower values for the correlation () resulted from literature examination (Helminiak2009; Meibom2009).

The situation is much more complex in the case of radii determination. The light curve analysis poses independent constraints on the sum of radii and on their ratio , causing a dependence between the two evaluated stellar radii; however, the magnitude and the sign of the correlation between the estimated radii depend on the actual values of and and of their errors. Direct computations from the values quoted in the literature showed correlations varying from to (Popper1986; Grundahl2008; Torres2009; Brogaard2011; Vos2012). Luckily, such high variability is not a serious problem because, as shown below, the actual magnitude of the correlation between radii has very little impact on the results. Thus, for the reference scenario we assumed uncorrelated radii.

To quantify the impact on the age estimate uncertainty of accounting for the aforementioned correlations in the procedure followed to build the synthetic datasets, we applied the SCEPtER pipeline on datasets produced by relying on different assumptions about correlations. As a first test, we provided a dataset of artificial binary stars with uncorrelated observables; i.e., the Gaussian noise was added assuming a diagonal covariance matrix. Although such a choice is quite extreme and not realistic, the comparison between the outcomes of the age estimate procedure on this dataset with those on the standard one will quantify the maximum effect of the correlations. Figure 2 shows the impact of the correlation between observables on the relative errors in the estimated age of the system. The neglect of correlation on the corresponding observables of the two stars causes a mean relative shrinkage of the envelope of about 20%. Direct simulations (not shown) proved that the shrinkage is due to dropping the correlation between the two masses, while the correlations between the effective temperatures and between the metallicities account for modest variations.

As a second test, we built two synthetic datasets by assuming the same correlations of the standard case, with the exception of the correlation coefficient between the radii of the binary components, namely . The comparison of the performances of the recovery procedure on these datasets with those on the standard one will assess the impact of the correlation amongst the radii. Table 2 reports, as a function of the mass of the primary star, the results obtained accounting for correlations between the radii. The mean variations on the envelope boundaries are about 0.2%, so it is safe to neglect the correlation among radii.

Primary star mass () | |||||||||
---|---|---|---|---|---|---|---|---|---|

0.8 | 0.9 | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 | 1.6 | |

-7.1 | -6.9 | -7.3 | -7.4 | -6.9 | -7.2 | -7.0 | -6.8 | -5.2 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

6.9 | 6.6 | 6.6 | 7.3 | 7.0 | 7.1 | 6.9 | 7.3 | 7.4 | |

-6.4 | -6.6 | -7.0 | -7.2 | -6.8 | -6.8 | -6.8 | -6.9 | -5.3 | |

0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |

7.1 | 7.4 | 6.7 | 7.4 | 7.5 | 7.0 | 7.0 | 7.1 | 7.0 |

### 3.2 Recovery algorithms

Figure 2 also shows the impact of the adopted age-recovery technique. The simple average of the individually independent age estimates of the two stars appears to be as good as the joint estimate in the case of a primary star of low mass ( ), while it shows a notable larger variance in binary systems composed of a high-mass primary and a low-mass secondary. For these combinations of objects, the secondary stars is sampled in the first stages of the evolution, since the maximum time is set by the faster evolution of the more massive star. The age estimate of a star in the early evolutionary stages is subject to larger relative uncertainty, and this uncertainty is propagated in the final mean. Such an effect is not present in the joint evaluation, since in this case the age is mainly determined by the more evolved star. In fact the age spread of the models in the box around the primary star is generally lower than the spread the secondary star, since the primary is sampled at a higher relative age than the secondary. Therefore the joint likelihood estimate is limited by the range of age spanned by the primary star.

The increase in the variance in age estimates when coevality is not explicitly assumed is shown further in Fig. 3. The top
row of the figure shows, as
a function of the masses and of the binary system components, the
position of the 2D lower (left panel) and upper (right panel)
envelopes^{8}

### 3.3 Contribution of the secondary star to age estimates

In Sect. 3.2 we discussed the dominant impact of the primary star on joint likelihood age determination. It is therefore interesting to explore the contribution of the secondary star to the whole system age determination. We performed this test by comparing the envelopes for the joint likelihood age estimate and for the primary single-star estimates. The results are presented in Fig. 4. Computed as a function of the primary mass, the joint estimate’s relative error envelope shows negligible differences for both the upper and the lower boundaries with respect to the envelope computed only for primary stars. The contribution of the secondary star is only visible for young systems, where the age of the primary star is less constrained. In fact, the right-hand panel of the figure shows that the relative age envelope that takes the secondary contribution into account is less prone to age underestimation at a relative age lower than 0.15. In summary, the secondary star plays a minor role with respect to the primary one because its age is generally estimated with a larger uncertainty being closer to the ZAMS than the more massive companion of the same age.

## 4 Stellar model uncertainty propagation

The accuracy of age estimates of real binary systems obviously depends on the reliability of the adopted grid of stellar models. In this section, following V14 and V15, we provide an analysis of the impact of the various uncertainty sources that affect the stellar computations. We focus here on some issues considered in V15: the convective core overshooting efficiency, the initial helium content, and the neglect of the microscopic diffusion. We constructed several artificial grids characterised by the variation in the single input under examination up to its extreme allowed values, while keeping all the others fixed to their reference values. We then built the corresponding synthetic datasets by sampling artificial binary systems from these non-standard grids of stellar models and by adding a Gaussian noise in all the observed quantities to simulate observational uncertainties following the same prescriptions for covariance matrix described in the standard case. Finally, the ages are estimated by means of the SCEPtER pipeline by relying on the standard grid of stellar models.

Although the use of a solar-calibrated mixing-length value for stars that differ from the Sun for mass, composition, or evolutionary phase could be inappropriate in the present paper; unlike V15, we do not consider this source of uncertainty. Such a choice is motivated by the fact that the dependence of the mixing-length value on the mass, metallicity, level of activity, and the evolutionary phase of the star is still under study and that no definitive conclusion has been reached (see, among many, Clausen2009; Deheuvels2011; Bonaca2012; Mathur2012; Tanner2014). In V15 we studied the impact on age estimates due to a change in the mixing-length value by adopting synthetic grids with non-solar mixing-length values, which did not vary inside a given grid regardless of the mass and metallicity of the star. Since in the present work, we consider systems of stars with different masses, hence in different evolutionary stages, we feel that the approach of V15 might not be appropriate. Besides the effect of a change in the mixing length of the sampling grid, a differential effect of the possible evolution of the mixing-length value with the mass of the stars could be important. Nevertheless, we performed a numerical experiment, by sampling from grids with with respect to our solar-calibrated value ( = 1.74), which showed a bias ranging from about to , the lowest values for mass ratio near 1.0. These simulations suggest that the bias due to the different choices of the mixing-length values is certainly important but not dominant. For a more rigorous investigation, better knowledge of the dependence of the external convection efficiency on stellar parameters, such as mass, chemical composition, evolutionary phase, and stellar activity, is mandatory.

The observables adopted here in the age reconstruction are different from those in V15, which made use of classical (i.e. the effective temperature of stars and their metallicity [Fe/H]) and asteroseismic (i.e. the average large frequency spacing and the frequency of maximum oscillation power ) observable constraints. We therefore expect that the impact of the uncertainty sources is different from what is found in our previous work. In particular, the availability of precise stellar masses is paramount since it severely restricts the number of stellar models that can actually contribute to the recovery procedure.

Mass () | |||||||||

0.8 | 0.9 | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 | 1.6 | |

Overshooting | |||||||||

0.0 | -0.7 | -3.6 | -6.2 | -5.7 | -6.1 | -4.9 | |||

-0.2 | -0.2 | -1.0 | -3.0 | -4.7 | -6.6 | -6.8 | -6.6 | -6.1 | |

Overshooting | |||||||||

0.0 | -3.3 | -10.3 | -10.7 | -10.2 | -10.7 | -10.5 | |||

-0.8 | -1.6 | -4.0 | -8.2 | -11.9 | -10.9 | -11.2 | -11.2 | -9.4 | |

-6.4 | -6.8 | -9.7 | -10.1 | -10.7 | -10.0 | -9.8 | -9.2 | -7.8 | |

-7.9 | -8.3 | -10.0 | -10.2 | -9.9 | -9.7 | -9.2 | -8.1 | -6.8 | |

6.9 | 10.1 | 10.5 | 11.4 | 11.4 | 10.5 | 10.1 | 9.7 | 10.2 | |

9.4 | 10.5 | 11.6 | 10.6 | 10.5 | 10.4 | 9.7 | 8.7 | 8.6 | |

No microscopic diffusion | |||||||||

8.0 | 6.2 | 5.1 | 4.0 | 3.8 | 3.3 | 2.8 | 2.6 | 2.9 | |

6.4 | 4.9 | 3.9 | 3.4 | 3.3 | 3.1 | 2.6 | 2.0 | 1.7 |

Mass ratio | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

Input | 0.50 | 0.55 | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95 | 1.00 |

-3.7 | -3.6 | -3.7 | -3.8 | -3.9 | -3.8 | -3.5 | -3.2 | -2.6 | -3.4 | -3.4 | |

-8.2 | -9.2 | -9.7 | -9.8 | -9.5 | -9.0 | -8.1 | -7.3 | -6.1 | -5.0 | -4.6 | |

-7.6 | -8.0 | -9.0 | -9.3 | -9.7 | -9.8 | -9.9 | -9.8 | -8.9 | -9.1 | -8.9 | |

10.0 | 10.5 | 11.2 | 11.3 | 11.1 | 11.3 | 11.3 | 11.2 | 10.2 | 9.9 | 9.7 | |

No diffusion | 4.0 | 3.8 | 3.4 | 3.3 | 3.4 | 3.7 | 3.9 | 4.1 | 3.8 | 3.4 | 3.2 |

### 4.1 Convective core overshooting

As in V15 we parametrised the extension of the extra-mixing region beyond the canonical border, as defined by the Schwarzschild criterion, in terms of the pressure scale height : , where is a free parameter. To quantify the impact of taking the convective core overshooting into account, we computed – only for models more massive than 1.1 – two additional grids with values of = 0.2 and 0.4, the latter representing a possible maximum value (see e.g. the discussion in cefeidi). We then extended these two grids of stellar models at masses lower than 1.1 with the standard one, which neglect core overshooting.

From these non-standard grids we built two synthetic data sets (for and ) by sampling artificial binary systems each and adding a Gaussian noise. The age estimate is then performed by adopting our standard grid of models (without overshooting) in the recovery procedure. The results are shown in Tables 3 and 4 and in Figs. 5 and 6. To evaluate the magnitude of the obtained biases, they can be compared with the standard envelope due to observational uncertainties in Table 1. It turns out that the median bias due to the mild (i.e. ) core overshooting scenario – computed as function of the mass ratio – is about 50% of the half width of the standard envelope uncertainties. The corresponding bias for the strong () core overshooting scenario is about 120% of the standard envelope half width. In both cases the relative importance of the bias is generally greater for low values of ; for the strong overshooting scenario, in particular, the bias reaches values of about 160% of the standard envelope half width at = 0.5, while it is only 60% at = 1.0. This trend is expected since part of the systems with near 1.0 are composed of stars with both masses lower than 1.1 , for which core overshooting is not taken into account so that the discussed effect is partially masked.

Figure 6 shows that the dominant effect for the bias is due to the mass of the primary star. As an example, the mean bias for mild overshooting scenario (left panel in the figure) is for a primary star that is less massive than 1.1 , while it is for a massive primary coupled to a light secondary star. The shift due to having a massive secondary star is , so it only accounts for an additional change. The contribution of the secondary star is even less important in the strong overshooting scenario (right panel in the figure), when it accounts for about one-fifth of the primary star one.

The results presented in this section have particular importance since binary systems are often adopted to calibrate the efficiency of the convective core overshooting by isochrone fine tuning (see e.g. Claret2007; Lacy2008; Clausen2010). Although our results do not directly allow estimation of the errors in the recovered best convective core-overshooting parameter that is consistent with the observations, we found that the bias induced by a mild convective core overshooting scenario on age estimate is about one-half of the expected error in age estimate. This implies a difficulty unambiguously identifying the overshooting effect, which can be masked by random fluctuations.

The present results suggest that a mild overshooting scenario is hardly distinguishable from a no overshooting one. This result agrees with the finding by Claret2007 in the lower mass range considered in that paper. Although a statistical calibration – obtained from several systems – can be meaningful, the large errors that propagate into the final calibrated overshooting parameter shed some doubt on the calibration of single systems. Owing to the widespread use of binary systems for calibration purposes, we feel that further research is needed to statistically quantify the errors in the calibrated parameters.

### 4.2 Initial helium abundance

We quantified the impact on the age estimate of the current uncertainty in the initial helium abundance following V15. We computed two additional grids of stellar models with the same metallicity values as in the standard grid, but by changing the helium-to-metal enrichment ratio to values of 1 and 3. Then, we built two synthetic datasets, each of artificial binary systems, by sampling the objects from these two non-standard grids and adding Gaussian noise. The age of the objects is reconstructed using the SCEPtER pipeline relying on the standard grid with .

The results are presented in Tables 3 and 4 and in Figs. 7 and 8. The mean impact of the initial helium abundance uncertainty is higher than for the convective core overshooting. The median bias for – computed as a function of the mass ratio (Table 4) – is about the of the half-width of the standard envelope due to observational uncertainties reported in Table 1. The corresponding bias for is about 160% of the standard envelope half width in Table 1. Figure 8 shows that both the changes in the bias due to the primary and secondary stars are modest, so, in the examined mass range, the initial helium induced bias can be considered almost constant at about for a change in by .

The results presented in this section on the initial helium uncertainty impact are very different from those reported in V15, where this uncertainty source have been found to be nearly negligible. This striking difference results from the observational constraints adopted in the estimation. The analysis in V15 in fact showed that the little influence of the helium content on age estimates came from two opposite effects that nearly cancel each other. The first obvious effect was due to the impact of the helium content on the evolutionary time scale, so that helium-rich models evolve faster. A second effect was due to the placement on the standard recovery grid of the non-standard helium models; helium-rich models mimicked more massive standard models, leading to a mass overestimation bias and therefore to reduced age estimates. The two effects compensate for each other.

This balancing effect cannot occur here since the mass is adopted as an observational constraint, and this leads to the strong influence of the initial helium content on the final age estimates.

### 4.3 Element microscopic diffusion

We evaluated the impact on the age estimates of neglecting microscopic diffusion in the stellar models adopted in the recovery procedure following, as in V15, a different approach with respect to the previous cases. Here we preferred to build the synthetic data set by sampling artificial binary systems from stellar models that take diffusion into account and by adding Gaussian noise. We then reconstructed their age by means of the SCEPtER pipeline but adopting the non-standard grid of models, which neglect diffusion.

The results are presented in Tables 3 and 4 and in Figs. 9 and 10. The bias due to the microscopic diffusion has a mean value of 3.6% (see Table 4), it is similar in absolute value to the one due to mild overshooting, and it ranges from about one-half to one-third of the bias due to initial helium change. Figure 10 confirms an expected result that the diffusion is most important in systems composed of two low-mass stars. For primary and secondary masses lower than 1.1 , the bias in age estimates is 5.6%, while it drops to 3.7% for more massive primary stars. The trend can be understood since the evolutionary time scale of these stars is faster than for the diffusion processes, which are therefore quite inefficient (see left panel in Fig. 9).

As already discussed in Sect. 4.2, this result confirms that the bias in age grid estimates, owing to the different input in the evolutionary stellar codes, and strongly depends on the assumed observational constraints. In fact V15 has shown that – if adopting asteroseismic constraints – the microscopic diffusion-induced bias dominated the other examined biases, while it is not the case here.

## 5 The coevality problem: a statistical approach

In the literature there are several cases of real eclipsing binary systems whose components are estimated to be non-coeval, because they cannot be fitted simultaneously by a single isochrone (e.g. Clausen2009; Vos2012; Sandquist2013). Such an occurrence is often used to claim that the current generation of stellar models present some weaknesses and that additional physical processes should be included or re-calibrated. Given its considerable implications, it is worth discussing whether the estimates of non-coevality amongst the two components of real binary systems are statistically reliable or are simply a fluctuation. In this regard it is of particular interest to quantify the expected difference in the estimated ages of two genuine coeval stars caused simply by the observational uncertainties. To do that, we analysed our previously described synthetic datasets of coeval binary stars.

As detailed in Sect. 2, the artificial binary systems sampled from the grid of stellar models are subject to random perturbations to account for the observational errors. Therefore the ages of the single stellar members, which are equal by construction in the sampling stage, are generally estimated as different. It is interesting to estimate how large the expected differences are in age owing only to these random perturbations, since this effect also interests real world observations.

A detailed analysis of the problem has to deal with many technical aspects^{9}

We define and as the estimated ages of the two members, with . We focused our attention on the statistics , which then ranges from 0, when the stars are correctly estimated coeval, and 1, when one star is estimated to be much younger than the other. To develop a statistical test based on we have to estimate how large can become only because of the random fluctuations in age estimates and . Once the distribution of – by means of a Monte Carlo simulation has been empirically estimated, it is possible to choose a critical value by identifying the range of values of that are too extreme to be compatible with the coevality hypothesis. Usually the critical value is chosen as the 95th or the 99th quantile of the distribution of the statistics under consideration. The choice of the critical value to adopt (in the following, we refer to this as ) defines the ”level” of the test ( thus being the probability that an observed value of is larger than ). The set of values higher than the critical values are said to lie in the ”rejection region”.

We present the rejection regions at level = 0.05, hence the critical values identifying them, of the null hypothesis that the stars are coeval. Reconstructed values of higher than the critical value lead to the rejection of the null hypothesis, implying that standard stellar models cannot account for the coevality of the stars with the assumed input or parameters.

A first obvious complication comes from the value of depending on the mass of the two stars, on their metallicity, and on the (observationally unknown) relative age. Here we focus on the impact of the binary star masses.

To perform the critical-value estimation, we generated a sample of artificial binary systems, perturbed their observables as explained in Sect. 2, and then independently estimated the ages of the two stars. From these estimates, we computed the statistics . Then we binned the values according to the primary and secondary stellar masses, as for constructing the 2D envelope in Sect. 3.1. In each of these bins, we computed the empirical 95th quantile, which approximates the required critical value. A detailed discussion of the dependence of such a critical value on the actual uncertainties affecting the observational constraints and on the evolutionary phase of the two binary members will be presented in a forthcoming paper (Valle et al. 2015, in prep.).

The results of the analysis are displayed in Fig. 11. The figure shows that genuine coeval stars can be reconstructed as non-coeval with a sizeable relative age difference only because of the current uncertainty in the observational constraints. While the critical value is about 0.25 for systems of nearly equal masses, it can be greater than 0.6 for unbalanced binary masses.

This result should be carefully considered when determining the stellar ages of a binary system. Thus, a determination of two discrepant stellar ages, which leads to a value of , does not allow a statistically grounded rejection of the hypothesis that the stars are coeval since the difference might only be due to a random fluctuation in the observational constraints.

As discussed in Sect. 4.1, one should be aware of this effect whenever the extent of convective core overshooting is calibrated on eclipsing binary by requiring that both stellar components are simultaneously fitted by a single isochrone, i.e. by imposing their coevality. If the ages estimated by means of stellar models without core overshooting give a value of , one cannot conclude on statistical grounds that these models are not able to simultaneously fit both stars and that a convective core overshooting must be taken into account. Therefore, the reliability of the calibration method of the convective core overshooting consisting in the isochrones fine tuning to the observation has to be cautiously evaluated.

In contrast, when standard stellar models provide ages of the binary components such that , the possibility that non-coevality is the result of observational uncertainties should be considered as negligible. As an example, binary systems in which at least one star is suspected of supporting a strong surface activity, generally show a large discrepancy among the ages of the two components, such as V636 Cent and EF Aqr (Clausen2009; Vos2012). In this system, the active secondary stars are estimated to be much older than the primary, far above the limit due to random fluctuations. The coevality is recovered adopting a significantly lower mixing-length value for the secondary stars. This is an expected behaviour since it is recognised that in active stars the dynamo-induced magnetic fields can suppress convection and produce starspots, causing differences in temperature and radii with respect to standard stellar models, which can be mimicked by a lower mixing-length value (Gabriel1969; Cox1981; Clausen1999). It is, however, clear that further research is mandatory to determine the expected errors in the calibrated mixing length in a sound statistical way.

## 6 Conclusions

We performed a detailed theoretical investigation of the age estimation of detached double-lined eclipsing binaries by means of a maximum likelihood grid-based technique. We analysed the impact of the current uncertainty on several input of stellar model computations, such as the chemical composition, the efficiency of convective core overshooting, and the efficiency of microscopic processes.

We adopted the grid-based pipeline SCEPtER, which is extensively described in scepter1 and eta. As observational constraints we assumed the stellar effective temperature, the metallicity [Fe/H], the mass, and the radius of the stars. The grid of stellar models, computed for the evolutionary phases from ZAMS to the central hydrogen depletion covers the mass range [0.8; 1.6] .

We computed the statistical errors arising from the uncertainties in observational quantities. We found an overall relative uncertainty of about in age estimates of the system, which is nearly independent of the masses of the stars and of the mass ratio of the binary systems. However, as described in eta, the relative error in the age estimates sensitively depends on the evolutionary phase and becomes larger for models near those ZAMS for which it can reach values of about 90% for the upper boundary, with an envelope half width of about 50%. The large uncertainty and the bias toward higher ages found near the ZAMS are mainly an edge effect distortion, since the estimated age cannot be negative. The relative error envelope quickly shrinks as the stars evolve from the ZAMS, and for primary stars of relative age 0.3, the half width of the error envelope is about 15%.

We studied the impact on the age estimates of different choices of the sampling algorithms, of the correlation between the effective temperatures, the masses, the radii, and [Fe/H] of the two binary components, and of the method of estimating the joint age of the stars. It resulted that the first source of uncertainty plays a minor role since it accounts for a variation in the relative error half width of about 1%. The neglect of correlations between the observables of the two stars was shown to reduce the age error envelope width by about one-fifth. The age estimates of the binary system obtained by averaging the ages of the two stars computed independently – without explicitly assuming coevality – lead to negligible differences only for stars of nearly equal masses. Whenever an unbalanced system is considered ( and ), the envelope width of the age relative errors for the average of independent age estimates is larger by a factor of two, but it can be as high as a factor of five for systems with .

We evaluated the systematic biases on age estimates owing to the uncertainties in initial helium content, in the convective core overshooting, and in the microscopic diffusion efficiency. We found that the bias due to a mild ( = 0.2) and strong ( = 0.4) convective core overshooting scenarios are about 50% to 120% of the half width of the standard age relative error envelope due only to the observational uncertainties. The bias due to the initial helium content was explored assuming an uncertainty of on the helium-to-metal enrichment ratio . The effect is about and 160% of the random error, so slightly more than that of the strong overshooting scenario. Finally, the bias due to the neglect of the element microscopic diffusion was about 60% of the half width of the standard age-relative error envelope.

The comparison of these results to those reported in eta, which were computed by adopting asteroseismic constraints and without the knowledge of stellar mass and radius, showed relevant differences. In fact, eta showed that the helium-induced bias on age estimates was negligible, while the bias due to the microscopic diffusion was dominant. These differences arise from the corresponding bias on mass estimates that could occur in eta, which led to some compensations for the age estimates. These mechanisms cannot occur in the present study, since the stellar mass is assumed as an observational constraint. An important point is therefore that the bias in the grid-based age estimates – due to adopting different input in a stellar evolutionary code – strongly depends on the assumed observational constraints and should be evaluated case-by-case whenever a new grid is adopted for estimates.

We also introduced a statistical test on the expected differences in the estimated ages of two genuine coeval binary components due simply to the observational uncertainties affecting the observables used in the recovery procedure. We found that the relative difference between the reconstructed ages of two coeval binary members depends on the mass ratio . For systems with , the relative difference in the inferred ages can be about 20%, while it is higher than 60% for system with . This result should be carefully taken into account before using the apparent non-coevality of stars in binary systems – i.e. the inability to find a single isochrone that fits both components – to claim the weakness of current generation of stellar models and to calibrate the efficiency of some poorly known physical process, such as convective core overshooting, by isochrone fine tuning with the observations.

###### Acknowledgements.

We thank our anonymous referee for the stimulating comments that helped us to clarify and improve the paper. This work has been supported by PRIN-MIUR 2010-2011 (Chemical and dynamical evolution of the Milky Way and Local Group galaxies, PI F. Matteucci), and PRIN-INAF 2012 (The M4 Core Project with Hubble Space Telescope, PI L. Bedin).## Appendix A Synthetic dataset: sampling strategy

In constructing the sample adopted in the error estimation, we did not try to match any observed distribution of the mass of the binary stars or of the mass ratio . Since we primarily aimed to investigate the performances of grid-based age estimates of binary stars, we paid more attention when sampling the whole range in mass, metallicity, and evolutionary phases covered by the available grid of models.

In detail, the datasets of artificial stars are sampled by adopting the following scheme. First, a star is randomly sampled from the grid. No priors on mass, age, and metallicity are assumed in the sampling procedure. Then, a second star is coupled to the first one with two requirements: it must share with the first star both the same initial [Fe/H] and the age, with a tolerance of 10 Myr. The couple of stars is re-ordered to have the most luminous star as a first member.

To show which models are selected as primary and secondary stars, we estimated
the joint density of mass and relative age for primary and secondary stars^{10}^{11}

The effect of the sampling could be seen in Table 1, which shows a small shrink in the age relative error envelope at greater than 0.9. This is a direct consequence of the relative error envelope narrowing at high relative age. To show that this artificial trend could be removed with different sampling strategies, in Table 1 we also present two additional series of results. In fact, we first explored the impact of a different sampling. For this test, we sampled the set of primary stars, then we coupled each of those objects to the secondary star by sampling only from the models respecting the constraint in metallicity and age, which are also less massive than the primary object. The results are in the rows labelled “mass ratio (alternative sampling)” in Table 1. As a second test, to see the impact of the ad-hoc removal of the trend of primary relative ages versus shown in Fig. 13, we checked an alternative approach that rejects some objects in the subset , based on the standard method. In fact, since this is the region that is more populated, it was possible to draw from it a sample of size equal to that of the adjacent subset , with the requirement that the new sample’s primary relative age matches the distribution of the subset . The corresponding results are in the rows ”mass ratio (rejection step)” in Table 1. It appears that the shrinkage in the upper region is reduced in both two cases.

In contrast, the shrinkage of the error envelope at low mass ratio (see Table 1) is due to an edge effect since there are only models in this region with mass 1.6 coupled with 0.8 models. In fact the lower left-hand region of Fig. 13 is not populated because to produce a coeval combination, the more massive model cannot be too young. In fact, in this case the secondary is still in the pre-main sequence phase, which is not considered in the grid. Figure 13 shows that the minimum relative age for these massive models is about 0.2. Therefore the most critical region for age estimation – that of low relative age – is avoided, and the estimate relative error envelope is narrower.

## References

### Footnotes

- offprints: G. Valle, valle@df.unipi.it
- http://CRAN.R-project.org/package=SCEPtERbinary
- The chosen value allows a mean relative accuracy of 0.1% to be reached on the confidence interval.
- http://astro.df.unipi.it/stellar-models/
- The ZAMS and the central hydrogen depletion models are defined as the models whose central hydrogen abundance drops below 99% of the initial value and , respectively.
- The half width of the windows is chosen to maintain the error on the mass envelope due to Monte Carlo sampling at a level of about 0.2%, without introducing too much smoothing. The corresponding error on the envelope – relative to 2.5th and 97.5th quantiles – is about 1%.
- As usual, is defined as the ratio between the masses of the secondary and the primary stars.
- The bidimensional envelopes are computed with a similar technique to the 1D ones, with mass steps of 0.1 and moving-window half width of 0.08 . The mass binning choice is such that the envelope is affected by a typical random uncertainty of about 0.5%.
- Such as the choice of the sampling scheme for the binary systems, the evaluation of the precision and the accuracy on the estimated critical values, and of the sample size needed to attain a given level of precision.
- The joint density was obtained by a bi-dimensional kernel density estimation using the function kde2d in the R library MASS. Details on the technique can be found in venables2002modern.
- A LOESS (LOcal regrESSion) smoother is a non-parametric, locally weighted polynomial regression technique that is often used to show the underlying trend of scattered data (see e.g. Feigelson2012; venables2002modern).