Issue 
A&A
Volume 635, March 2020



Article Number  A26  
Number of page(s)  12  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201935565  
Published online  02 March 2020 
Uncertainties and biases in modelling 16 Cygni A and B
^{1}
Division of Sciences, New York University, Abu Dhabi, UAE
email: mb6215@nyu.edu
^{2}
Center for Space Science, NYUAD Institute, Abu Dhabi, UAE
Received:
28
March
2019
Accepted:
23
December
2019
Aims. In this study I assess how existing data for the solar analogues 16 Cyg A and B, in particular the asteroseismic measurements obtained from Kepler, constrain theoretical stellar models. The goal is twofold: first to use these stars as benchmarks to discuss which precisions can realistically be expected on the inferred stellar quantities; and second to determine how well “nonstandard” prescriptions, such as microscopic diffusion and overshoot, are constrained.
Methods. I used a Bayesian statistical model to infer the values of the stellar parameters of 16 Cyg A and B. I sampled the posterior density of the stellar parameters via a Markov chain Monte Carlo algorithm, tested different physical prescription, and examined the impact of using different seismic diagnostics.
Results. General good agreement is found with several recent modelling studies on these stars, even though some discrepancies subsist regarding the precise estimates of the uncertainties on the parameters. An age of 6.88 ± 0.12 Gyr is estimated for the binary system. The inferred masses, 1.07 ± 0.02 M_{⊙} for Cyg A and 1.05 ± 0.02 M_{⊙} for Cyg B, are shown to be stable with respect to changes in the physical prescriptions considered for the modelling. For both stars, microscopic diffusion has a significant effect on the estimates of the initial metallicity. Overshoot is confined to very small regions below the convective zone. I show that a proper treatment of the seismic constraints is necessary to avoid biases in the estimate of the mass.
Key words: stars: individual: 16 Cyg A / stars: individual: 16 Cyg B / stars: oscillations / stars: solartype / asteroseismology / methods: statistical
© 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 Sunlike stars were obtained from groundbased telescopes using highprecision 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 Sunlike stars and the ability to obtain precise estimates of their masses, ages, and initial chemical compositions (Silva Aguirre et al. 2017).
Among the Sunlike 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 nonstandard 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 nonstandard 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 nonstandard 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 nonlinearly on the stellar parameters of mass, age, initial chemical composition, and mixinglength 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 stellarparameter 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 maximumlikelihood perspective, optimising their criteria using a genetic algorithm (Metcalfe et al. 2009). Bellinger et al. (2016) opts for a Bayesian approach using neuralnetwork 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 nonstandard 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 setups 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
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
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 nonseismic measurements (effective temperature, surface luminosity, surface metallicity and radius) I used are the same as those adopted in Metcalfe et al. (2015). Nonseismic constraints are listed in Table 1. Nonseismic measurements are assumed independent, therefore their likelihood is simply the product of the individual likelihoods. These are considered Gaussian, , with μ_{k} the observed value of the observable k (for a given ordering of the nonseismic observations) and its variance.
Nonseismic observational properties of 16 Cyg A and B.
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_{01}, r_{02}, and r_{13} 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 nonseismic data are uncorrelated. This means that the likelihood can be written as f(θ) = f(θ)_{ns}f(θ)_{s}; in this equation, and in the following, ns and s stand for nonseismic 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_{01}, r_{02}, and r_{13}, the coefficient c_{i, j} of the covariance matrix is given by c_{i, j} = ∂_{m}x_{s, i}c_{f; m, n}∂_{n}x_{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 frequencyratio 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
The likelihood function departs from that considered in Metcalfe et al. (2015) in that the authors did not include the nondiagonal terms in Σ (Metcalfe, priv. comm.). This likelihood function also differs from Creevey et al. (2017) since the authors only considered nondiagonal terms in Σ for pairs of r_{01} individual ratios, but not for pairs of r_{02} ratios or r_{01}/r_{02} pairs. Furthermore, Metcalfe et al. (2015) use r_{010} ratios instead of r_{01} ratios, the former being obtained by the addition of the r_{10} ratios (Roxburgh & Vorontsov 2003) to the latter. This procedure has been criticised by Roxburgh (2018) on the grounds that in an r_{010} ratio 2N quantities are fitted, while stellar models provide only N frequency phase corresponding to both the r_{01} and r_{10} sequences. This may lead to an overfitting 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 μ_{k}s, the σ_{k}, and the coefficients of Σ, and therefore θ = θ_{⋆}.
The physical model I adopt is a spherically symmetric, nonrotating, nonmagnetic star. The corresponding equations for the stellar structure and evolution are solved numerically via ASTEC (ChristensenDalsgaard 2008a) and those for stellar pulsations via adipls (ChristensenDalsgaard 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 mixinglength formalism of BöhmVitense (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 ^{14}N(p, γ)^{15}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 hydrogenmass fraction, X_{0}, and metallicity, Z_{0}, and the mixinglength parameter, α. This latter is a proportionality coefficient between the meanfree path of a fluid parcel in the convective zone and the pressure scale height. The mixinglength 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.
Lower and upper bounds used for the prior uniform densities for each stellar parameter.
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_{01}, r_{02}, and r_{13} as seismic diagnostics for the first five runs. For the sixth run, r_{01} was replaced with r_{010} to test the potential effects of overfitting. For each simulation I provide three different estimates of the stellar parameters parameters M_{⋆}, t_{⋆}, X_{0}, Z_{0}, α and, when relevant, α_{ov}. I also provide estimates for the initial heliummass fraction, Y_{0}, to facilitate comparison with other studies, this parameter often being reported instead of X_{0}. For a given setup, 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.
Setup of the MCMC simulations and physical quantities inferred using the ASTEC stellar evolution code for 16 Cyg A.
The baseline case was chosen so that it corresponds to the setup used in Metcalfe et al. (2015), who also used ASTEC^{1}. This setup 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 nonseismic 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 r_{13}, r_{02}, and r_{01} observed ratios, respectively. The red lines show the best models for all runs in Tables 3 and 4. 
Fig. 2.
Nonseismic 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 twodimensional joint PDFs for this setup are shown in Fig. 3 alongside the corresponding correlation coefficients. The mass correlates most significantly with X_{0}. The age correlates extremely well with the mixinglength parameter and slightly less with X_{0}. It may be counterintuitive 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 onetoone 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_{0} (Pearson correlation coefficients ≳0.8) both correlate strongly with the luminosity and the radius on the zeroage main sequence (ZAMS). This is related to the massluminosity relation on the ZAMS, in which the meanmolecular 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 meanmolecular weight. This latter corresponds to an increase in the initial hydrogenmass fraction and a decrease of the initial heliummass fraction because the two quantities are perfectly anticorrelated. The interaction between the mass and the hydrogenmass 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_{⋆}, X_{0}, Z_{0}, 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 mixinglength parameter correlates very strongly with the effective temperature. It is well known that these two quantities are related Clayton (1968). Again, the mixinglength parameter allows for setting the initial conditions of the evolutionary sequence. The age also anticorrelates 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 metallicitytohydrogen ratio. It is much stronger (Pearson correlation coefficient >0.9) than the correlation between Z/X and X_{0}. This is because too large variations of X_{0} would also affect the radius and the luminosity, whereas variations of Z_{0} 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 onedimensional 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 setup 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^{2}. 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 neuralnetwork 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 heliummass 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_{0} and t_{⋆}, coupled to a small increase of Y_{0} (all below the 1σ uncertainties from the baseline case). The estimate of Z_{0} 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; ChristensenDalsgaard 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 hydrogenmass fraction and metallicity, at least for stars up until the first dredgeup, 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 (ChristensenDalsgaard 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 twodimensional and onedimensional 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 mixinglength parameter is lower. Within their 68.3% credible intervals, all estimates agree albeit marginally. The estimated X_{0} and Z_{0} 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_{0} and Z_{0} are similar, but deviate more from the baseline case. The estimate of X_{0}, which now decreases by more than 1σ and Z_{0} increases by ∼15% (1.8σ from baseline case) for 16 Cyg A and ∼11% (1.3σ from baseline case). This is a wellknown 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_{0} needs to increase. However, contrary to the nondiffusion case, the estimated age this time decreases and the mixinglength 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 mixinglength parameters. The initial metallicity and hydrogenmass 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 mixinglength 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 longstanding problem in stellar physics. Stellar models use mixinglength 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 semiquantitative analysis, which suggests using penetration to model overshoot below the convective envelope. The depth of the overshoot/penetration zone is set using
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_{⋆}, X_{0}, Z_{0}, α, 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 ChristensenDalsgaard 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 threedimensional. Twodimensional simulations (Hurlburt et al. 1986, 1994; Rogers & Glatzmaier 2005; Rogers et al. 2006) often exhibit penetrative overshoot. On the other hand, in threedimensional 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 plumelike structure in twodimensional than in threedimensional simulations. Our results, even though they do not address chemical mixing, are consistent with the picture of a very small to nonexistent penetration region.
3.3. Using r_{01}/r_{10} ratios
The final MCMC runs for both 16 Cyg A and B aim at assessing the effect of using r_{010} instead of r_{01}. 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_{010} also biases the results towards lower X_{0} and α values and, marginally, towards higher values of t_{⋆}. The uncertainties on the age are also biased when using r_{010} as seismic constraints. The estimates in this case are twice as small as those obtained using r_{01}.
Fig. 5.
Marginal densities for the stellar parameters M, t_{⋆}, X_{0}, Z_{0}, 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 r_{010} is used as a seismic constraint instead of r_{01}. 
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_{010} and not r_{01} as done in this case. This can be explained by the fact that they neglected the nondiagonal terms in Σ. To test this interpretation, I use MCMC simulations for both 16 Cyg A an B using r_{010} 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 nondiagonal terms of Σ. A good indicator of this behaviour is the value the argument of the exponential in Eq. (3), that is the χ^{2} values, of the best models for the baseline case and the case constrained using r_{010}. 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 nondiagonal 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 overfitting 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_{01} (or r_{10}) rather than r_{010} 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 mixinglength 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 subadiabatic 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 Sunlike stars. The present results aim at helping to prepare their modelling.
The full results of the runs can be found at https://amp.phys.au.dk/browse/simulation/767 for 16 Cyg A and https://amp.phys.au.dk/browse/simulation/768 for 16 Cyg B.
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
 Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Angulo, C., Champagne, A. E., & Trautvetter, H.P. 2005, Nucl. Phys. A, 758, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Antia, H. M., & Basu, S. 2011, J. Phys. Conf. Ser., 271, 012034 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., & Sauval, A. J. 2005, ASP Conf. Ser., 336, 25 [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2009, IAU Symp., 253, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., Pinsonneault, M. H., & Wasserburg, G. J. 1995, Rev. Mod. Phys., 67, 781 [NASA ADS] [CrossRef] [Google Scholar]
 Basu, S. 2016, Liv. Rev. Sol. Phys., 13, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Basu, S., & Antia, H. M. 2008, Phys. Rep., 457, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Bazot, M. 2013, EAS Pub. Ser., 63, 105 [CrossRef] [Google Scholar]
 Bazot, M., Bouchy, F., Kjeldsen, H., et al. 2007, A&A, 470, 295 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bazot, M., Bourguignon, S., & ChristensenDalsgaard, J. 2008, Mem. Soc. Astron. It., 79, 660 [NASA ADS] [Google Scholar]
 Bazot, M., Campante, T. L., Chaplin, W. J., et al. 2012, A&A, 544, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bazot, M., ChristensenDalsgaard, J., Gizon, L., & Benomar, O. 2016, MNRAS, 460, 1254 [NASA ADS] [CrossRef] [Google Scholar]
 Bazot, M., Creevey, O., ChristensenDalsgaard, J., & Meléndez, J. 2018, A&A, 619, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bazot, M., Benomar, O., ChristensenDalsgaard, J., et al. 2019, A&A, 623, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bedding, T. R., Butler, R. P., Kjeldsen, H., et al. 2001, ApJ, 549, L105 [NASA ADS] [CrossRef] [Google Scholar]
 Bedding, T. R., Kjeldsen, H., Arentoft, T., et al. 2007, ApJ, 663, 1315 [NASA ADS] [CrossRef] [Google Scholar]
 Bellinger, E. P., Angelou, G. C., Hekker, S., et al. 2016, ApJ, 830, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Bellinger, E. P., Basu, S., Hekker, S., & Ball, W. H. 2017, ApJ, 851, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, J. O., & Berger, J. O. 1985, Statistical Decision Theory and Bayesian Analysis (New York: SpringerVerlag) [CrossRef] [Google Scholar]
 Birnbaum, A. 1962, J. Am. Stat. Assoc., 57, 269 [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Bull. Am. Astron. Soc., 42, 215 [NASA ADS] [Google Scholar]
 Bouchy, F., & Carrier, F. 2001, A&A, 374, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bouchy, F., Bazot, M., Santos, N. C., Vauclair, S., & Sosnowska, D. 2005, A&A, 440, 609 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brooks, S. P., & Gelman, A. 1998, J. Comput. Graph. Stat., 7, 434 [Google Scholar]
 Brown, T. M., ChristensenDalsgaard, J., WeibelMihalas, B., & Gilliland, R. L. 1994, ApJ, 427, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 [Google Scholar]
 Castro, M., Vauclair, S., & Richard, O. 2007, A&A, 463, 755 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ChristensenDalsgaard, J. 2008a, Ap&SS, 316, 13 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2008b, Ap&SS, 316, 113 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Monteiro, M. J. P. F. G., & Thompson, M. J. 1995, MNRAS, 276, 283 [NASA ADS] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [CrossRef] [Google Scholar]
 Clayton, D. 1968, Principles of Stellar Evolution and Nucleosynthesis (Chicago: University of Chicago Press) [Google Scholar]
 Creevey, O. L., Monteiro, M. J. P. F. G., Metcalfe, T. S., et al. 2007, ApJ, 659, 616 [NASA ADS] [CrossRef] [Google Scholar]
 Creevey, O. L., Metcalfe, T. S., Schultheis, M., et al. 2017, A&A, 601, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cunha, M. S., Aerts, C., ChristensenDalsgaard, J., et al. 2007, A&ARv, 14, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Davies, G. R., Chaplin, W. J., Farr, W. M., et al. 2015, MNRAS, 446, 2959 [NASA ADS] [CrossRef] [Google Scholar]
 de Porto Mello, G. F., da Silva, R., da Silva, L., & de Nader, R. V. 2014, A&A, 563, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
 Gough, D. O. 2012, ASP Conf. Ser., 462, 429 [NASA ADS] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse, 15 [Google Scholar]
 Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Guzik, J. A., Watson, L. S., & Cox, A. N. 2006, Mem. Soc. Astron. It., 77, 389 [NASA ADS] [Google Scholar]
 Haario, H., Saksman, E., & Tamminen, J. 2001, Bernoulli, 7, 223 [CrossRef] [MathSciNet] [Google Scholar]
 Herwig, F. 2000, A&A, 360, 952 [NASA ADS] [Google Scholar]
 Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J.P. 1994, ApJ, 421, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Kjeldsen, H., Bedding, T. R., Butler, R. P., et al. 2005, ApJ, 635, 1281 [NASA ADS] [CrossRef] [Google Scholar]
 Kjeldsen, H., Bedding, T. R., & ChristensenDalsgaard, J. 2008, ApJ, 683, L175 [NASA ADS] [CrossRef] [Google Scholar]
 Korre, L., Garaud, P., & Brummell, N. H. 2019, MNRAS, 484, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalfe, T. S., Creevey, O. L., & ChristensenDalsgaard, J. 2009, ApJ, 699, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalfe, T. S., Chaplin, W. J., Appourchaux, T., et al. 2012, ApJ, 748, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalfe, T. S., Creevey, O. L., Doğan, G., et al. 2014, ApJS, 214, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalfe, T. S., Creevey, O. L., & Davies, G. R. 2015, ApJ, 811, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Michaud, G., & Proffitt, C. R. 1993, ASP Conf. Ser., 40, 246 [NASA ADS] [Google Scholar]
 Monteiro, M. J. P. F. G., ChristensenDalsgaard, J., & Thompson, M. J. 1994, A&A, 283, 247 [NASA ADS] [Google Scholar]
 Moore, D. W. 1967, IAU Symp., 28, 405 [NASA ADS] [Google Scholar]
 Richard, O., Vauclair, S., Charbonnel, C., & Dziembowski, W. A. 1996, A&A, 312, 1000 [NASA ADS] [Google Scholar]
 Richard, O., Michaud, G., & Richer, J. 2001, ApJ, 558, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, O., Michaud, G., Richer, J., et al. 2002, ApJ, 568, 979 [Google Scholar]
 Robert, C. 2007, The Bayesian Choice: From DecisionTheoretic Foundations to Computational Implementation, Springer Texts in Statistics (New York: Springer) [Google Scholar]
 Robert, C. P., & Casella, G. 2005, Monte Carlo Statistical Methods (Springer Texts in Statistics) (Secaucus, NJ, USA: SpringerVerlag, New York, Inc.) [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2005, ApJ, 620, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., Glatzmaier, G. A., & Jones, C. A. 2006, ApJ, 653, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W. 2017, A&A, 604, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roxburgh, I. W. 2018, ArXiv eprints [arXiv:1808.07556] [Google Scholar]
 Roxburgh, I. W., & Vorontsov, S. V. 2003, A&A, 411, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Silva Aguirre, V., Basu, S., Brandão, I. M., et al. 2013, ApJ, 769, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Silva Aguirre, V., Lund, M. N., Antia, H. M., et al. 2017, ApJ, 835, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Singh, H. P., Roxburgh, I. W., & Chan, K. L. 1995, A&A, 295, 703 [NASA ADS] [Google Scholar]
 Vauclair, S., Laymand, M., Bouchy, F., et al. 2008, A&A, 482, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Veronis, G. 1963, ApJ, 137, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Meakin, C., Prat, V., & Arnett, D. 2015, A&A, 580, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 von Toussaint, U. 2011, Rev. Mod. Phys., 83, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
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 pseudoscale (Gelman & Rubin 1992) and multivariate pseudoscale (Brooks & Gelman 1998) criteria. The general ruleofthumb 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 pseudoscale factor is always superior to pseudoscale 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 setup 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.
Pseudoscale (R, full lines) and multivariate pseudoscale factors (R_{p}, 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.
Pseudoscale (R, full lines) and multivariate pseudoscale factors (R_{p}, 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
Lower and upper bounds used for the prior uniform densities for each stellar parameter.
Setup of the MCMC simulations and physical quantities inferred using the ASTEC stellar evolution code for 16 Cyg A.
All Figures
Fig. 1.
Seismic diagnostic for 16 Cyg A (left panel) and B (right panel). The squares, triangles, and circles denote the observed r_{13}, r_{02}, and r_{01} observed ratios, respectively. The red lines show the best models for all runs in Tables 3 and 4. 

In the text 
Fig. 2.
Nonseismic 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_{⋆}, X_{0}, Z_{0}, 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_{⋆}, X_{0}, Z_{0}, α, 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_{⋆}, X_{0}, Z_{0}, 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 r_{010} is used as a seismic constraint instead of r_{01}. 

In the text 
Fig. A.1.
Pseudoscale (R, full lines) and multivariate pseudoscale factors (R_{p}, 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.
Pseudoscale (R, full lines) and multivariate pseudoscale factors (R_{p}, 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.