Solar-cycle–related Changes in the Helium Ionization Zones of the Sun

It has been known for about 35 years that solar oscillation frequencies increase with activity (Woodard & Noyes 1985; Elsworth et al. 1990; Libbrecht & Woodard 1990). Subsequent analyses showed that the frequency changes were correlated with the average line-of-sight magnetic field over the solar surface (Woodard et al. 1991), and this was confirmed later using much longer data sets (e.g., Howe et al. 2002). The correlation between frequency shifts and activity is also seen for other solar activity proxies (Jain et al. 2009; Howe et al. 2018, etc.). Both ground-based and space-based data have shown that solar dynamics changes dramatically over the solar cycle (see e.g., Howe 2009; Basu & Antia 2019; etc.). Changes in solar structure, however, are more difficult to determine. There have been some indications of changes in structure at the base of the convection zone (CZ; Eff-Darwich et al. 2002; Baldner & Basu 2008) and indications are that changes are very small. There is also evidence that the immediate sub-surface layers of the Sun change with activity, and that there is a latitudinal dependence of the changes (Antia et al. 2001; Basu et al. 2007).

Basu & Mandel (2004, hereafter BM04) examined the then-available data from the Global Oscillation Network Group (GONG; Hill et al. 1996) and the Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory spacecraft (Scherrer et al. 1995) to look for changes slightly deeper into the CZ, specifically the second helium ionization zone around 0.98R_⊙ based on work by Gough (2002), who suggested that there could be solar-cycle related variations in that region. The helium ionization zone is an "acoustic glitch," a region where the sound speed changes rapidly over a region shorter than the radial wavelength of the modes. A second, sharper glitch is the base of the CZ where the glitch is caused because of the transition of the temperature gradient from being radiative in the deeper layers to adiabatic in the CZ. In the ionization zone, the glitch is caused by the decrease in the adiabatic index, Γ₁, because of ionization. An acoustic glitch contributes a characteristic oscillatory signature in the frequencies (Gough 1990). The wavelength of the oscillations is the acoustic depth of the glitch while the amplitude of the oscillations is governed by how "severe" the glitch is. In the case of the He II ionization zone; normally, the amplitude would be a measure of the helium abundance—the higher the abundance of helium, the higher the amplitude. A change in Γ₁ can also occur in the presence of magnetic fields (Basu et al. 2004) because of changes in the local temperature. Thus a changing magnetic field in the layers around the He II ionization zone will change the amplitude of the signal. BM04 measured the amplitude and found that it decreases as the flux of the 10.7 cm flux, an activity indicator, increases.

The BM04 results were confirmed by Verner et al. (2006) and Ballot et al. (2006) using low-degree data from the Birmingham Solar Oscillation Network (BiSON) and the Global Oscillations at Low Frequencies (GOLF), respectively. BM04 interpreted the change in the ionization zone as one caused by a change in temperature. However, these results, and those of Verner et al. (2006) and Ballot et al. (2006), were challenged by Gough (2013). Gough (2013) analyzed the potential significance of the results under the assumption that the change in the amplitude of the He ii signal is caused by a dilution of the influence of the ionization zone on the wave propagation speed by a broadly distributed magnetic field. He concluded that if the variations of the He ii signal were a direct result of the presence of a temporally varying large-scale magnetic field, then the total solar cycle change of the spatial average of the magnetic field 〈B²〉 in the vicinity of the He ii ionization zone is greater than most estimates.

The BM04 analysis work was done with data for about half a solar cycle. We now have data for two solar cycles. In this work, we repeat the work of BM04 but use improved fitting methods. BM04 had not taken correlated errors into account; we do that now. Additionally, the previous work had assumed that the uncertainty of the result of every set is the same; we abandon that assumption and determine uncertainties for each set by boot-strapping. Another change that we make is that we use two different methods to fit the data—the model is highly non-linear and there is often a tendency for the solution to be stuck at a local χ² minimum.

We use helioseismic data from three sources: (I) the ground-based GONG (Hill et al. 1996), (II) the MDI on board the Solar and Heliospheric Observatory spacecraft (Scherrer et al. 1995), and (III) the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012) on board the Solar Dynamics Observatory.

We use 236 GONG data sets that cover a period from 1995 May 7 to 2018 March 18. The data are labeled by GONG "months," each "month" being 36 days long. Solar oscillation frequencies and splittings of sets starting Month 2 are obtained using 108 day (i.e., three GONG months) time series. There is an overlap of 72 days between different data sets. Data from MDI cover the period from 1996 May 1 to 2011 April 24, with frequencies in each set obtained from 72 day time series; the sets have no overlap in time. We use 48 sets of HMI data with the first set starting at 2010 April 30 and the last ending on 2019 October 16. Like MDI, HMI frequencies sets are obtained with non-overlapping 72 day time series. Note that we essentially use data covering two full solar cycles. We use the 10.7 cm radio flux as an indicator of solar activity (Tapping & Morton 2013 and Tapping 2013) ³

In order to compare our results with those of BM04, we use the same model to determine the amplitude of the glitch signature. As with that work, we first amplify the oscillatory signature by taking the fourth differences of the frequencies:

$$\begin{eqnarray}\begin{array}{rcl}{\delta }^{4}{\nu }_{n,{\ell }} & = & {\nu }_{n+2,{\ell }}-4{\nu }_{n+1,{\ell }}+6{\nu }_{n,{\ell }}\\ & & -4{\nu }_{n-1,{\ell }}+{\nu }_{n-2,{\ell }},\end{array}\end{eqnarray} \tag{ 1 }$$

where ν_n,ℓ is the frequency of a mode of radial order n and degree ℓ. The fourth differences are fitted to the form

$$\begin{eqnarray}\begin{array}{rcl}{\delta }^{4}\nu & = & \left[{a}_{1}+{a}_{2}\nu +({a}_{3}+{a}_{4}L)/{\nu }^{2}\right]\\ & & +\left({a}_{5}+\displaystyle \frac{{a}_{6}}{{\nu }_{m}^{2}}+\displaystyle \frac{{a}_{7}L}{{\nu }_{m}^{2}}+\displaystyle \frac{{a}_{8}L}{{\nu }_{m}^{4}}\right)\sin (2{\nu }_{m}{\tau }_{a}+{\phi }_{a})\\ & & +\left({b}_{1}+\displaystyle \frac{{b}_{2}}{{\nu }_{n}^{2}}+\displaystyle \frac{{b}_{3}L}{{\nu }_{n}^{2}}+\displaystyle \frac{{b}_{4}L}{{\nu }_{n}^{4}}\right)\sin (2{\nu }_{n}{\tau }_{b}+{\phi }_{b})\end{array}\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}L & = & {\ell }({\ell }+1),\quad {\nu }_{m}=\nu -{\gamma }_{a}L/2{\tau }_{a}\nu ,\,\mathrm{and}\\ {\nu }_{n} & = & \nu -{\gamma }_{b}L/2{\tau }_{b}\nu .\end{array}\end{eqnarray} \tag{ 3 }$$

The parameters a₁–a₄ define the overall smooth term, a₅–a₈, τ_a , and ϕ_a define the oscillatory component due to the He ii ionization zone, and the terms b₁–b₄, τ_b , and ϕ_b define the oscillatory component due to the base of the CZ. The frequencies of the two components, τ_a and τ_b , are approximately the acoustic depths of the He ii ionization zone and the CZ base, respectively; they also include a contribution from the frequency-dependent part of the phases ϕ_a and ϕ_b that is not taken into account explicitly. The terms γ_a and γ_b take into account the fact that the response to acoustic glitches has a degree dependence, and these two parameters are fitted along with the other parameters. The shifted frequencies ν_m and ν_n in Equation (3) are needed to account for the degree dependence of the response to acoustic glitches; the derivation of this term can be found in Basu et al. (1994).

The parameters are obtained with non-linear minimum-χ² fits, and this is where we differ from BM04. Equation (1) shows that the errors of neighboring fourth differences will be correlated. BM04 did not take error correlations into account and assumed that the diagonal terms of the full error-covariance matrix sufficed. In this work, we use the full error-covariance matrix, and hence, we define χ² as

$$\begin{eqnarray}&&{\chi }^{2}={\left(y-y^{\prime} \right)}^{T}{C}^{-1}(y-y^{\prime} ),\end{eqnarray} \tag{ 4 }$$

where y is a column matrix with the data, $y^{\prime}$ the result of the model fit, and C is the covariance matrix.

Like BM04, we fit the form in Equation (2) in two degree ranges 5 ≤ ℓ ≤ 25, and 30 ≤ ℓ ≤ 60. All modes in the first set have lower turning points below the CZ base. For the second set, we choose only those modes that have lower turning points above the CZ base. Modes in the first set (the "low-degree" set) sample both glitches, the He II ionization zone, as well as the CZ base and, consequently, all three terms of Equation (1) are fitted. Modes in the second set (the "high-degree" set) sample only the He II ionization zone and, consequently, only the first two terms of Equation (2) apply. Given that the amplitude of the two oscillatory terms in Equation (2) are a function of frequency, like BM04, we use an average amplitude between 2 and 3.5 mHz for our investigation. We define the average amplitude 〈A〉 such that it does not depend on the degree. Defining ${\nu }_{\min }$ and ${\nu }_{\max }$ to be the minimum and maximum frequency, we define 〈A〉 as

$$\begin{eqnarray}\begin{array}{rcl}\langle A\rangle & = & \displaystyle \frac{{\displaystyle \int }_{{\nu }_{\min }}^{{\nu }_{\min }}\left({a}_{5}+\tfrac{{a}_{6}}{{\nu }^{2}}\right)d\nu }{{\displaystyle \int }_{{\nu }_{\min }}^{{\nu }_{\min }}d\nu }\\ & = & \displaystyle \frac{1}{{\nu }_{\max }-{\nu }_{\min }}\left(\Space{0ex}{3.5ex}{0ex}{a}_{5}({\nu }_{\max }-{\nu }_{\min })\right.\\ & & +\left.{a}_{6}\left[\displaystyle \frac{1}{{\nu }_{\max }}-\displaystyle \frac{1}{{\nu }_{\min }}\right]\right).\end{array}\end{eqnarray} \tag{ 5 }$$

We use two different codes to do the fits. The first is the same as the one used by BM04 which uses the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm as implemented by Antia (2002). The second is one that uses the python code scipy.curve, which uses the Levenberg-Marquardt (LM) algorithm to fit the data. The highly non-linear nature of the model means that the χ² surface has many local minima; given the differences between the two algorithms, a similar result returned by the two different methods is unlikely to be that from a local minimum. The complicated model also means that whether or not we obtain a successful fit within the stipulated number of iteration depends on the initial guesses. In another departure from BM04, we have increased the number of initial guesses in order to obtain 250 sets of fitted parameters to find a global minimum. However, we go even further and determine the uncertainties using a traditional boot-strapping method by simulating 100 realizations of the observations, fitting each one of them in the same manner as the original data. We use the aggregate of results to determine the value of the mean amplitudes and their uncertainties.

We first compare the amplitudes obtained by BM04 with the new fits with the same code but using the full error-correlation matrix; these can be seen in Figure 1. It should be noted that BM04 only fitted non-overlapping sets of GONG data and there were no HMI data at that time. Note that while there are small differences in the result obtained—the differences are at most 1.5σ, there are systematic differences. In the case of GONG, there is a slope with respect to the 10.7 cm flux which implies that any correlation that we obtain will be somewhat different than that found by BM04. For MDI, the differences appear to be a constant. For both GONG and MDI, the new amplitudes are higher than the old ones for the higher-degree cases.

Zoom In Zoom Out Reset image size

In Figure 2 we show the difference in results obtained with the BFGS code and the LM code. We were unable to get stable fits with the LM code for the high-degree set sets. The number of converged results were too low to estimate the amplitude and their uncertainties. We are not completely sure why this is the case; the solutions for most parts seemed to get stuck in local minima and often the reduced χ² was greater than 100. Consequently, we concentrate on results obtained for the low-degree set. The results for the GONG sets obtained by the two codes are not very different, neither are those for the HMI sets. The difference in results for the MDI sets, on the other hand, show a somewhat significant slope with respect to the 10.7 cm flux. However, the quality of the LM fits to the space-based data are worse than the GONG case, most likely because each GONG set is obtained with 105 days of data, while MDI and HMI sets were obtained with only 72 days of data. The difference in the results obtained by the two codes thus has to be taken into account while judging if there is a change in the amplitude of the helium glitch signature with change in solar activity.

Zoom In Zoom Out Reset image size

The average amplitude of the glitch signature is plotted as a function of the 10.7 cm flux in Figure 3. The results shown are for the BFGS code. In all panels, we see that that the amplitude decreases when the activity index increases. The slope is much larger than those caused by the changes caused by using the full error-correlation matrix. The slope of the change as given by the LM fits to the low-degree sets corresponding to panels (a), (c), and (e) are, respectively,  − 2.71 × 10⁻⁴ ± 4.67 × 10⁻⁵,  − 2.88 × 10⁻⁴ ± 8.42 × 10⁻⁵ and  − 2.23 × 10⁻⁴ ± 9.61 × 10⁻⁵. Thus, the slopes are less steep than those returned by the BFGS code; nevertheless, they are significant. We therefore believe that, although the exact correlation between the glitch amplitude the 10.7 cm flux is uncertain, the trend is robust.

Zoom In Zoom Out Reset image size

The two solar cycles' worth of data allows us to plot the results against time to examine whether or not a change with solar activity is easily discernible. We show this in Figure 4. Note that the time variation and the anti-correlation with solar activity is very clear for the low-degree sets, though visually less clear for the high-degree ones.

Zoom In Zoom Out Reset image size

We have repeated the work of BM04 with data for two complete solar cycles to determine whether or not there are changes inside the Sun at depths corresponding to the He ii ionization zone. This was done by measuring the amplitude of the signature of the acoustic glitch caused by the second ionization of helium. We find that the amplitude of the glitch indeed changes with time and that it is anti-correlated with solar activity. This behavior is seen in both ground-based and space-based data. This shows, unambiguously, that there are solar-cycle related changes in the Sun at these depths.

The results, which are very similar to those obtained by BM04, lead to the question of why the Gough (2013) analysis led to magnetic-field estimates of nearly 300 Kilo Gauss, an unrealistically large field. We believe there are two reasons for this. The first is that the analysis was based on the common assumption that the oscillatory signature in the frequencies arises solely from the dip in Γ₁ caused by He ii ionization. Work that has been done subsequently (Broomhall et al. 2014; Verma et al. 2014) however, showed that the signature arises from a depth that corresponds to the peak in Γ₁ between the Hei and He ii ionization zones (shown in Figure 5 as the vertical dashed line). Shallower layers have lower gas pressure and hence need smaller magnetic fields to change their structure. Direct inversion results of solar-cycle–related changes in global modes are too noisy to determine how changes in solar magnetic activity modify the Γ₁ profile; consequently, we look instead to Γ₁ changes below sunspots of different strengths and these are shown in Figure 5. These results were obtained by taking 15° × 15° regions around sunspots. Using such a large region considerably dilutes their magnetic field, but also allows one to obtain results up to depths corresponding to the He ii ionization zone. Therefore, the magnetic field of this large region is lower than that of the active region itself. Inversion results for three such regions with surface magnetic field strengths of about 20, 68, and 80 Gauss are shown in Figure 5. The bottom panel shows that a larger change takes place in the boundary between the Hei and He ii ionization zones. The second reason is that the Gough (2013) analysis, based on Houdek & Gough (2007), assumed that the He ii can be modeled as a Gaussian with a finite width and that the magnetic fields change the amplitude of the Gaussian. However, as can be seen from Figure 5, while magnetic fields do make the He ii ionization-related dip in Γ₁ smaller, they also distort the shape of the dip. Thus the model used to associate the changes in amplitude to changes in Γ₁, and then the conversion of the change in Γ₁ needs to be changed. All these factors need to be accounted for to determine a change in the magnetic fields from the amplitude of the glitch signature.

Zoom In Zoom Out Reset image size

A detailed analysis carried out by Baldner & Basu (2008) revealed that between solar cycle 23 maximum and minimum, there were small but significant changes in the sound speed at the equator and at a latitude of 15°. The results were too noisy at higher latitudes to see anything, but at the low latitudes the changes were seen even below 0.9R_⊙. Our results indicate that changes also occur in a spherically symmetric sense.

We would like to thank the referee for comments that have helped us improve this Letter.