© ESO 2020

1. Introduction

In recent years, asteroseismology has become one of the most important tools to obtain precise estimates of stellar parameters. The first conclusive results for Sun-like stars were obtained from ground-based telescopes using high-precision spectrographs. The following bright stars were first observed: α Cen A (Bouchy & Carrier 2001; Bedding et al. 2001; Bazot et al. 2007), α Cen B (Kjeldsen et al. 2005), μ Ara (Bouchy et al. 2005), ι Hor (Vauclair et al. 2008), β Hyi (Bedding et al. 2007), and 18 Sco (Bazot et al. 2012). The space missions CoRoT (Baglin et al. 2009) and Kepler (Borucki et al. 2010) later provided much longer photometric time series for many more fainter stars. In particular, Kepler allowed for the study of an important number of Sun-like stars and the ability to obtain precise estimates of their masses, ages, and initial chemical compositions (Silva Aguirre et al. 2017).

Among the Sun-like stars observed by Kepler, the solar analogues 16 Cyg A and B stand out as particularly important. These are two of the brightest stars observed during the mission. Very precise measurements of their photometric flux were obtained over a period of 2.5 years. Spectral analysis of the resulting flux time series resulted in the detection of 54 and 56 global pulsation modes for 16 Cyg A and B, respectively (Metcalfe et al. 2012; Davies et al. 2015). For both stars, the uncertainties on the estimated eigenfrequencies are of the order of a few tenths of microhertz. Such precisions have allowed us to constrain stellar models tightly and obtain in turn a very good precision on the physical parameters of these stars, most notably their ages. Metcalfe et al. (2015) claim estimated uncertainties on the age of the order of 250 Myr. Modelled separately, 16 Cyg A and B appear to have ages 7.07 Gyr and 6.74 Gyr, respectively, favouring an age for the system around 7 Gyr. Based on these studies, the masses of these stars are expected to be slightly higher than the Sun, by a few percent, and their initial metallicities almost solar. These results confirmed previous spectroscopic studies that classified these stars, based on their atmospheric characteristics, as good solar analogues and sometimes even as solar twins (for a discussion, see for instance de Porto Mello et al. 2014).

It has often been argued that solar twins and analogues are good candidates to test non-standard stellar physics (see e.g. Bazot et al. 2018). This is because current stellar models have been mostly calibrated using the Sun, relying on the very good precisions of solar data. Consequently, they are usually well tested in regimes close to the solar regime, even for non-standard physics such as overshooting or microscopic diffusion (Basu 2016). The interest in using solar twins and/or analogues is thus that they can be modelled using the same physics without making further assumptions. This implies that we can confidently test these non-standard physics for stars that are not the Sun and therefore explore their behaviour when stellar physical characteristics, such as mass, age, or chemical composition, vary.

Another important aspect of studying 16 Cyg A and B is that the exquisite precision on their frequencies allow us to discuss the statistical reliability of the estimates of the physical parameters of 16 Cyg A and B. This is critical since observable stellar quantities, such as effective temperature, luminosity, seismic frequencies, radius, and metallicity, depend non-linearly on the stellar parameters of mass, age, initial chemical composition, and mixing-length parameter. Therefore good precisions on the data are needed to avoid complex behaviours of the underlying densities of the stellar parameters (e.g. Bazot et al. 2018). Even though some estimates of the stellar parameters may be found to vary significantly in some studies because of the variety of estimation strategies and numerical codes used for the modelling, most of these estimates show remarkable agreement. For instance, the estimated ages of 16 Cyg A span a range of ∼0.9 Gyr in Silva Aguirre et al. (2017), for an average value of the error bars of the order of 0.35 Gyr. In the framework of stellar-parameter estimation, the cases of 16 Cyg A and B the probability densities of the parameters are strongly constrained by the data. This allows for a careful discussion not only concerning the agreement of various estimates of their physical quantities, but also regarding the robustness of the associated uncertainties.

My first objective in this work is to evaluate the results published in previous studies that already show a good level of agreement (Metcalfe et al. 2012, 2015; Bellinger et al. 2016; Creevey et al. 2017). These works adopt various estimation strategies to obtain the stellar parameters. Metcalfe et al. (2015) and Creevey et al. (2017) consider the problem from a frequentist maximum-likelihood perspective, optimising their criteria using a genetic algorithm (Metcalfe et al. 2009). Bellinger et al. (2016) opts for a Bayesian approach using neural-network techniques to sample the probability densities of the stellar parameters. In this study, I also adopt the Bayesian approach and use Markov chain Monte Carlo (MCMC) sampling methods to approximate the probability density of the stellar parameters because the observations are fixed. This method has not been considered in previous studies.

The second objective is to explore the impact of some non-standard processes, namely microscopic diffusion and overshoot. The former is routinely used in solar models because it helps to reproduce properly the sound speed as measured by helioseismology (Basu 2016). In this work, the main concern is to compare the effect of including or not the diffusion of metals below the convective envelope.

The method itself and the algorithmic set-ups are described in Sect. 2. In Sect. 3, I present the results of the estimation problem and compare these with those of previous studies.

2. Method

2.1. Bayesian statistical model and sampling

This study was carried out within the Bayesian framework. It is particularly well suited to parameter estimation in the context of stellar physics. Indeed, contrary to laboratory experiments, for operational reasons it is usually not possible to repeat measurements of observable stellar quantities. The classical approach of frequentist statistics postulates that the data are random variables. These are described by probabilities distributions, which often depend on parameters. In these statistical models, any parameter is considered a deterministic quantity. Owing to the lack of repeated measurements, it is nevertheless difficult to estimate these parameters.

Bayesian statistics, thanks to the introduction of the prior distribution of the parameters, reverse this picture. The data become the fixed quantities and the parameters the random quantities. A significant practical advantage is that we can sample the space of parameters using various numerical strategies and, subsequently, use statistical tools to estimate the stellar parameters. This replaces advantageously the need for multiple measurements in classical frequentist statistics.

Bayesian statistics have been described in many monographs and review articles (see e.g. Berger & Berger 1985; Robert & Casella 2005; von Toussaint 2011), including in the context of stellar parameter estimation (see e.g. Bazot et al. 2008, 2012, 2016, and references therein). In this section, I simply recall the most important points of the approach and briefly describe the algorithmic strategy chosen for parameter estimation. Central to Bayesian statistics is the Bayes formula

$$\begin{array}{c}\hfill p(\mathit{\theta }|\mathit{x})\propto p(\mathit{x}|\mathit{\theta })p(\mathit{\theta }).\end{array}$$(1)

In this formula p(θ|x) is called the posterior density of θ conditional on x. This definition is perfectly general, but in the following I identify θ to the parameters of the model and x to the data. The likelihood is the quantity f(θ) = p(x|θ). Importantly, this is a probability density for the data (i.e. when p(x|θ) is seen as a function of x) but not for the parameters. It can be shown in frequentist statistics that a random experiment can be fully specified by the determination of a likelihood function (Birnbaum 1962; Robert & Casella 2005). Finally, p(θ) is the prior density on the parameters. This is the central concept of Bayesian statistics. The addition of a prior density is the key element that allows for shifting from the data to parameters as the random quantities. The prior density encodes the information one has on the parameters, θ, before their inference based on the measurements x.

Once the priors have been specified, we can make statistical statements based on the posterior density. The estimation of moments of p(θ|x) and credible intervals requires the integration of the density either over subsets of its domain. Therefore, I wish to obtain an approximation of p(θ|x), which does not have, in general, a closed form. In order to obtain it, I used an MCMC algorithm similar to that described in Bazot et al. (2019). The adaptive MCMC algorithm (Haario et al. 2001) was run independently on ten chains. Convergence of the MCMC simulations are discussed in Appendix A.

2.2. Data and likelihood

In order to set up a Bayesian statistical model, we need to define the likelihood. I first assume that the measurement errors are random and additive. This allows us to write

$$\begin{array}{c}\hfill \mathit{x}=\mathcal{S}({\mathit{\theta }}_{\mathbf{\star }})+\mathit{ϵ},\end{array}$$(2)

where 𝒮 is a theoretical model, which depends on the parameters θ_⋆; and ϵ is a random vector with zero mean, that is unbiased measurements are assumed, and the variance is determined using the observations. Or, depending on your meaning, “...mean, that is unbiased measurements are assumed and the variance is determined using the observations”.

For the sake of comparison, the non-seismic measurements (effective temperature, surface luminosity, surface metallicity and radius) I used are the same as those adopted in Metcalfe et al. (2015). Non-seismic constraints are listed in Table 1. Non-seismic measurements are assumed independent, therefore their likelihood is simply the product of the individual likelihoods. These are considered Gaussian, $\mathcal{N}({\mu }_k,{\sigma }_k^2)$, with μ_k the observed value of the observable k (for a given ordering of the non-seismic observations) and ${\sigma }_k^2$ its variance.

I chose not to use individual frequencies as a seismic diagnostic because of the need to estimate the surface effects that affect these frequencies (Kjeldsen et al. 2008), which are not included in the physical model 𝒮(θ_⋆). Taking these effects into account can be problematic (Bazot 2013) and I instead adopt the frequency ratios, r₀₁, r₀₂, and r₁₃ defined by Roxburgh & Vorontsov (2003). These can be straightforwardly computed from Table S2 of Davies et al. (2015). They are well understood theoretically and it has been emphasised, using empirical arguments, that they indeed seem largely independent of the surface effects (Silva Aguirre et al. 2013). I first assume that seismic and non-seismic data are uncorrelated. This means that the likelihood can be written as f(θ) = f(θ)_nsf(θ)_s; in this equation, and in the following, ns and s stand for non-seismic and seismic, respectively.

The assumption is made that the statistical noise on the measured frequencies is Gaussian and that their covariance matrix is diagonal. This is not the case for the individual ratios, which are in general correlated, since any given pair of separation ratios may involve one or several common eigenfrequencies. I assume that the seismic likelihood is a Gaussian random vector distributed as 𝒩(0, Σ). The covariance matrix Σ has to be evaluated from the variances of the individual frequencies. In this work, I use linear approximations as suggested by Roxburgh (2017). For given orderings of the individual frequency and the frequency ratios r₀₁, r₀₂, and r₁₃, the coefficient c_i, j of the covariance matrix is given by c_i, j = ∂_mx_s, ic_f; m, n∂_nx_s, j, where c_f; m, n is a coefficient of the covariance matrix, Σ_f, of the frequencies (for 16 Cyg A and B, c_f; m, n = 0 if m ≠ n), and the symbol ∂_m indicating derivation with respect to the mth frequency. I note that this covariance matrix could also be approximated using simulated eigenfrequencies sampled from 𝒩(0, Σ_f), simply computing the ratios for each realisation and then evaluating the covariance matrix Σ from the corresponding frequency-ratio sample. In general this method give results in fair agreement with the linear approximation, provided the number of realisations is large enough. I adopt the former method for the sake of simplicity and efficiency.

Given these considerations, the likelihood function in Eq. (1) can be written as

$$\begin{array}{cc}\hfill p(\mathit{x}|\mathit{\theta })\propto & exp(-\frac{1}{2}\sum _k\frac{{({\mathcal{S}}_{\mathrm{ns},k}({\mathit{\theta }}_{\mathbf{\star }})-x_{\mathrm{ns},k})}^2}{{\sigma }_k^2})\hfill \\ \hfill & \times exp(-\frac{1}{2}{({\mathit{x}}_{\mathrm{s}}-{\mathcal{S}}_{\mathrm{s}}({\mathit{\theta }}_{\mathbf{\star }}))}^T{\mathbf{\Sigma }}^{-1}({\mathit{x}}_{\mathrm{s}}-{\mathcal{S}}_{\mathrm{s}}({\mathit{\theta }}_{\mathbf{\star }}))).\hfill \end{array}$$(3)

The likelihood function departs from that considered in Metcalfe et al. (2015) in that the authors did not include the non-diagonal terms in Σ (Metcalfe, priv. comm.). This likelihood function also differs from Creevey et al. (2017) since the authors only considered non-diagonal terms in Σ for pairs of r₀₁ individual ratios, but not for pairs of r₀₂ ratios or r₀₁/r₀₂ pairs. Furthermore, Metcalfe et al. (2015) use r₀₁₀ ratios instead of r₀₁ ratios, the former being obtained by the addition of the r₁₀ ratios (Roxburgh & Vorontsov 2003) to the latter. This procedure has been criticised by Roxburgh (2018) on the grounds that in an r₀₁₀ ratio 2N quantities are fitted, while stellar models provide only N frequency phase corresponding to both the r₀₁ and r₁₀ sequences. This may lead to an over-fitting configuration, possibly introducing biases into the estimation process.

2.3. Physical models and priors

In Eq. (1), the parameters may encapsulate any relevant quantity of the Bayesian statistical model that ought to be estimated. In practice, all the parameters that do not enter as an argument of 𝒮 are fixed. This includes the μ_ks, the σ_k, and the coefficients of Σ, and therefore θ = θ_⋆.

The physical model I adopt is a spherically symmetric, non-rotating, non-magnetic star. The corresponding equations for the stellar structure and evolution are solved numerically via ASTEC (Christensen-Dalsgaard 2008a) and those for stellar pulsations via adipls (Christensen-Dalsgaard 2008b). The equation of state was set according to the OPAL prescription (Rogers & Nayfonov 2002). The opacities are also provided by the OPAL collaboration (Iglesias & Rogers 1996). Convection is treated following the mixing-length formalism of Böhm-Vitense (1958). Nuclear reaction rates are from the NACRE collaboration (Angulo et al. 1999) and supplemented by the values given in Angulo et al. (2005) for the ¹⁴N(p, γ)¹⁵O reaction.

Other physical prescriptions were varied to test their impact on parameter estimation. Microscopic diffusion of elements heavier than hydrogen is treated using the formalism of Michaud & Proffitt (1993). I either considered simple helium diffusion or helium and metal diffusion. The ratios of chemical element abundances are taken either from Grevesse & Noels (1993) or Grevesse & Sauval (1998). For the sake of simplicity and comparison with Metcalfe et al. (2015) and Creevey et al. (2017), I disregarded the abundance ratios provided by Asplund et al. (2009), which have raised many difficulties in the solar case (see e.g. Guzik et al. 2006; Castro et al. 2007; Basu & Antia 2008; Antia & Basu 2011; Gough 2012; Basu 2016). This case is left to future studies. Finally, I considered the possibility of overshoot below the convective envelope, that is transport process beyond the limit point predicted by a linear stability analysis.

Solving the equations for stellar structure and evolution under the above assumptions demand to set five free parameters, plus one in case overshoot is included. The basic parameters are the stellar mass, M_⋆, age, t_⋆, initial chemical composition, given by the initial hydrogen-mass fraction, X₀, and metallicity, Z₀, and the mixing-length parameter, α. This latter is a proportionality coefficient between the mean-free path of a fluid parcel in the convective zone and the pressure scale height. The mixing-length parameter sets the depth of the convective zone. When it is included, a new parameter α_ov must also be taken into account, whose signification is discussed in Sect. 3.2.

The last necessary step in the construction of the Bayesian statistical model is to specify the priors on the stellar parameters. Following Bazot et al. (2018) I assume uniform priors on all parameters. Their boundaries are given in Table 2.

3. Results and discussion

The results of the MCMC simulations are given in Tables 3 and 4 for 16 Cyg A and B, respectively. For both stars I used r₀₁, r₀₂, and r₁₃ as seismic diagnostics for the first five runs. For the sixth run, r₀₁ was replaced with r₀₁₀ to test the potential effects of over-fitting. For each simulation I provide three different estimates of the stellar parameters parameters M_⋆, t_⋆, X₀, Z₀, α and, when relevant, α_ov. I also provide estimates for the initial helium-mass fraction, Y₀, to facilitate comparison with other studies, this parameter often being reported instead of X₀. For a given set-up, the first line of estimates correspond to the maxima of the marginal densities. These maxima are given alongside credible intervals that are defined as the smallest interval with probability 0.683 that contains the maximum. The second line is the maximum a posteriori (MAP) obtained from the joint distribution p(θ_⋆|x). I do not provide uncertainties for this point estimate. In the third the posterior means (PM) and posterior standard deviations (PSD) are reported; these estimates are computed from the marginal densities.

The baseline case was chosen so that it corresponds to the set-up used in Metcalfe et al. (2015), who also used ASTEC¹. This set-up uses the Grevesse & Noels (1993) abundances ratios. Only helium is included in the diffusion equation of chemical elements. Overshoot is not taken into account.

The observed seismic diagnostics are indicated in Fig. 1. I also show the best model for each case in Tables 3 and 4 (red thin lines). These optimal theoretical values are barely distinguishable from one model to the other. This indicates that they may be equally valid in order to reproduce the data, and that statistical model comparison methods (see e.g. Robert 2007, Chap. 7) alone will not suffice to distinguish them. Further discussions on model comparison are given in Sects. 3.2 and 3.1. Visual inspection of Fig. 1 also indicates that the accuracy of the present results is similar to what was obtained by Metcalfe et al. (2015). Figure 2 shows the MAP estimates of all non-seismic observations. They are all reproduced with good accuracy, within 1σ of the observed values.

Fig. 1. Seismic diagnostic for 16 Cyg A (left panel) and B (right panel). The squares, triangles, and circles denote the observed r13, r02, and r01 observed ratios, respectively. The red lines show the best models for all runs in Tables 3 and 4.

Fig. 2. Non-seismic observations for 16 Cyg A (left panel) and B (right panel). The horizontal lines represent the observed values (full lines) and the corresponding uncertainties (dashed lines). The dots denote the estimated values (MAPs of the marginal densities) for the corresponding observation, for each run. The runs are labelled on the abscissa according to the numbering given in Tables 3 and 4.

The two-dimensional joint PDFs for this set-up are shown in Fig. 3 alongside the corresponding correlation coefficients. The mass correlates most significantly with X₀. The age correlates extremely well with the mixing-length parameter and slightly less with X₀. It may be counter-intuitive that the mass does not correlate strongly with the age. A likely explanation is that, due to the level of constraint imposed by the seismic data, these two parameters cannot vary enough to allow such a correlation to be observed. Interpreting these correlations is not straightforward because there is in general no one-to-one correspondence between a parameter and an observation. Thus elements of interpretation can be gathered by looking at the correlation between the stellar parameters and other quantities. The mass and X₀ (Pearson correlation coefficients ≳0.8) both correlate strongly with the luminosity and the radius on the zero-age main sequence (ZAMS). This is related to the mass-luminosity relation on the ZAMS, in which the mean-molecular weight enters (see e.g. Clayton 1968). To maintain a roughly constant luminosity, an increase in mass must be compensated by a decrease in the mean-molecular weight. This latter corresponds to an increase in the initial hydrogen-mass fraction and a decrease of the initial helium-mass fraction because the two quantities are perfectly anti-correlated. The interaction between the mass and the hydrogen-mass fraction thus sets the initial conditions of the stellar evolution sequence. Such a behaviour was already noted in Bazot et al. (2018), even though the trend is not as clear in this work.

Fig. 3. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg A (left) and B (right). The central panels show the joint marginal densities of the paired parameters. In the side panels are plotted the individual marginal densities. The red shaded areas in the central panels and the full lines in the side panels represent posterior densities for the baseline case, that is without diffusion of metals. The black contours in the central panels and the dashed lines in the side panels represent posterior densities for models with diffusion of metals. The numbers in each panel give the Pearson correlation coefficient for the two variables. The values in red correspond to the red shaded densities and the values in black to the black contours.

The mixing-length parameter correlates very strongly with the effective temperature. It is well known that these two quantities are related Clayton (1968). Again, the mixing-length parameter allows for setting the initial conditions of the evolutionary sequence. The age also anti-correlates well with the effective temperature (Pearson correlation coefficient ∼0.8), more than with any other observable. This is an evolutionary effect. Those two correlations explain the relation between the age and the mixing length.

Another very strong correlation of note is that between the initial metallicity and the surface metallicity-to-hydrogen ratio. It is much stronger (Pearson correlation coefficient >0.9) than the correlation between Z/X and X₀. This is because too large variations of X₀ would also affect the radius and the luminosity, whereas variations of Z₀ do not have the same effect. These are very general observations. In order to better understand the behaviour of these correlations with the observational constraints, it would be necessary to use simulated data in the spirit of Brown et al. (1994) or Creevey et al. (2007). In Fig. 3 also shows one-dimensional marginal densities, which are all Gaussian to a very good approximation. This means that their PSD can be interpreted as a 0.683 credible interval as defined above, the PM and maxima of the marginal being equal.

As is shown below, the stellar mass is of importance as it is a good marker of the robustness of the asteroseismic (and interferometric) constraints. For the baseline case I obtain 1.07 ± 0.02 M_⊙ and 1.05 ± 0.02 M_⊙ for 16 Cyg A and B. These are fairly close to the estimates of Metcalfe et al. (2015) and Creevey et al. (2017), both of which use ASTEC, albeit with slightly different set-up for the latter. The same is true for the other parameters with the possible exception of the age of 16 Cyg A quoted in Creevey et al. (2017), which is higher by 6.5% (i.e. roughly 2.4σ away from the PM value). The high precisions on all parameters can be traced back to the use of asteroseismology. This can be seen qualitatively by comparing the current results to those of Bazot et al. (2018) on 18 Sco. This latter star, a solar twin, has even more precise spectrophotometric and interferometric data than 16 Cyg A and B. However, the precision achieved on the estimates of its parameter using much poorer asteroseismic data is far worse. For the age, there is an order of magnitude difference in precision.

A difficulty arises when comparing the uncertainties on these parameters. Indeed both Metcalfe et al. (2015) and Creevey et al. (2017) use samples obtained from a genetic algorithm to estimate confidence intervals (Metcalfe et al. 2014). In the case of MCMC algorithms, asymptotic properties of Markov chains are used to show that the resulting sample is generated according to the target distribution, i.e. p(θ_⋆|x) in this work. The issue with genetic algorithms is that there exist no such theorem guaranteeing the convergence of the resulting sample to the joint probability of the parameters being optimised². Therefore, and even though the quoted uncertainties in these previous studies are often in good numerical agreement with those given in this work, we cannot draw strong conclusions by simply comparing their estimates with those of Tables 3 and 4.

Other studies have provided estimates for the stellar parameters of 16 Cyg A and B. In Silva Aguirre et al. (2017) results are reported from multiple groups using various codes and estimation strategies. The values quoted in this study bracket those represented in Tables 3 and 4. Finally, other works to be mentioned are those of Bellinger et al. (2016, 2017). In the first paper, the author used neural-network strategies to estimate stellar parameters in a Bayesian fashion. Their estimates are relatively close to those given in this study. In particular, the age estimates are exactly the same for both 16 Cyg A and B. An exception is the mass of 16 Cyg B, which is slightly lower than that seen in Table 4. As discussed in Sect. 3.3, this can potentially be explained by the treatment of the seismic constraints. Their estimated uncertainty on the age of 16 Cyg A is twice as large as those given in Table 4. Such a discrepancy is not straightforward to explain. A possible factor is the inclusion in their model of a parameter that controls the efficiency of diffusion. However, we would expect 16 Cyg B to be also be affected to some extent, which does not seem to be the case. This problem will require further work. Finally, their distribution for the initial helium-mass fraction of 16 Cyg A shows a bimodality that I do not observe. This stresses that, despite a satisfying general agreement between the two methods, the posterior densities used to derive the stellar parameters still show discrepancies that ought to be explained.

The estimates of the age for both stars are in extremely good agreement, 6.93 ± 0.19 Gyr and 6.83 ± 0.17 Gyr. Confirming previous studies, given the independent modelling of 16 Cyg A and B this corresponds to an age of 6.88 ± 0.12 Gyr for the system. This value is close to that claimed by Metcalfe et al. (2015), roughly departing by 1σ.

3.1. Chemical abundances and microscopic diffusion

The treatment of chemical abundances may have a strong impact on the physical characteristics of the stars. In this section, I study the potential biases induced by different prescriptions for the abundances ratios and treatments of microscopic diffusion. The results are listed in Tables 3 and 4. The first interesting result is that the mass estimate is remarkably stable with respect to such changes. This reflects the fact that the radius and density, through the seismic data, are well constrained, hence providing a robust mass estimate (see for instance Bazot et al. 2012, 2018 for another example of robust mass estimate using interferometry and asteroseismology; see also Creevey et al. 2007 and Cunha et al. 2007 for more general discussions on the interplay between asteroseismology and interferometry).

The main impact of the change from Grevesse & Noels (1993) to Grevesse & Sauval (1998) abundance ratios is a small decrease in X₀ and t_⋆, coupled to a small increase of Y₀ (all below the 1σ uncertainties from the baseline case). The estimate of Z₀ remains stable. Overall, the biases to be expected from such a change should remain small. This may not be the case if I had considered Asplund et al. (2005) abundances. These are known to produce discrepancies between solar models and helioseismic observations (see Basu 2016, for a review). Because of these discrepancies I did not include these abundances in this study; the models used in this work were satisfyingly validated in the solar case for the Grevesse & Noels (1993) and Grevesse & Sauval (1998) abundance ratios.

Microscopic diffusion of chemical abundances plays an important role in the characterisation of the solar sound speed profile, which in turn relates to the oscillation frequencies (Bahcall et al. 1995; Christensen-Dalsgaard et al. 1996; Richard et al. 1996). Evidence of gravitational settling of heavy elements has also been found in more massive stars (e.g Richard et al. 2001). The Kepler data may help to understand the magnitude of the bias induced by neglecting diffusion for all or some elements. In theory, if we consider particles settling in an hydrogen background, then all heavier elements should be included. However, treatments of diffusion in the literature are somewhat inconsistent. First of all, diffusion has not been systematically included. This has the advantage of helping to set the initial hydrogen-mass fraction and metallicity, at least for stars up until the first dredge-up, since their surface ratio does not evolve and correspond to the observed metallicity. Diffusion is sometimes also excluded from stellar models due to numerical issues. For ASTEC difficulties may arise when treating diffusion of elements heavier than helium when the star has a convective core (Christensen-Dalsgaard 2008a). This is why only helium diffusion was considered in Metcalfe et al. (2015) and Creevey et al. (2017). Sometimes diffusion is neglected above a certain stellar mass (Silva Aguirre et al. 2017; Creevey et al. 2017). The argument given is that the short diffusive timescales corresponding to shallow convective envelopes would deplete the upper layers of elements heavier than hydrogen too fast. This actually reflects the incompleteness of most stellar models that do not include “radiative levitation”, which becomes an efficient competing phenomena at intermediate masses (e.g. Richard et al. 2002). 16 Cyg A and B present the advantage of having low enough mass (even when accounting for the uncertainties on this parameter) so that no convective core ever develops (this is also due to their almost solar metallicity; see e.g. Bazot et al. 2012, 2016).

I ran two additional simulations for each star: the first does not include diffusion and the second takes into account settling of helium and heavy elements. The resulting two-dimensional and one-dimensional posterior marginal densities are shown in Fig. 3. It should first be noticed that this does not change the magnitude of the estimated uncertainties for any parameter. When diffusion is not included, the estimated age is higher than in the baseline case and the mixing-length parameter is lower. Within their 68.3% credible intervals, all estimates agree albeit marginally. The estimated X₀ and Z₀ are lower and higher, respectively, but they remain closer to those of the baseline case. When diffusion of metal is included, the trends for X₀ and Z₀ are similar, but deviate more from the baseline case. The estimate of X₀, which now decreases by more than 1σ and Z₀ increases by ∼15% (1.8σ from baseline case) for 16 Cyg A and ∼11% (1.3σ from baseline case). This is a well-known effect resulting from the depletion of the external convective envelope of its metals. In order to reproduce the observed surface metallicity, the average value Z₀ needs to increase. However, contrary to the non-diffusion case, the estimated age this time decreases and the mixing-length parameter does not change.

This gives us some insight into potential biases caused by the details of microscopic diffusion in stellar models. To that effect, I now consider that the most accurate physical model is that including microscopic diffusion for the metals. Then, provided that the asteroseismic data are similar in quality to the Kepler time series for 16 Cyg A and B and that the star is in a physical state close enough to the solar state, we could expect that neglecting diffusion implies higher estimated ages and lower mixing-length parameters. The initial metallicity and hydrogen-mass fraction may also be marginally affected, respectively, decreasing and increasing. On the other hand, including helium diffusion does not so strongly impact the estimates for the age and the mixing-length parameter, but induces higher deviation in the initial chemical composition. The strongest bias is on the age when diffusion is omitted. To conclude, it should be noted that even moderate bias may become important when attempting to perform statistical studies such as those seen in stellar open cluster analysis or Galactic archaeology. This may cause the averages used for the quantities of interest to converged towards a wrong value.

3.2. Overshooting

Overshooting is a long-standing problem in stellar physics. Stellar models use mixing-length theory to model convection. It is a local theory in which convective motions are expected to stop at the boundary of the unstable envelope defined by the Schwarzschild criteria, that is when the radiative temperature gradient equals the adiabatic temperature gradient. It has long been recognised that convective movements are likely to persist beyond this limit (Veronis 1963; Moore et al. 1967). It is not clear however whether downward convective flows are strong enough to thermalise the subadiabatically stratified layers below the convective boundary. This is the case when the Péclet is high enough; such a situation is often referred to as penetration (Zahn 1991). This is opposed to chemical mixing, when overshooting only homogenises abundances without affecting the temperature gradient. In this study I follow the recommendation of Viallet et al. (2015), based on a semi-quantitative analysis, which suggests using penetration to model overshoot below the convective envelope. The depth of the overshoot/penetration zone is set using

$$\begin{array}{c}\hfill {\ell }_{\mathrm{ov}}={\alpha }_{\mathrm{ov}}{H}_{\mathrm{p}},\end{array}$$(4)

where ℓ_ov is the size of the overshoot layer and H_p the pressure scale height at the boundary of the convective envelope.

There have been a few estimates of the overshoot parameter for 16 Cyg A and B. In Silva Aguirre et al. (2017), six stellar evolution codes out of seven include overshoot. However, the overshoot parameter has been fixed to calibrated values. A notable exception to this lack of published estimates is the work of Bellinger et al. (2016) who obtained values for the overshoot parameter of 0.07 ± 0.03 and 0.11 ± 0.03 for 16 Cyg A and B, respectively. Their prescription for overshooting differs from the present one in that chemical mixing (and not penetration) is achieved through a diffusive process (Herwig 2000).

I ran two additional simulations for each star to estimate α_ov. In the first simulation only helium settles. In the second, diffusion of metals is taken into account. The idea is that the diffusion of heavy elements affects the layers immediately below the convective zone and that it may impact overshoot, which is expected to occur in the same region. The corresponding PDFs of the stellar parameters are shown in Fig. 4 and their estimates are given in Tables 3 and 4. The first obvious result is that α_ov is poorly constrained. The marginal densities for α_ov, contrary to the other parameters, are not Gaussian, but closer to a gamma distribution. From the marginal MAP estimates I obtain uncertainties above 100%. This is in sharp contrast with the results from Bellinger et al. (2016) who estimated uncertainties of the order of 40%. It is not clear if this discrepancy is due to the difference in the estimation strategy or in the physical modelling. It is also noteworthy that α_ov does not correlate with any other stellar parameters. There is therefore no bias to be found in the other stellar parameters for the inclusion of penetrative overshoot. We see that α_ov is marginally larger than when only helium is diffusing. However, the uncertainties remain extremely large and the two prescriptions cannot be distinguished.

Fig. 4. Marginal densities for the stellar parameters M, t⋆, X0, Z0, α, and αov of 16 Cyg A (left) and B (right). The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels. The red shaded areas in the central panels and the full lines in the side panels represent posterior densities for models with penetrative overshoot and helium diffusion (#5 in Tables 3 and 4). The black contours in the central panels and the dashed lines in the side panels represent posterior densities for models with penetrative and metal diffusion (#6 in Tables 3 and 4). The numbers in each panel give the Pearson correlation coefficient for the two variables. The values in red correspond to the red shaded densities and the values in black to the black contours.

In order to estimate α_ov for both overshoot and penetration, I ran an additional simulation for each star. Besides the inclusion of penetrative convection, the physical models are identical to the baseline case. The corresponding PDFs of the stellar parameters are shown in Fig. 4 and their estimates are given in Tables 3 and 4. The first obvious result is that α_ov is poorly constrained. The marginal densities for α_ov, contrary to the other parameters, are not Gaussian, but closer to a gamma distribution. From the marginal MAP estimates I obtain uncertainties above 100%. This is in sharp contrast with the results from Bellinger et al. (2016) who estimated uncertainties of the order of 40%. It is not clear if this discrepancy is due to the difference in the estimation strategy or in the physical modelling. It is also noteworthy that α_ov does not correlate with any other stellar parameters. There is therefore no bias to be found in the other stellar parameters for the inclusion of penetrative overshoot.

The extent of the overshotting region in the Sun has been measured by several authors. For instance Monteiro et al. (1994) and Christensen-Dalsgaard et al. (1995) have used solar models including penetration to reproduce the signatures that sharp transitions in the sound speed profile induce in the oscillation frequencies. These authors estimate a depth for the solar overshoot region in an approximate range 0.07H_p–0.1H_p. When helium diffusion only is included, the upper limits of the 68.3% credible intervals for α_ov are below these solar values for both 16 Cyg and B. When metal diffusion is included, the 68.3% credible intervals for α_ov overlaps marginally with the range given for the solar value.

Numerical investigations have also been used to evaluate the characteristics of stellar overshooting. They are obviously limited in that they they assume Prandtl numbers that are too high and Rayleigh numbers that are too low with respect to a typical stellar interior. These numerical invesigations nevertheless provide interesting insight into the general behaviour of stellar convection. Many studies have been published considering either fluids in the Boussinesq approximation, the anelastic approximation, or fully compressible fluids. Broadly speaking two types of results emerge depending on whether they are the models are two- or three-dimensional. Two-dimensional simulations (Hurlburt et al. 1986, 1994; Rogers & Glatzmaier 2005; Rogers et al. 2006) often exhibit penetrative overshoot. On the other hand, in three-dimensional simulations (Singh et al. 1995; Brummell et al. 2002; Korre et al. 2019) penetration is not observed but rather chemical mixing occurs down to significant depth. This is most likely due to a larger density of plume-like structure in two-dimensional than in three-dimensional simulations. Our results, even though they do not address chemical mixing, are consistent with the picture of a very small to non-existent penetration region.

3.3. Using r₀₁/r₁₀ ratios

The final MCMC runs for both 16 Cyg A and B aim at assessing the effect of using r₀₁₀ instead of r₀₁. The obvious significant result is that the estimates of the masses depart from those of all other runs. The discrepancy is significant, i.e. of the order of 3% for both stars, which corresponds to ∼1.5σ from the baseline case. The corresponding posterior densities are shown in Fig. 5. In light of the criticism expressed in Roxburgh (2018), this discrepancy can be interpreted as a bias with respect to the “correct” value of the parameter, which according to the previous discussions are likely to be around 1.07 M_⊙ and 1.05 M_⊙ for 16 Cyg A and B, respectively. The use of r₀₁₀ also biases the results towards lower X₀ and α values and, marginally, towards higher values of t_⋆. The uncertainties on the age are also biased when using r₀₁₀ as seismic constraints. The estimates in this case are twice as small as those obtained using r₀₁.

Fig. 5. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg A (left) and B (right). The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels. The red shaded areas in the central panels and the full lines in the side panels represent posterior densities for the baseline case. The black contours in the central panels and the dashed lines in the side panels represent posterior densities when r010 is used as a seismic constraint instead of r01.

The effect uncovered in this work happens to be quite subtle. Indeed, it has already been pointed out that the results for the baseline case are in good agreement with those of Metcalfe et al. (2015) who used r₀₁₀ and not r₀₁ as done in this case. This can be explained by the fact that they neglected the non-diagonal terms in Σ. To test this interpretation, I use MCMC simulations for both 16 Cyg A an B using r₀₁₀ but computing the likelihood (3) using only the diagonal terms of Σ. I then find mass estimates in agreement with the baseline case.

The general picture that emerges from these results is that decent models for 16 Cyg A and B may be found at masses around 1.04 M_⊙ and 1.02 M_⊙, respectively. Their emergence as local or global minima strongly depends on the precise computation of the likelihood, in particular on the non-diagonal terms of Σ. A good indicator of this behaviour is the value the argument of the exponential in Eq. (3), that is the χ² values, of the best models for the baseline case and the case constrained using r₀₁₀. In the former case, it is indeed larger for the model at 1.07 M_⊙, while in the latter it is larger for the 1.04 M_⊙ model. This behaviour is reversed again if the non-diagonal terms in Σ are set to zero. This shows how delicate this problem can be. There are competing effects and, in the case of 16 Cyg A and B, the bias induced by over-fitting is compensated by the approximation that consists in neglecting correlations between frequency ratios. However, there is no indication that this behaviour can be straightforwardly generalised to other stars. Therefore, the present results strengthen the more general claim of Roxburgh (2018), outlining the need to properly take into account the correlations between seismic indicator and using r₀₁ (or r₁₀) rather than r₀₁₀ ratios.

4. Conclusions

In this study, I used Bayesian statistics to derive robust estimates and uncertainties for the physical parameters of the solar analogues 16 Cyg A and B. Using a statistical method independent from the previous studies on this stars, I obtain a precision on the mass of the order of 4% and on the age of the order of 6%. This is in fair agreement with most of the published estimates. I outline the need to use the proper seismic diagnostics to avoid biases on the mass, initial chemical composition, and mixing-length parameter (and, marginally, the age), and their uncertainties. I also pointed out the changes induced in the estimates of the initial chemical composition induced by changes in the abundances ratios and microscopic diffusion. The latter, limiting the comparison to the abundance ratio of Grevesse & Noels (1993) and Grevesse & Sauval (1998), are negligible. The effects of microscopic diffusion are much more important. They may cause biases in the estimates in the range 7–8%. They may also introduce some biases in the estimated initial metallicity. Considering penetrative overshoot below the convective envelope, I find that the sub-adiabatic region can be thermalised only in a very shallow layer. Including penetrative overshoot does not affect the estimates of the other parameters. The results presented in this study complete those previously obtained on 16 Cyg A and B, which are excellent benchmark of what can be achieved using precise asteroseismic data. The future PLATO mission will allow us to observe many more Sun-like stars. The present results aim at helping to prepare their modelling.

Acknowledgments

I thank the referee for his/her helpful comments, which really helped improving this article. I thank S. Hannestad for providing him access to the Grendel cluster at DCSC/AU of which important use has been made during this work. I also thank T. Metcalfe for answering many questions on his previous work, this was extremely helpful in constructing this analysis. This material is based upon work partly supported by the NYU Abu Dhabi Institute under grant G1502.

References

Appendix A: Convergence of the MCMC algorithm

In this appendix, I provide assessments on the convergence of the MCMC simulations. In Figs. A.1 and A.2, I show the pseudo-scale (Gelman & Rubin 1992) and multivariate pseudo-scale (Brooks & Gelman 1998) criteria. The general rule-of-thumb regarding convergence is to have these criteria below ∼1.2. In the limit of infinitely many iterations, they can be expected to tend to one. All the runs are well behaved with respect to these convergence diagnostics. The multivariate pseudo-scale factor is always superior to pseudo-scale factor, which is expected, the latter converging to one faster than the former.

In theory, a single chain run for a very long time, as well as several shorter chains, would sample the target density. However, given the large computational cost of stellar models, the current set-up represents a real gain in that samples generated from the same stationary distribution can be merged to get larger samples, Furthermore, such tests as those performed in this work are not available with one single chain, hence representing an operational advantage.

Fig. A.1. Pseudo-scale (R, full lines) and multivariate pseudo-scale factors (Rp, dashed line) for all cases presented in Table 3 16 Cyg A. The dashed horizontal red lines denote the 1 and 1.2 values.

Fig. A.2. Pseudo-scale (R, full lines) and multivariate pseudo-scale factors (Rp, dashed line) for all cases presented in Table 4 16 Cyg B. The dashed horizontal red lines indicate the 1 and 1.2 values.

All Tables

All Figures

	Fig. 1. Seismic diagnostic for 16 Cyg A (left panel) and B (right panel). The squares, triangles, and circles denote the observed r13, r02, and r01 observed ratios, respectively. The red lines show the best models for all runs in Tables 3 and 4.
In the text

Fig. 2. Non-seismic observations for 16 Cyg A (left panel) and B (right panel). The horizontal lines represent the observed values (full lines) and the corresponding uncertainties (dashed lines). The dots denote the estimated values (MAPs of the marginal densities) for the corresponding observation, for each run. The runs are labelled on the abscissa according to the numbering given in Tables 3 and 4.

In the text

Fig. 3. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg A (left) and B (right). The central panels show the joint marginal densities of the paired parameters. In the side panels are plotted the individual marginal densities. The red shaded areas in the central panels and the full lines in the side panels represent posterior densities for the baseline case, that is without diffusion of metals. The black contours in the central panels and the dashed lines in the side panels represent posterior densities for models with diffusion of metals. The numbers in each panel give the Pearson correlation coefficient for the two variables. The values in red correspond to the red shaded densities and the values in black to the black contours.

In the text

Fig. 4. Marginal densities for the stellar parameters M, t⋆, X0, Z0, α, and αov of 16 Cyg A (left) and B (right). The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels. The red shaded areas in the central panels and the full lines in the side panels represent posterior densities for models with penetrative overshoot and helium diffusion (#5 in Tables 3 and 4). The black contours in the central panels and the dashed lines in the side panels represent posterior densities for models with penetrative and metal diffusion (#6 in Tables 3 and 4). The numbers in each panel give the Pearson correlation coefficient for the two variables. The values in red correspond to the red shaded densities and the values in black to the black contours.

In the text

Fig. 5. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg A (left) and B (right). The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels. The red shaded areas in the central panels and the full lines in the side panels represent posterior densities for the baseline case. The black contours in the central panels and the dashed lines in the side panels represent posterior densities when r010 is used as a seismic constraint instead of r01.

In the text

	Fig. A.1. Pseudo-scale (R, full lines) and multivariate pseudo-scale factors (Rp, dashed line) for all cases presented in Table 3 16 Cyg A. The dashed horizontal red lines denote the 1 and 1.2 values.
In the text

	Fig. A.2. Pseudo-scale (R, full lines) and multivariate pseudo-scale factors (Rp, dashed line) for all cases presented in Table 4 16 Cyg B. The dashed horizontal red lines indicate the 1 and 1.2 values.
In the text