1 Introduction

The properties of near-surface convection are subject to slight changes over the course of the activity cycle (e.g., a decrease in granule size in phase with the solar cycle found by Macris et al., 1984; Muller, 1988, and a decrease of granular contrast with increasing level of magnetic activity, see Muller, Hanslmeier, and Saldaña-Muñoz, 2007). As solar p modes are stochastically driven by the acoustic noise generated by convective motion, the related p-mode parameters also vary over the solar activity cycle. However, most previous studies that investigated these activity-related changes of solar p-mode parameters, aside from mode frequencies, were either limited to unresolved solar time series and were thus restricted to low harmonic degrees \(l\lesssim3\) (e.g., Pallé, Régulo, and Roca Cortés, 1990; Elsworth et al., 1993; Chaplin et al., 1998; Howe et al., 2003; Jiménez-Reyes et al., 2004; Salabert et al., 2007; Broomhall, Pugh, and Nakariakov, 2015), or had a very short data baseline to work on (e.g., Jefferies et al., 1991 with 460 hr of data with a 54% duty cycle). Over the last years, little attention has been given to the activity-related changes of solar p-mode parameters for resolved data aside from mode frequencies and frequency splitting coefficients (the exceptions being Komm, Howe, and Hill, 2000b; Korzennik et al., 2013; Korzennik, 2017). Now, with two full solar sunspot cycles’ worth of data from the Global Oscillation Network Group (GONG) available, we examine the temporal variations of parameters of solar p modes of harmonic degrees \(l= 0\,\mbox{--}\,150\).

It is well established that mode widths, \(\Gamma\), which are inversely proportional to the lifetimes of the mode, vary in phase with the level of solar magnetic activity (e.g., Jefferies et al., 1991; Komm, Howe, and Hill, 2000a; Chaplin et al., 2000; Jiménez, Roca Cortés, and Jiménez-Reyes, 2002; Jiménez-Reyes et al., 2003, 2004; Salabert et al., 2007; Burtseva et al., 2009; Broomhall, Pugh, and Nakariakov, 2015). This is taken to indicate that the oscillations are damped by the presence of a magnetic field. The mode amplitudes, \(A\), are observed to be in antiphase with magnetic activity (e.g., Pallé, Régulo, and Roca Cortés, 1990; Elsworth et al., 1993; Chaplin et al., 2000; Komm, Howe, and Hill, 2000a; Jiménez, Roca Cortés, and Jiménez-Reyes, 2002; Jiménez-Reyes et al., 2003, 2004; Broomhall, Pugh, and Nakariakov, 2015). The magnitude of the change of \(\Gamma\) and \(A\) depends on both mode frequency and harmonic degree. Along with mode amplitudes, the mean squared velocity \(\langle v^{2}\rangle\) and energy \(E\) of the p modes change with the solar cycle, where highest mode velocities and energies are observed during solar minimum (Komm, Howe, and Hill, 2000b; Jiménez-Reyes et al., 2003; Salabert et al., 2007; Jiménez-Reyes et al., 2004). The energy supply rate to the modes, which is proportional to mode energy times mode width, has been measured not to change throughout the solar cycle for modes of low harmonic degree (Chaplin et al., 2000; Jiménez-Reyes et al., 2003, 2004; Broomhall, Pugh, and Nakariakov, 2015).

In this article, we do not consider mode frequencies. As these are more accessible than the parameters studied here, they have been subject of numerous studies in the past, and their behavior throughout the solar activity cycle is well documented for a wide range of harmonic degrees (e.g., Woodard and Noyes, 1985; Elsworth et al., 1990; Libbrecht and Woodard, 1990; Jiménez-Reyes et al., 1998; Chaplin et al., 2001; Salabert, García, and Turck-Chièze, 2015; Tripathy, Jain, and Hill, 2015; Broomhall, 2017). Mode frequencies are correlated with solar activity, being at their highest during times of strong activity. For very high frequency modes, around and above the acoustic cutoff frequency, this correlation turns into an anticorrelation: the frequencies of these modes are higher during times of weak activity (see, e.g., Woodard and Libbrecht, 1991; Howe et al., 2008; Rhodes et al., 2010). The magnitudes of the shifts of mode frequencies depend on both mode frequency and harmonic degree (see Basu, 2016 and references therein).

This article is structured as follows: The data and the corrections applied to them are described in Section 2. This includes correcting for the spatial masking of GONG Dopplergrams and azimuthal averaging (Section 2.2), and measures to account for the imperfect duty cycle (Section 2.3) and for spurious jumps in the mode parameters (Section 2.4). We present our results for the temporal variation of mode widths and amplitudes in Section 3 and focus on the temporal average of the mode amplitudes in Section 3.3. We proceed to consider the temporal variation of the quantities mean squared velocity, energy, and energy supply rate of the p modes (Section 4). A summary and discussion of our findings is given in Section 5.

2 Data and Method

2.1 Data Sets

In the present analysis, we use mode parameter data that were produced by the standard GONG pipeline (Anderson et al., 1990; Hill et al., 1996; Hill and Howe, 1998) from solar full-disk Dopplergrams.Footnote 1 The mode parameters were obtained by fitting symmetric Lorentzian profiles to the power spectra of 108-day-long datasets, where every third dataset was independent, i.e., they overlapped by 72 days. These time series were concatenations of three GONG months, with one GONG month being 36 days. This ensures a reasonable frequency resolution of the fitted power spectra and number of independent data points. The power spectra were computed and fit for all harmonic degrees \(l\) and azimuthal orders \(-l\le m\le l\) up to \(l=150\). The model that was fit to each peak in the power spectra consisted of a linear two-parameter background and a Lorentzian profile, which depends on the parameter width \(\Gamma\), amplitude \(A\), and frequency \(\nu\). Mode frequencies and the background parameters are not considered in the following. The GONG peak-fitting algorithm fits symmetric Lorentzian profiles to the spectra (Hill et al., 1996). Since the mode peaks are known to be asymmetric (e.g., Nigam et al., 1998), this might introduce small, temporally varying systematics to the mode parameter analyzed, as the asymmetry changes over the solar cycle (Jiménez-Reyes et al., 2007; Korzennik, 2013; Howe et al., 2015).

2.2 Azimuthal Averaging and Correction for Spatial Masking

In full-disk images of the Sun, pixels close to the solar limb are subject to several unwanted effects. First, because GONG measures line-of-sight velocities, projection effects increase toward the solar limb. Second, the resolution of each pixel, measured in distance on the solar surface, increases toward the limb. A spatial mask is applied in the GONG pipeline to cut away these unwanted pixels. This leads to a suppression of mode amplitudes depending on the ratio of azimuthal order and harmonic degree \(|m/l|\): While spherical harmonics with \(m=0\) oscillate over all latitudes (with a temporal dependence of \(\mbox{e}^{-\mathrm{i}\omega t}\), where \(\omega\) is the angular frequency), spherical harmonics with azimuthal order \(m\ne0\) are confined to latitudes closer to the equator. The confinement to low latitudes becomes more pronounced for higher harmonic degree and higher \(|m/l|\). Therefore, by masking out regions near the solar limb, more mode power is cut away from modes with low \(|m/l|\) than from those with higher \(|m/l|\). In contrast, the damping width of a mode in a power spectrum, representing the lifetime of the mode, is negligibly affected by the masking.

Figure 1 shows the outputs of the GONG pipeline for mode width \(\Gamma\), amplitude \(A\), and the products \(\Gamma\cdot A\) and \(\Gamma^{2}\cdot A\) of the mode \((n,l)=(10,50)\) from GONG month 5, where \(n\) is the radial order. In the top left panel, the mode widths of the azimuthal components are shown as a function of \(m/l\), but no visible dependence is observed. The solid red line is the error-weighted mean of the data. In the top right panel, the black data points are the measured values of the mode amplitudes. As discussed above, mode amplitudes and products that include mode amplitudes show a dependence on \(m/l\), which is introduced by the spatial mask. It is accounted for by fitting a polynomial proportional to \((m/l)^{2k}\) with \(k=0,1,2\) to the data. This polynomial function is empirically determined (see Komm, Howe, and Hill, 2000a). The resulting fit is shown as a solid black line in the top right and two lower panels of Figure 1. The dependence is then removed, where the values at \(|m/l|=1\) are kept constant. The error-weighted mean of the corrected data points is adopted as the representative value of the mode parameter (shown as a solid red line). In order to obtain meaningful weighted mean values, in the subsequent analysis, only those multiplets are included for which at least one-third of the azimuthal components are fit by the GONG pipeline. This minimum number of fit components is set to two-thirds for modes with low harmonic degree \(l\le10\), to account for their small number of azimuthal components. In the following, we avoid the term multiplet, as only azimuthal averages are considered, and we use the term mode instead.

Figure 1
figure 1

Correction for the effects of the spatial masking for the mode \((n,l)=(10,50)\) in GONG month 5. Black data points are measurement values. Solid black lines are polynomial fits to the data to account for the effect of the spatial masking. Red data points are the corrected values. Solid red lines indicate the weighted mean of the corrected data.

2.3 Correction for Window Function

The temporal fill factor of GONG is lower than unity. This leads to a redistribution of the power from the mode peak to neighboring frequency bins and into side lobes of the main peak, which in turn results in a diminished mode amplitude and an increased mode width. The side lobes are caused by repeating gaps in the data. These side lobes, especially those caused by daily gaps, are considerably suppressed for a network distributed around the globe as GONG is (Hill et al., 1996; Leibacher, 1999). A machine-readable table with the fill factor of each 36-day GONG month is available at https://gong2.nso.edu/fill.txt . The specific structure of the gaps does not have to be taken into account as the fill factor is rather high, with values between 69% and 93%. To account for the impact of gaps, we adopted the approach described in Komm, Howe, and Hill (2000b), who showed that a linear regression is sufficient to correct for the effect of the temporal window function of GONG data.

In Figure 2, the mode width \(\Gamma\), mode amplitude \(A\), mode area \(\Gamma\cdot A\), and mode width times mode area \(\Gamma ^{2}\cdot A\) of the mode \((n,l)=(10,50)\) are shown as a function of fill factor. The products \(\Gamma\cdot A\) and \(\Gamma^{2}\cdot A\) are used in Section 4 to calculate the quantities of mean squared velocity, mode energy, and energy supply rate of the modes. The data are normalized to the value of the linear regression at \(\mbox{fill}=1\). To obtain the corrected parameter values, the resulting fit of a linear function, shown as a red line in each panel of Figure 2, is first subtracted from the data and then its value at \(\mbox{fill}=1\) is added. As previously demonstrated by, e.g., Komm, Howe, and Hill (2000a), mode widths increase with decreasing fill factor and amplitudes decrease with decreasing fill factor. It should be noted that the background amplitude increases with decreasing fill factor. This indicates that power from the mode peaks is indeed distributed into the background.

Figure 2
figure 2

Correction of the effect of the temporal window function on mode parameters of mode \((n,l)=(10,50)\). Parameter values are shown as a function of the temporal fill factor of the individual GONG months. Linear fits to the data are shown as red lines. Data are normalized to the value of the linear fit at \(\mbox{fill}=1\).

2.4 Correction for Jumps in Sensitivity

There are two jumps in the mode amplitudes that we corrected for with an empirical correction factor. The first jump around GONG month 60 is due to a camera upgrade. The second jump around month 100 is of uncertain origin. The correction factor \(C\) is given in Equation 12 in Appendix A. More details about the applied corrections and the jumps are given in Appendix A.

2.5 Proxy for Solar Activity

Many quantities connected to the Sun are known to vary with the solar activity cycle, e.g., the sunspot number and the emission in the Ca ii H & K lines. Here, we use the solar radio flux \(F_{10.7}\) as a proxy for the level of solar activity.Footnote 2 The \(F_{10.7}\) is the total emission on the solar disk at a wavelength of 10.7 cm integrated over one hour. \(F_{10.7}\) is measured in solar flux units sfu, where \(1~\mbox{sfu} = 10^{-22}~\mbox{W}\,\mbox{m}^{-2}\,\mbox{Hz}^{-1}\). The \(F_{10.7}\) index is the averaged \(F_{10.7}\) flux scaled for 1 AU. It is a proxy for the level of activity in the upper chromosphere and the lower corona (Tapping, 2013). A comparison of different proxies of solar magnetic activity and their correlation with the activity-related frequency shifts of solar p modes over the last three solar cycles can be found in Broomhall and Nakariakov (2015).

3 Temporal Variation of Mode Parameters

3.1 Mode Widths

To obtain clear results on their temporal variation, mode parameters have to be averaged over many modes. From here on, only modes that are present at all time samples are considered. The analysis of both mode widths and mode amplitudes is restricted to modes with frequencies \(1500~\upmu\mbox{Hz}\le\nu\le4500~\upmu\mbox{Hz}\) and harmonic degree \(0\le l \le150\). This results in a set of 1275 mode widths and 1272 mode amplitudes. Different error flags for bad mode widths and amplitudes result in the small difference in the number of modes in the sets. As both mode widths and amplitudes are dependent on mode frequency and harmonic degree, they are normalized to the mean over all time samples for each mode individually. Next, the normalized parameters of all included modes are averaged. The result for the observed mode widths \(\Gamma_{\text{obs, mean}}\) is shown in the first panel of the top row of Figure 3. Here, the correction for the effect of the fill factor has not yet been applied. The need for such a correction can be seen from the middle panel of the top row, where the same data are plotted as a function of fill factor. The linear dependence on fill is apparent. In the right panel, the mode widths are shown as a function of the \(F_{10.7}\) index. The red data points in all panels of Figure 3 indicate time samples with higher than median level of activity as computed from the \(F_{10.7}\) flux.

Figure 3
figure 3

Normalized mode widths averaged over all modes as a function of time (left column), fill factor (middle column), and of the magnetic activity proxy \(F_{10.7}\) index (right column). The top row shows the raw data, the middle row is normalized for the effects of the fill factor, and the bottom row is corrected for fill and for change in apparent radius of the Sun.

In the middle row of Figure 3, the mode widths \(\Gamma_{\mathrm{fill}=1, \mathrm{mean}}\) are shown after the correction for the temporal window function. A second correction is required to account for the apparent change in solar radius throughout the year. This change affects the measured values of oscillation parameters as the spatial resolution of the Dopplergrams changes with it. We applied a linear correction to the mode parameters as a function of apparent solar radius to correct for this. The bottom row shows the mode widths \(\Gamma_{\mathrm{fill}=1, \mathrm{radius}, \mathrm{mean}}\) after both the effect of the fill factor and the change of the apparent solar radius have been removed.

For better visibility, the left panel of the bottom row of Figure 3 is shown again in the top panel of Figure 5. As can be seen, the uncertainties on the normalized variation of the mode widths are on the order of a few tenths of a percent. The error bars given here are computed as the standard error of the mean. The solid black line is the one-year running average. The red line is the \(F_{10.7}\) index. It is boxcar smoothed over one year and scaled to match the extrema of the one-year boxcar-smoothed variation of the mode widths. The correlation coefficient (Spearman’s rank correlation coefficient \(\rho\)) between the variation of the width and the \(F_{10.7}\) index (both unsmoothed) is \(\rho=0.62\) with \(p= 2\cdot10^{-9}\). Only independent data points were used to calculate this. The number of independent data points is 74, i.e., every third GONG 108-day dataset. We used Spearman’s \(\rho\) to asses the level of correlation of the two quantities as the relation between mode widths and activity is not strictly linear: As can be seen from the right panel in the bottom row in Figure 3, the variation of mode widths increases with the level of the \(F_{10.7}\) index. However, for values of \(F_{10.7}\gtrsim130\), mode widths appear to be in saturation. The largest discrepancy between the \(F_{10.7}\) index and the mode widths can be seen during the time of the camera upgrade in the GONG network around month 60 (in 2001). Averaged over all modes, the minimum-to-maximum variation of the one-year boxcar-smoothed variation of mode widths is \(11.5\pm0.2\%\).

In Figure 6, the normalized variation of averages of mode widths (after the corrections for fill and apparent solar radius variation have been applied) for modes of different frequency and harmonic degree ranges are shown as a function of time. Table 3 in Appendix B holds information on the number of modes used in each panel of Figure 6, the minima and maxima of the parameter variation, the mean uncertainty of each data point, and the correlation of mode widths with the \(F_{10.7}\) index.

The minimum-to-maximum variation of the mode widths over the two observed solar cycles is dependent on mode frequency and harmonic degree. The largest fractional variation is found for modes with \(2400~\upmu\mbox{Hz}\le\nu\le3300~\upmu\mbox{Hz}\) and \(31\le l \le60\) with a minimum-to-maximum variation of \(26.6\pm0.3\%\). The smallest variation of mode widths is found for the modes in the \(1500~\upmu\mbox{Hz}\le\nu\le2400~\upmu\mbox{Hz}\) and \(101\le l \le150\) parameter regime with a minimum-to-maximum variation of \(5.5\pm0.2\%\). To suppress the impact of outliers on these minimum-to-maximum variations, they were calculated for the one-year boxcar-smoothed values of the mode widths. The middle frequency range with \(2400~\upmu\mbox{Hz}\le\nu\le3300~\upmu\mbox{Hz}\) shows the largest minimum-to-maximum variation over the solar cycle for each subset of modes of the same range of harmonic degrees. The exception to this is the subset of modes with \(101\le l \le150\), for which the variation over the solar cycle is largest for modes in the high-frequency range.

The last two columns of Table 3 give the correlation between the variation of mode widths and the \(F_{10.7}\) index as well as the associated two-sided significance value. The highest level of correlation is found for modes with \(1500~\upmu\mbox{Hz}\le\nu \le2400~\upmu\mbox{Hz}\) and \(0\le l \le30\), with a Spearman’s rank correlation coefficient of \(\rho=0.84\) and a two-sided significance of \(p<10^{-10}\). The lowest correlation is found for modes with \(1500~\upmu\mbox{Hz}\le\nu\le2400~\upmu\mbox{Hz}\) and \(31\le l \le60\), with \(\rho=0.31\) and \(p\approx0.01\). Within each range of harmonic degrees, the correlation coefficient increases with mode frequency. The exception to this are modes with \(0\le l \le30\), where modes in the middle frequency range show the highest level of correlation with the \(F_{10.7}\) index.

3.2 Mode Amplitudes

In Figure 4 the mode amplitudes \(A\) are shown as functions of time, fill, and \(F_{10.7}\) index (columns) and for three levels of correction (rows), as was previously described for the mode widths (Figure 3). The effect of the correction for the change in apparent size of the Sun is more pronounced for mode amplitudes than it is for mode widths. As can be seen in the left panel of the middle row, there is a distinct yearly variation in the mode amplitudes \(A_{\mathrm{fill}=1, \mathrm{mean}}\) after the fill correction has been applied. This is largely removed in the left panel of the bottom row. The change to the mode amplitudes due to this correction is most obvious through months 110 – 150, which corresponds to the activity minimum between Cycles 23 and 24. In the right panel of the bottom row of Figure 4, mode amplitudes \(A_{\mathrm{fill}=1, \mathrm{radius}, \mathrm{mean}}\) are shown as a function of the \(F_{10.7}\) index. A clear anticorrelation between the two quantities is visible. For levels of activity up to \(F_{10.7}\approx130\), the relation appears linear. For higher levels of activity, there is little to no change in mode amplitudes, similar to the behavior seen for mode widths (Figure 3).

Figure 4
figure 4

Same as Figure 3, but for mode amplitudes.

The normalized, averaged, and corrected temporal variation for mode amplitudes is shown in the lower panel of Figure 5. The last row of Table 4 in Appendix B gives more detailed information about the number of modes within each parameter range and the error bars presented in this plot. The error bars in Figure 5 were computed as the standard error of the mean. The solid lines are the data (black) and \(F_{10.7}\) index (red) smoothed with a one-year boxcar. The \(F_{10.7}\) index is scaled to match the minimum-to-maximum variation of the mode amplitudes and flipped for better appreciation of the anticorrelation of the two quantities. Averaged over all modes, the minimum-to-maximum variation is \(17.3\pm 0.2\%\). The correlation coefficient (Spearman’s rank correlation coefficient) between amplitudes and the \(F_{10.7}\) index (both unsmoothed) is \(\rho=-0.91\) with \(p<10^{-10}\). Only independent data points were used in the computation of the correlation coefficient.

Figure 5
figure 5

Top panel: Normalized mode widths averaged over all modes (black data points). Error bars indicate the error of the mean. The solid black line is the one-year average of the mode widths. The red line is the \(F_{10.7}\) solar radio flux as a function of time smoothed by boxcar with a width of one year. It is scaled to match the variation of the mode widths. Bottom panel: Same as top panel, but for mode amplitudes. The \(F_{10.7}\) radio flux was flipped for better appreciation of the anticorrelation of the change in mode amplitudes and the level of activity.

The equivalent of Figure 6, where the normalized variation of mode widths is shown for different mode parameter regimes, is presented in Figure 7 for mode amplitudes. Detailed information on the individual panels is given in Table 4 in Appendix B. The largest variation is found for modes with \(2400~\upmu\mbox{Hz}\le\nu\le3300~\upmu\mbox{Hz}\) and \(61\le l \le100\) with a minimum-to-maximum variation of \(27.4\pm0.4\%\). Modes with frequencies \(1500~\upmu\mbox{Hz}\le\nu\le2400~\upmu\mbox{Hz}\) and harmonic degrees \(0\le l \le30\) exhibit the smallest variation over the 22 years of data, with a fractional change of \(11.6\pm0.5\%\). The anticorrelation between the level of activity and the change in mode amplitudes is highest for modes with \(2400~\upmu\mbox{Hz}\le\nu \le3300~\upmu\mbox{Hz}\) and \(0\le l \le30\) with a rank correlation of \(\rho=-0.94\) and \(p<10^{-10}\). For the three ranges of harmonic degrees (\(l\) between 0 – 30, 31 – 60, and 61 – 100) the largest variation is found in the middle frequency range (\(2400\,\mbox{--}\,3300~\upmu\mbox{Hz}\)). For modes with \(101\le l \le150\), this is observed for the high-frequency modes.

Figure 6
figure 6

Temporal variation of mode widths for different ranges of harmonic degrees (rows, harmonic degree range indicated to the right of the fourth column) and mode frequencies (columns, frequency range indicated above the first row). Widths are normalized to the mean for each mode and then averaged over all modes in the respective range of frequency and degree. Months with higher than median \(F_{10.7}\) solar radio flux are highlighted by red points. The one-year average is shown by the solid black line. Levels of 1.1 and 0.9 of the mean are indicated by gray dashed lines.

Figure 7
figure 7

Same as Figure 6, but for mode amplitudes.

3.3 Frequency Distribution of Mode Amplitudes

In the top panel of Figure 8, the amplitudes of all modes with \(2\le l \le150\), which are present in at least 50% of the GONG months, are shown as a function of mode frequency on a logarithmic ordinate divided by \(10^{4}\). Here, the amplitudes interpolated to \(\text{fill}=1\) were used. Any variation due to the residual change in apparent solar radius was averaged out as the data cover 22 years. By lowering the presence rate to 50%, more modes, especially of \(l\gtrsim 100\), were included and the mode amplitudes could be investigated for separate ranges of harmonic degrees (see the bottom panel of Figure 8 for an example). This increased the total number of included modes from 1272 to 1543.

Figure 8
figure 8

Mode amplitudes as a function of mode frequency. Amplitudes are averaged over all GONG months after correction for the fill factor. The red curve is an asymmetric Voigt profile fit to the data. Top panel: Amplitudes of modes with \(2\le l \le150\). Bottom panel: Amplitudes of modes with \(31\le l \le40\).

The solid red lines in Figure 8 are fits of asymmetric Voigt profiles to the mode amplitudes. This profile is described by

$$\begin{aligned} V(\nu) = \mathcal{A}\cdot \biggl(b+ a\cdot \int\limits _{-\infty} ^{\infty}G(x) L(\nu-x)\,\mathrm{d}x \biggr), \end{aligned}$$
(1)

where \(b\) is an offset, \(a\) is a factor to adjust the height of the profile, and the integral extends over the entire frequency axis with frequencies \(x\). The Gaussian function \(G(\nu)\) and the Lorentz function \(L(\nu)\) in Equation 1 are given by

$$\begin{aligned} G(\nu) &=\frac{\exp(-\nu/2\sigma^{2})}{\sigma\sqrt{2\pi}}, \end{aligned}$$
(2)
$$\begin{aligned} L(\nu) &=\frac{\gamma}{\pi (\nu^{2}+\gamma^{2} )}, \end{aligned}$$
(3)

with the standard deviation of the Gaussian \(\sigma\) and the half-width at half-maximum of the Lorentzian \(\gamma\). The asymmetry is introduced with the function

$$\begin{aligned} \mathcal{A} = \frac{1}{\pi} \bigl(\tan^{-1} \bigl(S \cdot (\nu- \nu_{\text{max}} )/\Sigma \bigr) + 0.5 \bigr), \end{aligned}$$
(4)

where \(S\) is the asymmetry parameter, \(\nu_{\text{max}}\) is the frequency of the maximum of the symmetric Voigt profile, and the full width at half-maximum of the Voigt profile, \(\Sigma\), can be approximated by (Whiting, 1968; Olivero and Longbothum, 1977)

$$\begin{aligned} \Sigma= 1.0692 \gamma+ \sqrt{0.86639 \gamma^{2} + 8 \ln(2)\sigma ^{2}}. \end{aligned}$$
(5)

The resulting fit parameters for amplitudes of modes with \(2\le l \le150\) are given in Table 1. For the frequency of maximum amplitude, we find \(\nu_{\text{max}}=3079.76\pm0.17~\upmu\mbox{Hz}\). The width of the Voigt profile is \(\Sigma= 611.8\pm0.5~\upmu\mbox{Hz}\). The parameter \(S=-0.100\pm0.002\) indicates that the distribution is slightly left-tailed (left-skewed, right-leaning). The \(\chi^{2}_{\text{red}}\) of the fit is 32.8. This rather high \(\chi^{2}_{\text{red}}\) value is due to the spread in the distribution of the mode amplitudes for modes of different harmonic degree and the small error bars. Other profiles were tested (pure Gaussian, pure Lorentzian, both Gaussian and Lorentzian including asymmetry), but the asymmetric Voigt yielded the best \(\chi^{2}_{ \text{red}}\). We excluded modes with \(l=0,1\) because their amplitudes are less well defined, and including them here would increase the \(\chi^{2}_{\text{red}}\) to \({\approx}\,135\).

Table 1 Fit parameters of the frequency distribution of mode amplitudes for modes with \(2\le l \le150\).

It should be noted that the maximum of the mode amplitudes \(\nu_{\text{max}}\) presented here is not exactly what is often referred to as \(\nu_{\text{max}}\) in the literature. It is usually measured from Sun-as-a-star data, which only include low harmonic degrees up to \(l\lesssim4\) (Broomhall et al., 2009). Here, \(\nu_{\text{max}}\) is the maximum of the mode amplitudes in the frequency spectrum of the spherical harmonic transform of the GONG Doppler velocity data for all harmonic degrees up to \(l=150\). Typical values for \(\nu_{\text{max}}\) include \(3050~\upmu\mbox{Hz}\) by Kjeldsen and Bedding (1995), \(3104~\upmu\mbox{Hz}\) as a calibrated value to obtain better results from asteroseismic scaling relations by Mosser et al. (2013), and \(3120~\upmu\mbox{Hz}\) from SOHO/VIRGO data by Kallinger et al. (2010). If \(\nu_{\text{max}}\) is measured from the velocity of the modes (as we discuss in Section 4, the mean velocity of the modes is proportional to mode width \(\Gamma\) times mode amplitude \(A\)), the damping widths \(\Gamma\) are included in the measured quantity. As mode widths change with mode frequency and harmonic degree (see next section), the maxima in mode amplitude \(A\) and mean velocity of the modes (proportional to \(\Gamma\cdot A\)) occur for slightly different frequencies.

We also investigated the mode amplitudes for smaller ranges of harmonic degrees. For this, we separated the modes into groups of typically ten harmonic degrees. As an example, we show the amplitudes of modes with \(31\le l \le40\) (black data points) and the fit of an asymmetric Voigt profile (solid red line) in the lower panel of Figure 8. As can be seen from Figure 8, the frequency dependence of the amplitudes of modes of similar harmonic degree is well captured by the asymmetric Voigt profile. The resulting fit parameters for this and more ranges of harmonic degrees are presented in Table 5 in Appendix C. There, we also show figures of the change of the fit parameters with harmonic degree. Except for the width of the Voigt profile \(\Sigma\), all parameters show systematic variations with mode degree. The frequency of maximum amplitude \(\nu_{\text{max}}\) has a maximum value for intermediate-degree modes and decreases substantially at low and high degrees (below \(l=30\) and above \(l=90\)). The skewness \(S\) is anticorrelated with \(\nu_{\text{max}}\). The width \(\Sigma\) is approximately constant, except for the band of highest harmonic degrees considered. Amplitudes first increase with harmonic degree, reach a maximum around \(l\approx40\), and then decrease again.

4 Results for Physical Quantities

The amplitude \(A\) is given as power per frequency bin \(\mbox{m}^{2}\,\mbox{s}^{-2}\,\mbox{Hz}^{-1}\). Hence, calculating the product of mode width and mode amplitude \(\Gamma_{nl}\cdot A_{nl}\) has the unit of squared velocity. The mean squared velocity of the modes can be calculated as (Goldreich, Murray, and Kumar, 1994)

$$\begin{aligned} \bigl\langle v^{2}_{nl}\bigr\rangle = \frac{\pi}{2} C_{\text{vis}} \Gamma_{nl} A_{nl}. \end{aligned}$$
(6)

Symmetric Lorentzian profiles are fitted to the peaks in the power spectrum by the GONG pipeline. This is taken into account by the scaling factor \(\pi/2\). The quantity \(\Gamma_{nl}\cdot A_{nl}\) is referred to as the mode area, i.e., the area below the fitted peak in the spectrum (Komm, Howe, and Hill, 2000b). The factor \(C_{\text{vis}}=3.33\) corrects for the reduced visibility of modes due to leakage effects (Hill and Howe, 1998).

Because of the different cavities in which modes of different harmonic degree and frequency propagate, the fraction of the mass of the Sun that is affected by different modes varies (see, e.g., Basu, 2016). The mode mass \(M_{nl}\) is given by

$$\begin{aligned} M_{nl} = M_{\odot} I_{nl}, \end{aligned}$$
(7)

where \(M_{\odot}\) is the solar mass and the mode inertia is calculated as (Christensen-Dalsgaard and Berthomieu, 1991)

$$\begin{aligned} I_{nl} = \frac{4\pi\int_{0}^{R_{\odot}} \rho(r) ({\xi^{r}_{nl}(r)} ^{2} + l(l+1){\xi^{h}_{nl}(r)}^{2} )r^{2}\,\mathrm{d}r}{M_{\odot } ({\xi^{r}_{nl}(R_{\odot})}^{2} + l(l+1){\xi^{h}_{nl}(R_{ \odot})}^{2} )} \,. \end{aligned}$$
(8)

Here, \(\xi^{r}_{nl}\) and \(\xi^{h}_{nl}\) are the radial and horizontal displacement eigenfunctions associated with the oscillation, \(\xi_{nl}^{r}(R_{\odot})\) and \(\xi_{nl}^{h}(R_{\odot})\) are their values at the photospheric radius, and \(\rho\) is the mass density. We calculated mode inertias from the eigenmodes of the standard solar model ‘Model S’ (Christensen-Dalsgaard et al., 1996).

The energy that is stored in the modes can be calculated as the product of mode mass (Equation 7) and mean squared velocity (Equation 6):

$$\begin{aligned} E_{nl} = M_{nl} \bigl\langle v^{2}_{nl} \bigr\rangle . \end{aligned}$$
(9)

This is the total mode energy, i.e., the sum of kinetic and potential energy (Goldreich, Murray, and Kumar, 1994). The rate at which energy is supplied to the modes can be calculated by the product of the energy in the modes and the radian damping rate of the modes \(2\pi\Gamma_{nl}\) (Kumar, Franklin, and Goldreich, 1988; Goldreich, Murray, and Kumar, 1994):

$$\begin{aligned} \frac{\mathrm{d}E_{nl}}{\mathrm{d}t} = 2\pi E_{nl} \Gamma_{nl} = \pi ^{2} C_{\text{vis}} M_{nl} A_{nl} \Gamma_{nl}^{2}, \end{aligned}$$
(10)

which is proportional to the product of squared mode width \(\Gamma _{nl}^{2}\) and mode amplitude \(A_{nl}\).

In Figure 9 the mode width (damping rate) \(\Gamma\), mean squared velocity \(\langle v^{2}_{nl}\rangle\), mode energy \(E\), and energy supply rate \(\mathrm{d}E/\mathrm{d}t\) of modes with harmonic degree \(l>10\) are shown on logarithmic scales as functions of mode frequency. Modes with \(l\le10\) are excluded here to further reduce the scatter in these plots. The parameters of these modes are less well defined owing to their small number of azimuthal components, see Section 2.2. For each mode, the values for \(\Gamma\), \(A\), \(\Gamma\cdot A\), and \(\Gamma^{2} \cdot A\) were calculated as the mean over all time samples after correction for the fill factor. Any effects from the activity cycle or residual yearly variations were averaged out as the dataset spans two complete Schwabe cycles. The conversion of these quantities into physical units was then made according to Equations 6, 9, and 10.

Figure 9
figure 9

Mode widths \(\Gamma\), mean squared velocity \(\langle v^{2}\rangle\), mode energies \(E\), and energy supply rates \(\mathrm {d}E/\mathrm{d}t\) as functions of frequency. All quantities are averaged over all GONG months.

The mode widths (upper left panel in Figure 9) increase from \(0.3\,\mbox{--}\,0.5~\upmu\mbox{Hz}\) at low mode frequencies \({\approx}\,1500~\upmu\mbox{Hz}\) to values between \(1\,\mbox{--}\,2~\upmu\mbox{Hz}\) for mode frequencies between \(2400\,\mbox{--}\,3000~\upmu\mbox{Hz}\). For higher mode frequencies, mode widths increase again, reaching \(\Gamma\approx 10~\upmu\mbox{Hz}\) for the highest mode frequencies. The mean squared velocity (upper right panel) increases monotonically with mode frequency until it reaches its maximum for modes with frequencies of \({\approx}\,3200~\upmu\mbox{Hz}\). The maximum mean velocity value is \(v\approx 37~\mbox{cm}\,\mbox{s}^{-1}\). The mode widths, mean squared velocity, and the energy supply rate (lower right panel) show ridges of modes with equal radial order. Different ridges are slightly offset from one another, but show the same behavior with frequency. The ridge structure is not apparent in the mode energies (lower left panel). The difference between mode energy \(E\) and the other three quantities is the inclusion of mode inertia, see Equations 6, 7, 9, and 10.

A normalized mode inertia can be calculated by computing the ratio of the mode inertia \(I_{nl}\) to the inertia of radial modes \(\overline{I _{n0}}\) interpolated to the frequency of the mode \(\nu_{nl}\) (Christensen-Dalsgaard et al., 1996; Aerts, Christensen-Dalsgaard, and Kurtz, 2010):

$$\begin{aligned} Q_{nl} = \frac{I_{nl}}{\overline{I_{n0}} (\nu_{nl} )}. \end{aligned}$$
(11)

The result of multiplying mode widths, mean squared velocities, and energy supply rates with \(Q_{nl}\) is shown in Figure 10. The overall frequency dependence of the quantities is unchanged. However, ridges of different radial orders are now largely collapsed into one. This can best be observed for \(\langle v^{2}_{nl}\rangle\) (top right panel), which is only composed of one thin line of data points without any residual ridge structure. The steep gradient in mode inertia at low frequencies is not perfectly represented by scaling \(Q_{nl}\) with \(\overline{I_{n0}}\). Some of the residual scatter in mode widths, energies, and energy supply rate at low mode frequencies is due to this. Different scalings, e.g., with the interpolated inertia of \(l=50\) modes, give somewhat different results (Komm, Howe, and Hill, 2000b). However, the scatter in the scaled quantities is never completely removed. The remaining dependence observed in mode widths across all frequencies may hold information on a degree dependence of the involved damping mechanisms.

Figure 10
figure 10

Same as Figure 9, but normalized for mode inertia with Equation 11. Mode energy \(E\) is unchanged.

To be able to appreciate the temporal change of the four quantities (mode width, mean squared velocity, mode energy, and mode energy supply rate), the parameters of many modes have to be averaged. As the quantities cover two to three orders of magnitude for the investigated set of modes, the frequency range of averaging has to be restricted. Otherwise, values that differ by orders of magnitude would contribute to the average. In Table 2, the frequency ranges of the modes that were taken for the averaging of the four quantities and the number of modes within these frequency ranges are given. The frequency ranges were chosen to select as many modes as possible with similar parameter values (within about a factor of two).

Table 2 Frequency ranges for averaging of physical quantities, number of modes in this frequency range, and correlation coefficients with the \(F_{10.7}\) index.

The result of this averaging is presented in Figure 11 as a function of time. The red data points indicate times of higher than median level of activity, and the solid black lines are the data smoothed with a one-year boxcar. The last two rows of Table 2 give the rank correlation of independent data points between the four quantities and the \(F_{10.7}\) index as well as the two-sided significance value. As before, the mode widths are correlated with the level of magnetic activity with \(\rho=0.69\). The temporal variations in mean squared velocity and mode energy are both highly anticorrelated with the level of solar activity with a value of \(\rho=-0.88\). The energy supply rate is not correlated with solar activity for the set of modes investigated here as \(\rho=-0.01\) and \(p=0.88\).

Figure 11
figure 11

Mode widths \(\Gamma\), mean squared velocity \(\langle v^{2}\rangle\), mode energies \(E\), and energy supply rates \(\mathrm {d}E/\mathrm{d}t\) as functions of time. Averaged over all modes in the frequency ranges given in Table 2.

5 Summary and Discussion

This study is the first to investigate the activity-related p-mode parameter changes of resolved data covering a complete solar magnetic cycle. We analyzed mode parameter data from the ground-based GONG network. Systematic effects due to the spatial masking of the GONG full-disk solar images, the imperfect temporal sampling, and a yearly modulation due to changes in the apparent radius of the Sun were corrected for. Measurements of mode widths and mode amplitudes are generally less robust than measurements of mode frequencies. Thus, averages over azimuthal orders as well as modes of different ranges of frequencies and harmonic degrees have to be performed in order to reach meaningful results for the temporal variation of these mode parameters.

Little attention has been given to mode widths and amplitudes (for resolved solar data; Komm, Howe, and Hill, 2000b; Korzennik et al., 2013; Korzennik, 2017 are the most notable exceptions), let alone their temporal variation over the solar cycle. The reason for this is that even without incorporating solar cycle variations, damping widths and amplitudes cannot be accurately modeled with the existing theory of convection and mode excitation (see, e.g., Houdek et al., 1999; Houdek and Dupret, 2015; Basu, 2016). Thanks to the long, uninterrupted GONG time series that is now available, we were now able to investigate the variation of the mode widths and amplitudes as functions of the level of solar activity for different subsets of mode frequencies and harmonic degrees. When the last studies similar to the present one were published, about three years’ (Komm, Howe, and Hill, 2000a) and four years’ (Komm, Howe, and Hill, 2000b) worth of GONG data had been recorded.

For the mode widths, which were averaged over 1275 modes, we found a variation of \(11.5\pm0.2\%\) between the minimum and maximum level of activity in the investigated time period. Mode amplitudes varied by \(17.3\pm0.2\%\) over the same time (averaged over 1272 modes). Overall, the results from previous analyses of the variation of mode widths and amplitudes with the level of solar activity by, e.g., Komm, Howe, and Hill (2000a, 2000b), and Broomhall, Pugh, and Nakariakov (2015) could be confirmed. However, we find a larger variation of mode widths and amplitudes than was expected by Komm, Howe, and Hill (2000a). By extrapolating the fractional change per Gauss of widths and amplitudes to the full minimum-to-maximum change of activity over a solar cycle, they expected widths to vary by 7% and amplitudes to vary by 16% averaged over all modes. Thus, the changes of mode widths we report here are about two-thirds larger than expected by Komm, Howe, and Hill (2000a), while the changes of mode amplitudes agree to within a few percent. The smaller estimated variation of mode widths by Komm, Howe, and Hill (2000a) was most likely due to the shorter time series and differences in which modes they were included.

Like Komm, Howe, and Hill (2000a), we also investigated the change of mode widths and amplitudes for subsets of mode frequencies and harmonic degrees. While they averaged over \(100~\upmu\mbox{Hz}\) or five harmonic degrees (their Figures 10 and 11), we did this for larger ranges in frequency and harmonic degree, see Figures 6 – 7 and Tables 3 – 4. We can confirm their findings that the change of mode widths and amplitudes with activity is largest for modes around the frequency of maximum amplitude. We also can confirm that the change is more strongly dependent on mode frequency than it is on harmonic degree: for example, the amplitude of the fractional variation of mode amplitudes for modes with \(61\le l \le100\) is 12.9%, 27.4%, and 19.6% for the low-, mid-, and high-frequency ranges, respectively. Keeping the frequency range the same, the fractional variation of mode amplitudes for modes with \(2400\le\nu\le3300\) is 24.8%, 23.7%, 27.4%, and 18.4% for the four ranges of harmonic degrees in ascending order.

Table 3 Results from the normalized and averaged variation of mode width for different parameter ranges shown in Figure 6. Columns 1 – 4: Parameter ranges of the individual panels. Column 5: Number of modes included. Columns 6 – 7: Extrema of the one-year smoothed parameter variation Column 8: Typical error on the individual unsmoothed data point. Columns 9 – 10: Spearman’s correlation coefficient between independent data points of the unsmoothed parameter variation and the \(F_{10.7}\) index and the associated p value.
Table 4 Results from the normalized and averaged variation of mode amplitudes for different parameter ranges shown in Figure 7. Columns 1 – 4: Parameter ranges of the individual panels. Column 5: Number of modes included. Columns 6 – 7: Extrema of the one-year smoothed parameter variation. Column 8: Typical error on the individual unsmoothed data point. Columns 9 – 10: Spearman’s correlation coefficient between independent data points of the unsmoothed parameter variation and the \(F_{10.7}\) index and the associated p value.

The amplitude of the shifts of mode frequencies with the level of solar activity is known to increase with mode frequency (Jiménez-Reyes et al., 1998; Salabert, García, and Turck-Chièze, 2015). In contrast, this behavior is observed neither for mode widths nor for mode amplitudes. For these two parameters, the largest variation over the solar cycle is observed around the frequency of maximum amplitude \(\nu_{\text{max}}\), see Tables 3 and 4, which has been measured by Komm, Howe, and Hill (2000b), see their Table 3. This frequency-dependent variation is in agreement with the theoretical calculations by Houdek et al. (2001), who found that mode damping widths for frequencies between \(2500\,\mbox{--}\,3000~\upmu\mbox{Hz}\) increase with decreasing characteristic size of granulation, i.e., with higher levels of magnetic activity.

Thus, the mechanisms through which the frequencies are changed by the presence of the magnetic field associated with the solar cycle differ from those physical mechanisms that perturb the mode damping widths and amplitudes. Mode widths and amplitudes are determined by the excitation and damping of the acoustic oscillations in very shallow layers, where convection is most vigorous (Balmforth, 1992; Rimmele et al., 1995; Houdek and Dupret, 2015). Hence, this is where the changing magnetic field over the solar cycle affects mode widths and amplitudes most strongly. In contrast to this, mode frequencies can be perturbed by magnetic fields that are located far more deeply in the Sun (see, e.g., Gough and Thompson, 1990). An analysis of the timing of the changes in mode frequencies, amplitudes, and widths may yield information on the evolution and possibly the gradual ascent of magnetic field concentrations through the outer part of the convection zone.

The mode amplitudes as a function of mode frequency were fit with an asymmetric Voigt profile. From this, we found a frequency of maximum amplitude at \(\nu_{\text{max}}=3079.76\pm0.17~\upmu\mbox{Hz}\). The frequency of maximum amplitude \(\nu_{\text{max}}\) is of great interest to asteroseismic studies because it can be used in scaling relations for stellar mass and radius. The value found here compares well to those found in the literature (see the remarks in Section 3.3). Mode amplitudes are a function of radial order and frequency, as can be seen in the top panel of Figure 8. The spread in the distribution of the mode amplitudes as a function of mode frequency could not be removed by multiplication with mode inertia as for, e.g., mean squared velocity of the modes. This indicates that the excitation of modes is intrinsically not only a function of mode frequency. We find that amplitudes of modes of similar harmonic degree follow the asymmetric Voigt profile very well (bottom panel of Figure 8) and that the fit parameters vary with harmonic degrees of the modes, see the discussion in Appendix C.

The physical quantities of mean squared velocity, mode energy, and energy supply rate were calculated from the mode widths and amplitudes. Ridges of the same radial order are evident in the mode widths, mean squared velocities, and the energy supply rates (Figure 9). This was removed by scaling them with their mode inertia \(Q_{nl}\) (Figure 10), indicating that \(Q_{nl}\) is the correct scaling for mode widths, mean squared velocities, and energy supply rates (Komm, Howe, and Hill, 2000b). The maximum of the velocity of the modes is approximately \(37~\mbox{cm}\,\mbox{s}^{-1}\). This value is somewhat higher than the value of \({\approx}\,27~\mbox{cm}\,\mbox{s}^{-1}\) found by Chaplin et al. (1998) for radial modes from BiSON data. As modes of very low harmonic degree (below \(l=10\)) are not included in this investigation, it may well be that the velocity of these modes would be closer to the value of Chaplin et al. (1998) than that of modes for higher degree. This is supported by the fact that, on closer inspection, modes of lower harmonic degree are concentrated at the lower edge of the distribution shown in the top right panel of Figure 9.

To obtain the temporal variation of these quantities (i.e., mean squared velocity, energy, and energy supply rate; in physical units, not normalized to their temporal mean), an average over modes within a frequency range specific to each quantity (see Table 2) was calculated for the inertia-corrected parameter values.

The average of the mode widths was found to vary between a maximum width of \(1.36\pm0.01~\upmu\mbox{Hz}\) at the maximum of Solar Cycle 23 and a minimum of \(1.19\pm0.01~\upmu\mbox{Hz}\) during the activity minimum between Cycles 23 and 24. The mean squared velocity varied between extrema of \(1065\pm2~\mbox{cm}^{2}\,\mbox{s}^{-2}\) and \(1249\pm4~\mbox{cm}^{2}\,\mbox{s}^{-2}\). The mean mode energy exhibits variation between \(1.64\pm0.01\times10^{28}~\mbox{erg}\) and \(2.06\pm0.01\times10^{28}~\mbox{erg}\). This highlights that detectability of solar-like oscillation is reduced in stars with high levels of magnetic activity, as has been investigated and shown by Chaplin et al. (2011). These fractional changes are considerably different for modes of different frequencies, as can be seen from Figures 6 and 7. Qualitatively, our results agree with those of, e.g., Komm, Howe, and Hill (2000b), Jiménez-Reyes et al. (2003), and Salabert and Jiménez-Reyes (2006). The differences in the included mode degrees and frequencies and in the length of the data make it somewhat difficult to quantitatively compare the measured fractional parameter changes.

The variation in mean energy supply rate, shown in the lower right panel of Figure 11, is not correlated with the level of solar activity for the investigated set of modes. This confirms the results of, e.g., Chaplin et al. (2000), Jiménez-Reyes et al. (2003), and Broomhall, Pugh, and Nakariakov (2015), who found that the energy supply rate to solar p modes of low harmonic degrees does not change with the level of activity. Overall, there is a decrease in supply rates over the observed time period of about \(8\%\). However, this decrease happened entirely before GONG month 60, hence, before the upgrade of the GONG network. Since then, it has remained at a rather constant level. A dedicated study of energy supply rates, focusing on mode sets from different frequency ranges and harmonic degree ranges, as we did in Sections 3.1 and 3.2 for mode widths and amplitudes, will give further insight into whether the energy supply rates are truly constant over the solar cycle.