© ESO 2018

1. Introduction

Amidst the labyrinthic zoology of stellar types and classes, one subset has gained considerable attention over the past few decades. These stars are called solar twins. Even though, with such a name, what they should be seems obvious, defining what they really are has so far been an ever-evolving process. Of course, one expects that a solar twin should have a physical state as close as possible to the Sun. This is actually the reason behind the relatively recent interest for these stars. They were first identified as a group by Cayrel de Strobel et al. (1981), based on spectroscopic arguments. Therefore, one can immediately see that for a “good” solar twin to be classified as such will depend on the precision one can reach to estimate its physical properties. This also explains that, before this pioneering work, this class of stars has largely remained ignored. Being G2V objects, they are, on average, relatively faint, hence demanding large telescopes and high-resolution spectrographs in order to obtain a good precision on atmospheric parameters. Conversely, alongside the on-going improvement of spectroscopic instrumental methods, the threshold for classification as solar twin has evolved considerably and many observational campaigns have been carried over in order to detect these stars (Porto de Mello & Da Silva 1997; King et al. 2005; Meléndez et al. 2006, 2010, 2014; Meléndez & Ramírez 2007; Takeda & Tajitsu 2009; Datson et al. 2012, 2014; Porto de Mello et al. 2014; Ramírez et al. 2014; Mahdi et al. 2016). As of today, roughly a hundred stars can be classified as solar twins, somewhat depending on the exact criterion retained for classification.

Studying solar twins offers multiple perspectives. On a statistical level, they offer a good benchmark for solar-like populations. Some studies have focused on the properties of the Sun itself, trying to determine if it was an outlier with respect to some solar-twin samples (Gustafsson 1998, 2008; Meléndez 2014; Dos Santos et al. 2016). Exploiting further this idea, some other studies have explored potential planet-star connections using samples of solar twins. It was suggested for instance that planet-hosting stars were deficient in refractory elements¹. This is also true for the Sun itself, which is refractory-deficient with respect to most of the known solar twins (Ramírez et al. 2010). They were also used to study the problem of Li depletion in planet-hosting stars (Israelian et al. 2009). Studying a sample of approximate twins (in the sense that their ages span a range that encompasses largely the solar age), Baumann et al. (2010) claimed that this trend is purely evolutionary and is not correlated to the presence of planets (see also Monroe et al. 2013; Carlos et al. 2016).

From a stellar-modelling point of view, solar twins are extremely interesting targets. We first recall that stellar evolution codes were largely developed using the Sun as a reference observational benchmark. Indeed it is the only star for which we have extremely precise measurements, independent of modelling, for age (Bouvier & Wadhwa 2010) and mass (Olive & Group 2014). Furthermore, the observational precision on its radius (Emilio et al. 2012) and luminosity (Fröhlich & Lean 2004) allows to calibrate solar models and assess their accuracy through helioseismic inversion (see e.g., Thompson 1991; Basu & Christensen-Dalsgaard 1997; Christensen-Dalsgaard 2002; Basu 2016).

As a result, we know with excellent confidence the main physical processes at play in the solar interior. However, these solar models are usually built on assumptions such as spherical symmetry, neglecting rotation or magnetic fields, that may quickly break down when the mass of the star varies. Likewise, the sizes of convective envelopes or cores depend strongly on the stellar mass and evolutionary phase. For low-mass stars, the outer convective envelope becomes deep enough that the radiative core disappears. This may change the stability properties of the star against small perturbations (Gabriel 1964, 1967; Rodríguez-López et al. 2012). For masses larger than the Sun, the convective envelope rapidly becomes very thin while convective cores start to appear. Their modelling is somewhat uncertain. This can be due to numerical issues in treating simultaneously microscopic diffusion and nuclear reactions in a convective region (see for instance Christensen-Dalsgaard 2008a). It might also be the consequence of other phenomena such as double-diffusive convection, in that case taking the form of semiconvection (Moore & Garaud 2016).

A good property of solar twins is that they shall not enter such regimes. At the same time, they might differ slightly from the Sun. Hence, by studying them, we can be confident that the general assumptions made for our stellar models hold. But at the same time, we can test this model by letting the stellar structure vary.

In this paper, we focus on 18 Sco, the brightest solar twin. It was the first observed by asteroseismology (Bazot et al. 2011, 2012). This allowed to estimate its global seismic properties such as the large separation, but also individual frequencies. In this paper we want to assess the impact of such measurements on our knowledge of the physical state of 18 Sco. In particular, we want to obtain statistically robust estimates of its global physical parameters but also of the related uncertainties.

Fig. 1. Examples of problematic PDFs of individual frequencies for n = 19, l = 0 (left panel) and n = 18, l = 2 (right panel). The black lines represent the histogram estimate and the red line a kernel-density estimation.

In Sect. 2 we present the statistical model we use to estimate these parameters. We also present the two estimation methods, namely Bayesian probability density estimation and frequentist optimisation, that we selected to carry out the estimation. Importantly, we reassess the seismic data given in Bazot et al. (2012), which needs to be treated carefully before being incorporated in a statistical model. In Sect. 3 we present our results and discuss them.

2. Statistical model

2.1. Generic formulation

By statistical model, we mean a mathematical formulation of the behaviour of our observations. A simple approach is to consider an additive model made of a deterministic and a random part, and hence having the form

$$\begin{array}{c}\hfill \mathit{X}=\mathcal{S}(\mathit{\theta })+\mathit{ϵ},\end{array}$$(1)

where X are the observations, $\mathcal{S}(\mathit{\theta })$ is the deterministic part (being a function of θ, the stellar parameters we want to estimate) and ϵ the realisation of a random variable that represents the observational noise.

In the following, we will make the (somewhat optimistic) assumption that our model is unbiased, that is that the noise component has zero mean and hence that the stellar model $\mathcal{S}$ is the expectation value of our observations. We assume that the observations are independent. The likelihood has thus the form

$$\begin{array}{c}\hfill \pi (\mathit{X}|\mathit{\theta })={\displaystyle \prod _i}\pi (X_i|\mathit{\theta }),\end{array}$$(2)

where the X_i are the components of the vector X.

2.2. Data and the physical model

2.2.1. Non-seismic data

One of the main advantages of studying solar twins is that stellar spectra can be analysed differentially relative to the solar one. The departures from the latter are small enough that they can be treated as a first-order perturbations. This assumes that a stellar spectrum and a solar one have been obtain from the same instrument (Meléndez et al. 2014). This translates in turn into smaller uncertainties on the atmospheric parameters T_eff, logg and [Fe/H] than those usually quoted for other stars (see e.g., Meléndez et al. 2014; Nissen 2015; Spina et al. 2016).

In this work we considered three recent estimates of the atmospheric parameters of 18 Sco from Meléndez et al. (2014), Nissen (2015) and Spina et al. (2016). They are given in Table 1. Assuming that data points have Gaussian parent distributions with the same mean but different standard deviations, a reasonable new estimate of these atmospheric parameters is the weighted average of the sample. The variance of the weighted mean was then used to compute the associated uncertainties. The resulting parameters are T_eff = 5817 ± 4 K, logg = 4.448 ± 0.012, [Fe/H] = 0.052 ± 0.005. For modelling purposes, we used only the effective temperature and surface metallicity, choosing to constrain our model with the luminosity rather than logg. The conversion between [Fe/H] and Z/X, which we effectively use as the output of the stellar code, was performed using the solar ratio (Z/X)_⊙ from Grevesse & Sauval (1998). We note that different values for this ratio have been derived since then (Asplund et al. 2005). However, these lower estimates raised the, still unsolved, solar-abundance problem, which put solar models and solar observations at odds (see e.g., Guzik et al. 2006; Castro et al. 2007; Basu & Antia 2008; Antia & Basu 2011; Gough 2012; Basu 2016). For the sake of conciseness and simplicity we chose to overlook this issue in the present study.

For the luminosity, we selected the value given by Boyajian et al. (2013). It is based on an estimate of the bolometric flux F_b = (17.34 ± 0.09)×10⁻⁸ erg s⁻¹ cm⁻². This was obtained using an aggregate of various photometric fluxes and compared to a library of stellar spectra (Pickles 1998). The corresponding luminosity is L = 1.0438 ± 0.0120 L_⊙.

A natural choice for the radius is the one derived by Bazot et al. (2012). It was obtained using the CHARA interferometric array and the PAVO interferometer. Its value is 1.010 ± 0.009 R_⊙.

The luminosity and radius were obtained using the HIPPARCOS parallax (van Leeuwen 2007). We note that a much higher value of 1.166 ± 0.026 R_⊙ has been derived by Boyajian et al. (2012). Their subsequent modelling leads however to far too large ages for 18 Sco. Consequently we decided not to take it into account through averaging (as we did for the effective temperature). Some recent analyses of angular diameters derived by Casagrande et al. (2014) and White et al. (2018) claim some possible systematic errors for stars observed by Boyajian et al. (2012) that are not very well resolved, justifying our decision not to use the radius measurement. Finally, note that using the independently-determined, but strongly-correlated observations effective temperature, radius and luminosity provides consistency check between the existing constraints. It is indeed not always easy to find a model that reproduce them all, as shown by the case of α Cen A (Miglio & Montalbán 2005).

For the non-seismic data, we always choose Gaussian densities. Therefore, the likelihood for these observations is proportional to

$$\begin{array}{c}\hfill exp[-\frac{1}{2}{\displaystyle \sum _{i=1}^{n_{\mathrm{ns}}}}\frac{{(X_i-{\mathcal{S}}_i(\mathit{\theta }))}^2)}{{\sigma }_i^2}],\end{array}$$(3)

where i labels the non-seismic observations and n_ns is their number.

2.2.2. A reassessment of the seismic data

The frequencies provided in Table 3 of Bazot et al. (2012) could potentially be used to derive seismic indicators such as the small separations (see below). However, one needs to clearly understand what these estimates mean. In order to do so, one must go back to the output of the Markov chain Monte Carlo (MCMC) simulation that was used to estimate the oscillation frequencies of 18 Sco. A thorough examination indicates that the marginal posterior density functions (PDFs) for the ν_n, l are extremely complex, mostly reflecting the noisy nature of the data and the difficult spectral window induced by the ground-based observations. Typical examples of such distributions are shown in Fig. 1. Therefore, even though the estimates given in Bazot et al. (2012) provide quantitative assessments on the oscillation frequencies, they are difficult to use in statistical models as such, that is to provide a likelihood for asteroseismic diagnostic. They are statistical summaries that only very partially capture the shape of their underlying parent distribution.

Fig. 2. PDF for the small separation δν0, 19. The black lines represent the histogram estimation, the red line the kernel-density estimation and the green one the result of mixture model fitting.

We chose not to derive the seismic indicators based on the estimates of the individual frequencies given in Bazot et al. (2012). This is motivated on one hand by the potential difficulty there is when facing likelihoods with multiple modes², such as shown in Fig. 1. These could lead to many degeneracies in the solution to the estimation problem. On the other hand, dealing with individual frequencies demands to take into account surface effects (see e.g. Kjeldsen et al. 2008). As a general rule, this could be problematic for any given star (Bazot et al. 2013).

A much more robust approach consists in analysing directly the posterior probability density of these seismic indicators themselves. This is extremely straightforward since one simply has to combine the MCMC samples for the individual frequencies to obtain the sample of a frequency separation or ratio. For instance, the two samples $\{{\nu }_{n,l}^{(1)},\cdots ,{\nu }_{n,l}^{(N)}\}\sim \pi ({\nu }_{n,l}|\mathit{y})$ (where y stands for the observed time series) and $\{{\nu }_{n-1,l+2}^{(1)},\cdots ,{\nu }_{n-1,l+2}^{(N)}\}\sim \pi ({\nu }_{n-1,l+2}|\mathit{y})$ allow us to obtain the sample $\{{\nu }_{n,l}^{(1)}-{\nu }_{n-1,l+2}^{(1)},\cdots ,{\nu }_{n,l}^{(N)}-{\nu }_{n-1,l+2}^{(1)}\}\sim \pi (\delta {\nu }_{n,l}|\mathit{y})$. Here we used the definition δν_n, l = ν_n, l − δν_{n − 1, l + 2} for the small separation. These were chosen because they are mostly sensitive to the innermost regions of the star (Tassoul 1980; Roxburgh & Vorontsov 1994).

It turns out that these samples are much easier to study than those of individual frequencies. An example of such a situation is shown in Fig. 2. A comparison with Fig. 1 shows that most of the multiple maxima found in the marginal PDFs for ν_19, 0 and ν_18, 2 are not found in the corresponding PDF for δν_19, 0 = ν_19, 0 − ν_18, 2. This is easy to understand if one considers that from one iteration of the MCMC algorithm to the other the value of these frequencies might vary greatly. However, their average small separation will have to remain centred around the same value in order to reproduce the seismic data.

This first step considerably simplifies our analysis. We also stress that any subsequent study of 18 Sco using the HARPS data should adopt this approach. However, our marginal PDFs for the small separations remain complicated. We noticed that some of them depart from normal densities. In most cases, they display important asymmetries and, sometimes, multiple modes. A very interesting approach in order to manipulate these distributions is to model them as mixtures of normal distributions (see e.g. Frühwirth-Schnatter 2006). Formally if we consider a random vector Z with density p_Z, we can model the latter as follows

$$\begin{array}{c}\hfill p_{\mathit{Z}}(\mathit{Z})={\displaystyle \sum _{j=1}^M}{P}_j\mathcal{N}({\mathit{\mu }}_{\mathit{j}},{\mathbf{\Sigma }}_{\mathit{j}}),\end{array}$$(4)

with $\mathcal{N}({\mathit{\mu }}_j,{\mathbf{\Sigma }}_j)$ a multivariate normal distribution with mean μ_j, a vector of size N = dim(Z), and covariance matrix Σ_j, of size N × N, and subject to the constraints

$$\begin{array}{c}\hfill {\displaystyle \sum _{j=1}^M}{P}_j=1\text{and}{P}_j\ge 0.\end{array}$$(5)

Equation (4) gives a Gaussian mixture model for a vector. Ideally, this is how one would treat any seismic indicator used to constrain a stellar model. For instance, when dealing with the so-called frequency ratios r₀₁ or r₁₀ (Roxburgh & Vorontsov 2003) one has to take into account that any values evaluated for order n₁ and n₂, such that n₁ ≠ n₂, are correlated through the frequencies that enter their computations (see e.g. Silva Aguirre et al. 2013). This is not true, however, of the individual small separations which are uncorrelated for different values of the couple (n, l). Lets call X^sis that vector that regroups all the seismic indicators we wish to reproduce. In the case its components are only small separations, we can model each one with a separate Gaussian mixture model. This is a very convenient simplification. Indeed, a popular approach to estimate the parameters of a mixture model is the expectation-minimisation (EM) algorithm. As it turns out estimating the parameters of a Gaussian mixture model of the form (4) using an EM algorithm becomes increasingly difficult when the dimension of the problem increases. Using small separation allows us to bypass this technical difficulty. We could then use a separate Gaussian mixture model to model each individual small-separation likelihood

$$\begin{array}{c}\hfill \pi (\delta {\nu }_{n,l}|\theta )={\displaystyle \sum _{j=1}^M}{P}_j\mathcal{N}({\mu }_j,{\sigma }_j^2).\end{array}$$(6)

If we map (bijectively) the (n, l) couples for which we have measured a small separation onto a single index i then the seismic data vector can be written X^sis = (δν₁, …, δν_i, …, δν_N). Here, N is the number of observed small separations. The corresponding likelihood is then

$$\begin{array}{c}\hfill \pi ({\mathit{X}}^{\mathrm{sis}}|\theta )={\displaystyle \prod _{i=1}^N}\left[{\displaystyle \sum _{j=1}^{M_i}}{P}_{j,i}\mathcal{N}({\mu }_{j,i},{\sigma }_{j,i}^2)\right].\end{array}$$(7)

The parameters of the Gaussian mixtures were estimated separately, for a given i, following the simple version of the EM algorithm given in Bishop (1995). In Table A.1 we give the main characteristics of the mixture model we used to describe the small separations.

A mixture model has the advantage, in the framework of MCMC sampling, of being easy to compute, since it has a simple analytic form. However, a word of caution is in order concerning the current implementation of this methodology in our analysis. One should note that we are using a rudimentary approach to mixture model fitting. In particular, we evaluate, by visual inspection, the number M of normal distributions to be included in the sum on the right-hand side of Eq. (4).

2.2.3. ASTEC

Our model for the evolution and oscillations of the stellar structure, $\mathcal{S}$, is composed of the Aarhus STellar Evolution Code (ASTEC) and adipls. Both have been extensively described in the literature (Christensen-Dalsgaard 1982, 2008a, b). Here we simply state the main settings we adopted.

We assumed a non-rotating, non-magnetic star. The opacities and equation-of-state tables in which we interpolate are taken from the OPAL collaboration, respectively from Iglesias & Rogers (1996) and Rogers & Nayfonov (2002). Nuclear reaction rates were taken from Angulo et al. (1999) with the additional inclusion of the values obtained by the LUNA collaboration for the ¹⁴N(p,γ)¹⁵O reaction (Formicola et al. 2004). Diffusion was included for He and heavy elements. These latter are treated as a block. It is fine to do so with ASTEC as long as we do not try to model stars with convective cores (Christensen-Dalsgaard 2008a).

Many parameters can be tuned in ASTEC. We only let a small subset vary, namely the mass, M, the age, t_⋆, the initial metallicity and hydrogen mass fraction, Z₀ and X₀. We also have the mixing-length parameter α as a free parameter. The latter is the proportionality coefficient between the mean-free path of a fluid element, in the mixing-length picture as described by Böhm-Vitense (1958), and the pressure scale height. Therefore, $\mathcal{S}$ is a mapping from a subspace of the parameter space to the observation space in which the vectors are respectively θ = (M, t_⋆, X₀, Z₀, α) and X = (T_eff, L, Z/X, R, {ν_n, l}).

2.3. Bayesian estimation

2.3.1. Bayesian statistical model

Bayesian density estimation is one of the two strategies we adopt to obtain values for the parameters of 18 Sco. To that effect, we shall supplement our parametric statistical model with a prior density for the parameters as per Bayes’ formula

$$\begin{array}{c}\hfill \pi (\mathit{\theta }|\mathit{X})\propto \pi (\mathit{\theta })\pi (\mathit{X}|\mathit{\theta }).\end{array}$$(8)

Here π(X|θ) is the likelihood, which is in our case computed as the product of Eqs. (3) and (7). π(θ|X) is the Posterior Density Function of the parameters conditional on the data, X, and π(θ) is the aforementioned prior density. The former is the object we are ultimately interested in, since its knowledge allows us to use the tools of statistics to provide estimates of the parameters. The latter is the fundamental feature on which Bayesian statistics rely and thanks to which one can, in practice, switch from the observations to the parameters being the random quantities in the problem.

The prior density encodes the information we possess on the parameter before carrying over the estimation. Its use has been the subject of many discussions and debates for many decades. These are far outside the scope of this paper. For the present study, suffice to say that one always needs to specify carefully the prior density considered. Indeed, its formulation will condition the final outcome of the estimation process. As general rule, using two different priors π₁(θ) and π₂(θ) in (8) shall ultimately result in different a posteriori estimates ${\widehat{\mathit{\theta }}}_1$ and ${\widehat{\mathit{\theta }}}_2$. That being said this does not mean that Bayesian Statistics are more “subjective” than a frequentist method. They simply provide a way to formalise assumptions one may have on the outcome of the estimation process, for instance forbidden regions in the parameter space or previous independent measurement on some of the parameters. In stellar physics a typical example are stellar mass estimates for members of close binaries (see for instance Bazot et al. 2016).

The prior information on the stellar parameters is sparse. In this study we only used uniform densities, their properties are given in Table 2. The only parameter for which reliable prior measurements exist is the mass. A previous study by Bazot et al. (2012) gives an estimate of the mass, namely 1.02 ± 0.03 M_⊙. This was obtained by combining density and radius estimates through homology relations. Therefore, this estimate is based on the same data we are using here and a much cruder physical model than ASTEC. For these reasons we decided not to use it as a prior, but simply compare it to the results inferred from the PDF.

Regarding the other stellar parameters, only two of them have clear cut upper limits. The age of 18 Sco ought to be smaller than ∼13 Gyr, the age of the Universe³. For X₀ helium measurements have shown that the earliest galaxies have an helium mass fraction ≤0.25 (Olive & Skillman 2004; Aver et al. 2013). If we neglect the metal abundance in these very old galaxies at the epoch of their formation, we can set an upper bound on the initial hydrogen mass fraction, X₀ = 0.75.

The other parameters are less-well constrained. In practice, upper and lower bounds can be obtained using test runs of our MCMC algorithm (see Sect. 2.3.2 below). After analysing their outcome, we can redefine the domain of definition of our prior by excluding regions of the space of parameters in which we are confident that models will not be accepted by the algorithm. This empirical procedure using MCMC test runs was used to set the priors for Z₀ and α. The resulting domains are sometimes very large with respect to the region in which the marginal PDF significantly differs from zero. This is largely due to the fact that if a parameter is only allowed to vary over a narrow region, then the approximate PDF could be artificially increased close to the boundaries. This is due to the MCMC algorithm trying to go past the upper or lower limits and thus getting “swamped” near the boundaries. Therefore, only sound physical arguments, such as the ones given for the upper bounds for t_⋆ and X₀, shall motivate strong restrictions in the space of parameters, and we chose to err on the safe side for the numerical setup.

One may wish to include an Y₀ − Z₀ relation in the prior. These have been observed previously. Some studies have shown, for instance, that there exists a linear relation between galactic abundances of helium and metals (Izotov & Thuan 2004; Fukugita & Kawasaki 2006; Balser 2006; Casagrande et al. 2007), and hence between X₀ and Z₀. Such relations have been used previously for stellar modelling (see e.g. Deal et al. 2017). However, we do not wish to incorporate a priori correlations between the parameters but would rather study them a posteriori. Thus we retain independent uniform priors for both Z₀ and X₀.

Finally, the case of α is a difficult one. The current formulation of the mixing length is somewhat heuristic, adopted in order to provide a convective flux in one-dimensional models. Numerical simulations far more precise than the one used here exist and have shown the mixing-length parameter to remain fairly constant across the HR diagram, at least for solar analogues (Trampedach et al. 2014). So far, these simulations have not been used directly to fit stellar observations. Some hybrid one-dimensional stellar codes that interpolate in the tables obtained from three-dimensional simulations have developed in order to obtain solutions in the upper stellar layers, and in particular the superadiabatic layers (Sonoi et al. 2015; Ball et al. 2016; Houdek et al. 2017; Jørgensen et al. 2017). These were mostly developed to account for surface effects on oscillation frequencies but could be of great interest to simply provide a more robust formulation of surface convection. Nevertheless, we do not have such a code implemented together with an MCMC interface. Moreover, the testing required to ensure good performance of a stellar code interpolating in a grid of three-dimensional atmosphere in the context of Bayesian estimation is outside the scope of this study. This implies that we have to let the mixing-length parameter vary significantly, since there does not exist a physically sound reason to limit it (Gough et al. 1977). Likewise, it would be poor practice to set the bounds of a uniform prior based on other numerical simulations. Thus we set the upper and lower bounds of α using the aforementioned empirical approach.

We finally assume that, besides X₀ and Z₀, all the parameters in the priors are uncorrelated, that is

$$\begin{array}{c}\hfill \pi (\mathit{\theta })=\pi (M)\pi (t_{\star })\pi (\alpha )\pi (Z_0)\pi (X_0)\pi ({\alpha }_{\text{ov}}).\end{array}$$(9)

2.3.2. Sampling method

Recovering the posterior density function π(θ|X) is the main technical issue of the estimation process. There is no closed-form solution to this inverse problem that gives θ as a function of X. We note that our use of the term “inverse problem” differs here from the more restricted scope encountered in solar physics (see e.g. Christensen-Dalsgaard 2002). We consider a much broader meaning as can be found for instance in Tarantola (2004) or Gregory (2005).

Numerical methods are thus in order and in this paper we adopted an MCMC algorithm to carry out the estimation. The details of the method are given in Appendix A.1. Their viability in the framework of stellar modelling has been discussed by Bazot et al. (2008, 2012, 2016). It has already been emphasised in these studies that fine tuning of an MCMC algorithm for stellar parameter estimation might turn out to be a subtle matter. These papers dealt with α Cen A, which is a component of a binary system. As such, we have a strong prior on its mass, which greatly facilitates the sampling, restraining significantly the relevant space of parameters.

2.4. Local optimisation

It is important to compare the resulting parameters and their uncertainties using classical local optimisation methods with a Bayesian one. The former are used widely in the literature and the shortcomings of such approaches need to be quantified and understood, and in particular, the (under-)determination of proper uncertainties.

Using the exact same code set-up as described in Sect. 2.2.3, we also used the Powell algorithm to find local solutions. This algorithm has the advantage that it can minimise any function where the uncertainties on the data are not necessarily described by simple Gaussian distributions. Using the baseline dataset X = (T_eff, L, [Fe/H],R, {δν}_n, l) (Sect. 2.2) we proceeded to find local solutions by optimising the likelihood that appears in Eq. (8). The one difference between the two parameter sets θ is the use of the initial helium abundance Y₀ instead of the initial hydrogen abundance X₀ in the optimisation. As X₀ + Y₀ + Z₀ = 1 this has no influence on the result.

Unlike the Bayesian approach, a local method suffers badly from correlations in parameters, that is it will find a solution close to the initial parameters if two of the free parameters are degenerate, whereas the Bayesian method will correctly extract all of these parameter correlations and additionally provide a better framework for interpreting the results.

Due to this local problem, we chose to work in a reduced 3-d parameter space to optimise M, t₀, Z₀ while fixing X_0, α. Such an approach is typical and necessary (Miglio & Montalbán 2005; Creevey et al. 2007, 2012; Stello et al. 2009; Dogan et al. 2013; Lebreton & Goupil 2014). The optimisations were repeated using many initial guesses of the parameters as well as using different combinations of the latter.

The estimation of our uncertainties is based on 1. generating a small grid around the optimal parameters, 2. perturbing the observations by their uncertainties and 3. finding the model from the grid that matches best to the perturbed observations. These simulations were repeated 10 000 times, and we used the resulting distributions of 1D parameters to describe their mean parameter and symmetric uncertainty.

3. Results and discussion

3.1. MCMC results

In Appendix A.1, we show that our MCMC simulations have converged to an acceptable level. We can thus merge the results from the independent chains we ran and obtain a posterior density for θ. The baseline case we analyse is for X = (T_eff, L, [Fe/H],R, {δν}_n, l) and the uniform prior on the mass as explained in Sect. 2.3.1. In Table 3 we show the estimates for the individual stellar parameters based on the corresponding marginal densities. We also display the 68.3% credible intervals intervals on the parameters. In the following, whenever a density is Gaussian the 68.3% credible interval will be given using the symbol “±”. In this case, the 68.3% credible interval then coincides with the 1σ credible interval. Otherwise, the intervals are summed up as asymmetric error bars. These are defined as the smallest intervals containing the Maximum A Posteriori (MAP) and for which the posterior density integrates to 0.683 (see Bazot et al. 2016, and references therein).

It can be noted immediately that the mass estimate is in very good agreement with the one given in Bazot et al. (2011). This is particularly interesting since, as already noted, these have not been obtained using the same assumptions for the stellar model. Of course, without a proper modelling of the star, one cannot make statements on the other parameters. Nevertheless, it is a nice a posteriori confirmation for homology techniques (Gough 1990).

Fig. 3. Two- and one-dimensional marginal densities for the stellar parameters of 18 Sco for X = (Teff, L, [Fe/H],R, {δν}n, l), with the uniform prior on the mass given in Sect. 2.3.1. The red shades in the central panels represent the two-dimensional marginal densities (normalised). The side panels represent the one-dimensional ones.

The two- and one-dimensional densities are shown in Fig. 3. Another notable feature in Fig. 3 is seen in the side panels representing the one-dimensional marginal densities. These appear to be far from Gaussian. In particular the one for the age is multimodal. Consequently, establishing the statistical summaries such as those given in Table 3 demands some care. For unimodal distributions, we report the MAP estimate. The uncertainties are given as the smallest interval containing the MAP estimates for which the parameters have a 68.3% probability to lie in. We define the level of precision as the estimated as the ratio of the length of the 68.3% credible interval to the MAP estimate. This is a global estimate that does not account for any asymmetry in the density. For our baseline case, we obtained levels of precision of the order 6% on the mass, 6% and 9% on X₀ and Z₀, and 35% on α. The precision on the age is discussed below in greater details.

Fig. 4. Marginal joint probability density for (M, X0) (magnified from Fig. 3). The full lines mark locii of constant ZAMS luminosity. The dashed line shows the slope as obtained by Pearson’s correlation coefficient.

Examining the correlation coefficients for each pair of parameters shows that they all correlate to some degree. However, we can identify two groups that display significant correlations. As a rule of thumb, we consider as significant a correlation for which the Pearson’s correlation coefficient (Pearson 1895) is ≳0.7. First, the mass correlates strongly with X₀. Second, t_⋆, α, X₀ and Z₀ are all tightly correlated. The correlation between the mass and X₀ can be explained by the adverse effects these two parameters have on the luminosity on the zero-age main sequence, L_ZAMS. Mass-luminosity relations obtained from simple homology considerations (see e.g. Clayton 1968) establish clearly that an increase of M induces an increase of L_ZAMS. Conversely, the dependence of opacity on X₀, often assumed being a power law, implies that an increase in hydrogen-mass fraction corresponds to a luminosity decrease. In Fig. 4, we display the posterior joint density for the couple (M, X₀). Overplotted are lines of constant luminosity on the ZAMS and the slope obtained using Pearson’s correlation coefficient. For this latter we retain its classical interpretation as the geometric average of the two regression slopes of M by X₀ and X₀ by M (Rodgers & Nicewander 1988). We see that all of these have the same direction. For that reason, we can associate the correlation between these two variables to be caused solely by the need to balance their effects on the initial luminosity. Since the mass does not correlate strongly with the other parameters, we can assume that once L_ZAMS is determined through the values of Z₀, α and X₀, then a value of M is imposed.

The correlation between X₀ and Z₀ is intuitive and corresponds to setting the initial metal-to-hydrogen ratio so that the evolved star can reproduce the observed value of [Fe/H]. On a side note, we can notice that the behaviour of these two parameters is opposite to what is seen in Bazot et al. (2016) in the case of α Cen A. This is due to the fact that there exists a strong prior on the mass for this latter star. Therefore, it is not possible to set an adequate initial luminosity by varying M and X₀ and this is instead achieved through relative adjustments of X₀ and Z₀. This shows how complex can be the dependence of the final posterior estimates on the precise functional form of the statistical model.

The correlations of α with Z₀ and X₀ can be partly explained by considering the effects of varying the mixing-length parameter. This affects mostly two characteristics of a stellar model: the depth of its convective zone and its effective temperature (and consequently its radius). The correlation between α and Z₀ thus sets the initial metal density in the convective zone, that is the ratio of Z₀ to the size of the convective zone. This is an important quantity because, together with the initial Z₀/X₀ it characterises the amount of metals that diffusive process ought to deplete the external convective zone in order to reproduce the observations. The correlation between α and X₀ is partly governed by the initial effective temperature. Because in the framework of the mixing-length theory an increase of α decreases the temperature gradient, then, all other things remaining equal, it will increase T_eff. Decreases of T_eff induced from increases of X₀ can be deduced from homology relations which predict that the effective temperature is proportional to some power of the molecular weight. However, the initial effective temperature is not as strongly linked to the (α, X₀) correlation as is the ZAMS luminosity to the (M, X₀). This indicates that other parameters, such as Z₀, have an influence on the initial T_eff. Such intricate interplay are difficult to disentangle.

The same can be said for the correlation between the stellar age and α, Z₀ and X₀. The age of the star is mostly controlled by the need to reproduce the internal layers of 18 Sco and thus its seismic characteristics. In that sense, it is the tight constraints on the age that imply the aforementioned correlations for the other parameters, as shall be discussed below. On that point suffice to say that Fig. 3 shows that the age of the star increases with α, Z₀ and X₀. The two first correlations are diffusion effects. Indeed, reaching the required amount of metal takes longer when Z₀ increases and diffusion becomes slower when, everything otherwise equal, the depth of the convective envelope becomes larger. The correlation with X₀ is related to the energetics of the star since the ZAMS luminosity is on average lower and that more time is necessary to reach the observed luminosity.

The case of the age demands a more careful discussion. It has been established that its marginal density is bimodal. To account for this we used a Gaussian Mixture Model. We found that a reasonable agreement is obtained for four Gaussian components (this allows to accommodate reasonably well for the bimodality but also for clear asymmetries in the densities). One mode peaks at 4.46 Gyr and the other at 6.71 Gyr. We can describe both modes satisfactorily using two of the components given in Table 3. The peak at 4.46 Gyr is well-described with the two modes with the smallest means and the one at 6.71 Gyr by the other two. Taking this into account, the weights of the lower- and upper-age modes are respectively 0.54 and 0.46. Therefore, one cannot conclude clearly on whether one of these two solutions is more likely than the other. Coming back to the problem of estimating credible intervals, one can separate the two peaks using the results of the Gaussian Mixture model fitting. In that case, it is necessary to renormalise the weights obtained, since only two components are used for each mode. Using such a procedure, we obtain as a MAP estimate $4.{67}_{-1.29}^{+0.87}$ Gyr. Likewise, for the upper-age mode, we derive a credible interval $6.{95}_{-0.89}^{+1.81}$ Gyr. The relative precisions are 46% for the former mode and 39% for the latter. We note that if we estimate the age in the sense of the Posterior Mean, then the two solutions are even further apart (4.18 ± 1.10 Gyr and 7.60 ± 1.17 Gyr).

Fig. 5. Examples of observed (full lines) and posterior (dashed lines) probability densities for the small separations δν15, 0 (left panel) and δν19, 0 (right panel).

From what preceded, it seems obvious that the bimodality observed in the age marginal density of our main result stems from the very nature of the seismic data we used. If we look at the joint probability of the age and the theoretical individual small separations (not displayed here), we indeed see that they are strongly (anti-)correlated. In order to understand how the small-separation measurements affect the age, we need to examine the adequacy between the theoretical and observational density of the individual small separations. In short, we want to assess whether or not we could reproduce the seismic data.

Looking at these densities, one sees immediately that it is difficult to model them properly. Two examples are shown in Fig. 5. The MCMC-simulated densities always show a bimodality that maps the age bimodality. In order to get a feeling of the closeness between the observed and theoretical densities, we can compute the Kullback–Leibler distance, which is defined as the distance between two probability densities p and q

$$\begin{array}{c}\hfill D_{\mathrm{KL}}(p,q)={\displaystyle {\int }_{-\infty }^{+\infty }p(x)log\left(\frac{p(x)}{q(x)}\right)\mathrm{d}x.}\end{array}$$(10)

The theoretical densities shown in Fig. 5 are those with the smallest (l = 0,  n = 15) and largest (l = 0,  n = 19) D_KL. None of these reproduce perfectly the observed ones. However, for the case l = 0,  n = 15 the absolute distance between the mean is 0.34 and the variance ratio is 2.08. Those values are 1.36 and 0.02 for the l = 0,  n = 19 small separation. This indicates that the model reproduces much better the former than the latter.

It is in fact those small separations with the highest D_KL that cause the bimodality of the age density. We sketch an explanation in Fig. 6. In there we plot the small separations obtained for the best models found in the MCMC sample for t_⋆ >  5 Gyr (full line) and t_⋆ <  5 Gyr (dashed line). For the sake of readability, we did not represent the full distribution of the small separations. Each corresponds to a different peak in the bimodal age density. We also plot the observational means and variance to provide an idea of the agreement between these local best models and the observed densities of the small separations. We see there that only five observed small separations are compatible with both models. Of the 13 remaining small separations, two are far closer to the older model, in particular δν_19, 0, which is likely to control the old-age solution. The other ones, are mostly compatible with the lower-age solution, even though they have larger variances than the distribution of δν_19, 0 and thus only partially compensate the impact of this latter.

From this discussion, we can conclude that the main source of error in our result does not come from the modelling of the observational errors. Rather, it is the data themselves that impose limitations on our analysis, either because of their intrinsic properties, that is the noise, or because their modelling in Bazot et al. (2012) was not accurate enough. The density of the small separations as we could derive them from MCMC samples in Sect. 2.2.2 suggests larger variations with the mode order (or the frequency) than what the stellar models can accommodate. Therefore, using only the data at hand cannot allow us to make a choice on which solution for the stellar age is the preferred. Either new data or a re-analysis of the existing time series with more adequate techniques is needed to go further.

3.2. Effect of observational constraints

In order to test the effect of the observations as constraints, we also ran three MCMC simulations with different observation vectors: X = (T_eff, L, [Fe/H],{δν}_n, l), X = (T_eff, L, [Fe/H],R), X = (T_eff, L, [Fe/H]). The resulting estimates of the parameters are given in Table 4. In Fig. 7 we show the two- and one-dimensional marginal densities of the stellar parameters for the case X = (T_eff, L, [Fe/H],{δν}_n, l). In Fig. 8 we display the same graphs for the two cases in which the seismic data were not included. We did not split up the spectro-photometric data since we considered that without these basic observations together, it is not possible to get any useful estimate of the stellar parameters. As discussed below, this assumption seems justified a posteriori.

Fig. 6. Observed and estimated small separations. The dashed lines show the small separations of the best model with t⋆ <  5 Gyr. The full lines show the small separations of the best model with t⋆ >  5 Gyr. The red dots show the mean of the observed small separations and their standard deviations. The full dots mark small separations with l = 0 and the open dots those with l = 1.

The case X = (T_eff, L, [Fe/H],{δν}_n, l) shows the effect of adding a radius measurement to the basic spectro-photometric data. We first note that the bimodal structure of the age density is preserved. Second, we see a significant decrease of the precision on the mass, which now drops to 15%. This a consequence of the loss of information on the average density of the star (Creevey et al. 2007; Bazot et al. 2011). Given the existing correlation between M and X₀, it is unsurprising that precision also decreases for this latter parameter, down to 12%. However, the decrease in precision is not as large as it is for the stellar mass. An explanation is that the lack of constraint on the final average stellar density allows for correlation between the mass and other stellar quantities besides X₀. We indeed see in Fig. 7 that it now correlates marginally with t_⋆, α and Z₀. The correlation with the age is a well-known trend in stellar physics. On isochrones, luminosity and effective temperature increase with the stellar mass. Therefore, to reproduce a similar data set with a larger mass, one needs to decrease the age. In retrospect this behaviour sheds light on the very small correlation seen between M and t_⋆ seen in Sect. 3.1. It indicates that, in the regime defined by the observations of 18 Sco, the age–mass relation is not as steep as the mass-X₀ relation and that the effects of the former can only be seen when M is allowed to vary on wider ranges. The correlations with α and Z₀ reflects, as above, the need to set up the ZAMS model adequately, only this time with one more degree of freedom. Rather than being fixed by balancing M and X₀, the initial luminosity is now the result of an interplay between M, α, X₀ and Z₀. Consequently, we do not observe such a good alignment between lines of constant L_ZAMS and the geometric average of the regression lines in the (M, X₀) plane as seen in Fig. 4 for the case including the radius measurement.

The precisions on α and Z₀ do not change, even though the MAP estimates do. This stresses that the radius does not provide such an important constraint on these parameters. It was already noticed by Creevey et al. (2007) that understanding the relation between the error on the radius and the uncertainty on α is difficult and depends on the details of the model. In the case of 18 Sco, this relative independence can be understood by the fact that the critical quantity upon which α and Z₀ act is not the radius, but the density of metals in the convective zone. This is not constrained by the radius but rather by the surface [Fe/H] ratio. Looking at the results for X = (T_eff, L, [Fe/H],R) and X = (T_eff, L, [Fe/H]) in Fig. 8 seems to confirm this. First we see that the age distribution does not display any significant bimodality. The loss of precision on the age compared with cases that include seismic data is very significant. If the radius is included precision is of the order 96%, otherwise it is of the order of 118%. This implies a small loss of precision in α, but not of the same magnitude. There is, however, almost no loss of precision in Z₀ and X₀. In Fig. 8, the two-dimensional marginal PDF that preserves a structure relatively similar to those seen in Figs. 3 and 7 is the one for (α, Z₀). This sheds light on the the role of the very precise measurement of [Fe/H] that exists for 18 Sco. This true regardless of whether or not the radius has been included as a constraint.

3.3. Comparison with local optimization

After presenting the results from the Bayesian approach, we aim to compare the resulting model parameters with those obtained from local optimisation and attempt to quantify the differences in the uncertainties using the two approaches. Optimal models and uncertainties from the local method are given in Table 5. These have been obtained for the baseline case with X = (T_eff, L, Z/X, R, {ν_n, l}). As can be noted the uncertainties are much smaller compared to those presented in Table 3, being based on fixing X₀ and α. In all cases, the optimal parameters are in good agreement with the marginal distributions obtained from the MCMC simulations. Comparing these values directly with Fig. 3, we can see that by restricting 0.692 <  X₀ <  0.711 (0.287 >  Y₀ >  0.268) the local uncertainties that we obtain are not entirely underestimated. For X₀ = 0.692 for example, solutions are found between 1.00 and 1.05 M_⊙, in agreement with the values of M = 1.041 M_⊙, 1.035 M_⊙ and 1.028 M_⊙ proposed in Table 5.

Fig. 7. Same as Fig. 3 but with X = (Teff, L, [Fe/H],[δν]n, l).

Fig. 8. Same as Fig. 3 but with X = (Teff, L, [Fe/H],R) and X = (Teff, L, [Fe/H]) in the left and right panel, respectively.

It is interesting to visualise the real discrepancy between them knowing that our assumptions on the local analysis are indeed unrealistic. There is no possibility to constrain the values of α and X₀ given the current (and most likely future) observations used in this work. The upper row in Fig. 9 illustrates the 1D marginal distributions for M, t_⋆ and Z₀ as inferred from the Bayesian analysis (black curves), along with the inferred parameters and uncertainties from local optimisation. The parameters and uncertainties from the local analysis are represented by the red Gaussian distributions for the three solutions provided in Table 5. We note that in particular, we entirely fail to obtain a solution in the second age range proposed by the Bayesian approach.

Some of this discrepancy can be explained by the need to fix some parameters in the local optimisation procedure. To test this we restricted the MCMC sample to ranges α ± 0.01 for all cases. We see in the bottom row of of Fig. 9 the corresponding distributions for all three values of the mixing-length parameter given in Table 5. A significant agreement is then reinstated for the age, although this means that the large-age solution observed in the original sample has been filtered out by our cut on α. Likewise, the mass distribution is also much closer to the one estimated from local optimisation. However, the metallicity density, whose estimates obtained from the local optimisation and MCMC strategies looked fairly similar, now becomes much narrower in the former case.

This highlights the difficulty there is to find a proper agreement between the two approaches. It is extremely difficult to rule out any similarity between the two outcomes as not being incidental. This is due to the fact that, while the dimension of the problem increases, it becomes more difficult to keep track of correlations between parameters This is the so-called curse of dimensionality. To that issue, MCMC algorithms offer a better operational solution, due to their ability to explore stochastically the space of parameters. They could potentially be used to serve as a benchmark for less-time-consuming local optimisation strategies.

Of course, this explanation of the discrepancy in terms of fixed parameters does not account for the difficulty to identify a second mode in the marginal age PDF using local optimisation. In this case, optimisation algorithms, which provide, by definition, point estimates, naturally underperform.

3.4. Other stellar parameter estimates

In Sect. 3.3 we have estimated the stellar parameters using a local optimisation strategy. The estimates there are in agreement with all those presented in the previous section. However, the relevant quantities are not only the point estimates for the parameters but also the uncertainties one can associate to these values. Those quoted in Tables 3–5 are consistent for Z₀. For the other parameters, they differ much, sometimes close to an order of magnitude. It is noteworthy that the uncertainties obtained from optimisation do not vary much when the observation vector X changes. This indicates that a lot of information is factored in the assumptions made in Sect. 2.4 for the derivation of the uncertainties and that this may lead to underestimating them.

We have seen above that it is hard to reproduce perfectly the seismic data, and that its inclusion leads to a double solution to the estimation problem. One could extrapolate to a case in which more accurate seismic data would be available and expect lower uncertainties on the age, that could come closer to those of Table 5. Nevertheless, even in that case, a proper modelling of the uncertainties points towards uncertainties of the order of ∼1 Gyr, which remains more than twice those found with optimisation. The convergence results from Appendix A.1 show that the samples generated from the MCMC are reliable, therefore one should clearly study carefully the details of the estimation strategy chosen to obtain stellar parameters before trusting the uncertainties. Contrary to what has been done in previous so-called “hare-and-hounds” exercises, this comparison only focuses on the differences in the methodology used to obtain the parameters. This means that we have used the exact same data and code (and code setup), which is not always the case in other comparison studies (see e.g. Reese et al. 2016).

Comparison with previously derived stellar parameters for 18 Sco is difficult precisely for this reason. It becomes very hard to disentangle the effect of the estimation strategy, the data that constrain the model and the precise numerics of the stellar evolution code used. We can point out a few recent estimates given in the literature, limiting ourselves to the age, which is the parameter the most difficult to assess but that could be crucial, in particular for studies that focus on Li depletion on the main sequence (Israelian et al. 2009; Meléndez et al. 2010). Carlos et al. (2016) give an estimate of 3.8 ± 0.5 Gyr. This was obtained by comparison with the Yonsei-Yale (Kim et al. 2002) set of isochrones and using only spectrophotometric constraints. A notable difference is that logg was considered instead of the luminosity. The estimates of the atmospheric parameters and their associated uncertainties also differ slightly. However, it remains extremely unlikely that these changes could account for the difference with the uncertainties in the range 3–4 Gyr obtained with X = (T_eff, L, [Fe/H]). At any rate, this should not allow to obtain uncertainties lower than those obtained using seismic data. Other recent estimates can be found in Ramírez et al. (2014) and Spina et al. (2018), they give $t_{\star }=3.0_{-0.6}^{+0.3}$ Gyr and $t_{\star }=4.2_{-0.5}^{+0.3}$ Gyr respectively. These uncertainties are again much lower, by an order of magnitude, than what is found using our method. These results were also obtained using isochrone fitting procedures, together with a Bayesian Statistical model. Only spectrophotometric parameters were considered, Spina et al. (2018) using both logg and the luminosity. A likely explanation for such a discrepancy is the difficulty to sample properly the space of stellar parameters using pre-computed isochrones (Bazot et al. 2012). Therefore, some models are not taken into account either due to incomplete sampling or because some stellar parameters have been fixed, thus reducing the final variance. As a sanity check we notice that the models from the MCMC simulation reproduce satisfactorily the probability density of the atmospheric parameters, indicating that a wide range of stellar ages can indeed account for such a combination.

Fig. 9. Upper row: marginal distributions (black) of M, t, and Z0 from the results presented in Fig. 3 (Sect. 3), and the distribution for these same parameters as determined from the local analysis (red Gaussians, Sect. 2.4, Table 5). Bottom row: marginal distributions for the same parameters, but restricted to a small interval around the values of α given in Table 5. The red lines are similar to those in the upper row.

Other studies provide ages derived using stellar population statistics. Some have focused on the so-called chemical-clocking methods, which are based on the chemical evolution of our galaxy, that is how much the interstellar medium from which the star was formed was enriched in Y, Mg and Al. Tucci Maia et al. (2016) estimate the age to be 3.090 ± 0.391 Gyr based on an average age–[Y/Mg] relationship. Spina et al. (2018) found age estimate ranging from 3.2 ± 0.9 Gyr to 4.3 ± 0.5 Gyr depending on whether they use an age–[Y/Mg] or an age–[Y/Al] relationship and on the precise nature of their fit. Interestingly, Nissen et al. (2017) seem to confirm these relationships using the Kepler LEGACY database (Silva Aguirre et al. 2017).

Finally, some other studies focused on age–activity relationship to provide an estimate to t_⋆ for 18 Sco. Noteworthy are Mittag et al. (2016) and Lorenzo-Oliveira et al. (2018) which give respectively 5.1 ± 1.1 Gyr and 4.6 ± 0.9 Gyr. Similarly to the chemical-clocking estimates, the point estimates for the age are compatible within their error bars. However, the critical point is that those error bars differ significantly from the ones found using direct modelling and the Bayesian Statistics approach coupled to MCMC sampling. This raises the question of understanding how average-based estimates are representative of single objects. In other words how are these estimates affected by systematics. One also needs to understand if the current estimates for the age provided in this paper hint at the need to recalibrate these relations using more realistic uncertainties.

4. Conclusion

In this paper, we obtained estimates of the physical parameters of the solar twin 18 Sco using existing seismic data. Special care was taken to describe the asteroseismic diagnostic for the stars. A Bayesian Statistical model was used to relate the observations and the stellar parameters and statistical samples were obtained using an MCMC algorithm. A bimodal solution is obtained for the age, due to the difficulty to reproduce the seismic data. The most likely result gives an age that is roughly solar. This also points out the limitation of the current ground-based seismic data. This result may thus be used has benchmarks to evaluate in a near future the improvements made using the forthcoming TESS data or, potentially, SONG measurements.

Comparison of the resulting uncertainties with those obtained from local optimisation shows a discrepancy, the MCMC simulations leading to much larger uncertainties. The same conclusion applies when comparing these results to previous estimates in the literature. This motivates a more thorough investigation of the strategies used to estimate uncertainties on the physical characteristics, and in particular the age, of other solar twins.

Acknowledgments

MB would like to thank S. Hannestad for providing him access to the Grendel cluster at DCSC/AU of which important use has been made during this work. This material is based upon work supported by the NYU Abu Dhabi Institute under grant G1502. Part of this research was supported through the Laboratoire Lagrange BQR funding scheme. Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant DNRF106).

References

Appendix A: MCMC sampling

A.1. Algorithm

For the sake of completeness we describe here the MCMC algorithm used in this work. In order to explore efficiently the space of parameters, we combined a Simulated Annealing (SA; Liang et al. 2010) algorithm with an Adaptive Markov chain Monte Carlo algorithm with Global Adaptive Scaling (AMGAS; Andrieu & Thoms 2008, Algorithm 4). The pseudo-code is given in Algorithm A.1. In there, the quantity $q(.|{\theta }_m^{(n-1)})$ is the proposal density used by the MCMC algorithm to obtain a trial parameter, conditional on the current value of the Markov chain. Additionally, we ran several Markov chains in parallel.

The SA and AMGAS components of the algorithms are expected to improve the efficiency of the sampler. The former will improve convergence to the regions of highest probability density, even when the initial guess is chosen close to another local minimum. The latter improves the efficiency of the classical Metropolis–Hasting algorithm (Metropolis 1953; Hastings 1970) by improving iteratively the proposal density of the algorithm. The adaption of the scaling factor is related to the idea of of optimal scaling (Rosenthal 2008) and allows to adapt the scaling factor of the (Gaussian) proposal density so that it approaches a certain acceptance rate.

The main idea behind running multiple chains is to improve convergence diagnostics of MCMC algorithm (Gelman & Rubin 1992; Brooks & Gelman 1998). In principle, multiple chains do not sample the space of parameters better than one single long chain. In practice though, in the special case of stellar models, which take a long time to compute, it is beneficial to construct a sample from different chains, provided we can assess with some confidence that these have converged, since we are able to share the computational load on several processors. This helps the post-processing analysis of the posterior densities. In particular, the modelling of densities is improved when the size of the sample increases.

A.2. Convergence assessments

We show here some common convergence diagnostics for a multiple chain sampler. In Fig. A.1 we show the pooled and within variance and in Fig. A.2 the pseudo-scale reduction factor first introduced by Gelman & Rubin (1992). We use the proper correction for the degrees of freedom given in Brooks & Gelman (1998). We also show the posterior pooled and within variances. These three indicators taken together give a decent indication of convergence. First it appears that the two variances are fairly stabilised in all cases and converge towards each other (the within variance being lower as is expected). The pseudo-scale reduction factor is always lower than 1.2 which is a reasonable indicator for convergence (it is expected to converge to 1 when N → +∞).

These are indicators defined for univariate distributions. Brooks & Gelman (1998) also provide a convergence diagnostic for multivariate distributions. We also display it in Fig. A.2. It is an upper bound to the pseudo-scale reduction factor, as an MCMC algorithm converges slower towards the marginals than the joint posterior. Nevertheless it still decreases towards 1, giving another confirmation that our chains have converged towards the same stationary density.

Finally, we also used the cumulative mean of the sample of these chains to control the mixing. These indicate that the algorithm may not be performing optimally. This is further confirmed by the fact that the acceptance rates are usually lower than what is expected. They vary between roughly 3% and 10%. We attribute this to the strong correlations seen in the PDF.

Fig. A.1. Pooled (dashed lines) and within (full lines) variances for the stellar parameters.

Fig. A.2. Pseudo-scale reduction factor (full line) for all stellar parameters. The multivariate pseudo-scale reduction factor is shown in each panel as a black dashed line.

All Tables

All Figures

	Fig. 1. Examples of problematic PDFs of individual frequencies for n = 19, l = 0 (left panel) and n = 18, l = 2 (right panel). The black lines represent the histogram estimate and the red line a kernel-density estimation.
In the text

	Fig. 2. PDF for the small separation δν0, 19. The black lines represent the histogram estimation, the red line the kernel-density estimation and the green one the result of mixture model fitting.
In the text

	Fig. 3. Two- and one-dimensional marginal densities for the stellar parameters of 18 Sco for X = (Teff, L, [Fe/H],R, {δν}n, l), with the uniform prior on the mass given in Sect. 2.3.1. The red shades in the central panels represent the two-dimensional marginal densities (normalised). The side panels represent the one-dimensional ones.
In the text

	Fig. 4. Marginal joint probability density for (M, X0) (magnified from Fig. 3). The full lines mark locii of constant ZAMS luminosity. The dashed line shows the slope as obtained by Pearson’s correlation coefficient.
In the text

	Fig. 5. Examples of observed (full lines) and posterior (dashed lines) probability densities for the small separations δν15, 0 (left panel) and δν19, 0 (right panel).
In the text

	Fig. 6. Observed and estimated small separations. The dashed lines show the small separations of the best model with t⋆ <  5 Gyr. The full lines show the small separations of the best model with t⋆ >  5 Gyr. The red dots show the mean of the observed small separations and their standard deviations. The full dots mark small separations with l = 0 and the open dots those with l = 1.
In the text

	Fig. 7. Same as Fig. 3 but with X = (Teff, L, [Fe/H],[δν]n, l).
In the text

	Fig. 8. Same as Fig. 3 but with X = (Teff, L, [Fe/H],R) and X = (Teff, L, [Fe/H]) in the left and right panel, respectively.
In the text

Fig. 9. Upper row: marginal distributions (black) of M, t, and Z0 from the results presented in Fig. 3 (Sect. 3), and the distribution for these same parameters as determined from the local analysis (red Gaussians, Sect. 2.4, Table 5). Bottom row: marginal distributions for the same parameters, but restricted to a small interval around the values of α given in Table 5. The red lines are similar to those in the upper row.

In the text

	Fig. A.1. Pooled (dashed lines) and within (full lines) variances for the stellar parameters.
In the text

	Fig. A.2. Pseudo-scale reduction factor (full line) for all stellar parameters. The multivariate pseudo-scale reduction factor is shown in each panel as a black dashed line.
In the text