© ESO 2019

1. Introduction

Since the launch of CoRoT (Baglin et al. 2009) and Kepler (Borucki et al. 2010) missions, the scientific community has access to a tremendous amount of asteroseismic data of unprecedented quality. Such data are essential to better constrain stellar structure and evolution and, in turn, improve the characterisation of exoplanets and stellar populations. However, it is essential to develop techniques that are able to retrieve stellar parameters as accurately as possible in order to benefit from the quality of the data. A very complex problem in determining stellar parameters is the model dependency of the results. The results are intrinsically dependent on the input physics such as the equation of state as well as the opacity tables used. It therefore becomes of prime importance to develop techniques that are able to test the influence of the input physics on stellar parameters or even techniques that provide results as model independent as possible.

With such precision in the data, and the precision of the future missions TESS (Ricker et al. 2014) and PLATO (Rauer et al. 2014), studying the signature of acoustic glitches becomes a natural step towards better models. This idea was originally proposed by Vorontsov (1988) and Gough (1990). They both highlighted the effect of a sharp feature in the stellar structure on the frequencies, either directly or on the second differences. Such considerations have already been the subject of several studies. As an example of the several techniques used, Mazumdar et al. (2014) illustrate four techniques either using the second frequency differences or striving to isolate the glitch oscillation directly from the frequencies. Also, most of the current methods focus on the localisation of the helium second ionisation zone or the base of the envelope convection zone (e.g. Monteiro et al. 2000). This is a crucial first step on the way to a better understanding of stellar physics. Indeed, a characterisation of the convective envelope extension allows to get constraints on convection itself as well as on overshooting (both its efficiency and nature). Therefore, the glitch provides the necessary observational data needed to refine current convection theories. Also, the study of the helium glitch should provide information on the helium surface content in low-mass stars. This is essential as, in such stars, it cannot be derived from spectroscopic data. Therefore, a methodology taking advantage of most of the oscillation spectrum aspects is required.

Finally, it is sometimes the case that studies use correlated constraints or discard pieces of information in their studies. For example, in method C from Verma et al. (2014), the information that is not contained in the glitch (the smooth component) is not used directly as the glitch information is isolated to draw inferences. Therefore, the usual indicators (e.g. large separation, small separation ratios, ...) are computed separately from the glitch and no method is proposed to determine properly the correlation between those indicators and that of the glitch.

For those reasons, we propose a method: Whole Spectrum and Glitches Adjustment (WhoSGlAd) that takes as much of the available spectral information as possible into account. This method defines new seismic indicators in such a way that they are as independent as possible and significant from a statistical point of view. To do so, it relies on linear algebra via Gram-Schmidt’s algorithm.

The present paper is organised as follows. We first present the method in a very general and mathematical way in Sect. 2. Then Sect. 3 defines the seismic indicators and their diagnosis power will be used to study solar-like pulsators. In the following section, we demonstrate its ability to extract and analyse the glitches signal. We also show its limitations. We present a first application to the case of 16 Cygni A (HD 186408) in Sect. 5. Let us insist on the fact that we do not present here a thorough study of 16 Cygni A but we rather show an example of the ability of the method to provide constraints on stellar structure. We conclude the paper by discussing the results and detailing future perspectives.

2. Method

In the present section, we describe the method we developed. It aims at using as much as possible of the information available in the oscillation spectrum of a star. Therefore, both the oscillatory and smooth part of the spectrum are simultaneously analysed in a single adjustment. This avoids multiple usage of the same information to draw different inferences. The very strength of the proposed method is that the different parameters obtained will be independent of each other, meaning that their covariance matrix will be the identity matrix. This will allow to build indicators which also are independent of each other and draw statistically relevant inferences. The independence of the parameters will be ensured by using Gram-Schmidt’s (Gram 1883; Schmidt 1907) algorithm. Then, the defined seismic indicators will be used as constraints to provide improved models in the framework of forward seismic modelling (see for example Miglio & Montalbán 2005 for one of the first use of Levenberg-Marquardt’s algorithm to adjust a model to seismic and non-seismic observables). Finally, such models may be used as initial guesses for inverse seismic modelling. (see Roxburgh & Vorontsov (2002a,b) for the application of the inversion technique on an artificial target, which shows the feasibility of such techniques, and Buldgen et al. (2016a,b) for examples of inversions in the case of 16 Cygni A.)

2.1. Mathematical description

Non-radial pulsation frequencies can be mathematically defined by three integer numbers; the radial order n, the spherical degree l, and the azimuthal order m (in this paper, we do not consider the seismic probing of rotation and consider only m = 0). The method we developed – WhoSGlAd – is based on linear algebra in a Euclidean space. The vector space we consider is the set of N observed oscillation frequencies ν_i. The standard deviation for each frequency is written σ_i. Given two frequency vectors x and y we define their scalar product as

(1)

Often in asteroseismology, it is useful to compare two sets of frequencies (e.g. theoretical and observed frequencies) using a merit function defined as

(2)

with ν_th and ν_obs, the theoretical and observed¹ frequencies. Taking advantage of the scalar product defined above and the associated norm, this simply becomes

(3)

In the presence of a glitch, Houdek & Gough (2007) showed that the oscillatory component in frequencies due to the glitch can be isolated from the rest of the spectrum, called the smooth component. Thus, to represent observed frequencies, we define a vector sub-space that is typically a polynomial space – the smooth component – associated with an oscillating component – the glitch –. The analytical formulation of those two components will be given in the following sections. This is very similar to what has been done by Verma et al. (2014).

The method consists in the projection of the observed and theoretical frequencies over the vector sub-space. Then, we define seismic indicators from the projections. Their definitions are given in Sect. 3. To do so, it is useful to define an orthonormal basis over the vector sub-space. This is done via Gram-Schmidt’s orthogonalisation process associated with the definition of the scalar product (1). For more information about this process as well as its equivalent form as a QR decomposition, the interested reader may read Appendix D.

If we write j and j₀ the indices associated with the basis elements, p_j the former basis elements, q_j0 the orthonormal basis elements, and the transformation matrix, we have:

(4)

where the dependence in n and l translates that the basis elements are evaluated at each observed value of the radial order n and the spherical degree l.

It is essential that the projections be done in a specific order to obtain the lowest possible value of the merit function. This will be the subject of the following subsection. Finally, we write a_j = ⟨ν|q_j⟩ the projections of the frequencies over the basis elements. The fitted frequencies are then given by

(5)

Let us add that, thanks to the orthonormalisation, the standard deviations of the coefficients a_j are σ(a_j) = 1 and they are independent (their covariance matrix is the identity).

2.2. Smooth component

Now that the mathematical context is given, we may detail the vector sub-space we selected to fit the smooth component of the frequency spectrum. As the set of observed radial orders and σ_i are usually different for each spherical degree and the smooth component depends on l, the smooth component basis elements depend on l. For each value of l, the frequencies will be projected over the different powers considered. We also point out that, for each spherical degree, the method requires at least the same number of frequencies as of powers considered in the formulation. The polynomials are then of the general form:

(6)

where δ_ll′ is the Kronecker delta which compares two spherical degrees l and l′, p_k(n) is a polynomial in the radial order n and k represents its ordering. We note that the previously defined j is now separated into two indices, the spherical degree l and the ordered power k. For a better understanding, the spherical degree and ordering will be explicitly written for the smooth component transformation matrix as . And the orthonormal basis elements are

(7)

which yield:

(8)

where thanks to the introduction of δ_ll′ in Eqs. (6) and (7), the sum over l′ collapses over one fixed degree l. For the smooth component, we treat separately each spherical degree, the parameters associated to a given degree only depend on the frequencies of this degree.

According to the asymptotic theory of non-radial oscillations (Gough 1986), we have the following formulation of the expected frequencies as a function of n and l:

(9)

where is the asymptotic large frequency separation, c(r) is the adiabatic sound speed, and R_* is the radius at the photosphere of the star.

It follows that the first two polynomials in n (taken from the right hand side of Eq. (6)) used to depict the spectrum smooth component will be:

(10)

(11)

Then, to provide the best fit to the observed frequencies, we methodically tested several combinations of powers to find the set giving the best agreement with the observations – the observations actually referring to theoretical models taken as observed stars in a set of calibrations –, hence the lowest χ² value. We get:

(12)

At this point, it is of prime importance to note that the construction of the basis via Gram-Schmidt’s process will have to be done following the ordering of the degrees because it will allow us to associate the seismic indicators to the projection of the frequencies on the successive basis elements. For example, the projection of the frequencies on the 0 order polynomial corresponds to a fit to a constant value. This estimates the mean frequency value. Moreover, we did not include other degrees as the fit of the smooth component was already very good². Also, adding higher order polynomials to the smooth component might account for some of the glitch oscillating features. It is essential to avoid such a behaviour as the definitions of the seismic indicators, and the inferences we draw from them, will be impacted.

Furthermore, we could also include a regularisation parameter λ (as in Verma et al. 2014, method C) in order to prevent those behaviours. This requires a new definition of the vector sub-spaces. The vectors are now of dimension 2N. The first N components are defined as before while the components from N + 1 to 2N are 0 for both the frequencies and the glitch polynomial, and the second derivative of the polynomial for the smooth part. The vectors are now:

(13)

Then, we have to define the scalar product of x and y as

(14)

Therefore, using the definition (3) with the new scalar product gives another value of the merit function. The inclusion of the regularisation parameter allows to minimise the oscillation of the smooth component as the minimisation of the merit function will lead to a minimisation of the second derivatives of the smooth component. Let us note that, for the regularisation constant to have an influence on the results, it must at least be of the order of the inverse of the frequencies standard deviation. However, we performed some tests with and without regularisation terms and it appears that the method is very stable without it. Moreover, we observed in many cases that the results were degraded when including it (see also Appendix A for an illustrative example). Therefore, it is not necessary to include these regularisation terms to properly extract the glitch in our method. As a consequence, the results presented in this paper do not include such terms. The fact that we use fewer fitting parameters than in (Verma et al. 2014; they consider polynomials up to the fourth degree whereas we only reach the second) might explain that the regularisation constant is not necessary in our case. In addition, we note that using λ = 0 leads to the classical χ² fitting.

Finally, as we have three polynomials for each value of l, we built, for the smooth component, a vector sub-space of dimension 3 × l which equals 12 if we have four values for l (e.g. 0, 1, 2, 3). As hinted earlier, more than three observed frequencies of each degree are necessary to apply the developed method.

2.3. Glitch

The formulations used by Verma et al. (2014) and Houdek & Gough (2007) allow us to fit properly the helium and convection zone glitches but they are highly non linear with respect to the free parameters. Below is the expression from Verma et al. (2014), which we adapted:

(15)

where the first term takes the helium glitch into account and the second, the convection zone glitch. The quantities τ_He, ϕ_He, and 𝒜_He represent respectively the acoustic depth of the second ionisation zone of helium, the phase of the helium glitch, and its amplitude. The same goes for the quantities τ_CZ, ϕ_CZ, and 𝒜_CZ in the case of the base of the envelope convective zone. Finally, c₂ is the rate of decrease in amplitude of the helium glitch with the squared frequency.

Moreover, Eq. (15) is implicit since the frequency appears on the right-hand side. To adjust at best the frequencies, it is therefore necessary to use non-linear least square fitting algorithms (e.g. Levenberg-Marquardt’s method, genetic algorithms, ...) which can be unstable and are very sensitive to the initial guesses on the optimal parameters.

For this reason, we decided to adopt the following linearised functions, expressed as a function of for the helium glitch:

(16)

(17)

with k = (4, 5); and, for the convection zone glitch:

(18)

(19)

To obtain the above expressions we replaced the value of the frequency ν by its first order approximation from the asymptotic formulation (9): . We also defined T_He = τ_HeΔ and T_CZ = τ_CZΔ.

Moreover, we approximated the exponential decrease in frequency by the combination of two polynomials in . The degrees −4 and −5 have been chosen to reproduce at best the decrease of the glitch amplitude towards high frequencies (which was described using a gaussian by Houdek & Gough 2007). To do so, we compared the polynomial formulation with the exponential one (Eq. (15)). Let us add that we adjusted the helium glitch using several sets of degrees and the method proved to be very stable and the results remained good. Both the fitted coefficients and defined indicators values were quasi unaffected by the choice of degrees. The glitch then writes

(20)

(21)

Furthermore, such a formulation allows us to move on with Gram-Schmidt’s process and generate orthonormal vectors to append the smooth component basis. As the glitch is generated in the superficial layers, it should not depend on l. Therefore, we defined the coefficients c and s to be independent of l. By doing so, the vector sub-space associated with the glitch is only of dimension 6 and this sub-space is used to complete the orthonormal basis over which frequencies are projected. We are thus able to write the glitch contribution to the frequencies as

(22)

At this point, one should be aware to the fact that the basis functions depend on the acoustic depths of the second ionisation zone of helium and the base of the envelope convective zone, respectively – through T_He and T_CZ –. We observe that the functions depend non-linearly on the values of T_He and T_CZ. To preserve the linearity of the method, it is necessary to provide values for τ_He and τ_CZ and leave them unchanged to generate the basis and project the observed frequencies over it. In the case of a theoretical model, this estimation is done using the definition of the acoustic depth, that is:

(23)

where R_* is the radius at the photosphere, r_He/CZ represents the radius of the helium second ionisation zone or of the base of the envelope convection zone (in practice we take the corresponding local maximum between the two local minima of Γ₁ due to the partial ionisation of He and H for r_He and the last point below the surface for which ∇ <  ∇_rad for r_CZ).

For observed data, we first generate the optimal model that does not take the glitches into account. We then fit this model for the glitches and retrieve the model values of the acoustic depths as estimators of the optimal values. Therefore, we do not provide any new means to estimate the involved acoustic depths.

Eventually, we may optimise over the values of τ_He and τ_CZ to get the best results. This is done through the use of Brent (1973)’s minimisation algorithm. However, it makes the problem non linear again. Furthermore, the optimised estimations of τ_He and τ_CZ always remain very close to the theoretical value and do not decrease significantly the χ² value. Also, we observed that by using an initial value of τ_He different from that at the Γ₁ maximum (e.g. at the minimum of the helium second ionisation zone) and adjusting it, we found back the value at the maximum. Moreover, we tried to find the value of Y_f giving the best agreement with the helium amplitude (defined in Sect. 3.3.1) observed in the case of 16 Cyg A at fixed values of τ_He corresponding either to the second local minimum or the local maximum of the Γ₁ profile. We noted that the difference between both values of Y_f is smaller than the standard deviation. In addition, freeing τ_He does not impact the surface helium abundance retrieved in a significant way. This stems from the fact that the acoustic depth value changes at most of 10% and has little influence on the measured helium glitch amplitude. Therefore, the influence on the calculated Y_f is negligible as well. Thus we finally decided to give up that last non-linear minimisation for theoretical models.

By keeping a fully linear implementation of the spectrum fitting we guarantee the stability of the algorithm, the independence³ of the parameters obtained via the projection on the orthonormal basis as well as small computation times. This is essential for it has to be included into a non-linear routine that searches for a stellar model accounting at best for the seismic and non-seismic observables.

3. Seismic indicators

The main advantage of the developed method is that it provides – via Gram-Schmidt’s process – fitted coefficients which are independent of each other. It therefore allows us to derive seismic indicators as uncorrelated as possible. We will define the ones that we explored in the current section. To characterise those indicators, we computed their evolution along the grid of models presented in the next subsection and using the set of modes observed for 16 Cyg A.

3.1. Models

The grid of models we used was computed using CLES (Scuflaire et al. 2008a) combined with LOSC (Scuflaire et al. 2008b) stellar evolution and oscillation codes. The models used the FreeEOS software (Cassisi et al. 2003) to generate the equation of state table, the reaction rates prescribed by Adelberger et al. (2011), the metal mixture of AGSS09 (Asplund et al. 2009), and the OPAL opacity table (Iglesias & Rogers 1996) combined with that of Ferguson et al. (2005) at low temperatures. Moreover, the mixing inside convective regions was computed according to the mixing length theory (Cox & Giuli 1968) and using the value α_MLT = l/H_p = 1.82 (where l is the mixing length and H_p the pressure scale height) that we obtained via a solar calibration. Microscopic diffusion was taken into account in the computation by using Thoul et al. (1994)’s routine. For each model, the temperature at the photosphere and the conditions above the photosphere are determined by using an Eddington T(τ) relationship. The models have masses ranging from 0.90 M_⊙ to 1.30 M_⊙ by steps of 0.01 M_⊙ and are in the main sequence phase. Moreover, each model has an initial composition of Y₀ = 0.25 and Z₀ = 0.016 to remain close to the solar case. Finally, unless specified otherwise, the observed frequencies have been corrected for the surface effects using Kjeldsen et al. (2008)’s prescription of which the coefficients a and b have been calibrated by Sonoi et al. (2015).

3.2. Smooth component indicators

3.2.1. Large separation

A commonly used indicator is the large separation which holds a local (i.e. based on the individual frequencies) and an asymptotic definition. To construct an estimator of the large separation, we will take inspiration in the asymptotic definition. In the asymptotic regime (n ≫ l), Eq. (9) is satisfied. We notice that, in this formulation, Δ represents the slope in n of the straight line fitting at best the frequencies. Moreover, to fit the spectrum smooth component, we project the frequencies over the basis in a specific order given by the sequence of degrees (0,1,2). This means that keeping only the expression of order 0 will give an adjustment of the frequencies by a constant term, thus estimating the mean value. Furthermore, if we now keep the expression of first order, we adjust the frequencies to a straight line of which the slope, if we rely on Eq. (9), is Δ. This is the most common way to define the mean large separation in seismic analyses. However, we must note that we have different basis vectors depending on the spherical degree considered. This means that, for each value of l, we will have a different regression to a straight line, therefore a different estimate of the large separation Δ_l. Using expressions (4), (5), (7), and (8), we isolate this slope to write:

(24)

We may finally average these indicators over l to estimate at best the large separation. Knowing that the standard deviation of a_l, 1 is 1, is the variance of Δ_l. The weighted mean of the large separations thus yields

(25)

Finally, we expect from Ulrich (1986) that Δ should be an estimator of the mean stellar density.

3.2.2. Normalised small separation

Two commonly used indicators are the small separations d₀₁(n) and d₀₂(n) of which the definitions are:

(26)

They allow a measurement of the spacing between the observations and the asymptotic relation (9). However, they happen to be sensitive to the surface effects. Therefore, Roxburgh & Vorontsov (2003) suggested to divide these expressions by the large separation in order to minimise such effects. Indeed, they showed that those ratios are almost independent of the structure of the outer layers of the star.

We thus introduce estimators of these ratios. Such ratios represent the spacing between ridges of spherical degrees 0 and 1 for Eq. (26) and degrees 0 and 2 for Eq. (27) in the échelle diagram (Grec et al. 1983). In a more general way, we approximate the mean difference between the ridges of spherical degrees 0 and l by comparing the mean values of the frequencies for those degrees. That is . Assuming expression (9) to be exact, this difference is . Then, we added to the expression to make its value come close to ϵ₀ − ϵ_l. We then obtained the following expression:

(27)

where and are respectively the weighted mean values of ν(n, l) and of n for the spherical degree l in accordance with the definition of the scalar product. In addition, the mean value equals as it is the fitting of the frequencies of degree l to a constant value. Finally, we stress that the above expression is slightly different from Eqs. (26) and (27) as they represent the local spacing between ridges in the échelle diagram and Eq. (28) corresponds to the mean spacing. Figure 1 shows the evolution of those indicators along the grid presented in Sect. 3.1⁴ for the set of modes observed in 16 Cygni A. The x-axis is the large separation of spherical degree 0 as we defined above and the y-axis the considered indicator. We also display the observed values for 16 Cygni A (HD 186408) using the frequencies determined by Davies et al. (2015). In blue is the observed value and, in red, the value corrected for the surface effects according to Kjeldsen et al. (2008)’s prescription. We note that these indicators are almost insensitive to surface effects excepted for the case of , which value is changed by about 1σ. We only show the standard deviation for the estimators of the small separation as the one for the large separation is too small to be visible on the plot. Indeed, we computed a standard deviation for Δ₀ of σ(Δ₀) = 5 × 10⁻³ μHz.

Fig. 1. Seismic HR diagram defined with the new indicators and computed along the grid presented in Sect. 3.1. The masses increase from right to left. The blue marker shows the observed value for 16 Cyg A while the red one shows the value corrected for the surface effects following Kjeldsen et al. (2008)’s prescription. The thick line represents the track for 1.06 M⊙. The blue line has been computed for Y0 = 0.27, the orange one for αMLT = 1.5, and the pink one for (Z/X)0 = 0.018. All the coloured lines have been computed for 1.06 M⊙.

To provide a comparison, we computed the evolution of the “usual” indicators along the same tracks as in Fig. 1 and display it in Fig. B.1. We observe that the new indicators exhibit the same behaviour as the usual ones and provide smaller standard deviations, therefore, tighter constraints.

We note on Fig. 1 that is a very good indicator of the core conditions and should hold information about the evolutionary stage on the main sequence as its evolution is almost monotonic. It is therefore very similar to the small separation that has been shown to carry information on the evolution (Christensen-Dalsgaard 1988) as it is sensitive to the sound speed gradient which in turn is sensitive to the chemical composition changes. On the other hand, is not a good indicator of the evolution but carries additional information. For example, de Meulenaer et al. (2010) showed that, for stars with masses and metallicities close to that of α Centauri A (HD 128620), it should provide an upper limit on the amount of convective-core overshooting. We also draw attention to the fact that shows a turn off for evolved stars. This means that we have to be cautious when fitting models to the observations as, for specific sets of input physics, there exists inaccessible regions. This should allow to constrain the input physics.

Finally, does not provide new information. In addition, we may observe from the comparison between the observations and the theoretical tracks that the expected masses retrieved from the different indicators are in agreement. Indeed, from , we should expect masses between 1.06 M_⊙ and 1.11 M_⊙. Then, from , the expected values are in between 1.06 M_⊙ and 1.07 M_⊙. Again, does not add some information as the values range from 1.05 M_⊙ to 1.06 M_⊙. As a consequence, we would expect the mass of 16 Cygni A to be around 1.06 M_⊙. This value has been highlighted on the figure by using a thick line. However, we must observe that this does not provide a precise estimate of 16 Cyg A mass as we only tested a specific chemical composition – Y₀ = 0.25, Z₀ = 0.016 – as well as given choice of α_MLT = 1.82. A proper adjustment is needed to draw conclusions. Those values are only given to illustrate the compatibility between the different indicators. We also show the influence of the composition and α_MLT on the coloured tracks. The blue line is for a higher initial helium abundance. Then, the orange line depicts the influence of a lower value of α_MLT. Finally, the pink line shows how a higher metallicity modifies the results. We observe in all cases that the inferred mass should be lower.

3.2.3. Δ_0l indicators

As it has been shown in several studies, the combination of the small separation ratios r₀₁ and r₁₀ first introduced by Roxburgh & Vorontsov (2003) into r₀₁₀ allows to provide inferences about the stellar central mixed region extension (Popielski & Dziembowski 2005; Deheuvels et al. 2010; Silva Aguirre et al. 2011). Indeed, the mean value and slope of this indicator occupy very specific regions in the parameter space according to the extent of the central mixed region. This should therefore provide constraint on the amount of overshooting necessary to reproduce observations. As an example, an extensive study of several Kepler targets has been realised by Deheuvels et al. (2016) who have been able to provide constraints on the overshooting parameter α_ov for eight of those targets. In the framework of this paper, the indicator represents an estimator of the mean value of r₀₁₀. We may therefore build an indicator for its slope as follows:

(28)

with Δ_l the large separation for modes of spherical degree l defined by Eq. (24). It is straightforward to show that, in the asymptotic regime, Δ₀₁ indeed represents the slope of the frequency ratio r₀₁.

As for the small separation indicators, we computed the evolution of such indicators along the grid of models presented in Sect. 3.1 to demonstrate their regularity and validity. This is shown in Fig. 2. We also display the observed value for the case of 16 Cygni A. This value is corrected for the surface effect using Kjeldsen et al. (2008)’s prescription. As in Fig. 1, the error on Δ₀ is too small to be visible. Moreover, we show the influence of a change in the composition and of α_MLT with the coloured tracks. The colours are the same as in Fig. 1.

Fig. 2. Seismic HR diagram defined with the new indicators Δ0l and computed along the grid presented in Sect. 3.1. The masses increase from right to left. The red marker shows the observed value corrected for the surface effects following Kjeldsen et al. (2008)’s prescription. The thick line represents the track for 1.06 M⊙. The blue line has been computed for Y0 = 0.27, the orange one for αMLT = 1.5, and the pink one for (Z/X)0 = 0.018. All the coloured lines have been computed for 1.06 M⊙.

Again, it is possible to get an estimate of the mass value for the given composition and physics. From Δ₀₁, we expect masses ranging from 1.06 M_⊙ to 1.07 M_⊙. Then, from Δ₀₂, we expect that they lie between 1.05 M_⊙ and 1.07 M_⊙. Finally, from Δ₀₃, the mass should be between 1.00 M_⊙ and 1.05 M_⊙. This time, we observe a slight incompatibility between the first two indicators and the last one. Let us add that we highlighted the value of 1.06 M_⊙ as in Fig. 1. Moreover, we note that, as opposed to the small separation ratio indicators, the relative behaviours of the coloured tracks are different for the three indicators. This could allow to discriminate various choices in the physics of the models considered as they represent different values of α_MLT, Y₀, and Z/X₀ and also to solve the slight mass discrepancy observed.

Finally, as detailed above, the simultaneous use of both and Δ₀₁ allows to provide estimations of the extent of stellar mixed cores. To illustrate this, we plotted main sequence evolutionary tracks for a 1.2 M_⊙, Y₀ = 0.25 and X₀ = 0.734 with several overshooting parameter values ranging from 0.005 to 0.3. The overshooting parameter gives the extent of the mixed core above the Schwarzschild limit through d = α_ovmin(H_p,h) where H_p is the pressure scale height and h the thickness of the convection zone. This is shown in Fig. 3. The ZAMS is at the converging point of the tracks and the TAMS is at their end. We observe that the various tracks occupy very specific regions in the (Δ₀₁,) diagram. Moreover, we note the striking resemblance of Fig. 3 and Fig. 3 of Deheuvels et al. (2016). This should allow us to constrain the amount of overshooting. However, Deheuvels et al. (2016) noted that for this diagnostic tool to be efficient, the mean large separation of the target should not exceed ∼110 μHz and mixed modes should not be present in the oscillation spectrum.

Fig. 3. MS evolutionary tracks in the (Δ01,) plane for models of 1.2 M⊙ with Y0 = 0.25 and X0 = 0.734 for several values of the overshooting parameter αov shown in the legend.

3.2.4.

To provide an estimate of ϵ, we define the following vector sub-space, where frequencies are described as

(29)

where and K are free parameters.

Then, we define an orthonormal basis over this sub-space: and . Finally, by projection of the frequencies over this basis and identification of the several coefficients with the asymptotic formulation, we can retrieve an expression for , the estimator of ϵ. We note that this projection also provides an expression to estimate the large separation that is different from Eq. (25).

Figure 4 shows the evolution of the indicator along the grid presented in Sect. 3.1 but with a step of 0.02 M_⊙. We may observe that this indicator is almost insensitive to the mass for the early stages of the main sequence. The influence of the mass only becomes visible when the stars become older. Moreover, the red marker in Fig. 4 shows the observed value for 16 Cygni A. We note a disagreement between theoretical and observed data. This disagreement can however be tackled by correcting the observed frequencies for the surface effects. This is what the blue and green markers represent. For the blue one, we have computed a correction to the surface effects following Kjeldsen et al. (2008)’s prescription. Then, for the green one, we have computed the correction prescribed by Sonoi et al. (2015). Therefore, it seems reasonable that the indicator we defined could be of some use to constrain the surface effects. It could provide a complementary method to that of Roxburgh (2016). Indeed the method presented uses differences between observed and model ϵ values – under the hypothesis that both the model and observed star have the same inner structure – to isolate only the surface contribution to the measured frequencies. This allows to account for surface effects without the need of empirical corrections. On the other hand, the present indicator should allow to discriminate several surface effects corrections without the need of any physical assumption.

Fig. 4. Evolution of along the grid presented in Sect. 3.1. The masses increase from top to bottom. The step is here of 0.02 M⊙. The red, green, and blue markers respectively represent the observation for 16 Cygni A, the observation for 16 Cygni A corrected for the surface effects using Kjeldsen et al. (2008)’s prescription and the one corrected using Sonoi et al. (2015)’s prescription.

3.3. Glitch indicators

3.3.1. Helium amplitude

With the aim of retrieving the photospheric helium abundance, we built an indicator of the helium glitch amplitude. Verma et al. (2014) obtain their indicator via an integration of the glitch amplitude over the spectrum. We prefer taking advantage of the scalar product and define the indicator as the norm of the helium glitch component, thus generating the following expression:

(30)

Thanks to the orthonormalisation, it is independent of the other indicators and its standard deviation equals 1. We calculated the evolution of this indicator with respect to the surface helium mass fraction Y_f and the surface mass fraction of metals Z_f. This is shown in Figs. 5 and 6. Figure 5 has been computed for stars with a fixed (Z/X)₀ ratio of 0.022 for several surface helium mass fractions, displayed on the abscissa, and for the three values of the mass shown in the legend. To have a reference, we imposed the models to have a fixed value of the large separation – the one observed for 16 Cyg A. In a similar way, Fig. 6 has been computed for a fixed value of Y₀ = 0.24 and for several masses shown in the legend. The helium glitch amplitude has then been computed for the values of Z_f displayed on the abscissa. We insist on the fact that these tracks do not represent the evolution of the helium amplitude with the surface composition along the evolution of a given model. Instead, each point corresponds to a given stellar model that fits the observed 16 Cygni A large separation for a given surface composition. This means that those models were not selected from the grid presented in Sect. 3.1.

Fig. 5. Evolution of the helium glitch amplitude AHe with the surface helium abundance Yf. Each track corresponds to a given mass, written in the legend. The dot-dashed line represents the amplitude for a 1.052 M⊙ model of which the frequencies have been corrected for surface effects as in Kjeldsen et al. (2008). Every model has an initial heavy elements abundance of (Z/X)0 = 0.022. Each point has been computed with the same large separation to remain at the same evolutionary stage.

Fig. 6. Evolution of the helium glitch amplitude AHe with the surface heavy elements abundance Zf. Each track corresponds to a given mass labeled in the legend.

We observe in Fig. 5 an increasing trend in the helium glitch amplitude with the helium mass fraction as well as with the mass. This has to be expected as a larger quantity of helium inside the star would lead to a more important depression of the first adiabatic index Γ₁ at the second ionisation zone of helium, and, therefore, a glitch of greater amplitude. Moreover, we also show the influence of the surface effects by computing the amplitude evolution for a 1.052 M_⊙ star for which the surface effects have been taken into account via Kjeldsen et al. (2008)’s prescription (dot-dashed line). It is apparent that they have little influence on the amplitude as the values remain in the 1σ error bars of the uncorrected models. This was expected as the glitch is of greater amplitude in the low frequencies regime while the surface effects corrections are greater in the high frequencies regime.

Furthermore, as shown in Fig. 6, the glitch amplitude and the metallicity are anti-correlated. This corroborates Basu et al. (2004)’s observations. Therefore, we are facing a degeneracy and the glitch amplitude alone will not be sufficient to estimate properly the surface helium mass fraction. Thus, the smooth component indicators defined above will be of great help. This clearly shows that the A_He–Y_f relation is model dependent which should never be forgotten.

Figure 7 illustrates this degeneracy. It represents the profile of Γ₁ as a function of the reduced radius in the superficial layers of stars of a fixed large separation but with several chemical compositions. We immediately notice that both an increase of the surface helium abundance and a decrease of the surface heavy elements abundance lead to a minimum that increases in magnitude. Therefore, the helium glitch amplitude becomes greater as well. We provide an interpretation of this phenomenon in Sect. 4.2.

Fig. 7. Evolution of Γ1 as a function of the reduced radius in the superficial layers of stars of fixed large separation. Every star has a mass of 1.052 M⊙. Two families of curves are displayed. The black ones have a fixed value of the ratio between the initial hydrogen and heavy elements abundances of (Z/X)0 = 0.022 and variable initial helium abundance. The red ones have a fixed initial helium abundance of Y0 = 0.24 and a variable initial heavy elements abundance. The different values are displayed on the figure.

3.3.2. Convection zone amplitude

The definition of the envelope convective zone glitch amplitude we provide is very similar to that of the helium glitch and is the following:

(31)

We expect this indicator to be a proxy of the sharpness of the transition between the envelope convective zone and the radiative zone. Again, thanks to the orthonormalisation, it is independent of the other indicators and its standard deviation is equal to 1. We present in Fig. 8 its evolution with the importance of the undershooting, characterised by the coefficient α_under. This coefficient determines the size d of the undershooting region at the bottom of the convective envelope. This size is given by d = α_undermin(H_p,h) where H_p is the pressure scale height and h the thickness of the convection zone. In the undershoot region, the temperature gradient is set to the adiabatic one and the mixing is assumed to be instantaneous.

Fig. 8. Evolution of the convection zone glitch amplitude ACZ with αunder for a star of mass 1.052 M⊙, X0 = 0.744, Z0 = 0.016 and at a constant value of the large separation.

We expect an increase of A_CZ with α_under as the introduction of undershooting in a stellar model will create a discontinuity of the temperature gradient. The temperature gradient jump increases with the value of α_under. Thus, the glitch amplitude increases as well. This is what we show in Fig. 8 where the computed models have a mass of 1.052 M_⊙, an initial hydrogen mass fraction of X₀ = 0.744, and an initial metal mass fraction of Z₀ = 0.016. As for the case of the helium amplitude evolution, we kept a fixed large separation – which is that of 16 Cyg A – for each model. Moreover, we noted that each computed model was at the same evolutionary stage (constant central hydrogen mass fraction). Therefore, the observed effect is not evolutionary but rather the effect of the temperature gradient discontinuity as expected. Finally, we may add that, when setting the temperature gradient to the radiative one in the undershoot region, we do not observe any significant trend in the amplitude with α_under. We may thus conclude that we observe the effect of the temperature gradient and not of the chemical composition.

4. Method characterisation

4.1. Capabilities

For the observed data, we have chosen the frequencies computed by Davies et al. (2015) for the component A of the binary system 16 Cygni (HD 186408). We come back to the particular case of 16 Cyg A in Sect. 5. Figure 9 shows the difference between the observed frequencies and the smooth component of the fitted frequencies, the fitted helium glitch alone and the fitted helium and convection zone glitches. Only the l = 0 fitted curves are displayed. We observe that the fit is good and that the helium glitch has been properly isolated. However, the convection zone glitch is of very low amplitude, compared to the helium glitch, and has a negligible contribution. This is visible in the negligible improvement of the χ² value from the results without including the convective zone glitch (about 10% variation). These results are similar to those of Verma et al. (2014).

Fig. 9. Fitted glitch to 16CygA (HD 186408) data (Davies et al. 2015). Only the l = 0 fitted curves are displayed.

4.2. Limitations

The presented method has been developed for the study of solar-like pulsators and to provide a comprehensive analysis of their oscillation spectra. However, it is not yet adapted to study evolved stars which exhibit mixed modes.

In addition, from masses around 1.25 M_⊙ and above as well as for the highest values of the helium abundances considered (from Y_f ∼ 0.195 and above), the evolution of the helium glitch amplitude with the surface helium mass fraction is less monotonic and inferences become unreliable. Figure 10 shows this limitation. Indeed, for 1.25 M_⊙ we observe a very sharp increase of the helium amplitude as we defined in Sect. 3.3.1 for a quasi constant surface helium mass fraction. This corresponds to a decrease of the surface Z/X ratio which is shown in Sect. 3.3.1 to lead to a higher amplitude.

Fig. 10. Evolution of the helium glitch amplitude AHe with the surface helium abundance Yf for high masses.

Also, we observe for 1.3 M_⊙ that the last point (Y_f ∼ 0.197 and A_He ∼ 96) moves backwards. This is due to microscopic diffusion. We observe that the size of the convective envelope decreases with the initial helium abundance. This leads to a more efficient gravitational settling as the diffusion velocities are greater close to the surface. Therefore, the surface abundance becomes smaller than for the previous point that had a lower initial helium abundance. We observe the same decrease in the surface metallicity which explains the increase of the amplitude even though the surface helium mass fraction remained constant.

We also noted in Sect. 3.3.1 that the surface abundance of metals has a relevant influence on the helium glitch amplitude and separating its contribution from that of the helium is not an easy task. To investigate such a behaviour, we developed a toy model for the first adiabatic index Γ₁ which is thought to be the main contributor to the helium glitch amplitude (Gough 1990; Houdek & Gough 2007). In this model, we trace back the influence of the helium and metals abundances on the dip of Γ₁ in the helium second ionisation zone. More information about the construction of this model is given in Appendix C. Using it, we were able to test the influence of the chemical composition decoupled from its evolutionary effect on the temperature and density profiles. To do so, we artificially modified the chemical composition profile of a reference model (black curve in Fig. 11), without changing its temperature and density profiles, to match the surface abundance, in either metals or helium, of a second reference model. This way, we were able to isolate the contribution of the chemical composition alone. Figure 11 shows the comparison between the helium second ionisation zone toy models with modified chemical composition profiles only (blue curves) and toy models with the same composition but for which its effect on temperature and density have been taken into account (red curves). We observe in the top panel that the effect of the helium abundance dominates over the effect of temperature. In the bottom panel, we observe that the metal abundance alone does not modify by a significant amount the dip in the helium second ionisation zone. However, when taking its influence on the temperature and density profiles into account, the effect becomes significant. This allows to better understand the degeneracy in composition on the helium amplitude. One effect, that of the helium abundance is direct on the shape of the Γ₁ profile in the helium second ionisition zone while the second, that of the metal abundances, is indirect as it influences the temperature and density profiles which in turn modify the Γ₁ profile.

Fig. 11. Comparison with a reference model (black curve) of the toy model Γ1 profiles in cases for which the temperature and density profiles have been decoupled from the composition profiles (blue curves) and coupled cases (red curves). The dashed lines have a common decreased abundance of the considered element and the dotted lines have an increased abundance. Top panel: effect of the helium abundance; bottom panel: effect of the metals abundance.

Finally, we observe that, in a lnP − lnT diagram, the curves for the different models in the temperature region of the Γ₁ dip are parallel to each other but at various height. A higher lnP − lnT curve corresponds to a shallower Γ₁ dip. We also observe this behaviour with a fixed composition and a variable mass. Higher masses models have lower curves and deeper depressions in Γ₁. Therefore, the toy model allows to understand the influence of the mass on the helium glitch amplitude as well.

5. Illustration with 16 Cygni A (HD 186408) observations

16 Cyg A (HD 186408) is one of the brightest stars in the Kepler field of view. It belongs to a binary system of solar analogs, both exhibiting solar-like pulsations. The quality and length of the collected time series make it the ideal subject to test the method. It should be noted that we do not provide a detailed study of 16 Cygni A. Indeed, we only present here the capability of the method to provide structural constraints. A detailed study will be the object of a future paper.

5.1. Methodology

To obtain constraints on 16 Cyg A, we compute the value of the seismic indicators using the frequencies determined by Davies et al. (2015). We have corrected the surface effects for the observed frequencies by using the power law prescribed by Kjeldsen et al. (2008) and the a and b coefficients fitted by Sonoi et al. (2015) as a function of T_eff and g. The authors have done this coefficient adjustment by comparing the adiabatic frequencies of patched models based on 3D simulations and that of unpatched standard 1D models. We then fit the observed values of Δ, , , and A_He with the age, mass, initial mass fraction of hydrogen, and the initial heavy elements over hydrogen abundance ratio as free parameters. To do so, we select an initial guess value of X₀ and derive the best fit values of Δ, , and . This results in a set of values for the mass, the age, and (Z/X)₀. Then, the adjustment is done by applying the secant method to find the value of X₀ in best agreement with the target A_He. At each step of the secant algorithm – at each value of X₀ –, Levenberg-Marqardt’s algorithm (L-M) computes the optimal set of stellar parameters giving the best fit of Δ, , and . Finally, when we have a good estimate for A_He, we use L-M’s algorithm one last time to derive the complete set of values for the mass, the age, X₀, and (Z/X)₀ fitting the observed parameters. We used the secant method in order to diminish the computational time needed to converge towards the solution.

5.2. Results

The values of the relevant seismic indicators as well as the associated standard deviations computed from the observed frequencies are given in Table 1. The column labeled Kjeldsen refers to the observed indicators corrected for the surface effects using Sonoi et al. (2015)’s coefficient fitted to Kjeldsen et al. (2008)’s prescription as explained in Sect. 3.1. The Sonoi column corresponds to the correction of surface effects using a Lorentzian profile as in Sonoi et al. (2015). We also derived the seismic indicators while adding smooth basis elements with n⁻¹. As the frequencies are projected in a specific order to build the seismic indicators, only the glitches amplitudes are affected and we have: A_He = 27.6 and A_CZ = 3.3. Using the methodology described above, we managed to derive the stellar parameters given in Table 2 for the different sets of observed seismic indicators. We note that changing the treatment of surface effects or adding the n⁻¹ basis elements have an impact on the fitted parameters comparable to that of the frequencies uncertainties. Thus, we will only discuss here the parameters adjusted to the data corrected for the surface effects using Kjeldsen et al. (2008)’s prescription and without adding n⁻¹ basis elements. The parameters we derived do not constitute a detailed characterisation of the target 16 Cyg A. Rather, they illustrate the ability of the method to provide constraints on a solar analog. Indeed, only one single set of input physics was tested. Therefore, the standard deviations tend to be underestimated as they are the ones intrinsic to the method. The abundances used for the computations were the solar ones determined by Asplund et al. (2009). Let us add that we obtain a surface helium abundance of Y_f = 0.242 ± 0.028 which lies in the interval obtained by Verma et al. (2014), Y_f, V ∈ [0.231,0.251]. This is comforting us in the idea that the developed method is efficient in isolating the glitches and drawing inferences from their signatures. In Fig. 12, we show the evolution of the helium glitch amplitude resulting from the three parameters adjustment of Δ, , and as a function of X₀. We also show the evolution of the value of (Z/X)₀. The very linear trend justifies that we used the secant method to provide successive estimates of A_He in order to lessen the computational charge. In addition, we illustrate both the observed value for A_He and the corrected value under Kjeldsen et al. (2008)’s prescription. We note that the corrected value is of about 0.79σ lower than the uncorrected one. Using Sonoi et al. (2015)’s prescription only leads to a 0.55σ

Fig. 12. Helium amplitude (solid line) and (Z/X)0 ratio (dotted line) versus the initial hydrogen abundance. The horizontal solid line represents the target value for the amplitude – corrected for the surface effects following Kjeldsen et al. (2008) – while the dashed lines represent the 1σ interval. The horizontal dotted line is the value of the amplitude without surface effects correction.

variation of the measured amplitude. This demonstrates that the surface effects have a small influence on the amplitude we derive using the method. Also, we retrieved a value of 970.97 s for τ_He and 3042.32 s for τ_CZ. As they were fixed to model values and, as a consequence, were neither parameters nor constraints of the adjustment, we do not provide uncertainties. However, we may compare their values to the ones adjusted by Verma et al. (2014). They obtained τ_He ∈ [868,944] s and τ_CZ ∈ [2992,3234] s. We observe that τ_CZ lies in the interval calculated by Verma et al. (2014) while the value of τ_He is slightly above the upper limit. However, as we show in Sect. 2.3 this does not impact the inferences drawn from the helium glitch amplitude in a significant way. To illustrate this statement, we freed the value of τ_He (but not that of τ_CZ as the convection zone glitch is of negligible amplitude compared to that of the helium glitch). The relative change between the optimised τ_He and its estimator is of only 6%. Moreover, the observed value of the helium glitch amplitude remains unchanged compared to its standard deviation (we observe a change of 0.06σ). As expected, the best fit Y_f also remains untouched. At this point, one should be reminded that we focus on the glitch amplitude to draw our inferences. And, as we showed in Sect. 2.3, the exact location of the glitch does not have a significant impact on its amplitude as well as on the derived surface helium abundance. Finally, Fig. 13 shows the fitted helium glitches⁵ for the observed data and the 4 parameters best fit model. It is visible that the observed and fitted glitches are close in amplitude and period. This demonstrates the ability of the method to both isolate the glitch and the related parameters and to provide a model reproducing at best those parameters. We also show in Fig. 14 the comparison in an échelle diagram of the observed frequencies and those of the best fit model. We observe in this last figure that the smooth component of the spectrum (not visible in the previous figure) is properly adjusted. Thus, by fitting the set of indicators Δ₀, , , and A_He we obtain a good representation of the observed frequencies. We also show in Appendix F the adjustment to frequencies corrected for surface effects using a Lorentzian profile as in Sonoi et al. (2015). We observe in Fig. F.1 that the agreement between the model and observed glitches is better than in Fig. 13. However, we also note in the corresponding échelle diagram (Fig. F.2) that there is an offset between observed and model ridges. This corresponds to different values of . Figure 4 also illustrates this behaviour as it is clearly visible that only Kjeldsen et al. (2008)’s prescription allow to reproduce the observed value for .

Fig. 13. Comparison between the observed helium glitch δνo (solid line) and the one resulting from the best fit model δνm (dot-dashed line), for l = 0. We also display the observed glitch as a function of the frequencies (errorbars) as well as the best model glitch associated with the theoretical frequencies (diamond).

Fig. 14. Comparison between the observed frequencies (diamonds) and the best model frequencies (circles) in an échelle diagram.

6. Discussion and conclusions

6.1. Principle

In the present paper, we provide a new method that uses as much of the available seismic information as possible. Indeed, we take advantage of the information contained in both the glitches and the smooth component (usually discarded in glitches analyses but then used separately for forward seismic modelling, see Verma et al. 2014 for example) of the spectrum to define indicators that are as independent as possible of each other. To do so, we take advantage of Gram-Schmidt’s algorithm to create an orthonormal basis over which we project the frequencies. The obtained coefficients are therefore independent of each other and consist in a linear combination of the frequencies. Thus, using the appropriate combination of those allows us to define uncorrelated indicators. Such indicators are constructed in order to reproduce the behaviour of “usual” indicators such as the large frequency separation. Up to this day and to our knowledge, no method, has been proposed to provide proper correlations between the smooth and glitch indicators as they are built separately. Therefore, our method provides the asteroseismologists with new diagnosis means.

6.2. Advantages

As the method only relies on linear algebra, it is very stable and the computation times are negligible – of the order of a fraction of a second –. Thus, it could easily be implemented in stellar model fitting algorithms, which are non-linear, without impacting the total computational time. This can be done with any algorithm as the method only focuses on the definition of new seismic indicators and does not rely on the physics of the model itself. Let us add that the defined indicators should be used in combination with non-seismic constraints through a single merit function while searching for a stellar model in order to obtain proper covariances between the inferred quantities – which is not often the case in seismic analyses. The illustration for 16 Cyg A truly demonstrates the possibility to use the new indicators to provide further constraints on stellar structure in the framework of forward seismic modelling.

In addition, the usual indicators often hold a local definition and may use correlated information while the new ones are built in a way that the information used is averaged over the whole frequency range and that it is not redundant (each observed frequency is used only once for each indicator). This results in smaller standard deviations and smoother behaviours. Appendix B shows the evolution of the classical indicators r₀₁(n) and r₀₂(n) as defined by Roxburgh & Vorontsov (2003). We indeed observe that the new definitions of the indicators give smaller error bars while preserving the expected trends, providing tighter constraints. Finally, the method has the advantage that it can be implemented even in cases where some modes are missing. It is not the case of the “usual” methods as they need successive modes to define some of their indicators. For example, the classical local definitions of the large separation and second differences require at least two consecutive frequencies. Therefore, whenever some modes are missing, pieces of information might be discarded.

6.3. Information carried by the indicators

We have defined indicators that provide estimates for the classical indicators that are the large separation, Δ and Δ_l, and the small separation ratios, . As expected, Δ provides an estimation of the stellar mean density (Ulrich 1986). Also, we defined the indicators Δ_0l which combine the large separations associated with the spherical degrees 0 and l which are known to give an estimate of the resonant cavity of the l degree modes (Monteiro 2002). Those indicators provide an estimator of the slope of r₀₁₀ which, combined with can be used to constrain the overshooting parameter as shown in Fig. 3.

Moreover, we observe in Fig. 1 that the indicator , defined to estimate the small separation ratio between the spherical degrees 0 and 2, is a proper indicator of the evolution of the star, as expected from Christensen-Dalsgaard (1988). Also, we observe that the indicator presents a degeneracy, as a turn-off occurs in the Δ₀– plane, and there is an inaccessible region. This is a very interesting observation as it should provide tight constraints on the stellar structure. Indeed, we show that a change in composition or in the value of α_MLT allows to modify that region. It will therefore be necessary to use such parameters to reproduce observed values.

Furthermore, Fig. 4 shows that we might get constraints on the surface effects from the indicator we defined. It should allow to discriminate from several empirical formulations meant to account for the surface effects. We have tested both formulations from Kjeldsen et al. (2008) and Sonoi et al. (2015) as an illustration of the diagnosis power of the indicator . However, we are aware of the existence of the formulation from Ball & Gizon (2014) and it should also be inspected in further studies. Besides, we showed that the helium amplitude indicator is almost unaffected by the surface effects as it has to be expected (see for example Fig. 12). Indeed, empirical surface effects corrections are important for the high frequencies compared to ν_max while the glitch is of great amplitude only for the low frequencies. Also, the helium glitch amplitude should allow to draw inferences on the surface helium abundance as shown by Fig. 5. However, attention has to be paid as it is also anti-correlated with the metallicity (see Fig. 6). We have demonstrated via a toy model for the first adiabatic index that both effects on the amplitude stem from the position of the adiabat which in turn determines the amplitude of the second local Γ₁ minimum due to helium partial ionisation.

In addition, we observe that the convective zone glitch amplitude has a significantly lower amplitude than that of the helium glitch and is correlated with the amount of undershooting at the base of the envelope convection zone (in agreement with Verma et al. 2014).

Finally, let us add that it is possible to define other indicators than those presented in this paper. This should therefore be carefully studied to take advantage of as much of the available information as possible.

6.4. Limitations

However, we show that the method is only fit to draw inferences about solar-like stars, that is low-mass stars on the main-sequence. Indeed, Fig. 10 illustrates that from masses around 1.3 M_⊙ and above, diffusion plays an important role and the relation between helium surface abundance and the helium glitch amplitude is not monotonic anymore. Therefore, the method will have to be adapted for massive and evolved stars.

Moreover, one could argue that using model values as estimators of the acoustic depths of the glitches is a major drawback of the method. However, the proposed method does not focus on those quantities. Rather, we focus on the information that the amplitude of the glitch (not the period of the signal) carries. This means that we only need proper, but not exact, estimators for the acoustic depths in order to draw inferences. Indeed, we showed in Sect. 2.3, in the case of the helium glitch, that a small excursion from the estimated value (either by manually setting another value, that of the second minimum in the first adiabatic index for the helium glitch, or by optimising over its value) does not lead to a significant change in the measured amplitude. Therefore, the inferences drawn remained unchanged. However, one could still regard this as a flaw of the method as we do not provide a new way of retrieving the acoustic depths of glitches. This could be explored in future studies by, for example, finding the global opitmum for the acoustic depths of the glitches and then using the measured value as a constraint for the best fit model. Nevertheless, this would make the calculations more time consuming and annihilate the benefit of the orthonormalisation that is the independence of the fitted parameters.

Also, the indicators should be used to complement non-seismic data as we have shown that, for example, there exists a degeneracy between the helium glitch amplitude, the helium surface abundance and the metallicity⁶ (see Figs. 5 and 6). Therefore, the helium amplitude alone is not sufficient and additional information are needed to lift such a degeneracy – these information may be contained in other indicators (seismic and non-seismic).

6.5. Future perspective

The next step will be the detailed study with the new method of several stars from the Kepler legacy sample (Lund et al. 2017). This sample consists of 66 main sequence stars for which at least one year of continuous observations has been made. Having such long time series provides the necessary precision to study glitches. Let us add that this future study is intended to be independent of what has been done in Verma et al. (2017).

Also, another important step will be the improvement of the method by enabling the study of glitches present in red subgiants spectra. The peculiarity of such pulsators is that they exhibit mixed modes. This will consist in a new challenge as their analytical formulation is not a simple task. However, this is a necessary step to improve our knowledge of the evolution of a star such as the Sun. Moreover, this will allow a better description of such stars and will provide a deeper understanding of their properties via, for example, the characterisation of the mixing processes.

Acknowledgments

The authors would like to thank the referee for their very constructive remarks that improved the paper in a significant way. M.F. is supported by the FRIA (Fond pour la Recherche en Industrie et Agriculture) – FNRS PhD grant. S.J.A.J.S. is funded by ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. GB acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, G.A. no. 772293).

References

Appendix A: Effect of the regularisation constant

Fig. A.1. Comparison between the glitch adjustment for several values of the regularisation constant λ for 16 Cyg A best model determined in Sect. 5.2. The values are given in the legend. The lines correspond to the fitted glitch and the markers to the observed glitch. Only the fitted curves for l = 0 are displayed.

Figure A.1 shows the influence of the regularisation constant λ on the quality of the fitted glitch. We fitted the glitch of the best model obtained in Sect. 5.2 using the values of λ shown in the legend. We observe that, when using a value of 10² which is already significant as discussed in Sect. 2.2, the fit (red dashed curve) remains close to the one obtained without including regularisation terms (solid black curve). There is only a slightly higher dispersion in the model frequencies subtracted from the smooth part of the spectrum. To show the degradation caused by a regularisation constant that dominates the adjustment, we used λ = 10³. We observe both a discrepancy between results without and with (dotted black curve) regularisation terms and a higher dispersion. This translates in a higher value of the merit function (from χ² ∼ 1 without regularisation to χ² ∼ 500 with λ = 10³).

Appendix B: Comparison with usual indicators

We present in Fig. B.1 the evolution of the “usual” indicators r₀₁(n) and r₀₂(n) (Roxburgh & Vorontsov 2003) used in asteroseismology in order to compare them with the new indicators. They are evaluated at the value of n = 21 which corresponds to the measured value of n_max for 16 Cyg A (n_max being the value of n at l = 0 of closest frequency to the ν_max value). First let us draw the attention to the fact that the usual indicators hold a local value, as Eqs. (26) and (27) show, while the new indicators are averaged over the complete set of available modes. It results from this averaging a lower standard deviation. Then, we notice that the behaviour of both the usual indicators and the new ones follow similar trends. This means that the definition is consistent with what has been done up to now and that it should hold the same information. Therefore, it should be able to provide similar diagnostics, but with higher precisions.

Fig. B.1. Evolution of the normalised small separation as defined by Roxburgh & Vorontsov (2003) and evaluated at n = 21 along the grid presented in Sect. 3.1. The red marker shows the observed value corrected for the surface effects following Kjeldsen et al. (2008)’s prescription. The thick line represents the track for 1.06 M⊙. Top panel: three points normalised small separation between spherical degrees 0 and 1; middle panel: five points normalised small separation between spherical degrees 0 and 1; bottom panel: normalised small separation between spherical degrees 0 and 2.

Appendix C: Γ₁ toy model

To build a toy model for Γ₁ that replicates at best its behaviour in the helium second ionisation zone we use the following hypotheses:

• In the helium second ionisation zone, hydrogen is fully ionised;

• Metals are in their atomic form;

• We consider a perfect gas.

We then define:

• The once ionised helium number of particles per unit volume He⁺ and twice ionised helium number of particles per unit volume He⁺⁺ such that:

  (C.1)

  where m_u is the atomic mass unit;

• The helium ionisation fraction:

  (C.2)

  which equals 0 if the helium is ionised only once and 1 if it is fully ionised;

• The number of electrons per unit volume:

  (C.3)

• The total number of particles:

  (C.4)

  where Z_i is the mass fraction of the metal labeled i and A_i its mass number.

Using Saha’s equation (Saha 1920):

(C.5)

where g is the statistical weight of helium at its fundamental state, h is Planck’s constant, k_B is Boltzmann’s constant, m_e the electron mass, and χ the helium second ionisation energy, we may obtain a second order equation for x_He:

(C.6)

which we solve to obtain the evolution of x_He with temperature and density. We may also derive this expression with respect to temperature at constant density. This yields

(C.7)

Then we have the following expressions:

(C.8)

(C.9)

(C.10)

where μ is the mean molecular weight and

(C.11)

Finally, we may insert their values in the relation linking Γ₁, P_T, P_ρ, and c_ν:

(C.12)

Appendix D: Gram-Schmidt’s process and QR decomposition

As a reminder, Gram-Schmidt’s algorithm consists in the construction of orthonormal basis elements from a set of non-orthonormal basis elements. Let us consider the element of index j₀ : p_j0. Let us also assume that we have already built the set of orthonormal basis functions up to index j₀ − 1, that is the set (q₁,…,q_j0 − 1). To build element q_j0, we first subtract to p_j0 its projection over the successive previous basis elements. We thus have

(D.1)

where u_j0 is the basis element orthogonal to the set (q₁,…,q_j0 − 1). Finally, it is normalised to obtain

(D.2)

It is also possible to express this process as a QR decomposition. To do so we call 𝒫_l the matrix of initial polynomials and 𝒬_l the matrix of orthonormal polynomials for a given spherical degree:

(D.3)

(D.4)

where n_min is the lowest observed radial order and n_max the highest one for the spherical degree considered.

We may then express Gram-Schmidt’s procedure in a matrix form as a QR decomposition as follows:

(D.5)

where, using our definition of the scalar product

(D.6)

is an upper triangular matrix.

We may now generalise by considering several spherical degrees. For the smooth part, the matrices become block matrices which are the combination of the matrices for the several spherical degrees considered. This block disposition illustrates the fact that we have different basis functions for the different spherical degrees – this is not the case for the glitch part –. This is represented by a Kronecker delta in Eq. (6). In these matrices, each row corresponds to a given mode and each column corresponds to a given element of the basis. Finally, we append the glitch part to those block matrices. They take the following form (we remind that the independent variable for the glitch basis is ):

(D.7)

(D.8)

And the R matrix becomes

(D.9)

where ⟨q_k|p_g⟩_l = 0 is a column vector of which the rows are the successive scalar products of the basis elements with the glitch function – denoted with the g index – such that ⟨q_k|p_g⟩_l = 0 = (⟨q₀|p_g⟩_l = 0⋯⟨q₂|p_g⟩_l = 0)^T. We may again write the QR decomposition as

(D.10)

Appendix E: A numerical example

We generated a set of frequencies for spherical degrees from 0 to 2 (listed in Table E.1) and applied our method. We show in Figs. E.1 and E.2 the successive adjustments of the basis functions. This shows the validity of using the set of functions described in Sect. 2 and allows oneself to compare their results with ours. As a reminder, we project the frequencies over the basis in the specific order explained in Sect. 2. Therefore, taken in the correct order, these plots provide intermediary results for the glitch adjustment. We also show in Table E.2 the values of the fitted coefficients. In this example, we also fitted n⁻¹ polynomials – inspired by the second order form of the asymptotic expansion – to show that it is not necessary, nor relevant, to add supplementary basis elements to our method. Indeed, the fitted coefficient values become comparable to the standard deviation that is equal to 1, through the orthonormalisation. Such values and plots were obtained as follows:

1. Considering the set of standard deviations from Table E.1, we use Gram-Schmidt procedure (Eq. (4)) associated with the definition of the scalar product (Eq. (1)) to produce the ordered orthonormal basis functions for each value of the spherical degree l. Thus, we successively project the former basis elements p_j(n, l) (i.e. the ordered set of polynomials n⁰,n¹,n² and the glitch functions for each spherical degree) on the already defined orthonormal basis elements q_j0(n, l). Then, we normalise those projections. This provides us with the orthonormal basis elements q_j0(n, l). as well as the transformation matrix .

2. For the smooth part and one spherical degree at a time, the frequencies from Table E.1 are projected on the orthonormal basis elements following the proper order to produce the fitted frequencies, ν_f(n, l), according to Eq. (5) where the fitted coefficients are given by Eq. (8). Then, for the glitch part, the frequencies are projected simultaneously for every spherical degree on the glitch basis elements following the same procedure. This is due to the fact that the glitch coefficients should not depend on l. This produces the coefficients from Table E.2.

Fig. E.1. Comparison between the successive adjustments and the observed radial modes frequencies listed in Table E.1. Upper panel: observed frequencies compared to the first order adjustment; lower panel: residual of the first order adjustment compared to the second order fit to those residuals.

Fig. E.2. Separated glitches adjustments to the frequencies in Table E.1 of both glithes shown for the spherical degree l = 0. Upper panel: helium glitch; lower panel: convection zone glitch.

Appendix F: Supplementary adjustments

Fig. F.1. Comparison between the observed helium glitch δνo (solid line) and the one resulting from the best fit model δνm (dot-dashed line) for l = 0. We also display the observed glitch as a function of the frequencies (errorbars) as well as the best model glitch associated with the theoretical frequencies (diamond). The observed frequencies have been corrected for surface effects using Sonoi et al. (2015)’s prescription.

Fig. F.2. Comparison between the observed frequencies (diamonds) and the best model frequencies (circles) in an échelle diagram. The observed frequencies have been corrected for surface effects using Sonoi et al. (2015)’s prescription.

Fig. F.3. Comparison between the observed helium glitch δνo (solid line) and the one resulting from the best fit model δνm (dot-dashed line) for l = 0. We also display the observed glitch as a function of the frequencies (errorbars) as well as the best model glitch associated with the theoretical frequencies (diamond). We include polynomials in n−1 to the basis functions.

We show in the present section the best fit model frequencies adjusted to the observed 16 Cygni A frequencies for several cases. First we show the adjustment to the frequencies corrected for surface effects using a Lorentzian profile as in Sonoi et al. (2015). Figure F.1 shows that both glitches are in good agreement. However, the frequencies are systematically shifted, as Fig. F.2 illustrates, as a consequence of the difference between the observed and theoretical values. Figure 4 illustrates such a discrepancy. This shows that, even though Sonoi et al. (2015) showed that for high frequency regimes Kjeldsen et al. (2008)’s prescription is not able to reproduce frequency differences between patched and unpatched model, it is the only tested empirical correction that allowed us to reproduce the observed value for . In a further study, it would be appropriate to try out a scaled formulation of Ball & Gizon (2014)’s correction such as presented in Manchon et al. (2018).

We also provide in Fig. F.3 the best fit model including terms in n⁻¹ in the basis functions for the smooth part of the spectrum. By eye, the glitch adjustment seems better than in Fig. 13. However, as shown in Table E.2, the improvement is not significant as the fitted parameters values are comparable to their standard deviation. Moreover, Table 2 demonstrates that the effect of including such terms in the adjustment is comparable to a variation of 1σ in the frequencies.

All Tables

All Figures

Fig. 1. Seismic HR diagram defined with the new indicators and computed along the grid presented in Sect. 3.1. The masses increase from right to left. The blue marker shows the observed value for 16 Cyg A while the red one shows the value corrected for the surface effects following Kjeldsen et al. (2008)’s prescription. The thick line represents the track for 1.06 M⊙. The blue line has been computed for Y0 = 0.27, the orange one for αMLT = 1.5, and the pink one for (Z/X)0 = 0.018. All the coloured lines have been computed for 1.06 M⊙.

In the text

Fig. 2. Seismic HR diagram defined with the new indicators Δ0l and computed along the grid presented in Sect. 3.1. The masses increase from right to left. The red marker shows the observed value corrected for the surface effects following Kjeldsen et al. (2008)’s prescription. The thick line represents the track for 1.06 M⊙. The blue line has been computed for Y0 = 0.27, the orange one for αMLT = 1.5, and the pink one for (Z/X)0 = 0.018. All the coloured lines have been computed for 1.06 M⊙.

In the text

	Fig. 3. MS evolutionary tracks in the (Δ01,) plane for models of 1.2 M⊙ with Y0 = 0.25 and X0 = 0.734 for several values of the overshooting parameter αov shown in the legend.
In the text

	Fig. 4. Evolution of along the grid presented in Sect. 3.1. The masses increase from top to bottom. The step is here of 0.02 M⊙. The red, green, and blue markers respectively represent the observation for 16 Cygni A, the observation for 16 Cygni A corrected for the surface effects using Kjeldsen et al. (2008)’s prescription and the one corrected using Sonoi et al. (2015)’s prescription.
In the text

Fig. 5. Evolution of the helium glitch amplitude AHe with the surface helium abundance Yf. Each track corresponds to a given mass, written in the legend. The dot-dashed line represents the amplitude for a 1.052 M⊙ model of which the frequencies have been corrected for surface effects as in Kjeldsen et al. (2008). Every model has an initial heavy elements abundance of (Z/X)0 = 0.022. Each point has been computed with the same large separation to remain at the same evolutionary stage.

In the text

	Fig. 6. Evolution of the helium glitch amplitude AHe with the surface heavy elements abundance Zf. Each track corresponds to a given mass labeled in the legend.
In the text

Fig. 7. Evolution of Γ1 as a function of the reduced radius in the superficial layers of stars of fixed large separation. Every star has a mass of 1.052 M⊙. Two families of curves are displayed. The black ones have a fixed value of the ratio between the initial hydrogen and heavy elements abundances of (Z/X)0 = 0.022 and variable initial helium abundance. The red ones have a fixed initial helium abundance of Y0 = 0.24 and a variable initial heavy elements abundance. The different values are displayed on the figure.

In the text

	Fig. 8. Evolution of the convection zone glitch amplitude ACZ with αunder for a star of mass 1.052 M⊙, X0 = 0.744, Z0 = 0.016 and at a constant value of the large separation.
In the text

	Fig. 9. Fitted glitch to 16CygA (HD 186408) data (Davies et al. 2015). Only the l = 0 fitted curves are displayed.
In the text

	Fig. 10. Evolution of the helium glitch amplitude AHe with the surface helium abundance Yf for high masses.
In the text

Fig. 11. Comparison with a reference model (black curve) of the toy model Γ1 profiles in cases for which the temperature and density profiles have been decoupled from the composition profiles (blue curves) and coupled cases (red curves). The dashed lines have a common decreased abundance of the considered element and the dotted lines have an increased abundance. Top panel: effect of the helium abundance; bottom panel: effect of the metals abundance.

In the text

	Fig. 12. Helium amplitude (solid line) and (Z/X)0 ratio (dotted line) versus the initial hydrogen abundance. The horizontal solid line represents the target value for the amplitude – corrected for the surface effects following Kjeldsen et al. (2008) – while the dashed lines represent the 1σ interval. The horizontal dotted line is the value of the amplitude without surface effects correction.
In the text

	Fig. 13. Comparison between the observed helium glitch δνo (solid line) and the one resulting from the best fit model δνm (dot-dashed line), for l = 0. We also display the observed glitch as a function of the frequencies (errorbars) as well as the best model glitch associated with the theoretical frequencies (diamond).
In the text

	Fig. 14. Comparison between the observed frequencies (diamonds) and the best model frequencies (circles) in an échelle diagram.
In the text

	Fig. A.1. Comparison between the glitch adjustment for several values of the regularisation constant λ for 16 Cyg A best model determined in Sect. 5.2. The values are given in the legend. The lines correspond to the fitted glitch and the markers to the observed glitch. Only the fitted curves for l = 0 are displayed.
In the text

Fig. B.1. Evolution of the normalised small separation as defined by Roxburgh & Vorontsov (2003) and evaluated at n = 21 along the grid presented in Sect. 3.1. The red marker shows the observed value corrected for the surface effects following Kjeldsen et al. (2008)’s prescription. The thick line represents the track for 1.06 M⊙. Top panel: three points normalised small separation between spherical degrees 0 and 1; middle panel: five points normalised small separation between spherical degrees 0 and 1; bottom panel: normalised small separation between spherical degrees 0 and 2.

In the text

	Fig. E.1. Comparison between the successive adjustments and the observed radial modes frequencies listed in Table E.1. Upper panel: observed frequencies compared to the first order adjustment; lower panel: residual of the first order adjustment compared to the second order fit to those residuals.
In the text

	Fig. E.2. Separated glitches adjustments to the frequencies in Table E.1 of both glithes shown for the spherical degree l = 0. Upper panel: helium glitch; lower panel: convection zone glitch.
In the text

Fig. F.1. Comparison between the observed helium glitch δνo (solid line) and the one resulting from the best fit model δνm (dot-dashed line) for l = 0. We also display the observed glitch as a function of the frequencies (errorbars) as well as the best model glitch associated with the theoretical frequencies (diamond). The observed frequencies have been corrected for surface effects using Sonoi et al. (2015)’s prescription.

In the text

	Fig. F.2. Comparison between the observed frequencies (diamonds) and the best model frequencies (circles) in an échelle diagram. The observed frequencies have been corrected for surface effects using Sonoi et al. (2015)’s prescription.
In the text

	Fig. F.3. Comparison between the observed helium glitch δνo (solid line) and the one resulting from the best fit model δνm (dot-dashed line) for l = 0. We also display the observed glitch as a function of the frequencies (errorbars) as well as the best model glitch associated with the theoretical frequencies (diamond). We include polynomials in n−1 to the basis functions.
In the text