© P. S. Houdayer et al. 2021

1. Introduction

Asteroseismology, which is the study of resonant modes in stars, reveals information on the physical properties of the layers that a wave passes through on its way to the surface. Coupled with a physical model of the star, it thus allows us to constrain the various internal processes involved better than any other method. Undoubtedly, the constraints thus obtained depend on what is chosen to be modelled and what is assumed to be known. This choice is highly complex and largely dependent on the quality of the information provided (i.e., the available precision on the observables). Although high-precision photometric data provided by CoRoT¹ (Baglin et al. 2006), Kepler (Gilliland et al. 2010; Lund et al. 2017), and now TESS² (Ricker et al. 2015; Stassun et al. 2019) allow ever more precise estimates of oscillation frequencies, the evaluation of certain stellar parameters remains uncertain, however. This illustrates that accuracy rather than precision becomes the limiting factor in this case. Although the value of seismic inference in the determination of stellar parameters is undeniable, it also opens the door to potential biases on these very parameters (some of them are illustrated in Fig. 1).

Fig. 1. Schematic illustration of several potential biases involved in the inference problem. Applied to the abundance determination problem, the observable space represents (notably) the oscillation frequency sets, while the unobservable space would represent the model abundances. Case a shows an ideal scenario, where for given constraints (blue circle) our modelling (in grey) provides unbiased abundances (red circle). Panel b: possible expression of physical processes that are missing or are poorly taken into account in the model. It thus results in physical biases although it has ideal observables. Case c provides an example of degeneracy, that is, parameters that can vary widely because they are too weakly constrained by the observables. The last panel d is a typical illustration of the solar problem, where no mixture seems to satisfy the constraints. In this case, inconsistencies in the model physics or much weaker constraints than actually provided by the oscillation frequencies must be considered in order to provide a solution.

In this respect, the determination of abundances using realistic stellar models constitutes a textbook case that is potentially susceptible to a wide range of biases. Identified sources of difficulties undoubtedly reside in the complexity of the physics that must be considered. In particular, the equation of state, opacities (Kosovichev et al. 1992), or transport processes such as diffusion, turbulent mixing (Christensen-Dalsgaard et al. 1993), or radiative acceleration (Deal et al. 2018) can be evoked, whose consideration or exclusion may result in physical biases. Recent studies using model grids also suggest a strong anti-correlation between the mass and initial helium abundance estimates (Lebreton & Goupil 2014; Noll et al. 2021). This degeneracy makes the inference more complex by resulting in a high volatility of both parameter estimates. Moreover, when a large number of frequencies are available, as for the Sun, inversion techniques highlight significant discrepancies, and solar models are forced to choose between inconsistent abundances, densities, or convective zone (CZ) depths (Basu & Antia 2004; Asplund et al. 2009; Serenelli et al. 2009). In addition to these processes that directly involve the composition, additional uncertainties surrounding the near-surface region are to be faced, namely surface effects (Christensen-Dalsgaard et al. 1988). As a result of theoretical developments (Canuto 1997; Belkacem et al. 2021) and 3D hydrodynamical simulations (Belkacem et al. 2019; Schou & Birch 2020), our understanding of the mechanisms involved in these effects (as well as their contributions) has definitely improved. However, the more costly procedures (Mosumgaard et al. 2020) are generally avoided in favour of ad hoc corrections of the oscillation frequencies (Kjeldsen et al. 2008; Ball & Gizon 2014; Sonoi et al. 2015), which may affect stellar parameter estimates (Nsamba et al. 2018).

Methods allowing independent measurements of abundances (in particular helium) are therefore the subject of much investigation. The challenge is in fact twofold because a reliable estimate of abundances in turn would provide a constraint on the internal physics that depends on it.

The specific connection between composition and seismic properties of a star has long been known and can easily be understood through the structural change caused by ionisation if ever the abundances were to vary. The fact that this change is localised (in this case in the near-surface region) causes what is called a seismic glitch and manifests itself through an oscillatory deviation of the observed frequencies from a chosen reference (Gough 1990). Although this effect is not specific to the ionisation regions, these regions have benefited from numerous treatments, being both the most pronounced glitch and a potential marker of the helium abundance (Perez Hernandez & Christensen-Dalsgaard 1994; Lopes et al. 1997). From this point on, many studies (Monteiro et al. 1994; Basu et al. 1994; Monteiro & Thompson 1998, 2005; Gough 2002; Houdek & Gough 2007) considered various shapes of structural perturbations to analytically derive the expected frequency shifts. Most procedures exploit inversion formulae (Dziembowski et al. 1990; Gough & Thompson 1991; Antia & Basu 1994) and rely on an analytical modelling of the glitch to give a parameterised form of the oscillation. The expression thus obtained only needs to be fitted to data in order to provide information on the introduced parameters. In order to overcome dependences on irrelevant disturbances in the frequencies (e.g., surface effects in our case), it has also been proposed that the second differences, d²ν, be studied, rather than the frequencies directly (Gough 1990),

$$\begin{array}{c}\hfill d^2{\nu }_{n,\ell }={\nu }_{n+1,\ell }-2{\nu }_{n,\ell }+{\nu }_{n-1,\ell },\end{array}$$(1)

where n and ℓ are the associated oscillation radial order and degree, respectively. The seismic diagnostic defined in this way is less sensitive to surface effects and core perturbations while highlighting signal that would come from an intermediate acoustic depth (Ballot et al. 2004).

The main benefit of this (now usual) approach is to avoid complex modelling and therefore part of the issues mentioned above. Being based on a model with far fewer parameters, which reflects a local rather than global structure, the method seeks to make the best use of the information residing only in the low-degree frequencies. Additionally, this procedure is lighter than a direct minimisation or inversion of the frequency differences, which makes it convenient for application to large samples of stars (Verma et al. 2019). On the other hand, it would be a mistake to think that the procedure is “model independent” as it depends on the form of the perturbation that is considered. This procedure may therefore well lead to nonphysical or inaccurate frequency shifts if it relies on a nonphysical or inaccurate glitch modelling. Moreover, if the parameter to be estimated does not appear directly in the model, there is no alternative to calibration on stellar models (Houdek & Gough 2007, 2011; Verma et al. 2014a, 2019), which adds its own uncertainty on the internal physics. This method is therefore largely dependent on the work done beforehand, and this first paper in the series is intended to cover this perspective. A second paper will focus on the seismic effects of the structural perturbation and on the information these effects provide about the abundances.

In the present paper, we propose a physical model of the ionisation region that allows us to derive a semi-analytic description of the structural perturbation caused by a change in abundances. In this respect, we first briefly describe previous modelling in Sect. 2 and provide more detailed picture of the most commonly used model of the ionisation glitch, which is that of Houdek & Gough (2007), hereafter HG07. The formalism of our own model of the ionisation region is introduced in Sect. 3 and its structural behaviour is studied in Sect. 4. In Sect. 5 we propose an analysis of the structural perturbation induced by the model and determine how it is related to analytical models of the ionisation glitch. The last part of the article is dedicated to our conclusions.

2. Previous seismic glitch frameworks

As we briefly mentioned, a twofold approach in the modelling is required to study seismic glitches. First, an expression for the structural perturbation must be considered, and then the shape of the associated frequency shift must be inferred.

2.1. Structural perturbation modelling

The modelling of the perturbation usually boils down to choosing a localised and parametrisable shape for the induced structural perturbation. Depending on the region and the phenomenon under study, various forms can be considered. For the transition between radiative and convective regions, very sharp shapes such as a Dirac or a step function are generally used to reflect a discontinuity in the density derivatives (Monteiro et al. 1994; Houdek & Gough 2007). However, because the ionisation region is more spread out, the shapes that are used to model the structural perturbations usually involve an additional dispersion parameter (Monteiro & Thompson 2005; Houdek & Gough 2007). In Fig. 2 we represent shapes introduced in previous papers with their associated parameterisation. In what follows, we focus on the most commonly used model (i.e. the one presented in HG07) for estimating helium abundances and on the derived expressions for the frequencies.

Fig. 2. Various shapes along with their associated parameters that are used to describe a structural perturbation. (a) Dirac function used in Monteiro et al. (1994) to model the variations in the acoustic potential (in the case of an overshoot) passing from a convective to a radiative region. For this modelling, it is only necessary to specify the acoustic depth τd and amplitude Aδ of the perturbation. (b) Triangular shape of the Γ1 perturbation in the second helium ionisation region as used in Monteiro & Thompson (2005). An additional parameter controls the width (β) of the perturbation. (c) Gaussian shape of the Γ1 perturbation in the second helium ionisation region as used in Houdek & Gough (2007). ${\mathrm{\Delta }}_{\text{II}}^2$ and ΓII describe the variance and the area of the distribution, respectively.

In HG07, the structural change caused by a change in helium abundance is first modelled as a Gaussian perturbation of the first adiabatic exponent Γ₁ = (∂lnP/∂lnρ)_S,

$$\begin{array}{c}\hfill \frac{\delta {\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}=-\frac{1}{\sqrt{2\pi }}\frac{{\mathrm{\Gamma }}_{{}_{\mathrm{HeII}}}}{{\mathrm{\Delta }}_{{}_{\mathrm{HeII}}}}{e}^{-{(\tau -{\tau }_{{}_{\mathrm{HeII}}})}^2/2{\mathrm{\Delta }}_{{}_{\mathrm{HeII}}}^2}.\end{array}$$(2)

The idea is as follows: The transition from a pure-hydrogen model to one partially composed of helium causes the appearance of a well in the Γ₁ profile (cf. Fig 3). This results from the second helium ionisation to which the index HeII refers. While various analytical functions can be used to model a well, the well appears to be closely reproduced by a Gaussian expressed as a function of the acoustic depth $\tau ={\int }_r^R\text{d}{r}^{\prime }/c$; c(r), designating the sound speed at a given position. An important point is that if the depth, position, or width of the well depends on the helium abundance Y or on the thermodynamic conditions of the CZ, then a connection (though implicit) can be made between physical properties of the CZ and the parameters Γ_HeII, Δ_HeII, and τ_HeII.

Fig. 3. Representation of a typical Γ1 profile in the ionisation region as a function of the acoustic depth τ (the surface corresponds to τ = 0). The helium abundance in this plot is Y ∼ 0.25. The contributions of the three main ionisation zones have been distinguished, i.e., the hydrogen (H), the first (HeI), and second (HeII) helium ionisation zones. Each of them causes a deviation from the Γ1 reference value of 5/3.

Before further discussion, it may be useful here to clarify an ambiguity concerning the notation δ that symbolises the perturbation. Because it denotes a difference between two profiles (e.g., between profiles ${\mathrm{\Gamma }}_1^A$ and ${\Gamma }_1^B$), the variable over which they are expressed has to be chosen carefully. For instance, for arbitrary variables x and y,

$$\begin{array}{c}\hfill {\delta }_x{\mathrm{\Gamma }}_1\equiv {\mathrm{\Gamma }}_1^B(x)-{\mathrm{\Gamma }}_1^A(x)\ne {\mathrm{\Gamma }}_1^B(y)-{\mathrm{\Gamma }}_1^A(y)\equiv {\delta }_y{\mathrm{\Gamma }}_1.\end{array}$$(3)

As shown above, the notation δ_x allows us to overcome this ambiguity. This notation was introduced by Christensen-Dalsgaard & Thompson (1997), which developed this point extensively. Because the quantity x should vary in the same range for the two profiles A and B, it is natural to choose it as a normalised variable, for example, r/R, m/M, t = τ/τ₀, ...(R, M, and ${\tau }_0={\int }_0^R\text{d}r/c$ are the total radius, mass, and acoustic radius of the model). Moreover, the variable on which the perturbation is expressed can differ from the variable that was used to calculate the difference; both δ_xΓ₁(x) and δ_xΓ₁(y) have a meaning.

This point clarified, it appears that the perturbation used in HG07 for Eq. (2) is δ_τΓ₁ because the Gaussian shape has been established after an expansion at fixed τ (cf. Eq (31) of HG07). As mentioned, δ_τ is ambiguous, however, if ${\tau }_0^A\ne {\tau }_0^B$. It can be replaced by δ_t following

$$\begin{array}{c}\hfill \frac{{\delta }_t{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}=\frac{{\delta }_{\tau }{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}+\frac{{\delta }_t\tau }{\tau }\frac{\mathrm{d}ln{\mathrm{\Gamma }}_1}{\mathrm{d}ln\tau }=\frac{{\delta }_{\tau }{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}+\frac{\delta {\tau }_0}{{\tau }_0}\frac{\mathrm{d}ln{\mathrm{\Gamma }}_1}{\mathrm{d}ln\tau },\end{array}$$(4)

with $\delta {\tau }_0={\tau }_0^B-{\tau }_0^A$ simply being the perturbation of a constant.

2.2. Expected frequency shift

The resulting signature in the frequencies can be divided in two distinct parts: the ionisation component, and the smooth component. We describe them in turn below.

Ionisation component. The connection between Eq. (2) and the frequency shift $\delta {\nu }_{n,\mathcal{l}}={\nu }_{n,\mathcal{l}}^B-{\nu }_{n,\mathcal{l}}^A$ is made by considering the following asymptotic form (Gough 1990):

$$\begin{array}{c}\hfill \frac{\delta {\nu }_{n,\ell }}{{\nu }_{n,\ell }}={\displaystyle {\int }_0^1[{\mathcal{K}}_{\rho c^2}^{n,\ell }\frac{{\delta }_x\rho }{\rho }+{\mathcal{K}}_{c^2\rho }^{n,\ell }\frac{{\delta }_x{c}^2}{c^2}]\mathrm{d}x\sim {\int }_0^1{\mathcal{K}}_{c^2\rho }^{n,\ell }\frac{{\delta }_x{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}\mathrm{d}x,}\end{array}$$(5)

where ${\mathcal{K}}_{\rho c^2}^{ n,\mathcal{l}}$ and ${\mathcal{K}}_{c^2\rho }^{ n,\mathcal{l}}$ denote the usual structural kernels that can be found in Gough & Thompson (1991). This derivation, as well as the majority of the frequency shift treatment, will be fully discussed in a second paper. We emphasise that the perturbation used in Eq. (5) must be applied with caution. It can be shown that

$$\begin{array}{c}\hfill \frac{{\delta }_x{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}=\frac{{\delta }_t{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}-[\frac{\delta {\tau }_0}{{\tau }_0}-\frac{\delta R}{R}+\frac{1}{x}{\displaystyle {\int }_0^x}\frac{{\delta }_{t}c}{c}\mathrm{d}{x}^{\prime }]\frac{\mathrm{d}ln{\mathrm{\Gamma }}_1}{\mathrm{d}lnr}.\end{array}$$(6)

This nuance cannot be appreciated without the notation introduced here, however.

The frequency shift derived in HG07 is written as a continuous function of frequency,

$$\begin{array}{c}\hfill \frac{\delta \nu }{\nu }=A_{{}_{\mathrm{HeII}}}\nu e^{-8{\pi }^2{\mathrm{\Delta }}_{{}_{\mathrm{HeII}}}^2{\nu }^2}cos2(2\pi {\tau }_{{}_{\mathrm{HeII}}}\nu +{ϵ}_{{}_{\mathrm{HeII}}}),\end{array}$$(7)

as are the second differences,

$$\begin{array}{c}\hfill d^2{\nu }_{{}_{{\mathrm{\Gamma }}_1}}=F_{{}_{\mathrm{HeII}}}{A}_{{}_{\mathrm{HeII}}}\nu e^{-8{\pi }^2{\mathrm{\Delta }}_{{}_{\mathrm{HeII}}}^2{\nu }^2}cos[2(2\pi {\tau }_{{}_{\mathrm{HeII}}}\nu +{ϵ}_{{}_{\mathrm{HeII}}})-{\delta }_{{}_{\mathrm{HeII}}}],\end{array}$$(8)

with A_HeII = πΓ_HeII/τ₀. F_HeII and δ_HeII are functions of the other parameters and the frequency (although it is assumed that they only fluctuate slowly with it). The total number of parameters is therefore only four.

Smooth component. The helium component is not the only perturbation expected in the second differences. In addition to the signature of the transition between radiative and convective regions, a smooth component d²ν_s can also be considered and modelled as

$$\begin{array}{c}\hfill d^2{\nu }_{\mathrm{s}}={\displaystyle \sum _{k=0}^3{a}_k{\nu }^{-k}.}\end{array}$$(9)

The idea is to first consider a star that would not contain any glitch. Without any structural perturbation (thus simplifying the problem), it can be assumed that its frequencies follow the asymptotic expansion provided by Tassoul (1980), that is,

$$\begin{array}{c}\hfill {\nu }_{n,\ell }=[\frac{1}{2}(2n+\ell +\epsilon )+\frac{2V_{\ell }}{2n+\ell +\epsilon }]\mathrm{\Delta }\nu +\mathcal{O}\left(\frac{1}{{{\nu }_{n,\ell }}^2}\right),\end{array}$$(10)

where ε and V_ℓ are two dimensionless parameters that are independent of the radial order n, and Δν ≡ 1/2τ₀ is the asymptotic large frequency separation. The last term of Eq. (9), which is proportional to ν⁻³, can then be interpreted as the leading term of the second differences applied to Eq. (10),

$$\begin{array}{c}\hfill d^2{\nu }_{\mathrm{asym}}\sim \frac{2V_{\ell }}{{\nu }^3}{(\mathrm{\Delta }\nu )}^4.\end{array}$$(11)

However, it is difficult to give as convincing a justification for the other terms, which hide a much more complex deviation to the asymptotic expansion. Houdek & Gough (2011) in particular listed hydrogen ionisation, a sharp stratification of the upper layers of the convection zone, non-adiabatic processes, and turbulent perturbations caused by the oscillations. In the end, this component adds four free parameters to those identified so far. It should be noted that the combination of positive powers of ν introduced in Verma et al. (2014a) to describe the smooth component is incompatible with the asymptotic expansion (10).

2.3. Improvements and limits of the approach

Although it is not apparent from this short summary, the derivation of Eqs. (5)–(8) is as complex as it is intricate, and an analytical derivation of second differences therefore requires making use of many approximations. The first consequence is that a fit of the second differences with Eq. (8) only allows a very approximate retrieval of the parameters introduced in Eq. (2). HG07 revealed typical discrepancies of 50, 33, and 25% on parameters Γ_HeII, Δ_HeII, and τ_HeII, respectively. These differences can be reduced to 5, 15, and 15% by introducing a frequency dependence in ϵ_HeII that is induced by the cut-off frequency and a more complete model by adding the first helium ionisation contribution in Eq. (2) (see Fig. 3). Having now two Gaussians (indexed by HeI and HeII), the expected second differences become

$$\begin{array}{c}\hfill d^2{\nu }_{{}_{{\mathrm{\Gamma }}_1}}={\displaystyle \sum _{i=\{\mathrm{HeI},\mathrm{HeII}\}}{F}_i{A}_i\nu e^{-8{\pi }^2{\mathrm{\Delta }}_i^2{\nu }^2}cos[2(2\pi {\tau }_i\nu +{ϵ}_i(\nu ))-{\delta }_i]}\end{array}$$(12)

instead of Eq. (8). Because the parameters of the two Gaussians are not independent, three of them are fixed in HG07 by the empirical relations A_HeI/A_HeII = 0.5, Δ_HeI/Δ_HeII = 0.9, and τ_HeI/τ_HeII = 0.7. The number of independent parameters in the expression is still four (only one is needed to determine ϵ_HeI(ν) and ϵ_HeII(ν)). However, despite the substantial improvement in the parameter retrieval, Eq. (8) remains the currently most frequently used equation for studying ionisation glitches (Verma et al. 2017, 2019; Farnir et al. 2019).

More importantly, the abundances, and in particular, the helium abundance Y, do not appear directly as parameters in Eq. (12). As mentioned above, without further theoretical work to establish a dependence between these parameters and Y, calibration based on realistic models therefore seems to be necessary in order to make it appear. This is probably the greatest criticism that can be made of the method. To illustrate this point, we place ourselves in a broader framework and consider two models 𝒫 and ℳ_⋆ that represent the structural perturbation caused by a change in helium δY and the structure of a realistic stellar model, respectively. The first is parameterised by a set of values ⃗θ_p and the second by the set ⃗θ_⋆ from which we explicitly distinguish Y for its particular role in this context. As an example, if the model 𝒫 is chosen to be a Gaussian, as done by HG07 (Eq. (2)), it would then involve three parameters ⃗θ_p = (Γ_HeII,Δ_HeII,τ_HeII). In contrast, ⃗θ_⋆ instead contains quantities needed to calculate a realistic model, such as fundamental parameters (mass M_⋆, radius R_⋆, age A_⋆, etc.), but also all the quantities required to model the physical processes involved in the star (overshoot, mixing length, parameters implied in diffusion or rotation, etc.; cf. Lebreton & Goupil 2014, which provides a broad idea of the various possible prescriptions). As a result, ⃗θ_⋆ generally contains many more components than ⃗θ_p. With this in mind, we would like to be able to relate a fitted set of parameters, ⃗θ_p, to a difference in helium δY from a chosen reference Y. For this purpose, the solution proposed by calibration is to assimilate the perturbation model 𝒫 to a model perturbation δℳ, and therefore to assume

$$\begin{array}{c}\hfill \mathcal{P}({\theta }_p)\simeq {\mathcal{M}}_{\star }(Y+\delta Y;{\theta }_{\star })-{\mathcal{M}}_{\star }(Y;{\theta }_{\star })\simeq {\frac{\partial {\mathcal{M}}_{\star }}{\partial Y}|}_{Y;{\theta }_{\star }}\delta Y\end{array}$$(13)

in order to numerically derive a relation ⃗θ_p(δY). We note that Eq. (13) is equivalent to Eq. (33) of HG07 when 𝒫 is chosen to be the sum of two Gaussians and Y = 0 as a reference (although the choice of ⃗θ_⋆ is not specified in HG07). In order to be able to associate a set of parameters ⃗θ_p with a unique δY, the calibration method therefore makes the assumption that ∂ℳ_⋆/∂Y does not depend on the reference point (Y; ⃗θ_⋆). However, despite some arguments that are given in favour of a relative independence regarding the choice of Y, we show in Sect. 5 that this assumption about ⃗θ_⋆ is inadequate. These dependences cannot be reflected in model 𝒫, which depends on too few parameters, namely ⃗θ_p (we note that this is part of the reasons for introducing 𝒫 in the first place). To minimise the error that is introduced, it is then necessary to determine as many ⃗θ_p(δY) relations as stars that are studied (given independent constraints on ⃗θ_⋆, see Verma et al. 2019), which can be quite costly. Moreover, by using such a relation anyway, we reintroduce the biases we wished to avoid when studying the ionisation glitch. The relation that is obtained will largely depend on the physics considered in ℳ_⋆, in particular, on mechanisms that are not necessarily relevant for modelling the ionisation region.

To summarise, even though it is based on a more sophisticated model than previous methods, the method described above faces challenges that are mainly due to its empirical introduction. In particular, the lack of an explicit dependence on Y greatly reduces its applicability. We would like to reverse the approach. The idea is firstly, to introduce a physics-based parameterisation allowing inferences concerning ionisation regions, and secondly, to give a mathematical description regarding the Γ₁ profile. The challenge is to provide such a modelling while keeping the number of parameters as low as possible. This problem is addressed in the two following sections.

3. First adiabatic exponent in an ionisation region

In order to model the structure of the ionisation zone, we first try to provide an analytical expression of Γ₁ in this specific region. As previously mentioned, we avoid introducing it through similarities with mathematical functions, but rather derive it from thermodynamic relations. Expressions like this have already been obtained (e.g. Cox & Giuli 1968; Kippenhahn et al. 2012), but are generally functions of the occupation number rather than state variables (such as T and V). Furthermore, because the chemical equilibrium resulting from Saha’s equations is complex, these equations are generally solved numerically, and analytical expressions are limited to the ionisation of single species. We try to present in the following the simplest model that can nevertheless explicitly involve chemical composition.

3.1. Free energy

Our starting point is the free energy of a multi-component ideal gas, that is,

$$\begin{array}{c}\hfill F(T,V,N_{\alpha })={\displaystyle \sum _{\alpha }}{N}_{\alpha }({\mu }_{\alpha }-kT),\end{array}$$(14)

where V designates the volume, kT is the temperature expressed in energy units, and N_α and μ_α(T, V, N_α) are the number and chemical potential of particles of type α. In the context of partial ionisation, α refers to (i, r), 1 ≤ i ≤ N_sp being an index on the chemical species and 0 ≤ r ≤ z_i the ionisation state (where z_i is the atomic number), or it may also correspond to e for an electron. The key point of Eq. (14) lies in its ability to provide us with the pressure P and entropy S through its first derivative. Moreover, its second derivative gives access to virtually any other thermodynamic quantity ranging from the speed of sound through the adiabatic exponents including Γ₁ to various compressibilities. Before the calculations, we would like to define the particular meaning of a derivative in stellar physics compared with pure thermodynamics in more detail. The first adiabatic exponent is often defined as follows:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1\equiv {\left(\frac{\partial lnP}{\partial ln\rho }\right)}_S=-{\left(\frac{\partial lnP}{\partial lnV}\right)}_S,\end{array}$$(15)

where the subscript S denotes a partial derivative taken at constant entropy. The above expression can easily be rewritten using the more common variables T and V,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathrm{\Gamma }}_1=\frac{V}{P}{{\left(\frac{\partial S}{\partial T}\right)}_V}^{-1}[{\left(\frac{\partial P}{\partial T}\right)}_V{\left(\frac{\partial S}{\partial V}\right)}_T-{\left(\frac{\partial P}{\partial V}\right)}_T{\left(\frac{\partial S}{\partial T}\right)}_V]\end{array}.\end{array}$$(16)

Although this equation is generally valid, its use in the present case requires some clarification because it only mentions the two state variables T and V and does not explicitly show the dependence on the various number of particles N_α. Assuming that species are indexed so that z_i = i ($N_i={\sum }_r{N}_i^r$ potentially being null), the free energy introduced in Eq. (14) is a function of N_sp(N_sp + 3)/2 + 3 variables. As a result, the partial derivatives of Eqs. (15) and (16) are generally ambiguous unless some indication is given concerning the N_sp(N_sp + 3)/2 + 1 implicit conservation equations. Partial derivatives are commonly considered at constant values of state variables, for example,

$$\begin{array}{c}\hfill S=-{\left(\frac{\partial F}{\partial T}\right)}_{V,N_{\alpha }},P=-{\left(\frac{\partial F}{\partial V}\right)}_{T,N_{\alpha }},\forall \alpha {\mu }_{\alpha }={\left(\frac{\partial F}{\partial N_{\alpha }}\right)}_{T,V,N_{\beta \ne \alpha }}.\end{array}$$(17)

However, it is clear that the derivatives appearing in Eq. (16) (and second derivatives of F considering Eq. (17)) are not of this kind. As mentioned in the introductory section, the Γ₁ profile used in stellar physics shows clear signs of ionisation that derivatives taken at constant N_α do not. This choice imposes a particular ionisation state in the whole region (described by the constant N_α values), which instead should fluctuate widely with the thermodynamic conditions in the CZ. The solution is to consider the following conservation equations instead:

$$\begin{array}{cc}& \forall i,\forall r>0,{\mu }_i^r+{\mu }_e-{\mu }_i^{r-1}=0\left[\frac{N_{\mathrm{sp}}(N_{\mathrm{sp}}+1)}{2}\right]\hfill \end{array}$$(18)

$$\begin{array}{cc}& \forall i,{\displaystyle \sum _{r=0}^{z_i}}{N}_i^r=N_i=\mathrm{cnst}\left[N_{\mathrm{sp}}\right]\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\displaystyle \sum _{i=1}^{N_{\mathrm{sp}}}\sum _{r=1}^{z_i}r{N}_i^r=N_e\left[1\right].}\hfill \end{array}$$(20)

The first equation corresponds to the chemical equilibrium of the ionisation reaction, $A_i^{r-1}\rightleftharpoons A_i^r+e^{-}$. The second equation is merely the conservation of each atom number and can also be considered as a chemical equilibrium in the absence of reactions and homogeneous mixing of the CZ. The last expression corresponds to the overall charge equilibrium, that is, electroneutrality. This set of relations can be interpreted as a local equilibrium for given temperature and volume conditions that should be verified at each point of the CZ, and partial derivatives taken with respect to these constraints are hereafter be denoted with the subscript EQ. The resulting number of constraints is indicated in square brackets, and these constraints add up to N_sp(N_sp + 3)/2 + 1 as required.

Hereafter, we use the notation ${\partial }_{\alpha \beta }^{2}F$ to designate a second derivative taken with respect to β (subject to N_α being constant) and then α (subject to EQ). For example, Eq. (16) becomes

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1\equiv -{\left(\frac{\partial lnP}{\partial lnV}\right)}_{S,\mathrm{EQ}}=V\frac{({\partial }_{\mathit{VT}}^{2}F)({\partial }_{\mathit{TV}}^{2}F)-({\partial }_{\mathit{VV}}^{2}F)({\partial }_{\mathit{TT}}^{2}F)}{({\partial }_{V}F)({\partial }_{\mathit{TT}}^{2}F)}\end{array}$$(21)

by using relations (17). Considering Eq. (14), the evaluation of the Eq. (21) now requires finding $N_i^r(T,V)$ at a given EQ.

3.2. Approximate local equilibrium

To solve the system (18)–(20), the expression of the chemical potential of each particle must first be explicitly be written out. We have

$$\begin{array}{c}\hfill {\mu }_{\alpha }(T,V,N_{\alpha })=kTln\left(\frac{N_{\alpha }{\lambda }_{\alpha }^3}{{\mathcal{Z}}_{\alpha }V}\right),\end{array}$$(22)

with λ_α(T) and 𝒵_α(T) the thermal De Broglie wavelength and the canonical partition function of a particle α, respectively. Conditions (18) thus become

$$\begin{array}{c}\hfill \forall i,\forall r>0,\frac{x_i^r}{x_i^{r-1}}=\frac{{\mathcal{Z}}_i^r{\mathcal{Z}}_e}{{\mathcal{Z}}_i^{r-1}}\frac{V}{N_e{{\lambda }_e}^3},\end{array}$$(23)

with $x_i^r=N_i^r/N_i$ the occupation number of state (i,r). In our effort to obtain the simplest model that can be used to handle chemical composition, the ideal gas assumption remains relevant. We note, however, that in order to take the ionisation process in a more accurate way into account, it would be necessary to consider the Coulomb effects (Rogers 1981; Rogers et al. 1996). In this ideal case, ${\mathcal{Z}}_i^r$ and 𝒵_e are often approximated by (e.g., Kippenhahn et al. 2012)

$$\begin{array}{cc}\hfill {\mathcal{Z}}_i^r(T)& \simeq u_i^r(T)exp\left({\displaystyle {\displaystyle \sum _{s=1}^r}{\chi }_i^s/kT}\right)\simeq g_i^{r}exp\left({\displaystyle {\displaystyle \sum _{s=1}^r}{\chi }_i^s/kT}\right)\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\mathcal{Z}}_e& \simeq 2,\hfill \end{array}$$(25)

with $u_i^r$ describing the fine structure of state (i,r) approached by the degeneracy of the ground state $g_i^r$, and ${\chi }_i^r$ the ionisation energy separating state (i,r) from (i,r − 1) ($\sum _s^r{\chi }_i^s/kT$ is then the energy from the reference state (i,0), which is assumed here to have a null internal energy). The second approximation (25) directly results from neglecting the electron spin energy contribution.

Thus, conditions (23) simply become Saha’s equations,

$$\begin{array}{c}\hfill \forall i,\forall r>0,\frac{x_i^r}{x_i^{r-1}}=\frac{2g_i^r}{g_i^{r-1}}\frac{V}{N_e{{\lambda }_e}^3}{e}^{-{\chi }_i^r/kT}.\end{array}$$(26)

These equations are highly coupled, especially because of the term N_e, which can be expressed from (20) as

$$\begin{array}{c}\hfill N_e=N{\displaystyle \sum _{i=1}^{N_{\mathrm{sp}}}\sum _{r=1}^{z_i}r{x}_i{x}_i^r}\end{array}$$(27)

with N = ∑_iN_i and x_i = N_i/N the total number of atoms and the number abundance of element i, respectively, which are two constants from (19). Because of this, Saha’s equations are generally solved numerically by iterating over estimates of the $x_i^r$. In the following, we try to give analytical approximations by considering a few simplifying assumptions. First, we introduce the cumulative occupation number $y_i^r={\sum }_{s\ge r}{x}_i^s$, which is more suited to the study of Eqs. (26) (Baker & Kippenhahn 1962). We also use from this point the convention

$$\begin{array}{c}\hfill {\displaystyle \sum _{\mathit{ir}}\equiv \sum _{i=1}^{N_{\mathrm{sp}}}\sum _{r=1}^{z_i}}\end{array}$$

In this way, we can write $N_e=N\overline{z}$ with

$$\begin{array}{c}\hfill \overline{z}={\displaystyle \sum _{\mathit{ir}}{x}_i{y}_i^r,}\end{array}$$(28)

the mean number of electrons per atom. Saha’s equations can be rewritten as

$$\begin{array}{c}\hfill \forall i,\forall r>0,\frac{y_i^r-y_i^{r+1}}{y_i^{r-1}-y_i^r}=\frac{2g_i^r}{g_i^{r-1}}\frac{V}{N\overline{z}{{\lambda }_e}^3}{e}^{-{\chi }_i^r/kT},\end{array}$$(29)

with the convention $y_i^{z_i+1}=0$ (and by definition $y_i^0=1$). We start by noting that it is rare to find more than two ionisation states simultaneously for a given element. It is thus advantageous to make the following approximation (Baker & Kippenhahn 1962):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \forall i,\forall r>0,\mathrm{if}{y}_i^r\ne 0,\mathrm{then},\forall s<r,y_i^s& =1\hfill \\ \hfill \forall s>r,y_i^s& =0\hfill \end{array}\end{array}$$(30)

Then, in the single domain of interest ($y_i^r\ne 0$), (29) is approached by

$$\begin{array}{c}\hfill \forall i,\forall r>0,\frac{y_i^r}{1-y_i^r}=\frac{2g_i^r}{g_i^{r-1}}\frac{V}{N\overline{z}{{\lambda }_e}^3}{e}^{-{\chi }_i^r/kT}.\end{array}$$(31)

A final assumption is needed to express $\overline{z}$ as a function of $y_i^r$ only and to completely decouple the system. One idea may be to exploit the fact that stars are predominantly composed of hydrogen in number: For typical Y values (i.e. 1/4 < Y < 1/3), x₁ ∼ 0.89 − 0.92, and it can be seen using (28) that after a complete ionisation of the hydrogen, $\overline{z}>x_1$. Another bound can then be obtained on $\overline{z}$ because the latter tends asymptotically towards ∑_iz_ix_i = 2 − x₁ + ∑_i > 2(z_i − 2)x_i. Because for heavy elements (i > 2), m_i ≃ 2z_im_u with m_i the mass of atom i and m_u the atomic mass unit, and because the mean mass m₀ < 2m_u in a stellar mixture, we have m₀/m_i < 1/z_i. Finally, the number and mass abundances are related by x_i = m₀X_i/m_i, meaning that $\overline{z}$ should remain below the upper boundary 2 − x₁ + Z ∼ 1.1 with Z = ∑_i > 2X_i the metal abundance in mass. This simple reasoning allows us to convince ourselves that $\overline{z}$ must remain close to 1 when the ionisation of hydrogen is completed, and in a way relatively independent of the mixture under consideration.

The pragmatic assumption can therefore be made that the ionisation of hydrogen occurs first and is completed before any other ionisation begins. Thus we can use the fact that hydrogen is the dominant component to approximate $\overline{z}$ by the mean electron number in a pure-hydrogen model. Although this is not entirely true (particularly not for a low ionisation state), this assumption greatly simplifies the problem. For all elements i > 1, $\overline{z}\sim 1$, which, as discussed above, is a reasonable approximation. For the ionisation equation corresponding to hydrogen (i = 1), we instead approximate $\overline{z}\sim {\mathit{y}}_1^1$. Our version of Saha’s equations can thus finally be written as

$$\begin{array}{c}\hfill \forall i,\forall r>0,\frac{{(y_i^r)}^n}{1-y_i^r}=\frac{2g_i^r}{g_i^{r-1}}\frac{V}{N{{\lambda }_e}^3}{e}^{-{\chi }_i^r/kT},\end{array}$$(32)

with n = 2 if i = 1 and n = 1 in any other case. We define $K_i^r(T,V)$ as the term on the right-hand side (RHS) of (32). We then obtain the desired relations.

$$\begin{array}{cc}\hfill y_1^1(T,V)& =\frac{1}{2}[\sqrt{K_1^1(K_1^1+4)}-K_1^1]\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill \forall i>1,\forall r>0,y_i^r(T,V)& =\frac{K_i^r}{1+K_i^r}\hfill \end{array}$$(34)

It is easy to retrieve $x_i^r(V,T)$ from there using $x_i^r = y_i^r-y_i^{r+1}$ (as well as $N_i^r(V,T)$ using $N_i^r=N_i{x}_i^r$). Figure 4 shows a comparison of hydrogen and helium ionisation fractions $x_1^1,x_2^1$, and $x_2^2$ obtained analytically from (33) and (34) for given T and V values with the numerical solution of Saha’s Eqs. (26). The result is shown as a function of the normalised acoustic depth t so that it can be reduced to a comparison of profiles rather than a (V, T) mapping. The solution presented above appears to be satisfying considering the number of assumptions made. However, $x_2^1$ is slightly underestimated because $\overline{z}$ is slightly lower than 1 in the first helium ionisation region. The opposite reasoning holds for the second helium ionisation region.

Fig. 4. Comparison of approximations (33) and (34) in the case i = 1, 2 (opaque curves) with the numerical solution of (26) (dashed light curves) for solar near-surface conditions of (T, V, Xi > 1). These quantities are represented as a function of the normalised acoustic depth t = τ/τ0 (the surface is then located at t = 0), along which all future plots are represented.

3.3. First adiabatic exponent

Approximations (33) and (34) allows us to derive all derivatives appearing in (21). First of all, we have

$$\begin{array}{c}\hfill {\partial }_{V}F=-P=\frac{-kT}{V}{\displaystyle \sum _{\alpha }{N}_{\alpha }=\frac{-NkT}{V}(1+\sum _{\mathit{ir}}{x}_i{y}_i^r),}\end{array}$$(35)

where we explicitly made the $y_i^r$ apparent to fully exploit relations (33) and (34). It is straightforward to show that we simply retrieve the classical expression for the pressure in an ideal gas. Second derivatives taken at EQ ensue,

$$\begin{array}{cc}\hfill {\partial }_{\mathit{VV}}^{2}F& =\frac{\mathit{NkT}}{V^2}[1+{\displaystyle \sum _{\mathit{ir}}{x}_i{y}_i^r(1-\frac{\partial ln{y}_i^r}{\partial lnV})}]\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill {\partial }_{\mathit{TV}}^{2}F& =\frac{-Nk}{V}[1+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(1+\frac{\partial ln{y}_i^r}{\partial lnT})].\hfill \end{array}$$(37)

More explicit expressions can be deduced by writing ${\displaystyle \frac{\partial ln{\mathit{y}}_i^r}{\partial ln\alpha }=\frac{\mathrm{d}ln{\mathit{y}}_i^r}{\mathrm{d}ln{K}_i^r}\frac{\partial ln{K}_i^r}{\partial ln\alpha }}$ and exploiting Eq. (32),

$$\begin{array}{cc}\hfill \forall i,\forall r>0,\frac{\mathrm{d}ln{y}_i^r}{\mathrm{d}ln{K}_i^r}& =\frac{1-y_i^r}{1+{\delta }_i^1(1-y_i^r)}\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill \frac{\partial ln{K}_i^r}{\partial lnV}& =1\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill \frac{\partial ln{K}_i^r}{\partial lnT}& =\frac{3}{2}+\frac{{\chi }_i^r}{\mathit{kT}}\equiv {\varphi }_i^r\hfill \end{array}$$(40)

${\delta }_i^j$ standing for Kronecker’s symbol. Resulting expressions of the derivatives are

$$\begin{array}{cc}\hfill {\partial }_{\mathit{VV}}^{2}F& =\frac{\mathit{NkT}}{V^2}[1+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(1-\frac{1-y_i^r}{1+{\delta }_i^1(1-y_i^r)})]\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill {\partial }_{\mathit{TV}}^{2}F& =\frac{-Nk}{V}[1+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(1+\frac{(1-y_i^r){\varphi }_i^r}{1+{\delta }_i^1(1-y_i^r)})].\hfill \end{array}$$(42)

The calculation of the two remaining derivatives is more subtle, however. We propose here to start with a determination of the energy E,

$$\begin{array}{c}\hfill E=-T^2{\left(\frac{\partial }{\partial T}\left(\frac{F}{T}\right)\right)}_{V,N_{\alpha }}=NkT(\frac{3}{2}+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r{\varphi }_i^r).\end{array}$$(43)

Because ${\partial }_{VT}^{2}F=-{(\partial S/\partial V)}_{T,\text{EQ}}$ and ${\partial }_{TT}^{2}F=-{(\partial S/\partial T)}_{V,\text{EQ}}$, it might be tempting to take directly the derivatives of (43) using dE = TdS − PdV. However, in the present framework,

$$\begin{array}{cc}\hfill \mathrm{d}E& =T\mathrm{d}S-P\mathrm{d}V+{\displaystyle \sum _{\alpha }{\mu }_{\alpha }\mathrm{d}{N}_{\alpha }}\hfill \\ \hfill & =T\mathrm{d}S-P\mathrm{d}V+N{\displaystyle \sum _{\mathit{ir}}{x}_i({\mu }_i^r+r{\mu }_e)\mathrm{d}{x}_i^r+N\sum _i{x}_i{\mu }_i^0\mathrm{d}{x}_i^0}\hfill \\ \hfill & =T\mathrm{d}S-P\mathrm{d}V+N{\displaystyle \sum _{\mathit{ir}}{x}_i({\mu }_i^r+{\mu }_e-{\mu }_i^{r-1})\mathrm{d}{y}_i^r}\hfill \\ \hfill & +N{\displaystyle \sum _{\mathit{ir}}{x}_i{\mu }_i^0\mathrm{d}{x}_i^r+N\sum _i{x}_i{\mu }_i^0\mathrm{d}{x}_i^0.}\hfill \end{array}$$(44)

Clearly, the last two sums cancel each other because $\forall i, {\sum }_{r>0}\text{d}{x}_i^r=d(1-x_i^0)=-\text{d}{x}_i^0$, which is verified in our derivation as long as we define $x_i^0\equiv 1-y_i^1$. Electroneutrality, used in line 2 of Eq. (44), is also verified by imposing $N_e\equiv N{\sum }_{ir}{x}_i{y}_i^r$. According to Eq. (18), ${\mu }_i^r+{\mu }_e-{\mu }_i^{r-1}=0$, and the third term must also cancel exactly. In practice, Eq. (18) is too complex to be perfectly solved, and numerical calculations as well as analytical approximations can lead to some small departures from exact cancellation. These may lead to slight thermodynamic inconsistencies (which could be corrected through departures from equality of partial mixed second derivatives of state functions such as the free energy F). However, these residuals remain small, and we considered that chemical equilibrium is perfectly satisfied. The only part that therefore remains from (44) is the following well-known identity:

$$\begin{array}{c}\hfill \mathrm{d}E=T\mathrm{d}S-P\mathrm{d}V.\end{array}$$(45)

This means that ${\partial }_{VT}^{2}F=-{(\partial E/\partial V)}_{T,\text{EQ}}-P$ and $T{\partial }_{TT}^{2}F = -{(\partial E/\partial T)}_{V,\text{EQ}}$. It follows from Eq. (43) that

$$\begin{array}{cc}\hfill {\partial }_{\mathit{VT}}^{2}F& =\frac{-Nk}{V}[1+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(1+\frac{(1-y_i^r){\varphi }_i^r}{1+{\delta }_i^1(1-y_i^r)})]\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill {\partial }_{\mathit{TT}}^{2}F& =\frac{-Nk}{T}[\frac{3}{2}+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(\frac{3}{2}+\frac{(1-y_i^r){({\varphi }_i^r)}^2}{1+{\delta }_i^1(1-y_i^r)})].\hfill \end{array}$$(47)

Denoting the bracket part of ${\partial }_{\alpha \beta }^{2}F$ as ${\partial }_{\alpha \beta }^{2}f$, we obtain

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1=\frac{({\partial }_{\mathit{VT}}^{2}f)({\partial }_{\mathit{TV}}^{2}f)+({\partial }_{\mathit{VV}}^{2}f)({\partial }_{\mathit{TT}}^{2}f)}{({\partial }_{V}f)({\partial }_{\mathit{TT}}^{2}f)},\end{array}$$(48)

where

$$\begin{array}{cc}\hfill {\partial }_{V}f& =1+{\displaystyle \sum _{\mathit{ir}}{x}_i{y}_i^r,}\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill {\partial }_{\mathit{VV}}^{2}f& =1+{\displaystyle \sum _{\mathit{ir}}{x}_i{y}_i^r(1-\frac{1-y_i^r}{1+{\delta }_i^1(1-y_i^r)}),}\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill {\partial }_{\mathit{TV}}^{2}f={\partial }_{\mathit{VT}}^{2}f& =1+{\displaystyle \sum _{\mathit{ir}}{x}_i{y}_i^r(1+\frac{(1-y_i^r){\varphi }_i^r}{1+{\delta }_i^1(1-y_i^r)}),}\hfill \end{array}$$(51)

$$\begin{array}{cc}\hfill {\partial }_{\mathit{TT}}^{2}f& =\frac{3}{2}+{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(\frac{3}{2}+\frac{(1-y_i^r){({\varphi }_i^r)}^2}{1+{\delta }_i^1(1-y_i^r)}).\hfill \end{array}$$(52)

The various symmetries present in these equations may give the impression that expression (48) can be further simplified. We show in Appendix A that Γ₁ can finally be written as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_1=\frac{5}{3}-\frac{2}{3}{\gamma }_1,\end{array}$$(53)

with 0 ≤ γ₁ < 1 simply expressed as

$$\begin{array}{c}\hfill {\gamma }_1=\frac{1}{{\partial }_{\mathit{TT}}^{2}f}{\displaystyle \sum _{\mathit{ir}}}{x}_i{y}_i^r(1-y_i^r)\frac{{({\chi }_i^r/kT)}^2}{1+{\delta }_i^1(1-y_i^r)}.\end{array}$$(54)

Together with Eqs. (32)–(34), Eqs. (52)–(54) provide an analytical approximation of Γ₁(T, V, x_i) in the CZ. Even though this expression is applicable to any kind of mixture (but is assumed to be homogeneous in the CZ!), the following examples are based on a hydrogen-helium mixture of mass fraction 1 − Y and Y (Y = 4x₂/(1 + 3x₂)) to facilitate the study. We recall that N/V = ρ/m₀, and represent the resulting Γ₁(ρ, T, Y) in Fig. 5 for Y = 0.25 as well as for extreme values, Y = 0 and Y = 1. Each panel also shows the relation ρ(T) extracted from a CESTAM solar model (Morel & Lebreton 2008; Marques et al. 2013). The latter is obviously only indicative as it depends on the composition, but it still provides an idea of what part of the map is visible in a Γ₁ profile. The figure clearly reflects the Γ₁ depressions caused by (from left to right) the hydrogen and first and second helium ionisations. Each element contribution is enhanced in the left panels. In the top right panel, the superposition of the hydrogen and helium first ionisation zones in typical solar conditions is presented, whereas the distinction becomes apparent at lower densities. In contrast, this signature appears to be more diffuse at the highest densities. This behaviour appears in each panel and is further discussed in Sect. 4. The last frame shows the variation δ_ρ, TΓ₁ caused by a change in helium abundance, δY = 0.1, from a reference value Y = 0.25. Although the variables are not normalised, the notation δ_ρ, T is still meaningful because Γ₁ does not refer to a profile here, but to a function defined for all (ρ, T) at a given Y. As expected, this variation in abundance causes the two helium wells to become lower while elevating the hydrogen well. However, the variations appear to be fairly disproportionate between the two elements. The change in the hydrogen ionisation region would be barely noticeable considering the typical solar relation ρ(T), shown in black. Appendix A provides clues that help to understand this difference in behaviour as well as the particular shape of the perturbation in the hydrogen region compared with that of helium. In particular, the amplitude of the variation seems more related to the relative change δx_i/x_i than δY. While we have δx₁ = −δx₂ for the hydrogen-helium case, we see that δx₁/x₁ = −δx₂/x₁ = −(x₂/x₁)δx₂/x₂ with (x₂/x₁) = 1/12 for Y = 0.25. This explains why the change in hydrogen is expected to cause a variation that is an order of magnitude lower than the change in helium.

Fig. 5. 2D maps of the first adiabatic exponent obtained from Eq. (53). Panels a–c: Γ1(ρ, T, Y) as a function of ρ and T for different values of Y. Panel d: variation δρ, TΓ1 caused by a perturbation δY = 0.1 from a reference value Y = 0.25. The dashed line in each panel shows the relation ρ(T) extracted from a CESTAM solar model.

Although Eqs. (53) and (54) appear to be promising, they are merely functional forms and do not allow us at this point to model a parameterised structure of the ionisation region. However, Γ₁, ρ and T are closely related in stellar interiors.

4. Model structure and properties

In this section we obtain the stellar structure associated with the Γ₁ expression derived in the previous section (see Eqs. (53) and (54)). To do this, we assume the conditions of a CZ that correspond to the assumptions for deriving the first adiabatic exponent. These are (1) an isentropic region, (2) chemical equilibrium at any position, (3) uniform abundances (due to homogeneous mixing), and (4) electroneutrality at any position. The last three conditions correspond to applying EQ conditions at any point in the region. Moreover, to avoid getting lost in the many possibilities offered by these relations in terms of the mixture, we again consider a relation of the form Γ₁(ρ, T, Y) in this part. It is possible o generalise all the observations made here to any (reasonable) mixture by replacing Y by X_i > 1, however.

4.1. Ionisation region structure

As mentioned above, the quantities Γ₁, ρ and T are not independent in stellar interiors, and although it may be interesting to study Γ₁ for any condition of temperature and density, relevant profiles in these environments are much more constrained. The appropriate equations to express these constraints are the well-known hydrostatic equilibrium and Poisson’s equation assuming spherical symmetry,

$$\begin{array}{cc}\hfill \frac{\mathrm{d}P}{\mathrm{d}r}& =-\rho g\hfill \end{array}$$(55)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}g}{\mathrm{d}r}+\frac{2g}{r}& =4\pi G\rho ,\hfill \end{array}$$(56)

where we have introduced the additional intermediate variable g, the norm of the gravity field. At this point, it is important to recall that we try to account for the variation of T and ρ in the region rather than that of P. This can be achieved by noting that for any quantity α, ${\displaystyle \frac{\mathrm{d}\alpha }{\mathrm{d}r}=\frac{-\rho g\alpha }{P}{\left(\frac{\mathrm{d}lnP}{\mathrm{d}ln\alpha }\right)}^{-1}}$. Then, crucially, we note that derivatives taken with respect to the radius r such as dlnP/dlnα = (dlnP/dr) (dlnα/dr)⁻¹ exactly correspond to derivatives taken at constant entropy and EQ (cf. assumptions 1–4 at the beginning of the section). Using this for α = T, ρ, Eqs. (55) and (56) can be changed into a differential system of three equations, the solution of which is (T, ρ, g)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}T}{\mathrm{d}r}& =-\frac{\rho Tg}{P}\frac{{\mathrm{\Gamma }}_2-1}{{\mathrm{\Gamma }}_2}\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\rho }{\mathrm{d}r}& =-\frac{{\rho }^{2}g}{P}\frac{1}{{\mathrm{\Gamma }}_1}\hfill \end{array}$$(58)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}g}{\mathrm{d}r}& =4\pi G\rho -\frac{2g}{r},\hfill \end{array}$$(59)

where we introduced the second adiabatic exponent Γ₂ such that ${\displaystyle \frac{{\mathrm{\Gamma }}_2}{{\mathrm{\Gamma }}_2-1}\equiv {\left(\frac{\partial lnP}{\partial lnT}\right)}_{S,\mathrm{EQ}}}$. Its expression can easily be deduced from Γ₁ by noting that their definitions are symmetrical with respect to a change V ↔ T,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{{\mathrm{\Gamma }}_2}{{\mathrm{\Gamma }}_2-1}& =\frac{T}{P}{{\left(\frac{\partial S}{\partial V}\right)}_T}^{-1}[{\left(\frac{\partial P}{\partial T}\right)}_V{\left(\frac{\partial S}{\partial V}\right)}_T-{\left(\frac{\partial P}{\partial V}\right)}_T{\left(\frac{\partial S}{\partial T}\right)}_V]\hfill \end{array}.\end{array}$$(60)

An analogous expression to Eq. (53) can be provided for Γ₂,

$$\begin{array}{c}\hfill \frac{{\mathrm{\Gamma }}_2}{{\mathrm{\Gamma }}_2-1}=\frac{5}{2}+{\gamma }_2,\end{array}$$(61)

by introducing

$$\begin{array}{c}\hfill {\gamma }_2=\frac{1}{{\partial }_{\mathit{VT}}^{2}f}{\displaystyle \sum _{\mathit{ir}}{x}_i{y}_i^r(1-y_i^r)\frac{({\chi }_i^r/kT)(5/2+{\chi }_i^r/kT)}{1+{\delta }_i^1(1-y_i^r)}.}\end{array}$$(62)

The pressure can also be removed using the ideal gas law (35), which leads to the following system:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}T}{\mathrm{d}r}& =-\frac{m_{0}g}{k}\frac{1}{{\partial }_{V}f}\frac{{\mathrm{\Gamma }}_2-1}{{\mathrm{\Gamma }}_2}\hfill \end{array}$$(63)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\rho }{\mathrm{d}r}& =-\frac{\rho m_{0}g}{\mathit{kT}}\frac{1}{{\partial }_{V}f}\frac{1}{{\mathrm{\Gamma }}_1}\hfill \end{array}$$(64)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}g}{\mathrm{d}r}& =4\pi G\rho -\frac{2g}{r}.\hfill \end{array}$$(65)

All quantities on the RHS are functions of ρ, T, and g only (for a given value of Y). By imposing the central values T_c and ρ_c (g_c must be null in any case), we now have a complete system, the solution of which is ${\displaystyle \left\{\begin{array}{c}T(r;Y,T_c,{\rho }_c)\\ \rho (r;Y,T_c,{\rho }_c)\\ g(r;Y,T_c,{\rho }_c)\end{array}\right\}}$. The resulting model is fully convective and contains an ionisation region whose properties are controlled by the triplet (Y, T_c, ρ_c). Although this triplet seems the most intuitive, it is possible to parameterise the model by considering other triplets derived from these three quantities. Hereafter, we instead consider the following dimensionless parameterisation:

$$\begin{array}{cc}\hfill Y_s& \equiv Y,\hfill \end{array}$$(66)

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{CZ}}& \equiv ln\left[\frac{{\rho }_c{h}^3}{2m_u{(2\pi m_{e}k{T}_c)}^{3/2}}\right],\hfill \end{array}$$(67)

$$\begin{array}{cc}\hfill \epsilon & \equiv \frac{{\chi }_{\mathrm{H}}}{k{T}_c}\hfill \end{array}$$(68)

where we introduced the Planck constant h, the electron mass m_e, and the ionisation potential of hydrogen χ_H. This particular choice and the associated notations will become clearer later on. At this point, this choice entirely constrains the Γ₁ value at the centre, the latter being only a function of Y, $2m_u/\rho {\lambda }_e^3$ (we recall that ${\lambda }_e^{-2}=2\pi m_{e}kT/h^2$), and χ_H/kT (cf. Eqs. (32)–(34) and Eqs. (52)–(54)). However, their impact on our model is much more significant than this.

Physical interpretation of the parameters. Because Y(r) is constant (see assumption 3), the helium amount Y_s is clearly imposed at any location of the model up to the surface symbolised by the index s.

The second parameter derives its notation from its close link to the electron degeneracy parameter $\psi ={\mu }_e/kT=\ln (n_e{\lambda }_e^3/2)$ in the ideal gas limit. We can show that

$$\begin{array}{c}\hfill \psi \underset{\mathrm{at}\mathrm{centre}}{=}{\psi }_{\mathrm{CZ}}-\frac{Y}{2}+\mathcal{O}(Y^2).\end{array}$$(69)

However, as long as Γ₁ = 5/3 (i.e., outside the ionisation region), ρ ∝ T^3/2, which means that ψ is preserved. Therefore ψ_CZ approximates the electron degeneracy parameter from the centre up to the ionisation region.

The last parameter ε results from the ratio of two energy scales in the model. The first is the ionisation energy for the hydrogen (13.6 eV), and the second is the order of magnitude of the thermal energy at the centre of the model (∼keV), making this ratio generally rather small (hence the notation ε). Through conservations, we just saw that Y_s and ψ_CZ reflect properties of the ionisation region; this also holds for its helium abundance and electron degeneracy. It is complicated to apply this to ε because χ_H/kT varies by several orders of magnitude in the model. We try to give an intuitive understanding of its impact on the model structure, however. We first note that if Γ₁ did not depend on χ_H/kT, the latter would be equal to 5/3 in the whole model. Because Γ₁ would depend only on Y (which is constant) and $2m_u/\rho {\lambda }_e^3$, the value of which is fixed as long as ρ ∝ T^3/2, the first adiabatic exponent would keep its central value, that is, 5/3. The value of χ_H/kT at a given position then alone determines whether Γ₁ deviates from the value 5/3, that is, if the ionisation region is entered starting from the model centre. It can be shown that the ionisation of the element i in state r takes place when ${\chi }_i^r/kT~-\psi$. The role of ε can now be clearly identified: When ε is to be large (e.g., of the order of −ψ), the ionisation will be complete close to the centre. In contrast, if ε is small before −ψ, it will take until T ≪ T_c for the ionisation to be completed and the region will be shifted closer to the surface. This shows that although it is not explicitly linked to any quantity of the ionisation zone, ε controls the extent of the region.

After this physical interpretation for the parameters, we now try to provide a geometrical interpretation of the impact of each parameter regarding the Γ₁ profile. We rely on Fig. 6, which shows the Γ₁(t) (with t = τ/τ₀, the normalised acoustic depth) profiles resulting from multiple parameter sets (Y_s, ψ_CZ, ε) by varying them one by one around a reference value.

Fig. 6. Γ1(t) profiles obtained by considering various sets (Ys, ψCZ, ε) around a reference $(Y_s^{\odot },{\psi }_{\text{CZ}}{}^{\odot },{\epsilon }^{\odot })=(0.255,-5,0.015)$ that corresponds to a Γ1 profile approaching that of the Sun (see Sect. 4.2). (a) Dependence of the profile on 0 ≤ Ys ≤ 1. (b) Dependence of the profile on −12 ≤ ψCZ ≤ −2. (c) Dependence of the profile on 0.01 ≤ ε ≤ 0.05.

Geometrical interpretation of the parameters regarding the Γ₁ profile. We gave a representation of Γ₁(t) for all possible values of Y_s in panel a of Fig. 6, including such extreme cases as pure-hydrogen or helium stars. As expected, the profile changes from one that only shows the hydrogen well (Y_s = 0) to one that only shows the two helium wells (Y_s = 1), passing through profiles that account for the three components (0 < Y_s < 1). Clearly, increasing Y_s has the effect of extending both helium wells while reducing the hydrogen well.

The effects of modifying ψ_CZ are shown in panel b of Fig. 6. The range under consideration corresponds to typical solar degeneracy values ranging from the surface (−12: very low degeneracy) to the centre (−2: high degeneracy). The impact is again clearly visible: A high degeneracy tends to spread out the ionisation wells, making them indistinguishable when ψ_CZ → −2. In the opposite situation, low degeneracies in the ionisation region result in much deeper and localised wells. It is even possible to distinguish the contributions of the hydrogen and the first helium ionisations when ψ_CZ approaches −12. This last point is consistent with what we showed in Fig. 5 : high degeneracies correspond to relations of the form lnρ(T) = C + (3/2)lnT with high C values. Thus, structures with high ψ_CZ correspond to ρ(T) that are located in the upper part of Fig. 5, where the ionisation regions tend to overlap. The opposite reasoning holds just as well.

The impact of ε is illustrated in panel c of Fig. 6. As expected from the previous paragraph, varying its value from 0.01 to 0.05 results in a shift of the ionisation region. However, this change differs between each well, the HeII well is more strongly affected than the hydrogen well. In order of magnitude, this effect can be explained by the dependence of the normalised acoustic depth on temperature close to the surface: $t\propto \sqrt{T/T_c}\propto \sqrt{\epsilon T}$. Thus, for a given ionisation temperature T_ion, the impact of ε becomes more noticeable on the associated depth t_ion as T_ion is high. This reasoning also explains why the shift does not appear to be linear with respect to ε for a given well (this is particularly visible for the HeII well); this shift varies as $\sqrt{\epsilon }$.

Figure 6 shows an interesting point of the model defined with (Y_s, ψ_CZ, ε): the relative independence of the impact of each parameter. When each well is treated as a distribution, all parameters appear to impact a specific moment. For instance, Y_s appears to be directly related to the distribution area and does not impact the shapes of the wells or the position much (for reasonable values of Y_s). ψ_CZ is clearly related to the dispersion of the distribution and ε to its position. These distinct impacts on the structure are encouraging for subsequent parameter retrieval as we can expect to avoid degeneracy biases.

4.2. Comparison with a realistic stellar model

The first and most important point to address is to verify the extent to which the model described in this article can approximate the ionisation region structure of a realistic stellar model. To do this, we compare in panel a of Fig. 7Γ₁(t) profiles extracted from a CESTAM solar model and from our model for a particular parameter set $(Y_s^{\odot },{\psi }_{\text{CZ}}^{\odot },{\epsilon }^{\odot })$. These values were chosen so that

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^1{\left(\frac{{\delta }_t{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}\right)}^2\mathrm{d}t={\int }_0^1\frac{{({\mathrm{\Gamma }}_1^{\mathrm{CESTAM}}(t)-{\mathrm{\Gamma }}_1(t;Y_s,{\psi }_{\mathrm{CZ}},\epsilon ))}^2}{{\mathrm{\Gamma }}_1^2}\mathrm{d}t}\end{array}$$(70)

Fig. 7. Correspondence between the CESTAM profiles and their respective analogues in our model. (a) Comparison of ${\mathrm{\Gamma }}_1^{ \text{CESTAM}}(t)$ extracted from a CESTAM solar model (in grey) and ${\mathrm{\Gamma }}_1(t;Y_s^{\odot },{\psi }_{\text{CZ}}{}^{\odot },{\epsilon }^{\odot })$ obtained from our model (in red). (b) Comparison of YCESTAM(t) extracted from the same model with $Y_s^{\odot }$. (c) Comparison of ψCESTAM(t) with ${\psi }_{\text{CZ}}^{\odot }$. In the last two panels, the regions in which the curves are expected to correspond are shaded in red.

is minimal. The behaviour of the realistic model can clearly be very well reproduced, the curves differ mainly in only two regions. A first difference is located in the first helium ionisation region (0.12 < t < 0.18), where the depth of the well is slightly overestimated. This is mainly caused by the assumption made in Sect. 3.2, where we assumed that $\overline{z}\sim 1$ in the helium ionisation region and thus underestimated the number of electrons in the HeI region. This causes the ionisation to start slightly deeper in the model and therefore to be more localised. By reducing the dispersion of the HeI well, the assumption slightly enlarged its depth. As expected, another difference is visible in the atmospheric region (t < 0.03). Because the ionisation has not yet started, our Γ₁ still has the value 5/3, which is obviously not the case in a realistic model containing an atmospheric modelling. The overall profile shape is nevertheless well reproduced in the ionisation region, especially considering that only three parameters need to be adjusted.

A second point to be clarified is the exact meaning of the adjusted values $(Y_s^{\odot },{\psi }_{\text{CZ}}^{\odot },{\epsilon }^{\odot })$ regarding the realistic model. When the minimal differences (δ_tΓ₁/Γ₁)² are assumed, the model parameterised by $(Y_s^{\odot },{\psi }_{\text{CZ}}^{\odot },{\epsilon }^{\odot })$ is interpreted as being the closest to the CESTAM model from the perspective of a particular structural aspect. This aspect will be discussed in much more detail in a forthcoming paper, but we can already relate it to the “visible” part of the structure through the frequency shift considering Eq. (5). In this sense, if our modelling and its interpretation are correct, we can expect our parameters to approach the quantities they refer to in the realistic model. However, in previous paragraphs we gave an understanding of these parameters in our model, but a realistic model only locally verifies the assumptions used to derive our structure. In particular, the relation between P and ρ is described by ${\displaystyle \gamma \equiv \frac{\mathrm{d}lnP}{\mathrm{d}ln\rho }}$ rather than Γ₁, the two differing outside the isentropic region. Accordingly, we expect the adjusted values to approximate the CESTAM model quantities where assumptions 1 to 4 are verified, that is, in the CZ. For the helium abundance Y_s, because Y_s = Y(t) at any position of our model, the value $Y_s^{\odot }$ should approximate Y^CESTAM(t) in the whole CZ up to the atmosphere. In the case of ψ_CZ, we saw that the approximation remains valid as long as Γ₁ = 5/3, therefore we expect ${\psi }_{\text{CZ}}^{\odot }$ to approach ψ^CESTAM(t) in the CZ part that is located below the ionisation region. As mentioned, because ε does not relate to any particular quantity of the CZ, it is difficult to deduce more than the position of the ionisation region. These approximations can be summarised as follows:

$$\begin{array}{cc}\hfill Y_s^{\odot }& \simeq Y^{\mathrm{CESTAM}}(t)\mathrm{for}\mathrm{BCZ}<t<\mathrm{ATM}\hfill \end{array}$$(71)

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{CZ}}^{\odot }& \simeq {\psi }^{\mathrm{CESTAM}}(t)\mathrm{for}\mathrm{BCZ}<t<\mathrm{ION},\hfill \end{array}$$(72)

with BCZ, ATM, and ION designating the base of the convective zone, the base of the atmosphere, and the end of the ionisation region, respectively. We illustrate this point in panels b and c of Fig. 7 by representing both Y^CESTAM(t) and ψ^CESTAM(t), as well as the estimates $Y_s^{\odot }$ and ${\psi }_{\text{CZ}}^{\odot }$ obtained through the minimisation of Eq. (70). We also highlight the regions in which the quantities are expected to correspond to each other based on Eqs. (71) and (72). The values thus obtained $(Y_s^{\odot },{\psi }_{\text{CZ}}^{\odot },{\epsilon }^{\odot })$ = (0.255, −5, 0.015) are consistent with the CESTAM quantities in the associated regions (see panels b and c of Fig. 7 and Table 1).

It is risky to be quantitative about this result because there are no obvious values with which the difference can be compared; in particular the relative difference does not appear to be relevant. Nevertheless, this example strongly suggests that the modelling defined here allows us to define structures that are close to realistic models (in the sense of criteria (70)) and to also obtain similar helium abundances and electron degeneracies in the regions of interest.

5. Analysis of first adiabatic exponent perturbations

So far, we focused on introducing a physical model of the ionisation region that depends on only a few parameters in order to study its properties. However, Eq. (5) shows that the analysis of frequency shifts relies more on the modelling of structural perturbations than on the structure itself. We therefore investigate what our model predicts as a perturbation caused by a change in surface helium abundance compared with the ad hoc profiles used in previous papers (see Fig. 2). The advantage of having modelled the structure is that we can easily reduce it to the study of a perturbation of helium abundance by considering the following difference between profiles:

$$\begin{array}{c}\hfill \frac{{\delta }_t{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}=\frac{{\mathrm{\Gamma }}_1(t;Y_s+\delta Y_s,{\psi }_{\mathrm{CZ}},\epsilon )-{\mathrm{\Gamma }}_1(t;Y_s,{\psi }_{\mathrm{CZ}},\epsilon )}{{\mathrm{\Gamma }}_1}.\end{array}$$(73)

The analysis of this perturbation is carried out in Fig. 8. An important aspect to bear in mind is the number of variables on which the above function depends. The perturbation naturally depends on the helium difference δY_s but also on the point (Y_s, ψ_CZ, ε) around which the differences are calculated, thus leading to four dimensions to explore. Our analysis of the 3D space required for the structure was obviously already incomplete in Fig. 6, and it follows that the perturbation analysis is necessarily even more superficial. Figure 8 provides plots for various reference points, (Y_s, ψ_CZ, ε), but for a fixed difference δY_s = 0.05. Nevertheless, it is straightforward to see by linearising Eq. (73) that a difference caused by $\delta Y_s^{\prime }=\alpha \times \delta Y_s$ would be relatively similar to the profile α × δ_tΓ₁/Γ₁. In this sense, Fig. 8 still provides a good representation of what the change in helium abundance might cause. Panel a shows that all the perturbed profiles roughly overlap, which reflects a relative independence of the helium amount chosen as a reference (the range considered, 0.2 < Y_s < 0.35, is representative of most realistic stellar helium abundances). The shape is quite similar to what is expected from a helium variation: Each helium ionisation region contributes to a Gaussian-like well (t ∼ 0.13 & t ∼ 0.19), and the second well is more pronounced. The hydrogen contribution is less intuitive, resulting in a very localised peak (t ∼ 0.03, present in every panel) at the beginning of ionisation. This peak reflects more of a shift of the hydrogen ionisation region (visible in panel a of Fig. 6) caused by the helium change than a real drop in the well. Nevertheless, this component is clearly noticeable in terms of amplitude.

Fig. 8. Profile differences δtΓ1/Γ1 = [Γ1(t;Ys+δYs,ψCZ,ε)−Γ1(t;Ys,ψCZ,ε)]/Γ1 obtained for a fixed value of δYs = 0.05 and various sets (Ys, ψCZ, ε) (the values used here are not the same as those in Fig. 6) around the reference $(Y_s^{\odot },{\psi }_{\text{CZ}}^{\odot },{\epsilon }^{\odot })$ = (0.255, −5, 0.015) as a function of the normalised acoustic depth t. (a) Dependence on 0.2 ≤ Ys ≤ 0.35. (b) Dependence on −9 ≤ ψCZ ≤ −3. (c) Dependence on 0.01 ≤ ε ≤ 0.05.

The situation in panel b is less intuitive. In order to provide realistic values for this rather uncommon parameter, we considered −9 < ψ_CZ < −3 as the electron degeneracy range, which corresponds to typical values inside solar-like oscillators (a justification of this point is provided later on). This time, the shape of the perturbation caused by a fixed change in the helium amount is very sensitive to the electron degeneracy value taken as a reference, and this highlights how counter-intuitive differences of profiles can be. At higher electron degeneracy levels, the perturbation looks like a more spread-out version of the perturbation visible in panel a, but its behaviour becomes more complex as the degeneracy decreases. Both helium wells lose their Gaussian aspect: the first well approaches a “heartbeat” shape, that is, it becomes clearly positive (probably under the influence of the hydrogen well), while the second well becomes progressively bi-modal. In contrast, the results in panel c could have been guessed from Fig. 6. The ε value shifts both reference and perturbed profiles, thus resulting in a scaled version of the perturbation shown in panel a.

As a consequence, Fig. 8 shows the issues inherent to calibration methods that we mentioned in Sect. 2. Even if it can be established that the choice of the helium amount Y_s taken as reference does not ultimately matter greatly, Fig. 8 illustrates the diversity of δ_tΓ₁/Γ₁ profiles corresponding to the same helium difference by simply considering different values of ψ_CZ and ε. Reusing the notations introduced in Sect. 2, we clearly find that any change in a component of ⃗θ_⋆ that may impact the best-fitting pair (ψ_CZ, ε)⊂⃗θ_p will lead to substantially different estimates of δY. In this respect, a parameter with a particularly strong impact is given in the forthcoming paragraphs, thus further distancing us from a unique relation of the form ⃗θ_p(δY) as assumed in Eq. (13). Furthermore, it should be added that the particular case of a perturbation around the reference value Y_s = 0 (as considered in HG07) is no exception. Although we only varied the two remaining parameters ψ_CZ and ε, Fig. 9 shows the many possibilities for δ_tΓ₁/Γ₁ even when an identical helium difference δY_s = 0.25 is considered from a pure hydrogen model Y_s = 0. The spread of these randomly drown curves provides an idea of the potential dispersion of the perturbed profile at any point of the structure even under these “favourable” conditions.

Fig. 9. Profile differences δtΓ1/Γ1 again obtained for a fixed value of δYs = 0.25, but from a pure-hydrogen reference model Ys = 0. The 500 profile differences in this plot were obtained by randomly drawing the remaining two reference quantities (ψCZ, ε) in the 2D box [ − 9, −3]×[0.015, 0.02] (the ε interval has been reduced compared to Fig. 8 for clarity). Two opposite corners of this box are highlighted in red and blue.

Secondly, Figs. 8 and 9 underline how accustomed we are to conceive the structural perturbations caused by a variation in helium under solar conditions of electron degeneracy. We may note here the added value of having introduced a model physically; it would have been complicated to imagine and then parameterise these types of perturbations by simply observing a Γ₁ profile in a realistic model. Accordingly, we may question the validity of using ad hoc profiles to approximate the perturbation caused by a change in helium abundance, and this particularly at low electron degeneracy (whose domain of relevance is discussed in the following). It appears to be likely that adjusting functions whose form is inappropriate for these complex profiles may result in inconsistencies regarding their parameterisations.

Up to this point, we only considered a given helium difference, in particular, to determine how the model presented in this paper agrees or contrasts with the types of perturbations for ionisation regions from previous studies. However, Eq. (73) clearly extends to the much more general difference

$$\begin{array}{c}\hfill \frac{{\delta }_t{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_1}=\frac{{\mathrm{\Gamma }}_1(t;Y_s+\delta Y_s,{\psi }_{\mathrm{CZ}}+\delta {\psi }_{\mathrm{CZ}},\epsilon +\delta \epsilon )-{\mathrm{\Gamma }}_1(t;Y_s,{\psi }_{\mathrm{CZ}},\epsilon )}{{\mathrm{\Gamma }}_1}.\end{array}$$(74)