MEASUREMENT OF ACOUSTIC GLITCHES IN SOLAR-TYPE STARS FROM OSCILLATION FREQUENCIES OBSERVED BY KEPLER

In this approach, we considered a broad family of theoretical models which mimic the global properties of the stars and their seismic properties, as listed in Table 1, but cannot be claimed to necessarily have frequencies that very closely match the observed ones. We deliberately spanned a very broad range in each of the global properties. This is because our aim here was to only determine the possible range in the locations of the acoustic glitches in stellar models similar to the target star, and compare the values estimated from the Kepler data. Further, we repeated the exact procedure of one of the methods (B) to obtain the acoustic depths τ_BCZ and τ_HeIIZ from the theoretical frequencies of the stellar models. This allowed us also to investigate possible systematic shifts between the theoretical acoustic locations of the glitches from the models and their estimated values from the oscillatory signal in the model frequencies themselves.

We used the Yale Stellar Evolution Code (YREC; Demarque et al. 2008) to model the stars. A detailed description of the models can be found in Appendix B. The first step of our modeling was to use the average large separation Δ₀ and the frequency of maximum power, ν_max, along with T_eff and metallicity to determine the masses of the stars using a grid-based Yale-Birmingham pipeline (Basu et al. 2010; Gai et al. 2011). Mathur et al. (2012) have shown that grid-based estimates of stellar masses and radii agree very well with those obtained from more detailed modeling of the oscillation frequencies of stars, though with slightly lower precision (see also Silva Aguirre et al. 2011).

For each star, we specified 8–12 initial masses scanning a 2σ range on either side of the mass obtained by the grid modeling. For each initial mass, we modeled the star using three values of the mixing length parameter, α (1.826, 1.7, and 1.5; note that α = 1.826 is the solar calibrated value of α for YREC). For each value of mass and α, we assumed at least four different values of the initial helium abundance, Y₀. In general, Y₀ ranged from 0.25 to 0.30. The models were evolved from zero-age main sequence and the properties were output at short intervals to allow us to calculate oscillation frequencies for the model as it evolved.

All models satisfying the observed constraints on large separation, small separation, T_eff, log g and [Fe/H] were selected for comparison. In order to get a reasonably large number of models, the uncertainty margin was assumed to be 2.0 μHz, 2.0 μHz, 200 K, 0.1 dex, and 0.1 dex, respectively. The number of models selected with these criteria ranged from about 80 to 800. For each of these models, the acoustic depths τ_BCZ and τ_HeIIZ were determined by a procedure exactly similar to method B adopted for the real data. The theoretically computed value of the frequency was taken as the mean value, and the uncertainty was adopted as that of the corresponding mode frequency in the Kepler data of the concerned star. Thus, the model data set mimicked the observed data set in terms of number of modes, specific modes used in the fitting, and the error bars on the frequencies.

The comparison of the fractional acoustic radii of the glitches in the models and our estimates from the Kepler data are shown in Figure 7 for two of the stars in our sample. One of them (KIC010963065) is a main-sequence star, while the other (KIC011244118) is in the subgiant phase. In this figure, we show the fractional acoustic radii of the glitches estimated from the Kepler data by method B along with the values obtained from the frequencies of the models by the same method, as described above. For the latter, the typical uncertainty (not shown in the graphs for the sake of clarity) would be similar to that of the values obtained from the Kepler data, since we have assumed the same uncertainties on the theoretical frequencies as the data. The acoustic depths obtained from the oscillatory signals have been converted to acoustic radii through the relation t = T₀ − τ ≈ (2Δ₀)⁻¹ − τ. The error bars on the values from these methods include the uncertainty propagated in the calculation of T₀ from Δ₀, as listed in Table 2. The figure also shows the theoretical values of the fractional acoustic radii of the glitches from the models calculated using Equations (3) and (5). The two illustrated stars are typical of the sample. Results for other stars are similar to these with the discrepancies between the models and observations being either smaller or higher in a couple of stars.

Zoom In Zoom Out Reset image size

We also compared our estimated acoustic locations of the glitches to those in an optimally fitted model of each star. The fitted model was obtained through the Asteroseismic Modeling Portal (AMP; Metcalfe et al. 2009). The details of the stars studied here can be found in Mathur et al. (2012). Briefly, the models were obtained by optimizing the match between the observed frequencies and the modeled ones. To overcome the issue of the surface effects, we applied the empirical formula of Kjeldsen et al. (2008). The model acoustic radii were calculated as per the definitions given in Equation (3). The comparison of AMP model values of the fractional acoustic radii and those obtained from the four methods is shown in Figure 8.

Zoom In Zoom Out Reset image size

For three of the stars, KIC008006161, KIC006106415, and KIC012009504, independent optimized models were made by fitting frequency ratios as described by Silva Aguirre et al. (2011, 2013). These models were constructed with the GARSTEC code (Weiss & Schlattl 2008), and are also shown in Figure 8.

Overall, we find remarkably good agreement between the values of the acoustic depths of the BCZ determined by different methods. In most cases, the τ_BCZ (or the corresponding T_BCZ) agree with each other well within the quoted 1σ error bars (cf. Table 2 and Figure 6). In terms of relative uncertainties, the τ_BCZ and τ_HeIIZ values from different methods agree pairwise within 10% and 5%, respectively, for most of the stars.

The general agreement between τ_HeIIZ values from different methods is not as good as those for τ_BCZ, although it is mostly within 1σ uncertainties, and the values always match within 5% of each other. While comparing our estimates of the acoustic depths τ_BCZ and τ_HeIIZ from methods A, B, and C to the acoustic radii from method D, or from typical models, we invoked the relationship between the total acoustic radius T₀ and the average large separation Δ₀, i.e., T₀ ≈ (2Δ₀)⁻¹. This relationship is an approximate one. Furthermore the observed average large separation, Δ₀, depends on the adopted frequency range over which it is measured and on the adopted method.

Particular attention should be paid to the systematic biases involved in treating the outermost layers of the star by different seismic diagnostics dealing with the determination of acoustic depths of glitches. Although one may hope that the different methods (seismic diagnostics) will provide similar stellar radii, r, for the acoustic glitches of both the ionization zones of helium and the base of the surface convection zone, the seismically measured acoustic depths, τ(r), of these glitches will in general be different, for they depend on the very details of the adopted seismic diagnostic. For example, method C approximates the outer stellar layers by a polytrope with a polytropic index m = 3.5, allowing for some account of the location of the outer turning point of an incident acoustic wave by approximating the acoustic cutoff frequency, ν_ac ≃ (m + 1)/τ (Houdek & Gough 2007). The location of the upper turning point affects the phase of an incident mode and consequently, also affects the location of the acoustic glitches, typically increasing their acoustic depths when the outer layers are approximated by a polytrope (Houdek & Gough 2007). A (mathematically) convenient way for estimating the phase at the upper turning point is to relate it to a fiducial location in the evanescent region far above the upper turning point, where the mode in the propagating region would "feel" the adiabatic sound speed, c, vanish, i.e., where c = 0. This location, lying well inside the evanescent zone of most of the acoustic modes, defines the acoustic radius or acoustic surface. The squared adiabatic sound speed, c², decreases in the adiabatically stratified region of an outer convection zone nearly linearly with radius (e.g., Balmforth & Gough 1990; Gough 2013). The location of the acoustic radius can therefore be estimated where the linearly outward extrapolated c², with respect to radius r, vanishes. In a solar model, the so-determined location of the seismic surface is about 111 s in acoustic height above the temperature minimum (Lopes & Gough 2001), or above 200 s in acoustic height above the photosphere (Monteiro et al. 1994; Houdek & Gough 2007). This shift may in part explain the differences between the location of the acoustic glitches between the stellar models and our seismically determined values.

There is also a noticeable systematic shift between the acoustic depth of the HeIIZ in a stellar model and its value estimated from the oscillatory signal in the model frequencies (see Figure 7). The fitted value of τ_HeIIZ, as obtained with methods A and B, is typically smaller than the model value by about 100 s, on average. However, we do not find such a systematic shift between the τ_BCZ values of the models and the fits. This is a general feature of most of the stars and might be because the depth of the He ii ionization zone cannot be uniquely defined since the ionization zone covers a finite range of radius. On the other hand, the base of the convection zone is a clearly defined layer within the star. Including the glitch contribution from the first stage of helium ionization in the seismic diagnostic, in addition to the contribution from the second stage of helium ionization, perhaps leads to a more accurate definition of the acoustic depth of helium ionization, as considered in method C, resulting in larger values for the fitted acoustic depth τ_HeIIZ as indicated in Table 2.

In Table 2 and Figures 7 and 8, we did not take account of these differences in the acoustic depths between the methods A–C and the calculated stellar models. It is also worth noticing that the T₀ values of AMP models are systematically smaller (by about 300 s, or, by about 4%, on average) than those estimated from Δ₀ of the Kepler data used in this paper (see Table 3). This may be due to a combination of factors, including imperfect modeling of the stars, especially of the surface layers. The AMP fitting has been subjected to the surface correction technique (Kjeldsen et al. 2008), while this was not included in the analysis of the present data. However, even with stellar models, a difference of up to 2% between the T₀ values derived from the asymptotic relation (Equation (5)) and that from the large separation is expected (Hekker et al. 2013). Nevertheless, this contributes partially to the differences in the fractional acoustic radii of the glitches between AMP models and the estimates from the data. We have not estimated the uncertainties in the model values of the acoustic radii. Thus, it is difficult to make any quantitative comparison between these values and those obtained from the data.

We further discuss the cases of the eight selected stars in detail below (ref. Figures 2–5). The cases of the 11 remaining stars are somewhat similar to these and a general understanding of the issues involved may be obtained from the selected subsample itself.

5.2.1. KIC008006161

5.2.2. KIC006603624

5.2.3. KIC010454113

5.2.4. KIC010963065

5.2.5. KIC004914923

5.2.6. KIC012009504

5.2.7. KIC006933899

5.2.8. KIC011244118

5.2.9. Other Discrepant Cases

In this method, the strategy is to isolate the signature of the acoustic glitches in the frequencies themselves. Compared to the methods involving second differences, described in Appendices A.2 and A.3, this has a different sensitivity to uncertainties (for a discussion on how the errors propagate in comparison with the uncertainties when using frequency combinations, please see the works by Ballot et al. 2004; Houdek & Gough 2007), does not require having frequencies of consecutive order and avoids the additional terms that need to be considered when fitting the expression to frequency differences. However, it is less robust on the convergence to a valid solution (the starting guess must be sufficiently close to the solution), since we do not put any constraints on what we consider to be a smooth component of the frequencies.

Here, we assume that only very low-degree data are available so that any dependence of the oscillatory signal on mode degree can be ignored. The expression for the signature due to the BCZ is taken from Monteiro et al. (1994), Christensen-Dalsgaard et al. (1995), and Monteiro et al. (2000). The expression for the signature due to the HeIIZ is from Monteiro & Thompson (1998) and Monteiro & Thompson (2005), with some minor adaptations and/or simplifications.

In the following, ν_r is a reference frequency, introduced for normalizing the amplitude. One good option would be to use ν_r = ν_max (frequency of maximum power) but the actual value selected is not relevant unless we want to compare the amplitude of the signal of different stars.

The signal from the BCZ, after removing a smooth component from the frequencies, is written as (assuming there is no discontinuity in the temperature gradient at the BCZ)

$$\begin{equation} \delta \nu _\mathrm{BCZ} \simeq A_\mathrm{BCZ} \left({\nu _r \over \nu }\right)^2 \sin (4\pi \tau _\mathrm{BCZ}\nu + 2\phi _\mathrm{BCZ}). \end{equation} \tag{ A1 }$$

The parameters to be determined in this expression are A_BCZ, τ_BCZ, ϕ_BCZ for amplitude, acoustic depth, and phase, respectively. The typical values to expect for a star like the Sun are A_BCZ ∼ 0.1 μHz, τ_BCZ ∼ 2300 s and ϕ_BCZ ∼ π/4.

The signal from HeIIZ, after removing a smooth component from the frequencies, is described by

$$\begin{eqnarray} \delta \nu _\mathrm{\mathrm{HeIIZ}} & \simeq & A_\mathrm{\mathrm{HeIIZ}} \left({\nu _r \over \nu }\right) \sin ^2 (2\pi \beta _\mathrm{\mathrm{HeIIZ}}\, \nu) \nonumber \\ & & \times \cos (4\pi \tau _\mathrm{HeIIZ}\nu + 2\phi _\mathrm{HeIIZ}). \end{eqnarray} \tag{ A2 }$$

The free parameters in this expression are $A_\mathrm{\mathrm{HeIIZ}}$, $\beta _\mathrm{\mathrm{HeIIZ}}$, τ_HeIIZ, and ϕ_HeIIZ corresponding, respectively, to amplitude, acoustic width, acoustic depth, and phase. The typical values to expect for a solar-like star for these parameters are $A_\mathrm{\mathrm{HeIIZ}}\sim 1.0$ μHz, $\beta _\mathrm{\mathrm{HeIIZ}} \sim 130$ s, τ_HeIIZ ∼ 700 s, and ϕ_HeIIZ ∼ π/4. When performing the fitting of Equation (A2), leading to these parameters, it must be noted that the signature described by Equation (A1) is ignored and is treated as noise in the residuals.

If low-degree frequencies are used for a solar-type star, then the signals present in the data correspond to

$$\begin{equation} \nu = \nu _\mathrm{s1} + \delta \nu _\mathrm{BCZ} \end{equation} \tag{ A3 }$$

or

$$\begin{equation} \nu = \nu _\mathrm{s2} + \delta \nu _\mathrm{\mathrm{HeIIZ}}. \end{equation} \tag{ A4 }$$

Here, ν_si (i = 1, 2) represents a "smooth" component of the mode frequency to be removed in the fitting. The parameters to fit these expressions are for Equation (A1) (BCZ): A_BCZ, τ_BCZ, ϕ_BCZ; and for Equation (A2) (HeIIZ): $A_\mathrm{\mathrm{HeIIZ}}$, $\beta _\mathrm{\mathrm{HeIIZ}}$, τ_HeIIZ, ϕ_HeIIZ. In this case, we treat δν_BCZ as noise, since ν_s2 is treated as including ν_s1 and δν_BCZ.

The fitting procedure used is the same method as described by Monteiro et al. (2000). This method uses an iterative process in order to remove the slowly varying trend of the frequencies, leaving the required signature given in either Equation (A1) or Equation (A2). A polynomial in n was fitted to the frequencies of all modes of given degree l separately using a regularized least-squares fit with third-derivative smoothing through a parameter, λ₀ (see Monteiro et al. 1994 for the details). The residuals to those fits were then fitted for all degrees simultaneously. The latter fit is the "signal" either from the BCZ or from the HeIIZ. This procedure was then iterated (by decreasing the smoothing); at each iteration we removed from the frequencies the previously fitted signal and recalculated the smooth component of the frequencies, this time using the smaller value of λ_j (j = 1, 2, 3, ...). The iteration converged when the relative changes to the smooth component of the signal fell below 10⁻⁶ (typically this happened for j ≃ 3).

The initial smoothing, selected through the parameter λ₀, defines the range of wavelengths whose variation we want to isolate. To isolate the signature from the BCZ, a smaller initial value of λ₀ was required in order to ensure that the smooth component also extracts the signature from the HeIIZ. However, when isolating the signature from the HeIIZ, the shorter wavelength signature from the BCZ was also retained. We show the fits to Equations (A1) and (A2) in Figure 2 for six of the Kepler stars that we studied.

In order to estimate the impact on the fitting of the observational uncertainties, we used Monte Carlo simulations. We did so by producing sets of frequencies calculated from the observed values with added random values calculated from a standard normal distribution multiplied by the quoted observational uncertainty. In the fit, we removed frequencies with large uncertainties (>0.5 μHz for the BCZ and >1.0 μHz for the HeIIZ). For each star, we produced 500 sets of frequencies, determining the parameters as the standard deviation of the results for the parameters. Only valid fits were used as long as the number of valid fits represented more than 90% of the simulations. Otherwise, we considered that the signature had not been successfully fit (even if a solution could be found for the observations).

This method involves the second differences of the frequencies to determine τ_BCZ and τ_HeIIZ simultaneously by fitting a functional form to the oscillatory signals.

The oscillatory signal in the frequencies due to an acoustic glitch is quite small and is embedded in the frequencies together with a smooth trend arising from the regular variation of the sound speed in the stellar interior. It can be enhanced by using the second differences

$$\begin{equation} \Delta _2 \nu (n,l) = \nu (n-1,l) - 2\nu (n,l) + \nu (n+1,l), \end{equation} \tag{ A5 }$$

instead of the frequencies ν(n, l) themselves (see, e.g., Gough 1990; Basu et al. 1994, 2004; Mazumdar & Antia 2001).

We fitted the second differences to a suitable function representing the oscillatory signals from the two acoustic glitches (Mazumdar & Antia 2001). We used the following functional form which has been adapted from Houdek & Gough (2007) (Equation (22) therein):

$$\begin{eqnarray} \Delta _2 \nu & = & a_0 \nonumber\\ && + {b_2\over \nu ^2} \; \sin (4\pi \nu \tau _\mathrm{BCZ}+ 2\phi _\mathrm{BCZ}) \nonumber \\ && + c_0\nu \, e^{-c_2\nu ^2} \sin (4\pi \nu \tau _\mathrm{HeIIZ}+ 2\phi _\mathrm{HeIIZ}), \end{eqnarray} \tag{ A6 }$$

where a₀, b₂, c₀, c₂, τ_BCZ, ϕ_BCZ, τ_HeIIZ and ϕ_HeIIZ are eight free parameters of fitting. For one of the stars in our set (KIC010018963), we needed to use a slightly different function in which the constant term representing the smooth trend, a₀, was replaced by a parabolic form: (a₀ + a₁ν + a₂ν²). This was necessitated by the sensitivity of the fitted τ values to small perturbations to the input frequencies. The two τ values for the BCZ and the HeIIZ had uncertainties about 10 times larger (and overlapping) with the constant form, as compared to those with the parabolic form. Although this parabolic form could, in principle, be used for all stars as well, it actually interferes with the slowly varying periodic HeIIZ component in the limited range of observed frequencies and makes it difficult to determine τ_HeIIZ. Therefore, in adopting Equation (A6), we essentially assumed that the smooth trend in the frequencies had been reduced to a constant shift in the process of taking the second differences. We ignored frequencies which have uncertainties of more than 1 μHz.

While Equation (A6) is not the exact form prescribed by Houdek & Gough (2007), it captures the essential elements of that form while keeping the number of free parameters relatively small. The ignored terms can be shown to have relatively smaller contributions. We note that Basu et al. (2004) have shown that the exact form of the amplitudes of the oscillatory signal does not significantly affect the results (however, see also Houdek & Gough 2006). The fits of Equation (A6) to the second differences of some selected stars in our sample are shown in the left panels of Figure 3.

The fitting procedure is similar to the one described by Mazumdar et al. (2012). The fit was carried out through a nonlinear χ² minimization, weighted by the uncertainties in the data. The correlation of uncertainties in the second differences was accounted for by defining the χ² using a covariance matrix. The effects of the uncertainties were considered by repeating the fit for 1000 realizations of the data, produced by perturbing the frequencies by random uncertainties corresponding to a normal distribution with a standard deviation equal to the quoted 1σ uncertainty in the frequencies. The successful convergence of such a nonlinear fitting procedure is somewhat dependent on the choice of reasonable initial guesses. To remove the effect of initial guesses affecting the final fitted parameters, we carried out the fit for multiple combinations of starting values. For each realization, the fitting was repeated for 100 random combinations of initial guesses of the free parameters and the fit which produced the minimum value of χ² was accepted.

The median value of each parameter for 1000 realizations was taken as its fitted value. The ±1σ uncertainty in the parameter was estimated from the range of values covering a 34% area about the median in the histogram of fitted values, assuming the error distribution to be Gaussian. Thus, the quoted uncertainties in these parameters reflect the width of these histograms on two sides of the median value. The histograms for τ_BCZ and τ_HeIIZ and the ranges of initial guesses of the corresponding parameters are shown in the right panels of Figure 3.

In some cases, the fitting of the τ_BCZ parameter suffers from the aliasing problem (Mazumdar & Antia 2001), where a significant fraction of the realizations are fitted with τ_BCZ equal to $\tilde{\tau }_\mathrm{BCZ} \equiv T_0- \tau _\mathrm{BCZ}$. This becomes apparent from the histogram of τ_BCZ which appears bimodal with a reflection around T₀/2 (e.g., for KIC010454113, shown in Figure 3). In such cases, we chose the higher peak in the histogram to represent the true value of τ_BCZ, and the uncertainty in the parameter was calculated after "folding" the histogram about the acoustic mid-point T₀/2, where the acoustic radius was estimated from the mean large separation, Δ₀. For a few stars, there are multiple peaks in the histogram for τ_BCZ, not all of which can be associated with the true depth of the BCZ or its aliased value (e.g., for KIC012009504, shown in the bottom right panel of Figure 3). In such cases, we chose only the most prominent peak to determine the median, and the uncertainty was estimated from the width of that peak. The remaining realizations were also neglected for estimating the other parameters such as τ_HeIIZ.

Approximate expressions for the frequency contributions, δν_i, arising from acoustic glitches in solar-type stars were recently presented by Houdek & Gough (2007, 2011), which we adopt here for producing the results presented in Figure 4 and Table 2. A detailed discussion of the method can be found in those references; therefore, we present only a summary here. The complete expression for δν_i is given by

$$\begin{equation} \delta \nu _i=\delta _\gamma \nu _i+\delta _{\rm c}\nu _i+\delta _{\rm u}\nu _i, \end{equation} \tag{ A7 }$$

where the terms on the right-hand side are the individual acoustic glitch components located at different depths inside the star.

The first component,

$$\begin{eqnarray} \delta _\gamma \nu &=&-\sqrt{2\pi }A_{\rm II}\Delta ^{-1}_{\rm II} \left[\nu +\textstyle \frac{1}{2}(m+1)\Delta _0\right]\nonumber\\ && \times \left [\tilde{\beta }\int _0^{T_0}\kappa ^{-1}{e}^{-{(\tau -\tilde{\eta }\tau _{\rm II})^2 \over 2\tilde{\mu }^2\Delta ^2_{\rm II}}}|x|^{1/2} |{\rm Ai}(-x)|^2\,{d}\tau\right. \nonumber\\ &&\left.+\int _0^{T_0}\kappa ^{-1}{e}^{-{(\tau -\tau _{\rm II})^2\over 2\Delta ^2_{\rm II}}}|x|^{1/2} |{\rm Ai}(-x)|^2\,{d}\tau \right ], \end{eqnarray} \tag{ A8 }$$

arises from the variation in γ₁ induced by helium ionization. The constant m = 3.5 is a representative polytropic index in the expression for the approximate effective phase, ψ, appearing in the argument x = sgn(ψ)|3ψ/2|^2/3 of the Airy function Ai(− x). We approximate ψ as

$$\begin{equation} \psi (\tau)=\kappa \omega \tilde{\tau }-(m+1)\cos ^{-1}[(m+1)/\omega \tilde{\tau }], \end{equation} \tag{ A9 }$$

if $\tilde{\tau }>\tau _{\rm t}$, and

$$\begin{equation} \psi (\tau)=|\kappa |\omega \tilde{\tau }-(m+1)\ln [(m+1)/\omega \tilde{\tau }+|\kappa |], \end{equation} \tag{ A10 }$$

if $\tilde{\tau }{\le }\tau _{\rm t}$, in which $\tilde{\tau }{=}\tau {+}\omega ^{-1}\epsilon _{\rm II}$, with _II (_I = _II) being a phase constant, and τ_t is the location of the upper turning point of the mode. The location of the upper turning point of the oscillation mode is determined approximately from the polytropic representation of the acoustic cutoff frequency, leading to the expressions $\kappa (\tau)=[1-(m+1)^2/\omega ^2\tilde{\tau }^2]^{1/2}$.

The coefficients associated with the glitch contribution from the first stage of helium ionization (He i) (first integral expression in Equation (A8)) are related to the coefficients of the second stage of helium ionization (He ii) by the constant ratios $\tilde{\beta }:=A_{\rm I}\Delta _{\rm II}/A_{\rm II}\Delta _{\rm I}$, $\tilde{\eta }:=\tau _{\rm I}/\tau _{\rm II}$ and $\tilde{\mu }:=\Delta _{\rm I}/\Delta _{\rm II}$, where A_II is the amplitude factor of the oscillatory He ii glitch component, Δ_II is the acoustic width of the glitch, and τ_II ≡ τ_HeIIZ is its acoustic depth beneath the seismic surface. For the constant ratios $\tilde{\beta }$, $\tilde{\eta }$, and $\tilde{\mu }$, we set the values at 0.45, 0.70, and 0.90, respectively, as did Houdek & Gough (2007, 2011). With this approach, the addition of the He i contribution does not introduce additional fitting coefficients in the seismic diagnostic (A8).

The second component in Equation (A7),

$$\begin{eqnarray} \delta _{\rm c}\nu &\simeq & A_{\rm c}\Delta _0^3\nu ^{-2} \left(1+1/16\pi ^2\tau _0^2\nu ^2\right)^{-1/2}\nonumber\\ && \times \left \lbrace \cos [2\psi _{\rm c}+\tan ^{-1}(4\pi \tau _0\nu)]\right. \nonumber \\ && -\left(16\pi ^2\tilde{\tau }_\mathrm{BCZ}^2\nu ^2\!+\!1\right)^{1/2} \rbrace, \, \end{eqnarray} \tag{ A11 }$$

arises from the acoustic glitch at the base of the convection zone resulting from a near discontinuity (a true discontinuity in theoretical models using local mixing length theory with a non-zero mixing length at the lower boundary of the convection zone) in the second derivative of density. We model this acoustic glitch with a discontinuity in the squared acoustic cutoff frequency $\omega _{\rm c}^2$ at τ_BCZ, with A_c being proportional to the jump in $\omega _{\rm c}^2$, coupled with an exponential relaxation to a putative, glitch-free, model in the radiative zone beneath, with a relaxation time scale τ₀ = 80 s, as did Houdek & Gough (2007, 2011). This leads to

$$\begin{eqnarray} \psi _{\rm c}&=&\kappa _{\rm c}\omega \tilde{\tau }_\mathrm{BCZ} \nonumber\\ && - (m+1)\cos ^{-1}[(m+1)/\tilde{\tau }_\mathrm{BCZ}\omega ] \nonumber\\ && +\pi /4, \end{eqnarray} \tag{ A12 }$$

where κ_c = κ(τ_BCZ) and $\tilde{\tau }_\mathrm{BCZ}=\tau _\mathrm{BCZ}+\omega ^{-1}\epsilon _{\rm c}$ with _c being a constant phase.

The additional upper-glitch component δ_uν_i (i enumerates individual frequencies), which is produced, in part, by wave refraction in the stellar core, by the ionization of hydrogen and by the upper superadiabatic boundary layer of the envelope convection zone, is difficult to model. We approximate its contribution to Δ₂ν as a series of inverse powers of ν, truncated at the cubic order:

$$\begin{equation} \Delta _2\delta _{\rm u}\nu _i=\sum _{k=0}^3a_k\nu _i^{-k}. \end{equation} \tag{ A13 }$$

The 11 coefficients η_α = (A_II, Δ_II, τ_HeIIZ, _II, A_c, τ_BCZ, _c, a₀, a₁, a₂, a₃), α = 1, ..., 11, were found by fitting the second differences (cf. Equation (A5)),

$$\begin{eqnarray} \Delta _{2i}\nu (n,l)\; &{:=}&\nu (n-1,l)-2\nu (n,l)+\nu (n+1,l)\nonumber \\ &\simeq & \;\Delta _{2i}(\delta _\gamma \nu +\delta _{\rm c}\nu +\delta _{\rm u}\nu) =: g_i(\nu _j;\eta _\alpha) \end{eqnarray} \tag{ A14 }$$

to the corresponding observations by minimizing

$$\begin{equation} E_{\rm g}=\sum _{i,j}(\Delta _{2i}\nu -g_i)C^{-1}_{\Delta ij}(\Delta _{2j}\nu -g_j), \end{equation} \tag{ A15 }$$

where $C^{-1}_{\Delta ij}$ is the (i, j) element of the inverse of the covariance matrix C$_\Delta$ of the observational uncertainties in Δ_2iν, computed, perforce, under the assumption that the uncertainties in the frequency data ν_i are independent. The covariance matrix C_ηαγ of the uncertainties in the fitting coefficients η_α was established by Monte Carlo simulation.

Results for six of the stars in our sample are shown in Figure 4 in which the individual acoustic glitch contributions are illustrated in the lower panels.

Determinations of acoustic depths of glitches are biased by surface effects (e.g., Christensen-Dalsgaard et al. 1995). A way to remove such biases is to consider acoustic radii of the signatures. First, it is possible to perform posterior determinations of radii by comparing depths to the total acoustic radius of the star derived from the average large separation Δ₀ (see Ballot et al. 2004); it is also possible to directly measure acoustic radii by considering the small separations d₀₁ and d₁₀ or the frequency ratios r₀₁ and r₁₀ as shown by Roxburgh & Vorontsov (2003).

We use the three-point differences:

$$\begin{equation} d_{01,n}=\frac{1}{2}(2\nu _{n,0}-\nu _{n-1,1}-\nu _{n,1}), \end{equation} \tag{ A16 }$$

$$\begin{equation} d_{10,n}=-\frac{1}{2}(2\nu _{n,1}-\nu _{n,0}-\nu _{n+1,0}) \end{equation} \tag{ A17 }$$

and the corresponding ratios:

$$\begin{eqnarray} r_{01,n}=\frac{d_{01,n}}{\Delta \nu _{1,n}},& r_{10,n} = \frac{d_{10,n}}{\Delta \nu _{0,n+1}}. \end{eqnarray} \tag{ A18 }$$

We denote the sets {d₀₁, d₁₀} and {r₀₁, r₁₀} by d₀₁₀ and r₀₁₀, respectively.

Using these variables, the main contributions of outer layers are removed (see Roxburgh 2005). The global trend of these variables then gives information on the core of the star (e.g., Silva Aguirre et al. 2011; Cunha & Brandão 2011). Nevertheless, the most internal glitches, such as the BCZ, also imprint their signatures over the global trend. Using solar data, Roxburgh (2009) showed that we can recover the acoustic radius of the BCZ (T_BCZ) by the use of a Fourier transform on the residuals obtained after removing the global trend. As a consequence of this approach, information about surface layers, including He i and He ii ionization zones, are lost.

We used an approach similar to Roxburgh (2009) and develop a semi-automatic pipeline which extracts glitches from a frequency table of l = 0 and 1 modes. Instead of fitting a background first, then making a Fourier transform, we do both simultaneously by fitting the variable y = ν*r₀₁₀ (or y = ν*d₀₁₀) with the following expression:

$$\begin{equation} f(\nu)= \sum _{k=0}^m \frac{c_k}{(\nu +\nu _\mathrm{r})^k} + A\sin (4\pi T \nu + \phi), \end{equation} \tag{ A19 }$$

where ν* = ν/ν_r and ν_r is a reference frequency. We have m + 3 free parameters: {c_k}, A, T, and ϕ.

To be able to estimate reliable uncertainties for T_BCZ, we performed a Markov Chain Monte Carlo (MCMC) to fit the data. Our MCMC fitting algorithm is close to the one described, for example, in Benomar et al. (2009) or Handberg & Campante (2011), but without parallel tempering. Moreover, in our case, the noise is not multiplicative and exponential but additive and normal. We fixed m = 2 and ν_r = 0.8ν_max, similar to the value adopted for the Sun by Roxburgh (2009). We used uniform priors for A, T, and ϕ. The prior for A was very broad (between 0 and 10max (|y|)); ϕ was 2π-periodic with a uniform prior over [0, 2π], and T was restricted to [δT, T₀ − δT], where δT = (max (ν) − min (ν))⁻¹. To determine priors for {c_k}, we first performed a standard linear least-square fitting without including the oscillatory component. Once we obtained the fitted values, $\tilde{c}_k$, and associated uncertainties, σ_k, we used for c_k priors which are uniform over $[\tilde{c}_k- 3\sigma _k,\tilde{c}_k+ 3\sigma _k]$ with Gaussian decay around this interval. After a 200,000-iteration learning phase (the learning method we used is based on Benomar et al. 2009), MCMC was run for 10⁷ iterations. We tested the convergence by verifying the consistence of results obtained by using the first and second half of the chain.

Posterior PDFs for T are plotted for different stars in Figure 5. Some PDFs exhibit a unique and clear mode that we attributed to the BCZ. T_BCZ was then estimated as the median of the distribution and 1σ error bars were taken as a 68% confidence limit of the PDF. Nevertheless, some other PDFs exhibit other small peaks. They can be due to noise but also can be due to other glitch signatures present in the data, especially tiny residuals from the outer layers. To take care of such situations, we also computed T_BCZ by isolating the highest peak and fitting it with a Gaussian profile. When the secondary peaks are small enough, the values and associated uncertainties obtained through the two methods are very similar. We have been successful mainly on main-sequence stars, while subgiant stars presenting mixed modes strongly disturb the analysis.