© ESO, 2015

1. Introduction

Studies of extrasolar planets rely on a good understanding of the stars that they orbit. To estimate the mass and radius of an extrasolar planet that transits its host star we require an estimate for the mass of the star. The mass of the star will also strongly influence the planet’s environment, e.g., the spectrum and intensity of the stellar flux intercepted by the planet and the nature and strength of the tidal interaction between the star and the planet. The ages of planet host stars are used to investigate the lifetimes of planets and the time scales for tidal interactions between the planet and the star (e.g., Matsumura et al. 2010; Lanza 2010; Brown et al. 2011).

The analysis of the light curve produced by a planetary transit yields an accurate estimate for the radius of the star relative to the semi-major axis of the planet’s orbit, R_⋆/a, provided that the eccentricity of the orbit is known. This estimate can be combined with Kepler’s laws to estimate the density of the host star. The density can be combined with estimates for the effective temperature and metallicity of the star to infer a mass and age for the star by comparison with stellar models or an empirical calibration of stellar mass. In general, the comparison with stellar models is done using a maximum likelihood method, i.e., taking the mass and age of the evolution track and isochrone that give the best fit to the observed density (or R_⋆/a) and effective temperature, estimated either by a least-squares fit to the observed properties or “by-eye”.

Maximum-likelihood estimates can be strongly biased in cases where the mapping between the observed parameters and the parameters of interest is non-linear. This is certainly the case for stellar ages because the observed parameters change very little during the main-sequence phase, but there are large changes to the observed properties during the rapid evolution of a star away from the main sequence. This can produce a “terminal age bias”, where the “best-fit” method applied to a sample of stars produces a distribution of ages that is a priori very unlikely, i.e., too few main-sequence stars and many stars in regions of the model paramater space corresponding to short-lived evolutionary phases. This problem is particularly acute for cases where the uncertainties on the observed parameters are large compared to the change in observed properties during the main sequence phase. This is often the case for the age estimates of single stars based mainly on surface gravity measured from the stellar spectrum or the absolute magnitude derived from the measured parallax and stellar flux. Bayesian methods that account for the a priori distribution of stellar ages have been developed to deal with this problem in spectroscopic stellar surveys (Jørgensen & Lindegren 2005; Pont & Eyer 2004; Serenelli et al. 2013; Schneider et al. 2014). Pont et al. (2009) applied a similar Bayesian approach in their study of the HD 80606 planetary system, but there have been few other examples of this approach in exoplanet studies. This may be because there is currently no software available to the exoplanet research community that can be used to apply these Bayesian methods and that is straightforward to use.

In general, the terminal age bias is expected to be less severe for planet host stars than for single stars because stellar density measurements based on the analysis of a planetary transit are usually more precise than surface gravity estimates based on the analysis of a stellar spectrum or the luminosity derived from parallax and flux measurements. The mean stellar density is also more sensitive to the change in radius of a star as it evolves away from the main sequence. However, precise stellar densities can also cause problems because the broad sampling in mass, age and metallicity used for many grids of stellar models can produce poor sampling of the observed parameter spacing, i.e., the typical difference in stellar density between adjacent model grid points can be much larger than the uncertainty on the observed value. This can produce systematic errors due to interpolation, particularly for stars near the end of their main-sequence evolution where the evolution tracks have a complex behaviour that is very sensitive to age, mass and composition. This also makes it difficult to make reliable estimates of the uncertainties on the mass and age. One method sometimes employed to estimate the uncertainties is to look for the mass and age range of all the models that pass within the estimated errors on the observed values, e.g., all the models within the 1σ error bars. This approach can be misleading because the errors on the mass and age are often strongly non-Gaussian and highly correlated . It also possible to miss some combinations of mass, age and composition that provide a reasonable match to the observed properties of the star but that are not sampled by the stellar model grid or that fall just outside the 1σ error bars, particularly when the fitting is done by-eye.

The recently-launched European Space Agency mission Gaia (Perryman et al. 2001) will measure parallaxes and optical fluxes for many stars that are transited by planets, brown-dwarfs and low-mass stars. The luminosity measurement derived from these observations can be an additional useful constraint on the mass and age of the star. If the density, effective temperature and composition of the star are also known then finding the best-fit mass and age of the star becomes an over-determined problem for stellar models with fixed helium abundance and mixing-length parameter, i.e., there are more observables than unknowns. If the best-fit to the observed parameters is poor then this opens up the possibility of using these stars to explore whether the assumed values of the helium abundance and mixing-length parameter, or some other factor, can be adjusted to improve the agreement between the stellar models and real stars.

To deal with these issues we have developed a Markov-chain Monte Carlo (MCMC) method that calculates the posterior probability distribution for the mass and age of a star from its observed mean density and other observable quantities using a grid of stellar models that densely samples the relevant parameter space. We have validated our method by applying it to data derived from different stellar models and by applying it to stars in eclipsing binary stars with precisely measured masses and radii. We have also quantified the systematic error in the estimated mass and age due to the variations in the assumed helium abundance and convective mixing length parameter. The method has been implemented as a program called bagemass that we have made available for general use.

2. Method

2.1. Stellar models and grid interpolation

Our method uses a grid of models for single stars produced with the garstec stellar evolution code (Weiss & Schlattl 2008). The methods used to calculate the stellar model grid are described in Serenelli et al. (2013) so we only summarise the main features of the models and some differences to that description here.

garstec uses the FreeEOS¹ equation of state (Cassisi et al. 2003) and standard mixing length theory for convection (Kippenhahn & Weigert 1990). The mixing length parameter used to calculate the model grid is α_MLT = 1.78. With this value of α_MLTgarstec reproduces the observed properties of the present day Sun assuming that the composition is that given by Grevesse & Sauval (1998), the overall initial solar metallicity is Z_⊙ = 0.01826, and the initial solar helium abundance is Y_⊙ = 0.26646. These are slightly different to the value in Serenelli et al. (2013) because we have included additional mixing below the convective zone in order reduce the effect of gravitational settling and so to better match the properties of metal-poor stars. Due to the effects of microscopic diffusion, the initial solar composition corresponds to an initial iron abundance [ Fe/H ] _i = +0.06.

The stellar model grid covers the mass range 0.6 M_⊙ to 2.0 M_⊙ in steps of 0.02 M_⊙. The grid of initial metallicity values covers the range [ Fe/H ] _i = −0.75 to −0.05 in steps of 0.1 dex and the range [ Fe/H ] _i = −0.05 to +0.55 in 0.05 dex steps. The initial composition of the models is computed assuming a cosmic helium-to-metal enrichment ΔY/ ΔZ = (Y_⊙ − Y_BBN) /Z_⊙, where Y_BBN = 0.2485 is the primordial helium abundance due to big-bang nucleosynthesis (Steigman 2010), so ΔY/ ΔZ = 0.984.

The initial abundance of all elements is scaled according to the value of [ Fe/H ] _i. For each value of initial mass and [ Fe/H ] _i we extracted the output from garstec at 999 ages from the end of the pre-main-sequence phase up to an age of 17.5 Gyr or a maximum radius of 3 R_⊙, whichever occurs first. We define the zero-age main sequence (ZAMS) to be the time at which the star reaches its minimum luminosity and measure all ages relative to this time.

To obtain the properties of a star from our model grid at arbitrary mass, [ Fe/H ] _i and age we use the pspline implementation of the cubic spline interpolation algorithm². We interpolated the stellar evolution tracks for stars that reach the terminal-age main sequence (TAMS) on to two grids, one from the start of the evolution track to the TAMS, and one that covers the post-main sequence evolution. The dividing line between the two is set by the age at which the central hydrogen abundance drops to 0. We use 999 grid points evenly distributed in age for each model grid. For stars on the main sequence we interpolate between models as a function of age in units of the main sequence lifetime at the specified values of mass and [ Fe/H ] _i. The interpolating variable for the post-main sequence properties is the age since the TAMS measured in units of the total time covered by the model grid. Splitting the grid of stellar models in this way improves the accuracy of the interpolation near the terminal-age main sequence.

We used garstec to calculate some models for solar metallicity at masses half-way between those used for our method and then compared the interpolated values to those calculated by garstec. For the masses that we checked (1.01 M_⊙ and 1.29 M_⊙) we find that the maximum error in the density is 0.01 ρ_⊙, the maximum error in T_eff is 25 K and the maximum error in log (L/L_⊙) is 0.015. The worst agreement between the calculated and interpolated models occurs at very young ages (<0.01 Gyr) and, for the 1.29 M_⊙ star, near the “blue hook” as the star evolves off the main sequence. Away from these evoutionary phases the error in the interpolated values is at least an order-of-magnitude smaller that these “worst case” values.

2.2. Input data

It is much easier to calculate the probability distribution functions of a star’s mass and age if the observed quantities can be assumed to be independent. To enable us to make this assumption we define the data to be a vector of observed quantities d = (T_eff,L_⋆, [ Fe/H ] _s,ρ_⋆), where T_eff is the effective temperature of the star, L_⋆ its observed luminosity, [ Fe/H ] _s characterises the surface metal abundance and ρ_⋆ is the stellar density.

The density of stars that host a transiting extrasolar planet can be determined directly from the analysis of the light curve if the eccentricity is known. From Kepler’s third law it follows that $\begin{array}{ccc}{\mathit{\rho }}_{\mathit{\star }}\mathrm{=}\frac{\mathrm{3}{\mathit{M}}_{\mathit{\star }}}{\mathrm{4}\mathit{\pi }{{\mathit{R}}_{\mathit{\star }}}^{\mathrm{3}}}\mathrm{=}\frac{\mathrm{3}\mathit{\pi }}{\mathit{G}{\mathit{P}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{q}\mathrm{)}}{\left(\frac{\mathit{a}}{{\mathit{R}}_{\mathit{\star }}}\right)}^{\mathrm{3}}\mathit{,}& & \end{array}$where P and a are the period and semi-major axis of the Keplerian orbit and q = M_c/M_⋆ is the mass ratio for a companion with mass M_c. The quantity (a/R_⋆) can be determined with good precision from a high quality light curve alone if the orbit is known to be circular since it depends only on the depth, width and shape of the transit (Seager & Mallén-Ornelas 2003). For non-circular orbits the same argument applies but in this case the transit yields (d_t/R_⋆), where d_t is the separation of the stars during the transit, so the eccentricity and the longitude of periastron must be measured from the spectroscopic orbit or some other method so that the ratio d_t/a can be determined. The error in the mass ratio will give a negligible contribution to the uncertainty on ρ_⋆ for any star where the presence of an extrasolar planet has been confirmed by radial velocity observations. The same argument applies to brown-dwarf or low mass stellar companions, but some care is needed to make an accurate estimate of the mass ratio, e.g., via the mass function using an initial estimate of the stellar mass, and the additional uncertainty in ρ_⋆ should be accounted for.

Most of the stars currently known to host exoplanets do not have accurately measured trigonometrical parallaxes. For these stars we set log (L_⋆/L_⊙) = 0 ± 5 so that this parameter has a negligible influence on the results. For stars with measured parallaxes the observed luminosity is L_⋆ = 4πd²f_⊕ where f_⊕, the flux from the star at the top of the Earth’s atmosphere corrected for reddening, and d, the distance to the star measured from its trigonometrical parallax. The values of T_eff and [ Fe/H ] _s can be derived from the analysis of a high quality stellar spectrum. These may be less accurate than the quoted precision because of the approximations used in the stellar model atmospheres and the uncertainties in atomic data used in the analysis. For that reason we set lower limits of 80 K on the standard error for T_eff and 0.07 dex for the standard error on [ Fe/H ] _s(Bruntt et al. 2010). For stars similar to the Sun a change of 80 K in the assumed value of T_eff results in a change of about 0.02 dex in the value of [ Fe/H ] _s derived from the analysis of the spectrum (Doyle et al. 2013). This is at least a factor of 3 lower than the minimum standard error that we have assumed for [ Fe/H ] _s so we ignore this weak correlation between T_eff and [ Fe/H ] _s. For a few stars the effective temperature can be determined directly from f_⊕ and the directly measured angular diameter. The quoted uncertainties on these directly-measured T_eff values should not be used if they are much smaller than about 80 K because our method does not account for the strong correlation between T_eff and L_⋆ in these cases and some allowance must be made for the uncertainties in the effective temperature scale of the stellar models.

2.3. Calculation of the Bayesian mass and age estimates

The vector of model parameters that are used to predict the observed data is m = (τ_⋆,M_⋆, [ Fe/H ] _i), where τ_⋆, M_⋆ and [ Fe/H ] _i are the age, mass and initial metal abundance of the star, respectively. The observed surface metal abundance [ Fe/H ] _s differs from the initial metal abundance [ Fe/H ] _i because of diffusion and mixing processes in the star during its evolution.

We use the MCMC method to determine the probability distribution function p(m | d) ∝ ℒ(d | m)p(m). To estimate the likelihood of observing the data d for a given model m we use ℒ(d | m) = exp(−χ²/ 2), where $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}& \mathrm{=}& \frac{{\left({\mathit{T}}_{\mathrm{eff}}\mathrm{-}{\mathit{T}}_{\mathrm{eff}\mathit{,}\mathrm{obs}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{T}}^{\mathrm{2}}}\mathrm{+}\frac{{\left(\log \mathrm{\left(}{\mathit{L}}_{\mathit{\star }}\mathrm{\right)}\mathrm{-}\log \mathrm{\left(}{\mathit{L}}_{\mathit{\star }\mathit{,}\mathrm{obs}}\mathrm{\right)}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\log \mathit{L}}^{\mathrm{2}}}\\ & & \mathrm{+}\frac{{\mathrm{(}{\mathrm{[}\mathrm{Fe}\mathit{/}\mathrm{H}\mathrm{]}}_{\mathrm{s}}\mathrm{-}{\mathrm{[}\mathrm{Fe}\mathit{/}\mathrm{H}\mathrm{]}}_{\mathrm{s}\mathit{,}\mathrm{obs}}\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{[}\mathrm{Fe}\mathit{/}\mathrm{H}\mathrm{]}}^{\mathrm{2}}}\mathrm{+}\frac{{\left({\mathit{\rho }}_{\mathit{\star }}\mathrm{-}{\mathit{\rho }}_{\mathit{\star }\mathit{,}\mathrm{obs}}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{\rho }}^{\mathrm{2}}}\mathrm{·}\end{array}$In this expression for χ² observed quantities are denoted by the subscript “obs”, their standard errors are σ_T, σ_log L, etc., and other quantities are derived from the model, as described above. In cases where asymmetric error bars are quoted on values we use either the upper or lower error bar, as appropriate, depending on whether the model value is greater than or less than the observed value. The probability distribution function p(m) = p(τ_⋆)p(M_⋆)p( [ Fe/H ] _i) is the product of the individual priors on each of the model parameters. The assumed prior on [ Fe/H ] _i normally has little effect because this parameter is well constrained by the observed value of [ Fe/H ] _s so we generally use a “flat” prior on [ Fe/H ] _i, i.e., a uniform distribution over the model grid range. It is possible in our method to set a prior on M_⋆ of the form p(M_⋆) = e^αM_⋆. The prior on M_⋆ is the product of the present day mass function for planet host stars, the mass distribution of the target stars in the surveys that discovered these planets and the sensitivity of these surveys as a function of stellar mass, but since none of these factors is well determined we normally use a flat prior for M_⋆, i.e, α = 0. We also use a flat prior on τ_⋆ over the range 0−17.5 Gyr. In general, and for all the results below unless stated otherwise, we simply set set p(m) = 1 for all models where the parameters are within the range of valid values. However, the implementation of our method does include the option to set priors on [Fe/H] or τ_⋆ of the form $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \end{array}& & \end{array}$We generate a set of points m_i (a Markov chain) with the probability distribution p(m | d) using a jump probability distribution f(Δm) that specifies how to generate a trial point m′ = m_i + Δm. The trial point is always rejected if any of the model parameters are outside their valid range. Otherwise, a point is always included in the chain if ℒ(m′ | d) > ℒ(m_i | d) and may be included in the chain with probability ℒ(m′ | d)/ℒ(m_i | d) if ℒ(m′ | d) < ℒ(m_i | d) (Metropolis-Hastings algorithm). If the trial point is accepted in the chain then m_{i + 1} = m′, otherwise m_{i + 1} = m_i(Tegmark et al. 2004).

We randomly sample 65 536 points uniformly distributed over the model grid range and take the point with the lowest value of χ² as the first point in the chain. From this starting point, m₀, we find the step size for each parameter such that $\mathrm{\left|}\ln \mathrm{\right(}\mathrm{ℒ}\mathrm{\left(}{{m}}_{\mathrm{0}}\mathrm{\right|}{d}\mathrm{\left)}\mathrm{\right)}\mathrm{-}\ln \mathrm{\left(}\mathrm{ℒ}\mathrm{\right(}{{m}}_{\mathrm{0}}^{\mathrm{\left(}\mathit{j}\mathrm{\right)}}\mathrm{\left|}{d}\mathrm{\right)}\mathrm{\left)}\mathrm{\right|}\mathrm{\approx }\mathrm{0.5}$, where ${{m}}_{\mathrm{0}}^{\mathrm{\left(}\mathit{j}\mathrm{\right)}}$ denotes a set of model parameters that differs from m₀ only in the value of one parameter, j. We then produce a Markov chain with 10 000 steps using a multi-variate Gaussian distribution for f(Δm) with the standard deviation of each parameter set to this step size. This first Markov chain is used to find an improved set of model parameters and the second half of this chain is used to calculate the covariance matrix of the model parameters. We then calculate the eigenvalues and eigenvectors of the covariance matrix, i.e., the principal components of the chain. This enables us to determine a set of transformed parameters, q = (q₁,q₂,q₃), that are uncorrelated and where each of the q_i have unit variance, as well as the transformation from q to m. We then produce a Markov chain with 50 000 steps using a multi-variate Gaussian distribution for f(Δq) with unit standard deviation for each of the transformed parameters. The first point of this second Markov chain is the set of model parameters with the highest value of ℒ(m_i | d) from the first Markov chain. This second Markov chain is the one used to estimate the probability distribution function p(m | d). We use visual inspection of the chains and the Gelman-Rubin statistic (Gelman & Rubin 1992) to ensure that the chains are well mixed, i.e., that they properly sample the parameter space. The number of steps used in the two Markov chains can be varied, e.g., longer chains can be used to ensure that the parameter space is properly sampled in difficult cases.

Our algorithm is implemented as a fortran program called bagemass that accompanies this article (available at the CDS) and that is also available as an open source software project³.

3. Validation of the method

3.1. Comparison to other models

In Table 1 we compare the predicted values of T_eff and ρ_⋆ from our garstec models for a star at an age of 4 Gyr with an initial mass of 0.9 M_⊙ and solar composition to those of three other grids of stellar models. These are all grids of “standard stellar models” in the sense that they assume a linear relation between helium enrichment and metallicity, and the mixing length parameter is calibrated using a model of the Sun. The Dartmouth Stellar Evolution Program (DSEP) model grid is described in Dotter et al. (2008). The “VRSS 2006” model grid is described in VandenBerg et al. (2006) and the Geneva 2012 models are described by Mowlavi et al. (2012). The VRSS 2006 models do not include diffusion or gravitational settling of elements.

The range in T_eff values is 65 K, with the garstec models being at the top end of this range. The values of ρ_⋆ vary by about 3.5%, with the garstec models predicting the lowest density while the VRSS 2006 models predict the highest density. The differences are mainly due to differences in the assumed solar metallicity that is used for the zero-point of the [Fe/H] scale and the assumed value of α_MLT in each model grid.

We used bagemass to find the best-fitting (maximum-likelihood) mass and age estimates for these model stars based on the values given in Table 1. We assigned standard errors of 80 K to T_eff and 0.07 dex to [ Fe/H ] _i, i.e., the minimum standard errors on these values that we use for the analysis of real stars. We assumed that the error on ρ_⋆ is ± 0.001ρ_⊙. The results are also shown in Table 1. The results for the garstec models show that our method is self-consistent to better than the 1-per cent level in recovering the stellar mass and age of the star. Comparing our results to those of other models, we see that the systematic error due to differences in the stellar models in the recovered stellar mass is less than about 2%. The systematic error in the recovered age is larger (≈25%) for this example. Note that these figures may not apply to stars at different masses or ages, e.g., differences in the treatment of convective overshooting can lead to large differences in the predicted age for more massive stars.

3.2. Eclipsing binary stars

We used DEBCat⁴ to identify 39 stars in 24 detached eclipsing binary systems that are suitable for testing the accuracy of the mass estimated using our method when applied in the mass range typical for planet-host stars. The masses and radii of the stars in this catalogue have been measured directly to a precision of better than 2%. We have restricted our comparison to stars in binary systems with orbital periods P_orb> 3 d in an attempt to avoid complications that might arise due to the stars having strong tidal interactions or very rapid rotation. We also exclude stars larger than 2 R_⊙, stars without a measured value of [ Fe/H ] _s that include an estimate of the standard error, and stars with [ Fe/H ] _s< −0.4. The limits on R_⋆ and [ Fe/H ] _s were chosen so that the stars are not close to the edge of the of the stellar model grid. For the majority of planet host stars currently known the luminosity of the star has not been measured directly, so we do not include the luminosity of the eclipsing binary stars as a constraint in this analysis. The properties of the stars in this sample and the masses we derive using our method are given in Table 2. Kepler-34 and Kepler-35 are also planet host stars, with planets on circumbinary orbits.

The observed and predicted masses are shown in Fig. 1. There is a clear tendency for our grid of standard stellar models to underpredict the mass of some stars by about 0.15 M_⊙. We imposed an arbitrary limit on the orbital period of the binary systems we have used so we have investigated whether the orbital period is a factor in this analysis. The approximate orbital period of each star is indicated by the plotting symbol/colour in Fig. 1. It is clear that all the stars with discrepant masses have short orbital periods (≲6 d), but that not all stars with short orbital periods have discrepant masses. This is more directly seen in Fig. 2, where we plot the mass discrepancy as a function of orbital period. All the stars with a mass underpredicted by more than 2 standard deviations have orbital periods P_orb< 7 d, and there is a clear tendency for the mass to be underpredicted for other stars with orbital periods in this range. For low mass stars with large convective envelopes in short-period binaries such as these, the rotational period of the star and the orbital period are expected to be equal (or nearly so for eccentric binaries) due to strong tidal interactions between the stars. This has been confirmed by direct observation of detached eclipsing binary stars, including some of the stars in this sample (Torres et al. 2010). If we assume that rotation is the reason why the predicted masses of some stars are discrepant, then we can conclude that our method is able to accurately predict the mass within the stated uncertainties for stars with rotation periods P_rot ≳ 15 d and in the mass range 0.8−1.3 M_⊙. There are no suitable data for stars in long-period binaries to test whether this conclusion holds in the mass range 1.3−1.6 M_⊙. There are also no suitable data to test our method in the mass range 0.6−0.8 M_⊙, or for stars with masses 0.8−1.3 M_⊙ and rotation periods 7 d <P_rot< 14 d.

The binary star WOCS 40007 (2MASS J19413393+ 4013003) is a member of the star cluster NGC 6819. Jeffries et al. (2013) find the age of NGC 6819 to be about 2.4 Gyr from coulor−magnitude diagram isochrone fitting and an age estimate for the binary system of 3.1 ± 0.4 Gyr based on stellar model fits to their mass and radius measurements of the two stars. We obtain ages of 3.1 ± 0.6 Gyr and 3.3 ± 1.8 Gyr for WOCS 40007 A and WOCS 40007 B, consistent with both the age estimate for the binary system and the cluster obtained by Jeffries et al.

Fig. 1 Bayesian mass estimates for stars in detached eclipsing binary star systems. Symbols/colours denote the following orbital period ranges: Porb< 6 d – filled circles/green; 6 <P< 12 d – open circles/black; P> 12 d – crosses/blue.

Fig. 2 Error in the Bayesian mass estimate for stars in detached eclipsing binary star systems as a function of orbital period. Errors in excess of 2 standard deviation are plotted with filled circles.

3.3. Systematic errors – helium abundance and mixing length

The results derived using our grid of stellar models with fixed values of α_MLT and Y calibrated on the Sun are subject to some level of systematic error due to the variations in these values between different stars. Fortunately, observed constraints on the variation in both of these factors are now available thanks to Kepler light curves of sufficient quality to enable the study of solar-like oscillations for a large sample of late-type stars. The Kepler data for 42 solar-type stars have been analysed by Metcalfe et al. (2014) including α_MLT and Y as free parameters in the fit of the stellar models to observed frequency spectrum. The values of Y and α_MLT derived from the asteroseismology are shown in Fig. 3. It can be seen that there is no strong correlation between Y and α_MLT, and that there is additional scatter in these values beyond the quoted standard errors. If we add 0.02 in quadrature to the standard errors on Y then a least-squares fit of a constant to the values of Y gives a reduced chi-squared value ${\mathit{\chi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1}$, so we use 0.02 as an estimate of the scatter in Y. Similarly, for α_MLT we assume that the scatter around the solar-calibrated value is 0.2. It is debatable whether these derived values of Y are reliable measurements of the actual helium abundance of these stars or whether the derived values of α_MLT are an accurate reflection of the properties of convection in their atmospheres. However, the frequency spectrum of solar-type stars is very sensitive to the mean density of the star so it is reasonable to use the observed scatter in Y and α_MLT from the study of Metcalfe et al. as a way to estimate the systematic error due to the uncertainties in these values in the masses and ages derived using our method.

Fig. 3 Helium abundance Y and mixing length parameter αMLT for 42 solar-like stars from Metcalfe et al. (2014).

In principle it may be possible to estimate an appropriate value of α_MLT to use for a given star based on an empirical or theoretical calibration of α_MLT against properties such as mass, surface gravity, etc. We have decided not to attempt this and instead to treat the scatter in α_MLT and Y as a sources of possible systematic error in the masses and ages. The systematic error in the mass due to the uncertainty on Y is given by the quantity σ_M,Y, which is the change in in the best-fit mass, M_b, produced by applying our method using the grid of stellar models in which the helium abundance is increased by + 0.02 compared to the value used in our grid of standard stellar models. Similarly, σ_τ,Y is the systematic error in the age due to the uncertainty on Y calculated in the same way. For the systematic errors in the mass and age due to the uncertainty in α_MLT we calculate the best-fitting mass and age for a grid of stellar models calculated with garstec for α_MLT = 1.50 and multiply the resulting change in mass or age by the factor 0.2/(1.50−1.78), i.e, σ_M,α is the change in the best-fitting mass due to an increase of 0.2 in α_MLT and similarly for the age and σ_τ,α. The two grids of stellar models necessary for these calculations are included in the version of bagemass that accompanies this article so the user is free to decide whether they should account for this potential source of additional uncertainty for the star they are analysing.

4. Results

We have applied our method to over 200 stars that host transiting planet or brown dwarf companions and have found that the software runs without problems in all these cases and that our results are generally in good agreement with published results. Here we only report the results for a selection of stars to illustrate the main features of our method and to highlight some interesting results.

The input data used for our analysis is given in Table 3. For stars that have a trigonometrical parallax in van Leeuwen (2007) with precision σ_π/π ≲ 0.1 we have calculated L_⋆ using a value of f_⊕ estimated by integrating a synthetic stellar spectrum fit by least-squares to the observed fluxes of the star. Optical photometry is obtained from the Naval Observatory Merged Astrometric Dataset (NOMAD) catalogue⁵(Zacharias et al. 2004), the Tycho-2 catalogue (Høg et al. 2000) and the Carlsberg Meridian Catalog 14 (Copenhagen University et al. 2006). Near-infrared photometry was obtained from the Two Micron All Sky Survey (2MASS)⁶ and Deep Near Infrared Survey of the Southern Sky (DENIS)⁷ catalogues (The DENIS Consortium 2005; Skrutskie et al. 2006). The synthetic stellar spectra used for the numerical integration of the fluxes are from Kurucz (1993). Reddening can be neglected for these nearby stars given the accuracy of the measured fluxes and parallaxes. Standard errors are estimated using a simple Monte Carlo method in which we generate 65 536 pairs of π and f_⊕ values from Gaussian distributions and then find the 68.3% confidence interval of the resulting log (L/L_⊙) values.

The results of the analysis for our selection of planet host stars are given in Table 4, where we provide the values of the mass, age and initial metallicity that provide the best fit to the observed properties of the star, M_b, τ_b and [ Fe/H ] _b, respectively, the value of χ² for this solution, and the mean and standard deviation for each of the posterior distributions of age and mass. The likely evolutionary state of the star is indicated by the quantity p_MS, which is the fraction of points in the chain for which the central hydrogen abundance is not 0, i.e., p_MS is the probability that the star is still on the main-sequence.

Fig. 4 Left panel: results of our MCMC analysis for HAT-P-13 in the Hertzsprung-Russell diagram. Black dots are individual steps in the Markov chain. The dotted line (blue) is the ZAMS. Solid lines (red) are isochrones for ages 6.6 ± 1.2 Gyr. Dashed lines (green) are evolutionary tracks for masses 1.21 ± 0.06 M⊙. All isochrones and tracks are for [ Fe/H ] i = 0.47. Right panel: joint posterior distribution for the mass and age and HAT-P-13 from our MCMC analysis.

5. Discussion

An example of applying bagemass to one star (HAT-P-13) is shown in Fig. 4. We have chosen this star as an example because it shows the difficulties that can arise in the analysis of stars close to the end of their main sequence evolution. The complexity of the evolution tracks and isochrones in the Hertzsprung-Russell diagram is clear, even though we have restricted the models plotted to one value of [ Fe/H ] _i (the best-fit value). The joint posterior distribution for the mass and age is clearly bimodal, with peaks in the distribution near (M,τ_⋆) = (1.20 M_⊙,6.7 Gyr) and (1.37 M_⊙, 4.1 Gyr). It can be difficult to ensure that the Markov chain has converged for distributions of this type where there are two solutions separated in the parameter space by a region of low likelihood. We used another MCMC analysis for this star with 500 000 steps in the chain to verify that the default chain length of 50 000 is sufficient even in this difficult case to produce an accurate estimate of the posterior probability distribution. The maximum-likelihood solution in this case has a mass and age (M_b,τ_b) = (1.21 M_⊙,6.5 Gyr). The mean and standard deviation of the joint posterior distribution for the mass and age are (⟨ M_⋆ ⟩ , ⟨ τ_⋆ ⟩ ) = (1.21 ± 0.06 M_⊙,6.7 ± 1.2 Gyr). The good agreement between the maximum-likelihood solution and Bayesian analysis is a consequence of the low uncertainty in the value of ρ_⋆. In logarithmic terms the standard error on ρ_⋆ for this star is 0.023 dex. This is much lower than the typical uncertainty for estimates of log g based on the spectroscopic analysis of a late-type star (≈0.1 dex, Doyle et al. 2013). In general, we should expect that the terminal age bias for planet host stars for isochrone fitting using the density is much less severe than for single stars using log g. Even so, the full Bayesian analysis is worthwhile to obtain an objective estimate of the true range of possible masses and ages for a planet host star, particularly for stars near the main sequence turn off or if the error on ρ_⋆ is large. We have also found that the joint probability distribution for mass and age from the full Bayesian analysis can also be useful for improving the power of statistical tests based on these quantities, e.g., comparing the ages derived from stellar models to gyrochronological ages (Maxted et al. 2015).

HD 209458 was the first transiting exoplanet discovered (Charbonneau et al. 2000) and is also one of the brightest and best-studied. The properties of this planetary system are also typical for hot Jupiter exoplanets. The fit to the observed properties of this star including the luminosity constraint is good (χ² = 0.19 for one degree of freedom). It is interesting to compare the results for HD 209458 to those for WASP-32, which is also a typical hot Jupiter system but one for which there is currently no accurate parallax measurement. If we compare the values of ⟨ τ_⋆ ⟩ and ⟨ M_⋆ ⟩ for these stars in Table 4 we see that the adding parallax measurement with the best precision currently available does not lead to a significant improvement in the precision of the mass and age estimates – the slightly better precision of the age estimate for HD 209458 is mostly due to the higher precision of the density estimate. A similar argument applies in the case of the stars HD 189733 and WASP-52. These stars are less massive than HD 209458 and close to the limit at which age estimates from stellar models become impossible because there is no significant evolution of the star during the lifetime of the Galaxy.

For stars without a parallax measurement the number of observables is the same as the number of model parameters, but it is still possible that no models in the standard model grid give a good fit to the observed properties of the star. This is the case for Qatar 2, which is a good example of a star that appears to be older than the Galactic disc (10 Gyr, Cojocaru et al. 2014). We used the Markov chain calculated for this star to estimate an upper limit to the probability that the age of this star derived from our stellar models is less than 10 Gyr and find P(τ_⋆< 10 Gyr) < 0.002. This probability is an upper limit because the best-fit to the observed parameters occurs at the edge of the model grid (τ_b = 17.5 Gyr). This also means that the values of σ_τ,Y and σ_τ,α are not reliable in this case.

It has long been known that some K-dwarfs appear to be larger by 5% or more than the radius predicted by standard stellar models (Hoxie 1973; Popper 1997). This “radius anomaly” is correlated with the rotation rate of the star, but also shows some dependence on the mass and metallicity of the star (López-Morales 2007; Spada et al. 2013). The dependence on rotation is thought to be the result of the increase in magnetic activity for rapidly rotating stars. Magnetic activity can affect the structure of a star by producing a high coverage of starspots, which changes the boundary conditions at the surface of the star, or by reducing the efficiency of energy transport by convection. Whatever the cause of the radius anomaly in K-dwarfs, the existence of inflated K-dwarfs is likely to be the reason why Qatar 2 has an apparent age >10 Gyr. One method that has been proposed to deal with the radius anomaly is to simulate the magnetic inhibition of convection by reducing the mixing length parameter (Chabrier et al. 2007). For Qatar 2, we found that models with α_MLT< 1.4 can match the properties of this star for ages less than 10 Gyr. This phenomenological approach has some support from stellar models that incorporate magnetic fields in a self-consistent way (Feiden & Chaboyer 2013). There is currently no well-established method to predict the correct value of α_MLT to use for a given star affected by the radius anomaly and so no way to estimate the ages of such stars.

Clausen et al. (2010b) found a reasonable match to the observed masses and radii of both stars in the eclipsing binary system BK Peg using both the VRSS stellar evolution models (VandenBerg et al. 2006) and the Yale-Yonsei “Y²” models (Demarque et al. 2004), so the difference between the observed mass of BK Peg A and the predicted mass with our grid of standard stellar models needs some explanation. Both model grids used by Clausen et al. (2010b) use old estimates for the reaction rate ¹⁴N(p,γ)¹⁵O. This reaction is the bottle-neck in the CNO cycle and the relevant reaction rate has been redetermined experimentally and theoretically in the last decade, resulting in a reduction by a factor of 2 compared to previous values (Adelberger et al. 2011). We are able to reproduce this good fit to the mass and radius of BK Peg A with garstec if we use the old (inaccurate) reaction rate. With the new reaction rate and the correct value of the star’s mass the “red-hook” in the stellar evolution track as the star leaves the main sequence does not reach values of T_eff low enough to match the observed T_eff value. We also experimented with variations in the assumed convective overshooting parameter, α_ov, but this does not help to produce a good fit.

Decreasing Y or, to a lesser extent, α_MLT leads to an increase in the estimated mass for stars with mass 0.9 M_⊙. This suggests that a comprehensive analysis similar to the one presented here but applied to eclipsing binary stars will yield useful insights into stellar models for such stars. Indeed, Fernandes et al. (2012) have undertaken just such an analysis using a maximum likelihood method to estimate the helium abundance and mixing length parameter for 14 stars in seven eclipsing binary stars with masses from 0.85 M_⊙ to 1.27 M_⊙. They find a weak trend for stars with large rotation velocities to be fit better by models with lower values of α_MLT, in qualitative agreement with the results that we have found here. We have assumed that rotation is also the reason why the mass esimated using our method is too low by about 0.1 M_⊙ for some stars with masses ≳1 M_⊙ in detached eclipsing binaries. This issue deserves further investigation, ideally using Bayesian methods similar to those developed here. Given the number of model parameters and other factors such as rotation and magnetic activity that should be considered, there is a clear need for accurate mass, radius, T_eff and [Fe/H] measurements for more eclipsing binary stars. In particular, there is a lack of data for stars in long period eclipsing binaries with masses 1.3−1.6 M_⊙, i.e, at the upper end of the mass distribution for planet host stars.

6. Conclusion

The software used to produce the results in this paper is written in standard fortran, is freely available and the installation requires only one, widely-used additional software library (fitsio, Pence 1998). The method has been validated against other models and tested using the observed properties of detached eclipsing binary stars. The method and its limitations are fully described herein. As a result, it is now straightforward to determine accurately the joint posterior probability distribution for the mass and age of a star eclipsed by a planet or other dark body based on its observed properties and a state-of-the art set of stellar models.

Acknowledgments

J.S. acknowledges financial support from the Science and Technology Facilities Council (STFC) in the form of an Advanced Fellowship. A.M.S. is supported by the MICINN ESP2013-41268-R grant and the Generalitat of Catalunya program SGR-1458. The pspline package is provided by the National Transport Code Collaboration (NTCC), a US Department of Energy (DOE) supported activity designed to facilitate research towards development of fusion energy as a potential energy source for the future. The NTCC website is supported by personnel at the Princeton Plasma Physics Laboratory (PPPL), under grants and contracts funded by the US Department of Energy.

References

All Tables

All Figures

	Fig. 1 Bayesian mass estimates for stars in detached eclipsing binary star systems. Symbols/colours denote the following orbital period ranges: Porb< 6 d – filled circles/green; 6 <P< 12 d – open circles/black; P> 12 d – crosses/blue.
In the text

	Fig. 2 Error in the Bayesian mass estimate for stars in detached eclipsing binary star systems as a function of orbital period. Errors in excess of 2 standard deviation are plotted with filled circles.
In the text

	Fig. 3 Helium abundance Y and mixing length parameter αMLT for 42 solar-like stars from Metcalfe et al. (2014).
In the text

Fig. 4 Left panel: results of our MCMC analysis for HAT-P-13 in the Hertzsprung-Russell diagram. Black dots are individual steps in the Markov chain. The dotted line (blue) is the ZAMS. Solid lines (red) are isochrones for ages 6.6 ± 1.2 Gyr. Dashed lines (green) are evolutionary tracks for masses 1.21 ± 0.06 M⊙. All isochrones and tracks are for [ Fe/H ] i = 0.47. Right panel: joint posterior distribution for the mass and age and HAT-P-13 from our MCMC analysis.

In the text