© ESO, 2015

1. Introduction

Knowledge of the ages of stars with planets (SWP) is important for studying several aspects of the evolution of planetary systems, such as dynamical interactions among planets (see e.g. Laughlin & Chambers 2001) and tidal effects induced by SWP (see e.g. Pätzold et al. 2004; Barker & Ogilvie 2009). See Hut (1980, 1981) for a theoretical approach regarding tides. Moreover, SWP ages enable assessment of the evolution of the atmosphere of the hosted planets caused by chemical reactions occurring on the planets themselves and by the consequences of tidal stripping or other atmospheric loss processes. Knowledge of the stellar ages is also useful for selecting candidates for planet detections.

Most of the known SWP are main sequence G-K type stars that belong to the nearby disk field population. It is very well known that determining the ages of these stars is difficult because of the degeneracy of parameters and the slow evolution of the observational quantities. The current uncertainties are higher than the accuracy needed for these studies.

The age is not a direct observable, so its computation should use models or a combination of models and empirical relations. Methods based on isochrones from stellar evolutionary models are often used, but other methods are also applied, based on empirical relations, such as gyrochronology and activity indices. Asteroseismology will be a very promising technique when more specific data is available, and in specific cases, the chemical analysis (the so-called chemiochronology) can be applied as well. See Soderblom (2010) for a broad review on this topic.

Most of the ages of SWP come from individual sources and different methodologies. A recent discussion of the ages of nearby field stars is presented in Haywood et al. (2013), while specific analyses of SWP have been published by Saffe et al. (2005) and Takeda et al. (2007). Finally Brown (2014) presents a new study based on a new geometrical approach for interpolating the grids of isochrones, and he shows that the results from gyrochronology give systematically younger ages.

In this paper we focus on the ages derived by using isochrones. To determine the age of a stellar cluster using isochrones, it is necessary to put its stars on the Hertzsprung-Russel diagram (HRD) and evaluate – among the isochrones having the metallicity of the cluster – which isochrone best fits the layout of the stars on the diagram (isochrone fitting). Instead, determining the age of a field star, in particular a MS star, is much more complex. The statistical treatment of the data plays a crucial role in the analysis, and it is necessary to face the problem of the degeneracy of parameters. Two different methods are typically applied: the isochrone placement and the Bayesian estimation (Jørgensen & Lindegren 2005; da Silva et al. 2006; Pont & Eyer 2004).

The isochrone placement technique consists in putting a star on the HRD together with its error bars in log T_eff and log L and in properly selecting the best isochrone to account for the error box. Instead, the Bayesian estimation technique requires getting a posterior probability density function (PDF) of the age of a star, assuming an a priori star formation rate (SFR), an a priori metallicity distribution, and an a priori initial mass function (IMF).

We note that sometimes statistical instruments such as Markov chain Monte Carlo (MCMC) have been applied in the literature. MCMC tools are able to sample PDFs (without obtaining an explicit expression for them), in order to infer parameters in Bayesian ambit. It can happen that useful parameters for computating stellar ages are recovered via MCMC tools and then used to compute the age through χ²-minimization-based methods. Algorithms where a preliminary Bayesian approach is followed by a frequentistic one are not strictly Bayesian.

The paper is organized in the following way. In Sect. 2 the input data and the isochrones we used are presented; in Sect. 3 the implementation of isochrone placement and Bayesian estimation techniques is described; in Sect. 4 the results obtained are discussed; while Sect. 5 reports a summary of our work.

2. The data

First of all, this section presents the input data that were used to test the reliability of the developed algorithms. The last two subsections are dedicated to describing the sample of SWP and of the theoretical models used to characterize the SWP. From here on, all the photometric parameters are expressed in the Johnson system.

2.1. The 3.2 Gyr synthetic stars catalogue

We built a catalogue of 927 synthetic stars located on an isochrone with initial [Fe/H] = 0 and log t = 9.5 (t ≈ 3.2 Gyr). To each of them we attributed a distance d, its uncertainty Δd, and an uncertainty on log gΔlog g, generating random numbers from normal distributions with means pc, $\overline{\mathrm{\Delta }\mathit{d}}\mathrm{=}\mathrm{0.1}\mathit{d̅}$, and $\overline{\mathrm{\Delta }\log \mathit{g}}\mathrm{=}\mathrm{0.1}$ dex, respectively. The uncertainty on [Fe/H] is Δ [ Fe / H ] = 0.05. All these values are typical of the stars belonging to the SWP catalogue that is described below.

2.2. SWP Catalogue

All the planetary and stellar parameters of SWP were collected from The Site of California and Carnegie Program for Extra-solar Planet Search: Exoplanets Data Explorer¹ (Wright et al. 2011). After having discarded binaries and stars without measurements of the apparent V magnitude, B − V colour index, parallax π from Hipparcos, [Fe/H], and log g, we considered the remaining 326 stars that constitute our SWP catalogue. These stars are represented on the CMD in Fig. 1 with the two 1-Gyr-isochrones corresponding to the minimum (Z = 0.00318) and the maximum (Z = 0.054) metallicity of the sample. Increasing the metallicity, the isochrones go towards redder colours. Some stars are located to the left of the Z = 0.00318 1-Gyr-isochrone. These are anomalous because they are expected on the right-hand side, considering their higher metallicity.

2.3. The isochrones

The theoretical models employed to determine the ages of the stars are the PARSEC² isochrones (version 1.0) by Bressan et al. (2012), corresponding to the solar parameters listed in Table 1.

The different sequences of isochrones are identified by the metallicity of a star at the moment of its birth: Z_ini. We used sequences spaced by 0.05 in log t (with t in years) starting from log t = 6 up to log t = 10.1.

Considering the solar sequence of isochrones identified by Z_ini = Z_{ini, ⊙} = 0.01774 as reported by Bressan et al. (2012) (note that this value is different from the present one, which is Z_⊙ = 0.01524), the interpolation between the log L/L_⊙ and M_bol tabulated values yields the correspondence 1 L_⊙ → M_bol = 4.770. We also computed the differences between the absolute bolometric and V magnitudes given by the grid (which correspond to the bolometric corrections in the V band BC = M_bol − M_V adopted by the authors), and finally, we looked for the BC value that gave T_eff = T_{eff, ⊙} and L = L_⊙. The interpolation gives a bolometric correction for the Sun BC_⊙ = −0.063. Finally, log g_⊙ and B − V_⊙ come from the interpolation in the solar isochrone grid, as well.

As already said, the parameter that identifies a given sequence of isochrones is the metallicity Z, which is linked to [Fe/H] by an exponential relation³, which assumes the form $\mathit{Z}\mathrm{=}{\mathrm{10}}^{\mathrm{[}\mathrm{Fe}\mathit{/}\mathrm{H}\mathrm{]}\mathrm{-}\mathrm{1.817}}\mathit{.}$(1)

3. Age determination methods

3.1. Preliminary considerations

Deriving the ages of stars by making use of theoretical models requires knowledge of the stellar mass M. In fact, since isochrones corresponding to different age values can be very close on the HRD, M can help in the selection of the one that fits the input data best. In some cases, there is even a degeneracy between the age t and the mass M of a given star. Figure 2, which shows the solar metallicity evolutionary tracks of 1 and 1.05 M_⊙, clarifies the situation: the stellar parameters (log T_eff, log L) = (3.7662, 0.0839) corresponding to the intersection of the two tracks are representative either of a 1 M_⊙ star with an age of 6.34 Gyr or of a 1.05 M_⊙ with an age of 26.8 Myr, so the knowledge of the stellar mass is fundamental to correctly establishing the evolutionary stage of any given star.

In the particular case where the comparison is between ages of a pre-main sequence (PMS) star and a MS star, it is also possible to remove the degeneracy considering that PMS stars are particularly active if compared with MS stars, and this implies that they have very high rotational velocities and cromospheric activity indices. Considering $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ as reference index for the cromospheric activity, we took the age-activity relation by Mamajek & Hillenbrand (2008) and we slightly shifted it, so that solar values adopted by the isochrones of Padova match the model. According to this relation $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathrm{4.48}$ for an age t = 500 Myr and $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathrm{4.27}$ for t = 100 Myr. The typical variation between the highest and lowest peaks in activity and the average level is ~0.2 dex for a solar type star; in fact the present mean solar $\log {\mathit{R}}_{\mathit{HK,}\mathrm{\odot }}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathrm{4.91}$, while it was $\log {\mathit{R}}_{\mathit{HK,}\mathrm{\odot }\mathit{,}\mathrm{Maunder}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathrm{5.105}$ during the Maunder minimum. To be conservative we assume that if $\log {\mathit{R}}_{\mathit{HK,}\mathit{\star }}^{\mathrm{\prime }}\mathit{<}\mathrm{-}\mathrm{4.65}\mathit{,}$(2)then the star has an age t> 500 Myr, while if $\log {\mathit{R}}_{\mathit{HK,}\mathit{\star }}^{\mathrm{\prime }}\mathit{<}\mathrm{-}\mathrm{4.47}\mathit{,}$(3)then the star has an age t> 100 Myr.

As mentioned before, another indicator of the activity of a star is its rotational velocity. According to the study of Denissenkov (2010), a star younger than 500 Myr at the very least has an angular velocity Ω_⋆ ≳ 2.65Ω_{e, ⊙}, while a star younger than 100 Myr has, at least, Ω_⋆ ≳ 4Ω_{e, ⊙}, where Ω_{e, ⊙} = 2.86 × 10^-6 rad/s is the present angular velocity of the solar envelope. In the absence of the $\log {\mathit{R}}_{\mathit{HK,}\mathit{\star }}^{\mathrm{\prime }}$ value, and assuming $\frac{\mathrm{4}}{\mathit{\pi }}\mathit{v}\sin \mathit{i}$ as the mean probable rotational velocity of a star, we conclude that a star has an age t> 500 Myr, if ${\mathrm{\Omega }}_{\mathit{\star }}\mathit{<}\mathrm{7.579}\mathrm{\times }{\mathrm{10}}^{-6} \mathrm{rad}\mathit{/}\mathrm{s}\mathit{,}$(4)while its age is greater than 100 Myr, if ${\mathrm{\Omega }}_{\mathit{\star }}\mathit{<}\mathrm{1.144}\mathrm{\times }{\mathrm{10}}^{-5} \mathrm{rad}\mathit{/}\mathrm{s}\mathit{.}$(5)Besides many other parameters, the databases we examined often reported the stellar mass M, but all these parameters were derived by different authors following different calibration procedures. To make as few assumptions as possible and to produce input parameters that enter a picture that is self-consistent with the theoretical values reported by the isochrones, we decided to start from the values of

• visual magnitude V,

• colour index B − V,

• parallactic distance d,

• metallicity [Fe/H], and

• spectroscopic log g,

which are available in the literature, and then to compute all the other needed parameters (i.e. the stellar mass M using the stellar effective temperature T_eff, luminosity L, and radius R), according to the calibrations that can be inferred from the values tabulated in the isochrones.

The results are sensitive to the bolometric corrections (BCs). Several published tables of bolometric corrections are reported in the literature, but – as pointed out by Torres (2010) – values given by an author can differ noticeably from the ones given by another author depending on the arbitrary zero point (traditionally set using the Sun as reference) that each author adopts. Moreover, there is sometimes no internal consistency between V_⊙, M_{bol, ⊙}, and BC_⊙.

For these reasons we obtained the correspondence between M_V and log L by interpolating the values of the isochrone grids and inferring the BCs from the photometric values tabulated in the isochrones. Assuming an internal uncertainty on the apparent bolometric magnitude equal to 0.03 mag, the uncertainty ΔL is associated to L through error propagation. The correspondence between B − V and log T_eff is obtained from the isochrones as well. An internal uncertainty of 1% is attributed to the resulting T_eff value.

From the parameters just derived, we can now compute the estimates of R and M that are used as input data in the isochrone placement technique: $\begin{array}{ccc}\mathit{R}& \mathrm{=}& \sqrt{\frac{\mathit{L}}{{\mathrm{(}{\frac{{\mathit{T}}_{\mathrm{eff}}}{\mathrm{5778}}}^{\mathrm{)}}}^{\mathrm{4}}}}\\ \mathit{M}& \mathrm{=}& \frac{\mathit{g}}{{\mathrm{10}}^{\mathrm{4.43}}}{\mathit{R}}^{\mathrm{2}}\mathit{.}\end{array}$In these equations L and M are in solar unities, T_eff is in K, and g is in cm/s². After applying the error propagation, all the data are accompanied by the respective uncertainties.

Another aspect that has been investigated is the temporal evolution of the surface stellar metallicity Z due to atomic diffusion. The interaction between different chemical species leads to a surface depletion of elements heavier than hydrogen, which sink downwards. The characteristic timescale for the diffusion of an element is (Chaboyer et al. 2001) ${\mathit{\tau }}_{\mathrm{diff}}\mathrm{\simeq }\mathit{K}\frac{{\mathit{M}}_{\mathrm{CZ}}}{\mathit{M}{\mathit{T}}_{\mathrm{CZ}}^{\mathrm{3}\mathit{/}\mathrm{2}}}$(8)where M is the stellar mass, M_CZ the mass of the surface convective zone, T_CZ the temperature at the base of the convective zone, and K a constant referring to the chemical element being taken into account. The direction of this process is parallel to the temperature and pressure gradients, while it is antiparallel to the chemical concentration gradient. For further information see Burgers (1969) and Chapman & Cowling (1970).

We considered several stellar evolutionary tracks, kindly provided by Leo Girardi, which illustrate the evolution of stars of different masses and initial metallicities on the HRD. Such tracks give the surface Z in correspondence of any given age, so we built up different Z_k,l = Z_k,l(t) functions depending on the stellar initial metallicity and mass (identified by the subscripts k and l, respectively). We observed that atomic diffusion is not negligible for stars having masses between 0.5 and 2 M_⊙. In this mass range, the surface metallicity decreases with the time during the MS phase and then increases as the envelope convection deepens, becoming constant once the initial metallicity value Z_k,l(0) is reached. As a consequence, the initial decrease followed by the increase in Z_k,l(t) spans a shorter time scale for higher masses, where the evolution is faster. An example of the temporal evolution of Z for stars of different masses, but characterized by the same initial metallicity, is shown in Fig. 3.

The evolution of the atmospheric chemical composition has a non-negligible effect. For instance, if we take a solar star (i.e. with the current luminosity and temperature of the Sun) with a present metallicity equal to Z_⊙ = 0.01524 and employ the sequence of isochrones characterized by such metallicity without taking into account that the metallicity value that identifies the isochrones is the one that a star had when it was born and not the present one, then we get an age t = 5.1 ± 2.8 Gyr instead of the t = 4.5 ± 0.1 Gyr that we obtain if we adopt the initial metallicity value Z_{ini, ⊙} = 0.01774 to select the isochrones.

Adopting the present-day atmospheric chemical composition in selecting the isochrones for the age computation of field stars generally produces a result that is slightly biased towards older ages, especially for intermediate-age stars, as shown in Fig. 4. In fact, higher metallicity isochrones are redder than lower ones for every age, and in the MS the older a star, the redder it is. Since considering the initial chemical composition of a star implies selection of an equally or higher metallic grid of isochrones, we expect t_noZ − t ≥ 0 for MS stars. Among pre-MS isochrones, instead, older isochrones are bluer for a given metallicity: in this case, we expect negative differences. The three stars in Fig. 4 with evident t_noZ − t< 0 are peculiar because their colours are too red for their luminosity and metallicity. They are located in a region where there are both old and pre-MS isochrones. Furthermore, they do not have activity indices to disentangle between young and old ages; in these conditions, the code may give t_noZ<t. This result could be eliminated by iterating the isochrone placement twice and using in the second iteration only those isochrones that differ ± Δt from the age value t recovered by the first iteration. We are considering implementating this in the near future.

3.2. Isochrone placement

The input data characterizing each star for which we want to establish the age are listed below:

Z V BC B − V d Δd T_eff ΔT_eff L ΔL g Δg M ΔM.

To determine the age and the other parameters, such as effective temperature, luminosity, gravity, and mass, of a given star according to the theoretical models, we first considered the sequence of isochrones Is characterized by the present metallicity of the star Z_⋆. If $\log {\mathit{R}}_{\mathit{HK}}^{\mathrm{\prime }}$ or vsini of a star was available (which happens for the 94% of the stars that belong to the SWP catalogue) and relations (2) or (4) were satisfied, we discarded all the rows reporting ages lower than 500 Myr from the sequence Is. If, instead, relations (2) and (4) did not hold, but relations (3) or (5) were satisfied, we discarded all the age values available in the isochrone grids lower than 100 Myr (activity cleaning).

Then, for each isochrone in the Is sequence, we considered the point on the CMD with the minimum distance from the given star and computed all its corresponding theoretical values by interpolation through the isochrone grid values. We thus built our reduced grid of isochrones Is’. After that, we developed the following procedure:

• 1.

  For each row of Is’ (from here on, the subscript iwill always represent the row index of Is’), multiply the the-oretical value reported in the ith row of Is’ (we refer to it as thegeneric variable X_i that indicates, time by time, T_eff, L, M, or g) by thebidimensional Gaussian distribution (window function):$\mathrm{𝒢}\mathrm{\left(}{\mathit{L}}_{\mathit{i}}\mathit{,}{\mathit{T}}_{\mathrm{eff}\mathit{,i}}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\left[{\left(\frac{\log {\mathit{T}}_{\mathrm{eff}\mathit{,i}}\mathrm{-}\log {\mathit{T}}_{\mathrm{eff}}}{\mathrm{\Delta }\log {\mathit{T}}_{\mathrm{eff}}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\log {\mathit{L}}_{\mathit{i}}\mathrm{-}\log \mathit{L}}{\mathrm{\Delta }\log \mathit{L}}\right)}^{\mathrm{2}}\right]}}{\mathrm{2}\mathit{\pi }\mathrm{\Delta }\log {\mathit{T}}_{\mathrm{eff}}\mathrm{\Delta }\log \mathit{L}}\mathrm{·}$(9)In this way we can consider that the probability that a given stellar parameter corresponds to the value reported by a certain row of Is’ decreases with increasing distance on the HRD between the star and the isochrone itself, but – at the same time – the isochrones falling out of the error bars in log T_eff and log L are not discarded definitely.

• 2.

  Compute the weight $\begin{array}{ccc}{\mathit{p}}_{\mathit{i}}& \mathrm{=}& [{\left(\frac{\log \mathit{L}\mathrm{-}\log {\mathit{L}}_{\mathit{i}}}{\mathrm{\Delta }\log \mathit{L}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\log {\mathit{T}}_{\mathrm{eff}}\mathrm{-}\log {\mathit{T}}_{\mathrm{eff}\mathit{,}\mathit{i}}}{\mathrm{\Delta }\log {\mathit{T}}_{\mathrm{eff}}}\right)}^{\mathrm{2}}\\ & & { \mathrm{+}{\left(\frac{\mathit{M}\mathrm{-}{\mathit{M}}_{\mathit{i}}}{\mathrm{\Delta }\mathit{M}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\log \mathit{g}\mathrm{-}\log {\mathit{g}}_{\mathit{i}}}{\mathrm{\Delta }\log \mathit{g}}\right)}^{\mathrm{2}}]}^{-1}\end{array}$(10)that must be attributed to X_i, so that the similarity between the stellar and the theoretical M and log g values is also taken into account. The greater the likeness between the stellar and the theoretical data, the bigger p_i.

Clearly, in the previous equations all the data with the subscript i are taken from the isochrones, while the others are the input stellar parameters.

Through weighted means, it is now possible to compute the age of a given star and its temperature, luminosity, mass, and gravity according to the Padova evolutionary models. According to what has just been described, the generic stellar parameter expressed by X_⋆ results in: ${\mathit{X}}_{\mathit{\star }}\mathrm{=}\frac{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n̅}}{\mathit{X}}_{\mathit{i}}\mathrm{𝒢}\mathrm{\left(}{\mathit{L}}_{\mathit{i}}\mathit{,}{\mathit{T}}_{\mathrm{eff}\mathit{,i}}\mathrm{\right)}{\mathit{p}}_{\mathit{i}}}{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n̅}}\mathrm{𝒢}\mathrm{\left(}{\mathit{L}}_{\mathit{i}}\mathit{,}{\mathit{T}}_{\mathrm{eff}\mathit{,i}}\mathrm{\right)}{\mathit{p}}_{\mathit{i}}}\mathrm{·}$(11)The corresponding uncertainty is given by $\mathrm{\Delta }{\mathit{X}}_{\mathit{\star }}\mathrm{=}\sqrt{\frac{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n̅}}\mathrm{(}{\mathit{X}}_{\mathit{i}}\mathrm{-}{\mathit{X}}_{\mathit{\star }}{\mathrm{)}}^{\mathrm{2}}\mathrm{𝒢}\mathrm{(}{\mathit{L}}_{\mathit{i}}\mathit{,}{\mathit{T}}_{\mathrm{eff}\mathit{,i}}\mathrm{)}{\mathit{p}}_{\mathit{i}}}{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n̅}}\mathrm{𝒢}\mathrm{\left(}{\mathit{L}}_{\mathit{i}}\mathit{,}{\mathit{T}}_{\mathrm{eff}\mathit{,i}}\mathrm{\right)}{\mathit{p}}_{\mathit{i}}}}\mathrm{·}$(12)Making use of these first guesses for t_⋆ and M_⋆, among the Z_k,M⋆ functions describing the evolution in metallicity of a star with M = M_⋆, we looked for the one where Z(t_⋆) = Z_⋆, from which we recovered the initial metallicity Z_ini,⋆ that a star of age t_⋆ had at its birth if the present metallicity is Z_⋆. After that, we considered the sequence of isochrones corresponding to the just estimated metallicity Z_ini,⋆ and we iterated all these operations until a convergence in the stellar age value was reached.

3.3. Bayesian estimation

As presented by Jørgensen & Lindegren (2005), determining stellar ages from isochrones requires a comparison between observational and theoretical data, according to the stellar evolutionary model adopted. If the model-relevant parameters are collected in a vector p, while the observational data are collected in a vector q, the theoretical model gives a map from the parameter space p to the data space q. Determining stellar ages represents the inverse problem, i.e. finding a map from q to p; here we considered p = (τ,Z,m) (where τ is the age and m is the mass) and q = ( [ Fe / H ] ,log T_eff,log L).

In Bayesian statistics, the parameters that have to be estimated (in our case τ, Z and m) are treated as random variables, and their posterior (joint) PDF is $\mathit{f}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}\mathrm{\propto }{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}\mathrm{ℒ}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}$(13)where f₀ is the prior PDF and ℒ the likelihood function. The value given by f(τ,Z,m)dτdZdm represents the fraction of stars with age inside [τ, τ + dτ], metallicity inside [Z, Z + dZ] and mass inside [m, m + dm]. The constant of proportionality must be chosen so that ${}^{\mathrm{\int }}\mathrm{\int }^{\mathrm{\int }}{}^{\mathrm{\int }}\mathrm{\int }^{\mathrm{\int }}{}^{\mathrm{\int }}\mathit{f}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}$dτdZdm = 1. The integration of f with respect to Z and m gives f(τ), which is the posterior PDF that a star has the age τ. Assuming the mode as the statistical index that synthesizes the function, the best estimate for the age is the value that maximizes f(τ). Other plausible choices are those referring to the most probable age of a star considering the mean of the PDF (corresponding to the centroid of the area under f(τ)) or the median, which, instead, is the value that bisects the area under f(τ). As already pointed out by Jørgensen & Lindegren (2005) – who considered a sample of 2968 synthetic stars in order to evaluate the best criterion to assess age in Bayesian statistics – the mean and the median suffer the bias of attributing an age that is in the centre of the sequence of age values reported by the isochrone grid employed. On the other hand, the mode tends to assign the extreme age values of the isochrone grid to the stars: in particular, selection of extreme ages arises for 2031 stars over the 2968 of their entire sample, which corresponds to a frequency of ~70%. In Sect. 4 we confirm the behaviour of these three statistical indices by applying the Bayesian statistics to the stars of both the 3.2 Gyr catalogue and the SWP catalogue.

Assuming independent Gaussian observational errors ${\mathit{\sigma }}_{\mathit{n}}^{\mathrm{obs}}$ for each ${\mathit{q}}_{\mathit{n}}^{\mathrm{obs}}$, the likelihood function is given by $\mathrm{ℒ}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}\mathrm{=}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{\mathit{n}}^{\mathrm{obs}}}\mathrm{·}\exp \mathrm{(}\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}$(14)where ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{3}}{\left(\frac{{\mathit{q}}_{\mathit{n}}^{\mathrm{obs}}\mathrm{-}{\mathit{q}}_{\mathit{n}}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{n}}^{\mathrm{obs}}}\right)}^{\mathrm{2}}\mathrm{·}$(15)Following the suggestion of Jørgensen & Lindegren (2005), we assumed ${\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}\mathrm{=}\mathit{\psi }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathit{\phi }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathit{\xi }\mathrm{\left(}\mathit{m}\mathrm{\right)}$(16)where ψ(τ) is the a priori SFR, φ(Z) the a priori metallicity distribution, and ξ(m) the a priori IMF. As one of our purposes is that of studying possible evolutionary peculiarities of our stellar samples, we assume that the prior SFRψ(τ) is flat. Regarding the metallicity prior φ(Z), the spectroscopically determined [Fe/H] values we employed are very reliable (for example, ~85% of the stars belonging to the SWP catalogue has σ_{[ Fe / H ]}< 0.08 dex, with the majority of them that having σ_{[ Fe / H ]} = 0.03 dex) and, since for very high accuracy metallicities (say σ_{[ Fe / H ]} ~ 0.05 dex) the likelihood function acts as a narrow window function that substantially suppresses the contribution of φ(Z), we consider a flat φ(Z) (see Pont & Eyer 2004, for a detailed discussion about the choice of the prior distributions). Assuming also that ξ(m) = m^-2.7, which is representative of the empirical IMF at around 1 M_⊙ (Kroupa et al. 1993), it is possible to obtain f(τ) as $\mathit{f}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{\propto }\mathit{G}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{\int }\mathrm{ℒ}\mathrm{\left(}\mathit{\tau ,Z,m}\mathrm{\right)}\hspace{0.17em}\mathit{\xi }\mathrm{\left(}\mathit{m}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{m}\hspace{0.17em}\mathrm{d}\mathit{Z}\mathit{.}$(17)Implementing the algorithm of the Bayesian determination of age, we first check that the condition $\begin{array}{c}\min \end{array}{\mathit{\chi }}_{\mathit{i}}^{\mathrm{2}}\mathit{<}{\mathit{\chi }}_{\mathrm{0.99}}^{\mathrm{2}}$(18)is satisfied; ${\mathit{\chi }}_{\mathrm{0.99}}^{\mathrm{2}}$ is the 99th percentile value of the chi-square distribution that in our case is 11.345, because we have three degrees of freedom (d.o.f.). Since the probability that χ² ≤ 11.345 for a 3 d.o.f. chi-square distribution is 99%, relation (18) states that we do not evaluate ages for stars whose input parameters have less than 1% probability (according to chi-square distribution) to actually be those measured. These stars are characterized by data points that are far away from any isochrone, so their PDFs turn out to be meaningless.

If the preliminary condition given by (18) is satisfied, then we numerically evaluate G(τ) considering a set of sequences of isochrones taken at constant steps of 0.05 dex in [Fe/H] within the interval whose bounds are empirically fixed at ± 3.5σ_{[ Fe / H ]} from the stellar metallicity, where σ_{[ Fe / H ]} is the uncertainty on [Fe/H]. Let m_jkl be the initial mass value read at line l of the isochrones grid of age τ_j and metallicity Z_k, then $\mathit{G}\mathrm{\left(}{\mathit{\tau }}_{\mathit{j}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}}\sum _{\mathit{l}}\mathrm{ℒ}\mathrm{\left(}{\mathit{\tau }}_{\mathit{j}}\mathit{,}{\mathit{Z}}_{\mathit{k}}\mathit{,}{\mathit{m}}_{\mathit{jkl}}\mathrm{\right)}\mathit{\xi }\mathrm{\left(}{\mathit{m}}_{\mathit{jkl}}\mathrm{\right)}\mathrm{(}{\mathit{m}}_{\mathit{jkl}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{m}}_{\mathit{jkl}}\mathrm{)}\mathit{.}$(19)Once we have obtained the vector G(τ) of components G(τ_j), we find the component that assumes the maximum value (say ) and divide each component by , obtaining the normalized function $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}$. After we have smoothed $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{G}\end{array}$ through a polynomial interpolation, the most probable age attributed to the star is the value that maximizes $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{G}\end{array}$ (modal value). We also compute the mean age as the age coordinate of the centroid of the area under $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{G}\end{array}$ and the median age as the age value that bisects the area under $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{G}\end{array}$.

As described by Jørgensen & Lindegren (2005), it is possible to prove that $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{G}}_{\lim }\end{array}\mathrm{=}\mathrm{0.61}$ sets the 68% confidence level of , so we provide a 68% confidence interval [ τ₁,τ₂ ] to be the shortest interval such that $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{G}\end{array}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathit{<}\mathrm{0.61}$ outside it⁴.

4. Discussion of the results

4.1. Isochrone placement vs. Bayesian estimation

The first aim of the paper is to compare the reliability of the developed algorithms in the computation of ages of field stars. To reach this goal, we applied these algorithms taking the 3.2 Gyr synthetic stars one by one. Then we checked the correspondence (B − V, M_V) → (log T_eff, log L) representing each star on both the colour magnitude diagram (CMD) and the HRD. The calibration is correct if a star has the same relative position with respect to the same reference isochrones on both the two diagrams. Finally, we checked to what extent the global age distribution of the entire sample is consistent with the expected age of 3.2 Gyr attributed to the stars a priori. The calibration between observational and theoretical parameters is perfect except for the RGB region, where the intersection between different isochrones makes it difficult; this will have consequences on the age determination.

In fact, the representation of the age of these stars determined with the isochrone placement technique versus their absolute magnitude M_V (Fig. 5) shows that – as expected – the majority of the stars fall in the horizontal region between 3 and 3.5 Gyr, but some deviations occur for very low and very high magnitudes. The ages of some very bright stars (in the RGB phase) can be imprecise because of the difficulties linked to the calibration between observational and theoretical stellar parameters owing to the particular shape of the isochrones in that region. Instead, the errors for the ages obtained for the faintest stars deals with the intrinsic difficulties in estimating the ages of low MS stars, even if the calibration from the CMD to the HRD is well done. Just to summarize, the algorithm implementing the isochrone placement technique gives reliable ages, except for

• some stars in the RGB phase, whose ages could also be completely wrong;

• low MS stars, with the obtained ages that can differ up to 50% from the correct value.

We decided to further analyse the ages of the synthetic stars belonging only to the MS (i.e. stars having 4 <M_V< 8), because they are the most common one amongst the SWP. In addition, since extremely young ages can be discarded considering the stellar activity, the ages are obtained by removing the theoretical isochrones with ages lower than 500 Myr from the fit procedure.

The histogram in Fig. 6 shows a comparison between the age distributions computed with the isochrone placement technique and the Bayesian estimation, which adopts the mode as synthesis index. It also reports the mean error bars associated to the results. Figures 7 and 8 show the age distributions of the same synthetic stars as were obtained through the Bayesian estimation technique that adopts the mean and the median as synthesis indices, respectively.

The isochrone placement turns out to be the most reliable technique, since all the stars fall in the bin [3, 3.5] Gyr (which contains the value of age of 3.2 Gyr a priori attributed to these stars) with a typical uncertainty in the age determination of ~0.15 Gyr, corresponding to an error of ~5%.

On the other hand, using the Bayesian statistics, the distribution derived from the modal age values presents the main peak in correspondence of [0.5, 1] Gyr. This is the outmost interval of the isochrone grid, and it is representative of an age that is lower than the correct one. Finally, the distribution obtained using the mean and the median as synthesis indices of the Bayesian PDFs are quite similar, essentially spanning an age range from 1 to 7 Gyr. The determination of ages through the mean or the median is therefore not very accurate, and confirming the conclusion by Jørgensen & Lindegren (2005) and by Takeda et al. (2007), these statistical indices tend to select ages in the middle of the sequence of age values reported by the isochrone grids, centralizing the distribution.

The mode is definitely the indicator to be preferred and its tendency to select the extremes in age values in the isochrone grids could be partially mitigated operating a proper numerical filtering. In fact, in some cases – an example of one of them is shown in Fig. 9 – the Bayesian PDF shows peaks corresponding to very low ages. Dealing with synthetic stars, in this context we can recognize such peaks as spurious, and we realize that they hide the presence of the peak centred at ~3 Gyr, which indicates the correct age. However, if no a priori indication is given about the ages of stars (the ordinary situation if the scientific aim is to determine of the ages of field stars), it is not possible to select the true peak.

Finally, we introduced a random Gaussian perturbation in the input M_V, B − V and log g of the synthetic stars, considering Gaussian distributions with 3σ = 1% of the unperturbed values. We show the age distributions deriving from the isochrone placement and the modal Bayesian age in Fig. 10. As expected, the distributions are broader and the mean error bars of the output results are higher. The mean age value of the isochrone placement does not change. The behaviour of the two techniques is similar to what we have just said about the ages derived from unperturbed input values.

In conclusion, the isochrone placement technique is the method chosen to compute the ages of the SWP.

4.2. Stars with planets ages

We analysed all the 326 stars belonging to the SWP Catalogue. Only nine of them have been removed because their observational parameters were not consistent with the theoretical ones. The age distribution of the remaining 317 stars determined using the isochrone placement technique is presented in Fig. 11. About 6% of the stars are younger than 0.5 Gyr, and then the distribution reaches a peak at [1.5, 2) Gyr and after that it generally decreases. There is a non-negligible number of stars (~7%) older than 11 Gyr. All the parameters of these stars are presented in Table 4.

The histogram of Fig. 11 shows that there are SWP with all possible ages, with a preponderance of stars with ages <~6 Gyr. The median value (~4 Gyr) appears slightly lower than the age of the Sun. Of course, older MS stars are fainter (making it more difficult to identify planets). It is not obvious that there is a selection bias for the ages in Fig. 11. Very old or very young stars can have large uncertainties in age because the oldest isochrones overlap the pre-MS region. The frequency of stars in the [0, 0.5) Gyr bin is still an open question. The majority of the youngest SWP are hot Jupiter (HJs) hosts, and they tend to have high rotational velocity, so the check of activity did not discard the youngest isochrones. Actually, it might happen that such high rotational velocity does not reflect the stellar activity, while it can be due to the spin-up induced by the HJ, as suggested for example by Poppenhaeger & Wolk (2014). Without definite indications from activity indices, another way to disentangle pre-MS from post-MS isochrones is to consider the stellar density ρ_⋆, which is observationally available in the case of SWP detected through the transit method (see e.g. Sozzetti et al. 2007). Here, ρ_⋆ gives indications of the evolutionary stage of a star, and it has been used, for instance, by Rouan et al. (2012) to reject an extremely low age for CoRoT-23. We are planning of inserting this kind of check in the future development of our algorithm.

The literature reports the ages of 189 SWP, estimated using isochrones, but with different theoretical models and techniques. The superimposition of the consequent literature age distribution on what is found here is shown in Fig. 12. The two distributions are quite in agreement for ages older than ~4.5 Gyr, while we found more stars in the domain of younger ages. In particular, in our sample, ~5% of the stars has an age between 0 and 0.5 Gyr, while no stars in the literature fall in this first bin.

The difference t_Lett − t between the ages found in the literature and the values computed here is represented in Fig. 13, which shows a wide spread in the age values. We provide new ages, which have the advantage of being derived using the same method and the same set of isochrones, and therefore useful for statistical investigations. Our age distribution in Fig. 12 is broader than in the literature. Our data show an overabundance of young stars. Some young stars in the literature are judged to be older by our technique. This may happen for stars located on the red side of the isochrone interval, where the locus of the old isochrones is very near to the pre-MS one.

Just for completeness, we also show the age distributions of 302 stars belonging to the SWP catalogue from the modal Bayesian estimation. Given that we find high frequencies in the first and in the last age bin of Fig. 14, it is again clear that the mode always tends to assign extreme age values to the stars. Figure 15 again suggests that the mean always produces a distribution centred on the middle of the age range. We avoid reporting the age distribution of SWP inferred from the median Bayesian age value, because it is very similar to Fig. 15 where the mean is employed as reference statistical index.

Figure 16 shows the difference t_BayesMode − t between the modal Bayesian age t_BayesMode and the age t computed through the isochrone placement represented versus t. Apart from the effect of the mode of attributing the most extreme age values available in the isochrone grids (visible from the locus of points in the upper right part of the figure), the other Bayesian-estimated age values appear slightly biased towards older ages. The median of the isochrone placement age distribution is ~4 Gyr and the mean error is ~1.15 Gyr, while the Bayesian modal age distribution has a median of ~4.25 Gyr, and the mean error is ~1.75 Gyr. This behaviour of the method contrasts with what is reported by Haywood et al. (2013), who used both a χ² minimization (somewhat comparable to our isochrone placement technique, even if they considered the χ² minimization simply on the (T_eff, M_V) plane) and the Bayesian technique described by Jørgensen & Lindegren (2005) to compute the ages of nearby field stars. In their paper, the Bayesian technique does give younger ages, so we again emphasize the crucial role played by the details within the specific method implemented to compute the ages of field stars.

In the upper left-hand side of Figs. 13 and 16 there is a single isolated point, which is the star CoRoT-23. The paper that discusses its detection (Rouan et al. 2012) reports V_⋆ = 15.63 mag and B_⋆ = 16.96 mag. V_⋆ also agrees with the magnitude reported by exoplanets.eu⁵. Instead, exoplanets.org reports V_⋆ = 16.96 mag. Since exoplanets.org reports B − V_⋆ = 1.33 (same value as inferred from Rouan et al. 2012), it is likely that B_⋆ = 16.96 ≠ V_⋆. However, after applying the isochrone placement with the new photometry, our code does not converge on any age value. This new photometry, in fact, moves the star farther out from the set of isochrones. This means that, first of all, there is a problem in the source of the photometry. Moreover, CoRoT-23 is a peculiar system considering its age and the orbital eccentricity of the hosted planet (see Rouan et al. 2012). As already said, taking the observational ρ_⋆ into account may help give a better answer to this problem.

Either way, using V_⋆ = 16.96 mag, we obtained an isochrone placement age t_⋆,Isoc = 3.3 Myr, while the Bayesian estimation gives t_⋆,Bayes = 11.3 Gyr. This is why CoRoT-23 appears as an outlier in Fig. 16. Thus t_⋆,Bayes ≫ t_⋆,Isoc because, even if the Bayesian PDF had a peak corresponding to a very low age, such a peak was judged as a spike by the smoothing through the polynomial interpolation, so it was erased. The next major peak was then at 11.3 Gyr.

Finally, we considered the work of Brown (2014), who used different isochrones and gyrochronological relations to assess the ages of a sample of SWP. We found 24 SWP in common with Brown (2014). The difference t_BrownYY − t between the ages t_BrownYY computed by Brown using the YY isochrones (Demarque et al. (2004)) and the ages t computed here through the isochrone placement are represented versus t in Fig. 17. Considering the small sample of stars, the spread is consistent with the uncertainties of the ages of MS stars. The outlier in the upper left-hand side in Fig. 17 is WASP-2. It is in the pre-MS region of the CMD, but we do not have any activity index, which could allow us to discard the youngest isochrones.

As stated by Brown (2014) in his paper, isochrones tend to give ages older than gyrochronology. If, instead, we employ our determination of the SWP ages t through the isochrone placement technique, we find that isochrones can also give ages younger than gyrochronology as clarified by Fig. 18, where the gyrochronological ages t_BrownGyro obtained from the relation of Barnes (2010) are used. Considering the typical age uncertainties, the agreement with Brown’s gyrochronological ages is very good with a lower dispersion if compared with the ages derived from YY models. It is likely that those stars, which turn out to be the oldest from the isochrone placement, do not have an accurate gyrochronological age considering that the age from the rotational velocity to age is very uncertain after some billion years.

4.3. Impact of the input parameters on the output ones

In this section we briefly discuss the sensitivity of the output parameters derived through the isochrone placement technique to the parameters assumed as input. We simulated the Sun, adopting the input parameters listed in the input column of Table 2. These are the same parameters as the isochrones we used to fit of the Sun. We attributed the typical uncertainties of the SWP catalogue stars to them, i.e. Δlog g = 0.1 dex, ΔL comes from the error propagation of Δm_bol = 0.03 mag and ; we have ΔL ~ 17%, ΔR ~ 10%, and ΔM ~ 40%. The output parameters obtained according to the isochrones are shown in the output column of Table 2. They are consistent with the expected solar parameters, and among them, the age is correct and characterized by very high precision.

To evaluate how changes in both the photometry and the spectroscopy propagate into the results, we perturbed the input V and B − V of the Sun alternatively by ± 0.01 mag and log g by ± 0.01 dex. Table 3 presents what we obtained, also listing the percentual variations Δ_oo of the output data assuming as reference values the output values obtained without any perturbation of the input parameters.

Table 3 shows that a variation of only 0.01 mag in V or B − V can lead to variations in the mean output values up to ~40% in the age and up to ~2% in the other parameters. This level of sensitivity for the derived age to the input photometry is understandable considering that we chose a MS star to perform these tests. In this region of the HRD, the isochrones are very close, so that even a small variation in magnitude and colours completely moves onto an isochrone that corresponds to a different age value. It is reasonable that if the star were in a different region of the HRD, where the isochrones are not so closed, the change in age induced by the perturbation of the input photometry would not be evident like this. Anyway, a very precise photometry is required in order to attribute the correct age to a star. On the other hand, even a slight perturbation in the input log g by only 0.01 dex induces a variation in the ages up to ~7%, while the variations in the other parameters can be considered almost negligible.

A final observation concerns the uncertainties accompanying the age values reported in Table 3, which also turn out to be an order of magnitude greater than the uncertainty of the output age of the unperturbed Sun. This is because the perturbations trigger so-called artificial stars, and the big uncertainties reflect the difficulty of properly matching all the input parameters in the isochrone grids, where parameters referring to stars that are supposed to exist are tabulated. On the other hand, since the physics underlying the isochrones is well performed and the Sun does exist, entering with consistent input stellar properties gives high precision output.

5. Conclusions

In this paper we uniformly derived the ages of 317 planet-hosting stars. We checked the reliability and accuracy of two techniques (isochrone placement and Bayesian estimation) from a sample of 3.2 Gyr synthetic stars. We found that the isochrone placement technique produces the expected age value. Instead, the estimation of age through the Bayesian statistics (using an explicit PDF) suffers the problem of primarily selecting the extreme age values available in the isochrones grid if we adopt the modal value of the PDF. It is worth emphasizing that this bias is present, although we used stars with well-defined ages. The mean or the median values of the PDF produce a bias towards the middle of the available age range and the age dispersion is large.

We found that ~6% of the stars with planets are younger than 0.5 Gyr: the reliability of this frequency will be subject to further investigations. The age distribution shows a peak in correspondence of the bin [1.5, 2) Gyr, then it decreases and ~7% of the stars are older than 11 Gyr. Approximately 60% of the stars in this sample are younger than 5 Gyr.

We found that a perturbation of the input V or B − V by only 0.01 mag can lead to variations up to 40% in the estimated age of a MS Sun-like star. These are the typical uncertainties that characterize the stars in our sample. Instead, a perturbation in the input log g by 0.01 dex can lead to variations up to 7% in the estimated age. The final accuracy depends on the combination of the input errors. Important uncertainties and biases can also be produced by the adopted stellar models, in particular by treating the element diffusion. We found that one of the most important consequences of the element diffusion is the choice of the correct isochrone metallicity because the initially metallicity can be considerably higher than the metallicity of the stellar atmosphere after some billion years for solar-type stars. If this effect is ignored, the ages can be overestimated up to about 25%.

Finally, we confirmed the presence of some stars with planets located on the left-hand side of the solar main sequence (as also noted by Brown 2014) and suggested further photometric observations to identify the source of the problem. Once the observational data is confirmed, we will investigate the nature using models that implement rotation. In fact, in the typical range of luminosities of MS stars with planets, MS isochrones that consider rotation are bluer than models that do not take rotation into account, as can be seen in Yang et al. (2013).

Online material

---

Acknowledgments

We thank the anonymous referee for the useful, thoughtful and deep comments that helped us in improving our paper. V.N. acknowledges partial support from INAF-OAPd through the grant “Analysis of HARPS-N data in the framework of GAPS project” (#19/2013) and “Studio preparatorio per le osservazioni della missione ESA/CHEOPS” (#42/2013). This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org.

References

All Tables

All Figures