Issue 
A&A
Volume 644, December 2020



Article Number  A37  
Number of page(s)  24  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202038522  
Published online  30 November 2020 
Thorough characterisation of the 16 Cygni system
I. Forward seismic modelling with WhoSGlAd
^{1}
Institut d’Astrophysique et Géophysique de l’Université de Liège, Allée du 6 août 17, 4000 Liège, Belgium
email: martin.farnir@uliege.be
^{2}
Observatoire de Genève, Université de Genève, 51 Ch. Des Maillettes, 1290 Sauverny, Switzerland
^{3}
LESIA, Observatoire de Paris, PSL Research University, CNRS, Université Pierre et Marie Curie, Université Paris Diderot, 92195 Meudon, France
Received:
28
May
2020
Accepted:
7
October
2020
Context. Being part of the brightest solarlike stars, and close solar analogues, the 16 Cygni system is of great interest to the scientific community and may provide insight into the past and future evolution of our Sun. It has been observed thoroughly by the Kepler satellite, which provided us with data of an unprecedented quality.
Aims. This paper is the first of a series aiming to extensively characterise the system. We test several choices of micro and macrophysics to highlight their effects on optimal stellar parameters and provide realistic stellar parameter ranges.
Methods. We used a recently developed method, WhoSGlAd, that takes the utmost advantage of the whole oscillation spectrum of solarlike stars by simultaneously adjusting the acoustic glitches and the smoothly varying trend. For each choice of input physics, we computed models which account, at best, for a set of seismic indicators that are representative of the stellar structure and are as uncorrelated as possible. The search for optimal models was carried out through a LevenbergMarquardt minimisation. First, we found individual optimal models for both stars. We then selected the best candidates to fit both stars while imposing a common age and composition.
Results. We computed realistic ranges of stellar parameters for individual stars. We also provide two models of the system regarded as a whole. We were not able to build binary models with the whole set of choices of input physics considered for individual stars as our constraints seem too stringent. We may need to include additional parameters to the optimal model search or invoke nonstandard physical processes.
Key words: asteroseismology / stars: oscillations / stars: solartype / stars: abundances
© ESO 2020
1. Introduction
In the past decade, the CoRoT (Baglin et al. 2009) and Kepler (Borucki et al. 2010) space missions have provided the stellar physics community with a wealth of data of unprecedented quality for solarlike stars. Such data allow stellar scientists, through the use of asteroseismology, to put their models to the test and to provide stringent constraints on the physical processes at hand, therefore highlighting the current shortcomings in the modelling (e.g. mixing processes, angular momentum transport, starplanet interaction). In addition, studying solarlike stars enables us to gather invaluable insight into the past and future of our Sun.
The 16 Cygni system is of great interest as it consists of binary solar twins which have been observed continuously for 928 days. Both stars are therefore among the solarlike pulsators with the best data available for seismic studies. Moreover, a great amount of information has yet to be accounted for. For example, differences in superficial lithium abundances remain unexplained (Friel et al. 1993 and King et al. 1997) observed that the B component is at least four times more Li depleted than its twin). The presence of a jovian companion to 16 Cygni B (Cochran et al. 1997) has been argued by Deal et al. (2015) to be the possible cause. This specific example illustrates that the system is an ideal testbench to constrain stellar models as well as to test nonstandard physical processes while taking advantage of asteroseismic techniques.
Solarlike oscillations, as both stars display, are stochastically excited by the outer convective layer. Such oscillation spectra may present the following two main features: a regular pattern, referred to as the smooth part of the spectrum, and an oscillating pattern of low amplitude, the glitch. An acoustic glitch is the oscillating signal observed in frequencies, which is caused by a sharp variation – compared to the wavelength of the oscillating mode – variation in the stellar structure. The first mention of the possible use of such signatures was by Vorontsov (1988) and Gough (1990) who theoretically demonstrated the effect of a sharp feature in the stellar structure on oscillation frequencies, either directly or in the second frequency differences. For example, in solarlike stars, we have the helium glitch, caused by a depression in the first adiabatic index^{1} in the second helium ionisation zone, and the convection zone glitch, due to the variation of the temperature gradient at the base of the convective envelope zone. They may constrain the surface helium content, inaccessible by other means for solarlike stars, (e.g. Basu et al. 2004; Verma et al. 2014) and the total extent of the envelope convective zone, as well as the mixing processes at hand that might explain such an extent (e.g. Monteiro et al. 2000).
In a previous paper, Farnir et al. (2019) described a new method to provide a robust analysis of the solarlike oscillation frequencies simultaneously accounting for the smooth and glitch parts, the WhoSGlAd (Whole Spectrum and Glitches Adjustment) method. This method relies on the GramSchmidt orthonormalisation process to define seismic indicators as uncorrelated as possible. It shows a great potential to provide precise, accurate and statistically relevant constraints on stellar physics. Compared to other seismic methods accounting for the glitches signature, the WhoSGlAd method has the advantage of decorrelating the information contained in both components of the oscillation spectrum (smooth and glitch parts) while accounting for them simultaneously. This leads to constraints which are, in turn, the least correlated possible and more stringent. Moreover, the measured frequencies are not fitted individually as it introduces large correlations. Rather, we use seismic indicators defined to be representative of the stellar structure and as little correlated as possible. This enables us to compute optimal models as accurate as possible.
This paper is part of a series of publications dedicated at providing the most accurate and complete picture of the 16Cyg binary system. In this first study, our goal is to establish a large sample of reliable structural models analysing the degeneracies stemming from variations in the micro and macrophysical prescriptions using our new consistent seismic modelling technique, WhoSGlAd. We model the system using asteroseismic, spectroscopic and interferometric constraints considering both 16CygA (KIC12069424) and 16CygB (KIC12069449) independently and as a joint system. We provide a suitable set of models for structural inversions to be studied in a second paper. Our thorough analysis also paves the way for an indepth description of potential traces of nonstandard processes acting (or having acted) during the history of the system. These include the effects of angular momentum transport processes (Eggenberger et al. 2010, 2019) as well as the effects of planetary formation and accretion on the lithium abundances of both stars (Deal et al. 2015; Thévenin et al. 2017).
The paper is structured as follows. First, we present in Sect. 2 the general methodology and recall the basics of the WhoSGlAd method. We then model the system. This is done in two steps. To take advantage of the great precision of the data for each star, we first provide, in Sect. 3, separate adjustments while testing different choices of micro and macrophysics. This allows to provide robust stellar parameter ranges accounting for the modelling uncertainties as well as to show discrepancies in the modelling for some cases. We select in Sect. 4 the models having consistent ages and initial compositions as initial guesses to compute models imposing a common age and initial composition, as those stars should have formed from a single molecular cloud. Even though no specific interaction between both stars is taken into account during their evolution, we refer to those models as binary models. We discuss the results in Sect. 5. Finally, we conclude our paper in Sect. 6.
2. Methodology
In the current section, we describe the optimisation scheme and the seismic and nonseismic constraints used. We then present the basic principle of the WhoSGlAd method. Finally, we describe the physics included in the models as well as the considered variations.
2.1. Bestfit model search
The search for bestfit models is undertaken by a LevenbergMarquardt (LM) algorithm. In doing so, we compare observed values of a set of constraints with model values, computed on the fly, through a χ^{2} function, to be minimised, defined as:
where C represents the N constraints, the obs (resp. mod) subscript the observed (resp. model) values and σ their associated standard deviations.
Except when mentioned otherwise, the set of constraints consists of the Δ, , , and A_{He} seismic indicators (Farnir et al. 2019) presented in Sect. 2.2 and the free parameters adopted in the modelling procedure are the mass (M), age (t), initial hydrogen abundance (X_{0}) and, metallicity (Z/X)_{0} of the considered star. Other nonseismic data, such as the effective temperature (T_{eff}), interferometric radius (R), or the metallicity ([Fe/H]), are used to discriminate between the several choices of input physics. In some cases, and when so stated, nonseismic data may be used as constraints to the model search while relaxing the mixing length parameter or including turbulent diffusive mixing with a free coefficient (see Sect. 2.3 for a description of the physics included in the models). Those data are gathered in Table 1. Finally, we do not use as a constraint the luminosity (L) from Metcalfe et al. (2012) as it results from asteroseismic modelling and would not be independent of our study. Instead, we compute it from the observed interferometric radius (White et al. 2013) and the definition of the effective temperature: , where s is the StefanBoltzmann constant.
Set of nonseismic data used througout this paper.
2.2. WhoSGlAd principle and seismic indicators
We recall here the set of WhoSGlAd seismic indicators used in the fitting procedure as well as the basics of the method. For a more detailed description, see Farnir et al. (2019).
Principle. The WhoSGlAd method relies on GramSchmidt’s orthogonalisation. To represent the observed frequencies, we define a Euclidean vector space of functions of the spherical degree, l, and radial order, n (only m = 0 modes are considered). The N observed frequencies at a given value of l are regarded as unknown vector functions of n and l which we write ν_{l} = (ν_{l, n1}, …, ν_{l, nN}). Two notable functions are the identity, 1, and linear function of the radial order, n_{l} = (n_{l, 1}, …, n_{l, N}).
Given two vector quantities, say the observed and theoretical vectors of frequencies, ν_{obs} and ν_{t}, we may define their scalar product as:
with σ_{i} the uncertainties associated with each component. From this scalar product is defined the norm of a vector ν_{obs}:
We may also define the weighted mean of a quantity, according to our scalar product and normalisation, as:
The functions used to describe the frequencies are separated into two contributions: a smooth part, represented by secondorder polynomials of n, and a glitch part, represented by oscillating functions of the frequency. The form of those functions is given in Appendix A. Using GramSchmidt’s orthogonalisation process, we build a basis of functions over that vector space. To provide an adjustment of the observed frequencies, those are projected over the basis elements. This provides completely independent coefficients which are combined into seismic indicators as little correlated as possible. One of the main advantages of this approach is that the glitch part of the adjustment is completely independent of the smooth part, even though both adjustments are carried out simultaneously. We define the seismic indicators as follows:
The large separation for modes of degree l. Corresponds to the value of the slope of the frequencies decorrelated from the contribution of the acoustic glitches and expressed as a linear function of the radial order n for each spherical degree l:
This is equivalent to its standard definition obtained through a linear regression (see e.g. Reese et al. 2012).
The large separation. Corresponds to the weighted mean value, , of the individual large separations for each spherical degree l
with σ(Δ_{l}) the uncertainty on the large separation of degree l.
The normalised small separations between degrees 0 and l.
Those indicators are analogous to the mean value of the local small separation ratios defined by Roxburgh & Vorontsov (2003) but are again completely independent from the contribution of acoustic glitches.
The helium glitch amplitude.
where δν_{He} is the helium glitch component.
Beside those indicators, we may define complementary seismic indicators which are presented in Appendix B. Those were not part of the constraints.
Table 2 gathers the values of the considered seismic indicators computed using the modes defined over the full length of the Kepler mission determined by Davies et al. (2015). We take out from those modes the ones with uncertainties above 1.5 μHz. This mostly corresponds to high frequency modes. A brief discussion of this choice is given in Appendix C. We have corrected the frequencies for the surface effects by using the power law prescribed by Kjeldsen et al. (2008) and the a and b coefficients fitted by Sonoi et al. (2015) as a function of T_{eff} and g. The authors have undertaken this coefficient adjustment by comparing the adiabatic frequencies of patched models based on 3D simulations and that of unpatched standard 1D models.
Observed seismic indicators.
2.3. Models
Unless specified otherwise all the models are computed using the CLES stellar evolution code (Scuflaire et al. 2008b) with the AGSS09 solar chemical mixture (Asplund et al. 2009), the OPAL opacity table (Iglesias & Rogers 1996) combined with that of Ferguson et al. (2005) at low temperatures, the FreeEOS software to generate the equation of state table (Cassisi et al. 2003), and the nuclear reactions rates prescribed by Adelberger et al. (2011). We also use the mixing length theory (Cox & Giuli 1968), with the solar calibrated value of α_{MLT} = l/H_{p} = 1.82 (where l is the mixing length and H_{p} the pressure scale height), to parametrise the mixing inside convective regions. This value is the result of a solar calibration that we carried out using the same set of input physics as described in the present section. The microscopic diffusion of elements is included and treated as in Thoul et al. (1994). The models do not include rotation and, therefore, rotationinduced mixing. Unless specified otherwise, models do not include overshooting at the boundary of convective layers. The temperature conditions above the photosphere are determined using an Eddington T(τ) relation, τ being the optical depth. We choose such a relation to remain consistent with Sonoi et al. (2015) whose fitted coefficients are used to correct surface effects on the observed frequencies. From now on, to distinguish from the several physical variations, we refer to the models with this specific set of input physics as the reference models. Finally, we compute theoretical adiabatic oscillation frequencies for each model via the LOSC oscillation code (Scuflaire et al. 2008a).
2.4. Variations in the input physics
As mentioned earlier, to provide the most reliable set of stellar parameter ranges while accounting for the choice of micro and macrophysics, we test choices by changing one ingredient at the time from the reference models. Those variations are:

The GN93 solar reference mixture (Grevesse & Noels 1993), in light blue in the figures (see Sect. 3.1.2);

The opacities from the opacity project (Badnell et al. 2005), denoted OP in light brown, the Los Alamos opacities (Colgan et al. 2016), written OPLIB in beige (see Sect. 3.1.2);

The CEFF equation of state (ChristensenDalsgaard & Daeppen 1992), in dark brown, and the revised OPAL equation of state (Rogers & Nayfonov 2002), written OPAL05 in grey (see Sect. 3.1.2);

A different choice of mixing length coefficient (α_{MLT} = 1.7), in yellow (see Sects. 3.1.2 and 3.2);

The inclusion (or not) of microscopic diffusion, in light pink (see Sect. 3.1.3);

The inclusion of turbulent mixing of chemical elements following the relation for the turbulent mixing coefficient (in cm^{2} s^{−1}), where ρ is the density, ρ_{0} the density at the bottom of the convective envelope and D_{turb}, n and D_{ct} are fixed at 7500, −3 and 0 respectively (Proffitt & Michaud 1991), shown in purple (see Sects. 3.1.3 and 3.2);

The inclusion of overshooting extending outside convective regions over a distance d = α_{ov}min(H_{p},h) where α_{ov} is the overshooting parameter, H_{p} the local pressure scale height and h the thickness of the convective region. The temperature gradient in the overshooting region is set to the radiative one and the mixing is instantaneous. We either include overshooting above the convective core, denoted α_{ov} and shown in red, or below the convective envelope, written α_{un} in khaki and referred to as “undershoot”. Both values are set to 0.1 (see Sect. 3.1.4);

The effect of a different choice of temperature profile above the stellar photosphere, in orange. We use the model temperature profile of the quiet sun by Vernazza et al. (1981) for which an analytical formulation may be found in Paxton et al. (2013) (see Sect. 3.1.5);

The impact of the surface effects, computing a model fitting seismic indicators defined with stellar frequencies which are not corrected for surface effects in dark green. Their values are shown in Table D.1. See also Sect. 5.2 and Appendix D.
3. 16 Cygni A and B seen as separate stars
In the present section, we look for individual models of each star representative of the observed data and accounting at best for the modelling uncertainties. The stellar parameters for every best model estimates are displayed in Appendix E. To find individual models, we test several choices of input physics without any specific hypothesis about the binarity of the stars. This allows to take advantage of the unprecedented quality of the data. The first part of this study is subdivided in two steps. We start by only considering seismic constraints. Then, we add nonseismic constraints, in Sect. 3.2, to further improve the models. The advantage of first computing individual models for each star is that it allows to have the same amount of constraints as free parameters and to obtain an exact solution, from a statistical point of view, to the minimisation process.
3.1. Fitting seismic constraints only
In the present section, we present the results of the modelling considering only seismic indicators. This allows to show the impact of the seismic indicators alone on optimal results as well as the possible limitations of such an exclusive approach. Furthermore, we test several choices of micro and macrophysics. This enables us to highlight their influence on the set of optimal parameters we retrieve. The models are computed as described in Sect. 2.3 while changing only one physical ingredient at a time.
The individual models for both components are shown in Figs. 1 and 2 (see Table E.3 for individual parameter values.) We show, in the upper panel of the figures the age versus the mass of the optimal model for each variation in the input physics along with the associated uncertainties. The middle panel displays the position of those models in a HertzsprungRussell (HR) diagram. We also represent, as a black box, the observed effective temperature and luminosity computed from the interferometric radius (White et al. 2013). Finally, the lower panels represent the initial hydrogen versus initial metallicity and their uncertainties for each model.
Fig. 1. Summary of 16 Cyg A bestfit models represented in a Mass – Age diagram (top panel), HR diagram (middle panel) and initial hydrogen abundance versus metal composition diagram (bottom panel). The luminosity and effective temperature constraints from White et al. (2013) are represented in the HR diagram as a black box. 
In both figures, the reference model is represented in dark blue and denoted AGSS09. For 16 Cyg A, it actually corresponds to the model presented in Sect. 5 of Farnir et al. (2019). Their Fig. 14 illustrates that the use of the WhoSGlAd seismic indicators as constraints allows to provide a proper agreement between observed and model frequencies.
3.1.1. Influence of seismic constraints on stellar parameters
Farnir et al. (2019) showed that the WhoSGlAd helium glitch amplitude is a good proxy of the surface helium abundance. This means that, when requiring our models to reproduce the observed helium glitch amplitude, we require a specific helium abundance. To illustrate this statement, we plot in Figs. 3 and 4 the surface helium content as a function of the surface metallicity of each model. We also represent, as a black box, the spectroscopic metallicity from Ramírez et al. (2009) and the asteroseismic Y_{s} range computed by Verma et al. (2014), taking advantage of the information contained in the helium glitch, along with their associated uncertainties. We indeed observe that the surface helium abundance is well constrained and, in most cases, in agreement with the study of Verma et al. (2014). They computed ranges of [0.231,0.251] and [0.218,0.266] for the A and B components respectively, which encapsulate most of our values. However, we also note a small scatter in the values. Again, Farnir et al. (2019) showed that the helium glitch amplitude is both correlated to the surface helium and metal abundances, with the helium abundance being the dominant factor. They indeed observe that, at constant surface helium abundance, a lower surface metal abundance or, at constant surface metal abundance, a higher surface helium abundance may both lead to a greater glitch amplitude. This is in fact due to a shift in the position of the adiabat in the second helium ionisation zone, where the first adiabatic index, , presents a large depression. This allows us to account for the small scatter observed. A direct consequence of the helium glitch amplitude constraint is the anticorrelation between the initial metallicity and helium abundance that we observe in the lower panels of Figs. 1 and 2.
Fig. 3. Surface helium abundance versus metallicity for 16 Cygni A. The black box represents the spectroscopic metallicity value computed by Ramírez et al. (2009) and the surface helium value from Verma et al. (2014) along with the corresponding uncertainties. 
Fig. 4. Surface helium abundance versus metallicity for 16 Cygni B. The black box represents the spectroscopic metallicity value computed by Ramírez et al. (2009) and the surface helium value from Verma et al. (2014) along with the corresponding uncertainties. 
Moreover, we point out that most models do not account for the spectroscopic metallicity constraint. This is clearly visible in Figs. 3 and 4. This does not come as a surprise as the presented models do not yet include the metallicity constraint in the fitting procedure. This also shows that the information contained in the helium glitch amplitude and the surface metallicity are complementary and it comes as a necessity to take advantage of both to provide the most accurate model possible.
Now looking at the middle panels of Figs. 1 and 2, the first striking feature is the fact that most models lie on a straight line. Such line corresponds to the locus of models of constant radius. This almost constance of the radii stems from the Δ indicator which provides a constraint on the mean stellar density (Farnir et al. 2019; Ulrich 1986). Thus, the models of constant mean density have almost constant stellar radii, as long as the mass remains close to constant. Actually, this is what we observe as the mean radius values of our models are 1.22 R_{⊙} and 1.11 R_{⊙} for 16CygA and B respectively (Typical uncertainties are of 0.02 R_{⊙} and 0.01 R_{⊙} respectively. Individual values are shown in Table E.3). This is in good agreement with the values in Table 1. We note that some models do not lie on the straight line with the other models. Such models are those with masses values that differ significantly from the mean value of other models. Finally, we observe that many models do not fall in the effective temperature and luminosity observational boxes. This comes from the fact that those constraints are not yet part of the modelling procedure and shows that their inclusion is necessary to provide the most accurate picture of the system.
3.1.2. Effect of the metallicity reference, opacity table, and equation of state
As is clearly visible in Figs. 1 and 2 using either the GN93 solar reference mixture or a greater value of the mixing length parameter produces models for both stellar components which are more luminous and have a greater effective temperature than the reference models. Furthermore, looking at Figs. 3 and 4, we observe that the GN93 solar reference tends to produce more metallic models, directly stemming from the fact that this solar reference is indeed more metallic than the AGSS09 one. However, in terms of metallicity, the mixing length parameter has opposite impacts on the two stellar components. A decrease of its value leads, in the case of 16 Cyg A, to models which become less metallic, while it produces more metallic models for 16 Cyg B.
Now looking at the influence of the opacity tables, we note that the OPLIB table leads to a model for 16 Cyg A which has a greater effective temperature while the effect is barely visible for 16 Cyg B. Conversely, the OP opacity table leads to models which are cooler for both stars but the effect is not as pronounced. The effect of the opacity tables is not clear on the surface composition as both models react in different ways.
Finally, the use of a different equation of state table also produces differential effects on both stars. On the first hand, in the case of 16 Cygni A, using either the CEFF or OPAL05 tables lead to hotter, more luminous stars and with decreased surface helium and metal abundances. On the other hand, for 16 Cygni B, both tables have very little influence on the position of the star in the HR diagram (see Fig. 2). However, the use of the OPAL05 table leads to a model of the B component which is both richer in helium and metals at its surface. The impact of the CEFF equation of state is barely visible.
3.1.3. Impact of diffusion
We note that the models we compute without diffusion of chemical elements are both older, heavier and richer in hydrogen than the reference, as represented in light pink in Figs. 1 and 2. As showed by Farnir et al. (2019), at a specific composition, more massive models present a stronger helium glitch signature. Therefore, to reproduce the observed signature of both stars, the models need to be poorer in helium as more massive models are favoured. This is even reinforced by the fact that no diffusion is included and the initial helium abundance has to match that of the surface. Moreover, the difference in surface helium abundance between models with and without diffusion is systematically of about 0.02 as was noted by Verma et al. (2019). This confirms their observation of the importance of diffusion in low mass stars of solar metallicity.
We observe that nonseismic data, that will be considered in Sect. 3.2, are not accounted for in most cases. We note that one way to account for them is to reduce the impact of microscopic diffusion, either partially for 16 Cygni A – by involving additional mixing processes as turbulent mixing – or completely for 16 Cygni B. This contradicts the conclusions of Buldgen et al. (2016a) who determined that models with increased diffusion efficiency (with diffusion velocities higher of about 10%) could help reproduce inversion results for the A star. However, this agrees with their second study (Buldgen et al. 2016b) where they noted that reducing the efficiency of diffusion allowed the computation of models consistent with the inversion results. This, however, lead their study to inconsistencies in surface composition between the two stars.
We also show the impact of the inclusion of turbulent mixing on the modelling by computing models with a turbulent mixing coefficient fixed at a value of D_{turb} = 7500 cm^{2} s^{−1}. We note that those models tend, for both stars, to be more luminous, hotter and less metallic. The overall agreement with nonseismic data is thus improved. We show in Sect. 3.2.1 the influence of the value of the turbulent mixing coefficient on the optimal results.
3.1.4. Extension of convective layers
To analyse the effect of the extension of the convective core – during premain sequence – on the stellar evolution, we include instantaneous overshooting in some of our models. Those are displayed in red in the figures. We note that the effect on the optimal models is not obvious. Indeed, such models are almost indistinguishable from the reference ones and lie within one another uncertainties. The same goes for the surface compositions retrieved. Including a greater value of the overshooting parameter (α_{ov} = 0.2 instead of 0.1) leads to great differences in the behaviours of both stars. Indeed, while the model for 16 Cyg B remains rather similar to that with a lower value of overshooting, thus similar to the one without overshooting, the model for the A component becomes significantly less massive, older, metal poor and with a smaller radius (see Table E.3). The reason for such a difference is that only the A component, because of its slightly greater mass and smaller metallicity, is able to maintain a convective core during the main sequence with such an overshooting parameter value. Therefore, its structure becomes significantly different from that of a model without overshooting which has a radiative core during the main sequence. However, this model maintaining a convective core on the main sequence seems to be a curiosity of the minimisation method. The retroaction of the presence of the convective core ultimately leads to a significant decrease of the optimal mass as well as the metallicity, becoming significantly subsolar ([Fe/H] = −0.39 ± 0.01). This further indicates the necessity to use nonseismic constraints to obtain proper models. We can therefore safely discard this model. Additionally, we may note that the significantly lower mass of the A model leads to a significantly smaller radius because of the Δ constraint, which again contradicts observations. Moreover, this pair of model has now a less massive A component than the B. Although, by a simple argument of scaling relations, we obtain a hint that the A component should be heavier than the B. Indeed, from Kjeldsen & Bedding (1995) we have that . With ν_{max, A} = 2188 μHz, ν_{max, B} = 2561 μHz taken from Lund et al. (2017), and the values of the effective temperatures and large separations presented earlier, we expect the mass of the A component to be larger than the other (M_{A} ≃ 1.04M_{B}). This is what we observe for most of our models.
The inclusion of undershooting below the base of the convective zone has no significant impact on the optimal stellar parameters beside a reduction of the individual uncertainties. Initial compositions and ages for both stars now fall out of each others uncertainties and the models are dismissed as valid candidates to compute binary models in Sect. 4. This statement is clearly illustrated in Fig. 7 where the initial hydrogen, metal abundances and ages of both stars are plotted against one another. The straight line displays the locus of identical parameters for both stars. In the upper panel, the khaki cross, representing the model with undershooting, does not meet the line any more.
3.1.5. Effects of the atmosphere and mixing length coefficient
We show the influence of the atmosphere on the optimal stellar parameters by using a temperature profile above the photosphere as in Vernazza et al. (1981) with a specifically calibrated value of α_{MLT} = 2.02 (dark pink small circle in the figures). We observe that the optimal model is very similar to the reference model. Indeed, both pairs of models lie within the uncertainties of one another. The optimal parameters are given in Table E.4. However, the models become hotter and more luminous. They are thus closer to the observed luminosities and effective temperatures. This indicates that it allows to provide better models in terms of spectroscopic constraints while preserving rather similar stellar parameters compared to the case using the Eddington relation. We must also point out that the use of a different temperature profile leads to significant changes from the reference models but the calibration of the mixing length parameter compensates for it. To illustrate this, we compute models with the temperature profile of Vernazza et al. (1981) while using the reference value of α_{MLT} = 1.82 calibrated for an Eddington T − τ relation. We observe that our models, shown in orange in the figures, are very similar to models using a lower value of the mixing length parameter, displayed in yellow. This is especially true for 16 Cygni B for which both models lie within respective uncertainties. For both stars, the computed models are older and lighter.
3.2. Fitting nonseismic constraints
To further improve individual models for each star, we may include nonseismic constraints into the minimisation process. We indeed note that, in most cases, they are not accounted for. Therefore, we use those constraints and include additional free parameters to add degrees of freedom. The considered constraints are the effective temperature (White et al. 2013) and the spectroscopic metallicity (Ramírez et al. 2009). The additional free parameters are the mixing length parameter α_{MLT} and the turbulent mixing coefficient D_{turb}. In what follows, we specify in every case which of those constraints and free parameters are used.
3.2.1. Accounting for the effective temperature
From Figs. 1 and 2, we expect that it is possible to improve the agreement with the observed effective temperature by increasing the value of α_{MLT}. As a consequence of the Δ constraint, which we showed in Sect. 3.1.1 to be a proper constraint on the radius, we also expect to produce better model luminosities.
To demonstrate that a variation in the mixing length parameter is indeed responsible for the improvement of the model effective temperature, we compute several models with different mixing length parameters. We only include seismic constraints in the fitting procedure to show the influence of α_{MLT} alone on optimal parameters. This is shown in Figs. 5 and 6 for 16 Cyg A and B respectively and the stellar parameters are given in Table E.5. The values considered are 1.7 [light blue], the solar value of 1.82 [dark blue], and 2.0 [brown]. Those models are connected with a blue line to improve visibility. We also display models with several choices for the turbulent mixing coefficient: D_{turb} = 2000 cm^{2} s^{−1} [light pink], 5000 cm^{2}s^{−1} [brown], 7500 cm^{2} s^{−1} [grey], 10 000 cm^{2} s^{−1} [dark pink], and no turbulent mixing [dark blue]. These are connected in red and stellar parameters are gathered in Table E.6. We observe that an increase of α_{MLT} leads to a better agreement with the observed effective temperature for both stars. Moreover, we note that the inclusion of turbulent mixing improves the agreement in effective temperature and metallicity for both stars. Even so, both stars exhibit too large values of the metallicity compared to the observed ones. We also note that, for 16 Cyg B, turbulent mixing alone is not sufficient for the observed and model T_{eff} to match. In addition, we observe a clear effect of saturation of the turbulent mixing coefficient, which occurs above a threshold value that is already exceeded by the considered values. The results are almost indistinguishable no matter which value is chosen. Therefore, the turbulent mixing and mixing length coefficient both have a impact on the model effective temperature (and thus luminosity) but using the turbulent mixing coefficient as a free parameter would be meaningless. Furthermore, we also display in both figures models which did not include microscopic diffusion of the chemical elements. We observe that those two models differ greatly from the models including turbulent diffusion, both being highly hotter and less metallic. The model for 16 Cygni B even properly reproduces the observed metallicity and effective temperature, although they are not yet part of the constraints.
Fig. 5. Variation of the optimal metallicity and effective temperature for 16 Cyg A with the mixing length parameter and turbulent mixing. Models of different α_{MLT} values are connected in blue while models with various turbulent mixing are connected by a red line. 
Fig. 6. Variation of the optimal metallicity and effective temperature for 16 Cyg B with the mixing length parameter and turbulent mixing. Models of different α_{MLT} values are connected in blue while models with various turbulent mixing are connected by a red line. 
Fig. 7. Comparison of adjusted stellar parameters for both stars at a given physics. The straight line shows the locus of identical stellar parameters for both stars. 
Finally, we compute two models for each star with a free mixing length coefficient and including T_{eff} as a constraint, with and without turbulent mixing (D_{turb} = 7500 cm^{2} s^{−1}). The set of input physics is that of the reference model. Both are displayed in Figs. 5 and 6 as orange diamonds labelled “T_{eff} − D_{turb}” and as yellow pentagons labelled “T_{eff}”. For the A component, we see that including turbulent mixing improves the results. However the agreement with the nonseismic constraints is not improved compared to the model including only turbulent mixing and the opacity table OPLIB with only seismic constraints, that already accounted for the effective temperature, as displayed on Fig. 1. Conversely, for the B component, the improvement is significant and the effective temperature is now well adjusted. Nevertheless, the inclusion of turbulent mixing does not have a significant impact on the resulting agreement with the nonseismic data. Furthermore, in both cases, the metallicity is still not properly accounted for. The corresponding set of stellar parameters is presented in Tables 3 and 4.
Adjusted stellar parameters including the T_{eff} with the reference set of input physics.
Adjusted stellar parameters including the T_{eff} constraint and turbulent mixing with a coefficient of D_{turb} = 7500 cm^{2} s^{−1}.
Analysing our results, we first notice from Table 3 that the calibrated mixing length parameters are very different. Indeed, looking at Fig. 2 from Magic et al. (2015), who used the same solar reference mixture as we do, we would expect that both adjusted values would remain rather close to the solar calibrated value (1.82 in our case) while being smaller as both stars have higher effective temperatures and smaller surface gravities than the Sun. Thus, both differences should translate into a lower mixing length parameter. However, we observe that for 16 Cygni A, and within the error bars, α_{MLT} remains solar while the calibrated value for 16 Cygni B is significantly higher than the solar value. We may need to invoke a special physical process acting on any of the components while being inefficient for the second to produce such a differential effect. Some of the possible scenarios are discussed in Sect. 5.3.
We note that slightly more massive and less metallic models than previously are favoured. Such an effect stems from the fact that more massive models are hotter and thus in better agreement with the observed effective temperature, while keeping the same density, because of the Δ constraint.
Another way to better reproduce the observed position in the HR diagram is to include extra mixing counteracting the diffusion of chemical elements. Indeed, for 16 Cyg B, the model computed without diffusion already accounts for these constraints. Which is striking as those were not yet constraints of the fit. What is more striking is that it also reproduces the spectroscopic metallicity. This strongly suggests that additional mixing processes may be necessary to properly and accurately model this star. For its twin, we note that the inclusion of turbulent mixing, a different opacity table, the Los Alamos one, or the use of a different equation of state, either CEFF or OPAL05, could help us account for the observed luminosity and effective temperature. Therefore, we perform another fit using the OPLIB opacity table, including turbulent mixing with a coefficient of D_{turb} = 2000 cm^{2} s^{−1}, and adding the effective temperature to the set of constraints. We do not include the turbulent mixing coefficient into the free parameters as Figs. 5 and 6 clearly demonstrate that it saturates. The set of constraints is now composed of our seismic indicators and the effective temperature and the free parameters are the age, mass, and initial composition. We are now able to get a suitable model which, again, accounts for the position in the HR diagram but also for the spectroscopic metallicity, which was not required. This is illustrated as a cyan square in the figures, with the label “OPLIB − D_{turb}”. This shows that modified opacities could help model the 16 Cyg A star as accurately as possible and also reinforces the argument that additional mixing processes might be necessary to model both stars. The values of the several stellar parameters are gathered in Table E.7.
3.2.2. Accounting for the metallicity
Up to now, the only two models accounting for the spectroscopic metallicity constraint from Ramírez et al. (2009) are models with a reduced impact of diffusion. Those correspond to the one with the OPLIB opacity table and including turbulent mixing with a fixed coefficient of D_{turb} = 2000 cm^{2} s^{−1}, labelled “OPLIB − D_{turb}” in the figures, for 16 Cyg A and the one without diffusion for 16 Cyg B. Both models reproduce the complete set of seismic and nonseismic constraints (that is Δ, , , A_{He}, T_{eff}, L, and [Fe/H]). However, the spectroscopic metallicity was not yet part of the fitting constraints. Moreover, as Figs. 3 and 4 clearly illustrate, most of the computed models do not agree with the spectroscopic metallicities.
As Figs. 5 and 6 clearly show, the turbulent mixing coefficient saturates and freeing its value cannot enable us to produce models that reproduce the metallicity. We also note that the impact of the mixing length parameter is mostly focused on the effective temperature. As a consequence, we expect a large variation of this parameter will be necessary to reproduce the metallicity. In order to verify this hypothesis, we compute models with a free mixing length parameter and the metallicity as a constraint. The set of free parameters is made of the age, mass, composition and, mixing length of the star and the constraints are the set of seismic constraints and the metallicity. The results are given in Table 5 and shown in grey brown with the label [Fe/H] in Figs. 1 through 4. We indeed observe that the necessary variations in α_{MLT} are incompatible with those to reproduce T_{eff}. Indeed, now that the model and observed metallicities agree, the effective temperature values do not. Trying to include both nonseismic constraints, using the mixing length parameter as a free coefficient and either including turbulent mixing or not did not lead to a satisfactory adjustment (i.e. with a reduced χ^{2} value inferior to 1). This clearly shows that we are not able, with the current set of parameters, to reproduce the complete set of seismic and nonseismic constraints without invoking special physical processes. We also note that the model for 16 Cyg B is both too massive and young compared to other studies (e.g. Buldgen et al. 2016b; Verma et al. 2017). This illustrates that one has to proceed with caution when modelling data as it is possible to provide a model which is representative of these data while having no physical meaning.
Adjusted stellar parameters including the metallicity constraint.
3.3. Individual best models
In the present section, we summarise and further analyse the two best models we obtain while regarding both stars as separate, that is to say without imposing a common initial composition and age. Those are the only models which simultaneously account for seismic and non seismic constraints and are the ones denoted “OPLIB − D_{turb}” for 16 Cyg A and “NoDiff” for 16 Cyg B. Table 6 shows the choice of input physics as well as the set of optimal parameters of those models. As both models do not have the same set of input physics, they may not be regarded as valid candidates to study the system as a whole as is done in Sect. 4. We indeed expect from binary stars with close stellar parameters that their internal physics should overall be identical. The goal of the present section is only to analyse in more depth models which fitted the complete set of considered constraints and to investigate whether those models still may be further improved.
Best individual models for each star, accounting for both seismic and nonseismic constraints.
We display the échelle diagrams of each star in Figs. 8 and 9. We observe that the frequency trend for both stars is well accounted for. However, we note a drift at high frequencies. We expect this effect to mainly result from the surface effects. To illustrate this claim, we display the échelle diagram of optimal models for both stars computed with seismic indicators which are not corrected for surface effects in Figs. D.1 and D.2. We indeed observe that the high frequency drift is reinforced.
Fig. 8. Échelle diagram of 16 Cygni A bestfit model. The crosses are the observed frequencies and the diamonds the theoretical ones. 
Fig. 9. Échelle diagram of 16 Cygni B bestfit model. The crosses are the observed frequencies and the diamonds the theoretical ones. 
We also note that, in the case of 16 Cygni B, the theoretical ridges are shifted with respect to the observed ones. This effect should mainly be due to the seismic indicator defined in Farnir et al. (2019). It corresponds to an estimator of the constant term in n in the asymptotic expression of frequencies as in Gough (1986) and has been shown to be sensitive to the surface effects. Its value along with several other WhoSGlAd indicators for every computed model are displayed in Figs. 10 and 11 (their definitions and observed values are given in Appendix B). We observe that the is not properly accounted for. However, it is worse in the case of 16 Cyg B, which explains why this effect is much more visible.
Fig. 10. Values of complementary seismic indicators for 16 Cyg A. Observed values along with their uncertainties are shown as a black box. 
Fig. 11. Values of complementary seismic indicators for 16 Cyg B. Observed values along with their uncertainties are shown as a black box. 
In Appendix B, we define other seismic indicators that were not part of the constraints. Now looking at those indicators displayed in Figs. 10 and 11, we see that, in most cases, none of them are properly accounted for. In the lower panels, we display the values of the base of the convection zone glitch amplitude and note that only a few models are within the one σ uncertainty region. However, one should not be alarmed by this observation as its value is hardly significant in the case of the 16 Cygni system (A_{CZ} = 2 ± 1 for both stars) and, therefore, bears little information.
We also represent the values for both Δ_{01} and Δ_{02} which represent the slopes of the individual frequency ratios r_{01} and r_{02} expressed as a function of the radial order. Again, every model presents a value which is significantly different from the observed ones. Nonetheless, accounting for such data is a complex task and we note that, in the present situation, only the modification of diffusion seems to provide an improvement for both stars.
We do not include the indicator in the modelling procedure as it has been shown by Farnir et al. (2019, Fig. 4) to be sensitive mostly to the surface effects and highly degenerate in the stellar mass. Moreover, it is tightly correlated to the large separation (see Eq. (B.2) and the asymptotic formulation of the frequencies in Gough 1986). The Δ_{0l} indicators are not used as they mostly carry information about central overshooting (Farnir et al. 2019, Fig. 3) which we presumed would not happen as both stars are below the approximative limit of ∼1.1 M_{⊙} and are expected to have a radiative core.
Finally, we plot in Figs. 12 and 13 the observed individual frequency ratios defined in Roxburgh & Vorontsov (2003), which we recall are not used as constraints in our fits, as a function of frequency against the best model ones. We observe that, although the overall agreement is good, the oscillation which is present in the observed ratios is not properly accounted for in the model frequencies. This clearly indicates that some information remains to be exploited to model the system as comprehensively as we can and inversion techniques may be of great help in doing so. However, with the use of only our indicators, instead of individual ratios, the overall trend seems to be well respected in both cases.
Fig. 12. Individual ratios for 16 Cygni A. Observed values along with their uncertainties are shown as crosses, best model values are represented by diamonds. 
Fig. 13. Individual ratios for 16 Cygni B. Observed values along with their uncertainties are shown as crosses, best model values are represented by diamonds. 
4. The system as a whole
In this section, we select the set of individual models which respect the binarity constraint and try to provide models while imposing a common age and initial composition.
4.1. Accepted models
The models satisfying the binarity constraint are those that have, within one another uncertainties, identical ages and initial composition. Figure 7 provides a clear illustration. Only models represented by a cross that meets the line, corresponding to identical stellar parameters in the three panels, are kept. Those models are referred to as the accepted models and are: the reference models (AGSS09, dark blue), those with turbulent mixing (D_{turb}, purple), without diffusion (NoDiff, light pink), with overshooting (α_{ov}, red), the models with the mixing length coefficient adjusted for the effective temperature (T_{eff}, light green), and with a temperature profile above the photosphere as in Vernazza et al. (1981) with a calibrated α_{MLT} of 2.02 (Vernazza α_{MLT} dark pink). From those models, we define the range of accepted stellar parameters given in Table 7.
Accepted stellar range defined as the centroid of the extremum values for each parameter.
Even though models without diffusion are included in the set of accepted models, we emphasise that it does not mean that microscopic diffusion should not be included when modelling the system but that other mixing processes might occur to counteract it. Moreover, those models largely shift the accepted parameter ranges to heavier, older and hydrogen rich – thus metal and helium poor – models. This significant difference in composition is a clear illustration of the degeneracy between helium and metal abundances in the helium glitch amplitude.
We previously noted that the only models accounting simultaneously for the complete set of seismic and nonseismic constraints were the model with the OPLIB opacity table and turbulent mixing, for 16 Cygni A, and the one without diffusion for 16 Cygni B. This again hints at the necessity to include nonstandard physical processes. However, we must point out that they are incompatible for a joint analysis, as is carried out in the present section, as their ages and compositions are significantly different. Moreover, as we use different opacity tables for both stars, they may not be used simultaneously to analyse the system as a whole.
4.2. Binary models
We now use the individual accepted models to compute models imposing a common age and composition. The adjustment is carried out as in the previous section (that is, with the same set of free parameters and constraints as in Sect. 3) only ages and initial compositions are required to be identical, reducing the set of free parameters by three. We use the average values of the initial composition and ages as initial guesses. The set of free parameters is composed of: one value of the age, initial hydrogen and metal fraction for both stars, and a different value of the mass for each star. The set of constraints corresponds to the individual values of the 4 seismic constraints considered in this paper (Δ, , , and A_{He}).
We are not able to provide an exact adjustment of both stars simultaneously (“exact” meaning that the reduced χ^{2} value should not exceed a value of 1). This may result from the fact that the size of the parameters space is reduced by three. Effectively, even though the accepted ages and initial compositions agree within their uncertainties (see Table 7), they are not identical and our seismic constraints may be too stringent to allow for an exact simultaneous fit while imposing identical ages and compositions. To illustrate this statement, we have plotted the optimal ages and initial composition of one stellar component against the other for each choice of input physics in Fig. 7. We observe in this figure that the three common free parameters almost never simultaneously agree for a given choice of input physics. As further illustration, we may compute the relative difference in the initial metallicity between both stars normalised by the quadratic sum of their uncertainties (). In the most favourable case, we obtain 0.2, while we get 2.7 in the least favourable one. This shows that the difference in initial composition for individual models is sometimes significant and may impair the convergence of the LevenbergMarquardt procedure.
One could argue that our inability to provide an exact adjustment for both stars comes from the reduced number of free parameters. Therefore, we try to include the mixing length parameter into the fitting parameters, allowing it to vary freely and independently for each star. However, this does not improve the results. As a matter of fact, the mixing length parameter value does not significantly vary. We retrieve optimal values of α_{MLT, A} = 1.8 ± 0.2 and α_{MLT, B} = 1.8 ± 0.1 respectively, compared to the fixed solar value of α_{MLT} = 1.82. From Magic et al. (2015) and the solartwin character of both stars (see also the discussion in Sect. 3.2.1) one might expect that the mixing length parameter value should remain close to solar.
Although we do not obtain models that exactly reproduce the seismic constraints, in some cases, we may find a reasonable agreement with most of them (but not all, therefore having a reduced χ^{2} value greater than 1). We obtain two sets of convincing results: one for models without diffusion, the other for models with a temperature profile above the photosphere as in Vernazza et al. (1981) and a corresponding solar calibrated α_{MLT} = 2.02 value. The optimal model stellar parameters are gathered in Tables 8 and 9, respectively, and Tables 10 and 11 show the differences between the observed and model seismic indicators normalised to the observed uncertainties. We also show the complete set of “optimal” model parameters, for each choice of input physics, as well as their relative agreement with the seismic constraints in Tables E.8 and E.9. We observe for the models without diffusion in Table 10 that, out of the 8 seismic constraints, only the large separation of 16 Cygni A was not properly accounted for. All the other indicators are within the 1σ uncertainty. For the models with a temperature profile above the photosphere as in Vernazza et al. (1981), both the large separation and small frequency ratio between radial and dipolar modes are poorly reproduced.
Best set of adjusted stellar parameters, models without diffusion, imposing a common age and initial composition for both stars.
Best set of adjusted stellar parameters, models with a temperature profile above the photosphere as in Vernazza et al. (1981) and with α_{MLT} = 2.02, imposing a common age and initial composition for both stars.
Differences between theoretical, model without diffusion, and observed values for the seismic constraints defined as δ = I_{obs} − I_{th}, in the units of the constraint.
Same as in Table 10, but for models model with a temperature profile above the photosphere as in Vernazza et al. (1981) and α_{MLT} = 2.02.
For the other models, presented in Tables E.8 and E.9, we note that it is always the large separation of the A star which is poorly fitted. The other indicators are rather well adjusted but the small separation ratios often fall out, while remaining close, of the 1σ uncertainty box. This difference in fitting between the several indicators mainly stems from the difference in their relative uncertainties. Indeed, the helium glitch amplitudes have relative uncertainties of about 3% and, by far, are the least stringent constraints. Then come the ratios with relative uncertainties around 0.8%. Finally, the Δ constraint relative uncertainties are of approximately 0.004%. Moreover, it was shown by Farnir et al. (2019) that the ratio is mostly sensitive to the evolutionary state of the star, therefore at a given composition and mass, to the stellar age. Then, the large separation is a proxy of the mean stellar density and decreases along the main sequence. As both stars are required to have identical ages and compositions, and with such stringent constraints, we may understand that only one large separation may be fitted at a time and that all the other constraints, from the clear imbalance in the relative uncertainties will adjust to it. To find a simultaneous agreement for those models that did not reach satisfactory convergence, we either need to relax this assumption, allowing for example different values for the initial composition, or to relax the seismic constraints.
We also display in Fig. 14 the agreement of the models with the nonseismic data for each of the considered variations in input physics, represented by the different symbols. The observed values along with their uncertainties are shown as boxes. We display the results for 16 Cyg A in blue and for 16 Cyg B in red. We note, as in the previous section, that the models for each star are almost constant in radius. We also note that very few models account for the position of the stars in the HR diagram. Actually, those are the models which, individually accounted for these data. It does not come as a surprise as the minimisation aims at finding a compromise between all the seismic constraints for both stars. Then, looking at the lower panel, we note that no model for 16 Cyg A is representative of the surface composition. For 16 Cyg B, only the model without diffusion agrees with these data – again with no surprise as the individual model already agreed – Let us add that, as the models without diffusion for both stars must have the same initial composition, their surface compositions are identical and both markers are indistinguishable in the lower panel of the figure.
Fig. 14. Computed values of nonseismic constraints against the observed ones, symbolised by a box, for the system seen as a whole (i.e. with common ages and initial compositions). Each variation in input physics is represented by a different symbol. The colour represents the star, blue for 16 Cyg A and red for 16 Cyg B. The upper panel is a HR diagram. The lower panel shows the surface helium abundance versus the metallicity. 
In a nutshell, we are able to produce two pairs of binary models that are in reasonable agreement with our seismic constraints. However, no model accounts simultaneously for the seismic and nonseismic data of both stars. This may result from a differential effect between both stellar components which could, for example, create differences in compositions as has been discussed by Maia et al. (2019). This would require the inclusion of nonstandard physical processes in the modelling.
4.3. Relaxing the common composition hypothesis
We noted in the previous section that simultaneously providing models of both stars while imposing them to have identical initial compositions as well as ages is a difficult task. In most cases, this resulted in our inability to build such models. Therefore, we try to model the system requiring only an identical age for both stars. The set of free parameters is thus made of the age, the individual initial hydrogen, and metal abundances and the individual masses. This adds up to 7 free parameters. The constraints are the seismic indicators used throughout this paper which represent 8 constraints.
We compute models for the several sets of input physics considered in the previous section. The individual stellar parameters of those models are given in Table E.10. To quantify whether the improvement of the results is significant given the increased number of parameters, we introduce the Bayesian Information Criterion (BIC; Schwarz 1978). It allows to compare models of different dimensionalities and provides a criterion for making a selection. It has the advantage over the simple χ^{2} value to penalise over the number of fitting parameters and may pinpoint overfitting models. Under the assumption that model errors are independent and normally distributed, it takes the form:
with the χ^{2} value as defined previously, k the number of free parameters and N the number of constraints. When comparing different models, the key ingredient is not the BIC value itself but rather the difference between values – one must however keep in mind that the lower the value, the better – For the BIC difference to be significant, it has to exceed 2. BIC differences from 6 and above will be regarded as strongly significant (Kass & Raftery 1995). In this specific case, we compare models from the previous section that had 5 free parameters and 8 constraints, the added term to the χ^{2} value thus equals 10.4. In the current section, it is equal to 14.6. The difference is 4.2 which the χ^{2} improvement has to exceed so that the model may be regarded as improved.
Looking at the models in Table E.10, we observe that the relative difference in compositions between the two stars is very small (the maximum difference reaches up to 0.07 dex in X_{0} and 0.002 dex in (Z/X)_{0}). However, that variation alone is able to greatly improve the BIC values, as the comparison of Tables E.9 and E.11 shows. Only two models were not significantly improved. The model considering the mixing length parameters calibrated over the effective temperature of the stars could not be improved. Also, the one using the temperature profile of Vernazza et al. (1981) above the photosphere and a calibrated value of α_{MLT} is not improved, BICwise (the raw χ^{2} value did decrease of about 1.7). The BIC variation is of about 2 which makes this difference barely relevant. Therefore, in most cases, it is relevant to allow the composition to slightly vary between the two stars. This again points toward the necessity to include nonstandard physical processes. We also note that it is often the case that the uncertainties on the individual stellar parameters are degraded. Finally, we note that, in most cases, 16 Cyg B is initially richer in metals and poorer in hydrogen than its twin. This validates the trend we observe in the middle and lower panels of Fig. 7. Therefore, we observe that, seismically speaking, the B component of the system is more metallic than the A, which is in opposition with the spectroscopic observations (see Table 1).
5. Discussion
5.1. General considerations about the approach
Compared to other asteroseismic approaches aiming at the modelling of solarlike pulsators, the greatest difference of our approach is that we do not use directly the complete set of individual frequencies or of individual frequency ratios as constraints. Instead, we build seismic indicators, via the WhoSGlAd method, which are as little correlated as possible and relevant to the stellar structure. Moreover, the search for optimal models is carried out by minimising a single cost function comparing simultaneously theoretical seismic and nonseismic constraints to observed data. This has the advantage of avoiding unnecessary correlations.
Furthermore, the direct use of the helium glitch amplitude in our modelling also makes up the peculiarity of our approach. Indeed, in several other studies, it is not used directly as a constraint to the modelling. For example, Verma et al. (2014, 2019) calibrate the model helium glitch amplitude with respect to the surface helium abundance in a set of optimal models representative of other constraints (namely individual frequencies, ratios, effective temperature and metallicity) to provide an estimate. In the present case, including the helium glitch amplitude as a constraint to the fit acts as a constraint on the model helium abundance, with some correlation with the metal content as showed in Sect. 3.1.1. This means that we do not assume a specific relation between the two quantities and the resulting helium abundance stems from the best model search only.
We noted that providing models of both stars while requiring a common age and composition proves to be an arduous task. For most choices of micro and macrophysics, we are not able to produce satisfactory models. However, when completely inhibiting the microscopic diffusion of chemical elements or imposing a temperature relation above the photosphere following the prescription of Vernazza et al. (1981) while using a specifically calibrated value of the mixing length parameter, we obtained a reasonable agreement, but did not go below the value of 1 for the reduced χ^{2}. On the first hand, the fact that including extra mixing produces better results illustrates the need to include nonstandard physical processes to properly model such complex data. On the other hand, the improvement of the results while using another temperature profile above the photosphere demonstrates the impact of the surface effects on the seismic indicators. Even though those indicators were defined in such a way to lessen this effect at most. We must also add that nonseismic constraints are not simultaneously accounted for in both stars (see Fig. 14). This clearly illustrates the need to include such constraints in the fitting procedure as well as additional physical processes. Some of the possible nonstandard physical processes that could be included in the modelling are discussed in Sect. 5.3. Another way to improve models of the system as a whole is to relax the hypothesis of a common composition. This leads to small differences in composition between both stars (never exceeding 0.07 dex in X_{0} and 0.002 dex in (Z/X)_{0}). This is however sufficient in many cases to improve the results significantly as the BIC values testify.
5.2. Impact of the surface effects
The impact of the surface effects correction of the frequencies on the optimal model is not clear. We computed models that adjusted seismic indicators built over the uncorrected frequencies. Those are displayed in dark green in the figures throughout the paper. Although it is barely significant we note that those models are heavier and older than models with corrected frequencies.
Furthermore, we noted in Sect. 4.2 that the large separation we compute is the most stringent of our indicators with a relative uncertainty of approximately 0.004%. However, this specific value of the uncertainty, resulting from the error propagation of the individual frequencies, is unrealistic. To provide a more robust estimate, we can quantify the contribution of the surface effects by computing the difference between corrected and uncorrected values. We obtain σ(Δ_{A}) = 0.9 μHz and σ(Δ_{B}) = 1.0 μHz, for 16 Cyg A and B respectively, which amounts to 1.1% relative uncertainty. This shows that, even though we build seismic indicators in such way that makes them the less dependent on the surface effects as possible, they are still impacted by the surface conditions.
The impact of the surface effects correction is further illustrated in Appendix D where we give the values of the seismic indicators for uncorrected frequencies and display the bestfit models échelle diagrams of both stars. We also compute models with frequencies corrected according to the prescription of Ball & Gizon (2014). The optimal parameters are given in Table D.2. We observe, compared to the reference models, that both stars become older and hydrogen rich. We also note that, while the A component becomes more metallic, the B one is then less metallic. 16 Cygni B also becomes heavier while it is not the case for its twin. Nonetheless, the differences are such that we could include those models in the set of accepted ones.
5.3. Nonstandard physical processes and modelling improvements
As was shown by Verma et al. (2019) and as we illustrate in Sect. 3.1.3, diffusion is very important in the modelling of solarlike stars. However, the models computed with CLES currently consider only three chemical species, the hydrogen, the helium and the metals. A significant improvement would stem from the consideration of several subspecies in the metals. For example, we may retrieve invaluable information by following the lithium evolution with the star. Indeed, the discussion regarding the effects of potential nonstandard processes is tightly linked to the lithium abundance in both stars and its connection with the formation and orbital evolution of planetary systems. Deal et al. (2015) have proposed that, since the B component was orbited by a Jovian planet (Cochran et al. 1997), accretion of matter from the planetary disc in the envelope of 16CygB would have triggered fingering convection and thus led to a strong decrease in the lithium abundance. They determined that the accretion of 0.66M_{⊕} would be enough to reproduce the lithium abundance of 16CygB that is 4 times lower than that of 16CygA (Friel et al. 1993; King et al. 1997), despite both stars having very similar rotational and structural properties. It is also interesting to note that, in the broader context of the Li abundance of solartwins, 16CygA seems to be more Lirich than similar solartwins, while 16CygB seems to follow the trend observed with age for these stars (Carlos et al. 2016). This could suggest a reverse scenario that should explain the high lithium abundance of 16CygA and in particular the possibility of an increase in lithium abundance related to planet engulfment (e.g. Montalbán & Rebolo 2002; Carlos et al. 2016). In this context, analysing both scenarios in light of potential traces of such events in seismic indicators may lead to new synergies between asteroseismology and exoplanetology, namely in the analysis of planetary formation and material accretion onto the surface of planethost stars.
Furthermore, our models do not include rotation, nor rotationinduced mixing. As a matter of fact, Davies et al. (2015) and Bazot et al. (2019) showed that rotation indeed occurs in both stars taking advantage of the rotational splitting present in stellar oscillation spectra of rotating stars. Bazot et al. (2019) even showed that differential rotation occurs in both stars, in a similar way to our Sun. Such a process could significantly affect the diffusion of chemical elements. This could be argued to be a flaw in our models. However, the additional mixing induced may be approximated by the turbulent mixing of elements as was performed for several models in this study. This prescription consists of an approximation and the models may still be improved. Therefore, we may, in the future, use improved models which account for the rotation. This will be discussed in further papers of the series.
Finally, we noted that choosing a different opacity profile, that of the Los Alamos project, while including turbulent mixing of elements, which counteracts diffusion, allowed to reproduce both seismic and non seismic constraints for 16 Cygni A. This provides clues that we may need to modify the opacity profile of the star to properly account for all the observed constraints. Therefore, inversion techniques could help us to further improve our models. However, the OPLIB opacity table has two different effects on the stellar structure. On the one hand, it modifies slightly the size of the convective envelope. On the other hand, it changes the temperature gradient in the central regions. According to which effects dominates this could be the illustration of the need of nonstandard mixing processes as well.
6. Conclusion
With the aim of characterising the 16 Cygni system as thoroughly as possible, we took advantage of the seismic indicators defined via the WhoSGlAd method to provide stringent constraints on the stellar structure and test several choices of micro and macrophysics. We built those indicators using the frequencies computed over the full length of the Kepler data by Davies et al. (2015) and corrected for the surface effects according to Kjeldsen et al. (2008)’s power law adjusted by Sonoi et al. (2015). The several choices of micro and macrophysics used in stellar models we tested are: the solar reference mixture, opacity and equation of state tables as well as the inclusion of turbulent mixing or of diffusion of chemical elements, a different choice for the mixing length parameter, the inclusion of overshooting outside of convective regions, a different choice of temperature profile above the photosphere or the effect of the correction of the frequencies for the surface effects.
Overall, our results agree with previous studies with slight differences according to the choice of physics included in the models. We showed that the use of the WhoSGlAd indicators allows to discriminate between several of those choices. However, we also note that those indicators alone do not suffice to provide a complete adjustment of the stars as, in most cases, the nonseismic constraints (i.e. luminosity, effective temperature, and metallicity) are not satisfied. Therefore, they need to be included in the fitting process to provide the most representative model. We show in Sect. 3.2.1 that the mixing length parameter has a clear effect on the modelled effective temperature of the star. Indeed, a value greater than the solar one allows to greatly improve the agreement with the observed value in the case of 16 Cyg B. Moreover, we also study the impact of the inclusion of turbulent mixing and show that it leads to a better fit of the effective temperature and metallicity. However, we observe that the turbulent mixing coefficient saturates and using it as a free parameter of the modelling procedure would be meaningless. Therefore, using a free α_{MLT} and T_{eff} as a constraint, we are able to produce models in agreement with this constraint. However, the observed metallicities are not reproduced by those models. We also demonstrate in Sect. 3.2.2 that varying the mixing length parameter while using the metallicity constraint allows to better reproduce its value for both stars but at the cost of the agreement with the observed effective temperatures. This illustrates the necessity to build more complex models in order to reproduce both seismic and nonseismic constraints.
Indeed, we show that, to reproduce the nonseismic constraints, we have to select only specific choices of input physics or even include nonstandard physical processes. Indeed, for 16 Cygni A, only a model with a modified opacity profile, from the Los Alamos Opacity Project, and including turbulent mixing of the chemical elements reproduced all constraints. For 16 Cygni B, it was a model that did not include diffusion that was able to account for these constraints. This illustrates in both cases that nonstandard physical processes may be necessary to inhibit diffusion and to properly models those stars. Such processes could be the accretion of planetary matter or rotationinduced mixing.
Adjusting both stars simultaneously while imposing a common age and composition proved to be a difficult task. In most cases, with models that were consistent within mutual uncertainties for both stars as initial guesses, we could not obtain a satisfactory adjustment. The large separation of 16 Cyg B, because of the high precision of the Δ constrain defined in our study, was dominant and the free parameters adapted at best to provide a compromise for the other constraints. This resulted in models that did not exactly fit the large separation of 16 Cygni A while reasonably fitting other constraints (only a few σ difference from the required value). Nonetheless, for models without diffusion or with a different temperature profile above the photosphere (that of Vernazza et al. 1981) and a calibrated value of α_{MLT}, we were able to account for the seismic constraints of both stars with a reduced χ^{2} of 3.4 and 2.9 respectively. The difficulty to provide satisfactory models of both stars with other choices of input physics indicates that it can be necessary to either relax the common initial composition assumption, the seismic constraints or to invoke special physical processes. For example, we showed in Sect. 4.2 that the differences in the initial metallicity between both stellar components of optimal individual models may sometimes be significant. Therefore, we computed models relaxing relaxing the common composition hypothesis in Sect. 4.3 and were able to significantly improve the results. We observe that a small difference between the initial compositions of both stars is sufficient.
With the aim of providing a broad sample of reliable models of the system, the extensive analysis of the degeneracies carried out by combining seismic and nonseismic constraints is of prime importance to fully grasp the uncertainties of inverse analysis but also to the extent in which we can constrain physical processes not implemented in standard stellar models linked for example to the effects of accretion of planetary matter, angular momentum transport and their link to both seismic indices and the lithium and beryllium abundances of both stars.
We showed that, even for our models that reproduced both seismic and nonseismic constraints, information remain to be analysed as we observe that other indicators defined in Farnir et al. (2019) are not properly represented (i.e. , Δ_{01}, Δ_{02}, and A_{CZ}). Further studies could focus on those other constraints. Finally, as the WhoSGlAd method also proves to provide very stringent seismic constraints, we will, in future studies, undertake the adjustment of the Kepler LEGACY sample (Lund et al. 2017) to try and retrieve global trends in solarlike oscillators. This data set contains the best set of solarlike oscillation spectra available to the community to this day as it is composed of 66 solarlike stars which have been continuously observed from space for at least one year. Therefore, we would be able to realise an ensemble study of stellar parameters. For example, we could study the evolution of the amount of central overshooting with the stellar mass.
Acknowledgments
The authors would like to thank the referee for their careful reading of the paper as well as for their very constructive remarks. M. F. is supported by the FRIA (Fond pour la Recherche en Industrie et Agriculture) – FNRS PhD Grant. G. B. acknowledges fundings from the SNF AMBIZIONE Grant No. 185805 (Seismic inversions and modelling of transport processes in stars). C.Pinçon is supported by the F. R. S – FNRS as a Chargé de Recherche. P. E. and S. J. A. J. S. have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 833925, project STAREX). C.Pezzotti is sponsored by the Swiss National Science Foundation (project number 200020− 172505).
References
 Adelberger, E. G., García, A., Robertson, R. G. H., et al. 2011, Rev. Mod. Phys., 83, 195 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Badnell, N. R., Bautista, M. A., Butler, K., et al. 2005, MNRAS, 360, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2009, in Transiting Planets, eds. F. Pont, D. Sasselov, M. J. Holman, et al., IAU Symp., 253, 71 [Google Scholar]
 Ball, W. H., & Gizon, L. 2014, A&A, 568, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Basu, S., Mazumdar, A., Antia, H. M., & Demarque, P. 2004, MNRAS, 350, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Bazot, M., Benomar, O., ChristensenDalsgaard, J., et al. 2019, A&A, 623, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, in American Astronomical Society Meeting Abstracts #215, Bull. Am. Astron. Soc., 42, 215 [Google Scholar]
 Buldgen, G., Salmon, S. J. A. J., Reese, D. R., & Dupret, M. A. 2016a, A&A, 596, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Buldgen, G., Reese, D. R., & Dupret, M. A. 2016b, A&A, 585, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carlos, M., Nissen, P. E., & Meléndez, J. 2016, A&A, 587, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cassisi, S., Salaris, M., & Irwin, A. W. 2003, ApJ, 588, 862 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., & Daeppen, W. 1992, A&A Rev., 4, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Cochran, W. D., Hatzes, A. P., Butler, R. P., & Marcy, G. W. 1997, ApJ, 483, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Colgan, J., Kilcrease, D. P., Magee, N. H., et al. 2016, ApJ, 817, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York: Gordon and Breach) [Google Scholar]
 Davies, G. R., Chaplin, W. J., Farr, W. M., et al. 2015, MNRAS, 446, 2959 [NASA ADS] [CrossRef] [Google Scholar]
 Deal, M., Richard, O., & Vauclair, S. 2015, A&A, 584, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., Miglio, A., Montalban, J., et al. 2010, A&A, 509, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., Buldgen, G., & Salmon, S. J. A. J. 2019, A&A, 626, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Farnir, M., Dupret, M. A., Salmon, S. J. A. J., Noels, A., & Buldgen, G. 2019, A&A, 622, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Friel, E., Cayrel de Strobel, G., Chmielewski, Y., et al. 1993, A&A, 274, 825 [NASA ADS] [Google Scholar]
 Gough, D. O. 1986, Highlights Astron., 7, 283 [CrossRef] [Google Scholar]
 Gough, D. O. 1990, in Progress of Seismology of the Sun and Stars, eds. Y. Osaki, & H. Shibahashi (Berlin Springer Verlag), Lect. Notes Phys., 367, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse, 15 [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [CrossRef] [Google Scholar]
 King, J. R., Deliyannis, C. P., Hiltgen, D. D., et al. 1997, AJ, 113, 1871 [NASA ADS] [CrossRef] [Google Scholar]
 Kjeldsen, H., & Bedding, T. R. 1995, A&A, 293, 87 [NASA ADS] [Google Scholar]
 Kjeldsen, H., Bedding, T. R., & ChristensenDalsgaard, J. 2008, ApJ, 683, L175 [NASA ADS] [CrossRef] [Google Scholar]
 Lund, M. N., Silva Aguirre, V., Davies, G. R., et al. 2017, ApJ, 835, 172 [CrossRef] [Google Scholar]
 Magic, Z., Weiss, A., & Asplund, M. 2015, A&A, 573, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maia, M. T., Meléndez, J., LorenzoOliveira, D., Spina, L., & Jofré, P. 2019, A&A, 628, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manchon, L., Belkacem, K., Samadi, R., et al. 2018, A&A, 620, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Metcalfe, T. S., Chaplin, W. J., Appourchaux, T., et al. 2012, ApJ, 748, L10 [Google Scholar]
 Montalbán, J., & Rebolo, R. 2002, A&A, 386, 1039 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Monteiro, M. J. P. F. G., ChristensenDalsgaard, J., & Thompson, M. J. 2000, MNRAS, 316, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Proffitt, C. R., & Michaud, G. 1991, ApJ, 380, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Ramírez, I., Meléndez, J., & Asplund, M. 2009, A&A, 508, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reese, D. R., Marques, J. P., Goupil, M. J., Thompson, M. J., & Deheuvels, S. 2012, A&A, 539, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
 Roxburgh, I. W., & Vorontsov, S. V. 2003, A&A, 411, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
 Scuflaire, R., Théado, S., Montalbán, J., et al. 2008a, Ap&SS, 316, 83 [Google Scholar]
 Scuflaire, R., Montalbán, J., Théado, S., et al. 2008b, Ap&SS, 316, 149 [Google Scholar]
 Sonoi, T., Samadi, R., Belkacem, K., et al. 2015, A&A, 583, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tassoul, M. 1980, ApJS, 43, 469 [Google Scholar]
 Thévenin, F., Oreshina, A. V., Baturin, V. A., et al. 2017, A&A, 598, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thoul, A. A., Bahcall, J. N., & Loeb, A. 1994, ApJ, 421, 828 [NASA ADS] [CrossRef] [Google Scholar]
 Ulrich, R. K. 1986, ApJ, 306, L37 [Google Scholar]
 Verma, K., Faria, J. P., Antia, H. M., et al. 2014, ApJ, 790, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Verma, K., Raodeo, K., Antia, H. M., et al. 2017, ApJ, 837, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Verma, K., Raodeo, K., Basu, S., et al. 2019, MNRAS, 483, 4678 [NASA ADS] [CrossRef] [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [Google Scholar]
 Vorontsov, S. V. 1988, in Advances in Helio and Asteroseismology, eds. J. ChristensenDalsgaard, & S. Frandsen, IAU Symp., 123, 151 [NASA ADS] [CrossRef] [Google Scholar]
 White, T. R., Huber, D., Maestro, V., et al. 2013, MNRAS, 433, 1262 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: WhoSGlAd decomposition
We describe here the basis of functions which are used to represent the oscillation spectrum of a solarlike pulsator. First, the smooth part of the spectrum is represented by a secondorder polynomial of n, the radial order. We thus have the following succession of polynomials:
with p_{k}(n) = n^{k}, k = 0, 1, 2 and δ_{ll′} the Kronecker delta comparing two values of the spherical degree l and l′.
Then, the helium glitch is described by the following oscillating functions:
where τ is the acoustic depth of the glitch, Δν the asymptotic large frequency separation and . is actually the first order approximation of ν_{l, n}. The asymptotic large frequency separation is defined as (Tassoul 1980), with the local radius r, the local sound speed c(r) and R_{*} the radius of the star at the photosphere.
We must add that the values of τ and Δν are estimated via a model that is representative of the seismic indicators of the smooth part, namely Δ, and . Farnir et al. (2019) showed that the exact value of τΔν has a negligible impact on the amplitude of the glitch. A 10% percent excursion from the optimal value is of negligible impact.
Finally, orthonormalisation of the basis function is carried out via GramSchmidt’s process. This produces the orthonormal elements over which we may project the frequencies to represent them. We thus retrieve completely independent coefficients.
Appendix B: Additional seismic indicators
In the present section, we describe supplementary seismic indicators defined in the WhoSGlAd method but that are not used directly in our optimisation for a model representative of the 16 Cygni system. The values of those indicators are given in Table B.1.
Additional observed seismic indicators.
: In taking inspiration in the asymptotic formulation of the frequencies (Gough 1986)
we may construct a vector subspace over which frequencies are represented by the function:
where and K are free parameters. By defining an orthonormal basis over this subspace, projecting the frequencies and identifying the several coefficients with the asymptotic expression we may get an expression for .
A_{CZ}.
where δν_{CZ} is the base of the convection zone glitch component.
Δ_{0l}. Corresponds to the slope of the individual frequency ratios r_{0l} as a function of the radial order n and is defined as:
Appendix C: Impact of high uncertainties modes
From the modes computed by Davies et al. (2015) we select those with uncertainties below 1.5 μHz. Those high uncertainty modes have a limited and negligible impact on the results as our indicators are averaged over the whole spectrum. Furthermore, as the high frequency modes are the ones which are the most affected by the surface effects this may render our results more robust. In Table C.1 we show the optimal set of stellar parameters retrieved when using the full set of frequencies. We observe that the results do not vary significantly from the case considering modes with uncertainties below the 1.5 μHz threshold, presented in Table E.3 under the label AGSS09. Only the uncertainties on the individual parameters are affected which should not be of any concern as the dominant factor remains the choice of input physics. As further validation of this choice, we display in Figs. C.1 and C.2 a comparison of the échelle diagram of optimal models of both stars using the complete set of frequencies (blue circle) with the ones with the reduced set (red diamond) and the observations (black crosses). We observe that the results do not significantly differ, only that the high frequency drift is more visible as more frequencies are displayed.
Fig. C.1. Échelle diagram of 16 Cygni A comparing optimal reference model with full set of frequencies (blue circles) from Davies et al. (2015; black crosses) with the optimal model with a set restricted to frequencies with uncertainties lower than 1.5 μHz (red diamonds). 
Fig. C.2. Échelle diagram of 16 Cygni B comparing optimal reference model with full set of frequencies (blue circles) from (Davies et al. 2015; black crosses) with the optimal model with a set restricted to frequencies with uncertainties lower than 1.5 μHz (red diamonds). 
Stellar parameters retrieved with the reference set of input physics and the complete set of frequencies.
Appendix D: Influence of the surface effects
We computed both models for frequencies which were not corrected for the surface effects or with theoretical frequencies corrected as in Ball & Gizon (2014) with the adjusted relation in large separation, effective temperature, surface gravity and opacity of Manchon et al. (2018). Table D.1 displays seismic indicators computed with frequencies which are not corrected for the surface effects. The échelle diagram for the models computed with these set of indicators are displayed in Figs. D.1 and D.2. The most striking feature is the large shift between theoretical and observed ridges in both figures which shows that the indicator is not well accounted for as Figs. 10 and 11 show. Moreover, the overall shape of the individual ridges is well represented by the theoretical frequencies. Table D.2 displays the set of optimal parameters for models using frequencies corrected as in Ball & Gizon (2014).
Fig. D.1. Échelle diagram of 16 Cygni A optimal model calculated with seismic indicators defined over frequencies which are not corrected for the surface effects. The crosses are the observed frequencies and the diamonds the theoretical ones. 
Fig. D.2. Échelle diagram of 16 Cygni B optimal model calculated with seismic indicators defined over frequencies which are not corrected for the surface effects. The crosses are the observed frequencies and the diamonds the theoretical ones. 
Observed seismic indicators with frequencies uncorrected for surface effects.
Adjusted stellar parameters with theoretical frequencies corrected as in Ball & Gizon (2014) with the adjusted relation in Manchon et al. (2018).
Appendix E: Individual models
In this section, we summarise the set of input physics used in the reference models and gather the individual stellar parameters as well as the uncertainties propagated during the LevenbergMarquardt adjustment for each model presented in this paper. Table E.1 presents the set of input physics used in the reference model while Table E.2 summarises the several variations of input physics considered throughout the paper. In the latter, the first column gives the label given to the models considering that specific choice of input physics, the second column is the physical ingredient which is varied upon, the third column is the corresponding value and columns 4 and 5 display the χ^{2} values obtained for the optimal models of both stars in each case. Finally, Tables E.3 through E.11 give the complete set of stellar parameters obtained for every case considered in the present paper.
Summary of the physical ingredients included in the reference models, denoted AGSS09.
Variations in the input physics, corresponding name and reduced χ^{2} values.
Summary of the fitted models with only the seismic constraints.
Adjusted stellar parameters using the temperature profile from Vernazza et al. (1981) with a solar calibrated value of α_{MLT} = 2.02.
Results of the modelling considering only seismic constraints with different values of the mixing length coefficient.
Results of the modelling considering only seismic constraints and including turbulent mixing with different values of the turbulent mixing coefficient.
Adjusted stellar parameters including the effective temperature constraint with the OPLIB opacity table and with a fixed turbulent mixing coefficient.
Adjusted stellar parameters imposing a common age and initial composition for both stars for the different variations in physics.
Differences between theoretical and observed values for the seismic constraints, for the several variations in physics, defined as δ = I_{obs} − I_{th}, in the units of the constraint.
Same as Table E.8, but without imposing a common composition for the two stars.
Same as Table E.9, but without imposing a common composition for the two stars.
All Tables
Adjusted stellar parameters including the T_{eff} with the reference set of input physics.
Adjusted stellar parameters including the T_{eff} constraint and turbulent mixing with a coefficient of D_{turb} = 7500 cm^{2} s^{−1}.
Best individual models for each star, accounting for both seismic and nonseismic constraints.
Accepted stellar range defined as the centroid of the extremum values for each parameter.
Best set of adjusted stellar parameters, models without diffusion, imposing a common age and initial composition for both stars.
Best set of adjusted stellar parameters, models with a temperature profile above the photosphere as in Vernazza et al. (1981) and with α_{MLT} = 2.02, imposing a common age and initial composition for both stars.
Differences between theoretical, model without diffusion, and observed values for the seismic constraints defined as δ = I_{obs} − I_{th}, in the units of the constraint.
Same as in Table 10, but for models model with a temperature profile above the photosphere as in Vernazza et al. (1981) and α_{MLT} = 2.02.
Stellar parameters retrieved with the reference set of input physics and the complete set of frequencies.
Adjusted stellar parameters with theoretical frequencies corrected as in Ball & Gizon (2014) with the adjusted relation in Manchon et al. (2018).
Summary of the physical ingredients included in the reference models, denoted AGSS09.
Adjusted stellar parameters using the temperature profile from Vernazza et al. (1981) with a solar calibrated value of α_{MLT} = 2.02.
Results of the modelling considering only seismic constraints with different values of the mixing length coefficient.
Results of the modelling considering only seismic constraints and including turbulent mixing with different values of the turbulent mixing coefficient.
Adjusted stellar parameters including the effective temperature constraint with the OPLIB opacity table and with a fixed turbulent mixing coefficient.
Adjusted stellar parameters imposing a common age and initial composition for both stars for the different variations in physics.
Differences between theoretical and observed values for the seismic constraints, for the several variations in physics, defined as δ = I_{obs} − I_{th}, in the units of the constraint.
Same as Table E.8, but without imposing a common composition for the two stars.
Same as Table E.9, but without imposing a common composition for the two stars.
All Figures
Fig. 1. Summary of 16 Cyg A bestfit models represented in a Mass – Age diagram (top panel), HR diagram (middle panel) and initial hydrogen abundance versus metal composition diagram (bottom panel). The luminosity and effective temperature constraints from White et al. (2013) are represented in the HR diagram as a black box. 

In the text 
Fig. 2. Summary of 16 Cyg B bestfit models. The colours are the same as in Fig. 1. 

In the text 
Fig. 3. Surface helium abundance versus metallicity for 16 Cygni A. The black box represents the spectroscopic metallicity value computed by Ramírez et al. (2009) and the surface helium value from Verma et al. (2014) along with the corresponding uncertainties. 

In the text 
Fig. 4. Surface helium abundance versus metallicity for 16 Cygni B. The black box represents the spectroscopic metallicity value computed by Ramírez et al. (2009) and the surface helium value from Verma et al. (2014) along with the corresponding uncertainties. 

In the text 
Fig. 5. Variation of the optimal metallicity and effective temperature for 16 Cyg A with the mixing length parameter and turbulent mixing. Models of different α_{MLT} values are connected in blue while models with various turbulent mixing are connected by a red line. 

In the text 
Fig. 6. Variation of the optimal metallicity and effective temperature for 16 Cyg B with the mixing length parameter and turbulent mixing. Models of different α_{MLT} values are connected in blue while models with various turbulent mixing are connected by a red line. 

In the text 
Fig. 7. Comparison of adjusted stellar parameters for both stars at a given physics. The straight line shows the locus of identical stellar parameters for both stars. 

In the text 
Fig. 8. Échelle diagram of 16 Cygni A bestfit model. The crosses are the observed frequencies and the diamonds the theoretical ones. 

In the text 
Fig. 9. Échelle diagram of 16 Cygni B bestfit model. The crosses are the observed frequencies and the diamonds the theoretical ones. 

In the text 
Fig. 10. Values of complementary seismic indicators for 16 Cyg A. Observed values along with their uncertainties are shown as a black box. 

In the text 
Fig. 11. Values of complementary seismic indicators for 16 Cyg B. Observed values along with their uncertainties are shown as a black box. 

In the text 
Fig. 12. Individual ratios for 16 Cygni A. Observed values along with their uncertainties are shown as crosses, best model values are represented by diamonds. 

In the text 
Fig. 13. Individual ratios for 16 Cygni B. Observed values along with their uncertainties are shown as crosses, best model values are represented by diamonds. 

In the text 
Fig. 14. Computed values of nonseismic constraints against the observed ones, symbolised by a box, for the system seen as a whole (i.e. with common ages and initial compositions). Each variation in input physics is represented by a different symbol. The colour represents the star, blue for 16 Cyg A and red for 16 Cyg B. The upper panel is a HR diagram. The lower panel shows the surface helium abundance versus the metallicity. 

In the text 
Fig. C.1. Échelle diagram of 16 Cygni A comparing optimal reference model with full set of frequencies (blue circles) from Davies et al. (2015; black crosses) with the optimal model with a set restricted to frequencies with uncertainties lower than 1.5 μHz (red diamonds). 

In the text 
Fig. C.2. Échelle diagram of 16 Cygni B comparing optimal reference model with full set of frequencies (blue circles) from (Davies et al. 2015; black crosses) with the optimal model with a set restricted to frequencies with uncertainties lower than 1.5 μHz (red diamonds). 

In the text 
Fig. D.1. Échelle diagram of 16 Cygni A optimal model calculated with seismic indicators defined over frequencies which are not corrected for the surface effects. The crosses are the observed frequencies and the diamonds the theoretical ones. 

In the text 
Fig. D.2. Échelle diagram of 16 Cygni B optimal model calculated with seismic indicators defined over frequencies which are not corrected for the surface effects. The crosses are the observed frequencies and the diamonds the theoretical ones. 

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.