{'fraction_non_alphanumeric': 0.03082000507743082, 'fraction_numerical': 0.021934501142421934, 'mean_word_length': 4.259279038718291, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solar Cycle Related Changes in the Helium Ionization Zones of the Sun

Helioseismic data for solar cycles 23 and 24 have shown unequivocally that solar dynamics changes with solar activity. Changes in solar structure have been more difficult to detect. Basu&Mandel (2004) had claimed that the then available data revealed changes in the HeII ionization zone of the Sun. The amount of change, however, indicated the need for larger than expected changes in the magnetic fields. Now that helioseismic data spanning two solar cycles are available, we have redone the analysis using improved fitting techniques. We find that there is indeed a change in the region around the HeII ionization zone that is correlated with activity. Since the data sets now cover two solar cycles, the time variation is easily discernible.


INTRODUCTION
It has been known for about thirty-five 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), henceforth 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, Γ 1 , 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 HeII 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 Γ 1 can also occur in the presence of magnetic fields  because of changes in the local temperature. Thus a changing magnetic field in the in the layers around the HeII ionization zone will change the amplitude of the cbwatson@bu.edu 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. These results, and those of Verner et al. (2006) and Ballot et al. (2006), were however challenged by Gough (2013). Gough (2013) analyzed the potential significance of the results under the assumption that the change in the amplitude of the Heii 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 Heii 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 2 > in the vicinity of the Heii 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 on the result of every set is the same; we abandon that assumption and determine uncertainties for each set by boot-strapping. Another change 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 χ 2 minimum.

DATA, MODEL, AND FITTING TECHNIQUES
We use helioseismic data from three sources: (I) The ground-based Global Oscillation Network Group (GONG) (Hill et al. 1996), (II) the Michelson Doppler Imager (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., 3 GONG months) time series. There is an overlap of 72 days between different data sets. Data from MDI cover the period from May 1, 1996 to April 24, 2011, 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 andTapping 2013) 1 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: where ν n,ℓ is the frequency of a mode of radial order n and degree ℓ. The fourth differences are fitted to the form The parameters a 1 -a 4 define the overall smooth term, a 5 -a 8 , τ a and φ a define the oscillatory component due to the Heii ionization zone, and the terms b 1 -b 4 , τ b , and φ b define the oscillatory component due to the base of the convection zone. The frequencies of the two components, τ a and τ b , are approximately the acoustic depths of the Heii 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 Eq. 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-χ 2 fits, and this is where we differ from BM04. Eq. 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 χ 2 as where y is a column matrix with the data, y ′ the result of the model fit, and C is the covariance matrix. Like BM04, we fit the form in Eq. 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 HeII ionization zone, as well as the CZ base and consequently all three terms of Eq. 1 are fitted. Modes in the second set (the 'high-degree' set) sample only the HeII ionization zone and consequently only the first two terms of Eq. 2 apply. Given that the amplitude of the two oscillatory terms in Eq. 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 ν min and ν max to be the minimum and maximum frequency, we define < A > as 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 χ 2 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 errorcorrelation matrix, these can be seen in Fig. 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 obtainedthe 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.
In Fig. 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 χ 2 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. Figure 1. The difference between the average amplitude of the glitch signature between 2 and 3.5 mHz obtained in this work and the work of BM04 that did not take the complete error-correlation matrix into account. The gray line shows the no-difference case, while the blue line is a weighted least-squares fit to the differences. Note that for GONG data the differences show a slope with respect to the 10.7 cm flux in solar flux units (1 SFU = 10 4 Jansky = 10 −22 W m −2 Hz −1 ). This indicates that the slope obtained in this work will be different from those obtained earlier even if we restricted ourselves to the same data set. For MDI, there seems to be a constant difference, rather than a difference in slope. Figure 2. The difference between the average amplitude of the glitch signature between 2 and 3.5 mHz obtained between the BFGS and LM codes for the low-degree mode set for (a) GONG and (b) MDI and HMI data. The grey line in each panel marks the no-difference case. There is a small, but statistically insignificant differences in the slopes with respect to the 10.7 cm flux obtained with GONG data (blue line in Panel a) and HMI data (green line in Panel b). The slope for the MDI sets is slightly more significant. Note that if MDI and HMI are used together, as we do later, then the differences have a slope of 1.64 × 10 −4 ± 9.61 × 10 −5 . Figure 3. The average amplitude of the helium glitch signature plotted as a function of the 10.7 cm flux. The results were obtained with the BFGS code. In each panel, the cyan horizontal like is the weighted average of the points. The blue solid line is a weighted least-squares fit to the points with the 95% confidence limit plotted as the blue dashed lines. Panels (a), (c), and (e) are low-degree results, while panels (b), (d), and (f) show high-degree results. We show all the sets for GONG in panels (a) and (b) and only non-overlapping sets in panels (c) and (d).

RESULTS
The average amplitude of the glitch signature is plotted as a function of the 10.7 cm flux in Fig. 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 ± 4.67 × 10 −5 , −2.88 × 10 −4 ± 8.42 × 10 −5 and −2.23 × 10 −4 ± 9.61 × 10 −5 . 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.
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 Fig. 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.

DISCUSSION AND CONCLUSIONS
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 Heii 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- high-degree sets. The magenta line is a weighted average of the blue points which show results for non-overlapping GONG sets. The results are from the BFGS codes. The grey curve in the background is the 10.7 cm flux whose value can be read from the axis on the right. 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 Γ 1 caused by Heii 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 Γ 1 between the Hei and Heii ionization zones (shown in Fig. 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 Γ 1 profile, consequently we look instead to Γ 1 changes below sunspots of different strengths and these are shown in Fig. 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 Heii ionization zone. The magnetic field of this large region is therefore 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 the figure. The lower panel shows that a larger change takes place in the boundary between the Hei and Heii ionization zones. The second reason is that the Gough (2013) analysis, based on Houdek & Gough (2007), assumed that the Heii 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 Fig. 5, while magnetic fields do make the He ii ionization-related dip in Γ 1 smaller, they difference between the active regions and the quiet Sun. The active region inversion results are from . The magnetic activity index defining the magnetic field of the region studied is also indicated for each active region. The dotted vertical line marks the position of the minimum in Γ1 in the He ii ionization zone for the model, the dashed line marks the peak in Γ1 between the Hei and He ii ionization zones. also distort the shape of the dip. Thus the model used to associate the changes in amplitude to changes in Γ 1 , and then the conversion of the change in Γ 1 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.
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 the paper.