11institutetext: Department of Computer Science, Aalto University, PO Box
15400, FI-00076 Aalto, Finland 22institutetext: Max Planck Institute for Solar
System Research, Justus-von-Liebig-Weg 3, D-37077 Göttingen, Germany
33institutetext: Nordita, KTH Royal Institute of Technology and Stockholm
University, Hannes Alfvéns väg 12, SE-10691 Stockholm, Sweden

# Solar-Cycle Variation of quiet-Sun Magnetism and Surface Gravity Oscillation
Mode

M. J. Korpi-Lagg 1133 A. Korpi-Lagg 1122 N. Olspert 11 H.-L. Truong 11

###### Abstract

Context. The origin of the quiet Sun magnetism is under debate. Investigating
the solar cycle variation observationally in more detail can give us clues
about how to resolve the controversies.

Aims. We investigate the solar cycle variation of the most magnetically quiet
regions and their surface gravity oscillation ($f$-) mode integrated energy
$E_{f}$.

Methods. We use 12 years of HMI data and apply a stringent selection criteria,
based on spatial and temporal quietness, to avoid any influence of active
regions (ARs). We develop an automated high-throughput pipeline to go through
all available magnetogram data and to compute $E_{f}$ for the selected quiet
regions.

Results. We observe a clear solar cycle dependence of the magnetic field
strength in the most quiet regions containing several supergranular cells. For
patch sizes smaller than a supergranular cell, no significant cycle dependence
is detected. The $E_{f}$ at the supergranular scale is not constant over time.
During the late ascending phase of Cycle 24 (SC24, 2011-2012), it is roughly
constant, but starts diminishing in 2013, as the maximum of SC24 is
approached. This trend continues until mid-2017, when hints of strengthening
at higher southern latitudes are seen. Slow strengthening continues, stronger
at higher latitudes than at the equatorial regions, but $E_{f}$ never returns
back to the values seen in 2011-2012. Also, the strengthening trend continues
past the solar minimum, to the years when SC25 is already clearly ascending.
Hence $E_{f}$ behavior is not in phase with the solar cycle.

Conclusions. The solar cycle dependence at the supergranular scale is
indicative for the fluctuating magnetic field being replenished by tangling
from the large-scale magnetic field, and not solely due to the action of a
fluctuation dynamo process in the surface regions. The absence of variation at
smaller scales might be an effect of the limited spatial resolution and
magnetic sensitivity of HMI. The anticorrelation of $E_{f}$ with the solar
cycle in gross terms is expected, but the phase shift of several years
indicates a connection to the large-scale poloidal magnetic field component
rather than the toroidal one. Calibrating AR signals with the QS $E_{f}$ does
not reveal significant enhancement of the $f$-mode prior to AR emergence.

###### Key Words.:

the Sun: helioseismology; magnetic fields; activity

##  1 Introduction

Localized regions of intense bipolar magnetic structures, called active
regions (ARs), are seen on the solar surface. Their numbers vary periodically
in an 11-year cycle and trace the butterfly diagram, which reveals a cyclic
magnetic activity of the Sun in a latitude–time domain. Such diagrams have
proven to be useful and reveal some properties of the solar large-scale
magnetic field. The large-scale magnetic field reverses its polarity every
second such cycle, which constitutes the solar magnetic cycle of 22 years. The
origin of this global variability is not yet fully understood (see, e.g.,
Charbonneau, 2010).

In addition to the global solar magnetic field, it is observationally known
that there are ubiquitous small-scale magnetic fluctuations, the origin of
which is equally debated - they could arise through the action of a small-
scale dynamo instability or tangling of the large-scale magnetic field due to
turbulence driven by convection, but it is not self-evident how the former
mechanism could work in the solar plasma, where the conditions are in general
unfavourable for it (see, e.g. Brandenburg & Subramanian, 2005). These
magnetic fluctuations are usually studied by investigating the magnetically
quiet regions on the solar surface (hereafter referred to as the quiet Sun,
QS). Several observational studies of QS magnetism have been conducted, and it
has been argued that the QS magnetic field is independent of the solar cycle
(see, e.g., Kleint et al., 2010; Buehler et al., 2013; Faurobert & Ricort,
2015; Jin & Wang, 2015a, b). On the other hand, some other studies have
proposed that some dependence should exist (e.g. Lites et al., 2014; Meunier,
2018; Faurobert & Ricort, 2021). Taken the theoretical and observational
controversies, further studies of QS magnetism are required, and long-term
investigations are now enabled by instruments like Helioseismic and Magnetic
Imager (HMI, Scherrer et al., 2012; Schou et al., 2012) on board the Solar
Dynamics Observatory (SDO, Pesnell et al., 2012) This instrument provides high
sensitivity, high spatial resolution, long-term stability, and constant
conditions. Currently, the data covers 12 years of observations, enabling
solar-cycle-scale studies.

Another motivation to study the QS comes from recent observational
investigations Singh et al. (2016); Waidele et al. (2022) that have reported
strengthening of solar surface or the fundamental $f$-mode about one to two
days before the formation of ARs using different kinds of local helioseismic
techniques. Accompanied with numerical simulations that have given similar
indications Singh et al. (2014, 2015, 2020), this seems as a very promising
avenue of studying the origin of solar sub-surface magnetism and the mechanism
of active region formation. The observational studies, however, suffer from a
lack of proper calibration method against the QS. Singh et al. (2016), for
example, rely on QS regions on the opposite hemisphere to compare with the AR
$f$-mode. This method requires that a quiet patch exists on the other
hemisphere, and hence limits the number of ARs that can be included in the
hindcasting procedure. It is also prone to be affected by the probable
fluctuations in the QS $f$-mode level. Although the results look promising,
proper calibration with statistically sound QS level, not just a comparison of
a random QS patch on the other hemisphere, is necessary to prove the
robustness of these findings. Also, such a calibration procedure is necessary
for increasing the sample size. Building such a QS calibration data product is
one of the main aims of this study: we carefully identify the quietest regions
on the solar surface based on the level of magnetic activity observed in line-
of-sight (LOS) magnetograms that are readily available from HMI, and compute
the $f$-mode energy at the central meridian as function of latitude and time
with suitable averaging. Building such a data product for the $f$-mode,
requiring us to identify the most inactive regions on the solar surface,
allows us to extract statistics of the QS magnetism as well, which is the
second main aim of this study.

The paper is organized as follows: in Sect. 2 we describe the data, the
necessary steps to clean it, and the automated pipeline we built for
harvesting the data and compiling the end products, namely the QS magnetism
data products, and the QS and AR $f$-mode data. In Sect. 3 we discuss our
findings for the QS magnetism, $f$-mode, and the lastly present some AR data
with the QS calibration applied.

##  2 Observations

Our analysis is based on data from HMI@SDO. We use two standard data products:
(i) full-disk line-of-sight (LOS) magnetograms, computed every 720 s by
combining filtergrams obtained over a time interval of 1260 s (hmi.M_720s),
and (ii) full-disk LOS dopplergrams, computed every 45 s from six positions
across the nominal 6173.3 Å spectral line (hmi.V_45s). We processed the two
data sets in a semi-automatic pipeline (see Fig. 2), optimized for obtaining
reliable information about the magnetic field in the QS regions and for a
robust computation of the $f$-mode power from the dopplergrams. The left tree
in Fig. 2 describes the pipeline used for the magnetograms, the right tree for
the dopplergrams. In the following we use the following notation and
definitions. We denote the solar latitude with $\lambda$, longitude with
$\varphi$, both in the Stonyhurst coordinate system Thompson (2006), and time
with $t$.

Figure 1: Data pipeline. The leftmost tree illustrates the pipeline to collect
the magnetograms in the most quiet patches, while the rightmost one the
pipeline to track data from the quiet patches from the HMI database, and
compute the $f$-mode energy. The central path of AR processing is otherwise
equivalent to the QS pipeline, but there AR coordinates are sent for tracking
in the MPS (in Germany) cluster environment, and an additional step, namely QS
calibration, is performed in the end. Rectangular boxes represent analysis
functions, and ellipsoids the derived data products. The rectangles with
yellow frames stand for the corrective functions applied to the data.

###  2.1 Magnetograms

The first data product, full-disk line-of-sight (LOS) magnetograms, provides a
direct measurement of the variability of the QS magnetism during a solar
cycle. To enhance the signal-to-noise ratio in the LOS magnetograms we
performed a newly developed algorithm for spatial and temporal averaging: The
full-disk HMI magnetograms starting from 27-Apr-2010 and ending at 04-May-2022
were downloaded from the Joint Science Operation Center (JSOC) hosted at
Stanford University (http://jsoc.stanford.edu) to a temporary storage (see
Fig. 2, ’Mahti storage’) and tracked at full spatial resolution for 8 hours to
compensate for the solar rotation (’Full disk tracking’). We neglected
differential rotation due to the short tracking time. The step between tracked
sequences was 4 hours, so that a total of 6 tracked sequences were gathered
per one day, resulting in more than $25765$ tracked sequences (as of
04-May-2022).

From each tracked sequence, we extracted two data products by dividing the
visible solar disk between latitudes and longitudes from $-80\degree$ to
+$80\degree$ into (i) $64\times 64$ overlapping patches of $15\degree$ (in
solar latitude and longitude) and (ii) $180\times 180$ patches of $1\degree$.
Every of these patches therefore contains a space-time cube of LOS
magnetograms at full spatial resolution at a 12-minute cadence, from which we
compute the root-mean-square (rms) of the magnetic field strength, averaged
over the full field-of-view of the tracked cube and over the full 8 hour
period ($B_{\rm rms}$=$\sqrt{\langle B^{2}\rangle_{\rm LOS}}$). The
$15\degree$ patches (i) are large enough to cover several supergranulation
cells containing network and internetwork fields Rieutord & Rincon (2010) with
a typical size of 30–35 Mm (we refer to them as Nw cubes), and the $1\degree$
patches (ii) are small enough that some of them lie completely in the
internetwork (INw cubes). The statistics for each patch were stored as data
products (see also Fig. 2) including the information about latitudinal and
longitudinal position as well as the Carrington longitude for network and
internetwork. We refer to these maps as the Nw and the INw statistical maps.

###  2.2 Quiet region selection

From the statistics computed from the cubes, $B_{\rm rms}$ turned out to be
the best tracer for determining the magnetic activity level. It could clearly
distinguish between $15\degree$ patches containing active regions, plage,
enhanced network and quiet network. Also, it depicted very well the low-field
internetwork regions.

The analysis of the solar-cycle variation of the QS magnetism required a
careful selection of the most quiet regions, defined as being free of enhanced
solar activity. We therefore searched for the minimum value of $B_{\rm rms}$
in both, the Nw and the INw statistical maps, on a latitudinal grid with a
$10\degree$ spacing fulfilling the following additional criteria:

  1. (i)

the most quiet pixel must be within $\pm 10\degree$ around the central
meridian,

  2. (ii)

this pixel must belong to the 10% most quiet pixels of the month,

  3. (iii)

this pixel must be the most quiet pixel within a 4-day interval.

Criterion (i) was chosen to get the strongest magnetic field signal along the
central meridian. Criterion (ii) guarantees an equal distribution of quiet
pixels over the 12-year period of available HMI measurements, and criterion
(iii) ensures that the quiet pixels for the 1-month period do not originate
from the same supergranular structure, since the dynamical evolution time of
the supergranulation lies between 24 and 48 h Rieutord & Rincon (2010). As a
result of applying these criteria we obtained two time series of the $B_{\rm
rms}$ for the most quiet patches in the network and the internetwork regions.
Note that this selection also efficiently removes the 24 h modulation present
in the HMI magnetograms.

###  2.3 Correction for HMI sensitivity change

The temporal evolution of INw $B_{\rm rms}$ value revealed clearly a change in
the HMI observing mode, performed on 13-Apr-2016. On this day, HMI switched to
a more efficient observing mode (see Hoeksema et al., 2018, 2014; Couvidat et
al., 2016). By combining both HMI cameras to determine the vector-field
observables the cadence for full-disk magnetograms could be reduced from 135 s
(observational mode MOD-C) to 90 s (MOD-L). This reduced the noise level for
Stokes $V$ measurements by 17%, resulting in a decrease of the noise level in
the ${\rm LOS}$ magnetograms by 5%.

For the long-term study presented in this paper, we need to correct for this
sensitivity change. A very accurate correction method can be derived from the
INw time series: since it contains only the most quiet pixels over a certain
latitude region and time, the sensitivity change results in a step function.
The value of the step was determined by fitting a polynomial to the $B_{\rm
rms}$ values determined from the INw time series plus a Heaviside step
function, centered at the date of the mode change. We used the Bayesian
information criterion (BIC, Stoica & Selen, 2004) to determine the degree of
the polynomial, which lies between 1 and 6 for the various latitudes. We want
to note that the retrieved amplitude of the Heaviside step function is only
weakly dependent on the degree of the polynomial. This fitting is exemplified
in Fig. 2 for the solar latitude 0∘, where the minimum value for the BIC was
reached for a fit with a polynomial of degree 3. The so-determined amplitude
of the Heaviside step function is added to the $B_{\rm rms}$ data points after
13-Apr-2016 for all data presented in this paper.

Figure 2: Determination of the correction for the HMI sensitivity change: the
observing mode change on April 13 2016 causes a discontinuity in the level of
$B_{\rm rms}$ values of the internetwork data. The original data are displayed
in with the dark red and black dots, the corrected data with the light red
dots. The dashed line indicates the polynomial fit of degree 3 used to obtain
the offset.

###  2.4 Dopplergrams

The second HMI data product used in this paper are the ${\rm LOS}$ velocity
maps. The goal is to compute the energy contained in the surface gravity mode,
the so-called $f$-mode, and investigate whether variability during the solar
cycle is present. Such variations could be caused by the presence of non-
emerging sub-surface magnetic fields, as the $f$-mode is known to be strongly
affected by the presence of magnetic fields (see, e.g., Cally et al., 1994;
Cally & Bogdan, 1997; Singh et al., 2016). It is also very important to study
this question in more detail, for establishing a reliable calibration method
for a measurement of the claimed $f$-mode enhancement prior to AR emergence
Singh et al. (2016); Waidele et al. (2022).

The dopplergram data is hosted in the German Data Center for SDO (GDC-SDO) on
a server at the Max Planck Institute for Solar System Research (MPS Göttingen,
Germany) whereas the analysis is executed in the CSC supercomputing
environment (Finland). We have utilized a function-as-a-service client based
on funcX Chard et al. (2020) for accessing the required data in the database
server. We have developed functions based on funcX API and deployed the
functions in the MPS environment. These functions leverage the
mtrack111http://hmi.stanford.edu/teams/rings/modules/mtrack/v25.html command
to prepare the dopplergram cubes within the GDC-SDO environment and
subsequently transfer them to the CSC environment. The coordinates of the
selected QS regions were sent to the funcX service, which invokes suitable
functions to automate the data retrieval and movement. The data processing
pipeline used Astropy222http://www.astropy.org, a community-developed core
Python package for Astronomy (Astropy Collaboration et al., 2013, 2018) and
version 3.1.6 of the SunPy333https://sunpy.org open source software package
The SunPy Community et al. (2020).

###  2.5 Computation of the $f$-mode energy

Subsequent processing involved calculating the three-dimensional power spectra
for each dopplergram cube and integrating in $(k_{x},k_{y})$-plane, for each
angular frequency $\omega$, over all wavenumbers,
$k=\sqrt{k_{x}^{2}+k_{y}^{2}}$, to obtain a collapsed power spectrum,
$P(\omega,k)$. An illustrative case for a frequency value at which $f$-mode is
strong, is shown in Fig. 3, upper panel. Such a procedure significantly
reduced the noise level leading to smooth one-dimensional $k-\omega$ spectra,
an example being shown in the lower panel of Fig. 3; in this kind of a plot,
the $f$-mode is the rightmost peak at the highest $k$ values. The adopted
procedure is justified for the QS spectra, as the ring diagrams are radially
symmetrical w.r.t. $\omega$ axis. To obtain the total energy contained in the
$f$-mode, $E_{f}$, we performed another integration over the separated
$f$-mode signal. In contrast to earlier studies (Singh et al., 2016; Waidele
et al., 2022), we did not perform fitting to the $f$-mode for its extraction.
Instead, we only determined the background signal, and subtracted this
contribution. The background determination used the fact that the position of
the $f$-mode peak for low values for $\omega$ (between 3 and 6 rad s-1) lies
outside the $k$ range used for the computation of the $f$-mode power. Since
the background level does not change significantly with $\omega$, the average
spectrum over this $\omega$ range provided a good estimate for the background
signal. Then, an integration range in $k$ space was selected for each constant
$\omega$ in a following way: first the maximum of $f$-mode $k_{\rm max}$ and
the minimum between $f$-mode and first $p$-mode $k_{\rm start}$ were detected;
$k_{\rm start}$ was chosen as the start of integration range and the end of
the integration range $k_{\rm end}$ was set as $k_{\rm end}=k_{\rm
max}+2(k_{\rm max}-k_{\rm start})$; $k_{\rm end}$ chosen in such a way
guarantees that the integration range is sufficiently wide to cover
significant part of the $f$-mode signal in $k$-space. Integration range in
$\omega$ space was chosen between values starting from 14.45 and ending at
29.03, where the significant part of the $f$-mode power resides. As the end
result, we have computed the $f$-mode energy

$E_{f}=\sum_{\mbox{$k_{\rm start}$}}^{\mbox{$k_{\rm
end}$}}\sum_{\mbox{$\omega_{\rm min}$}}^{\mbox{$\omega_{\rm
max}$}}P(\omega,k)$ (1)

After completing the described steps in the CSC supercomputing environment, we
have collected a total of 22680 $f$-mode energies for patches close to central
meridian over all latitudes covering most of SC24 and the ascending phase of
SC25.

Figure 3: Top: example of quiet-sun ring diagram at the angular frequency
$\omega=22.40$ rad s-1 for one quadrant of the full ring (logarithmic color
scale), computed from the data cube presented in Fig. 4 (left panel). Bottom:
collapsed spectrum obtained by integrating over the ring diagram. The red
vertical line marks the position of the maximum of the $f$-mode and the black
vertical lines the range of integration. The wavenumber $k$ is expressed in
units of the circumference of the Sun ($C_{\rm Sun}$).

#### Orbital correction of the $f$-mode energy

Since the $f$-mode energy is computed from the ${\rm LOS}$ dopplergrams, its
value depends strongly on the viewing geometry. This dependence, roughly
following the cosine of the solar latitude for data taken at the central
meridian, is additionally modulated by the orbital motion of the Earth around
the Sun, which changes the viewing angle at any given solar latitude by
$\approx\pm 7^{\circ}$ during one year. We compensated for this periodic
variation by fitting the parameters ($x_{0}(\lambda),...,x_{3}(\lambda)$) of
the following function to all observations of a given latitude $\lambda$:

$\displaystyle A_{\mbox{corr}}(\lambda)$ $\displaystyle=$ $\displaystyle
x_{0}(\lambda)(2((1+\cos(\Phi+x_{1}(\lambda)))/2)^{x_{3}(\lambda)}-1)$ (2)
$\displaystyle+$ $\displaystyle x_{2}(\lambda),$

with $\Phi$ being the phase angle of the Earth defined as the cotangens
computed from the $x,y$ barycentric position of the Earth. This correction
$A_{\mbox{corr}}(\lambda)$ is then subtracted from the computed $E_{f}$
values, which efficiently removes any yearly variation.

##  3 Results

After applying the selection criteria described in Sec. 2.2 and the subsequent
correction for the HMI’s sensitivity change (Sec. 2.3) we obtain the
dependency of the $B_{\rm rms}$ values of the most quiet region at the central
longitude as a function of heliographic latitude and time. Fig. 5 presents the
$B_{\rm rms}$ values from May 2010 until May 2022 for the Stonyhurst latitude
and longitude $(\lambda,\varphi)=(0,0)^{\circ}$. The individual data points
(red dots) are computed over an area of 15∘ in latitude and longitude and over
a time of 8 hours. At disk center, this corresponds to an area of $\approx
180\times 180$ Mm2, and therefore contains 30–40 supergranular cells. The
$B_{\rm rms}$ value therefore contains network and internetwork fields.
Typical magnetograms at disk center at solar minimum (July 2010) and maximum
(Nov 2014) are presented in Fig. 4. There is no obvious difference between the
two maps: both show similar minimum and maximum field strengths, and the
visible impression of the supergranular cells with the strong field patches of
both polarities surrounding the internetwork regions is indistinguishable.

Fig. 5 clearly indicates the solar cycle dependence of the QS magnetic fields,
consisting of network and internetwork. The $B_{\rm rms}$ values peak in the
second half of 2014, around half a year later compared to the declared maximum
of solar cycle #24 (April 2014, source: WDC-SILSO, Royal Observatory of
Belgium, Brussels). The variation is statistically significant, as indicated
by the the mode and the 97% percentiles, computed by fitting a log-normal
distribution over a 1-year sliding window.

A similar plot is presented in Fig. 6, but now the $B_{\rm rms}$ value was
computed only for a $1\degree$ window in latitude and longitude, corresponding
to an area of only $12\times 12$ Mm2 at disk center. The quiet-region
selection criteria (Sec. 2.2) guarantees, that this patch lies fully within
the internetwork, and is not ‘contaminated’ by network fields. Clearly, there
is no dependence of $B_{\rm rms}$ with solar cycle detectable. The slight
increase over the 12-year period might be an effect of the aging of the HMI
instrument and the resulting decrease in the sensitivity for magnetic field
measurements.

The average $B_{\rm rms}$ value in the internetwork is in $\approx 6$ G. Using
the area of one HMI pixel this corresponds to a magnetic flux of $3\times
10^{16}$ Mx. This value is certainly affected by the sensitivity of the HMI
instrument, and is reduced to $\lesssim 5$ G after the HMI sensitivity change
described in Sect. 2.3.

Figure 4: Comparison of two disk-center magnetograms at solar minimum
(2010-07-21) and around the maximum seen in the disk-center $B_{\rm rms}$
values (2014-11-09). The maps show the first frame of the 8 h data cube.
Figure 5: $B_{\rm rms}$ determined from the most quiet patches from 2010 to
2022 at disk center. The patches of 15∘ size in longitude and latitude contain
network and internetwork fields. The solid red lines display the 97%
percentile level, the dashed, gray line the mode of a log-normal fit to
yearly-binned data moved in a sliding window of 100 days length. For all
values of $\omega$ the magnitude of $P$ is always an order of magnitude higher
than the noise level outside the considered location of the ring.  Figure 6:
Same as Fig. 5, but for $B_{\rm rms}$ determined from patch size of only of 1∘
in longitude and latitude, corresponding to internetwork patches.

###  3.1 Quiet Sun $f$-mode butterfly diagram

To investigate whether the $f$-mode energy has any significant solar cycle
dependence, we construct butterfly-diagram-like plots by computing and
removing the temporal mean, and plotting the so-obtained data in a binned
latitude-time diagram. We compute the temporal mean over the whole time series
as

$\overline{E_{f}}(\lambda)=1/n_{t}\sum_{t}E_{f}(\lambda,t),$ (3)

where $n_{t}$ is the total number of time points at a certain $\lambda$, and
denote it with an overbar. Next, we compute the variation of the $f$-mode
energy around this level, defined as

$E_{f}^{\prime}(\lambda,t)=E_{f}(\lambda,t)-\overline{E_{f}}(\lambda),$ (4)

and plot this quantity for the central meridian as function of latitude and
time as in Figure 7, in patches with a 3-month binning in time and roughly 6
degrees in latitude.

To estimate the overall variability level of the QS, and compare it to the
cycle variation level defined above, we will use the standard deviation of the
data in a patch. We take $N_{t}$ bins in time, $N_{\lambda}$ in latitude, and
denote the number of $E_{f}$measurements contained in such a patch by $N_{P}$.
We compute the average $f$-mode energy in a patch as

$\langle E_{f}\rangle=1/N_{P}\sum_{N_{P}}E_{f}(\lambda,t).$ (5)

Then we can straightforwardly calculate the standard deviation as

$\sigma=\sqrt{\frac{1}{N_{P}}\sum_{\mbox{patch}}(E_{f}(\lambda,t)-\left<{E_{f}}\right>)^{2}}.$
(6)

From this we form the average error at each latitude by taking an average over
time and denote it as $\overline{\sigma}$. This quantity is plotted as the
error bar in the right panel of Fig. 7. By inspecting the level of variability
around the temporal mean and the overall variability of the QS, we must
conclude that these signals are of the same order of magnitude. In this
respect, the solar cycle dependent signal is very weak, but if we compare that
to the noise level in the data, both the QS variability level and cycle-
dependent signal exceed the noise level by an order of magnitude, hence this
is not an issue concerned with the instrument sensitivity.

What could be the cause of the strong QS variability level is most likely
related to the strongly varying sub-surface magnetic fields that affect the
$f$-mode, invisible in the surface $B_{\rm rms}$ which was used as the
selection criterion for really quiet patches. This is illustrated in Fig. 9,
where we show the measured $f$-mode energy versus the $B_{\rm rms}$ in the
patches. The $f$-mode energy is not correlated with the $B_{\rm rms}$ seen at
the surface, and strong, nearly constant, variability is seen at different
epochs of the solar cycle. This poses clear limitations of the calibration
method proposed here to detect weak transient signals in the $f$-mode
evolution.

Figure 7: Butterfly diagram of $f$-mode energy variations,
$E^{\prime}_{f}(\lambda,t)$, ToDo: we should still describe the binning.
computed from the collapsed ring diagrams (see bottom panel of Fig. 3). The
temporal average, $\overline{E_{f}}(\lambda)$, has been subtracted from the
data, and is shown as the black line in the right panel. Its error bars
represent the time-averaged standard deviation of the fluctuations of the
$f$-mode energy around the signal for every bin
($\overline{\sigma}(\lambda)$).  Figure 8: Top: Cut through butterfly diagram
of $f$-mode energy variation for three different latitudes (see legend). The
solid lines represent a fit using a sine function. Bottom: NOAA sunspot
number.

Figure 7 displays the butterfly diagram of the QS $f$-mode energy determined
from $15\degree\times 15\degree$ patches that showed only weak, small-scale
magnetic structures. To compute the $f$-mode energy, we have used the whole
ring diagram data. Temporally averaged profile in latitude is removed for
highlighting the variation in the $f$-mode power. The period covered is
somewhat larger than the 11 year SC24, which began to show some magnetic
activity from early 2010; this cycle ended by around 2020 (see also the lower
panel of Figure 8, showing the sunspot number of SC24 and SC25 up to the
writing of this manuscript). Broadly, we find an anti-correlation of the
$f$-mode energy with the activity cycle of the Sun — during the solar minimum
and early rising phase ($\sim 2010-12$), the $f$-mode energies are larger, and
around cycle maximum to the declining phase, the $f$-mode energy is
suppressed. This correlation is strongest near the solar equator, within about
$\pm 20\degree$ in latitude, and it decreases with higher latitudes.
Interestingly, we observe a temporal shift in the $f$-mode butterfly diagram
with respect to the solar cycle, such that the minimum of the $f$-mode occurs
later, i.e., in the declining phase of SC24, roughly around the end of the
year 2016.

During mid-2017, we see weak strengthening of the $f$-mode strength to
commence, first occurring at higher Southern latitudes. During the minimum
(2018-2019) between the SC24 and 25, the strengthening appears stronger in the
higher latitudes. It also continues after the minimum, when SC25 has already
started its ascending phase. What is also noteworthy is that the $f$-mode
energy never rises as high as the values observed during the ascending phase
of SC24 (2011-2012).

In Figure 8 we present cuts through the butterfly diagram at three different
latitudes, equator, and $\pm 29\degree$, and sinusoidal fits to the data
points. The variability is the strongest at the equator, but also contains the
largest scatter. The minimum $f$-mode energy is obtained there around the
beginning of the year 2018. There, the difference in magnitude during SC24 and
25 is the largest. Southern higher latitudes show very weak variation in
contrast to the equator, while the Northern one with an intermediate
magnitude. The shift w.r.t. the solar cycle is somewhat latitude dependent. At
the higher northern latitude, the $f$-mode minimum occurs roughly a year
earlier than at the equator. The shift in the $f$-mode mode phase w.r.t. the
sunspot cycle is of the order of three to four years, hence it is more in
phase with the poloidal phase of the magnetic field than the toroidal one.

###  3.2 AR $f$-mode with QS calibration

In this section we compare the $f$-mode time evolution before and after the
emergence of two ARs (11130 in Figure 10 and 11105 in Figure 11), that were
reported to show an enhancement both in Singh et al. (2016) and Waidele et al.
(2022) w.r.t. a QS control patch in the opposite hemisphere several days
before their emergence. The main difference in our analysis is that the
$f$-mode energy is now normalized to the actually measured average QS level
around the specific time and latitude. To minimize the QS fluctuations,
Gaussian smoothing to the QS $E_{f}$ is applied, but the smoothing kernel
width used is kept at $5^{\circ}$ in longitudinal and latitudinal direction so
that the fit to the data can be considered accurate. $f$-mode energy of 1 is
hence equivalent to that of the QS, and values stronger/weaker indicate
enhancement/quenching of the $f$-mode w.r.t. QS. Another difference is that we
use the full ring diagram when computing the $f$-mode energy, while the other
papers used only $k_{x}=0$ cuts. We have, however, performed analysis with
$k_{x}=0$ and $k_{y}=0$ cuts, and see no significant difference between the
two, nor w.r.t. the full ring data, except for the increased noise level in
$E_{f}$ derived from single cuts. As the measurement is done exactly in the
same manner for the QS and ARs, we anyway anticipate that this difference
should not influence the results significantly. Also, as more integration is
performed, the less noisy is the signal, hence improving the data quality.

The results from our AR analysis are shown with green symbols and lines in
Figure 10 and Figure 11: in the top panel, each orange point shows the
normalized $f$-mode energy, $\tilde{E_{f}}$, with four-hour cadence, smoothed
with a 24-hour boxcar to remove a daily fluctuation otherwise prominently
present in the HMI data. The shaded areas show the QS variance,
$E^{\prime}_{f}$, determined from the standard deviation of the variation
presented in Fig. 9, is roughly constant (around 4–5 percent) at all latitudes
and longitudes. When it is used as normalization of the $f$-mode, then the
uncertainty at large longitudes becomes very high, as is indicated by the
rapidly fanning shaded areas toward the limbs. The $B_{\rm rms}$ measured
simultaneously at the surface is shown in the lowest panels of the figures
with black symbols.

In the middle panels of Figure 10 and Figure 11, we compare the different
calibration curve candidates. The calibration curve used by Singh et al.
(2016), based on more simple cosinus dependence of the form
$\cos{\alpha}\left[q+\left(1-q\right)\cos{\alpha}\right]$ with $q=0.5$ and
$\cos\alpha=\cos\lambda\cos\varphi$, is shown with blue symbols and lines in
the aforementioned panels. As this specific curve does not fit the QS data
very well, we also performed a fit and determined the optimal value of $q$,
that turns out to be $-0.045$; this curve and calibration with it are plotted
with orange symbols and lines in the top and middle panel of the figures. The
differences between these calibration curves are illustrated in the middle
panels: The Singh et al. (2016) calibration curve does not fit the data,
overplotted with gray dots, well at larger longitudes. There, it is decreasing
less strongly than the data. Near the disk center the amplitude is too low to
present the data well. The new cosinus fit and the Gaussian convolution
applied to the QS data represent it better, with the consequence that the
calibration curves decrease more steeply toward the limbs. The fitted cosinus
curve has somewhat larger amplitude than the Gaussian convolution. This has
the effect that the $\tilde{E_{f}}$ calibrated with the cosinus fit tends to
be always below the QS level. The Gaussian convolution calibration results at
$\tilde{E_{f}}$ values maximally at the QS level before the AR emergence.

As can be seen from Figure 10 and Figure 11 top panels, applying the new QS
calibration curve flattens the signal profile in comparison to the earlier
$q=0.5$ calibration used. For both of the selected ARs, the $f$-mode energy
close to the limb is below the QS level to start off with, which could
indicate that both of these active regions occur in an environment already
affected by a sub-surface magnetic field. Applying the earlier calibration
curve, a clear enhancement above the QS level is visible for AR11130, while
for AR11105 the profile is similar, but the enhancement amplitude remains
within the level of QS variability (shaded area). The new calibration shows a
more modest enhancement, reaching its peak value at similar times as with the
earlier calibration for AR11130, while somewhat earlier for AR11105, but its
amplitude is not exceeding QS level significantly. Unfortunately, the possible
enhancement of $f$-mode energy of the selected ARs occurs so close to the limb
that it is not completely certain whether it is an artefact of the rather
uncertain limb calibration data, or a true effect. Obviously, a larger sample
size should be studied, including also ARs that emerge closer to the west limb
to completely prove/refute this scenario.

Figure 9: $f$-mode energy as a function of $B_{\rm rms}$ for the most quiet
patches at disk center. The colors indicate the grouping according to the
years. The standard deviation for the $f$-mode variation is 87. Figure 10:
Calibrated $f$-mode strength ($\tilde{E_{f}}$, top), QS calibration data
($E_{f}$, gray dots) and fits (middle) and $B_{\rm rms}$ (bottom) for AR11130.
The line-plot colors in the top and the middle panel indicate
$\tilde{E_{f}}$and $E_{f}$for various calibration methods as indicated in the
legend. The shaded area in the top panel represents the standard deviation of
the QS calibration derived from Fig. 9 and scaled with the fitted calibration
functions. Figure 11: Same as Fig. 10, but for AR11105.

##  4 Discussion and conclusions

The first part of the paper presented an analysis of QS magnetism, where the
main findings were the solar cycle independence of the internetwork magnetic
fluctuations, studied by using $1\degree$ patches, while network and
internetwork fields, studied by using $15\degree$ patches, were found to show
a statistically significant solar cycle dependence. The maximum of the
magnetic fluctuations at the disk center was observed to exhibit a phase shift
of roughly half a year w.r.t. to the sunspot-number-defined cycle. We
interpret these findings as follows. The independence of the internetwork
magnetic fields is either reflecting the fact that these fields, clearly on
smaller scales than the supergranular cells, are generated by the fluctuation
(or small-scale) dynamo alone, or then the sensitivity of the HMI instrument
prevents one to see any variation. The solar cycle dependence of the network
and internetwork fields reflects the fact that two separate processes are
responsible for generating these fields. Turbulence is tangling the solar-
cycle dependent large-scale sub-surface field resulting in magnetic
fluctuations that have the same dependence, and these are then expelled to the
edges of the supergranular cells. The magnetic fluctuations generated by the
fluctuation dynamo also become expelled in a similar fashion, but as the
growth of these fluctuations is exponential due to the dynamo instability, the
fluctuation dynamo can quickly replenish the fluctuations also in the
internetwork regions. The tangling process can be envisioned to have an
amplification time scale that is slower, only linear over time.

As is extensively reviewed in Brandenburg & Subramanian (2005), in turbulent
dynamo theory, the magnetic fields grow under the constraint of magnetic
helicity conservation, which leads to a bihelical spectrum of magnetic
helicity. Signs of such bihelical nature of the solar surface magnetic field
have already been reported by Brandenburg et al. (2017); Singh et al. (2018).
Accumulation of small-scale magnetic helicity leads to the quenching of large-
scale dynamo. In order to further grow its large-scale magnetic field, the
system must shed its small-scale helicity. Being an open system, the Sun may
have fluxes of magnetic helicity where ARs could play a vital role in removing
the magnetic helicity from small-scales, thus leading to a rejuvenation of the
large-scale dynamo. This mechanism was proposed to lead to a phase shift
between the cycles for sunspot and magnetic helicity, that was seen to peak
half a year later than the sunspot cycle in Singh et al. (2018). We see
similar evidence here for the solar-cycle dependent part of the magnetic
fluctuations in the network and internetwork - the phase shift by half a year
is consistent with the large-scale field being rejuvenated later, when ARs
have already peaked, and shed some small-scale magnetic helicity out from the
system to allow for the growth of the global field.

The second part of our study built a novel method for calibrating the $f$-mode
energy using the very quiet patches harvested in the first part. Even though
more than 22 000 patches were analysed to obtain the statistics over the solar
cycle, the $f$-mode energy variations of the QS remained large, and they show
no clear correlation with the rms magnetic field strength measured at the
surface. We interpret this as an indication of the $f$-mode energy being
affected by sub-surface magnetic fields, and the surface rms magnetic field
being an inadequate indicator to find the QS level. The measured QS $f$-mode
energy, indeed, shows a solar-cycle variation, which is, in broad terms, anti-
correlated with the solar cycle, the $f$-mode energy being lower/higher when
the rms magnetic field is strong/weak. Such a behavior is expected in the
light of many previous investigations, which have suggested $f$-mode damping
to occur as a result of absorption of the $f$-mode mode by the magnetic fields
in ARs (Cally et al., 1994; Cally & Bogdan, 1997; Singh et al., 2016). This
effect is the dominant one over the possible enhancement prior to AR
emergence, as can be seen also in the hindcasted AR data in this study.
Combined with the heavy averaging, any possible enhancement will be washed out
from a statistic measuring the global behavior.

The phase shift of several years between the sunspot cycle and the $f$-mode
modulation is too large to be explicable by the scenario suggested above for
the half a year phase shift of the QS magnetism maximum w.r.t. the sunspot
cycle. The $f$-mode energy modulation appears to follow more closely the
poloidal component of the global field than the toroidal part that is commonly
thought to give rise to the ARs. These two field components are known to have
a systematic phase difference of roughly $\pi$ (see, e.g., Charbonneau, 2010).

The general trend of $f$-mode energy being reduced during SC25 in comparison
to SC24 could be a result of the instrument aging, but the clear differences
seen over different latitudes and hemispheres is speaking against of an
overall instrumental degradation. If it were a physical effect, it would
suggest that the sub-surface magnetic field is stronger during the ascending
phase of SC25 than the one of SC24, and that could lead to SC25 being stronger
than SC24. Further investigations and fine tuning of the method are required
to verify the results and conclusions reached in this study, but the QS
$f$-mode energy seems to have potential applications in diagnosing the sub-
surface magnetic fields in the Sun.

The third and last part of our study presented two hindcasted AR $f$-mode
signals, that have been proposed to show a clear enhancement of the $f$-mode
energy in the earlier studies (Singh et al., 2016; Waidele et al., 2022), with
our novel QS calibration method. In summary, our analysis shows that the
transient and weak enhancement signal is very sensitive to the calibration
method used. The new method for the QS calibration developed here tends to
show a more flat signal with mild enhancement only in comparison to the
earlier calibration method of Singh et al. (2016). Although in general a
signature of an enhancement is seen, it remains unclear whether it is an
effect caused by the uncertainties close to the limb or a real effect. In the
light of the QS data, the magnitude of the enhancement is weaker than the QS
variability level, and hence its detection and usage as an AR emergence
predictor is not possible with the QS calibration method.

###### Acknowledgements.

We are indebted to Dr. Harsha Raichur, the developer of the original $f$-mode
analysis pipeline, wherefrom this work stems. We acknowledge the fruitful
discussions with Prof. Nishant Singh on the theoretical interpretation of the
results and practical insights to the analysis. All SDO data used are publicly
available from the Joint Science Operations Center (JSOC) at Stanford
University supported by NASA Contract NAS5- 02139 (HMI), see
http://jsoc.stanford.edu/. The data analysis has been carried out on
supercomputers in the facilities hosted by the CSC—IT Center for Science in
Espoo, Finland, which are financed by the Finnish ministry of education. The
data were also processed at the German Data Center for SDO (GDC-SDO), funded
by the German Aerospace Center (DLR), and hosted by the Max Planck Institute
for Solar System Research (Göttingen, Germany). This project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (Project UniSDyn, grant agreement
n:o 818665).

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Brandenburg et al. (2017) Brandenburg, A., Petrie, G. J. D., & Singh, N. K. 2017, ApJ, 836, 21, doi: 10.3847/1538-4357/836/1/21
  * Brandenburg & Subramanian (2005) Brandenburg, A., & Subramanian, K. 2005, Phys. Rep, 417, 1, doi: 10.1016/j.physrep.2005.06.005
  * Buehler et al. (2013) Buehler, D., Lagg, A., & Solanki, S. K. 2013, A&A, 555, A33, doi: 10.1051/0004-6361/201321152
  * Cally & Bogdan (1997) Cally, P. S., & Bogdan, T. J. 1997, ApJ, 486, L67, doi: 10.1086/310833
  * Cally et al. (1994) Cally, P. S., Bogdan, T. J., & Zweibel, E. G. 1994, ApJ, 437, 505, doi: 10.1086/175014
  * Charbonneau (2010) Charbonneau, P. 2010, Living Reviews in Solar Physics, 7, 3, doi: 10.12942/lrsp-2010-3
  * Chard et al. (2020) Chard, R., Babuji, Y., Li, Z., et al. 2020, in (ACM), doi: 10.1145/3369583.3392683
  * Couvidat et al. (2016) Couvidat, S., Schou, J., Hoeksema, J. T., et al. 2016, Sol. Phys., 291, 1887, doi: 10.1007/s11207-016-0957-3
  * Faurobert & Ricort (2015) Faurobert, M., & Ricort, G. 2015, A&A, 582, A95, doi: 10.1051/0004-6361/201526298
  * Faurobert & Ricort (2021) —. 2021, A&A, 651, A21, doi: 10.1051/0004-6361/202140705
  * Hoeksema et al. (2018) Hoeksema, J. T., Baldner, C. S., Bush, R. I., Schou, J., & Scherrer, P. H. 2018, Sol. Phys., 293, 45, doi: 10.1007/s11207-018-1259-8
  * Hoeksema et al. (2014) Hoeksema, J. T., Liu, Y., Hayashi, K., et al. 2014, Sol. Phys., 289, 3483, doi: 10.1007/s11207-014-0516-8
  * Jin & Wang (2015a) Jin, C., & Wang, J. 2015a, ApJ, 806, 174, doi: 10.1088/0004-637X/806/2/174
  * Jin & Wang (2015b) Jin, C. L., & Wang, J. X. 2015b, ApJ, 807, 70, doi: 10.1088/0004-637X/807/1/70
  * Kleint et al. (2010) Kleint, L., Berdyugina, S. V., Shapiro, A. I., & Bianda, M. 2010, A&A, 524, A37, doi: 10.1051/0004-6361/201015285
  * Lites et al. (2014) Lites, B. W., Centeno, R., & McIntosh, S. W. 2014, PASJ, 66, S4, doi: 10.1093/pasj/psu082
  * Meunier (2018) Meunier, N. 2018, A&A, 615, A87, doi: 10.1051/0004-6361/201730817
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Rieutord & Rincon (2010) Rieutord, M., & Rincon, F. 2010, Living Reviews in Solar Physics, 7, 2, doi: 10.12942/lrsp-2010-2
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207, doi: 10.1007/s11207-011-9834-2
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Singh et al. (2015) Singh, N. K., Brandenburg, A., Chitre, S. M., & Rheinhardt, M. 2015, MNRAS, 447, 3708, doi: 10.1093/mnras/stu2540
  * Singh et al. (2014) Singh, N. K., Brandenburg, A., & Rheinhardt, M. 2014, ApJ, 795, L8, doi: 10.1088/2041-8205/795/1/L8
  * Singh et al. (2018) Singh, N. K., Käpylä, M. J., Brandenburg, A., et al. 2018, ApJ, 863, 182, doi: 10.3847/1538-4357/aad0f2
  * Singh et al. (2016) Singh, N. K., Raichur, H., & Brandenburg, A. 2016, ApJ, 832, 120, doi: 10.3847/0004-637X/832/2/120
  * Singh et al. (2020) Singh, N. K., Raichur, H., Käpylä, M. J., et al. 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 196, doi: 10.1080/03091929.2019.1653461
  * Stoica & Selen (2004) Stoica, P., & Selen, Y. 2004, IEEE Signal Processing Magazine, 21, 36, doi: 10.1109/MSP.2004.1311138
  * The SunPy Community et al. (2020) The SunPy Community, Barnes, W. T., Bobra, M. G., et al. 2020, The Astrophysical Journal, 890, 68, doi: 10.3847/1538-4357/ab4f7a
  * Thompson (2006) Thompson, W. T. 2006, A&A, 449, 791, doi: 10.1051/0004-6361:20054262
  * Waidele et al. (2022) Waidele, M., Roth, M., Singh, N., & Käpylä, P. 2022, arXiv e-prints, arXiv:2202.11236.  https://arxiv.org/abs/2202.11236