Properties of the ionisation glitch I. Modelling the ionisation region

Context. Determining the properties of solar-like oscillating stars can be subject to many biases. A particularly important example is the helium-mass degeneracy, where the uncertainties regarding the internal physics can cause a poor determination of both the mass and surface helium content. Accordingly, an independent helium estimate is needed to overcome this degeneracy. A promising way to obtain such an estimate is to exploit the so-called ionisation glitch, that is, the deviation from the asymptotic oscillation frequency pattern caused by the rapid structural variation in the He ionisation zones. Aims. Although it is progressively becoming more sophisticated, the glitch-based approach faces problems inherent to its current modelling such as the need for calibration using realistic stellar models. This requires a physical model of the ionisation region that explicitly involves the parameters of interest, such as the surface helium abundance, Y s . Methods. Through a thermodynamic treatment of the ionisation region, an analytical approximation for the ﬁrst adiabatic exponent Γ 1 is presented. Results. The induced stellar structure is found to depend on only three parameters, including the surface helium abundance Y s and the electron degeneracy ψ CZ in the convective region. The model thus deﬁned allows a wide variety of structures to be described, and it is in particular able to approximate a realistic model in the ionisation region. The modelling work we conducted enables us to study the structural perturbations causing the glitch. More elaborate forms of perturbations than those that are usually assumed are found. It is also suggested that there might be a stronger dependence of the structure on the electron degeneracy in the convection zone and on the position of the ionisation region rather than on the amount of helium itself. Conclusions. When analysing the ionisation glitch signature, we emphasise the importance of having a relation that can take these additional dependences into account.


Introduction
Asteroseismology, i.e. the study of resonant modes in stars, reveals information on the physical properties of the layers that the 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 more than any other method.Undoubtedly, the constraints thus obtained depend on what one chooses to model or not and what one assumes to be known or not.This choice is highly complex, and largely dependent on the quality of the information provided (i.e. the precision available on the observables).But, 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) allows ever more precise estimates of oscillation frequencies, the evaluation of certain stellar parameters remains uncertain.This illustrates that accuracy, rather than precision, becomes in this case the limiting factor.Indeed, and 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 being illustrated in Fig. 1).
In this respect, the determination of abundances using realistic stellar models constitutes a textbook case, the latter being 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, one can evoke 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) whose consideration or not may result in physical biases.Recent studies using model grids also suggest that there is 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 high volatility of both parameter estimates.Moreover, in the case of a large number of frequencies being available like that of the Sun, inversion techniques highlight significant discrepancies and solar models are forced to compromise 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 directly involving the composition, one must face additional uncertainties surrounding the near-surface region, namely surface effects (Christensen-Dalsgaard et al. 1988).Due to theoretical 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) highlights what can be the expression of physical processes missing or poorly taken into account in the model, thus resulting in physical biases although having ideal observables.Case (c) provides an example of degeneracy, i.e. parameters that can vary widely because they are too little constrained by the observables.The last panel (d) typically illustrates the solar problem, where no mixture seems to satisfy the constraints.In this case, in order to provide a solution, one must consider inconsistencies in the model's physics or much weaker constraints than actually provided by the oscillation frequencies.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 impact 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, since a reliable estimate of abundances would in turn provide a constraint on the internal physics that depends on it.
The specific influence between composition and seismic properties of a star is long 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, the latter 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 parameters introduced.Also, in order to overcome dependencies on irrelevant disturbances in the frequencies (e.g.surface effects in our case) it has been proposed to study the second differences, d 2 ν, rather than the frequencies directly (Gough 1990): where n and are the associated oscillation radial order and degree respectively.The seismic diagnostic thus defined 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 thus making it convenient to apply 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 considered.Therefore, such a procedure may well lead to nonphysical or inaccurate frequency shifts if it relies on an 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. 2014aVerma et al. , 2019) ) which adds its own uncertainty on the internal physics.This method is therefore largely dependent on the work done beforehand and it is in this perspective that this first paper is intended.A second paper will focus on the seismic effects of the structural perturbation and what information these 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 that respect, we will first give a short description of previous modelling in Section 2, as well as a more in-depth picture of the most commonly used model of the ionisation glitch, namely Houdek & Gough (2007), hereafter HG07.The formalism of our own model of the ionisation region will be introduced in Section  et al. (1994) to model the variations of the acoustic potential (in the case of an overshoot) passing from a convective to a radiative region.For such 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 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 in Houdek & Gough (2007).∆ 2 II and Γ II describe in this case respectively the variance and the area of the distribution.
3 and its structural behaviour studied in Section 4. In Section 5, we propose an analysis of the structural perturbation induced by the model and how it is related to analytical models used for the ionisation glitch modelling.The last part of the article will be dedicated to our conclusions.

Previous seismic glitch frameworks
As briefly mentioned, studying seismic glitches requires a twofold approach in the modelling.One must first consider an expression for the structural perturbation and then infer the shape of the associated frequency shift.

Structural perturbation modelling
The modelling of the perturbation usually boils down to choosing a localised and parametrisable shape for the induced structural change.Thus, depending on the region and the phenomenon under study, various forms can be considered.For the transition between radiative and convective regions, one will generally use very sharp shapes like a Dirac or a step function to reflect a discontinuity in the density derivatives (Monteiro et al. 1994;Houdek & Gough 2007).However, the ionisation region being more spread out, the shapes 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 will 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.
In HG07, the structural change caused by a change in helium abundance is first modelled as a Gaussian perturbation of the first adiabatic exponent Γ 1 = (∂ ln P/∂ ln ρ) S : 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 Γ 1 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 latter seems closely reproduced by a Gaussian expressed as a function of the acoustic depth τ = R r dr /c; c(r) designating the sound speed at a given position.An important point is that, if the well depth, position or width 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 .
Before further discussion, it may be useful here to clarify an ambiguity concerning the notation " δ " which symbolises the perturbation.Indeed, since it denotes a difference between two profiles (say between profiles Γ A 1 and Γ B 1 ), the variable over which these are expressed has to be chosen carefully.For instance, for arbitrary variables x and y: As shown above, the notation " δ x " allows us to overcome this ambiguity.The latter has been introduced by Christensen-Dalsgaard & Thompson (1997) which develop this point extensively.Since the quantity x should vary in the same range for both profiles A and B, it is natural to choose it as a normalised variable e.g.r/R, m/M, t = τ/τ 0 , ... (R, M, and τ 0 = R 0 dr/c being the total radius, mass, and acoustic radius of the model).Also, one may note that the variable on which the perturbation is expressed can differ from the one that has been used to calculate the difference; both δ x Γ 1 (x) and δ x Γ 1 (y) have a meaning.
This point clarified, it appears that the perturbation used in HG07 for Eq. ( 2) is δ τ Γ 1 since the Gaussian shape has been established after an expansion at fixed τ (cf.Eq (31) of HG07).As mentioned, δ τ is however ambiguous if τ A 0 τ B 0 .It can be replaced by δ t following: with δτ 0 = τ B 0 − τ A 0 simply being the perturbation of a constant.

Expected frequency shift
The resulting signature in the frequencies can be divided in two distinct parts: Ionisation component The connection between Eq. ( 2) and the frequency shift δν n, = ν B n, − ν A n, is made by considering the following asymptotic form (Gough 1990): where K n, ρc 2 and K n, c 2 ρ 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.Once more, one must pay attention to the perturbation used in (5); it can be shown that: This nuance cannot, however, be appreciated without the notation introduced here.The frequency shift derived in HG07 is written as a continuous function of frequency: as are the second differences: HeII ν 2 cos 2 2πτ HeII ν + HeII − δ HeII (8) with A HeII = πΓ HeII /τ 0 .F HeII and δ HeII are functions of the other parameters and the frequency (though it is assumed that they only fluctuate slowly with it), thus making the total number of parameters only 4.

Smooth component
The helium component is not the only perturbation expected in the second differences.Apart from the signature of the transition between radiative and convective regions, a smooth component d 2 ν s can also be considered and modelled as: 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), i.e.: with ε and V two dimensionless parameters independent of the radial order n and ∆ν ≡ 1/2τ 0 the asymptotic large frequency separation.The last term of (9), which is proportional to ν −3 , can then be interpreted as the leading term of the second differences applied to (10): 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) list in particular 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 4 free parameters to the ones identified so far.1 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 intricate, and so an analytical derivation of second differences 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 reveal 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 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 instead of Eq. ( 8).Since the parameters of the two Gaussians are not independent, 3 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 still amounts to 4 (only one is needed to determine both HeI (ν) and HeII (ν)).However, despite the substantial improvement in the parameter retrieval, (8) remains currently the most used equation for studying ionisation glitches (Verma et al. 2017(Verma et al. , 2019;;Farnir et al. 2019).
More importantly, the abundances, and in particular the helium abundance Y, do not appear directly as parameters in (12).As mentioned above, without further theoretical work to establish a dependency between these parameters and Y, calibration based on realistic models therefore seems necessary in order to make it appear.This is, in our opinion, the greatest criticism that can be made of the method.To illustrate this point, we will place ourselves in a broader framework and consider two models P and M representing respectively the structural perturbation caused by a change in helium δY and the structure of a realistic stellar model.The first is parameterised by a set of values θ p and the second by the set θ from which we will explicitly distinguish Y for its particular role in this context.As an example, if the model P is chosen to be a Gaussian as done by HG07 (Eq.( 2)), it would then involve three parameters θ p = Γ HeII , ∆ HeII , τ HeII .In contrast, θ will contain instead 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 P to a model perturbation δM and therefore to assume: in order to numerically derive a relation θ p (δY).It can be noticed that Eq. ( 13) is equivalent to Eq. ( 33) of HG07 when P is chosen to be the sum of two Gaussians and Y = 0 as a reference (although the choice of θ is not specified in HG07).Thus, in order to be able to associate a set of parameters θ p to a unique δY, the calibration method makes the assumption that ∂M /∂Y does not depend on the reference point (Y; θ ).However, despite some arguments that will be given in favour of a relative independence regarding the choice of Y, we will show in Section 5 that such an assumption about θ is inadequate.Such dependencies cannot be reflected in model P which depends on too few parameters, namely θ p (and it should be noted that this is actually part of the reasons for introducing P in the first place).To minimise the error introduced, it is then necessary to determine as many θ p (δY) relations as stars studied (given independent constraints on θ , see Verma et al. (2019)), which can be quite costly.Moreover, by using such a relationship anyway, we end up reintroducing the biases we wanted to avoid when looking at the ionisation glitch.Indeed, the relationship obtained will largely depend on the physics considered in M and in particular on mechanisms not necessarily relevant for the ionisation region modelling.
To summarise, despite being based on a more sophisticated model than the previous ones, the method described above faces challenges 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 physically based parameterisation allowing inferences concerning ionisation regions and secondly to give a mathematical description regarding the Γ 1 profile.The challenge will be to provide such a modelling while keeping the number of parameters as low as possible.This whole issue will be addressed in the two following sections.

First adiabatic exponent in an ionisation region
In order to model the structure of the ionisation zone, we will first try to provide an analytical expression of Γ 1 in this specific region.As previously mentioned, we will avoid introducing it via similarities with mathematical functions but rather derive it from thermodynamic relationships.Such expressions 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, due to the complexity of the chemical equilibrium resulting from Saha's equations, these equations are generally solved numerically, and analytical expressions are limited to the ionisation of single species.We will try to present in the following the simplest model that can nevertheless explicitly involve chemical composition.

Free energy
Our starting point is the free energy of a multi-component ideal gas, i.e.: where V designates the volume, kT the temperature expressed in energy units, and N α and µ α (T, V, N α ) 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 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 via 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 Γ 1 , to various compressibilities.Before diving into the calculations, we would like to take a closer look at the particular meaning of a derivative in stellar physics compared with pure thermodynamics.The first adiabatic exponent is often defined as follows2 : 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: Although this equation is generally valid, its use in the present case will require some clarification as 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 = r N r i 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) & ( 16) are generally ambiguous unless some indication is given concerning the N sp (N sp + 3)/2 + 1 implicit conservation equations.It is usual to consider partial derivatives at constant values of state variables, e.g.: 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 Γ 1 profile used in stellar physics shows telltale signs of ionisation that derivatives taken at constant N α do not.Indeed, 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 instead consider the following conservation equations: The first one corresponds to the chemical equilibrium of the ionisation reaction: A r−1 i A r i + e − .The second 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, i.e. 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 will hereafter be denoted with the subscript "EQ".The resulting number of constraints has been indicated in square brackets and add up to N sp (N sp + 3)/2 + 1 as required.
Hereafter, we will use the notation ∂ 2 αβ F to designate a second derivative taken with respect to β (subject to the N α being constant) and then α (subject to EQ).For example, Eq. ( 16) becomes: by using relations (17).Considering Eq. ( 14), the evaluation of the Eq. ( 21) now boils down to finding the N r i (T, V) at a given EQ.

Approximate local equilibrium
To solve the system (18),( 19), ( 20), one must first explicitly write out the expression of the chemical potential of each particle.We have: with λ α (T ) and Z α (T ) the thermal De Broglie wavelength and the canonical partition function of a particle α, respectively.Conditions (18) thus become: ).In the ideal case3 , Z r i and Z e are often approximated by (e.g.Kippenhahn et al. 2012): with u r i describing the fine structure of state (i,r) approached by the degeneracy of the ground state g r i and χ r i the ionisation energy separating state (i,r) from (i,r − 1) ( r s χ s i /kT is then the energy from the reference state (i,0) 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: These equations are highly coupled, especially because of the term N e that can be expressed from (20) as: with N = i N i and x i = N i /N respectively the total number of atoms and the number abundance of element i, which are two constants from ( 19).Because of this, Saha's equations are generally solved numerically by iterating over estimates of the x r i .In the following, we will try to give analytical approximations by considering a few simplifying assumptions.First, let us introduce the cumulative occupation number y r i = s≥r x s i , more suited to the study of Eqs. ( 26) (Baker & Kippenhahn 1962).Also, we will use from this point the convention: with the convention y z i +1 i = 0 (and by definition y 0 i = 1).Let us start by noting that it is rare to find more than 2 ionisation states simultaneously for a given element.It is thus advantageous to make the following approximation (Baker & Kippenhahn 1962): Then, in the single domain of interest (y r i 0), (29) will be approached by: ∀i, ∀r > 0, A final assumption is needed to express z as a function of y r i 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 1 ∼ 0.89 − 0.92 and it can be seen using ( 28) that after a complete ionisation of the hydrogen z > x 1 .Another bound can then be obtained on z since the latter tends asymptotically towards  33) and ( 34) in the case i = 1, 2 (opaque curves) with numerical solution of (26) (dashed light curves) for solar near-surface conditions of (T, V, X i>1 ).Note that these quantities are represented as a function of the normalised acoustic depth t = τ/τ 0 (making the surface located at t = 0), along which all future plots will be represented.
mass unit, and since the mean mass m 0 < 2m u in a stellar mixture, we have m 0 /m i < 1/z i .Finally, the number and mass abundances are related by x i = m 0 X i /m i meaning that z should remain below the upper boundary 2 − x 1 + Z ∼ 1.1 with Z = i>2 X i the metal abundance in mass.This simple reasoning allows us to convince ourselves that z must remain close to 1 once the ionisation of hydrogen completed, and in a way relatively independent of the mixture under consideration.
A pragmatic assumption that can therefore be made is that the ionisation of hydrogen happens first and finishes before any other ionisation begins.Thus we can use the fact that hydrogen is the dominant component to approximate z by the mean electron number in a full hydrogen model.Despite not being entirely true (and particularly for a low state of ionisation), this assumption greatly simplifies the problem.Indeed, for all elements i > 1, z ∼ 1 which is, as discussed above, a reasonable approximation.For the ionisation equation corresponding to hydrogen (i = 1), we will instead approximate z ∼ y 1 1 .Our version of Saha's equations can thus finally be written as: with n = 2 if i = 1 and n = 1 in any other case.Let us define K r i (T, V) as the term on the right hand side (RHS) of (32).We then obtain the desired relationships: One can easily retrieve the x r i (V, T ) from there using x r i = y r i − y r+1 i (as well as the N r i (V, T ) using N r i = N i x r i ).Fig. 4 shows a comparison of hydrogen and helium ionisation fractions x 1 1 , x 1 2 , 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 seems satisfying considering the number of assumptions made.One may note, however, that x 1 2 is slightly underestimated because z is actually a little lower than 1 in the first helium ionisation region.The opposite reasoning holds for the second helium ionisation region.

First adiabatic exponent
Approximations ( 33) and (34) allows us to derive all derivatives appearing in (21).First of all, we have: where we explicitly made the y r i apparent to fully exploit relationships ( 33) and ( 34).It is straightforward to see that we simply retrieve the classical expression for the pressure in an ideal gas.Second derivatives taken at EQ ensue: More explicit expressions can be deduced by writing ∂ ln α and exploiting Eq. ( 32) δ j i standing for Kronecker's symbol.Resulting expressions of the derivatives are: The calculation of the two remaining derivatives is however more subtle.We propose here to start with a determination of the energy E: , one might be tempted to take directly the derivatives of (43) using dE = T dS − PdV.However, in the present framework: Clearly, the last two sums cancel each others since ∀i, r>0 dx r i = d(1 − x 0 i ) = −dx 0 i , which is verified in our derivation as long as we define x 0 i ≡ 1 − y 1 i .Electroneutrality, used in line 2 of Eq. ( 44), is also verified by imposing N e ≡ N ir x i y r i .According to Eq. ( 18), µ r i + µ e − µ r−1 i = 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 thanks to departures from equality of partial mixed second derivatives of state functions like the free energy F).However these residuals remain small and we will consider here that chemical equilibrium is perfectly satisfied.Therefore, the only part that remains from (44) is the following well known identity: This means that T ∂ 2 VT F = −(∂E/∂V) T,EQ − P and T ∂ 2 T T F = −(∂E/∂T ) V,EQ .It follows from Eq. ( 43): Denoting the bracket part of ∂ 2 αβ F as ∂ 2 αβ f , one obtains: where The various symmetries present in these equations may give the impression that expression (48) can be further simplified.In fact, we show in appendix A that Γ 1 can finally be written as: with 0 ≤ γ 1 < 1 simply expressed as: Together with Eqs. ( 32)-(34), Eqs. ( 52)-( 54) provide an analytical approximation of Γ 1 (T, V, x i ) in the CZ.Even if this expression is applicable to any kind of mixture (though assumed to be homogeneous in the CZ!), the following examples will be based on a hydrogen-helium mixture of mass fraction 1 − Y and Y (Y = 4x 2 /(1 + 3x 2 )) to facilitate the study.Remembering that N/V = ρ/m 0 , we gave a representation of the resulting Γ 1 (ρ, 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 nevertheless provides an idea of what part of the map is visible on a Γ 1 profile.The figure clearly reflects the Γ 1 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, one can notice the superposition of the hydrogen and helium first ionisation zones in typical solar conditions whereas the distinction becomes apparent at lower densities.In contrast, this signature seems to be more diffuse at the highest densities.This behaviour appears in each panel and will be further discussed in Section 4. The last frame shows the variation δ ρ,T Γ 1 caused by a change of helium abundance, δY = 0.1, from a reference value Y = 0.25.One may note that, even if the variables are not normalised, the notation "δ ρ,T " is still meaningful since Γ 1 does not refer here to a profile but a function defined for all (ρ, T ) at a given Y.As expected, such a variation in abundance results in a lowering of the two helium wells while elevating the hydrogen one.However, the variations seem fairly disproportionate between the two elements.In fact, 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 the helium one.In particular, the amplitude of the variation seems more related to the relative change δx i /x i than δY.While we have δx 1 = −δx 2 for the hydrogen-helium case, we see that δx 1 /x 1 = −δx 2 /x 1 = −(x 2 /x 1 )δx 2 /x 2 with (x 2 /x 1 ) = 1/12 for Y = 0.25.This explains why the change in hydrogen is expected to be an order of magnitude lower than the helium one.
Although Eqs. ( 53) and (54) seem promising, they are merely functional forms and do not allow us at this point to model a parameterised structure of the ionisation region.However, Γ 1 , ρ and T are closely related in stellar interiors.

Model structure and properties
In this section we will obtain the stellar structure associated with the Γ 1 expression derived in the previous section (see Eqs. ( 53) and ( 54)).For that, we place ourselves in the conditions of a CZ, which correspond to the assumptions while deriving the first adiabatic exponent: 1. Isentropic region 2. Chemical equilibrium at any position 3. Uniform abundances (due to the homogeneous mixing) 4. Electroneutrality at any position The last 3 correspond to applying EQ conditions at any point in the region.Also, to avoid getting lost in the many possibilities offered by these relations in terms of the mixture, we will once more consider a relation of the form Γ 1 (ρ, T, Y) in this part.It is however possible to generalise all the observations made here to any (reasonable) mixture by replacing Y by X i>1 .

Ionisation region structure
As mentioned above, the quantities Γ 1 , ρ and T are not independent in stellar interiors and, although it may be interesting where we have introduced the additional intermediate variable g, the norm of the gravity field.At this point, it is important to keep in mind here that we are trying to account for the variation of T and ρ in the region rather than that of P.This can be achieved by noticing that, for any quantity α, Then, crucially, we note that derivatives taken with respect to the radius r such as d ln P/d ln α = (d ln P/dr) (d ln α/dr) −1 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 3 equations, the solution of which is (T, ρ, g): where we introduced the second adiabatic exponent Γ 2 such that . Its expression can easily be deduced from Γ 1 by noticing that their definitions are symmetrical with respect to a change V ↔ T : Note that an analogous expression to Eq. ( 53) can be provided for Γ 2 : by introducing One can also get rid of the pressure using the ideal gas law (35) which leads to the following system: 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 . 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 actually possible to parameterise the model by considering other triplets derived from these 3 quantities.Hereafter, we will consider instead the following dimensionless parameterisation: 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, one can see that this choice entirely constrains the Γ 1 value at the centre, the latter being only a function of Y, 2m u /ρλ 3 e (we keep in mind that λ −2 e = 2πm e kT/h 2 ), and χ H /kT (cf.Eqs. ( 32)-( 34) & ( 52)-( 54)).However, their impact on our model is much more significant than that.

Physical interpretation of the parameters
-Since 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 gets its notation from its close link to the electron degeneracy parameter ψ = µ e /kT = ln(n e λ 3 e /2) in the ideal gas limit.In fact, one can show that: However, as long as Γ 1 = 5/3 (i.e.outside the ionisation region), ρ ∝ T 3/2 which means that ψ is preserved.Therefore, ψ CZ actually 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 actually reflect properties of the ionisation region; respectively its helium abundance and electron degeneracy.It is complicated to do so with ε since χ H /kT varies by several orders of magnitude in the model.We will however try to give an intuitive understanding of its impact on the model structure.The first thing to notice is that, if Γ 1 did not depend on χ H /kT , the latter would be equal to 5/3 in the whole model.Indeed, since Γ 1 would depend only on Y (which is constant) and 2m u /ρλ 3 e the value of which is fixed as long as ρ ∝ T 3/2 , the first adiabatic exponent would keep its central value, i.e. 5/3.One can then see that the value of χ H /kT at a given position alone determines whether Γ 1 deviates from the value 5/3, i.e. if one enters the ionisation region starting from the model centre.
In fact, it can be shown that the ionisation of the element i in state r takes place when: χ r i /kT ∼ −ψ.The role of ε can now be clearly identified: if one sets ε to be large (say 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 end and the region will be shifted closer to the surface.Thus, despite not being explicitly linked to any quantity of the ionisation zone, ε controls the position of the region.
Having given a physical interpretation for the parameters, we will now try to provide a geometrical interpretation of the impact of each parameter regarding the Γ 1 profile.We will rely on Fig. 6, which shows the Γ 1 (t) (with t = τ/τ 0 , the normalised acoustic depth) profiles resulting from multiples parameter sets (Y s , ψ CZ , ε) by varying them one by one around a reference value.
Geometrical interpretation of the parameters regarding the Γ 1 profile -We gave a representation of Γ 1 (t) for all possible values of Y s in panel (a) of Fig. 6, including such extreme cases as a full pure hydrogen or helium stars.As expected, we go from a profile showing only the hydrogen well (Y s = 0) to one showing only the two helium wells (Y s = 1), passing through profiles that account for the 3 components (0 < Y s < 1).
Clearly, increasing Y s has the effect of extending both helium wells while reducing the hydrogen one.-The effects of modifying ψ CZ can be seen 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).Once more, the impact is 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 gets closer to −12.In fact, this last point is consistent with what has been seen in Fig. 5 by noticing that high degeneracies correspond to relations of the form ln ρ(T ) = C + (3/2) ln T 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 one of HeII being more affected than the hydrogen one.In order of magnitude, this effect can be explained by the dependence of the normalised acoustic depth on temperature close to the surface: t ∝ √ T/T c ∝ √ εT .Thus, for a given ionisation temperature T ion , the impact of ε will become more noticeable on the associated depth t ion as T ion is high.This reasoning also explains why the shift does not seem linear with respect to ε for a given well (this is particularly visible for the HeII well); the latter actually varies as √ ε.
Regarding Fig. 6, an interesting point of the model defined with (Y s , ψ CZ , ε) is the relative independence of the impact of each parameter.Indeed, by treating each well as a distribution, all parameters seem to impact a specific moment.For instance, Y s seems directly related to the distribution area, and does not impact much the shapes of the wells nor the position (for reasonable values of Y s ).ψ CZ is clearly related to the dispersion of  the distribution and ε to its position.Such distinct impacts on the structure are encouraging for subsequent parameter retrieval as one can expect to avoid degeneracy biases.

Comparison with a realistic stellar model
The first and most important point to address is to verify to what extent the model described in this article can approximate the ionisation region structure of a realistic stellar model.To that end, we compare in panel (a) of Fig. 7 Γ 1 (t) profiles extracted from a CESTAM solar model and from our model for a particular parameter set (Y s , ψ CZ , ε ).These values have been chosen so that is minimal.One can see that the behaviour of the realistic model can be very well reproduced, the curves differing mainly in only 2 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 mainly results from the assumption made in Sect.3.2, where we assumed that z ∼ 1 in the helium ionisation region and thus underestimated the number of electrons in the HeI region.This causes the ionisation to start a little deeper in the model and therefore to be more localised.By reducing the dispersion of the HeI well, the assumption made its depth slightly greater.As expected, another difference can be seen in the atmospheric region (t < 0.03).Because the ionisation has not yet started, our Γ 1 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 3 parameters need to be adjusted.
A second point to be clarified is the exact meaning of the adjusted values (Y s , ψ CZ , ε ) regarding the realistic model.Having the minimal differences (δ t Γ 1 /Γ 1 ) 2 , the model parameterised by (Y s , ψ CZ , ε ) is interpreted as being the closest to the CESTAM model from the perspective of a particular structural aspect4 .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 ρ will be described by γ ≡ d ln P d ln ρ rather than Γ 1 , 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, i.e. in the CZ.For the helium abundance Y s , since Y s = Y(t) at any position of our model, the value Y s 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 Γ 1 = 5/3, so we expect ψ CZ to approach ψ CESTAM (t) in the CZ part located below the ionisation region.As mentioned, since ε 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: with BCZ, ATM and ION designating respectively the base of the convective zone, the base of the atmosphere and the end of the ionisation region.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 and ψ CZ obtained through the minimisation of Eq. ( 70).We also highlight the regions were the quantities are expected to correspond to each other based on Eqs. ( 71) and ( 72).The values thus obtained (Y s , ψ CZ , ε ) = (0.255, −5, 0.015) are consistent with the CESTAM quantities in the associated regions (see panels (b) & (c) of Fig. 7 and Table 1).
It is risky to be quantitative about this result since there are no obvious values with which the difference can be compared; in particular the relative difference does not seem 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.

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, as can be seen from Eq. ( 5), the analysis of frequency shifts relies more on the modelling of structural perturbations than on the structure itself.We will 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: 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 will naturally depend on the helium difference δY s but also on the point, (Y s , ψ CZ , ε), around which the differences are calculated, thus leading to 4 dimensions to explore.Our analysis of the threedimensional space required for the structure was obviously already incomplete in Fig. 6, and it follows that the perturbation analysis will necessarily be 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 δY s = α × δY s would be relatively similar to the profile α × δ t Γ 1 /Γ 1 .In this sense, Fig. 8 still provides a good representation of what a change in helium abundance might possibly cause.We see in panel (a) that all the perturbed profiles roughly overlap, which reflects a relative independence from the helium amount chosen as a reference (the range considered, 0.2 < Y s < 0.35, being 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), the second one being more pronounced.The hydrogen contribution is less intuitive, resulting in a very localised peak (t ∼ 0.03, present on every panel) at the beginning of ionisation.The latter 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, it should be noted that this component is clearly noticeable in terms of amplitude.
Things get less intuitive in panel (b).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 in Fig. 11).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 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 one visible in panel (a), however its behaviour becomes more complex as the degeneracy decreases.Both helium wells lose their Gaussian aspect: the first one gets closer to a "heartbeat" shape thus becoming clearly positive (probably under the influence of the hydrogen well) while the second 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 first warns us about the issues inherent to calibration methods already mentioned in Section 2.Even if it can be established that the choice of the helium amount Y s taken as reference does not ultimately matter much, Fig. 8 illustrates the diversity of δ t Γ 1 /Γ 1 profiles corresponding to the same helium difference by simply considering different values of ψ CZ and ε.Reusing the notations introduced in Section 2, it is clear that any change in a component of θ that may impact the best-fitting couple (ψ CZ , ε) ⊂ θ p will lead to substantially different estimates of δY.In this respect, a parameter with a particularly strong impact will be 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 Γ 1 /Γ 1 even when considering an identical helium difference δY s = 0.25 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.
Secondly, Figs. 8 & 9 underline how accustomed we are in conceiving 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 such types of perturbations by simply observing a Γ 1 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 will be discussed in the following).Indeed, it seems likely that adjusting functions whose form is inappropriate for these complex profiles may result in inconsistencies regarding their parameterisations.
Up to this point, we have only dealt with a given helium difference, in particular to see how the model presented in this paper agrees or contrasts with the types of perturbations for ionisation regions from previous studies.However, it is clear that Eq. ( 73) extends to the much more general difference: Indeed, when applied to frequency shift analysis, it seems unlikely that the model differs from its target only in its helium abundance.Therefore, in the following, we propose to study the structural impact of a δψ CZ or a δε difference.It is obvious that the problems mentioned above in terms of the number of dependencies will only be further amplified.For this reason, we fixed the reference set of parameters by taking for example the values (Y s , ψ CZ , ε ) and varied one by one the differences δY s , δψ CZ and δε.The resulting profiles are shown in Fig. 10.As expected, the variation observed in panel (a) is very similar to the one in panel (a) of Fig. 8 scaled by a factor α = δY s /0.05.The other two perturbations, however, take different shapes and notably make the wells unrecognisable.A change in the electron degeneracy (panel (b)) has a strong impact on the hydrogen ionisation structure, especially at a point (t ∼ 0.09) where a change in helium has almost no effect.The ionisation region shift is also enhanced while the second helium ionisation contributes to a weaker and more localised well.The variation in the position δε of the ionisation region is also fairly counterintuitive.Indeed, panel (c) illustrates a contribution that reaches its maximum between the two helium wells (t ∼ 0.15), the one corresponding to the peak in the Γ 1 (t) profile.
A first observation to highlight from Fig. 10 is how elaborate the ionisation perturbation profile becomes in the general case.It should be kept in mind that a structural perturbation from a reference (Y s , ψ CZ , ε ) can be approximated as a combination of profiles composing Fig. 10.Indeed, from Eq. ( 74): with e.g., being the structural perturbation studied in panel (a) within a multiplicative factor.With regard to Fig. 10, the result of the linear combination (75) can become highly complex and is unlikely to have anything in common with the functions generally used to fit δ t Γ 1 /Γ 1 .Combined with the issues already raised about calibration, this tends to suggest that the study of the ionisation glitch based on ad hoc perturbation profiles in order to calibrate them afterwards might lead to fairly inconsistent results.Also, we have tried to represent in the figure structural disturbances that are similar in order of magnitude.This allows us to qualitatively relate the values of the different parameters.In this sense, we see that perturbation with δY s = 0.1, δψ CZ = 1 or δε = 0.002 have a comparable structural impact in amplitude and hence on the frequencies.However, although the first one nearly spans the entire range of relevant values (the difference in helium between two stars can hardly exceed 0.1), this is not the case for the two others.In fact, regarding the ε parameter, this change is actually very small: two consecutive curves in panel (c) of Fig. 6 are separated by exactly δε = 0.002 which is the entire span for the perturbation represented in Fig. 10!This last point suggests a very high dependence of the perturbation on a shift of the ionisation region.As mentioned above, it is interesting to note that the glitch signature would come from the Γ 1 peak rather than from a well if this contribution were to dominate.This is consistent with observations already made in previous studies (Broomhall et al. 2014;Verma et al. 2014b), although it provides an alternative explanation.
To put into perspective the significance of a difference in the ψ CZ parameter, it can be interesting to relate it to more intuitive quantities like fundamental parameters.Since the electron degeneracy ψ is known to fluctuate with the mass of the star (e.g.Hayashi et al. 1962), we represented in Fig. 11 various ψ(t) profiles (along with the corresponding Γ 1 (t) profiles) obtained from CESTAM models of different masses.The latter cover a wide range of values from 0.7M to 1.35M appropriate for solar-like pulsators.All the models also share the same composition (which is the solar one) and subsequently the same helium abundance as well as the same evolutionary state (X c = 0.6).Clearly, despite having an identical composition, these models show fairly distinct Γ 1 (t) profiles in the ionisation region.In fact, these curves show a fairly similar behaviour to those shown in panel (b) of Fig. 6 (which can be extrapolated to the wider range 0.45M < M < 1.5M presented in dashes as an indication).Thus, at small masses, the wells seem to merge and then split at higher masses making in particular the HeI well visible.In addition to this effect, an expansion and contraction of the profiles is visible depending on the mass range considered (< 1M or > 1M , cf. arrows shown in panels (a) & (b)).Looking at the similarities with Fig. 6, we interpret these changes as manifestations of variations in ε and the electron degeneracy just below the ionisation region.Although the variations of ε can only be evaluated qualitatively (remember that this quantity is only defined in entirely convective models), this is not the case for the electron degeneracy.Indeed, ψ(t) profiles typically undergo little to no variation in convection zones as shown by the plateaus in the curves in panel (c), thus making it easy to identify ψ CZ .The range of ψ CZ values covered by these plateaus corresponds to the range studied in Fig. 8 (∼ −3 for 0.7M stars and ∼ −9 for the 1.35M model).It is then possible to seek to characterise the significance of the change δψ CZ = 1 mentioned above.To provide an order of magnitude, we represented in panel (c) the ψ span corresponding to the mass uncertainty (1.14 < M * < 1.32) of HD 52265 taken from Lebreton & Goupil (2014).The latter provides an example of a seismically determined mass uncertainty taking into account a large number of biases as discussed in Section 1.The corresponding extent in terms of electron degeneracy is about δψ CZ ∼ 2, which is 2 times larger than the range considered in the panel (b) of Fig. 10 and thus larger than   (c)) profiles extracted from CESTAM early main sequence models with masses ranging from 0.7M to 1.35M (and from 0.45M to 1.5M as an indication in dashes).All models share the same composition which is the solar one as well as the same evolutionary state: X c = 0.6.The qualitative variations of ε with mass are represented as well as the ψ span corresponding to the mass uncertainty of HD 52265 taken from Lebreton & Goupil (2014).
any possible variation in helium abundance.Note that this range in ψ CZ depends on the mass; such a mass difference would correspond to a lower δψ CZ at 0.7M and would instead increase at higher masses.In this regard, a mass difference between two models5 could contribute to a perturbation of the structure of the ionisation region as large as (if not larger than) that of helium.
Finally, to further echo the point made in Fig. 8 about calibration, the mass taken as a reference, M , provides a representative example of a component of θ that impacts the (ψ CZ , ε) pair as already observed by Verma et al. (2014b).In this article, the authors proposed to use it by considering different calibrations to determine the helium abundance of stars of distinct masses.However, besides the fact that this presupposes a relatively good knowledge of the star in question before being able to estimate its helium abundance, we have seen that stars whose masses are already seismically constrained can possess significantly different degeneracy levels in the convective region.Using the example of HD 52265, the question then arises whether it is possible to obtain consistent results using calibration and being able to consider both models with ψ CZ ∼ 6 or ψ CZ ∼ 8, which implies perturbation profiles that differ in both form and magnitude for a fixed difference in helium amount (see panel (b)) of Fig. 8).This statement does not apply exclusively to the mass, M , but can reasonably be extended to any fundamental parameters or physical processes (i.e. the other components of θ ) that may affect the average density in the convection zone (related to ψ CZ ) or the relative acoustic depth of the ionisation region (linked to ε).
Besides this last remark, by comparing the structural perturbations in Fig. 10, it seems likely that the component caused by a helium change between two stars is not the only one to contribute.In fact, it could even be dominated by a change in electron degeneracy or in the ionisation region position.From the perspective of frequency shifts, it thus seems incorrect to assimilate the ionisation glitch as a repercussion of a helium difference only in order to infer its value.In this light, it appears essential to look at forms of frequency shifts that can take into account these additional dependencies, thus relying on physical models as proposed in this paper.

Conclusion
Determining abundances can be subject to many biases in its classical approach, which we intend to overcome by exploiting the ionisation glitch.However, although progressively becoming more sophisticated, the glitch-based approach faces problems inherent to its current modelling such as the need for calibration by realistic stellar models.
To address these problems, a physical model of the ionisation region is proposed here, explicitly involving the parameters of interest, such as the surface helium abundance Y s , and which can be generalised to more elaborate compositions.In the case of a hydrogen-helium mixture, the model involves 3 parameters and highlights the importance of a characterisation of the electron degeneracy state in the convective zone ψ CZ as well as the position of the ionisation region here controlled by ε.While it is well known that the abundance of surface helium contributes greatly to the appearance of the first adiabatic exponent (Γ 1 ) profile by shaping the size of the helium wells, the state of degeneracy seems to affect the profile just as much by altering their dispersion.Taking this into consideration, the model is thus able to describe a wide variety of ionisation structures while providing them with a physical meaning.By comparing them to a CES-TAM model, we also verified that it was able to approach realistic Γ 1 profiles for consistent parameter values, conferring it a predictive capacity.
The modelling work conducted allowed us to study the shape of structural perturbations by analysing differences of profile with distinct parameterisations.In particular, the form of the perturbations caused by a helium difference was addressed.Expected shapes were found such as a Gaussian-like contribution for the helium ionisation regions with in addition a sharper com-ponent at the surface, caused by a shift in the hydrogen ionisation region.However, we have observed that this form is only valid for a restricted case, namely a helium difference under solar electron degeneracy conditions.This perturbation seems, however, highly variable with respect to these conditions, and can easily become much more complex at lower electron degeneracies.This point stresses, in particular, the major dependence of the perturbation on the choice of realistic models chosen for calibration purposes.We then went on to study more general parameter differences.More elaborate forms of perturbations than the ones usually assumed are found.It is also suggested that there is a stronger dependence on the electron degeneracy in the CZ or on the position of the ionisation region than on the helium amount itself.Also, the connection between electron degeneracy and stellar mass therefore enables us to clarify the strong dependence of the helium glitch amplitude on the stellar mass already observed by Verma et al. (2014bVerma et al. ( , 2019)).Moreover, the fluctuation in the ionisation region position thus induced seems compatible with a glitch signature coming from the peak of the first adiabatic exponent as reported by Broomhall et al. (2014) or Verma et al. (2014b).
When analysing the oscillation frequencies, we therefore emphasise the importance of having a relationship that can take into account these additional dependencies.In this sense, a second paper based on the introduced modelling is planned in order to derive more general forms of frequency shifts.The objective should be to interpret the ionisation glitch as a combination of multiple contributions for which a variation of helium abundance is one component among others.
i δx i = 0).Considering Eqs. ( 32)-(34), the y r i thus defined are functions of T and V but not of x i .It can thus be established that δ V,T (y r i ) = 0.The induced perturbations of χ 2 and ∂ 2 T T f will therefore be given by: δ V,T (χ 2 ) = ir δx i x i x i y r i (1 − y r i ) .12) thus leading to the following impact on γ 1 : For the same reasons as mentioned above, we will consider that these two terms are negligible in the following.Since δ V,T (Γ 1 ) = −2/3δ V,T (γ 1 ), an analytical expression of the perturbation is given by: .
To simplify the derivation, we will just present the parts that do not result in terms containing x i x j y r i y s j (1 − y r i )(1 − y s j ) which will be neglected: The resulting perturbation of the numerator is approximated by: And finally the overall perturbation on δ ρ,T (Γ 1 ) will thus be given by the two sums:  5. We have also represented the perturbation given by (A.17) which seems closer to the shape that can be expected for the perturbation, at least based on the panel (b) of Figure 5.The fact that the perturbation moves towards the borders of each ionisation region when passing from δ V,T (Γ 1 ) to δ ρ,T (Γ 1 ) is then interpreted as the result of a shift in the y r i caused by perturbing at constant ρ.In other words, at a given density, a perturbation on the abundances makes the ionisation regions move and this results in a less intuitive map.This effect can however be seen as "artificial" since it is obtained by fixing a variable that mixes the specific effects of various state variables (i.e.V and each N i ).where the "∼ ir " notation stands for the dominant term in the ionisation region (i, r).Comparing this with the numerator of Eq. (A.17) gives an intuition of why the order of magnitude to keep in mind for the perturbation is δx i /x i rather than δx i .Thus, while we have δx 1 = −δx 2 for the hydrogen-helium case, we see that δx 1 /x 1 = −δx 2 /x 1 = −(x 2 /x 1 )δx 2 /x 2 .However for Y = 0.25, (x 2 /x 1 ) = 1/12, resulting in a perturbation about 12 times lower in the hydrogen region.This estimate is consistent with

Fig. 2 .
Fig. 2. Various shapes along with their associated parameters used to describe a structural perturbation.(a) Dirac function used in Monteiroet al. (1994)  to model the variations of the acoustic potential (in the case of an overshoot) passing from a convective to a radiative region.For such 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 inMonteiro & 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 inHoudek & Gough (2007).∆ 2 II and Γ II describe in this case respectively the variance and the area of the distribution.

Fig. 3 .
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 these causes a deviation from the Γ 1 reference value of 5/3.
one can write N e = N z with: of electrons per atom.Saha's equations can be rewritten: ∀i, ∀r > 0,

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

Fig. 5 .
Fig. 5. Panels (a), (b) and (c): Γ 1 (ρ, T, Y) as a function of ρ and T for different values of Y. Panel (d) represents the variation δ ρ,T Γ 1 caused by a perturbation δY = 0.1 from a reference value Y = 0.25.The dashed line present on each panel shows the relation ρ(T ) extracted from a CESTAM solar model.

Fig. 7 .
Fig. 7. (a) Comparison of Γ CESTAM 1 (t) extracted from a CESTAM solar model (in grey) and Γ 1 (t; Y s , ψ CZ , ε ) obtained from our model (in red).(b) Comparison of Y CESTAM (t) extracted from the same model with Y s .(c) Comparison of ψ CESTAM (t) with ψ CZ .In the last two panels, regions where the curves are expected to correspond are shaded in red.

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

Fig. 11 .
Fig. 11.Γ 1 (t) (both panels (a) & (b) for clarity) and ψ(t) (panel (c)) profiles extracted from CESTAM early main sequence models with masses ranging from 0.7M to 1.35M (and from 0.45M to 1.5M as an indication in dashes).All models share the same composition which is the solar one as well as the same evolutionary state: X c = 0.6.The qualitative variations of ε with mass are represented as well as the ψ span corresponding to the mass uncertainty of HD 52265 taken fromLebreton & Goupil (2014).
17) Appendix A.3: Perturbing Γ 1 at given ρ, T valuesThe perturbation at given ρ and T values is more complex since δ ρ,T y r i 0 because of the variable change ρ/m 0 = N/V.For any quantity α, the resulting perturbation change can be found by applying:δ ρ,T (α) = δ V,T (α) − δ V,T(ρ) can now find the perturbation of any quantity expressed in terms of y r i using∂ ln y r i ∂ ln ρ = d ln y r i d ln K r i ∂ ln K r i ∂ ln ρ = − 1 − y r i 1 + δ 1 i (1 − y r i ) neglected many terms to obtain Eq. (A.21), this equation still gives a useful approximation.For instance, the difference between Γ 1 (ρ, T, 0.35) − Γ 1 (ρ, T, 0.25) (shown in Fig. 5) and the relation (A.21) is bounded by 10 −2 at any position of the (ρ, T ) plane and one can verify that the panel (b) of Fig. A.1 is actually very similar to what has been found in the panel (d) of Fig.

Table 1 .
Inferred electron degeneracy parameter and helium abundance and actual CESTAM quantities Note that our ψ estimate results from the combination ψ CZ − Y s /2 in accordance with Eq. (69).
13) Article number, page 18 of 20 Pierre S. Houdayer et al.: Properties of the ionisation glitch By calculating the numerator explicitly, one can see the apparition of 4 sums: (∂ 2 T T f ) 2 δ V,T (γ 1 ) =