They Do Change After All: 25 Years of GONG Data Reveal Variation of p-Mode Energy Supply Rates

It has been shown over and over again that the parameters of solar p modes vary through the solar activity cycle: frequencies, amplitudes, lifetimes, energies. However, so far, the rates at which energy is supplied to the p modes have not been detected to be sensitive to the level of magnetic activity. We set out to re-inspect their temporal behaviour over the course of the last two Schwabe cycles. For this, we use Global Oscillation Network Group (GONG) p-mode parameter tables. We analyse the energy supply rates for modes of harmonic degrees $l=0\text{-}150$ and average over the azimuthal orders and, subsequently, over modes in different parameter ranges. This averaging greatly helps in reducing the noise in the data.We find that energy supply rates are anti-correlated with the level of solar activity, for which we use the $F_{10.7}$ index as a proxy. Modes of different mode frequency and harmonic degrees show varying strengths of anti-correlation with the $F_{10.7}$ index, reaching as low as $r=-0.82$ for low frequency modes with $l=101\text{-}150$. In this first dedicated study of solar p-mode energy supply rates in GONG data, we find that they do indeed vary through the solar cycle. Earlier investigations with data from other instruments were hindered by being limited to low harmonic degrees or by the data sets being too short. We provide tables of time-averaged energy supply rates for individual modes as well as for averages over disjunct frequency bins.


INTRODUCTION
Solar acoustic oscillations (p modes) are stochastically excited in the turbulent, outer convective layers of the Sun (e.g., Balmforth 1992;Rimmele et al. 1995;Houdek & Dupret 2015). They are well described by a damped and stochastically forced harmonic oscillator in the time domain (e.g., Anderson et al. 1990). In turn, in the Fourier domain a single p-mode peak is well characterised by a Lorentzian profile with a frequency , width Γ, and height (neglecting the observed mode asymmetry, e.g., Nigam et al. 1998;Korzennik 2017;Philidet et al. 2020). Measurements of these mode parameters in the Fourier domain can be used to learn about the energetics of the oscillator in the time domain: the peak width Γ is inversely proportional to the mode's lifetime and thus holds information about the damping. The product of peak height -or amplitude -and Γ, after proper scaling with the corresponding mode mass, encodes the energy that is stored in the respective p-mode oscillation (e.g., Goldreich et al. 1994;Baudin et al. 2005). Ultimately, the energy that is supplied to a mode per unit time / is proportional to the product Γ 2 · (e.g., Goldreich et al. 1994). This quantity is thus directly connected to the forcing of the mode (e.g., Chaplin et al. 2003). Possible temporal changes in the forcing function of the os-★ E-mail: r.kiefer@warwick.ac.uk cillator can thus be measured via time-resolved measurement of the quantity Γ 2 · .
Magnetic activity affects p modes in various ways: frequencies of modes below the acoustic cut-off (Balmforth & Gough 1990), i.e., modes that are essentially trapped in the solar interior, are correlated with the level of magnetic activity: frequencies are highest during times of strong solar activity (e.g., Woodard & Noyes 1985;Elsworth et al. 1990;Libbrecht & Woodard 1990;Jiménez-Reyes et al. 1998;Chaplin et al. 2001;Salabert et al. 2015;Tripathy et al. 2015;Broomhall 2017). For oscillations above the acoustic cut-off -these are waves that escape into the solar atmosphere (Balmforth & Gough 1990;Fossat et al. 1992;Vorontsov et al. 1998) -frequencies turn anti-correlated with the activity cycle (see, e.g., Woodard & Libbrecht 1991;Howe et al. 2008;Rhodes et al. 2010). The magnitude of activity-related shifts of mode frequencies depend on both mode frequency and harmonic degree (see Basu 2016 and references therein).
On the energetic side, p modes are increasingly damped with more activity on the surface, which leads to broader peak widths, as the two are inversely related to each other (e.g., Jefferies et al. 1991;Komm et al. 2000b;Chaplin et al. 2000;Jiménez et al. 2002;Jiménez-Reyes et al. 2003;Jiménez-Reyes et al. 2004;Salabert et al. 2007;Burtseva et al. 2009;Kiefer et al. 2018 and see Broomhall et al. 2014 for a more ex-tensive list of references on the matter). The observed changes along the solar cycle are measured to be between 10-20 % depending on mode frequency and harmonic degree (e.g., Kiefer et al. 2018). Mode amplitudes are anti-correlated with activity (e.g., Pallé et al. 1990;Elsworth et al. 1993;Chaplin et al. 2000;Komm et al. 2000b;Jiménez et al. 2002;Jiménez-Reyes et al. 2003;Jiménez-Reyes et al. 2004;Kiefer et al. 2018). Again, the magnitude of these changes is somewhat dependent on mode frequency and harmonic degree (e.g., Kiefer et al. 2018) but are typically between 10-25 % over the solar cycle. As mode amplitudes can be converted into mode energies, these are also anti-correlated with the level of activity.
The energy supply rate (or excitation rate, forcing rate, energy dissipation rate through the oscillator) of the p modes has thus far been assumed to be constant: This has been explained on the simple reason that changes in damping are of the same magnitude as the accompanying changes in mode energies. Thus, changes in mode damping are sufficient to explain the entire observed phenomenology (Chaplin et al. 2000): modes get broader and their amplitudes get smaller. In turn, the increased damping explains the lower mode energies during phases of higher activity without having to propose a variation to the forcing function. For more details on how mode damping, excitation, and forcing are related under the assumption of a harmonic oscillator, we refer to Chaplin et al. (2000), in particular their Section 2.
In addition, measurements of the energy supply rates -via the compound quantity Γ 2 · -have repeatedly returned constancy through the solar activity cycle (Chaplin et al. 2000;Jiménez-Reyes et al. 2004;Jiménez-Reyes et al. 2003;Salabert et al. 2007;. However, all of these studies investigated only low harmonic degrees ≤ 3 from diskintegrated Sun-as-a-star data. The energy supply rate of p modes is intimately connected with the spatial and temporal properties of convection (e.g., Samadi et al. 2003). The temporal correlation of convective eddies is crucial to the theoretical understanding and modelling of mode energy supply rates (e.g., Houdek 2010). Different analytical descriptions of the eddy-time-correlation can lead to vastly different values in the predicted energy supply rates (e.g., Belkacem et al. 2010). Given the observed changes of near-surface convection properties through the solar cycle, e.g., a decrease in the size of granules (Macris et al. 1984;Muller 1988) and decreasing granular contrast with increasing level of magnetic activity (Muller et al. 2007), it can be expected that the energy supply rates are in fact not entirely constant through the solar activity cycle.
Short term departures from constancy have already been detected for disk-integrated helioseismic observations before: Chaplin et al. (2003) report on a roughly 100 d long augmentation of the mode forcing during 1998 in Birmingham Solar Oscillation Network (BiSON) data.
A hint in GONG data towards activity-related changes of pmode energy supply rates was found by Komm et al. (2000a) but has never been substantiated. Also, Komm et al. (2000c) noted that energy supply rates decrease with increasing activity, but the detected change amounted to a mere −4.4 ± 4.3 % for the average over modes with 15 ≤ ≤ 95 and 1.6 mHz ≤ ≤ 3.1 mHz. Thus, no long-term variation of the mode energy supply rates that is (anti-)correlated with the level of magnetic activity has been conclusively confirmed as yet.
In the following we describe the data we used in this study (Section 2) and how we prepared and corrected these data (Section 3). The results from our analyses are presented in Section 4. A detailed  discussion follows in Section 5, closing with the conclusions we draw from this study and its findings in Section 6.

GONG mode parameters
At the time of writing, the publicly available GONG mode parameter files cover more than one full cycle of the solar magnetic field, i.e., more than two 11-year Schwabe cycles: the start date of the time series is May 7, 1995, the last available segment ends on April 12, 2020. We use the p-mode parameter tables, which were produced by the standard GONG pipeline from solar full-disk Dopplergrams (Anderson et al. 1990;Hill et al. 1996;Hill & Howe 1998) and are available online 1 . The parameters in these tables are obtained by fitting the power spectra of 108-day long GONG Doppler-velocity time series. Consecutive datasets overlap by 72 days, i.e., each third dataset is independent. Every 36 days a new dataset is added, which defines the unit of one GONG month. Thus, over one year about ten GONG months are accumulated.
In each GONG dataset's power spectrum, the oscillation multiplets are fitted with one symmetrical Lorentzian profile per azimuthal order − ≤ ≤ and a linear background for all harmonic degrees up to = 150. The free parameters of the fit model are mode frequency , mode width Γ (full width at half maximum), mode amplitude , background offset 0, , and background slope 1, , where the triple ( ) indicates the radial order , harmonic degree , and azimuthal order of a mode. In this article, we concentrate on the product Γ 2 · , which is, as we will see in further detail in Section 3.2, proportional to the mode energy supply rate / .
2.2 Solar radio flux 10.7 As a proxy for solar magnetic activity, we used the solar radio flux given by the 10.7 index (Tapping 2013), which is a good proxy for the level of activity in the upper chromosphere and the lower corona (Tapping 1987;. Its time series is available online 2 . For a measurement of 10.7 , the total emission on the solar disk at a wavelength of 10.7 cm is integrated over one hour. It is given in solar flux units sfu, where 1 sfu = 10 −22 W m −2 Hz −1 . The 10.7 index is the averaged 10.7 flux scaled for a distance of 1 AU. There is no continuous time series that covers the entire time period that GONG has been in operation. Thus, we concatenate the tables F107_1947_1996.txt, F107_1996_2007.txt, and fluxtable.txt (data since 2004), which are available at the given URL, at the end date of each preceding table, ignoring the overlapping data. Measurements which are taken on the same day are averaged. These daily values were then averaged over the 108day periods of the GONG datasets to yield 10.7 .
In Figure 1, the resulting 10.7 index is shown as a function of time, limited to the period covered by the available GONG data. The blue horizontal line indicates the median value of the 10.7 index through the used time frame, 10.7, median = 95.6 sfu. Data points above this value are coloured red, points below are black. We will use this colour code in later figures to ease the interpretation of the data.
The start and end dates of the periods of solar activity minima and maxima through the GONG observations are listed in Table 1. These periods are shaded grey in Figure 1 and the respective identifier is printed at the top of the panel. Again, we will reuse this background colouring in later figures to ease identification of activity minima and maxima. The dates for the first three extrema are taken from Broomhall (2017). The dates of Max 24 were chosen such as to have a similar number of GONG datasets as the preceding maximum at times of high activity. The time when 10.7 has reached a comparable level of low activity as Min 24 was chosen as the starting date of the currently ongoing minimum. The number of GONG datasets during each of these periods are given in the last column of Table 1.

DATA PREPARATION
The data processing steps and the corrections that were applied to the quantity Γ 2 · are the same as described in detail by Komm et al. (2000b,c), and Kiefer et al. (2018) for its constituent quantities Γ and . Here, we will briefly lay out the rationale of these corrections and refer to Kiefer et al. (2018) for more details.

Correction for spatial masking and azimuthal averaging
Due to projection effects that increase towards the solar limb -affecting pixel resolution as well as measured radial velocities -mode amplitudes of azimuthal orders with lower | / | are suppressed. This suppression follows an empirically determined polynomial in ( / ) with = 0, 2, 4. To correct for this, we subtracted the fit from each multiplet and added the value of | / | = 1.
This fit, as all linear regressions in this article, was performed with 's weighted linear regression routine (Seabold & Perktold 2010). In this particular case, the fit was done multilinearly in ( / ) 2 and ( / ) 4 . The variance-weighted mean of the corrected data points is then adopted as the representative value of Γ 2 · of each multiplet. For a multiplet to be considered at all, a minimum of one third of azimuthal orders was required to be present for multiplets with > 10 and two thirds for multiplets with ≤ 10.
In the following, the index is dropped from all quantities, as the azimuthal orders of a multiplet have been averaged over. The frequency of the = 0 singlet is taken as the multiplet frequency. If the = 0 was not fitted, then the mean of the = ±1 singlet frequencies is used. As the azimuthal orders have been averaged over, from here on out the phrase multiplet will be avoided where possible and the term mode will be used instead.

Correction for temporal window
GONG is a network of six observing stations that are distributed around the Earth to maximise continuous viewing of the Sun. Still, the duty cycle (henceforth called fill) of the time series is lower than unity for all GONG datasets. 3 The fill of the 108-day time series varies between a minimum of 69.7 % and a maximum of 93.7 % with a median value of 86.7 %. This temporal window function artificially increases mode widths and decreases mode amplitudes.
As the fill is rather high and the gaps are typically not ordered in a repeating pattern, it is sufficient to account for its impact by a linear regression (Komm et al. 2000b): Γ 2 · of each mode are fitted with a linear function as a function of fill, which is then subtracted and the extrapolated value at fill = 1 is added.

Apparent solar radius
Just as Kiefer et al. (2018), we corrected for residual annual changes of the solar angular radius. For this, we performed a linear fit of Γ 2 · as a function of the solar apparent radius and adjusted Γ 2 · to its extrapolated value at the maximum of the apparent solar radius during the time series. Kiefer et al. (2018) discovered two jumps in the mode amplitudes, which had to be corrected for with an empirical correction factor. The first jump occurs around the year 2001 and is due to an upgrade of the GONG hardware. The second jump, around 4 years later, is of uncertain origin. We continue using the correction factor of Kiefer et al. (2018), which they give in their appendix A, equation (12).

Energy supply rates
Subsequently, we will consider either the energy supply rate parameter of individual modes, which, for brevity, we baptise Θ = Γ 2 · , or the physical quantity of energy supply rate / , see Eq.
(1). If the average over modes in certain ranges of mode frequency and harmonic degrees is taken, we drop the indices ( ). Figure 2 shows the normalised variation of Θ as a function of time and mode frequency (left panel) and harmonic degree (right panel).

Maps of the temporal variation of Θ
Here, the modes' Θ are first normalised to the mean over the entire time series. In the left panel, modes are grouped in independent bins of 50 µHz and averaged. In the right panel, the average is taken for all modes of the same harmonic degree. The colour map is capped at [0.9, 1.1]. It can be seen that there is considerable temporal variation in Θ, with bluer, i.e., larger values during activity minima as stated in Table 1. Red values, which indicate smaller Θ, can predominantly be found at times of high magnetic activity.
The rather sharp blue vertical feature at around the year 2001 is due to the GONG hardware upgrade at the time and can be considered an instrument artefact.
From these two panels, it can already be estimated that any cyclic variation in the mode energy supply rates is likely less apparent in modes with frequencies 3500 µHz and harmonic degrees 120. Also, even though the variation with apparent solar radius has been accounted for, a seemingly repeating, quasi-yearly pattern can be identified. We will discuss this further in Section 5.

Conversion to energy per unit time
From the measured quantity Θ = Γ 2 · , the energy that is supplied to the p-mode oscillations per mode can be calculated by (see, e.g., Goldreich et al. 1994) where Γ is in units of Hz, has units cm 2 s −2 Hz −1 . The numerical factor vis = 3.33 corrects for the GONG specific reduced visibility of modes due to leakage effects (Hill & Howe 1998). We will comment on the use of this factor in Section 5.3. The factor 2 stems from the Lorentzian profile that is used to fit the mode peaks and the assumption of scattering being the main contribution to mode damping. The mass of the mode ( ) is calculated as with the solar mass ⊙ in g and the mode inertia of solar Model S (Christensen- Dalsgaard et al. 1996), which is defined as where is density, and ℎ are the mode's radial and horizontal displacement eigenfunctions, ⊙ is the solar radius, is the radial coordinate, and ′ is the radial position at which the normalisation is done. We normalise the mode inertiae at either the photospheric radius ′ = ⊙ or the GONG observation height at ′ = ⊙ + 240 km ≡ obs (Baudin et al. 2005). Unless stated otherwise, mode inertiae are normalised at the GONG observation height obs . Values normalised at ⊙ are provided in an online-supplement table.
We use cgs units, i.e., / has units erg s −1 . More detailed discussions of how Eq. (1) comes about can be found in, e.g., Goldreich et al. (1994) and Komm et al. (2000c). We note that normalised variations of / as a function of time are equivalent to normalised variations of Θ in time, as the factors 2 vis then cancel. Figure 3 shows an --diagram of the 1580 modes which are present in at least 50 % of the GONG datasets. We excluded modes with = 0, 1 as they appear as outliers in the distribution of supply rates (see Figure 4). There are 26 radial and dipole modes, i.e., the total number of modes with a presence rate of at least 50 % is 1606. The colour of the points indicates log 10 / , where the variance-weighted average over all time samples (indicated by the overline) is taken as the representative value for each mode. The vertical and horizontal dashed grey lines indicate the boundaries in mode frequencies and harmonic degrees that separate mode sets which we use in later sections.
The top left panel of Figure 4 shows the energy supply rates of all modes as a function of mode frequency. Again, the varianceweighted average over all time samples is taken. Squares highlight radial modes, triangles are used to accent dipole modes. Here it can clearly be seen that radial and dipole modes are outliers from the general bulk of data. This problem arises from the non-inclusion of the instrument-specific leakage matrix in the GONG fitting pipeline. We do not adopt individual correction factors for different harmonic degrees in this article (consistently with Komm et al. 2000b). From here on out we shall not consider radial and dipole modes any more.
In the bulk of the data points, ridges can be identified. These are made up of modes of like radial order , as is highlighted by the discrete colour map. These ridges can be seen more clearly when / is plotted as a function of harmonic degree , see the bottom left panel of Figure 4.
In Table A1 we give the energy supply rates for averages over 32 modes consecutive in mode frequency. The second and third columns of this table are energy supply rates without and with rescaling by averaged over the entire time series. From the second column we find that energy supply rates peak at around 3460 µHz with a value of 13.394 ± 0.004 × 10 22 erg s −1 . The entries in the following columns are averaged over periods of activity extrema, i.e., all activity minima or maxima as defined in Table 1, and scaled by in the last two columns. This table is also available in machinereadable form online.

Accounting for different mode inertiae
As was shown by, e.g., Komm et al. (2002), when plotted as a function of mode frequency the ridges of individual radial orders are collapsed almost completely into one ridge by multiplication of / with the scaling factor Here, the radial mode inertia 0 is interpolated to the mode frequency to obtain 0 ( ). The right column of panels in Figure 4 shows the time-averaged energy supply rate values as in the left column, but scaled by multiplication of . Now, as can be seen in the top right panel of Figure 4, the supply rates can well be described by a function which depends only on mode frequency . We used a non-linear least squares optimization to fit the exponential to the base-10 logarithm of the resulting energy supply rates. This fit is shown by the solid red line. In Eq. (5), the parameter tuple comprises: offset , magnitude , central frequency 0 , and width of the exponential . The best-fit parameters and their standard uncertainties are given in Table 2. Though present in the panel, the radial and dipole modes are excluded from the fit. We also tested an asymmetric and a symmetric Voigt profile, both of which had a comparable goodness-of-fit as the exponential of Eq. (5). Larger deviations from the exponential behaviour occur at low mode frequencies 1800 µHz and at high mode frequencies 3600 µHz. At low frequencies, this is at least partially caused by the limited frequency resolution of the 108-day time series which is ≈ 0.107 µHz. Thus, only a few frequency bins are available per mode width which become smaller with decreasing mode frequency, see, e.g., Komm et al. (2000c) and Kiefer et al. (2018). This can lead to an imprecise estimation of mode widths and thus also of / . At high frequencies 3600 µHz some radial orders depart from the monotonic increase with mode frequency. This can be appreciated in the bottom panels of Fig. 4, where the data points with 20 show a dip in / . In a machine-readable online-only table, we provide the energy supply rates of all 1606 modes which are present in at least 50 % of the GONG samples. There, we provide mode inertiae and scaling factors normalised at both the photospheric radius ⊙ and the GONG observation height in the solar atmosphere obs . The energy supply rates are then provided for 12 different averages: averaged over all time samples with and without scaling by , averaged over times of the activity minima and maxima, both with and without -scaling, and all of these at both normalisation heights ⊙ and obs .

Averages over parameter ranges
As can already be appreciated from Figure 2, there is a temporal variation in the energy supply rates. In order to decrease the uncertainties, we take variance-weighted averages over four mode frequency ranges and five ranges in harmonic degree. These ranges are given in the first four columns of Table 3 and are indicated by grey dashed lines in Figure 3. As we set the presence rate of modes to 50 %, not all modes are present for all time samples. The range of number of modes within each parameter range that is averaged over is given in the fifth column of Table 3.
To select the modes which contribute to the variance-weighted average within a parameter range we do the following: First, we select all the modes, which are within the minmax and minmax ranges. The energy supply rates of each mode are then normalised to the variance-weighted mean over the entire time series. Modes which do not satisfy the 50 % quota are eliminated from the set. At each time sample we then remove all modes whose normalised (1) as function of mode frequency (top panel) and harmonic degree (bottom panel). Right column: As left column, but supply rates are inertia-scaled by multiplication with Eq. (4). Top right panel: The solid red line is a fit of the exponential function Eq. (5) to the supply rates. All panels: Radial orders are indicated by the colours of the data points. Radial and dipole modes are highlighted by red-bordered triangles and squares, respectively. variation is not within the interval [0.5, 1.5], as very large deviations from 1 are not expected to occur. Further, we remove all modes that deviate by more than four times the standard deviation of the contributing modes at each time sample. These two outlier rejections only remove a handful of modes per time sample, if any.
The variance-weighted averages of the normalised variation of / over these mode sets are presented in Figure 5. Parameter ranges of the mode set entering each panel are indicated to the top of each column and the right of each row. The error bars indicate 5 , i.e., uncertainties on the individual data points are very small due to the averaging over azimuthal orders and the modes in each range.
The minimum and maximum variations of each panel are given in the sixth and seventh column of Table 3. To reduce short-term variation caused by seasonal variations in data quality that are not captured by the fill factor correction, these values are calculated from the one-year boxcar smoothed data. The median uncertainty of the unsmoothed data is given in the last column of Table 3. The temporal variation of / is statistically significant in all 20 parameter ranges.

Low harmonic degrees
To emulate the observations of disk-integrated helioseismic observations, we also averaged the lowest available harmonic degrees = 2-4 over the complete frequency range. Figure 6 shows the resulting normalised variation as a function of time. Due to the smaller number of azimuthal orders per multiplet and the smaller number of multiplets overall (between 28 and 43, depending on the time sample), the averages are less well constrained than in the wider harmonic degree ranges shown in Figure 5. The solid red line shows the running one-year variance-weighted mean with its uncertainty in the red shaded band. As in Fig. 5 the error bars indicate 5 uncertainties.
The extrema of these smoothed data are 92.2 ± 0.3 % and 106.4 ± 0.3 %, where the median uncertainty on the un-smoothed data point is 0.9 %. We will discuss this further in Section 5.

Cross-correlating magnetic activity and
To assess the correlation between the variation observed in / and the level of solar magnetic activity, as measured by the 10.7 index, we computed the cross-correlation function (CCF) for the independent time samples, i.e., every third time step. We again split the modes into the subsets we used before: the modes included in the calculation of the functions shown in the panels of Figure 7 correspond to those of Figure 5. In Figure 7, the cross-correlation functions are shown by solid blue lines with blue dots. The first of the three possible sets of independent data points was chosen for analysis. Fourth-order polynomial fits to the CCFs are shown by the solid red curves. The dotted vertical blue line indicates the global minimum of each panel's CCF, while the dashed vertical red line indicates the minimum of the fitted function. These fits are performed with 's least-square polynomial fit routine numpy.polyfit (Oliphant 2006). The blue shaded bands indicate the minimum-to-maximum range of the three possible independent data sets. The correlation values typically vary only by a few percent between the three data sets.
In Table 4 the Pearson correlation coefficient and its -value for the different mode sets are given for different lag values: for Figure 5. Temporal variation of energy supply rates 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). Mode energy supply rates are normalised to the mean for each mode and then the variance-weighted average over all modes in the respective range of frequency and degree is taken. Data with higher than median 10.7 solar radio flux are highlighted by red points. Levels of 1.1 and 0.9 of the mean are indicated by dashed lines. Grey shaded times corresponds to the activity minima and maxima laid out in Figure 1 and Table 1. un-lagged data in columns 5 and 6; for the lag of the minimum of the cross-correlation function, in columns 7-9; for the lag of the minimum of the quadratic fit to the cross-correlation function, in columns 10-12. The last two columns give Spearman's rank correlation for un-lagged data and its corresponding -value.
For the un-lagged data the -values are < 0.05 for 15 out of the 20 mode sets. For low harmonic degrees = 2-30 only the low frequency modes have < 0.05. In addition, for the high frequency modes, the ranges = 61-100 and = 101-150 have > 0.05. This essentially also holds true for the other lag values (global minimum of the CCF and minimum of the fitted function) and the Spearman : For the low-, and mid-frequency sets, as well as for the full mode set, the -value of the correlation coefficients is < 0.05 for all harmonic degree ranges. Indeed, overall the correlation coefficients indicate a moderate to strong anti-correlation between the level of solar magnetic activity and the p-mode energy supply rates.
The strongest anti-correlation with = −0.82 is found for the low frequency mode set with = 101-150 at zero lag. Except for the high harmonic degree sets, the minima of all (statistically significant) CCFs are found at positive lag values. Table 3. Normalised and averaged variation of energy supply rates for modes sets with different parameter ranges (defined in columns 1-4) that are presented in the panels of Fig. 5. Column 5 gives the range of number of modes included in the GONG datasets. Columns 6 and 7 give the extrema of the one-year smoothed variation and their uncertainty. The last column gives the median error of the individual un-smoothed energy supply rate in this parameter range.  Table 4. Correlation coefficients between mode energy supply rates of different parameter ranges (defined in columns 1-4) and the 10.7 solar radio flux for different lags: Pearson and corresponding -value for un-lagged data in columns 5 and 6. Lag of the minimum of the cross-correlation function, this lag's Pearson and -value in columns 7-9. The lag of the minimum of the polynomial fit to the cross-correlation function, this lag's Pearson and -value in columns 10-12. Spearman's rank correlation for un-lagged data and corresponding -value in columns 13 and 14.

Comparison of activity extrema
We compare the energy supply rates during the different activity extrema listed in Table 1 with each other. For this, we average the energy supply rates of each mode over the time samples during the periods of Table 1, take the difference between two extremal periods, and normalise by the mean over the entire time series.
As activity and energy supply rates are anti-correlated, we subtract successive periods of high activity from periods of low activity for better comparability of the results. Also, to assess the relative depth or strength of minima and maxima, we compare the minima of cycles 25 and 24 as well as the maxima of cycles 24 and 23. Only energy supply rates of modes which are present in two periods can be subtracted from each other. Thus, if two rows of Table 5 are subtracted from each other, e.g. row three from row two to seemingly yield Max 24 -Max 23 , there is a small difference to the values given in the last row. In this example, the last row only necessitates the modes to be present in Max 24 and Max 23 , whereas taking the difference between rows three and two demands the modes to be present in all of these rows' four extrema.
In Figure 8, the normalised differences between the different extrema are shown for all modes as a function of mode frequency (black data points). The extrema that are compared in the respective panel are indicated to the left of each panel. In order to smooth out the scatter in the normalised differences, we calculated the rolling variance-weighted average over 100 modes consecutive in mode frequency. In each panel, this is shown by the solid coloured curve and its 10 confidence interval is given by the coloured band. The 1 uncertainty is calculated as the standard error of the weighted mean. Due to the heavy averaging that has been applied up to this point, the uncertainties are rather small. Figure 9 shows the same as Figure 8 but as a function of harmonic degree. Here, the weighted averaging is first done over the modes of each individual degree and then the rolling variance-weighted average over five consecutive degrees is calculated.
For better comparability, Figure 10 collects the smoothed curves from Figures 8 and 9. Here, the left panel shows the normalised differences as a function of mode frequency and the right panel as a function of harmonic degree.
In Tables 5 and 6, the percental energy supply rate change of modes in the four mode frequency ranges and five harmonic degree ranges are listed. We will discuss these results further in the next section.

Dependence on mode frequency and harmonic degree
The energy supply rates of solar p modes vary over the solar activity cycle. The magnitude of this variation and the correlation with the 10.7 index varies with mode frequency and harmonic degree.
From Fig. 5 and Table 3 we find that energy supply rates of modes in the range = 61-100 with frequencies < 2400 µHz vary most strongly with a minimum of 87.5 ± 0.1 % during times of high activity and a maximum of 105.9 ± 0.1 % at low activity levels. The anti-correlation of this mode set with activity is strong with (0) = −0.75 and < 10 −3 . The strongest correlation is found for the mid-frequency mode set with = 101-150, for which (0) = −0.82 and < 10 −3 . The amplitude of the variation over the activity cycle is smaller for higher frequency modes, which can be appreciated from Fig. 5, Table 3, and particularly from Fig. 10. The anti-correlation at zero lag between these high frequency modes and 10.7 has < 0.05 only for the range = 31-60 and the full range of harmonic degrees.
Interestingly, the strongest -and most significant -anticorrelation values for all mode sets except the highest harmonic degrees = 101-150 are found at positive lag values. I.e., it takes the energy supply rates a certain amount of time to react to changes in the level of magnetic activity.

Why has this not been detected before?
As the detected variation of / over the activity cycle is rather large and very significant for the majority of mode sets, it is worthwhile to ask why this variation has not been detected until now.
Previous studies of the solar-cycle variation of p-mode parameters often focused on mode frequencies (see Section 1 for references), as these hold the most accessible information about the deep layers of the Sun. When studies did include the other mode parameters -damping width Γ, mode amplitude , and quantities compound out of these two -they were almost exclusively limited to harmonic degrees ≤ 2. The total number of azimuthal orders per radial order of the multiplets = 0, 1, 2 is 9. What is more, this effectively reduces to 6, as only modes for which + is even are clearly seen in Sun-as-a-star data. Thus, even in the case if all harmonic azimuthal orders are fitted equally precisely (which is not the case for unresolved observations e.g., Chaplin et al. 2004), this is equivalent to using only the = 3 modes from our study and averaging over their azimuthal orders. Mind that Figure 6 in fact includes = 2-4 from resolved observations, i.e., 21 azimuthal orders per radial order.
In addition, it can be seen from the right panel of Figure 10 that lower harmonic degrees have smaller normalised differences between activity extrema. Thus, these modes' / appear to be less sensitive to activity compared to modes of higher degree. The larger deviations from unity around the end of Max 23 and just before the beginning of Min 24 seen in Figure 6 can also be found in the larger mode sets in Figure 5: low-degree modes do appear to have at least some response to the changing level of activity. Overall, the Pearson correlation-coefficient at zero lag between the normalised differences of the = 2-4 mode set plotted in Figure 6 and the 10.7 index is only (0) = −0.18 with a -value of 0.1.
From this we conclude that some combination of a smaller number of multiplets in total, the small number of azimuthal orders to average over, the shorter length of the observations, and the fact that lower harmonic degrees appear to be less sensitive to magnetic activity, lead to the null results of earlier studies on the matter of activity-related variations of / . Figure 11 shows the normalised variation of / for the full mode set (average over all frequencies and harmonic degrees) as a function of 10.7 . From a linear fit (blue line) to this we find that

Sensitivity of / to magnetic activity
[ / ] = (1.050 ± 0.007) − (4.8 ± 0.6) · 10 −4 sfu −1 10.7 , where the square brackets indicate that the energy supply rates are normalised to the mean over the full cycle and the entire mode set is included. Again, the linear fit shows the anti-correlation of the energy supply rates with activity: The quiet Sun supply rates are > 1 and the slope is < 0. Many of the data points above the linear fit belong to the rising activity phase going from minimum 23 to maximum 23. Whether this is part of a hysteresis pattern in / over a full magnetic cycle or an artefact from the GONG data and the GONG hardware upgrade just at the end of the maximum of cycle 23 should be tested with data from other instruments. The much larger percental change in / between the minimum and maximum of cycle 23 compared to the later extrema can also be seen in Tables 5 and 6. When comparing Min 24 to Min 23 directly, we find that the supply rates of the full mode set are 5.51 ± 0.02 % lower for Min 24 , which is obviously in contrast to the surface tracers of activity, e.g., the sunspot cycle. It is also conceivable that, due to the different strength of cycle 24, the energy supply rates have settled at a slightly different level when compared to the earlier stronger cycles which had larger amplitudes between their activity minima and maxima.
As Figure 8 and Table 5 show, modes around and above the Sun's frequency of maximum oscillation amplitude max (≈ 3080 µHz, e.g., Kiefer et al. 2018) are less sensitive to magnetic activity than modes below max . Interestingly, this is in contrast to mode frequencies for which the amplitude of change along the cycle increases with mode frequency (e.g., Broomhall 2017). Furthermore, let us consider the results of Kiefer et al. (2018), in particular their figures 6 and 7 which show the analogons of our Figure 5 but for mode widths Γ and mode amplitudes . While Kiefer et al. (2018) did not compare activity extrema separately, it can be appreciated that the percental variation in the low frequency mode set going from one activity extremum to the next is larger for / than it is for Γ or . Also, the variation is larger for the low frequency mode set than the mid frequency set for / , but for Γ and separately, the variation is larger in the mid frequency set. This is most apparent for the damping widths Γ in figure 6 of Kiefer et al. (2018). All this shows that it is indeed not an increase in the damping alone that is causing the variations in Γ and as has been assumed thus far. Indeed these observations call for an activity related variation in the forcing function of solar p modes.
Lag between / and 10.7 As we have established, the high frequency mode set is least reliable for the purpose of investigating / ; so let us exclude this set from the considerations of this subsection and focus on the other three frequency ranges.
Looking at the minima of the CCFs between / and 10.7 and the minima of the fits to them, we find that the lag value for these decreases with increasing harmonic degree. For example, in the mid frequency range, the position of the global minimum of the CCF decreases from 864 d for modes in the range = 2-30, to 756 d for = 31-60, to 108 d for = 61-100, to finally 0 d for modes with = 101-150. This is interesting as the higher harmonic degree modes are confined to more shallow layers than modes of lower harmonic degrees (Basu 2016). It thus appears that surface activity affects those modes earlier which are confined closer to the surface and modes that penetrate deeper into the Sun are fully affected later. As p modes are assumed to be excited very close to the  Table 1 as function of mode frequency. The variance-weighted average of each mode's energy supply rate is calculated over the extent of each extremal period. Differences are normalised by the mean over the entire GONG time series. The variance-weighted moving mean over 100 modes consecutive in frequency and its 10 uncertainty are shown by the coloured lines and the shaded bands. Table 5. Differences in energy supply rates between the activity extrema stated in the first column for modes with frequencies stated in the top row. All values are given in %. The dates of the activity extrema are taken from Table 1. Values are normalised by the mean over the complete time series, as in Figure 8. Then, the variance-weighted average over all modes in the parameter range is taken. surface (e.g., Chaplin & Appourchaux 1999) it is not immediately apparent why low degree, high inertia modes should take longer to react to changes in surface activity. We speculate that the higher mode inertia of lower frequency and lower harmonic degree modes acts as a sort of buffer against changes in the mode forcing, i.e., a stronger or longer-lasting perturbation is needed to induce a change in their energy supply rates.
We further note the curious behaviour in all mode sets' CCF just after the zero lag: compared to the polynomial fit, the CCF dips to smaller correlation values at around a lag of 108 d. It appears to oscillate back above the fit only to decrease again slightly. We speculate that this (quasi-)oscillatory behaviour in the CCF might be connected to the Sun's quasi-biennial variations that have been observed in, e.g., p-mode frequencies (Broomhall et al. 2012). This dip persists even if all time samples are included in the calculation  of the cross-correlation functions. This feature is thus likely not caused by the time resolution of 108 d used in Fig. 7.

Quasi-yearly variation of Θ
Even after removing the residual change with apparent solar radius, the energy supply parameter Θ in Figure 2 shows a quasi-yearly variation that needs to be further investigated in future studies. We attempted to remove this variation by fitting Θ as a function of the solar position angle and the 0 angle. However, whatever small variation with these quantities were removed by the linear fit, the observed variation in Θ persisted. We thus refrained from considering this matter further and call for a more detailed investigation of this (quasi-)yearly variation. We have not investigated periodograms of the variation of / . It is therefore entirely possible that the short term variation that can so clearly be seen in Figure 2 has a period that is quite different from 1 yr and might be connected to the above mentioned quasi-biennal oscillations. Part of the residual seasonal variation may be connected to changes in the GONG leakage matrix (Hill & Howe 1998).

Impact on asteroseismology
The accurate prediction of the amplitudes of stellar p-mode oscillations is crucial in the target selection of photometric space missions. Such a target selection has been done for the asteroseismic targets of the NASA Transiting Exoplanet Survey Satellite mission (TESS, Ricker et al. 2014) based on the results of Campante et al. (2016) who calculated a detection probability for solar-like oscillations. Their equation for the prediction of the oscillation amplitudes in the power spectrum contains factors that account for the length of the data strings, sampling and instrumental effects, the pattern of p-mode peaks in the power spectrum and a factor that gives the maximum oscillation amplitude of radial modes. The calculation of this amplitude factor is based on scaling relations based on the work of Kjeldsen & Bedding (1995) and Chaplin et al. (2011). Campante et al. (2016) mention the need for accounting for the effect of magnetic activity in the estimation of mode amplitudes. During the Kepler mission (Koch et al. 2010;Borucki et al. 2010), over 2000 stars were observed in the survey phase. Of these, only 540 were detected to exhibit solar-like oscillations (Mathur et al. 2019). It is known that -as in the solar case -stellar oscillations are suppressed by magnetic activity (e.g., García et al. 2010;Salabert et al. 2017;Kiefer et al. 2017;Santos et al. 2018). Indeed, Mathur et al. (2019) were able to attribute 32 % of nondetections to magnetic activity. The remaining 68 % however are still to be accounted for.
Based on the results we present in this article, we conjecture that it is not only increased damping of oscillations and the ensuing attrition of mode amplitudes that impede the detection of stellar p-modes and need to be accounted for, but in fact the oscillations are fed with less energy per unit time even for moderately active stars like the Sun.
In the target selection of future asteroseismic surveys this should be taken into account for stars where a reliable measure of activity is available. If the detection probability of oscillations is to be maximised, more active stars ought to be deprioritised. More work on this has to be done, e.g., on the extent to which the lowest harmonic degrees are affected by activity. This problem could be tackled by careful analyses of the available photometric time series on the one hand (Salabert et al. 2017) and simulations of mode excitation in solar and stellar convection zones (e.g., Belkacem et al. 2006;Samadi et al. 2007;Zhou et al. 2019;Zhou et al. 2020) that investigate the dependence of the mode excitation rate in these simulations on the magnetic field strength on the other hand.

Shortcomings
We have corrected the jumps in mode amplitudes that were noticed by Kiefer et al. (2018) in the same empirical way they used: we applied a correction factor to remove the two obvious jumps (see Section 3.1). This is unsatisfactory, as neither the origin of the second jump is explained, nor can we be certain that the direction of the adjustment is correct. It is entirely possible that the absolute values of the mode energy supply rates as presented in Figures 3, 4 and Table A1  We have neglected possible contributions from time-dependent mode asymmetry in this study. We are aware that some fraction of the detected activity-related change of supply rates can potentially be attributed to changes in mode asymmetry. With the publicly available official GONG data such analyses cannot be done, as the profiles that are used in the GONG mode fitting pipeline are symmetrical. It is therefore expedient to repeat similar analyses as we have done here with data that include a mode asymmetry parameter (e.g., Korzennik 2017) and, in particular, with data from the SOHO-MDI (Domingo et al. 1995;Scherrer et al. 1995) and SDO-MDI (Pesnell et al. 2012;Schou et al. 2012) instruments.
The GONG leakage matrix is not included in the fitting of the mode peaks. We account for part of this shortcoming by removing the bowl-shape of the energy supply rates as a function of | / |, see Section 3.1. This effectively takes care of the -dependence of the self-leakage. However, the overall level of self-leakage is only naively corrected for by the inclusion on the factor vis in Eq. 1. We stress that using this factor for all harmonic degrees is consistent with Komm et al. (2000b). In reality however, it can vary for different harmonic degrees (see, e.g., Baudin et al. 2005). We do not anticipate vis to vary over the course of the Sun's activity cycle, thus it cannot contribute to the detected cycle variation -it only affects the absolute values of energy supply rates. In a numerical estimation of the geometrical visibilities of different harmonic degrees in resolved Doppler observations, we found that the selfleakage levels off at relatively low harmonic degree ≈ 4. Thus, the reported absolute averages of energy supply rates are most likely only off by a small constant factor, as the majority of modes have > 4.

SUMMARY AND CONCLUSIONS
We have, for the first time, analysed the solar p-mode energy supply rates up to = 150 through two solar activity cycles. Contrary to the long-held opinion that the energy that is fed into the p-modes per unit time is constant, we have found that they are indeed sensitive to the level of magnetic activity. Averaged over the full mode set, the minimum supply rates are found during activity maxima, with a minimum of 91.2 ± 0.1 % normalised to the entire time series. At activity minimum they reach as high as 107.6 ± 0.1 %, i.e., supply rates are in anti-phase with the activity cycle with a correlation of = −0.50 ( < 10 −3 ), where the 10.7 index was used as a proxy for solar activity.
The fractional change through the activity cycle depends somewhat on mode frequency and only weakly on harmonic degree: supply rates of lower frequency modes < 3000 µHz are more strongly affected than higher frequency modes. It appears that the supply rates of lower harmonic degree modes 20 tend to be least affected by activity.
According to the mode energy supply rates, the currently ongoing minimum is as deep as the long minimum between cycles 23 and 24. Averaged over the full mode set, there is very little variation between these two minima with a difference of only 0.33 ± 0.02 %. Also, it is reflected in the energy supply rates that the maximum of cycle 23 was stronger than that of cycle 24: the energy supply rates were 2.43 ± 0.02 % higher during the latter.
As the energy supply rates are a crucial ingredient in simulations and prediction of the solar and stellar p-mode amplitudes, it is worthwhile to establish their ground-state, i.e., quiet Sun levels. We provide a table with the energy supply rates as a function of mode frequency for the average over the entire time series, corrected for mode inertia, averaged over the periods of activity minima as well as maxima, and with the mode inertiae normalised at both the photospheric radius ⊙ as well as the GONG observation height obs . Table A1 as well as a larger table with the energy supply rates of 1606 modes are provided online.
Our findings could -and should -be corroborated by analysing the last two activity cycles with data from multiple instruments. Optimally, such a study is done with results obtained with both symmetrical and asymmetrical mode profiles. This way, it could be determined to what extent the here reported change in supply rates over the cycle is in fact due to variations in mode asymmetry. Finally, we note that the convective forcing of acoustic mode in stars with different metallicity and surface gravity than the Sun might respond more weakly or strongly to magnetic activity. Table A1. Energy supply rates averaged over 32 modes consecutive in mode frequency; the last row only contains 12 modes. The first column gives the mean mode frequency of each bin. The second and third columns are energy supply rates without and with rescaling by averaged over the entire time series with variance-weighting over the 32 modes. The values in the following columns are averaged over only the indicated periods of activity extrema, i.e., all activity minima or maxima as defined in Table 1, and are scaled by as indicated. To maintain readability, uncertainties are rounded up to 0.001. More digits are given in the online version of this