# Analysis of the linear approximation of seismic inversions for various structural pairs

G. Buldgen, D. Reese, M. Dupret

► To cite this version:

G. Buldgen, D. Reese, M. Dupret. Analysis of the linear approximation of seismic inversions for various structural pairs. Astronomy & Astrophysics - A&A, 2017, 598, pp.A21. 10.1051/0004-6361/201629485. hal-02266197

HAL Id: hal-02266197

https://hal.science/hal-02266197v1

Submitted on 14 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Analysis of the linear approximation of seismic inversions for various structural pairs

G. Buldgen<fnref n="1" />, D. R. Reese<fnref n="2" />, and M. A. Dupret<fnref n="1" />

Received 5 August 2016 / Accepted 1 October 2016

## ABSTRACT

*Context.* Thanks to the space-based photometry missions CoRoT and *Kepler*, we now benefit from a wealth of seismic data for stars other than the sun. In the future, K2, Tess, and *Plato* will complement this data and provide observations in addition to those already at hand. The availability of this data leads to questions on how it is feasible to extend kernel-based, linear structural inversion techniques to stars other than the sun. Linked to the inversion problem is the question of the validity of the linear assumption. In this study, we analyse the limitations of this assumption with respect to changes of structural variables.

*Aims.* We wish to provide a more extended theoretical background to structural linear inversions by doing a study of the validity of the linear assumption for various structural variables. We thus point towards limitations in inversion techniques in the asteroseismic and helioseismic cases.

*Methods.* First, we recall the origins of the linear assumption for structural stellar inversions and explain its importance for asteroseismic studies. We also briefly recall the impact of unknown structural quantities such as the mass and the radius of the star on structural inversion results. We then explain how kernels for new structural variables can be derived using two methods, one suited to asteroseismic targets, the other to helioseismic targets. For this second method, we present a new structural pair, namely the (A, Y) structural kernels. The kernels are then tested in various numerical experiments that enable us to evaluate the weaknesses of different structural pairs and the domains of validity of their respective linear regime.

*Results.* The numerical tests we carry out allow us to disentangle the impact of various uncertainties in stellar models on the verification of the linear integral relations. We show the importance of metallicity, the impact of the equation of state, extra-mixing, and inaccuracies in the microphysics on the verification of these relations. We also discuss the limitations of the method of conjugated functions due to the lack of extremely precise determinations of masses or radii in the asteroseismic context.

**Key words.** stars: interiors – stars: oscillations – stars: fundamental parameters – asteroseismology

## 1. Introduction

Asteroseismology is now considered the golden path to the study of stellar structure. This young research field benefits from high quality data for a large sample of stars stemming from the successes of the CoRoT, *Kepler*, and K2 missions (Baglin et al. 2009; Borucki et al. 2010; Chaplin et al. 2015). More specifically, the detection of solar-like oscillations in a large sample of stars now allows a more accurate study of stellar structure. In the future, the Tess and *Plato* missions (Rauer et al. 2014) will carry on what is now called the space-photometry revolution.

Historically, the successes of asteroseismology were preceded by the successes of helioseismology, the study of solar pulsations. Indeed, the quality of seismic data of the sun is still far beyond what is achievable for other stars, even in the era of the space missions. In the particular field of helioseismology, numerous methods were developed to obtain constraints on the solar structure. Amongst them, inversion techniques lead to the successful determination of the position of the base of the convective envelope, the helium abundance in this region, and the rotational profile of the sun (Schou et al. 1998). The determination of the sound speed and density profiles also demonstrated the importance of microscopic diffusion for solar models (e.g. Basu & Antia 2008).

In the context of asteroseismology, the use of these inversion techniques can now be considered for a limited number of targets under the conditions of validity of all the hypotheses hiding behind the basic equations defining their applicability domain. The most constraining of these hypotheses is to assume a linear relation between frequency differences and structural differences. In this paper, we propose to analyse in depth the issues surrounding the linearity of these relations for various structural pairs and more specifically for kernels of the convective parameter. To that extent, we derive new structural kernels for the convective parameter and helium abundance and compare their linear behaviour to other, pre-existing kernels. Differences in the verification of the linearity of the frequency-structure relation mean that care should be taken when combining seismic diagnostics of various kernels, even in the solar case. Differences can stem from the intrinsic non-linear behaviour of the structural variables considered, but can also be the results of inaccuracies in terms of numerical quality of the models and/or kernels.

## 2. The variational principle and linear frequency-structure equations

The variational principle is a well-known property of adiabatic stellar oscillation equations. In fact, it can be extended to more

<footnote n="1">¹ Institut d'Astrophysique et Géophysique de l'Université de Liège, Allée du 6 Août 17, 4000 Liège, Belgium
e-mail: gbuldgen@student.ulg.ac.be</footnote>

<footnote n="2">² LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon Cedex, France</footnote>
general objects than stars and generalized beyond the classical case presented in helio- and asteroseismology. The history of the variational principle can be traced back to stability analysis in structural mechanics, but its application in seismology stems from the pioneering work of Chandrasekhar (Chandrasekhar 1964) and the generalization of his study by other authors the following years (Clement 1964; Lynden-Bell & Ostriker 1967).

Far beyond the historical interest of the discovery of this mathematical property, the hypotheses that lay behind this principle are still important since they are at the heart of intrinsic limitations of the frequency-structure relation. From the mathematical point of view, the variational principle is a consequence of the symmetry of the operator associated with adiabatic stellar pulsations. Mathematically, this means that given two functions $\xi$ and $\psi$ and $\mathcal{D}$ the operator associated with adiabatic oscillations, we have the property

$$ \langle \xi, \mathcal{D}(\psi) \rangle = \langle \psi, \mathcal{D}(\xi) \rangle, \qquad (1) $$

where $\langle, \rangle$ denotes the scalar product associated with the functional space defined as

$$ \langle \xi, \psi \rangle = \int_V \rho \xi \psi^* dV, \qquad (2) $$

with the symbol $^*$ denoting the complex conjugate. The absence of symmetry in the non-adiabatic case is the reason why all inversions are carried out using the hypothesis of adiabaticity of stellar oscillations.

Moreover, the hypothesis of linearity of the frequency-structure relations is a strong restriction to the validity of seismic inversions and in this section we briefly discuss how this hypothesis influences structural diagnostics from inversion techniques. The relation between perturbations of the frequencies and small perturbations of the stellar structure can be obtained by perturbing the variational expression of the pulsation frequencies. The symmetry of the stellar pulsation operator is then used to eliminate perturbations of the eigenfunctions in the resultant expression. Other effects such as perturbations of the eigenfunctions associated with each pulsation frequency can be neglected to first order. This implies that a direct relation can be obtained between structural differences and pulsation frequencies only. This relation is formally written:

$$ \delta v = \frac{\langle \xi, \delta \mathcal{D}(\xi) \rangle}{\langle \xi, \xi \rangle}, \qquad (3) $$

with $\delta v$ the perturbation of an oscillation frequency, $\delta \mathcal{D}$ the associated perturbation to the operator of adiabatic oscillations. The eigenfunctions $\xi$ and the unperturbed operator $\mathcal{D}$ are known and defined from the reference model.

In practice, Eq. (3) implies that small differences in frequencies can be used to analyse the associated differences in the operator of adiabatic pulsations. The main problem is the scalar product which implies integral relations and thus an ill-posed problem. However, it should be noted that the validity of the variational expression is limited, since we are speaking of small perturbations, the term “small” being misleading because it is often retroactively defined. In other words, perturbations are small because the variational expression is satisfied, but the quantification of how small a perturbation can be and if all variables can be similarly perturbed remains uncertain.

The classical equation for inversion techniques is the result of further developments introduced in the variational expression, assuming spherical symmetry of the star, and carrying out integration by parts and permutation of integrals. This

ends leading to the following formally simple equation (see Gough & Thompson 1991, for a full demonstration of this expression and its hypotheses):

$$ \frac{\delta v^{n,l}}{v^{n,l}} = \int_{0}^{R} K_{\rho,c^2}^{n,l} \frac{\delta \rho}{\rho} dr + \int_{0}^{R} K_{c^2,\rho}^{n,l} \frac{\delta c^2}{c^2} dr + O(\delta^2), \qquad (4) $$

with the following definition:

$$ \frac{\delta x}{x} = \frac{x_{\text{obs}} - x_{\text{ref}}}{x_{\text{ref}}}. \qquad (5) $$

The quantity $x$ can be the oscillation frequency of a particular mode, $v^{n,l}$, the density, $\rho$, the squared adiabatic sound speed, $c^2$, or other quantities for which kernels can be derived.

First of all, we note that to this expression is usually added the surface effects term, which is an empirical correction that is added to Eq. (4) to take into account the improper modelling of surface layers in the computation of oscillation frequencies of stellar models. In this study, we do not consider this surface term since we will only compare theoretical models and the validity of the linear approximation for various test cases between these models.

In terms of seismic diagnostic, the linear hypothesis puts strong restrictions on the applicability of inversion techniques. In fact, in some regions of the Hertzsprung-Russell diagram, it seems obvious that differences in frequencies may not be linearly related to structural changes. For example, in evolved stars, changes in the mixed modes frequencies will have strong impact on the coupling of the $p$ and $g$ mode cavities. Thus a small change in frequency will imply a strong change in the eigenfunctions. In this particular case, the second order terms neglected in the variational analysis may well become dominant and have to be modelled to efficiently use kernel based inversions as a seismic diagnostic.

For the case of $p$ modes observed in solar-like stars, one could say that provided the model is good enough, the linear approximation may be used. However, the linear approximation as presented is usually for a slow-rotating, non-magnetic, isolated star<fnref n="1" />. The problem is not to carry out the inversion, since that can be done provided a sufficient number of frequencies is available, but to decide whether the inverted results can be trusted. The errors due to linearity are intrinsically not seen by the inversion technique. However, it is still possible to witness their effects on inversion results and indirectly assess the quality of the reference model. To do so, one simply has to start from various reference models and analyse the variation of the results with the model. This simple and straightforward method is well-adapted to global optimization techniques which generate a large sample of models. However, this does not mean that using a large number of models, one can go beyond the linear approximation of the variational principle, it only implies that one can analyse the errors coming from the non-linear effects and decide whether the results should be trusted or not.

An other important aspect of asteroseismic inversions which has been reported by Basu (2003) and described in Buldgen et al. (2015) is that the inversion scales its results implicitly. This scaling stems from the assumption that integral relations are defined on the same domain for the reference model as for the observed target. In other words, if we define $R_{\text{Ref}}$, the radius of the reference model and $R_{\text{Tar}}$, the radius of the observed target, the inversion will wrongly consider that both radii are equal. However,

<footnote n="1">¹ By isolated, we mean that it is not in a close binary system where the gravitational influence of the neighbouring companion would change the geometry of the star and its oscillation modes.</footnote>
<table><caption>**Table 1.** Summary of the properties of the differential equations for various structural pairs with the direct method.</caption><thead><tr><th>Pair</th><th>Order of o.d.e.</th><th>Integration by parts</th></tr></thead><tbody><tr><td>ρ, Γ<sub>1</sub> (or Y)</td><td>0-algebraic</td><td>No</td></tr><tr><td>g, Γ<sub>1</sub> (or Y)</td><td>1</td><td>No</td></tr><tr><td>P, Γ<sub>1</sub> (or Y)</td><td>2</td><td>No</td></tr><tr><td>c<sup>2</sup>, Γ<sub>1</sub> (or Y)</td><td>2</td><td>No</td></tr><tr><td>u, Γ<sub>1</sub> (or Y)</td><td>2</td><td>No</td></tr><tr><td>A, Γ<sub>1</sub> (or Y)</td><td>3</td><td>Yes</td></tr><tr><td>N<sup>2</sup>, c<sup>2</sup></td><td>3</td><td>Yes</td></tr></tbody></table>

since the inversion uses seismic information, the mean density of the observed target is known. Consequently, the mass of the scaled target, denoted here $\tilde{M}_{\text{Tar}}$, which is studied by the inversion is $\tilde{M}_{\text{Tar}} = R_{\text{Ref}}^3 \tilde{\rho}_{\text{Tar}}$, with $\tilde{\rho}_{\text{Tar}}$ the mean density of the observed target. This means that structural quantities such as the squared sound speed, or indicators defined by integrated quantities, are not determined for the observed target itself, but for the scaled one and are related to the observed quantities through an homology. While this does not reduce the diagnostic potential of inversion techniques, it should of course be taken into account when comparing results inverted from various reference models.

### 3. Changing the structural pair

The calculation of new kernels is particularly interesting in the context of asteroseismology, where the change of structural variables can significantly improve the ability to fit a certain target while reducing the contribution from the so-called “cross-term”. Additional kernels have also been used in helioseismology to test the equation of state used in solar models and to try to determine the helium abundance in the convective envelope. In this section, we present two methods to derive additional structural kernels from Eq. (4) and discuss in more details their implementation and respective strengths and weaknesses.

#### 3.1. Direct method

We call this approach direct because it consists in a direct change of variables within Eq. (4) (or any similar relation), from which a linear differential equation is derived (This equation can be of the first, second or third order depending on the variables involved). The resolution of this equation allows us to determine new kernels, provided the proper boundary conditions are applied.

In practice, this method gives access to any function of $\rho$, $c^2$, $\Gamma_1$ or their integrals (e.g. hydrostatic pressure, P or the gravitational acceleration, g), or combinations of these variables (e.g. the squared isothermal sound speed, $u = \frac{P}{\rho}$)$^2$. However, it should be noted that this method does not give access to any function of the derivative of the density without further integration by parts when deriving the differential equation. The kernels that can be obtained through the direct method are listed in Table 1. We mention that this list contains only kernels for which the equations have been derived, but one could be interested to define new thermodynamical variables and to obtain kernels for these new variables.

This method has been partially presented in a previous paper (Buldgen et al. 2015) and referred to as Masters’ method, because it was developed as an extension of an approach presented in Masters (1979) for geophysical applications that was mentioned in Gough & Thompson (1991) as a potential method for obtaining kernels for the Brunt-Väisälä frequency$^3$. Originally, Masters’ approach proposed to solve directly the integral relations between structural kernels used in geophysics. In asteroseismology, the method could have been similar. First, we start with Eq. (4) and consider for example the change from the $(\rho, c^2)$ structural pair to the $(g, c^2)$ structural pair, where g is the gravitational acceleration and is written:

$$g = \frac{Gm(r)}{r^2}, \qquad (6)$$

with $m(r)$ the mass of stellar material contained in a sphere of radius $r$ and being defined:

$$m(r) = \int_{0}^{r} 4\pi \tilde{r}^{2} \rho d\tilde{r}. \qquad (7)$$

If we consider the linear relative perturbation of the gravity acceleration, we obtain:

$$\frac{\delta g}{g} = \frac{\frac{G \partial m}{r^2}}{\frac{G m}{r^2}} = \frac{\delta m}{m} = \frac{1}{m} \int_{0}^{r} 4\pi \tilde{r}^{2} \delta \rho d\tilde{r}. \qquad (8)$$

This expression can be used directly in the integral relations for the structural kernels. Indeed, if kernels of the pair $(g, c^2)$ can be found, they must satisfy the following relation:

$$\begin{aligned} \frac{\delta v^{n,l}}{\nu^{n,l}} &= \int_{0}^{R} K_{g,c^2}^{n,l} \frac{\delta g}{g} dr + \int_{0}^{R} K_{c^2,g}^{n,l} \frac{\delta c^2}{c^2} dr \\ &= \int_{0}^{R} K_{\rho,c^2}^{n,l} \frac{\delta \rho}{\rho} dr + \int_{0}^{R} K_{c^2,\rho}^{n,l} \frac{\delta c^2}{c^2} dr. \end{aligned} \qquad (9)$$

From the second equality, we have the integral relation that we searched. One has only to introduce the perturbation of the gravitational acceleration and permute the integrals such that the perturbation of density is in the outermost integral. From there it is easy to obtain a simple relation between kernels. Indeed, Eq. (9) must be satisfied for any perturbation within the linear regime, since the kernels must be dependent on the reference model only. We then obtain simple relations for each kernel:

$$K_{c^2,g}^{n,l} = K_{c^2,\rho}^{n,l} \qquad (10)$$

$$4\pi r^2 \rho \int_{r}^{R} \frac{K_{g,c^2}^{n,l}}{m} = K_{\rho,c^2}^{n,l}. \qquad (11)$$

The proposition of Masters (1979) was to solve directly this integral relation, which can be done using, for example, an iterative relaxation method to solve the integral equation. In practice, we favour a more efficient approach by deriving a differential equation for these kernels, simply by taking the derivative of Eq. (11) after having divided it by $r^2\rho$. We then obtain the following very simple differential equation:

$$-m \frac{d}{dr} \left( \frac{K_{\rho,c^2}^{n,l}}{\rho r^2} \right) = K_{g,c^2}^{n,l}, \qquad (12)$$

Since this equation is extremely simple, the kernels are straightforward to obtain. However, this development was just for the

<footnote n="2">2 More generally, this function could be written $f(\rho, P, g, \Gamma_1)$ or $f(\rho, P, g, Y)$ with any f that can be written in terms of linear perturbations of these quantities.</footnote>

<footnote n="3">3 We describe how this can be done in Appendix A.</footnote>
sake of illustration and a good example of the difficulties associated with this method is illustrated in Appendix A.

A more elaborated case, which has already been involved in practical applications is that of the $(u, \Gamma_1)$ and $(u, Y)$ kernels. These kernels are obtained by solving a second order differential equation which is recalled here:

$$ -y \frac{d^2 \mathcal{K}'}{(dy)^2} + \left[ \frac{2\pi y^{3/2} \tilde{\rho}}{\tilde{m}} - 3 \right] \frac{d \mathcal{K}'}{dy} = y \frac{d^2 \mathcal{K}}{(dy)^2} \\ - \left[ \frac{2\pi y^{3/2} \tilde{\rho}}{\tilde{m}} - 3 + \frac{\tilde{m} \tilde{\rho}}{2y^{1/2} \tilde{P}} \right] \frac{d \mathcal{K}}{dy} \\ + \left[ \frac{\tilde{m} \tilde{\rho}}{4y \tilde{P}^2} \frac{d \tilde{P}}{dx} - \frac{\tilde{m}}{4y \tilde{P}} \frac{d \tilde{\rho}}{dx} - \frac{3}{4y^{1/2} \tilde{P}} \frac{d \tilde{P}}{dx} - \frac{\tilde{m} \tilde{\rho}}{2y^{3/2} \tilde{P}} \right] \mathcal{K}, \quad (13) $$

with $\mathcal{K} = \frac{K_{u,\Gamma_1}^{n,l}}{r^2 \rho}$ and $\mathcal{K}' = \frac{K_{\rho,\Gamma_1}^{n,l}}{r^2 \rho}$ in the case of the differential equation of the $(u, \Gamma_1)$ kernels or with $\mathcal{K} = \frac{K_{u,Y}^{n,l}}{r^2 \rho}$ and $\mathcal{K}' = \frac{K_{\rho,Y}^{n,l}}{r^2 \rho}$ for the equation of the $(u, Y)$ kernels. Furthermore, in Eq. (13), one has the following definitions: $x = \frac{r}{R}$, $y = x^2$, $\tilde{m} = \frac{m}{M}$, $\tilde{\rho} = \frac{R^3 \rho}{M}$, $\tilde{P} = \frac{R^4 P}{GM}$. We also recall here that using kernels such as the $(\rho, Y)$ or the $(u, Y)$ kernels requires to introduce the equation of state by using the following definition:

$$ \frac{\delta \Gamma_1}{\Gamma_1} = \left( \frac{\partial \ln \Gamma_1}{\partial \ln P} \right)_{Z,Y,\rho} \frac{\delta P}{P} + \left( \frac{\partial \ln \Gamma_1}{\partial \ln \rho} \right)_{Z,Y,P} \frac{\delta \rho}{\rho} + \left( \frac{\partial \ln \Gamma_1}{\partial Y} \right)_{Z,P,\rho} \delta Y \\ + \left( \frac{\partial \ln \Gamma_1}{\partial Z} \right)_{Y,P,\rho} \delta Z. \quad (14) $$

Two hypotheses are made when using helium kernels. First, one assumes that the equation of state of the reference model and that of the target model are the same. Secondly, one usually drops the last term in $\delta Z$ of Eq. (14). This is often considered to be a benign hypothesis and we will review its impact for various kernels in Sect. 4.

The problem of this method is that, when deriving the differential equation for the kernels, one may be faced with discontinuous terms within the equation coefficients. These discontinuities are due to the effects of the transition from radiative regions to convective regions and have to be treated correctly if one does not wish to introduce numerical errors in the resolution. Typically, these discontinuities appear when taking first or second derivatives of the density (or any quantity related to the density through an algebraic relation). For example, the second derivative of the adiabatic squared sound speed, $c^2$ shows a discontinuity at the base of the convective envelope. This also means that the differential equation must be solved on separated domains and that continuity conditions have to be applied for each sub-domain. These conditions typically serve as constraints to solve the differential equations of structural kernels. For example, for the $(u, \Gamma_1)$ and $(u, Y)$ kernels, the resolution of the second order differential equation uses one central boundary condition that is derived from the differential equation itself and one “boundary” condition that stems from the integral equation. Namely, one assumes that the kernels have to satisfy their integral equation at some point of the sub-domain. For the next sub-domain, a continuity relation on the kernels is derived, since they have to be continuous for continuous variables, and the integral relation is again used to obtain an additional condition for the sub-domain.

In practice, the use of the integral relation for the additional condition is not trivial, since sometimes one can be confronted

with integrals of the layers above the layer on which one wishes to solve the differential equation$^4$. The problem is even more complicated when facing separated domains. Thus, one has to find a workaround based on the linearity of the problem and ends up solving a system of two differential equations on each sub-domain, where the equations are simultaneously connected through continuity relations and integral equations. With a little algebra, this can be done using a simple direct solver and finite difference discretization (in our case, we used the prescriptions of Reese (2013) for the grid on which the equation is solved). This leads in practice to a good accuracy in the results when care is taken in the computation of the derivatives of the coefficients and of the already known kernels. Indeed, these derivatives can be a source of significant numerical noise when calculated on a reference model of poor quality or when the eigenfunctions have been computed with a poor accuracy.

## 3.2. Method of conjugated functions – application to the A–Y kernels

The method of conjugated functions is quite different from what is done in the direct method, although the starting point is still the equality of two variational expressions for different structural variables. This method was presented for the first time in a paper by Elliott (1996) in the context of helioseismology and a more thorough presentation of the method can be found in Kosovichev (1999). In this section, we recall the basis of the method and apply it to the derivation of new kernels associated with the $(A, Y)$ structural pair.

The quantity A is called the convective parameter and is closely related to the Brunt-Väisälä frequency. It is defined as follows:

$$ A = \frac{d \ln \rho}{d \ln r} - \frac{1}{\Gamma_1} \frac{d \ln P}{d \ln r}. \quad (15) $$

This quantity has the convenient property to be zero in adiabatically stratified convective regions. It is also very sensitive to changes in depth of the base of the convective zone and changes in upper regions of convective envelopes. The problem we will define is thus related to determining the change of structural variables from $(\rho, Y)$ kernels to the $(A, Y)$ kernels. The $(\rho, Y)$ pair is a convenient starting point but one could choose other starting variables. Thus, our goal is to find the functions $K_{A,Y}^{n,l}$ and $K_{Y,A}^{n,l}$ for a given stellar model such that for two models that are sufficiently close to each other, we have:

$$ \begin{aligned} \frac{\delta Y^{n,l}}{\nu^{n,l}} &= \int_0^R K_{\rho,Y}^{n,l} \frac{\delta \rho}{\rho} dr + \int_0^R K_{Y,\rho}^{n,l} \delta Y dr \\ &= \int_0^R K_{A,Y}^{n,l} \delta A dr + \int_0^R K_{Y,A}^{n,l} \delta Y dr. \end{aligned} \quad (16) $$

First, we have to relate the linear perturbation of A to the other structural variables. In the approach of conjugated functions, one starts by defining a system of differential equations between the model perturbations, where one wishes to relate the different perturbed structural quantities found in the starting and final integral relations. In this particular case, one has a system of 3 differential equations that relates all the quantities together. This system

<footnote>$^4$ Thus the arguments of the integral are unknown.</footnote>
is written:

$$r \frac{d}{dr} \left( \frac{\delta \rho}{\rho} \right) = \delta A + \frac{Gm}{rc^2} \frac{\partial \ln \Gamma_1}{\partial Y} \Big|_{P,\rho} \delta Y + \frac{Gm}{c^2 r} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln P} \Big|_{\rho,Y} + 1 \right] \frac{\delta P}{P}, \\ + \frac{Gm}{c^2 r} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln \rho} \Big|_{P,Y} - 1 \right] \frac{\delta \rho}{\rho} - \frac{Gm}{rc^2} \frac{\delta m}{m}, \quad (17)$$

$$r \frac{d}{dr} \left( \frac{\delta m}{m} \right) = \frac{4\pi r^3 \rho}{m} \left( \frac{\delta \rho}{\rho} - \frac{\delta m}{m} \right), \quad (18)$$

$$r \frac{d}{dr} \left( \frac{\delta P}{P} \right) = \frac{Gm\rho}{rP} \left[ \frac{\delta P}{P} - \frac{\delta \rho}{\rho} - \frac{\delta m}{m} \right]. \quad (19)$$

As we will show in this section, the method of conjugated functions uses equations closely related to the system presented above. A major advantage is that this approach leads to a system with simple coefficients, for which the problem of numerical derivatives will not be as important as for the direct method. However, this method uses more hypotheses than the direct method and is consequently less well-suited for asteroseismology. Typically, the problem stems from the boundary conditions that are used to close the system and select a unique solution. For the surface boundary conditions, we have to assume that the mass of the observed target is the same as the mass of the reference model. At first, we thought that only the mean density was required to be fitted to ensure a verification of the variational expression but we will see how we were mistaken in the next section. Indeed, it can be argued that kernels for structural pairs such as the $(A, \Gamma_1)$ pair or the $(A, Y)$ pair will never offer a good accuracy in the asteroseismic case as will be illustrated in Sect. 4.

The goal of the method of conjugated functions is to determine the unknown tridimensional vector $v = (v_1, v_2, v_3)$, which is a conjugated function linked to the structural kernels (see Eq. (33)), solution of the following system (related to the system of Eqs. (17) to (19)):

$$-r \frac{dv}{dr} - v = \mathcal{A}^T v + C^T \mathbf{K}_1, \quad (20)$$

where we have used the following definitions:

$$\mathcal{A} = \begin{pmatrix} \frac{Gm}{rc^2} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln \rho} \Big|_{P,Y} - 1 \right] & - \frac{Gm}{rc^2} & \frac{Gm}{rc^2} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln P} \Big|_{\rho,Y} + 1 \right] \\ \frac{4\pi r^3 \rho}{rP} & - \frac{Gm\rho}{rP} & 0 \\ - \frac{Gm\rho}{rP} & - \frac{Gm\rho}{rP} & \frac{Gm\rho}{rP} \end{pmatrix}, \quad (21)$$

$$\mathbf{K}_1 = (K_{\rho,Y}^{n,l}, K_{Y,\rho}^{n,l}), \quad (22)$$

$$C = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \quad (23)$$

We also introduce the following definitions:

$$\mathbf{x} = \begin{pmatrix} \frac{\delta \rho}{\rho} \\ \frac{\delta m}{m} \\ \frac{\delta P}{P} \end{pmatrix}, \quad \mathbf{s}_1 = \begin{pmatrix} \frac{\delta \rho}{\rho} \\ \delta Y \end{pmatrix}, \quad \mathbf{s}_2 = \begin{pmatrix} \delta A \\ \delta Y \end{pmatrix}, \quad (24)$$

$$\mathcal{B} = \begin{pmatrix} 1 & \frac{Gm}{rc^2} \frac{\partial \ln \Gamma_1}{\partial Y} \Big|_{P, \rho} \\ 0 & 0 \\ 0 & 0 \end{pmatrix}, \quad (25)$$

$$\mathbf{K}_2 = (K_{A,Y}^{n,l}, K_{Y,A}^{n,l}), \quad (26)$$

$$\mathcal{D} = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}. \quad (27)$$

We use the following boundary conditions for $r = 0$:

$$3\bar{v}_1(0) + 3\bar{v}_2(0) = -\frac{K_{\rho,Y}^{n,l}}{r^2 \rho}(0), \quad (28)$$

$$\bar{v}_3(0) = 0, \quad (29)$$

with $\bar{v}_i = \frac{v_i}{r^2 \rho}$. Using $\bar{v}_i$ as variables for the system is motivated by the central limit of the structural kernels, such as $K_{\rho,Y}^{n,l}$, which goes as $O(r^2)$ in central regions. These boundary conditions can be obtained from the limit as $r$ goes to 0 of Eq. (20) itself similarly to what is presented in Unno et al. (1989) for the pulsation equations. The last boundary condition of Eq. (20), at $r = R$ is defined as follows:

$$\frac{\delta \rho}{\rho}(R)v_1(R) + \frac{\delta m}{m}(R)v_2(R) + \frac{\delta P}{P}(R)v_3(R) = 0, \quad (30)$$

and results from the elimination of surface term in the integration by parts in Eq. (39) which can be changed using l'Hospital's theorem to avoid having to define $\frac{\delta P}{\rho}(R)$ with the equation of hydrostatic pressure, thus considering both $P(R)$ and $\delta P(R)$ to be 0:

$$\frac{\delta \rho}{\rho}(R)v_1(R) + \frac{\delta m}{m}(R)v_2(R) + \left( \frac{\delta \rho}{\rho}(R) + \frac{\delta m}{m}(R) \right) v_3(R) = 0. \quad (31)$$

The main problem with this equation is that both the $\delta m$ and $\delta \rho$ terms are unknown, it is thus impossible to derive a simple boundary conditions and the components of $v$ without an additional hypothesis. In helioseismology, one states that the mass of the observed target is known and one ends up with a simple equation only with $\delta \rho$. One then simplifies the $\delta \rho$ term and ends up with the following simple relation:

$$v_1(R) + v_3(R) = 0. \quad (32)$$

This problem is intrinsic to the method of conjugated functions since one defines the boundary conditions of the system using an expression containing structural perturbations.

Using these definitions, it can be proven that if $v$ is the solution to this problem, then the structural kernels $K_{A,Y}^{n,l}$ and $K_{Y,A}^{n,l}$, for each perturbation of the model defined by $x$ and $s_2$, can be determined using the following relations:

$$K_{A,Y}^{n,l} = v_1, \quad (33)$$

$$K_{Y,A}^{n,l} = K_{Y,\rho}^{n,l} + K_{A,Y}^{n,l} \frac{Gm}{rc^2} \frac{\partial \ln \Gamma_1}{\partial Y} \Big|_{\rho,\rho}, \quad (34)$$

To demonstrate this property, let us first rewrite the system of Eqs. (17) to (19) in its vector form using the definitions we have just introduced:

$$r \frac{dx}{dr} = \mathcal{A}x + \mathcal{B}s_2. \quad (35)$$

One can also write a trivial matrix relation between vectors $s_1$ and $s_2$:

$$s_1 = Cx + Ds_2. \quad (36)$$

We now apply the scalar product of Eqs. (35) and (36) with $v$, defining the scalar product on the functional space as:

$$\langle a, b \rangle = \int_{0}^{R} a(r)b(r)dr, \quad (37)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of various kernels for the (A, Y) structural pairs for various degrees and radial orders.</figcaption></figure>

which is done in this case for each component of **v** and **x**. We
then obtain:

$$
\left\langle \mathbf{v}, r \frac{d\mathbf{x}}{dr} \right\rangle = \langle \mathbf{v}, \mathcal{A}\mathbf{x} \rangle + \langle \mathbf{v}, \mathcal{B}s_2 \rangle, \qquad (38)
$$

$$
- \left\langle r \frac{dv}{dr} + \mathbf{v} \cdot \mathbf{x} \right\rangle + [rv \cdot \mathbf{x}]_{0}^{R} = \langle \mathbf{v}, \mathcal{A}\mathbf{x} \rangle + \langle \mathbf{v}, \mathcal{B}s_2 \rangle, \quad (39)
$$

where we have applied an integration by parts and thus obtained
a differential equation for **v**. If one considers both Eqs. (39)
and (20), we obtain:

$$
\langle \mathbf{K}_1, \mathbf{C}x \rangle = \langle \mathbf{v}, \mathcal{B}s_2 \rangle. \quad (40)
$$

The new kernels can then be determined using:

$$
\langle \mathbf{K}_2, s_2 \rangle = \langle \mathbf{K}_1, \mathbf{C}x \rangle + \langle \mathbf{K}_1, \mathcal{D}s_2 \rangle \quad (41)
$$

$$
= \langle \mathbf{v}, \mathcal{B}s_2 \rangle + \langle \mathbf{K}_1, \mathcal{D}s_2 \rangle, \quad (42)
$$

where we have used Eq. (40). If we develop the scalar products,
we obtain the following integral relations:

$$
\int_0^R K_{A,Y}^{n,l} \delta A dr + \int_0^R K_{Y,A}^{n,l} \delta Y dr = \int_0^R v_1 (\delta A + \frac{Gm}{rc^2} \frac{\partial \ln \Gamma_1}{\partial Y}) |_{P_\rho} \delta Y dr \\
+ \int_0^R K_{Y,\rho}^{n,l} \delta Y dr. \quad (43)
$$

From these relations, we directly obtain the relations (33)
and (34) and have thus demonstrated that determining the vec-
tor **v** satisfying Eq. (20) allowed us to determine the kernels of
the (A, Y) structural pair.

However, a few comments must be made on Eq. (31) since
it leads to a strong limitation in the use of the method of con-
jugated functions. As previously stated, the boundary conditions
applied are that the mass of the observed target and the reference
model are the same. In asteroseismology, this is not necessarily
the case. When using this method for other kernels, we could
avoid this limitation by using the relation $\frac{\delta m}{m} = \frac{\delta \bar{\rho}}{\bar{\rho}}$ if the radius
is fixed. Ultimately, one ends up with the same implicit scaling

presented before for the direct method. It is a considerable ad-
vantage of the direct method that it does not explicitly uses any
hypothesis on the mass of the observed target.

In this particular case, scaling the perturbations is impossi-
ble since the quantities δA and δY are adimensional and not ex-
pressed as relative perturbations, obviously because $\frac{\delta A}{A}$ would
be undetermined when A goes to zero for the reference model.
Consequently, the trick of the implicit scaling cannot be used
and we are limited by the accuracy of radii determinations for
asteroseismic targets.

However, even with scaled models, the problem can still be
present for the helium integrals. Indeed, for the density or the
sound speed, the link is quickly done since these variables are
explicitly part of what is called the acoustic structure of the stel-
lar model and are directly linked to the oscillation frequencies.
The question is more difficult when one thinks about the helium
mass fraction. The problem is to link the helium mass fraction
profile of the scaled target model to the helium profile of the real
target. As such, there is no clear link between both profiles and
helium cannot be directly related to the dynamical time since it is
not an explicit variable of the acoustic structure. Therefore, cau-
tion as to be taken when determining helium abundances from
inversion techniques when there is no strong constraints on the
radius<fnref n="5" />.

We illustrate the (A, Y) kernels in Fig. 1 for various degrees and radial orders. It should be noted that the kernels associated with the convective parameter A of the (A, Γ₁) structural pair are quite similar to the kernels associated with A for the (A, Y) structural pair and could thus be used to carry out inversions of similar indicators without the need to introduce the equation of state in the problem. The main problem is then to cope with the high amplitude of the cross-term kernels but ultimately, the presence of pairs of kernels with similar behaviours can be used to check the robustness of the inversion for observed data since it should lead to similar results if the cross-term is properly damped for

<footnote n="5">⁵ As such the mass of the model would not be a problem if one considers that the mean density can be very accurately determined using seismology. Thus, if one knows the radius accurately, an accurate estimate of the mass can be determined provided good seismic data.</footnote>
<table><caption>**Table 2.** Physical ingredients of the target models used for the hare-and-hounds exercises.</caption><thead><tr><th></th><th>Target model 1</th><th>Target model 2</th><th>Target model 3</th><th>Target model 4</th></tr></thead><tbody><tr><td>Mass (M<sub>☉</sub>)</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><td>Radius (R<sub>☉</sub>)</td><td>1.0712</td><td>1.0822</td><td>1.0394</td><td>1.0770</td></tr><tr><td>Age (Gyr)</td><td>5.0</td><td>5.0</td><td>5.0</td><td>4.5</td></tr><tr><td>EOS</td><td>CEFF</td><td>OPAL</td><td>CEFF</td><td>OPAL</td></tr><tr><td>Abundances</td><td>GN93</td><td>GN93</td><td>AGSS09</td><td>AGSS09</td></tr><tr><td>X<sub>0</sub></td><td>0.7</td><td>0.7</td><td>0.7</td><td>0.67</td></tr><tr><td>X<sub>0</sub></td><td>0.015</td><td>0.015</td><td>0.015</td><td>0.02</td></tr><tr><td>α<sub>MLT</sub></td><td>1.7</td><td>1.7</td><td>1.7</td><td>1.7</td></tr><tr><td>Mixing</td><td>—</td><td>—</td><td>—</td><td>Settling+turbulent diffusion</td></tr></tbody></table>

both structural pairs. One additional striking feature of the (A, Y) structural pair is the high amplitude of the helium kernels when compared to those of the convective parameter. It is pretty unusual since as was already noticed for the (ρ, Y) structural pairs and confirmed for other structural kernels we derived, the helium kernels tend to have very low amplitudes and are thus very well adapted as cross-terms of inversions<fnref n="6" />. This makes these kernels very interesting for inversions of helium abundance using appropriated indicators in the solar case, where the data is abundant and the radius of the observed target is very well constrained and used to build the standard solar model used as a reference for the inversion. In terms of numerical quality, the verification of the initial system of differential equations is done up to relative differences of the order of $10^{-14}$ on the average. Typically the resolution is more accurate ($10^{-16}$ or less) in central regions and less accurate at the surface ($10^{-13}$). It should be noted that the numerical quality of the results is naturally still subject to the number of points of the models and the variables and unknowns considered in the system of equations.

## 4. Numerical experiments

In this section, we describe a few numerical experiments carried out to analyse the importance of various hypotheses used to compute structural kernels. All models were computed using the Clés stellar evolution code (Scuflaire et al. 2008b) with the following ingredients: the CEFF equation of state (Christensen-Dalsgaard & Daeppen 1992), the OPAL opacities from Iglesias & Rogers (1996), supplemented at low temperature by the opacities of Ferguson et al. (2005) and the effects of conductivity from Potekhin et al. (1999) and Cassisi et al. (2007). The nuclear reaction rates are those from the NACRE project (Angulo et al. 1999), supplemented by the updated reaction rate from Formicola et al. (2004) and convection was implemented using the classical, local mixing-length theory (Böhm-Vitense 1958). We also used the implementation of microscopic diffusion from Thoul et al. (1994), for which three groups of elements are considered and treated separately: hydrogen, helium and the metals (all considered to have diffusion speeds of $^{56}$Fe). The oscillation frequencies and eigenfunctions were computed using the Liège adiabatic oscillation code (Scuflaire et al. 2008a). We took much care to analyse the numerical quality of the eigenfunctions and the models before computing structural kernels. Irregularities and poor quality of the computed eigenfunctions can bias the results and lead to wrong structural kernels and thus wrong inferences from inverted results. From our experience in hare-and-hounds exercises and

inversions, we have determined that adding seismic constraints to the model is very efficient at bringing the reference model into the linear regime thus validating the inversion process. In other words, fitting the average large and small frequency separations is already a big improvement in terms of linearity, although individual seismic constraints, such as individual frequency ratios and individual small frequency separations are the best way to maximise the chances of being in the linear regime. Individual large frequency separations can also be used, but due to their sensitivity to surface effects, they should not be used in observed cases. As such, since in this study we did not use very elaborate seismic fitting techniques, our tests serve the only purpose of isolating various contributions to the errors and to test various hypotheses usually done when carrying out structural inversions in the context of helio- and asteroseismology.

We started by computing 4 target models with different physical ingredients summarized in Table 2. Among these effects, we tested opacity changes, changes in the equation of state, the impact of the metallicity, the impact of individual abundance tables along with changes of typical parameters used for seismic fits such as the mixing-length parameter, $α_{MLT}$ and the hydrogen abundance. For each target, we computed reference models with the same mass and similar physical ingredients. To ensure that both target models and reference models had the same radius, we used a minimization algorithm to fit the mean density of the target model varying the age of the reference model. In other words, since on the main-sequence the radius is changed due to slight core contraction and envelope expansion, we could ensure with this simple method a straightforward fit of all targets. Of course, this approach is limited. For instance, a model which includes efficient microscopic diffusion or a completely different chemical composition will not be strongly constrained by the fit of the mean density and thus will surely not be lying in the linear regime. This should be kept in mind throughout the paper since it is not what is done in typical seismic studies where all the available information is used. The verification of the linear relations between frequencies and structural profiles is demonstrated by plotting the relative differences between the right-hand side and left-hand side of the linear integral relations, denoted $E_{s_1, s_2}^{n,l}$ defined as follows:

$$E_{s_1, s_2}^{n,l} = \frac{\frac{\delta y^{n,l}}{\gamma^{n,l}} - \left( \int_0^R K_{s_1, s_2}^{n,l} \frac{\delta s_1}{s_1} dr + \int_0^R K_{s_2, s_1}^{n,l} \frac{\delta s_2}{s_2} dr \right)}{\frac{\delta y^{n,l}}{\gamma^{n,l}}}, \quad (44)$$

with $s_1$ and $s_2$ being any of the structural variables for which structural kernels can be obtained. Using this approach offers a straightforward method to compare the validity of the linear relations for each kernel and each mode, pointing out possible weaknesses and inaccuracies. For each comparison, we used the modes with $l = 0, 1, 2$ and $3$ and $n$ between 6 and 41.

<footnote n="6">⁶ Although errors on the equation of state can be non-negligible at the levels of accuracy of helioseismology.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Left panel: verification of the linear integral relations between structure and frequencies for target 1 and a model with $\alpha_{\text{MLT}} = 1.5$ and $X_0 = 0.69$. Right panel: same as left panel but for a model with $\alpha_{\text{MLT}} = 1.9$ and $X_0 = 0.715$. We clearly see that different structural pairs do not satisfy the linear relations to within the same accuracy. Each curve corresponds to a given $\ell$ of the modes.</figcaption></figure>

## 4.1. Limits of the linear regime

First of all, we illustrate in Fig. 2 the verification of the linear relations between frequency differences and structural differences for various structural pairs for target 1 and two reference models, with slightly different $\alpha_{\text{MLT}}$ and $X_0$ values. The model associated with the left panel has $\alpha_{\text{MLT}} = 1.5$ and $X_0 = 0.69$, whereas the model used as a reference for the right panel has $\alpha_{\text{MLT}} = 1.9$ and $X_0 = 0.715$. We can also see that all structural pairs do not satisfy the linear relations to within the same accuracy. Typically, kernels for the convective parameter $A$ can be problematic, especially kernels of the $(A, Y)$ structural pair. This can mean that all perturbations of the quantities may not be in the linear regime, and that for some kernels, second order terms should be considered. Ultimately this can be the case for variables other than the convective parameter and the $(A, Y)$ kernels can sometimes satisfy the linear relations whilst the $(\rho, c^2)$ kernels do not. Two other problems of hare-and-hounds exercises using various kernels have to be mentioned: first, the insufficient numerical quality of the model and of the eigenfunction themselves; second, the changes of the parameters of the models can sometimes be inappropriate to test these relations and thus, inversion techniques. This means that we are intrinsically limited in our tests for robustness of inversions and that to some extent, other approaches could be sought to fully constrain the limitations of inversions in the context of asteroseismology. The first point is quite straightforward and linked to various problems that can be found in stellar evolution codes. For example, the quality of numerical derivatives, which is a function of both the derivation scheme that is used and the quality of the grid on which the model or the eigenfunction is computed. Another highly underestimated error is the finite accuracy with which a stellar evolutionary model satisfies hydrostatic equilibrium. In other words, the intrinsic consistency of thermodynamical quantities used to describe the acoustic structure of the model must be checked. To these two sources of errors, we must add the possible differences stemming from intrinsic methods used to compute the models in various stellar evolutionary codes.

Intrinsic non-linearity is a recurring problem when using the frequency structure relations. In Fig. 3, we illustrate the arguments of the structural integrals from the $(\rho, c^2)$ pair and $(A, Y)$ pair. The $\rho$ and $c^2$ arguments have very regular patterns naturally more concentrated towards the surface regions due to the higher amplitude of the kernels. Similarly, the amplitude of the $Y$ contribution in the lower right panel is only important in the surface regions. Although smaller than the other contributions, this helium integral is a factor 2 larger than the helium integral from the $(\rho, Y)$ structural pair. As we will see later, this has important implications for the limitations of the linear regime with the $(A, Y)$ kernels. In the lower-left panel of Fig. 3, we can see that the $A$ term is much more important in the surface regions, with a small contribution coming from the base of the convective envelope. This means that in practice, this structural pair might well be very sensitive to surface effects. From the numerical point of view, this means that to use the $(A, Y)$ or $(A, \Gamma_1)$ pair, a very good quality of the grid as well as of the structure equations in the uppermost regions of the model is necessary to avoid important numerical uncertainties. We emphasize here that being able to build the $(\rho, c^2)$ pair to within a good accuracy does not mean that numerical errors remain small when building new kernels from the existing ones. This is particularly true for the $(A, Y)$ and $(A, \Gamma_1)$ kernels but can also be seen for other pairs.

Another extreme is the case where the perturbation of certain thermodynamic quantities can be considered small and thus within the linear regime while other cannot. In this case, certain linear relations might be valid while others are not. The case can be illustrated with kernels related to helium. Let us take two models, with the same mass, radius, chemical composition and mixing-length parameter. In one of the models, we include microscopic diffusion but not with its full intensity by multiplying the diffusion speeds by a factor $D$ smaller than one (here for example, we chose 0.5). The surface helium abundance has significantly changed. We see a difference in mass fraction of the order of 0.025, in other words, nearly 10%. It is obvious that the changes cannot be considered small and it is then no surprise to see that these models are within the linear regime for the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Arguments of the integrals of the linear relations between frequency and structure. $\rho$ argument (upper left) and $c^2$ argument (upper right) from the $(\rho, c^2)$ pair. A argument (lower left) and $Y$ argument (lower right) from the $(A, Y)$ pair. Each curve corresponds to a given $\ell$ of the modes.</figcaption></figure>

$(A, \Gamma_1)$ kernels but not the $(A, Y)$ kernels. Again this means that caution is required when changing the structural pair in an inversion process and that usually, the validity of the linear regime can be assessed by using different reference models to carry out the inversion with one structural pair. Ultimately, if the inversion result, let us say, changes significantly with the structural pair that is used, then there is a problem with the inversion process. In the case of the $Y$ kernels, the problem can also arise due to the assumption that the equation of state is known, since it is used to derive the kernels. In these test cases, we always used the same equation of state for both target and reference model, except when it is specifically mentioned as in the other test cases below.

However, even when the equation of state is the same, we noticed that Eq. (17) is not always perfectly satisfied. If the same equation is written for the $(A, \Gamma_1)$ kernels, then the agreement is improved, meaning that some of the errors seen for the $(A, Y)$ kernels can be attributed to the verification of Eq. (14). This hypothesis has been tested and we clearly saw a disagreement between the left-hand side and the right-hand side of Eq. (14). This disagreement did not seem to arise from numerical uncertainties but rather from the intrinsic non-linearity of the equation, due for example to shifts in the ionisation zones that were not reproduced by the linear expansion with derivatives of $\Gamma_1$. Knowing this, it thus seems perfectly normal to see a stronger non-linear behaviour for the $(A, Y)$ kernels since they have a much higher helium contribution than other kernels. This of course implies limitations on direct helium determinations from kernel inversions and requires further investigation.

## 4.2. Effects of metallicity

The effect of metallicity are extremely important to quantify since they are often neglected when trying to assess the helium

abundance. Using a few models, we review the impact of small changes of metallicity on the verification of the linear relations. This impact is illustrated in Fig. 4. An important point to mention is that in asteroseismic observed cases, the metallicity is calculated through the spectroscopic observations of [Fe/H]. One then uses the sun as a reference but it should be emphasized that there is no agreement to this day on the solar metallicity and that this uncertainty as such has an impact on linear relations, especially when using $Y$ related kernels.

In Fig. 4, we can disentangle the impact of metallicity, since on the left-hand plot, both the target and reference models have the same $Z$, whereas on the right-hand plot, we changed the metallicity by 0.002. Of course, since the models do not have the same age, some changes can be seen due to intrinsic differences in the models, but it is still striking to see that the difference in $Z$ can have an impact in some contexts, which is in contradiction with what was previously believed. The case of the high frequency range of the $(A, Y)$ kernel is a very good illustration of how this can be a problem. However, we note that other kernels, such as the $(A, \Gamma_1)$ pair were affected by the changes in metallicity, but not as much, so the intrinsic differences coming from $Z$ is at least a few per cent. From inspection of the behaviour of the $(\rho, \Gamma_1)$ pair, we can say that the differences for the $(\rho, Y)$ pair also stem from intrinsic differences and not only from the term in $\delta Z$ neglected in Eq. (14). This means that metallicity can be extremely important for some fitting processes in terms of the validity of the linear structural relations due to the intrinsic differences that can be generated between the target and reference models.

We also plotted in purple the verification of the linear relations for a model with a metallicity change of 0.001. We can see that the errors are divided by approximately a factor 2. This means that the effect of the metallicity is rather global and goes beyond the neglect of the additional term in Eq. (14). This is
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Left panel: verification of the linear integral relation for a reference model which has the same metallicity as target 1. Right panel: effects of a 0.002 shift to the metallicity on the verification of the linear integral relations. The effects of a 0.001 shift to the metallicity for the (A, Y) kernels only is shown in purple and referenced with a *. Each curve corresponds to a given *l* of the modes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Left panel: verification of the linear integral relation for a reference model which has the same equation of state as target 1. Right panel: effects of a change from the CEFF equation of state on the OPAL equation of state to the verification of the linear integral relations. Each curve corresponds to a given *l* of the modes.</figcaption></figure>

further confirmed by the impact of the metallicity on the (A, Γ₁) structural pair, where the δZ contribution is not explicitly involved. Figure 4 also shows that the (ρ, c²) kernels are not affected by metallicity changes, as expected, thus leaving their diagnostic potential unaltered.

### 4.3. Effects of the equation of state

In Fig. 5, we illustrate the same plot as Fig. 4, but changing the equation of state of the target model to the Opal equation of state (Rogers & Nayfonov 2002). The reference models are built with the CEFF equation of state which is used to compute the derivatives of Γ₁ and derive consistent variational expressions. Figure 5 thus illustrates the impact that not knowing the equation of state of the target has on the verification of the linear structural relations. While there is some impact, it is not as large as expected.

Of course this does not mean in any case that the equation of state is not important for the linear integral relations used in inversions, but when compared to the impact of metallicity, it

seems that in this case Z has a larger impact. This is not to be generalized but means that we have to be careful with the approximations made and perhaps, in the case of the solar metallicity problem, both the uncertainties on metallicity and the equation of state have to be taken into account. One point worth mentioning about this test case is that both the CEFF and the OPAL equations of state are very similar for solar conditions, so the small impact is a result of similarities between theoretical equations of state and might not be representative of the differences between the true equation of state in the sun and one of the theoretical ones.

### 4.4. Effects of abundances and radii inaccuracies

To test the impact of microphysics, we computed the target model with the AGSS09 (Asplund et al. 2009) heavy elements mixture and the same metallicity as the reference model, computed with the GN93 abundances (Grevesse & Noels 1993). Indeed, this changes significantly the thermodynamic quantities
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Left panel: effect of mismatches in radii between the reference and target models. The relative values of the mismatches are respectively $10^{-4}$ for the blue dots, $10^3$ for the magenta dots and $3 \times 10^{-3}$ for the red dots. Right panel: effects of changes in the abundances of heavy elements for different structural pairs. Each curve corresponds to a given *l* of the modes.</figcaption></figure>

inside the star and affects significantly the opacity. As such, this test can be seen as one way to demonstrate that microphysics also has a large impact on the verification of the linear integral relation used to carry out inversions. From the right hand side panel of Fig. 6, we can see that all the kernels are affected by the microphysics. The effects are mostly seen for the ($\rho$, $c^2$), ($\rho$, $\Gamma_1$) and ($\rho$, $Y$) pairs, but the good results of the other pairs are likely due to chance since all profiles have been significantly affected by the modified microphysics. The problem of the radius fit that was discussed in Sect. 3.2 is illustrated on the left side panel of Fig. 6. The tendency is clearly seen since introducing progressively an error on the mean density produces an important error on the verification of the linear relations. The problem would be similar if one would consider the mean density to be known within an excellent accuracy but the mass to be unknown. In such case, a small error on the mass introduces an error of the order of $R^3$ thus an even larger departure from the linear integral relation. The problem is intrinsically due to the adimensional nature of A, meaning that it cannot be scaled to take into account our ignorance of the mass or radius of the target. Indeed, this effect is not seen for any kernel computed with the direct method if the proper scaling is applied to the structural variables. This leads to intrinsic limitations of the application of the ($A, \Gamma_1$) and (A, Y) kernels in asteroseismology. The problem may not be solved in this case by changing the fitting method since fitting seismic constraints may not always ensure a good fit of the radius of the observed target. Ultimately, these kernels are only limited to the very best asteroseismic targets for which excellent interferometric measures of the radii are available.

### 4.5. Effects of extra-mixing

The extra-mixing term is very common in stellar physics. It is used to introduce additional hydrodynamical processes not taken into account in standard stellar models. The problem of lithium abundances is a good illustration that some extra-mixing is actually taking place in real stars. Thus, it seems perfectly normal to ask the question whether the non-inclusion of additional mixing processes could affect the verification of linear integral relation

between frequencies and structural quantities. The answer to this question is illustrated in Fig. 7. To carry out this test case, we used target model 4 which includes turbulent diffusion in addition to microscopic diffusion. To test the robustness of the linear relations, we used reference models for which turbulent diffusion had been inhibited. For example, in the left panel of Fig. 7, the reference model had a slightly higher helium of 0.005 abundance and a less efficient turbulent diffusion. The verification of the linear relations is still good, but it seems that the ($A, \Gamma_1$) and the ($\rho, \Gamma_1$) structural pairs are strongly affected by the neglect of extra-mixing. This statement is confirmed when looking at the right panel of Fig. 7 for which the reference model has a nearly constant extra-mixing throughout all layers of the model. The ($\rho, c^2$) and ($u, \Gamma_1$) kernels seem not to be too much affected by extra-mixing. In general, the impact of extra-mixing is much reduced for the models we tested here. This statement, of course, only applies for physical conditions similar to solar and for the fitting process we use in these numerical tests. This does not mean that the problem could not reappear for models with convective cores, for which extra mixing could change significantly the evolutionary path and the acoustic structure.

## 5. Conclusion

In this paper, we analysed the verification of the linear integral relations between frequency and structural quantities frequently used in helioseismology for various structural pairs. In Sects. 3.1 and 3.2, we presented the two methods to change the structural variables in the linear integral relations. The direct method we present in this paper has the advantage of being more general than the method of conjugated functions, which explicitly uses the same radius for the observed target and the reference model. Although this is implied in the direct method, it is not used as a pre-requisite in the derivation of the equations leading to the kernels. Moreover, it has been shown that it is possible circumvent the problem by re-scaling the information provided by the inversion (see Buldgen et al. 2015; Basu 2003, for a discussion of this problem). In that sense, the method we propose offers a good alternative that is applicable in the asteroseismic case. However,
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Effect of extra mixing on the verification of linear structural relations for both models. The left panel is associated with a model with a slightly inhibited extra-mixing intensity, whereas the right panel has a nearly constant extra-mixing throughout the model layers, but with an even smaller intensity. Each curve corresponds to a given *ℓ* of the modes.</figcaption></figure>

we have also shown that the method we propose leads to some-
what complicated coefficients which can be difficult to derive if
the numerical quality of the model is not ensured.

Furthermore, in Sect. 3.2, we showed how the conjugated
functions approach could be used to derive (A, Y) kernels. These
kernels have the particularity of showing a very high sensitivity
to the chemical composition since they are the only structural
kernels for which the helium kernels have a higher amplitude
than those associated with the secondary variable of the struc-
tural pair. This property is extremely important in the context of
structural inversions since the amplitude of helium kernels was
the main motivation behind their use as a cross-term in helio-
seismology but also their main handicap for direct kernel-based
inferences of the helium profile using classical inversion tech-
niques such as the SOLA or the RLS method. However, by no
means would these inversions be independent of the equation of
state since it is introduced in the very equations leading to the
(A, Y) kernels.

In Sect. 4, we presented various experiments showing the in-
trinsic limitations in the linear regime of structural pairs. These
limitations can be due to numerical inaccuracies or to the in-
trinsic non-linear behaviour of different variables. The most
striking example is that of helium, for which extra-mixing can
change significantly the local abundance while hardly changing
the sound speed or density profile. In that sense, the numerical
experiments we presented, although intrinsically limited, show
that changing the structural pair in the integral relations is not
often innocent, especially at the verge of non-linearity<sup>7</sup>.

In addition, we analysed the importance of various structural
changes, such as the impact of metallicity changes. We showed
that small changes of metallicity could affect significantly the
linear structural relations, especially for the (A, Y) structural pair.
However, we stress here that the linear behaviour of the inte-
gral relations is strongly dependent on the fitting process. This
emphasizes again that in the asteroseismic case, all information
available should be used to ensure the verification of the linear
structural relations.

In addition to the effect of metallicity, we also analysed the
impact of the equation of state on the integral relations. Surpris-
ingly, we find them to be less important than previously stated
and even less important in some cases than the metallicity ef-
fects on the verification of the linear relations. This could have
an impact on the potential of inversion techniques in the solar
case. Changing the equation of state has little impact, and may
be the result of a bias from the fitting process we used or it may
be due to the similarity under solar conditions of most equations
of state.

We also analysed the importance of radius constraints and
constraints on the microphysics by changing the heavy elements
mixture. The test case on the radius inaccuracies shows the im-
portance of this additional constraint for kernels derived with
the method of conjugated functions, whereas the kernels derived
with the direct method are found to be more robust if the proper
scaling is applied when analysing the inversion results. How-
ever, we also emphasize that adimensional variables, like A or Y,
cannot be rescaled. The test case on the heavy elements mixture
showed the important sensitivity to microphysics in stellar mod-
els. As such, the reference and target models were quite different
and it is not surprising that linear structural relations are strongly
affected. However, from our experience in seismic modelling, we
know that these differences can be reduced by introducing addi-
tional constraints. Ultimately, models with different abundances
can be very similar in terms of thermodynamical quantities due
to compensations.

Finally, we also analysed the importance of additional mix-
ing acting during the evolution of the target model. We found
that for solar conditions, additional extra-mixing processes could
change slightly the verification of the linear relations, but that
these changes were not as significant as those obtained from in-
accuracies in metallicity, for example. We stress that this analy-
sis should be extended to other parts of the HR diagram, where
extra-mixing can have a more significant impact on the acous-
tic structure of stellar models, and thus on the verification of the
linear relations between frequencies and structure.

To conclude, the advent of the space photometry era and the
quality of data provided by past and upcoming space missions
will allow us to use new seismic approaches to extract efficiently
seismic information. However, it is still important to provide a

<footnote n="7">⁷ Which may well be the case in the context of asteroseismic inversions.</footnote>
theoretical framework for these methods, to test their limitations
and to determine what additional information (spectroscopic, in-
terferometric, ...) or methodological improvement are necessary
to enable the use of seismic inversions in asteroseismology. As
such, this study only gives answers to limited theoretical ques-
tions and is only one step towards the improvement of our use of
seismic information.

*Acknowledgements.* G.B. is supported by the FNRS (“Fonds National de la Recherche Scientifique”) through a FRIA (“Fonds pour la Formation à la Recherche dans l’Industrie et l’Agriculture”) doctoral fellowship. This article made use of an adapted version of InversionKit, a software developed in the context of the HELAS and SPACEINN networks, funded by the European Commissions’s Sixth and Seventh Framework Programmes.

References

Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Baglin, A., Auvergne, M., Barge, P., et al. 2009, in IAU Symp. 253, eds. F. Pont, D. Sasselov, & M. J. Holman, 71

Basu, S. 2003, Ap&SS, 284, 153

Basu, S., & Antia, H. M. 2008, Phys. Rep., 457, 217

Böhm-Vitense, E. 1958, Z. Astrophys, 46, 108

Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977

Buldgen, G., Reese, D. R., & Dupret, M. A. 2015, A&A, 583, A62

Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094

Chandrasekhar, S. 1964, ApJ, 139, 664

Chaplin, W. J., Lund, M. N., Handberg, R., et al. 2015, PASP, 127, 1038

Christensen-Dalsgaard, J., & Daeppen, W. 1992, A&ARv, 4, 267

Clement, M. J. 1964, ApJ, 140, 1045

Elliott, J. R. 1996, MNRAS, 280, 1244

Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585

Formicola, A., Imbriani, G., Costantini, H., et al. 2004, Phys. Lett. B, 591, 61

Gough, D. O., & Thompson, M. J. 1991, in Solar Interior and Atmosphere, eds. A. N. Cox, W. C. Livingston, & M. S. Matthews (Tucson: University of Arizona Press), 519

Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. Vangioni-Flam, & M. Casse, 15

Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943

Kosovichev, A. G. 1999, J. Comp. Appl. Math., 109, 1

Lynden-Bell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293

Masters, G. 1979, Geophys. J., 57, 507

Potekhin, A. Y., Baiko, D. A., Haensel, P., & Yakovlev, D. G. 1999, A&A, 346, 345

Rauer, H., Catala, C., Aerts, C., et al. 2014, Exper. Astron., 38, 249

Reese, D. R. 2013, A&A, 555, A148

Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064

Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390

Scuflaire, R., Montalbán, J., Théado, S., et al. 2008a, Ap&SS, 316, 149

Scuflaire, R., Théado, S., Montalbán, J., et al. 2008b, Ap&SS, 316, 83

Thoul, A. A., Bahcall, J. N., & Loeb, A. 1994, ApJ, 421, 828

Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars (Tokyo: University of Tokyo Press)
## Appendix A: Convective parameter kernels from the direct method

As stated in the core of the paper, kernels for the structural pair $(A, Y)$ or $(A, \Gamma_1)$ can be also obtained from the direct method. We give a few steps in the derivation of the third order differential equation that leads to these kernels and discuss a few problems regarding its numerical resolution. The first step is to introduce the helium and A perturbations in the linear integral equation.

$$ 
\begin{aligned}
\int_0^R K_{\rho,Y}^{n,l} \frac{\delta \rho}{\rho} dr + \int_0^R K_{Y,\rho}^{n,l} \delta Y dr &= \int_0^R K_{A,Y}^{n,l} \delta A dr + \int_0^R K_{Y,A}^{n,l} \delta Y dr \\
&= \int_0^R K_{A,Y}^{n,l} \left[ r \frac{d}{dr} \left( \frac{\delta \rho}{\rho} \right) + \frac{Gm}{rc^2} \frac{\delta m}{m} \right. \\
&\qquad \left. - \frac{Gm}{rc^2} \left( \frac{\partial \ln \Gamma_1}{\partial \ln \rho} \Big|_{\rho,Y} - 1 \right) \frac{\delta \rho}{\rho} \right. \\
&\qquad \left. - \frac{Gm}{rc^2} \left( \frac{\partial \ln \Gamma_1}{\partial \ln P} \Big|_{\rho,Y} + 1 \right) \frac{\delta P}{P} \right] dr \\
&\quad + \int_0^R \left[ K_{Y,A}^{n,l} - K_{A,Y}^{n,l} \frac{Gm}{rc^2} \frac{\partial \ln \Gamma_1}{\partial Y} \Big|_{\rho,P} \right] \delta Y dr.
\end{aligned}
\quad (A.1) $$

Then, we have to use the definition of hydrostatic pressure and mass and permute the integrals. We can already notice that the term with the derivative of density will be problematic and will require an integration by parts.

$$ 
\begin{aligned}
& \left[ K_{A,Y}^{n,l} r \frac{\delta \rho}{\rho} \right]_0^R - \int_0^R \frac{d(r K_{A,Y}^{n,l})}{dr} \frac{\delta \rho}{\rho} dr \\
& + \int_0^R 4\pi r^2 \rho \left[ \int_r^R \frac{G K_{A,Y}^{n,l}}{\bar{r} c^2} d\bar{r} \right] \frac{\delta \rho}{\rho} dr \\
& - \int_0^R \frac{G m \rho}{r^2} \left[ \int_r^R K_{A,Y}^{n,l} \frac{G m}{\bar{r} c^2} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln P} \Big|_{\rho,Y} + 1 \right] d\bar{r} \right] \frac{\delta \rho}{\rho} dr \\
& - \int_0^R 4\pi r^2 \rho \left[ \int_r^R \frac{G \rho}{\bar{r}^2} \left( \int_0^\bar{r} \frac{G K_{A,Y}^{n,l}}{\bar{r} c^2 P} \left( \frac{\partial \ln \Gamma_1}{\partial \ln P} \Big|_{\rho,Y} + 1 \right) d\bar{r} \right) d\bar{r} \right] \frac{\delta \rho}{\rho} dr \\
& - \int_0^R K_{A,Y}^{n,l} \frac{G m}{rc^2} \left[ \frac{\ln \Gamma_1}{\ln \rho} \Big|_{\rho,Y} - 1 \right] \frac{\delta \rho}{\rho} dr = \int_0^R K_{\rho,Y}^{n,l} \frac{\delta \rho}{\rho} dr.
\end{aligned}
\quad (A.2) $$

The first term of this equation is exactly zero for $r = 0$. However, it is not for $r = R$ and it is a strong hypothesis to consider that the surface relative density differences are exactly zero. Nevertheless, it is the only way to obtain an equation for the $(A, Y)$ kernels with the direct method. The problem is exactly the same for the $(N^2, c^2)$ kernels which can be easily derived but will also face the same problem due to the density derivative. Moreover, one could argue that the contribution of the surface term is negligible when compared to the integrals and that this term has no impact on the final result of the kernels<fnref n="8" />. This simplification could also be seen as a boundary condition, stating that the kernels we are searching for have to be exactly 0 at the surface boundary. Ultimately, after a few additional algebraic operations, we obtain a third order differential equation that we write here as a function of $K = \frac{K_{A,Y}^{n,l}}{r^2 \rho}$

and $K' = \frac{K_{A,Y}^{n,l}}{r^2 \rho}$.

$$ 
\begin{aligned}
& - \frac{r^6}{m} \frac{d^3 K}{dr^3} + \left[ \frac{4\pi r^8 \rho}{m^2} - \frac{10r^5}{m} - \frac{r^6}{m} \frac{d \ln \rho}{dr} - \frac{r^5 A_1}{m} \right] \frac{d^2 K}{dr^2} \\
& - \left[ \frac{20r^4}{m} - \frac{16\pi r^7 \rho}{m^2} + \frac{2r^6}{m} \frac{d^2 \ln \rho}{dr^2} + 2 \frac{r^5}{m} \frac{d A_1}{dr} + \frac{7r^5}{m} \frac{d \ln \rho}{dr} \right. \\
& \left. - \frac{4\pi r^8 \rho}{m^2} + \frac{5r^4 A_1}{m} - \frac{4\pi r^7 \rho A_1}{m^2} - \frac{4\pi r^6 \rho}{c^2 m} \right] \frac{d K}{dr} \\
& - \left[ \left[ \frac{5r^4}{m} - \frac{4\pi r^7 \rho}{m^2} \right] \left[ \frac{d \ln \rho}{dr} + \frac{d A_1}{dr} \right] + 4GA_2 \right. \\
& + \frac{r^5}{m} \frac{d^2 A_1}{dr^2} + \frac{r^6}{m} \frac{d^3 \ln \rho}{dr^3} + \frac{16\pi^2 r^8 \rho^2}{c^2 m^2} + \frac{8\pi r^6 \rho}{c^3 m} \frac{dc}{dr} \\
& - \frac{4\pi r^6}{c^2 m} \frac{d \rho}{dr} + \frac{24\pi r^5 \rho}{c^2 m} + \left[ \frac{7r^5}{m} - \frac{4\pi r^8 \rho}{m^2} \right] \frac{d^2 \ln \rho}{dr^2} \\
& \left. K' = \frac{r^5}{m} \frac{d^2 K'}{dr^2} + \left[ \frac{5r^4}{m} - \frac{4\pi r^7 \rho}{m} \right] \frac{d K'}{dr} \right]
\end{aligned}
\quad (A.3) $$

with the following additional definitions:

$$ A_1 = \frac{Gm}{rc^2} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln \rho} \Big|_{\rho,Y} - 1 \right] \quad (A.4) $$

$$ A_2 = \frac{Gm\rho r}{c^2 P} \left[ \frac{\partial \ln \Gamma_1}{\partial \ln P} \Big|_{\rho,Y} + 1 \right]. \quad (A.5) $$

Now the central boundary conditions are found using additional transformations. Typically, we solve the equation using $r^2$ as the independent variable for radial position. Taking the limit of the differential equation when $r^2$ goes to 0 and simplifying leads to simple central boundary conditions. Surface conditions are found by stating that Eq. (A.2) has to be satisfied and that the kernels must be 0 at the surface. One must also take care of the discontinuous coefficients, meaning that, again, the system must be solved in the radiative and convective regions independently and reconnected using proper continuity conditions. We face the same problem as for the $(u, Y)$ or $(u, \Gamma_1)$ kernels but can use the same algebraic manipulation to solve the system. Ultimately, we have to solve three discretized equations on two different domains (six systems in total) and reconnect those solutions. Now in addition to the numerical cost of such manipulations, we can see in Eq. (A.3) third derivatives of the density and second derivatives of the $A_1$ function. From numerical experiments not presented here, we have seen that these coefficients contain numerical noise due to the quality of the reference model. The noise can be reduced by increasing the quality of the model and of the finite difference scheme, but smoothing is still necessary to a certain extent. The concern with the smoothing process is that it could in some pathological cases change significantly the form of the kernels. In conclusion, despite the applicability of the method to the $(A, Y)$ kernels and the fact that it uses different hypotheses to obtain structural kernels, we state that for the $(A, Y)$ kernels, the direct method is not well suited in large scale automated studies. Moreover, due to the intrinsic problems of inversions with the adimensional variables mentioned before and the difficulties in determining accurate radii, the $(A, Y)$ pair might well be restricted to solar inversions for which the conjugated functions method is perfectly valid and should be preferred since it leads to simpler equations.

<footnote n="8">⁸ This could be done by analysing how accurate Eq. (A.2) is without the additional term.</footnote>