© ESO 2022

1. Introduction

The advent of space-based photometry missions such as CoRoT (Baglin et al. 2009), Kepler (Borucki et al. 2010; Chaplin et al. 2015) and TESS (Ricker et al. 2015) have made asteroseismic constraints standard tools for studying the internal structure and rotation of distant stars. Asteroseismology is also the golden path to determining precise and accurate stellar parameters, which is now an important goal and a necessary condition for the characterisation of exoplanetary systems. It even constitutes key requirements for the preparation of the PLATO mission (Rauer et al. 2014). Another major consequence of the availability of high-quality data was the generalisation of seismic analyses techniques until then only used in helioseismology, where they were very successful (see e.g. Christensen-Dalsgaard 2002, 2021; Basu & Antia 2008; Kosovichev 2011; Buldgen et al. 2019a and references therein for reviews). In this context, the 16Cygni binary system constitutes a prime target for their application because very high-quality seismic, spectroscopic, interferometric and astrometric data are available. This system has been extensively studied in the past years (see e.g. Mathur et al. 2012; Gruberbauer et al. 2013; Verma et al. 2014; Metcalfe et al. 2015; Bazot 2020), and is now considered a benchmark system for stellar modellers.

The generalisation of linear inversion techniques has already been foreseen and tested on synthetic data (see e.g. Gough & Kosovichev 1993b; Roxburgh et al. 1998; Thompson & Christensen-Dalsgaard 2002; Takata & Montgomery 2002; Basu 2003). An application to real data can be found a few years later in di Mauro et al. (2004) on Procyon A and more recently in Kosovichev & Kitiashvili (2020). Despite the successes of the space missions, the data quality in most cases remained far below that of helioseismic observations and the absence of high-degree modes, as a result of geometric cancellation, limited the potential of linear seismic inversions. However, applications of non-linear techniques can be found in Appourchaux et al. (2015), who attempted to carry out a full inversion of the internal structure of HIP 93511.

As a full scan of the internal structure of a distant star would not be possible with linear variational relations, Reese et al. (2012) focused on developing inversions for global quantities, starting with the mean density. The main goal was to exploit at best low-degree modes, for which some degree of localisation can only be achieved in the core, while they sometimes still show strong dependences on envelope properties. Thus, Reese et al. (2012) and subsequent works on indicator inversions studied target functions that could more easily be fitted and still provided relevant physical constraints on the stellar structure. Buldgen et al. (2015a,b, 2018) defined additional quantities, denoted indicators, that could be determined for main-sequence stars observed by Kepler (Buldgen et al. 2017b). They applied their technique to a few targets, including the 16Cyg binary system (Buldgen et al. 2016a,b). In their latest study of 16Cyg, they concluded that a full re-modelling was required, especially in light of the publication of new opacity tables (Mondet et al. 2015; Colgan et al. 2016; Le Pennec et al. 2015). While both 16Cyg A&B constitute excellent targets for testing non-standard physical processes such as turbulence at the base of the convective envelope and accretion of planetary material (Deal et al. 2015), these processes should be studied in a systematic way on a large set of models. In this approach, seismic inversions may be very attractive, as they may indicate the limitations of stellar evolution models. They are therefore highly valuable to improving the theory of stellar structure and evolution, but also to improve the determination of stellar fundamental parameters such as mass, radius, and age derived in a model-dependent way from stellar evolution computations. In turn, these models lead to more accurate planetary parameters that are crucial when exoplanetary systems are characterised in detail.

A full re-analysis with various evolutionary models was provided by the extensive modelling work of Farnir et al. (2020, hereafter Paper I). The authors carried out a detailed study of 16Cyg A&B using the WhoSGlAd oscillation spectrum modelling technique. They defined relevant structural indicators (Farnir et al. 2019), and varied physical ingredients, hypotheses of the modelling (stars seen as independent or joint), and observational constraints. This work provides a suitable set of evolutionary models that can be further tested using seismic inversions of indicators and local corrections of the structure. We start in Sect. 2 by recalling the principles of linear seismic inversions of stellar structure using the variational equations. We then present the evolutionary models we used in our study in Sect. 3. Section 4 presents our inversion results for the mean density, two core-condition indicators, and an envelope indicator for the 16Cyg binary system, and Sect. 5 shows results of local corrections of the squared isothermal sound speed, denoted u and defined as $u=\frac{P}{\rho }$, with P the local pressure and ρ the local density. Finally, in Sect. 6, we discuss our results of the inversions of indicators and local corrections and compare them to results reported in the literature (Buldgen et al. 2016a,b; Bellinger et al. 2017).

2. Formalism of variational structural seismic inversions

The concept of helioseismic or asteroseismic inversions is quite broad. It encompasses all approaches of inferences from seismic data of the internal structure of the Sun or of a distant star. For example, the determination of a set of optimal parameters for an evolutionary stellar model from the adjustment of seismic data already constitutes an inference of what its internal structure is, given a set of hypotheses regarding its evolution and the physical processes that govern it. In addition to such inferences based on evolutionary models, other approaches, such as the computation of static models, have been used for SdB stars or white dwarfs (see e.g. Charpinet et al. 2008; Giammichele et al. 2018) and also constitute examples of seismic inversions.

In the context of solar-like oscillators, structural inversions have taken various forms. Iterative inversions based on the reconnection of partial wave solutions have been developed (e.g. Roxburgh 2002, 2010) and applied to a few Kepler targets (Appourchaux et al. 2015; Roxburgh 2016, 2015) as well as to synthetic data (Roxburgh & Vorontsov 2002a,b). Applications of the variational inversions were also investigated early on by Gough & Kosovichev (1993a,b), Roxburgh et al. (1998), Thompson & Christensen-Dalsgaard (2002), Basu et al. (2002), and Basu (2003). In this context, the words “inversion techniques” most often referred to the application of the variational formulae to determine corrections to the internal structure of a star, thus only to a subcategory of seismic inversion techniques.

The basic equations for the linear inversion problem derive from the variational principle of adiabatic stellar oscillations (Chandrasekhar 1964; Clement 1964; Lynden-Bell & Ostriker 1967). In the case of a non-rotating, non-magnetic, isolated star, the problem can be written as a linear integral relation derived by Dziembowski et al. (1990),

$$\begin{array}{c}\hfill \frac{\delta {\nu }^{n,\ell }}{{\nu }^{n,\ell }}={\displaystyle {\int }_0^R{K}_{s_1,s_2}^{n,\ell }\frac{\delta s_1}{s_1}\mathrm{d}r+{\int }_0^R{K}_{s_2,s_1}^{n,\ell }\frac{\delta s_2}{s_2}\mathrm{d}r,}\end{array}$$(1)

where ν^n, ℓ is the oscillation frequency of radial order n and degree ℓ, s₁ and s₂ are variables of the adiabatic oscillation equations such as the density, ρ, squared adiabatic sound speed, c², pressure, P, the first adiabatic exponent ${\mathrm{\Gamma }}_1=\frac{\partial ln\rho }{\partial lnP}{|}_S$, etc. The functions $K_{s_i,s_j}^{n,\mathcal{l}}$ are the so-called kernel functions, which are related to the variable s_i in the structural pair (s_i, s_j). These kernel functions depend on the eigenfunctions of the oscillation modes and on the structure of the reference model. The expression δ denotes a linear perturbation of a given quantity such as the frequency or a structural variable, defined as

$$\begin{array}{c}\hfill \frac{\delta y}{y}=\frac{y_{\mathrm{Obs}}-y_{\mathrm{Ref}}}{y_{\mathrm{Ref}}},\end{array}$$(2)

where y is a given quantity. The suffixes Obs and Ref refer to the observed quantity and that of the reference model respectively. The reference model should be already close enough of the actual target so that the linear approximation defined in Eq. (1) is valid.

These relations can take various forms. The original expressions presented in Dziembowski et al. (1990) were derived for the density ρ and the squared adiabatic sound speed, $c^2=\frac{{\mathrm{\Gamma }}_{1}P}{\rho }$, but variable changes can be applied (see amongst others, Elliott 1996; Kosovichev 1999; Buldgen et al. 2017a) to derive expressions for any quantity of the adiabatic oscillation equations. Assuming also a linear development of the equation of state, we may define kernels related to chemical composition or temperature (see e.g. Gough & Kosovichev 1988), which were used for example to derive kernels for the helium abundance, denoted Y. While this will lead to non-negligible biases in helioseismic inversions (see Basu & Christensen-Dalsgaard 1997; Basu et al. 2009 for a discussion), the larger uncertainties in the context of asteroseismology allow us to consider that these biases are negligible. Basu et al. (2009) found variations of about 1% for density inversions and Basu & Christensen-Dalsgaard (1997) found variations of less than 0.1% for u inversions. These equations can also be written for frequency combinations or frequency ratios, as in Otí Floranes et al. (2005), and we can also define kernel functions related to Lagrangian perturbations instead of Eulerian ones, as in Christensen-Dalsgaard & Thompson (1997), or even perturbations at fixed acoustic radius, as noted in Pijpers & Thompson (1994).

Equation (1) suffers from a few caveats. First, it assumes a direct linear relation between the frequency differences and the structural corrections, which has only a limited validity that depends on the nature of the oscillation modes, as well as on the quality of the reference model, which will play a crucial role in asteroseismic inversions (Thompson & Christensen-Dalsgaard 2002). This impact justifies the construction of suitable sets of models to which the linear inversion is applied. In practice, this implies that the inference problem is separated into at least two steps, of which the linear inversion is the last. Second, a few hypotheses and simplifications regarding surface regions are made when deriving Eq. (1). The hypothesis of adiabaticity of the oscillations, the use of simplified boundary conditions for the adiabatic oscillation problem (neutrality condition on the Lagrangian pressure perturbation), the poor modelling of the upper convective layers by the mixing-length theory, as well as a few surface terms from integrations by parts are neglected in the derivations. Thus, the formalism does not apply in the upper layers of stars, and ad-hoc corrections for the so-called surface effects have to be included. These take the forms of polynomial corrections, for example following Rabello-Soares et al. (1999) or more recent empirical correction by Ball & Gizon (2014) or Sonoi et al. (2015). In practice, this implies that a third term needs to be added to the variational formula, written as

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{Surf}}={\displaystyle \sum _i}{a}_i{\nu }^i,\end{array}$$(3)

which can then be limited to take the form of the two terms of the Ball & Gizon (2014) correction, a linearised formulation of the Lorentzian correction of Sonoi et al. (2015) or the polynomial formulation applied in helioseismology, and where the coefficients a_i are related to the scaled inertias of the modes.

In addition to these shortcomings, Eq. (1) also suffers from the implicit definition of the boundary of the integral relation. In the Eulerian expression of the perturbations, the radius is assumed to be the same for the observed target and the reference model. If this is not the case, some caution has to be taken regarding the interpretation of the inferred values from the variational expressions. This issue was discussed in Buldgen et al. (2015b), who reviewed the validity of the variational formulations for kernels related to the squared isothermal sound speed $u=\frac{P}{\rho }$, as well as in Buldgen et al. (2017a).

Additional regularisation terms can be applied to the inversion to take care of the radii differences, as in Takata & Gough (2001) for the helioseismic case. In asteroseismology, the best scenario occurs when an independent radius determination is available from interferometry or from a combination of Gaia parallaxes and high-quality spectrocopy, although in the latter case, uncertainties regarding bolometric corrections and extinction have to be considered. A way to circumvent the issue is to add a term to the cost function of the inversion to scale the frequencies with respect to M/R³ (i.e. the dynamical time) and make them adimensional. This solution was suggested by Gough & Kosovichev (1993a) and also applied in Roxburgh et al. (1998).

In the case of 16Cyg, the availibility of high-quality spectroscopic, seismic, astrometric, and interferometric data makes the two stars prime targets for variational seismic inversions. This context indeed motivated previous analyses of their internal structure with variational techniques, such as those of Buldgen et al. (2016a,b), and Bellinger et al. (2017).

To summarise, the goal of the seismic inversion techniques is to determine meaningful corrections to stellar evolutionary models by solving Eq. (1) using dedicated numerical techniques. Hereafter, we use the substractive optimally localised averages technique (SOLA, Pijpers & Thompson 1994), implemented in the InversionKit software (Reese & Zharkov 2016).

3. Reference evolutionary models

We used the full set of models of Paper I, computed with the Liège stellar evolution code (CLES; Scuflaire et al. 2008a). The adiabatic oscillations were computed using the Liège oscillation code (LOSC; Scuflaire et al. 2008b). We show in Fig. 1 the different positions in the Hertzsprung-Russel (HR) diagram of the set of reference models, and in Fig. 2 we show their positions in a mass−age diagram.

Fig. 1. HR diagram showing the position of the reference models of Paper I with respect to both 16CygA&B. The green and red dots indicate the positions of 16CygA and 16CygB with their respective 1σ error bars.

Fig. 2. Age values of the various models determined from seismic evolutionary modelling in Paper I as a function of their mass.

All models were computed using the method of Farnir et al. (2019). We varied their physical ingredients, including or excluding classical constraints in the modelling and considered independent or joint modelling, where in the latter case, the same age and chemical composition is imposed to both components of the system. We count about 30 models for each star. We refer to Paper I for the details of the evolutionary modelling procedure and the detailed results and discussion associated with the considered variations in physical ingredients and constraints.

We briefly recall the variety of physical ingredients that we used to compute these models in Table 1. Multiple modifications were sometimes considered simultaneously, as mentioned in Paper I. Whenever turbulent diffusion was included, the following parametric mixing from Proffitt & Michaud (1991) was used

$$\begin{array}{c}\hfill D_{\mathrm{turb}}=D_T{\left(\frac{{\rho }_{\mathrm{bcz}}}{\rho (r)}\right)}^n,\end{array}$$(4)

with values of D_T = 1000, 2500, 10 000, n being fixed to 3, ρ_bcz denoting the density at the base of the convective envelope of the model, and ρ(r) being the local value of the density at a radial coordinate r in the model. For some models, microscopic and turbulent diffusion were entirely neglected. Envelope and core overshooting, when included, were treated as leading to an instantaneous mixing, with the radiative and adiabatic temperature gradient enforced, respectively, in the overshooting regions, following recommendations of Viallet et al. (2015). The value of the efficiency was fixed to 0.1H_P in both cases, with $H_{\mathrm{P}}=\frac{-\mathrm{d}r}{\mathrm{d}lnP}$ the local pressure scale height.

This set of stellar evolutionary models includes various physical ingredients and provides an extensive exploration of the possible solutions of detailed seismic modelling for 16Cyg A&B. Moreover, thanks to the optimisation procedure, we have an ensemble of models that is suitable for linear seismic inversions, built with a technique that takes the smooth component of the spectrum and the helium glitch into account simultaneously. As stated in Paper I, the convective envelope glitch could not be constrained in this case as its signal was too small with respect to the uncertainties. The main goal is to determine whether a subset of models stands out from the analysis of linear variational inversions that break free from the evolutionary framework, thus providing constraints on the physical processes acting inside the two stars.

The dataset we used to carry out the linear inversions is the same as we used for the evolutionary modelling in Paper I. It is the oscillation spectrum determined from the full duration of the Kepler mission for these stars presented by Davies et al. (2015). Non-seismic data were taken as in Paper I from Ramírez et al. (2009) for [Fe/H] and White et al. (2013) for T_eff and R.

4. Inversions of seismic indicators

The motivation for the inversion of seismic indicators has been briefly described above. It mainly results from the fact that the observations of only low-degree modes will limit the physical information that can be inferred from the frequencies. In other words, a full tomography of the internal structure of distant stars cannot be achieved, but a few meaningful physical quantities can be defined.

The SOLA inversion technique is well suited for this goal. Pijpers & Thompson (1994) already mentioned that it can be used to determine rotation gradients, for example, and that the target function can be adapted to more general forms than the classical Gaussian function used in helioseismic inversions. The original paper of Pijpers & Thompson (1994) started by giving a general form for the target function, the Gaussian target being a specific case aimed at improving results with respect to the original MOLA techniques regarding oscillatory wings of the averaging kernels.. The cost function of the SOLA method for indicator inversions is defined as

$$\begin{array}{cc}\hfill \mathcal{J}(c_i)=& {\displaystyle {\int }_0^R[K_{\mathrm{Avg}}(r)-\mathcal{T}(r)]\mathrm{d}r+\beta {\int }_0^R{K}_{\mathrm{Cross}}^2(r)\mathrm{d}r}\hfill \\ \hfill & +\lambda {\displaystyle [n-\sum _i{c}_i]+tan\theta \frac{{\sum }_i{(c_i{\sigma }_i)}^2}{⟨{\sigma }^2⟩}+\sum _k{a}_k\sum _i{c}_i{\psi }_k({\nu }_i),}\hfill \end{array}$$(5)

with 𝒯 the target function of the inversion, λ a Lagrange multiplier, c_i the inversion coefficients, n an integer linked to the indicator definition, θ and β the trade-off parameters, σ_i the uncertainties of the individual frequency differences, $⟨{\sigma }^2⟩=\frac{1}{N}\sum _{i=1}^N{\sigma }_i^2$, and ∑_ka_kψ_k(ν_i) is the polynomial expression of the surface correction (e.g. the correction of Ball & Gizon 2014). In addition to these quantities, we define in Eq. (5) two terms, K_Avg and K_Cross, the averaging and cross-term kernels. They are defined from the recombination of the structural kernels of Eq. (1) using the inversion coefficients.

In the most general case, the target function 𝒯 of an indicator A, related to the structural variable s₁, is defined using an integral relation of the form

$$\begin{array}{c}\hfill \frac{\delta A}{A}={\displaystyle {\int }_0^R\mathcal{T}(r)\frac{\delta s_1}{s_1}\mathrm{d}r.}\end{array}$$(6)

For this indicator, the averaging and cross-term kernels are defined from the inversion coefficients c_i as

$$\begin{array}{c}\hfill K_{\mathrm{Avg}}={\displaystyle \sum _i{c}_i{K}_{s_1,s_2}^i,}\end{array}$$(7)

$$\begin{array}{c}\hfill K_{\mathrm{Cross}}={\displaystyle \sum _i{c}_i{K}_{s_2,s_1}^i.}\end{array}$$(8)

The series of coefficients, c_i, that are determined from the minimisation of Eq. (5), allow us to estimate an inverted value for the indicator defined by Eq. (6) from the frequencies following

$$\begin{array}{c}\hfill {\left(\frac{\delta A}{A}\right)}_{\mathrm{Inv}}={\displaystyle \sum _i{c}_i\frac{\delta {\nu }_i}{{\nu }_i}.}\end{array}$$(9)

The accuracy of the inversion depends on the quality of the fit of the averaging kernel to the target function, on the damping of the contribution from the cross-term kernel, and on the accurate reproduction of the surface effects. Meanwhile, its precision is determined by the fourth term in Eq. (5), related to the propagation of the observational uncertainties of the individual frequencies. The determination of the trade-off parameters is made using a so-called L-curve analysis (see e.g. Backus & Gilbert 1970; Pijpers & Thompson 1994; Rabello-Soares et al. 1999).

The integer n is a form of normalisation, introduced for mean density inversions by Reese et al. (2012) and generalised in later works. It is an adaptation of the unimodularity constraint that is commonly used in helioseismic inversions, related to the link between the indicator A and the inverse of the dynamical time, going as ${\overline{\rho }}^{1/2}$, with $\overline{\rho }$ the mean density. In other words, n is defined from the relation $A\propto {\overline{\rho }}^{n/2}$ and can be used as an additional regularisation constraint on the inversion coefficients because it can be shown that their sum should be close to n (Reese et al. 2012).

In practice, each indicator will not be proportional to the dynamical time itself. This requires further discussion regarding the implicit scaling applied by the inversion technique, which has a direct impact on the verification of the integral relations between frequency and structure. The issue is for example seen when the inverted $\frac{\delta u}{u}$ values from synthetic data are compared to the actual differences between two models. If the radii of the two models are different, the relative differences of any quantity that do not scale with the mean density are automatically rescaled by the implicit hypothesis made on the radius when the variational equations are used. It can be shown that they behave as if the “observed” model had a mass defined as ${\stackrel{\sim }{M}}_{\mathrm{Obs}}=\frac{4\pi {\overline{\rho }}_{\mathrm{Obs}}{R}_{\mathrm{ref}}^3}{3}$, where the suffix “Obs” denotes the target model and “Ref” denotes the reference model.

The issue was already noted in early works, where it was argued that the seismic mean density and the spectroscopic surface gravity values would provide the ratio of M/R, which would solve the problem. Other works (Gough & Kosovichev 1993a; Roxburgh et al. 1998) also included a term in the SOLA cost function that is directly related to the minimisation of the mean density. While this certainly is an important effect, it does not disqualify the use of asteroseismic inversions, but shows that comparisons of inversion results from models with different radii should take this scaling into account.

As the inversion problem is defined as a trade-off between precision and accuracy, the evaluation of the optimal set of trade-off parameters from the sole use of the L-curve can sometimes be difficult. This is especially true for asteroseismic inversions, where the quality of the averaging kernels is not necessarily always good and the applicability of the variational equations is uncertain. In some cases, high accuracy and precision can be accidentally achieved through compensations, meaning that the inversion appears artificially robust. Consequently, when asteroseismic inversions are tested on artificial data, it is useful to compute the actual differences between the known solution and the inverted one using the following expressions:

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{Avg}}& ={\displaystyle {\int }_0^R(K_{\mathrm{Avg}}-{\mathcal{T}}_A)\frac{\delta s_1}{s_1}\mathrm{d}r,}\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{Cross}}& ={\displaystyle {\int }_0^R\left(K_{\mathrm{Cross}}\right)\frac{\delta s_2}{s_2}\mathrm{d}r,}\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{Res}}& =\frac{A_{\mathrm{Inv}}-A_{\mathrm{Obs}}}{A_{\mathrm{Ref}}}-{ϵ}_{\mathrm{Avg}}-{ϵ}_{\mathrm{Cross}},\hfill \end{array}$$(12)

with A_Inv the inverted value of the indicator, A_Obs the actual value of the indicator for the artificial target, and A_Ref the value of the indicator for the reference model of the inversion. For an artificial target, these quantities are trivial to compute because s₁ and s₂ are known. These expressions can be used to detect potential non-linearities, issues with surface effect corrections, or a compensation of errors within the computation of the integrals that would lead to an artificially high accuracy.

4.1. Mean density

Mean density inversions were defined in Reese et al. (2012) and then further studied by Buldgen et al. (2015a, 2019b) and were applied to various observed targets. They are applicable to a large number of stars due to their easily fitted target function.

The definition of the target function is given by

$$\begin{array}{c}\hfill {\mathcal{T}}_{\overline{\rho }}=\frac{4\pi \rho x^2}{{\rho }_R},\end{array}$$(13)

with $x=\frac{r}{R}$, ${\rho }_R=\frac{M}{R^3}$, and n is equal to 2 in Eq. (5).

The inversion results for both 16Cyg A&B are shown in Fig. 3 as a function of the mass of the reference models. The verifications of robustness we carried out are presented in Appendix A.1. We used the models of Paper I and computed the inversions using either the two-terms of Ball & Gizon (2014), the Sonoi et al. (2015) surface correction prescriptions, or no correction at all. The reference models are already in excellent agreement with the inversion, especially for 16CygA, where some of them are within the range provided by the SOLA inversion. In the case of 16CygB, a trend is observed where the reference models always show a slightly higher mean density than the inversion. The disagreement is, significantly reduced, however, if surface corrections are included. Because these corrections themselves have their own uncertainties and the variations between the reference values and the inverted ones are about 0.2%, we can consider that most of the evolutionary models already agree very well with the inversions.

Fig. 3. Inverted mean density as a function of mass for the set of reference models for 16Cyg A (left panel) and 16Cyg B (right panel). The red crosses are the reference values of the mean density, the orange, green, and blue crosses are the inversion results without any surface corrections, with the surface corrections of Ball & Gizon (2014) and Sonoi et al. (2015) respectively.

Overall, the precision of the inversion is very high. It reaches 0.3% even taking the dispersion due to model-dependency and surface effects into account. We can also note that the values remain consistent even though they span a relatively wide mass range between 0.1 and 0.2 M_⊙ for the two stars. The final values for the inverted mean densities of 16CygA&B determined from our study are ${\overline{\rho }}_{\mathrm{Inv},A}=0.827\pm 0.003$ g cm⁻³ and ${\overline{\rho }}_{\mathrm{Inv},B}=1.055\pm 0.003$ g cm⁻³. The already good agreement of the reference models with these inverted values gives confidence in the WhoSGlAd oscillation spectrum modelling technique that we used in Paper I to determine reliable fundamental stellar parameters.

4.2. Core condition indicators

Two core condition indicators were developed in Buldgen et al. (2015b, 2018), denoted t_u and S_Core. These inversions are more demanding in terms of the quality of the seismic data and are therefore only applicable to a more limited number of targets. They are also more prone to exhibit non-linear behaviours. The confirmations for these inversions are shown in Appendix A.2.

The target function for the t_u indicator is given by

$$\begin{array}{c}\hfill {\mathcal{T}}_{t_u}=-\frac{2u}{t_u}\frac{\mathrm{d}}{\mathrm{d}r}\left(f(r)\frac{\mathrm{d}u}{\mathrm{d}r}\right),\end{array}$$(14)

with f(r) defined by

$$\begin{array}{c}\hfill f(r)=r(r-R)exp(-7\frac{r^2}{R}),\end{array}$$(15)

with r the radial coordinate, R the stellar radius and $u=\frac{P}{\rho }$, the squared isothermal sound speed. In addition, t_u is defined by the integral

$$\begin{array}{c}\hfill t_u={\displaystyle {\int }_0^{R}f(r){\left(\frac{\mathrm{d}u}{\mathrm{d}r}\right)}^2\mathrm{d}r.}\end{array}$$(16)

From these definitions, Buldgen et al. (2015b) showed that n is equal to 4 in Eq. (5).

The target function for the S_Core indicator is given by

$$\begin{array}{c}\hfill {\mathcal{T}}_{S_{\mathrm{Core}}}=\frac{-g(r)}{S_{\mathrm{Core}}{S}_{5/3}},\end{array}$$(17)

with $S_{5/3}=\frac{P}{{\rho }^{5/3}}$ a proxy for the entropy of the stellar plasma. The weight function g(r) is defined by

$$\begin{array}{cc}\hfill g(r)=& r({\alpha }_{1}exp(-{\alpha }_2{(\frac{r}{R}-{\alpha }_3)}^2)+{\alpha }_{4}exp(-{\alpha }_5(\frac{r}{R}-{\alpha }_6)))\hfill \\ \hfill & \times tanh\left({\alpha }_7{(1-\frac{r}{R})}^4\right),\hfill \end{array}$$(18)

with α₁ = 16, α₂ = 26, α₃ = 0.06, α₄ = α₅ = 6.0, α₆ = 0.07, and α₇ = 50. These parameters may be varied depending on the star that is studied or the seismic dataset. The integral defition of the S_Core indicator is given by

$$\begin{array}{c}\hfill S_{\mathrm{Core}}={\displaystyle {\int }_0^R\frac{g(r)}{S_{5/3}}\mathrm{d}r.}\end{array}$$(19)

In this case, Buldgen et al. (2018) showed that n is equal to $\frac{-2}{3}$ in Eq. (5).

Both indicators are defined from the expected physical behaviours of core conditions and the rather intricate definitions of the target functions stem from the difficulties of accomodating a limited number of observed frequencies (compared to helioseismology) and the degeneracies linked to the availability of only low ℓ modes, often with larger uncertainties. Consequently, not all target functions are equivalent for various sets of observed modes and stars. The differences in mass and evolutionary stages directly affect the behaviour of the structural kernels and the target they can fit.

In Fig. 4, we present the results for the t_u indicator for 16CygA and 16CygB as a function of the mean density. In each case, we carried out the inversion without and with surface corrections, using the Sonoi et al. (2015) empirical formula in log g and T_eff. This approach allowed us to measure the impact of surface effects while keeping the same quality of fit of the target function by the averaging kernel, as adding two additional parameters to the cost function of the inversion reduces its accuracy, without affecting the conclusions of our study. We started with the case of 16CygB, for which the analysis is straightforward. All models constructed by the WhoSGlAd method agree with the inversion results, with or without surface corrections. A narrow domain of mean density is defined by the mean density inversion, as already noted in Fig. 3, but no significant correction in terms of t_u is observed. This indicates that the models built from the WhoSGlAd technique agree much better than the models of Buldgen et al. (2016a), which were built by fitting the individual small frequency separations and the inverted mean density or the inverted acoustic radius.

Fig. 4. Inverted core condition indicator tu as a function of the inverted mean density for the set of reference models for 16Cyg A (left panel) and 16Cyg B (right panel).

For 16CygA, the results are in line with those of Buldgen et al. (2016a) regarding the t_u values, with the inversion favouring a higher range of inverted indicator values. However, as already noted by these authors, the uncertainties of the inversion imply that the indicator cannot be directly used in a modelling procedure. The disagreement appears less important than the one found in Buldgen et al. (2015b). We discuss its potential origin below, but none of the models of Paper I appears to be in the higher range of values of the inversion results, despite the wide variety of tested physical ingredients. We also note that slightly increasing the trade-off parameter θ to reduce the error amplification leads to a smaller correction than is illustrated in Fig. 4, but still shows an increase in the t_u value for 16CygA. Similarly, including a surface correction also leads to slightly lower t_u values.

In Fig. 5, we illustrate the results for the two stars determined using the S_Core indicator. The conclusions are quite similar to those found for the t_u inversion. For 16CygA, a slight disagreement is found, with the inversion favouring models with a slightly lower mass or a larger radius, whereas for 16CygB, the models built from the evolutionary modelling procedure agree excellently with the inverted results. The slight disagreement for 16CygA almost completely disappears when a surface correction is applied to the frequencies. This implies that the corrections we have seen for the t_u inversion might not be the result of a disagreement in the internal structure of the models of 16CygA.

Fig. 5. Inverted core condition indicator SCore as a function of the inverted mean density for the reference models for 16Cyg A (left panel) and 16Cyg B (right panel).

These results emphasize the importance of the quality of the evolutionary modelling procedure to carry out the inversion. For 16CygB, all models are validated by the inversion, but not all models agree perfectly with the classical constraints or the interferometric radii. This contrasts with the results of Buldgen et al. (2016a), where by construction, the weighting between the classical and asteroseismic constraints in the cost function was made so that the seismic constraints did not dominate completely. In this case, the WhoSGlAd method allows us to reproduce very well the seismic data for 16CygB, but as already noted in Paper I, simultaneously reproducing all classical constraints would require additional free parameters. For 16CygA, the disagreement might be due to surface effects. Nevertheless, it seems that it remains overall less important than the disagreement found in Buldgen et al. (2016a), most likely as a result of the more efficient evolutionary modelling technique. We further discuss in Sect. 6 the implication of these results in terms of modifications of the physical inputs of the models.

4.3. Envelope indicator

The last seismic indicator we used is the one defined as a convective envelope indicator in Buldgen et al. (2018), denoted S_Env. This indicator aims at estimating the height of the entropy proxy plateau in the convective envelope, which indicates potential limitations of the temperature and mean molecular weight gradients in the upper radiative layers of a solar-type star. The target function of this indicator is defined as

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathrm{Env}}=\frac{h(r)S_{5/3}}{S_{\mathrm{Env}}},\end{array}$$(20)

with the weight function h(r) defined as

$$\begin{array}{cc}\hfill h(r)& =[{\alpha }_{1}exp(-{\alpha }_2{(\frac{r}{R}-{\alpha }_3)}^2)+{\alpha }_{4}exp(-{\alpha }_5{(\frac{r}{R}-{\alpha }_6)}^2)\hfill \\ \hfill & +\frac{0.78}{1+(exp((\frac{R}{r}-\frac{1}{{\alpha }_7})/{\alpha }_8))}]r^{{\alpha }_9}tanh\left({\alpha }_{10}(1-{\left(\frac{r}{R}\right)}^4)\right),\hfill \end{array}$$(21)

with α₁ = 30, α₂ = 120, α₃ = 0.31, α₄ = 7.3, α₅ = 26, α₆ = 0.33, α₇ = 1.7 α₈ = 1.2, α₉ = 1.5, and α₁₀ = 50. The S_Env indicator is defined as

$$\begin{array}{c}\hfill S_{\mathrm{Env}}={\displaystyle {\int }_0^{R}h(r)S_{5/3}\mathrm{d}r.}\end{array}$$(22)

For this indicator, Buldgen et al. (2018) showed that n is equal to $\frac{2}{3}$ in Eq. (5). As for the S_Core indicator, the α_i parameters will vary depending on the target and the available oscillation modes. In practice, the α_i parameters will vary with the observed star to determine the best target function. These will be mostly α₁, α₂, and α₃ as they define the amplitude, position, and localisation of the highest peak of the function (see e.g. the left panel of Fig. A.11) that leads to the largest variations of the indicator.

In Fig. 6 we illustrate the inversion results for both 16CygA&B for the S_Env indicator, and the tests of robustness are shown in Appendix A.3. We carried out the inversion with and without surface corrections and applied the same approach as in Sect. 4.2. The results show no significant corrections for the evolutionary models. In hindsight, this is expected because the method of Farnir et al. (2019) focused on reproducing a series of seismic indexes and the helium ionisation glitch signature simultaneously. The variations in the properties of the convective envelope are actually small amongst the models and do not exceed 2% for the two stars. This is beyond the resolution of the inversion technique.

Fig. 6. Inverted envelope condition indicator SEnv as a function of the inverted mean density for a subset of reference models for 16Cyg A (left panel) and 16Cyg B (right panel).

5. Localised inversions using a Gaussian target

In addition to studying global indicators, we attempted to carry out localised inversions of the profile of the squared isothermal sound speed, u. The most suitable structural kernels for this inversion are those of the (u,Y) pair (Basu et al. 2002; Thompson & Christensen-Dalsgaard 2002; Basu 2003). With these kernels, the cross-term of the inversion will be naturally reduced by the low amplitude of the Γ₁ derivatives with respect to Y and contributes to almost nothing to the cost function, almost transforming Eq. (1) into a single-integral relation as for rotation inversions. In addition, we can argue that the fitting technique used in Paper I has also intrinsically reduced the cross-term contribution.

We followed the definition of Rabello-Soares et al. (1999), used for the solar case

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathrm{Gauss}}=Arexp\left({(\frac{r-r_0}{R\mathrm{\Delta }}+\frac{R\mathrm{\Delta }}{2r_0})}^2\right),\end{array}$$(23)

with A a normalisation constant, r₀ the position at which one wishes to localise the averaging kernel and $\mathrm{\Delta }=\frac{{\mathrm{\Delta }}_{A}c(r_0)}{c_A}$ relates to the width of the Gaussian, with Δ_A a free parameter of the SOLA inversion and c(r₀) and c_A the adiabatic sound speed at the inverted point and at a reference radius of 0.2R, respectively.

The unimodularity term in this case is the one found in usual helioseismic inversions

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^R{K}_{\mathrm{Avg}}\mathrm{d}r=1.}\end{array}$$(24)

As for the previous cases, we started by carrying out inversions in specific test cases using various evolutionary models that are illustrated in Appendix A.4.

For the inversion of the observed data, we chose five models for 16CygA&B in our set of reference models. We note that results remained very similar for other cases, and we only present five results to avoid redundancy. The inversions were carried out using the two-term surface correction of Ball & Gizon (2014). Because the stellar radius and mass are unknown, the frequencies need to be adimensionalised (Gough & Kosovichev 1993a). In other words, this implies that the target and reference models have mean density values as close to each other as possible so that we can consider that relative differences in dimensional and adimensional frequencies are equal, such that

$$\begin{array}{c}\hfill \frac{{\nu }_{n,\ell }^{\mathrm{Obs}}-{\nu }_{n,\ell }^{\mathrm{Ref}}}{{\nu }_{n,\ell }^{\mathrm{Ref}}}\simeq \frac{{\stackrel{\sim }{\nu }}_{n,\ell }^{\mathrm{Obs}}-{\stackrel{\sim }{\nu }}_{n,\ell }^{\mathrm{Ref}}}{{\stackrel{\sim }{\nu }}_{n,\ell }^{\mathrm{Ref}}},\end{array}$$(25)

with $\stackrel{\sim }{\nu }={(\frac{\mathit{GM}}{R^3})}^{-1/2}\nu$ the adimensional frequency. We ensured that the reference model had the same mean density as the target well within 1% by using mean density inversions (see Sect. 4.1).

In Fig. 7, we illustrate our inversion results for 16CygA&B, using the surface correction Ball & Gizon (2014) in the cost function of the inversion technique as for the test cases. We rejected the models without an appropriate mean density value for both cases because we saw in the test cases that they could lead to spurious inversion results.

Fig. 7. Inversion of the adimensional squared isothermal sound speed $(\stackrel{\sim }{u})$ at various points of the internal structure as a function of the normalised radius for 16Cyg A (left panel) and 16Cyg B (right panel) for five models, including the Ball & Gizon (2014) correction in the cost function. The models and their associated colour code are provided in Tables 2 and 3.

The properties of the models summarised in Tables 2 and 3 for 16CygA and 16CygB respectively. They were chosen to present a variety of physical ingredients in their chemical mixture, opacities, and transport of chemicals they included but also in the constraints and free parameters used in their modelling in Paper I. The constraint “Same age” implies that the modelling has been carried out simultaneously for 16CygA&B imposing that the two stars have the same age.

Localised kernels could be obtained up to 0.25 stellar radii. Above this limit, we consider that the oscillatory wings (see Fig. A.14) of the kernels make the inversion unreliable for our set of trade-off parameters. Bellinger et al. (2017) showed kernels with a slightly less oscillatory behaviour for their kernels and trust their results up to 0.3 stellar radii. We were able achieve this with our software using different trade-off parameters, namely a larger width of the target function that allowed us to have kernels with an almost Gaussian form at a few additional locations. We also note that they do not show the amplitude of their kernels above the photosphere, where these have high amplitudes that change the results (see e.g. Fig. 2 of Basu 2003) and also strongly influence the choice of the trade-off parameters. Nevertheless, in our case, adding a few points between 0.25 and 0.3 stellar radii to the inversion does not affect our conclusions.

Figure 7 shows that the results for 16CygB indicate that there is no additional information to be extracted from the inversion. Variations of the trade-off parameters lead to similar conclusions, as the inverted values are essentially fully consistent with 0.

The case of 16CygA is slightly more interesting. The inversion seems to pinpoint slight differences at 0.05 and 0.23 stellar radii. It might even be argued that the points above 0.15 stellar radii seem to indicate a different slope in the u profile. Nevertheless, the differences are basically in agreement with 0 at 1.5σ for all models at all inverted points. They are also entirely suppressed for some values of the trade-off parameters.

The model that pecifically fits the interferometric radius and the effective temperature shows the best overall agreement with the inversion. In particular, it excellently agrees in effective temperature and luminosity, but has a higher metallicity by 2σ than the value that was determined from spectroscopy by Ramírez et al. (2009) and Tucci Maia et al. (2014), but it agrees slightly better with the values determined by Morel et al. (2021). From the inversion alone, it seems difficult to advocate that the models need strong corrections of their internal structure. A slight variation in radius of 0.006 R_⊙ between the light blue model results and those of the orange model changes from being consistent with 0 at 1σ regarding the inversion results. It thus appears that only slight modifications would likely bring most of the models in full agreement with the inversion, as we already devised from the indicator inversions. From a closer analysis, it also appears that non-standard models including additional transport at the base of the convective zone do not fare much better regarding the core inversions of isothermal sound speed. The orange model is indeed essentially a standard model while the red model includes an additional transport of chemicals in the form of turbulent diffusion. There is no clear argument to be made in favour of one or the other, implying that in this specific case, the SOLA inversion of the localised relative squared isothermal sound speed difference has reached its limits.

6. Discussion

The combination of various inversions, clearly shows that the models derived for the two stars do not show significant discrepancies in their inversion results. This can be expected for most models based on the dispersion of the results in terms of mass and age in Fig. 2 and on the results of Bazot (2020). Most of the models lie in a relatively narrow range, the tests of robustness of the inversion technique show that their internal structure in terms of indicators and adimensional isothermal sound speed are very similar; they differ by about 1% or less at most. This is well below the resolution limit of the inversion techniques and implies that if the models disagree with an inversion result for 16CygA&B, they should all do so in similar fashion.

To some extent, this is shown in the left panels of Figs. 4, 5 and 7. However, we should keep in mind that slight mismatches in the averaging kernels can influence the results and that the accuracy of the inversion itself is not perfect. Moreover, adding surface corrections causes the discrepancies to disappear completely for the S_Core indicator and reduces them slightly for the t_u indicator, implying that the cause of the differences might well not be the internal structure of the models but, the treatment of the surface effects in the cost function of the SOLA method. Furthermore, potential correlations between the determined individual frequencies (recently reported by O. Benomar for the specific case of 16CygA), could also lead to slightly larger uncertainties on the inverted results, bringing them in agreement with the evolutionary models used here. Slight departures from the solar-scaled chemical mixture could also have an impact on the agreement of the models with the inversions, without the need for invoking non-standard processes.

The fact that both standard and non-standard models fare equally well with respect to the inversions implies that to some extent, the internal structure of the star and its mean density are reproduced within a similar degree of agreement with various physical ingredients, as expected, and that the effects of the change in physics are compensated for by a change in chemical composition, mass and age. This result means that a certain degree of degeneracy exists in asteroseismic constraints. In other words, additional observational constraints on chemical abundances, interferometric radii, parallaxes, and so on will help lift some degeneracies and help provide more reliable models. For example, the disagreement of some of our models with spectroscopic constraints, already seen in Paper I, demonstrates some degree of divergence between the various observations, as illustrated in Fig. 8. Only a few models agree with the constraints for the surface chemical composition, especially for the case of 16CygB. This was already found in paper I and may be our best indicator, alongside lithium and beryllium depletion, that macroscopic mixing is acting inside these stars or that opacity modifications could be required. In this sense, both stars, being solar twins, are directly influenced by the revision of key ingredients of solar models.

Fig. 8. Surface helium abundance as a function of the surface metallicity for our reference models. The vertical lines indicate the limits derived from spectroscopic estimates by Ramírez et al. (2009). The helium abundance range is that given from the fit using the WhoSGlAd method.

Therefore, high-quality spectroscopic data cannot be neglected and the results of seismic inversions do not allow us to select a subset of models that agree significantly better with these non-seismic constraints. Most likely, the most efficient solution is to include them in the evolutionary modelling procedure, as was done for some models of Paper I.

Our models do not confirm the conclusions of Bellinger et al. (2017), who claimed that missing physical processes, improper applications of known processes or improper inputs in the model computations could be at the origin of their observed differences. We note that in another paper, Bellinger et al. (2019) reported that the uncertainties in their studies of 16Cygni are too large to conclude on differences with the predictions of stellar evolutionary models, as we find here. However, Bellinger et al. (2021) mentioned that the internal sound speed in the core of 16Cygni must be corrected and compared it to the case of HR7322, for which they invoked radiative accelerations. However, Buldgen et al. (2016a) were able to recover models with high inverted values of their t_u indicator by recomputing models that included more efficient microscopic diffusion, leading to significantly lower masses and younger ages for both stars. The results of Buldgen et al. (2016a) are not confirmed here either, as using a larger set of models built with a more robust seismic modelling technique leads to a much less obvious need for corrections in t_u and S_Core. While the trend remains, its origin is not entirely clear and might even potentially be spurious. In this respect we consider the results of the evolutionary modelling of Farnir et al. (2020) to be far more robust than the claimed accepted parameters of Buldgen et al. (2016a) for the 16Cyg system.

Regarding the comparison made in Bellinger et al. (2021) between HR7322 and 16Cyg, we can safely say that radiative accelerations can be directly ruled out, as they are not expected to drastically change the internal structure of a solar twin such as 16CygA (Deal et al. 2015). The inclusion of macroscopic mixing at the base of the convective zone is a promising candidate, however, not to explain the inversion results, as seen here, but rather to reproduce the lithium and beryllium depletions in both stars (King et al. 1997; Deliyannis et al. 2000). In this respect, the availability of rotation inversions from Bazot et al. (2019) allows us to test prescriptions of transport related to the effects of shear-induced turbulence and magnetic instabilities, as well as to place both stars in the context of the lithium and beryllium destruction in the Sun. The analysis of the impact of this transport on the intriguing discrepancies in light element abundances between 16CygA and 16CygB as well as the potential impact of accretion of planetary matter needs to be placed in perspective with formation, migration and engulfment scenarios for these stars. This requires a multidisciplinary approach as well as a robust seismic modelling procedure, here provided by the WhoSGlAd technique.

7. Conclusion

We studied the 16Cyg binary system using variational inversion techniques and the stellar evolutionary models from Farnir et al. (2020), to determine whether non-standard evolutionary models including macroscopic transport, or models including more recent opacity tables, can be differentiated from standard models by means of linear seismic inversions alone.

Our findings are similar to the results of Buldgen et al. (2016a) for 16CygA. Inversions of the core condition indicator, t_u, seem to indicate the need for a slightly lower mass and larger radius for 16CygA. The need for corrections is not as clear however, especially for inversions of the S_Core indicator and internal profile inversions, wich appear to agree with the internal structure of the models we used. Consequently, it seems difficult to advocate for a particular non-standard process that should be acting to reconcile the 16CygA models. Based on a larger set of models, we see that inversion results often agree with the reference models. When showing corrections, these are barely beyond the 1σ uncertainties. This implies that model dependency, uncertainty on the surface correction or on the individual frequencies could explain such small differences.

In the case of 16CygB, the models built using the WhoSGlAd technique agree excellently with all inversions, both standard and non-standard models including additional mixing.

Our results do not find the trends observed in Bellinger et al. (2017) and show the importance of a detailed and robust seismic modelling procedure before carrying out linear inversions of the stellar structure, coupled to models built including variations of their physical ingredients. In this particular case, the linear inversions validate the WhoSGlAd technique and show its efficiency and robustness for solar twins.

Nevertheless, the presence of non-standard mixing of chemical elements during the evolution of 16CygA&B can still be investigated using the lithium and beryllium abundances of the two stars, as in Deal et al. (2015). The hypothesis of the accretion of planetary matter needs to be placed in context with a detailed study of the two stars that takes the rotation inversion results of Bazot et al. (2019) and their implication for the depletion of light elements from macroscopic transport of chemicals and angular momentum into account. This will be the subject of a forthcoming paper in this series.

Acknowledgments

G.B and J.B. acknowledge fundings from the SNF AMBIZIONE grant No 185805 (Seismic inversions and modelling of transport processes in stars). M. F. is supported by the FRIA (Fond pour la Recherche en Industrie et Agriculture) – FNRS PhD Grant and acknowledges the support STFC consolidated grant ST/T000252/1. C.P. acknowledges fundings from the Swiss National Science Foundation (project Interacting Stars, number 200020-172505). 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.P. acknowledges support from F. R. S.-FNRS (Belgium) as a Chargé de Recherche. M.D. acknowledge the support of FCT/MCTES through the research grants UIDB/04434/2020, UIDP/04434/2020 and PTDC/FIS-AST/30389/2017, and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitividade e Internacionalização (grant: POCI-01-0145-FEDER-030389). M.D. is supported by national funds through FCT in the form of a work contract.

References

Appendix A: Inversion robustness tests

Because the linear approximation used to derive Eq. 1 has intrinsic limitations, the robustness and stability of the inversion must be tested thoroughly before it is applied to the actual observed data. For this study, we chose to carry out verifications for each inversion and for each star. We detail our approach in this section.

The reference and target models for the inversion were chosen from evolutionary models of Paper I. We considered that if the inversion technique is able to distinguish between various evolutionary models and the actual stars, it should be able to distinguish between evolutionary models built with various physical ingredients. Thus we chose the models by varying their reference solar abundances, opacity tables, including or excluding additional turbulence at the base of their convective envelope, as well as the inclusion or exclusion of non-seismic constraints in their modelling with the WhoSGlAd technique. The dataset for the inversion was chosen to be as similar as possible to that of the real targets, so that the trade-off problem for the verification is the same as for the inversion of the actual data. In other words, we used the exact same oscillation modes (in terms of n and ℓ) as were detected in 16CygA&B and assigned to each mode their actual observed uncertainties. We also provide the actual error contributions in Appendix. B for each indicator inversion verification presented below.

A.1. Mean density

The results of the mean density verifications are illustrated in Fig. A.1, and the averaging and cross-term kernels associated with some of these test cases are given in Figs. A.2 and A.3 together with the errors of the inversions in Tables B.1 and B.2.

Fig. A.1. Inversion results for the mean density, as a function of mass, using various stellar evolutionary models of 16CygA (left panel) and 16CygB (right panel) as artificial targets and the exact same dataset as the observed dataset in both cases. The results include error bars but they are almost invisible, and lead to an overestimated precision, especially in the case without surface corrections that we used here for the verification step.

Fig. A.2. Averaging kernel for the mean density inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the mean density inversion of 16CygA (right, blue) and its associated target function (green).

Fig. A.3. Averaging kernel for the mean density inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the mean density inversion of 16CygB (right, blue) and its associated target function (green).

The mean density inversion is robust and quite accurate. However, as noted in Reese et al. (2012) and Buldgen et al. (2015a), it overestimates its actual precision. The derivation of the uncertainties of inverted results in the SOLA technique is not always optimal. Thus, for mean density inversions, the actual precision of the method must consider some dispersion coming from model dependence and from the imperfect reproduction of the target function of the inversion. This emphasizes the need to carry out the inversion from multiple reference models using a wide variety of physical ingredients.

A.2. Core condition indicators

The results of the verifications are shown in Figs. A.4 and A.7. For both cases, the results are well within 1σ of the actual target values for both t_u and S_Core. The mean density values are taken from the mean density inversion carried out in Appendix A.1. The inversion also remains stable because no large spurious corrections are observed. However, given the small variations found between the models, we can already foresee that the inversion will be at its resolution limit to provide constraints on the internal structure of 16CygA&B. Illustrations of the averaging and cross-term kernels for specific test cases for the two stars are shown in Figs. A.5 and A.6.

Fig. A.4. Inversion results for the core condition indicator, tu, as a function of the mean density, using various stellar evolutionary models of 16CygA (left panel) and 16CygB (right panel) as artificial targets and the exact same dataset as the observed dataset in each case.

Fig. A.5. Averaging kernel for the tu inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the tu inversion of 16CygA (right, blue) and its associated target function (green).

Fig. A.6. Averaging kernel for the tu inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the tu inversion of 16CygB (right, blue) and its associated target function (green).

In the case of the S_Core inversion, the results are stable and can also lead to meaningful corrections, as illustrated by the crosses on the left in the two panels of Fig. A.7. Thus, the inversion appears to be accurate and stable. The averaging and cross-term kernels for some of these inversions are shown in Figs. A.8 and A.9.

Fig. A.7. Inversion results for the core condition indicator, SCore, as a function of the mean density, using various stellar evolutionary models of 16CygA (left panel) and B (right panel) as artificial targets and the exact same dataset as the observed dataset for both stars.

Fig. A.8. Averaging kernel for the SCore inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the SCore inversion of 16CygA (right, blue) and its associated target function (green).

Fig. A.9. Averaging kernel for the SCore inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the SCore inversion of 16CygB (right, blue) and its associated target function (green).

A.3. Envelope indicators

As shown in Fig. A.10, the inversion remains stable and does not provide corrections between the reference model and the target. This may again indicate that it might not be able to provide any significant corrections for this indicator, but confirms the reliability of the models. This inversion is particularly difficult because the properties of the convective envelope are difficult to extract using only low-degree modes. Nevertheless, a relatively good agreement between the averaging kernels and the target function can be achieved, as illustrated for both stars for a specific test case in Figs. A.11 and A.12.

Fig. A.10. Inversion results for the envelope condition indicator, SEnv, as a function of the mean density, using various stellar evolutionary models of 16CygA (left panel) and B (right panel) as artificial targets and the exact same dataset as the observed dataset in both cases.

Fig. A.11. Averaging kernel for the SEnv inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the SEnv inversion of 16CygA (right, blue) and its associated target function (green).

Fig. A.12. Averaging kernel for the SEnv inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the SEnv inversion of 16CygB (right, blue) and its associated target function (green).

A.4. Localised inversions

We illustrate the results of the robustness tests of localised inversions in Fig. A.13, where the left panel shows the inversions of 16CygA models and the right panel shows those for 16CygB models. The crosses are the corrections obtained using the SOLA technique, and the dashed and plain lines are the actual dimensional and adimensional differences between the two models in squared isothermal sound speed, defined following

$$\begin{array}{c}\hfill \frac{u^{\mathrm{Obs}}-u^{\mathrm{Ref}}}{u^{\mathrm{Ref}}}=\frac{(M_{\mathrm{Obs}}/R_{\mathrm{Obs}}){\stackrel{\sim }{u}}^{\mathrm{Obs}}-(M_{\mathrm{Ref}}/R_{\mathrm{Ref}}){\stackrel{\sim }{u}}^{\mathrm{Ref}}}{(M_{\mathrm{Ref}}/R_{\mathrm{Ref}}){\stackrel{\sim }{u}}^{\mathrm{Ref}}},\end{array}$$(A.1)

Fig. A.13. Inversion of relative differences of the squared isothermal sound speed as a function of the normalised radius using various evolutionary models as references and targets for both 16CygA (left panel) and 16CygB (right panel). The inverted differences are compared in terms of the dimensional, u, and adimensional squared adiabatic sound speed, $\stackrel{\sim }{u}$. The colours are associated with three different pairs of reference models and targets, all computed in Paper I. The continuous and dashed lines correspond to the actual differences between the target and reference models in adimensional and dimensional squared isothermal sound speed.

Fig. A.14. Averaging kernel for the localised inversion of adimensional squared isothermal sound speed of 16CygA for various positions (left). Same for 16CygB (right).

with M_Obs and R_Obs the mass and radius of the observed target and M_Ref and R_Ref, those of the reference model.

In both panels, green symbols are related to a target model with a mean density difference of about 2% with respect to the reference model. In these conditions neither the dimensional nor adimensional squared isothermal sound speed can be estimated reliably.

The blue and magenta symbols show results for models with mean density values well within 1% of each other. In the left panel, the results are very similar for the adimensional and the dimensional isothermal sound speed within 1σ. The actual differences between the models remain within 1%, except for the very deep core, that is not probed by pure p-modes. A similar situation is found in helioseismology, where the sound speed inversions usually stop at about 0.05R_⊙.

In the right panel, the magenta results are related to specific models with the same mean density, but very different M/R. In these conditions, no meaningful information is obtained on the dimensional isothermal sound speed. The inversion results also almost disagree with the actual differences of $\stackrel{\sim }{u}$ around 0.12R. The blue symbols are related to models with almost exactly the same mass and radius, thus u and $\stackrel{\sim }{u}$ are essentially the same. In this case, the inversion results agree excellently with the actual differences.

From our tests, we conclude that inversions of $\stackrel{\sim }{u}$ are only meaningful if the mean density is determined well below 1%. Similarly, the dimensional sound speed can only be determined in the presence of reliable mass and radius estimates (beyond the reaches of classical seismic modelling when model-independent inversions results are the goal). In these conditions, spectroscopic binary systems or targets with independent radius determination using interferometry are prime candidates.

However, these test cases also indicate a similar behaviour for most of the indicators, namely that the SOLA inversions will mostly validate the quality of the models built by Farnir et al. (2020).

Appendix B: Error contributions of the inversion robustness tests

We provide here some additional data regarding the verification inversions between the various models in Tables B.1 and B.2 including the errors as defined in Eqs. 10, 11, and 12.

All Tables

All Figures

	Fig. 1. HR diagram showing the position of the reference models of Paper I with respect to both 16CygA&B. The green and red dots indicate the positions of 16CygA and 16CygB with their respective 1σ error bars.
In the text

	Fig. 2. Age values of the various models determined from seismic evolutionary modelling in Paper I as a function of their mass.
In the text

	Fig. 3. Inverted mean density as a function of mass for the set of reference models for 16Cyg A (left panel) and 16Cyg B (right panel). The red crosses are the reference values of the mean density, the orange, green, and blue crosses are the inversion results without any surface corrections, with the surface corrections of Ball & Gizon (2014) and Sonoi et al. (2015) respectively.
In the text

	Fig. 4. Inverted core condition indicator tu as a function of the inverted mean density for the set of reference models for 16Cyg A (left panel) and 16Cyg B (right panel).
In the text

	Fig. 5. Inverted core condition indicator SCore as a function of the inverted mean density for the reference models for 16Cyg A (left panel) and 16Cyg B (right panel).
In the text

	Fig. 6. Inverted envelope condition indicator SEnv as a function of the inverted mean density for a subset of reference models for 16Cyg A (left panel) and 16Cyg B (right panel).
In the text

	Fig. 7. Inversion of the adimensional squared isothermal sound speed $(\stackrel{\sim }{u})$ at various points of the internal structure as a function of the normalised radius for 16Cyg A (left panel) and 16Cyg B (right panel) for five models, including the Ball & Gizon (2014) correction in the cost function. The models and their associated colour code are provided in Tables 2 and 3.
In the text

	Fig. 8. Surface helium abundance as a function of the surface metallicity for our reference models. The vertical lines indicate the limits derived from spectroscopic estimates by Ramírez et al. (2009). The helium abundance range is that given from the fit using the WhoSGlAd method.
In the text

Fig. A.1. Inversion results for the mean density, as a function of mass, using various stellar evolutionary models of 16CygA (left panel) and 16CygB (right panel) as artificial targets and the exact same dataset as the observed dataset in both cases. The results include error bars but they are almost invisible, and lead to an overestimated precision, especially in the case without surface corrections that we used here for the verification step.

In the text

	Fig. A.2. Averaging kernel for the mean density inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the mean density inversion of 16CygA (right, blue) and its associated target function (green).
In the text

	Fig. A.3. Averaging kernel for the mean density inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the mean density inversion of 16CygB (right, blue) and its associated target function (green).
In the text

	Fig. A.4. Inversion results for the core condition indicator, tu, as a function of the mean density, using various stellar evolutionary models of 16CygA (left panel) and 16CygB (right panel) as artificial targets and the exact same dataset as the observed dataset in each case.
In the text

	Fig. A.5. Averaging kernel for the tu inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the tu inversion of 16CygA (right, blue) and its associated target function (green).
In the text

	Fig. A.6. Averaging kernel for the tu inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the tu inversion of 16CygB (right, blue) and its associated target function (green).
In the text

	Fig. A.7. Inversion results for the core condition indicator, SCore, as a function of the mean density, using various stellar evolutionary models of 16CygA (left panel) and B (right panel) as artificial targets and the exact same dataset as the observed dataset for both stars.
In the text

	Fig. A.8. Averaging kernel for the SCore inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the SCore inversion of 16CygA (right, blue) and its associated target function (green).
In the text

	Fig. A.9. Averaging kernel for the SCore inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the SCore inversion of 16CygB (right, blue) and its associated target function (green).
In the text

	Fig. A.10. Inversion results for the envelope condition indicator, SEnv, as a function of the mean density, using various stellar evolutionary models of 16CygA (left panel) and B (right panel) as artificial targets and the exact same dataset as the observed dataset in both cases.
In the text

	Fig. A.11. Averaging kernel for the SEnv inversion of 16CygA (left, blue) and its associated target function (green). Cross-term kernel for the SEnv inversion of 16CygA (right, blue) and its associated target function (green).
In the text

	Fig. A.12. Averaging kernel for the SEnv inversion of 16CygB (left, blue) and its associated target function (green). Cross-term kernel for the SEnv inversion of 16CygB (right, blue) and its associated target function (green).
In the text

Fig. A.13. Inversion of relative differences of the squared isothermal sound speed as a function of the normalised radius using various evolutionary models as references and targets for both 16CygA (left panel) and 16CygB (right panel). The inverted differences are compared in terms of the dimensional, u, and adimensional squared adiabatic sound speed, $\stackrel{\sim }{u}$. The colours are associated with three different pairs of reference models and targets, all computed in Paper I. The continuous and dashed lines correspond to the actual differences between the target and reference models in adimensional and dimensional squared isothermal sound speed.

In the text

	Fig. A.14. Averaging kernel for the localised inversion of adimensional squared isothermal sound speed of 16CygA for various positions (left). Same for 16CygB (right).
In the text