{'fraction_non_alphanumeric': 0.030760603449454953, 'fraction_numerical': 0.016262497777413044, 'mean_word_length': 4.361443132653809, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comprehensive stellar seismic analysis : New method exploiting the glitches information in solar-like pulsators

Aims: We develop a method that provides a comprehensive analysis of the oscillation spectra of solar-like pulsators. We define new seismic indicators that should be as uncorrelated and as precise as possible and should hold detailed information about stellar interiors. This is essential to improve the quality of the results obtained from asteroseismology as it will provide better stellar models which in turn can be used to refine inferences made in exoplanetology and galactic archaeology. Methods: The presented method - WhoSGlAd - relies on Gram-Schmidt's orthogonalisation process. A Euclidean vector subspace of functions is defined and the oscillation frequencies are projected over an orthonormal basis in a specific order. This allows the obtention of independent coefficients that we combine to define independent seismic indicators. Results: The developed method has been shown to be stable and to converge efficiently for solar-like pulsators. Thus, detailed and precise inferences can be obtained on the mass, the age, the chemical composition and the undershooting in the interior of the studied stars. However, attention has to be paid when studying the helium glitch as there seems to be a degeneracy between the influence of the helium abundance and that of the heavy elements on the glitch amplitude. As an example, we analyse the 16CygA (HD 186408) oscillation spectrum to provide an illustration of the capabilities of the method.


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: WhoSGlAdfor Whole Spectrum and Glitches Adjustment -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.

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, i.e. 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.)

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 -WhoS-GlAd -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: with ν th and ν obs , the theoretical and observed 1 frequencies. Taking advantage of the scalar product defined above and the associated norm, this simply becomes: 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 subspace 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 subspace. 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 subspace. 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 App. D.
If we write j and j 0 the indices associated with the basis elements, p j the former basis elements, q j0 the orthonormal basis elements, and R −1 j,j0 the transformation matrix, we have: 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 to note that the projections will 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 will then be given by: 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).

Smooth component
Now that the mathematical context is given, we may detail the vector subspace 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: 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 R −1 l,k,k0 . And the orthonormal basis elements are: which yield: where thanks to the introduction of δ ll ′ in Eqs. 6 and 7, the sum over l ′ collapses over a 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: where ∆ = 2 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 : 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 χ 2 value. We get: 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 2 . 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 subspaces. 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: Then, we have to define the scalar product of x and y as: 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 on 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 Sect. 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 χ 2 fitting. Finally, as we have three polynomials for each value of l, we built, for the smooth component, a vector subspace 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 3 observed frequencies of each degree are necessary to apply the developed method.

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: where the first term takes the helium glitch into account and the second, the convection zone glitch. The quantities τ He , φ He , and A 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 A CZ in the case of the base of the envelope convective zone. Finally, c 2 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 n = n + l 2 for the helium glitch: with k = (4, 5); and, for the convection zone glitch: To obtain the above expressions we replaced the value of the frequency ν by its first order approximation from the asymptotic formulation (9): n∆. 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 n∆. 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: 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 subspace associated with the glitch is only of dimension 6 and this subspace 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: [C He,k q HeCk ( n, T He ) + S He,k q HeSk ( n, T He )] We draw the reader's attention 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, i.e.: 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 Γ 1 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 χ 2 value. Also, we observed that by using an initial value of τ He different from that at the Γ 1 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 Γ 1 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 3 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 nonlinear routine that searches for a stellar model accounting at best for the seismic and non-seismic observables.

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.

Models
The grid of models we used was computed using CLES (Scuflaire et al. 2008b) combined with LOSC (Scuflaire et al. 2008a) 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 ), 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 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.90M ⊙ to 1.30M ⊙ by steps of 0.01M ⊙ and are in the main sequence phase. Moreover, each model has an initial composition of Y 0 = 0.25 and Z 0 = 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).

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), equation (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: 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, (R −1 l,1,1 ) 2 is the variance of ∆ l . The weighted mean of the large separations thus yields: Finally, we expect from Ulrich (1986) that ∆ should be an estimator of the mean stellar density.

Normalised small separation
Two commonly used indicators are the small separations d 01 (n) and d 02 (n) of which the definitions are: 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 (ν 0 − ν l ) /∆ 0 . Assuming expression (9) to be exact, this difference is n 0 + ǫ 0 − (n l + ǫ l + l/2). Then, we added −n 0 + n l + l/2 to the expression to make its value come close to ǫ 0 − ǫ l . We then obtained the following expression: where ν l and n l 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 ν l equals a l,0 R −1 l,0,0 as it is the fitting of the frequencies of degree l to a constant value. Finally, it has to be stressed 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. Fig. 1 shows the evolution of those indicators along the grid presented in Sect. 3.1 4 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 ofr 02 , 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 ∆ 0 of σ (∆ 0 ) = 5 10 −3 µHz.
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 of Appendix B. 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 thatr 02 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, r 01 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 (HD128620), it should provide an upper limit on the amount of convective-core overshooting. We also draw attention to the fact thatr 01 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,r 03 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, fromr 01 , we should expect masses between 1.06M ⊙ and 1.11M ⊙ . Then, fromr 02 , the expected values are in between 1.06M ⊙ and 1.07M ⊙ . Again,r 03 does not add some information as the values range from 1.05M ⊙ to 1.06M ⊙ . As a consequence, we would expect the mass of 16 Cygni A to be around 1.06M ⊙ . 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 = 0.25, Z 0 = 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.

∆ 0l indicators
As it has been shown in several studies, the combination of the small separation ratios r 01 and r 10 first introduced by Roxburgh & Vorontsov (2003) into r 010 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 indicatorr 01 represents an estimator of the mean value of r 010 . We may therefore build an indicator for its slope as follows: 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, ∆ 01 indeed represents the slope of the frequency ratio r 01 .
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 ∆ 0 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.
Again, it is possible to get an estimate of the mass value for the given composition and physics. From ∆ 01 , we expect masses ranging from 1.06M ⊙ to 1.07M ⊙ .Then, from ∆ 02 , we expect that they lie between 1.05M ⊙ and 1.07M ⊙ . Finally, from ∆ 03 , the mass should be between 1.00M ⊙ and 1.05M ⊙ . 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.06M ⊙ 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 0 , and Z/X 0 and also to solve the slight mass discrepancy observed.
Finally, as detailed above, the simultaneous use of botĥ r 01 and ∆ 01 allows to provide estimations of the extent of stellar mixed cores. To illustrate this, we plotted main sequence evolutionary tracks for a 1.2M ⊙ , Y 0 = 0.25 and X 0 = 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 = α ov min (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 (∆ 01 ,r 01 ) 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.

3.2.4.ǫ
To provide an estimate of ǫ, we define the following vector subspace, where frequencies are described as: where∆ and K are free parameters. Then, we define an orthonormal basis over this subspace: q 0 and q 1 . 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.02M ⊙ . 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.

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: A He = C 2 He,5 + S 2 He,5 + C 2 He,4 + S 2 He,4 .
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. Fig. 5 has been computed for stars with a fixed (Z/X) 0 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 = 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. 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 Γ 1 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.052M ⊙ 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  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. be of great help. This clearly shows that the A He -Y f relation is model dependent which should never be forgotten. Fig. 7 illustrates this degeneracy. It represents the profile of Γ 1 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.

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:  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 = α under min (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.
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.052M ⊙ , an initial hydrogen mass fraction of X 0 = 0.744, and an initial metal mass fraction of Z 0 = 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.

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 χ 2 value from the results without including the convective zone glitch (about 10% variation). These results are similar to those of Verma et al. (2014).

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.25M ⊙ 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. Fig. 10 shows this limitation. Indeed, for 1.25M ⊙ 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. Also, we observe for 1.3M ⊙ 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 Γ 1 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 Γ 1 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. Fig. 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 Γ 1 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 Γ 1 profile.
Finally, we observe that, in a lnP − lnT diagram, the curves for the different models in the temperature region of the Γ 1 dip are parallel to each other but at various height. A higher lnP − lnT curve corresponds to a shallower Γ 1 dip. We also observe this behaviour with a fixed com- position and a variable mass. Higher masses models have lower curves and deeper depressions in Γ 1 . Therefore, the toy model allows to understand the influence of the mass on the helium glitch amplitude as well.

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 makes 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.

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 ∆,r 01 ,r 02 , 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 0 and derive the best fit values of ∆,r 01 , and r 02 . This results in a set of values for the mass, the age, and (Z/X) 0 . Then, the adjustment is done by applying the secant method to find the value of X 0 in best agreement with the target A He . At each step of the secant algorithmat each value of X 0 -, Levenberg-Marqardt's algorithm (L-M) computes the optimal set of stellar parameters giving the best fit of ∆,r 01 , andr 02 . 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 0 , and (Z/X) 0 fitting the observed parameters. We used the   secant method in order to diminish the computational time needed to converge towards the solution.

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 −1 . 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 −1 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 −1 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), 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 3 parameters adjustment of ∆,r 01 , andr 02 as a function of X 0 . We also show the evolution of the value of (Z/X) 0 . 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σ 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. 4.6 10 −3 3.8 10 −3 0.1 10 −3 ∆ 02 5.9 10 −3 4.8 10 −3 0.1 10 −3 ∆ 03 14.9 10 −3 10.6 10 −3 0.6 10 −3 Finally, Fig. 13 shows the fitted helium glitches 5 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 ∆ 0 ,r 01 ,r 02 , and A He we obtain a good representation of the observed frequencies. We also show in App. 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ǫ. Fig. 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ǫ.
6. Discussion and conclusions

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 δνo fit δνo with σ δνm fit δνm 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 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.

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. The Appendix B shows the evolution of the classical indicators r 01 (n) and r 02 (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.

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,r 0l . 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 one specific spherical degree l and 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 010 which, combined withr 01 can be used to constrain the overshooting parameter as shown in Fig. 3.
Moreover, we observe in Fig. 1 that the indicatorr 02 , 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 indicatorr 01 presents a degeneracy, as a turn-off occurs in the ∆ 0 -r 01 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 Γ 1 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.

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.3M ⊙ 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 of the 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 nonseismic data as we have shown that, for example, there exists a degeneracy between the helium glitch amplitude, the helium surface abundance and the metallicity 6 (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).

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.  λ 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 2 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 3 . 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 χ 2 ∼ 1 without regularisation to χ 2 ∼ 500 with λ = 10 3 ).

Appendix B: Comparison with usual indicators
We present in Fig. B.1 the evolution of the 'usual' indicators r 01 (n) and r 02 (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) 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.06M⊙. The top panel is the three points normalised small separation between spherical degrees 0 and 1, the middle one is the five points normalised small separation between spherical degrees 0 and 1 and the bottom panel is the normalised small separation between spherical degrees 0 and 2.
Appendix C: Γ 1 toy model To build a toy model for Γ 1 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 parti-cles per unit volume He ++ such that: where m u is the atomic mass unit; -The helium ionisation fraction: x He = He ++ He + + He ++ , (C.2) which equals 0 if the helium is ionised only once and 1 if it is fully ionised; -The electron number of particles per unit volume: -The total number of particles: where Z i is the mass fraction of the metal labeled i and A i its mass number.
Using Saha's equation (Saha 1920): 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 : 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: Then we have the following expressions: where µ is the mean molecular weight and: Finally, we may insert their values in the relation linking Γ 1 , P T , P ρ , and c v : 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 0 : p j0 . Let us also assume that we have already built the set of orthonormal basis functions up to index j 0 − 1, that is the set q 1 , · · · , 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: where u j0 is the basis element orthogonal to the set q 1 , · · · , q j0−1 . Finally, it is normalised to obtain: It is also possible to express this process as a QR decomposition. To do so we call P l the matrix of initial polynomials and Q l the matrix of orthonormal polynomials for a given spherical degree: q 0 (n min ) · · · q 2 (n min ) . . . . . . q 0 (n max ) · · · q 2 (n max )    = (q 0 · · · q 2 ) l , 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: where, using our definition of the scalar product: 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 n = n + l/2): Article number, page 17 of 20 p He,C,5 (n + 3/2) · · · p C,S (n + 3/2) q He,C,5 (n, l) · · · q C,S (n, l) . . . . . . .
where − −−−− → q k |p g l=0 is a column vector whose rows are the successive scalar products of the basis elements with the glitch function -denoted by the g index -such that − −−−− → q k |p g l=0 = q 0 |p g l=0 · · · q 2 |p g l=0 T . We may again write the QR decomposition as:

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 −1 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 0 ,n 1 ,n 2 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 R −1 j,j0 .  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  several cases. First we show the adjustment to the frequencies corrected for surface effects using a Lorentzian profile as in Sonoi et al. (2015). Fig. 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. Fig. 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 −1 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.