11institutetext: Max-Planck-Institut für Sonnensystemforschung, Justus-von-
Liebig-Weg 3, 37077 Göttingen, Germany
11email<EMAIL_ADDRESS>22institutetext: Georg-August-Universität,
Friedrich-Hund-Platz 1, 37077 Göttingen, Germany

# Contribution of flows around active regions to the north-south helioseismic
travel-time measurements

P.-L. Poulier 11 Z.-C. Liang 11 D. Fournier 11 L. Gizon 1122

(Received ; )

###### Abstract

Context. In local helioseismology, the travel times of acoustic waves
propagating in opposite directions along the same meridian inform us about
horizontal flows in the north-south direction. The longitudinal averages of
the north-south helioseismic travel-time shifts vary with the sunspot cycle.

Aims. We aim to study the contribution of inflows into solar active regions to
this solar-cycle variation.

Methods.  To do so, we identify the local flows around active regions in the
horizontal flow maps obtained from correlation tracking of granulation in
SDO/HMI continuum images. We compute the forward-modeled travel-time
perturbations caused by these inflows using 3D sensitivity kernels. In order
to compare with the observations, we average these forward-modeled travel-time
perturbations over longitude and time in the same way as the measured travel
times.

Results.  The forward-modeling approach shows that the inflows associated with
active regions may account for only a fraction of the solar-cycle variations
in the north-south travel-time measurements.

Conclusions. The travel-time perturbations caused by the large-scale inflows
surrounding the active regions do not explain in full the solar-cycle
variations seen in the helioseismic measurements of the meridional
circulation.

###### Key Words.:

Sun: activity – Sun: helioseismology

##  1 Introduction

The Sun’s meridional flow at the surface is poleward with a maximum amplitude
of about $15$ m/s (Duvall 1979). The meridional ciculation, both at the
surface and in the deep convection zone, is believed to be a key ingredient in
flux-transport dynamo models (e.g., Wang et al. 1991; Dikpati & Gilman 2006,
2009). Observationally, the meridional circulation is given as the
longitudinal average of the north-south flows. This longitudinal average is
not constant in time: its amplitude and latitudinal dependence change over the
solar cycle (e.g., Komm et al. 1993; Hathaway & Rightmire 2010, 2011). It has
been proposed that extended inflows around solar active regions (Gizon et al.
2001) modulate the meridional flow at the surface (Gizon 2004; Gizon et al.
2010). These inflows have amplitudes of up to $50$ m/s near the surface and
extend up to $10^{\circ}$ from the center of the active regions or further
(Gizon et al. 2001; Hindman et al. 2003; Haber et al. 2004; Hindman et al.
2009; Braun 2019; Gottschling et al. 2021).

Helioseismic travel-time shifts in the north-south direction are sensitive to
the meridional flow (e.g., Giles et al. 1997). The measurements made by Gizon
et al. (2020) show a solar-cycle modulation. We aim to determine how much of
this modulation may indeed be attributed to the near-surface active-region
flows. To this end, we isolate and measure the active-region flows in local
correlation tracking (LCT) flow maps. By assuming a depth dependence of these
flows, we compute forward helioseismic travel-time perturbations in the north-
south direction to estimate their contribution to the fluctuations seen in the
time-distance measurements. Section 2 presents the data and the method.
Section 3 shows the latitudinal and longitudinal components of the resulting
active-region flows. We compute in Section 4 the north-south forward
helioseismic travel-time perturbations associated with these flows, using 3D
Born sensitivity kernels. We attempt to model the inflows with a simple model
based on the latitudinal gradient of the unsigned magnetic field in Section 5,
in order to extend the analysis to the previous solar cycle. We compare our
results with helioseismic measurements in Section 6.

##  2 Horizontal flows from granulation tracking

###  2.1 Flow maps from LCT

We use the horizontal flow maps computed by Löptien et al. (2017). The
original data set covered the period from May 2010 to April 2016, and has
later been extended till April 2019 (Gottschling et al. 2021). The maps were
obtained by using the Fourier Local Correlation Tracking code (Welsch et al.
2004; Fisher & Welsch 2008). The code tracked pairs of consecutive full-
resolution intensity images from SDO/HMI (Schou et al. 2012), so that the
flows represent the surface motions of solar granulation. The cadence of the
flow maps is 30 minutes.

The data contained systematics like the orbital motions of SDO and the
shrinking-Sun effect (Lisle & Toomre 2004; Löptien et al. 2016). Löptien et
al. (2017) decomposed the flow maps into Zernike polynomials and filtered the
time series of the coefficient amplitudes to remove the mean of the time
series, the periods of 24 h, 1 year, and their corresponding harmonics, which
takes care of most of the systematics.

The resulting data product contains the time-varying part of the rotation
(torsional oscillations) and of the meridional circulation, plus potentially
residual systematics (Gottschling et al. 2021). The velocities are in CCD-
frame units (pixels per second). The size of these filtered images is
$1024\times 1024$ pixels; we perform a $5\times 5$ binning on these images,
yielding a spatial resolution of about 10 arcsec per pixel which corresponds
to 0.6 heliographic degree per pixel at disk center.

###  2.2 Construction of background in CCD frame

Gottschling et al. (2021) reported that the LCT flow maps contain residual
systematics. This means that large-scale background flows shall be removed
from the flow maps in order to measure the flows related to active regions.
These systematics likely depend on the position on the visible disk, in
particular if they are a residual of the shrinking-Sun effect. Therefore, the
background will be measured in the CCD frame; the procedure is described
below.

####  2.2.1 Contours around active-region flows

We define the background as the area located far away from the magnetic
activity and its related flows. In order to measure the background flows, we
need to identify the active regions first. We use the HMI magnetograms at the
same time steps as the LCT data. We assume that the magnetic field is purely
radial and thus divide the line-of-sight magnetic field $B_{\mathrm{LOS}}$ by
the cosine of the great-circle angular distance to disk center to obtain
$B_{r}$.

We aim to determine smooth contours around the active-region flows based on a
spatially-smoothed magnetic field. We use a 2D Gaussian kernel to smooth
$|B_{r}|$, with a full width at half maximum (FWHM) of $10\degr$ (here,
$1\degr$ is $1$ heliographic degree) because it is the typical horizontal
extent of the flows into active regions (see Fig. 10). We then define the
contours around the active regions as the lines along which the smoothed
$|B_{r}|$ is equal to a magnetic threshold $b$. We choose $b$ to be $3.5$ G
(see Appendix A). By aligning the centers of active regions with a total
magnetic flux above $10^{21}$ Mx in a way similar to Braun (2019), we indeed
find that this choice allows us to capture flows as far as $12^{\circ}$ away
from the active region centers. The contours capture not only the active-
region flows, but also the flows around the diffuse flux that is far less
concentrated, as seen in Fig. 1.

Figure 1: Contours around active-region flows (thick black lines) during solar
minimum (panel a) and during the rising phase (panel b), superimposed on
SDO/HMI magnetograms (line-of-sight magnetic field). The thin black lines
represent the contours of the magnetic mask discussed in Section 4.2. We only
consider pixels whose great-circle distance from disk center is less than
$50\degr$, and hence the contours do not encompass the magnetic features near
the limb.

####  2.2.2 Background flows

We consider the flows outside the contours as the background flows which
consist of the potential systematics. We generate the background flow maps by
computing monthly averages of the flows using only the pixels outside the
contours.

The monthly background flow maps can contain little to no data at mid-
latitudes during the solar maximum. Therefore, we smooth in time with a
Gaussian of FWHM of one year. This width is large enough that it smooths out
the regions containing no data for several months, while the background is not
expected to vary significantly on that time scale (see Appendix B). We further
smooth in space with a 2D Gaussian of FWHM of ten pixels, corresponding to an
angular distance of $6\degr$ at disk center.

These smoothed background flows exhibit a modulation on the time scale of a
solar cycle. This modulation is consistent with what is described by
Gottschling et al. (2021), although we compute the background flows in a
different way. They are studied in more detail in Appendix B.

###  2.3 Active-region flows

The flows inside the contours are the superposition of the background flows
and the active-region flows. We therefore subtract the smoothed background
flows from the 30-minute cadence flow maps; the remaining flows are the
active-region flows. The residual flows outside of the contours are not
discarded after the subtraction of the background flows. They are expected to
be random noise from convection.

We track these active-region flow maps at the Carrington rate on a daily basis
using noon as a reference, and remap them using the Plate Carree projection
into the heliographic coordinate system. The spatial resolution is
$0.6^{\circ}$ per pixel. Following the formulas in the appendix of Löptien et
al. (2017), we convert the velocities from pixels per second to meters per
second. Then we average the flows daily. The colatitudinal and longitudinal
components of these daily-averaged flows are denoted respectively by
$v_{\theta}^{\mathrm{AR}}(\theta,\phi,t)$ and
$v_{\phi}^{\mathrm{AR}}(\theta,\phi,t)$.

##  3 Temporal variation of active-region flows

Figure 2: Longitudinal average of $v_{\theta}^{\mathrm{AR}}$ (panel a;
positive values are southward) and $v_{\phi}^{\mathrm{AR}}$ (panel b; positive
values are prograde), as a function of time and latitude. They were further
smoothed in latitude with a Gaussian of FWHM of $3.6^{\circ}$ and in time with
a Gaussian of FWHM of 6 months. The black lines show the mean active
latitudes. The ticks on the horizontal axis indicate the beginning of each
year. The shaded areas indicate the times when the edge effects become visible
due to the one-year temporal smoothing of the background flows. Figure 3:
Longitudinally- and yearly-averaged $v_{\theta}^{\mathrm{AR}}$ and
$v_{\phi}^{\mathrm{AR}}$ (black lines), binned every $3.6\degr$ in latitude.
The standard error of the mean in each binning interval is about $1$ m/s for
$v_{\theta}^{\mathrm{AR}}$ and $0.4$ m/s for $v_{\phi}^{\mathrm{AR}}$, and is
not shown here. The magenta filling indicates the distribution of $|B_{r}|$
(normalized to the same arbitrary constant in all time periods), which was
averaged and binned in the same way as the flows.

To examine the temporal and latitudinal evolution of the active-region flows,
we average $v_{\theta}^{\mathrm{AR}}$ and $v_{\phi}^{\mathrm{AR}}$ in
longitude within $\pm 15^{\circ}$ of the central meridian. Figure 2 shows the
longitudinal averages of $v_{\theta}^{\mathrm{AR}}$ and
$v_{\phi}^{\mathrm{AR}}$ as a function of time and latitude (denoted as
$\lambda=90\degr-\theta$). Figure 3 shows the yearly-averaged flows as a
function of latitude, together with the magnetic activity.

Figure 2a shows an inflow pattern converging toward the mean latitude of
activity in both hemispheres. The inflows are stronger and cover a greater
latitudinal range on the poleward side of the active latitudes. We note that,
in Braun (2019) and Gottschling et al. (2021), the inflows are mostly
symmetric with respect to the center of the active regions. One possible
reason is that they use an ensemble averaging over active regions whose
centers have been aligned with each other, while we use a longitudinal long-
term averaging that makes the inflows on the equatorward side of the northern
and southern active latitudes partly cancel each other out. In the left column
of Fig. 3, the amplitude of the yearly-averaged $v_{\theta}^{\mathrm{AR}}$
reaches extrema of over $7$ m/s in the north in 2011 – 2012 and over $6$ m/s
in the south in 2014 – 2015. Likewise, the solar magnetic activity peaks in
2011 in the northern hemisphere and in 2014 in the southern hemisphere.
Therefore the amplitude of the inflows is clearly correlated with the strength
of the solar activity. The amplitude is consistent with that from, e.g., Gizon
(2003), González Hernández et al. (2008), Komm et al. (2020), who isolated
regions of magnetic activity and their surroundings in a way similar to what
we do here. We note that Spruit (2003) interpreted the inflows as a
consequence of the enhanced cooling in the magnetic regions and predicted
inflows of $\sim$6 m/s toward the activity belts. In the bottom two panels,
$v_{\theta}^{\mathrm{AR}}$ shows no active-region flows but rather random
oscillations that have the same magnitude at all latitudes, as expected toward
the end of the solar cycle.

Figure 2b shows that the toroidal component consists of a retrograde flow
(with respect to the background flow) on the poleward side of the active
latitudes, and a generally prograde flow on the equatorward side. This pattern
may be consistent with the results of Komm et al. (2020). It is known that the
torsional oscillations, which are the time-varying part of the solar rotation,
exhibit a shear flow around the active latitudes, with a faster-rotating band
on the equatorward side and a slower-rotating one on the poleward side (e.g.,
Howard & Labonte 1980). Therefore, the toroidal component of the active-region
flows seen here contributes to the torsional oscillations to some extent.

The right column of Fig. 3 shows that the amplitude of the toroidal component
is also correlated with the amplitude of the solar magnetic activity. The
retrograde flow on the poleward side of the active latitudes reaches over $8$
m/s in the north in 2011 – 2012 and over $6$ m/s in the south in 2013 – 2014.
We note that Braun (2019) and Gottschling et al. (2021) reported the presence
of a retrograde flow surrounding active regions preferentially on the poleward
side. In addition, there is a weaker prograde flow on the equatorward side of
the activity belts. The presence of this signal might be consistent with the
idea of a cyclonic circulation around active regions that derives from the
model of Spruit (2003) and that was also described by Hindman et al. (2009).

##  4 Forward-modeled helioseismic travel-time perturbations

###  4.1 Computation of travel-time perturbations

Figure 4: Forward travel-time perturbations $-\tau^{\mathrm{AR}}$, as a
function of latitude and separation distance. They have been averaged in
longitude around the central meridian and in time from January 2011 to
December 2014.

Figure 5: Panel a: Longitudinally-averaged active-region flows
$v_{\theta}^{\mathrm{AR}}$ as plotted in Fig. 2, shown again here for
convenience of comparison. Panel b: Longitudinally-averaged forward travel-
time perturbations $-\tau^{\mathrm{AR}}$, further averaged over the separation
distances $6\degr$ – $18\degr$. Panel c: same as panel b, but for the travel-
time perturbations with the magnetic mask applied,
$-\tau_{\mathrm{m}}^{\mathrm{AR}}$ (discussed in Section 4.2). We applied the
same smoothing as in Fig. 2. The black lines show the mean active latitudes.
Panels d–f: Averages of the top panels over the latitudes
$\lambda=20^{\circ}-35^{\circ}$ in each hemisphere. The typical errors in the
flows over the latitudinal range of interest are about $0.6$ m/s. For the
travel-time perturbations, the errors are in the range of $0.002$ – $0.008$ s,
with the largest error for $\tau_{\mathrm{m}}^{\mathrm{AR}}$ during the peaks
of activity. The errors are not shown here. In all panels, the shaded areas
indicate the times when the edge effects become visible due to the one-year
temporal smoothing of the background flows.

We compute the forward travel-time perturbations associated with the active-
region flows to compare with the solar-cycle variations of the helioseismic
measurements of the meridional circulation. To this end, we first construct a
3D flow profile associated with the active-region flows at the surface. We
neglect the contribution from the vertical flow to the travel times. For the
colatitudinal and longitudinal components $u_{\theta}$ and $u_{\phi}$, we
choose a constant profile with depth such that

$\displaystyle
u_{\theta}(\boldsymbol{r},t)=v_{\theta}^{\mathrm{AR}}(\theta,\phi,t),$ (1)
$\displaystyle
u_{\phi}(\boldsymbol{r},t)=v_{\phi}^{\mathrm{AR}}(\theta,\phi,t),$ (2)

where $\boldsymbol{r}=(r,\theta,\phi)$ with $r$ being the distance to the
solar center. Although it would be more realistic to choose a flow profile
where the amplitude decreases with depth, choosing a constant profile allows
us to estimate an upper bound for the amplitude of the travel-time
perturbations due to active-region flows.

In order to directly compare with Gizon et al. (2020), we use an arc-to-arc
geometry, with the travel-time perturbations being computed between pairs of
points placed on two opposite arcs. The paired arcs, each subtending an angle
of $30^{\circ}$, are aligned in the north-south direction as described by
Liang et al. (2017). Let’s denote $\psi$ the angle between a meridian and the
ray path connecting the paired points on the arcs, $(\theta_{0},\phi_{0})$ the
mid-point between the paired arcs, and $\Delta$ the separation distance
between the paired points. The travel-time perturbation is thus defined as
(e.g., Gizon et al. 2017; Fournier et al. 2018)

$\begin{split}\tau(\theta_{0},\phi_{0},\Delta,\psi,t)=\int_{\odot}&\Bigl{(}{K}_{\theta}({\boldsymbol{r}};\theta_{0},\phi_{0},\Delta,\psi)\;u_{\theta}(\boldsymbol{r},t)\\\
&+{K}_{\phi}({\boldsymbol{r}};\theta_{0},\phi_{0},\Delta,\psi)\;u_{\phi}(\boldsymbol{r},t)\Bigl{)}\;{\rm
d}\boldsymbol{r},\end{split}$ (3)

where the integral is taken over the whole volume of the Sun, and $K_{\theta}$
and $K_{\phi}$ are the colatitudinal and longitudinal components of the
sensitivity kernel (see Appendix C). Here, by convention, a northward flow
perturbation corresponds to a positive travel-time perturbation for short
separation distances.

The forward travel-time perturbation is averaged in the same way as in Gizon
et al. (2020); that is,

$\tau^{\mathrm{AR}}(\theta_{0},\Delta,t)=\left(N_{\phi_{0}}N_{\psi}\right)^{-1}\times\sum\limits_{\phi_{0},\psi}\;\tau(\theta_{0},\phi_{0},\Delta,\psi,t),$
(4)

where the sum over $\phi_{0}$ is taken within $\pm 15\degr$, the sum over
$\psi$ is taken within $\pm 15\degr$, $N_{\phi_{0}}$ is the number of points
used in the longitudinal average, and $N_{\psi}$ is the number of points on an
arc for each distance.

Figure 4 shows the forward travel-time perturbations as a function of latitude
and separation distance after we averaged over the active period from January
2011 to December 2014. Even though we chose a constant flow profile with
depth, the forward travel-time perturbations could still change sign with
increasing separation distance, for example, when the separation distance is
larger than the spatial scale of active region flows. Our $\tau^{\mathrm{AR}}$
decreases with separation distance but does not change sign up to at least
$\Delta=18^{\circ}$, which is similar to the modeled travel-time perturbations
from Liang et al. (2018).

We show the longitudinally-averaged $v_{\theta}^{\mathrm{AR}}$ (as plotted in
Fig. 2) in Fig. 5a, and the corresponding $\tau^{\mathrm{AR}}$ averaged over
the separation distances $6^{\circ}$ – $18^{\circ}$ in Fig. 5b. We can clearly
see the inflow pattern converging toward the active latitudes in Fig. 5b, with
an amplitude that varies throughout the solar cycle as in panel Fig. 5a. When
averaged over the latitudes $20^{\circ}-35^{\circ}$ in each hemisphere, the
amplitude of the inflows reaches $5$ m/s in 2011 in the north and $-5$ m/s in
2014 in the south (Fig. 5d). Similarly, the amplitude of $\tau^{\mathrm{AR}}$
reaches extrema of $0.3$ s during these active periods (Fig. 5e). For
comparison, we also computed $\tau^{\mathrm{AR}}$ using the radial flow
profile from the shallow model LC2 described in Liang et al. (2018); in this
case, the amplitude of $\tau^{\mathrm{AR}}$ reaches extrema of $0.15$ s.

###  4.2 Effect of masking travel-time perturbations inside magnetic regions

Gizon et al. (2020) excluded the travel-time perturbations measured within
magnetic regions from the averages, as those measurements introduce a
systematic error that resembles a divergent flow pattern (Liang & Chou 2015).
In order to compare with their results, we apply the same masking as they did.
The contours of the mask are shown as the thin black lines in Fig. 1. We
rewrite the forward travel-time perturbation averaging with a weighting
function $w$ that is equal to zero if the forward travel-time perturbation is
excluded and one elsewhere. Equation 4 then becomes

$\begin{split}\tau^{\mathrm{AR}}_{\mathrm{m}}(\theta_{0},\Delta,t)&=\left(\sum\limits_{\phi_{0},\psi}w(\theta_{0},\phi_{0},\Delta,\psi,t)\right)^{-1}\times\\\
&\sum\limits_{\phi_{0},\psi}\;w(\theta_{0},\phi_{0},\Delta,\psi,t)\;\tau(\theta_{0},\phi_{0},\Delta,\psi,t).\end{split}$
(5)

Figure 6: Longitudinally-averaged forward travel-time perturbations, further
averaged over the separation distances $6\degr$ – $18\degr$ and over the time
period from January 2011 to December 2014, as a function of latitude. They
have been binned every $3.6\degr$ in latitude. The typical standard error of
the mean, computed in each binning interval, is about $0.003$ s, and is not
shown here. The vertical gray shaded areas show the mean active latitudes.

The third column of Fig. 5 shows the forward travel-time perturbations after
the magnetic mask is applied, denoted by $\tau_{\mathrm{m}}^{\mathrm{AR}}$.
These perturbations are associated with the flows located outside the thin
lines and inside the thick lines in Fig. 1. The
$\tau_{\mathrm{m}}^{\mathrm{AR}}$ in Fig. 5c resembles the
$\tau^{\mathrm{AR}}$ in Fig. 5b but the amplitude of the inflow pattern is
strongly reduced during solar maximum. Figure 5f shows that, when averaged
over mid-latitudes, the amplitude reaches extrema in 2011 – 2012 in the north
and in 2014 – 2015 in the south; the extrema are three times smaller than that
of $\tau^{\mathrm{AR}}$ in Fig. 5e.

Figure 6 presents a comparison between $\tau^{\mathrm{AR}}$ and
$\tau_{\mathrm{m}}^{\mathrm{AR}}$, averaged over four active years from
January 2011 to December 2014. $\tau^{\mathrm{AR}}$ shows an inflow pattern
with an amplitude reaching extrema at a latitude of about $20\degr$ in both
hemispheres. With the masking, the inflow pattern in
$\tau_{\mathrm{m}}^{\mathrm{AR}}$ converging toward the active latitudes is
nearly gone. The difference
$\tau_{\mathrm{m}}^{\mathrm{AR}}-\tau^{\mathrm{AR}}$ presents a clear inflow
pattern toward the active latitudes, comparable to $\tau^{\mathrm{AR}}$.
Therefore, the masking of the magnetic regions removes most of the inflow
pattern from the averaged travel-time perturbations.

We note that the masking of magnetic pixels implemented in Gizon et al. (2020)
checks an area of $4\times 4$ pixels around either of the paired foot points
in the arc-to-arc geometry (Liang et al. 2017); if the field strength is
greater than a threshold in this area, the paired points are excluded from the
averages. We found that if only the nearest pixel to the foot point (instead
of a $4\times 4$-pixel area) is checked, the reduction of the inflow pattern
in the forward travel-time perturbations is not as strong as the
aforementioned results.

##  5 Extension to May 1996 – April 2019

Figure 7: Modeled inflows $\widetilde{v}_{\theta}^{\mathrm{\,AR}}$ (panel a),
and travel-time perturbations without magnetic mask
$-\widetilde{\tau}^{\mathrm{\,AR}}$ (panel b) and with magnetic mask
$-\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ (panel c) from May 1996 to
April 2019. The conversion constant for producing
$\widetilde{\tau}^{\mathrm{\,AR}}$ and
$\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ is derived from the HMI data
averaged over the separation distances $6\degr$ – $18\degr$. The same
smoothing as in Fig. 2 is applied. The black lines show the mean active
latitudes. Figure 8: Panels a–c: Averages of
$\widetilde{v}_{\theta}^{\mathrm{\,AR}}$, $-\widetilde{\tau}^{\mathrm{\,AR}}$,
and $-\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ from Figs. 7a–c over the
latitudes $\lambda=20^{\circ}-35^{\circ}$ in each hemisphere (thin lines). The
errors on the models, estimated from the misfits, are about $0.6$ m/s for
panel a, $0.04$ s for panel b, and $0.03$ s for panel c; they are not shown
here. For comparison, the corresponding $v_{\theta}^{\mathrm{AR}}$,
$-\tau^{\mathrm{AR}}$, and $-\tau_{\mathrm{m}}^{\mathrm{AR}}$ from Figs. 5d–f
are overplotted in thick lines; the data in the first and last six months are
not shown because the edge effects become visible due to the one-year temporal
smoothing of the background flows. Panel d: 13-month running mean of the
hemispheric monthly sunspot numbers. The red (respectively blue) filling
indicates an excess of sunspot numbers in the northern (respectively southern)
hemisphere.

We aim to extend the analysis to cover the time period 1996 – 2019 in order to
compare the results with Gizon et al. (2020). Although the LCT data are only
available after May 2010 as they were computed from full-resolution HMI
intensity images, we could use the magnetic field as a proxy for the active-
region inflows from May 1996 to April 2010. The correlation between the
magnetic field and the inflows is visible in Fig. 3: the greater the
latitudinal gradient of the magnetic field strength, the greater the inflows.
Several models already exist in the dynamo literature to generate artificial
active-region inflows based on the observations of the magnetic field (e.g.,
De Rosa & Schrijver 2006; Cameron & Schüssler 2010, 2012).

We minimize the cost function

$\sum_{\theta,t}\left|\langle v_{\theta}^{AR}(\theta,\phi,t)\rangle-
c_{0}\frac{\partial\langle|B_{r}(\theta,\phi,t)|\rangle}{\partial\theta}\right|^{2}$
(6)

to determine the constant of proportionality $c_{0}$. Here
$\langle\,\cdot\,\rangle$ denotes the averaging in longitude within $\pm
15\degr$ of the central meridian and the smoothing in latitude and time in the
way described in Fig. 2. We use only the data from January 2011 to December
2014 and in the latitudinal range $20\degr$ – $35\degr$, where the inflows are
the strongest. We find $c_{0}\simeq 0.085$ m s-1 G-1.

We compute the modeled inflows $\widetilde{v}_{\theta}^{\mathrm{\,AR}}$ using
the SOHO/MDI magnetograms (Scherrer et al. 1995) from May 1996 to April 2010
and the HMI magnetograms from May 2010 to April 2019. The line-of-sight
magnetic field inferred from HMI data is smaller than that from MDI data (Liu
et al. 2012); the scaling factor depends on the location on the disk and on
the field strength. We determine this scaling factor using $B_{r}$ at
latitudes $20\degr$ – $35\degr$ in each hemisphere and longitudes within $\pm
15\degr$ of the central meridian, during the time period from May 2010 to
April 2011 when both data sets are available. We find the scaling factor is
about $0.74$ which is then applied to the MDI data.

Finally, based on the similarity between the inflows and the travel-time
perturbations shown in Fig. 5, we use a conversion constant to convert from
the modeled flows $\widetilde{v}_{\theta}^{\mathrm{\,AR}}$ to modeled travel-
time perturbations from May 1996 to April 2019. For an average over the
separation distances $6^{\circ}$ – $18^{\circ}$, we find the conversion
constant to be $-0.056$ s2/m for the case without the magnetic mask
$\widetilde{\tau}^{\mathrm{\,AR}}$ and $-0.023$ s2/m for the case with the
magnetic mask $\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$, where the tilde
indicates the modeled quantities. We present the scatter plots of the data
used for the modeling in Appendix D.

Figure 7 shows the modeled inflows and the travel-time perturbations without
and with magnetic mask as a function of time and latitude, over Cycle 23 (1996
– 2008) and Cycle 24 (2008 – 2019). All the panels exhibit the inflow pattern
with the flows generally converging toward the mean active latitudes, which is
similar to that in Fig. 5. The amplitude of the models is comparable to that
of the observations for both the inflows and the travel-time perturbations on
the poleward side of the activity belts. However, on the equatorward side, the
amplitude of the models is greater because we determined the proportionality
constant $c_{0}$ using only the latitudes $20\degr$ – $35\degr$ in each
hemisphere.

Figure 8 shows how the models compare with the observations when averaged over
the latitudes $20^{\circ}$ – $35^{\circ}$ in each hemisphere.
$\widetilde{v}_{\theta}^{\mathrm{\,AR}}$ (Fig. 8a) and
$\widetilde{\tau}^{\mathrm{\,AR}}$ (Fig. 8b) generally match
$v_{\theta}^{\mathrm{AR}}$ and $\tau^{\mathrm{AR}}$ well in both hemispheres,
especially during the peaks of magnetic activity. When the magnetic mask is
applied (Fig. 8c), the model $\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$
retrieves the correct order of magnitude but does not fully capture
$\tau_{\mathrm{m}}^{\mathrm{AR}}$. Furthermore, the amplitude of the models is
greater during the Cycle 23 solar maximum than that during the Cycle 24 solar
maximum, which is consistent with the sunspot number presented in Fig. 8d.

##  6 Comparison with helioseismic measurements

We compare in Fig. 9 our forward travel-time perturbations with the
measurements from Gizon et al. (2020) spanning Cycles 23 and 24. Gizon et al.
(2020) used the data from MDI and GONG (Harvey et al. 1996) for the periods
from May 1996 to April 2003 and from May 2003 to April 2019, respectively; the
travel times associated with the magnetic pixels were excluded in their
measurements. For our travel-time perturbations, we used
$\widetilde{\tau}^{\mathrm{\,AR}}$ and
$\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ for the period from May 1996
to April 2010, and $\tau^{\mathrm{AR}}$ and $\tau_{\mathrm{m}}^{\mathrm{AR}}$
for the period from May 2010 to April 2019.

The measurements of Gizon et al. (2020) exhibit a modulation in amplitude over
the solar cycles that is correlated with the magnetic activity. The
$\tau^{\mathrm{AR}}$ and $\widetilde{\tau}^{\mathrm{\,AR}}$ also exhibit this
modulation. We remind the reader that we used a constant profile with depth
for the active-region flow model to place an upper limit; that is, the
magnitude of $\tau^{\mathrm{AR}}$ and $\widetilde{\tau}^{\mathrm{\,AR}}$ is
overestimated. When the masking is taken into account, the magnitude of
$\tau_{\mathrm{m}}^{\mathrm{AR}}$ and
$\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ is much smaller than that of
the measurements. These results suggest that the true travel-time
perturbations caused by the inflows do not fully explain the solar-cycle
variations in the travel-time measurements.

Figure 9:  Comparison between the north-south travel-time perturbations as
measured by Gizon et al. (2020) (thick solid lines) and the forward travel-
time perturbations without (thin solid line) and with (thick dashed line) the
magnetic mask. The measurements were averaged over the latitudes
$\lambda=20\degr$ – $35\degr$ in each hemisphere, over the separation
distances $\Delta=6\degr$ – $18\degr$, and over consecutive 4-year time
intervals. The average over the quiet-Sun period January 2007 to April 2009
was subtracted from the measurements, and the sign was reversed for the
purpose of comparison. The corresponding $-\widetilde{\tau}^{\mathrm{\,AR}}$
and $-\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ are used from May 1996 to
April 2010, and $-\tau^{\mathrm{AR}}$ and $-\tau_{\mathrm{m}}^{\mathrm{AR}}$
from May 2010 to April 2019. For the measurements, the error bars represent
the standard error of the mean computed over the latitudinal range of
interest. For the forward travel-time perturbations, the standard error of the
mean is in the range of $0.001$ – $0.003$ s, and is not shown here. For the
modeled travel-time perturbations, the errors, not shown here, are about
$0.013$ s for $-\widetilde{\tau}^{\mathrm{\,AR}}$ and $0.009$ s for
$-\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$.

##  7 Summary and discussion

We used LCT flow maps over nine years during Cycle 24. We removed large-scale
background flows and extracted the flows associated with active regions. We
averaged these flows in longitude to study their latitudinal profile and time
evolution. The latitudinal flow exhibits an inflow pattern converging toward
the active latitudes with a maximal peak-to-peak amplitude of $10$ m/s (yearly
average) during solar maximum. This amplitude is consistent with that found by
Gizon (2003), González Hernández et al. (2008) and Komm et al. (2020). The
longitudinal flow exhibits a pattern around the activity belts resembling that
of the torsional oscillations. This pattern may be consistent with that
observed by Komm et al. (2020). The maximum peak-to-peak amplitude is about
$10$ m/s. The amplitude and the structure of both flow components are
correlated with the strength and the distribution of the magnetic activity.

We computed the corresponding forward helioseismic travel-time perturbations
using 3D sensitivity kernels and using an arc-to-arc geometry in the north-
south direction. We assumed constant inflows with depth to place an upper
limit on the contribution of inflows to helioseismic measurements of the
meridional circulation. For separation distances $\Delta=6^{\circ}$ –
$18^{\circ}$, the extrema of the averaged forward travel-time perturbations
are $\pm 0.3$ s at mid-latitudes during the peaks of solar activity. We also
averaged the forward travel-time perturbations with the masking of magnetic
regions as in Gizon et al. (2020), and found that the masking significantly
reduced the amplitude of the travel-time perturbations, leading to extrema of
about $\pm 0.1$ s during solar maximum.

We extended the active-region flows and the forward travel-time perturbations
to cover two solar cycles, from May 1996 to April 2019, using a model based on
the latitudinal gradient of the magnetic field strength. We assumed that the
travel-time perturbations are roughly proportional to the flows. This simple
model reproduces the inflow patterns throughout the solar cycle and allows us
to compare our results with the measurements of Gizon et al. (2020) over two
solar cycles. We found that the near-surface active-region flows do not
explain in full the solar-cycle variations seen in the measurements of the
meridional circulation.

We note that the background flows, defined as the flows that are far from
active regions, are expected to represent the systematics, but they might also
contain global-scale flows in the quiet regions, if any. We see that they
exhibit a large-scale modulation throughout the solar cycle that seems
correlated with the magnetic activity (Appendix B). Since the LCT data were
filtered to remove only the periods of 24 hr, 1 yr, and the mean of the time
series (Löptien et al. 2017; Gottschling et al. 2021), long-term variations of
the global-scale meridional flow could still remain in the background flows.
However, we cannot exclude that there exist systematics that vary with the
solar cycle. Separating the true global-scale flows from the systematics in
the LCT data is beyond the scope of this paper.

The fact that the surface active-region flows do not fully account for the
observed solar-cycle variations of the meridional flow implies that there may
be time-varying flows far from activity. Other similar studies in which the
active-region flows and quiet-region flows are separated seemed to confirm
this. González Hernández et al. (2008) used ring-diagram analysis to infer the
subsurface meridional flow, and found that the inflows persist even after the
flows surrounding the active regions are excluded. They attributed it to the
fact that their masking may not remove the weaker magnetic regions and the
diffuse field, but they also did not exclude the possibility that inflows may
exist in quiet regions. Similarly, Komm et al. (2020) studied the time
variations of subsurface flows for quiet and active regions, separately, over
the past two solar cycles and found a solar-cycle modulation in the quiet-Sun
flows.

Hathaway & Rightmire (2010) found that the meridional flow is stronger during
solar minimum, although they did not disentangle the active-region flows from
the global circulation. Lin & Chou (2018) assumed that the meridional flows
are simply a linear combination of the local inflows and global-scale
meridional flows and found a similar result. A number of other studies, in
which the authors subtracted a time-averaged meridional flow profile, observed
residuals during solar minima (e.g., Hathaway & Rightmire 2011; Komm et al.
2015; Getling et al. 2021); however, the patterns of these residuals differ
from one another, depending on the time periods used to compute the reference.
To avoid this dependence, González Hernández et al. (2010) subtracted a low-
order polynomial fit and also found meridional flow residuals during solar
minimum. This might explain why the active-region flows only account for a
fraction of the solar-cycle variations of the meridional flow.

We note finally that other phenomena related to the surface magnetic activity
can contribute to the temporal variation of the travel-time perturbations
measured by Gizon et al. (2020). In particular, the Woodard effect (Woodard
1997) can add systematics via the localized absorption of acoustic waves by
sunspots. This effect has not been taken into account in our study.

###### Acknowledgements.

P.-L. P. is part of the International Max Planck Research School. D. F. and L.
G. acknowledge funding from the ERC Synergy Grant WHOLE SUN #810218. We thank
B. Löptien for kindly accepting to provide the LCT data series, and N.
Gottschling for fruitful discussions. The HMI data are courtesy of NASA/SDO
and the HMI Science Team. SOHO is a project of international cooperation
between ESA and NASA. The computational resources were provided by the German
Data Center for SDO through grant 50OL1701 from the German Aerospace Center
(DLR). The sunspot numbers are from WDC-SILSO, Royal Observatory of Belgium,
Brussels.

## References

  * Braun (2019) Braun, D. C. 2019, ApJ, 873, 94
  * Cameron & Schüssler (2010) Cameron, R. H. & Schüssler, M. 2010, ApJ, 720, 1030
  * Cameron & Schüssler (2012) Cameron, R. H. & Schüssler, M. 2012, A&A, 548, A57
  * De Rosa & Schrijver (2006) De Rosa, M. L. & Schrijver, C. J. 2006, in ESA Special Publication, Vol. 624, Proceedings of SOHO 18/GONG 2006/HELAS I, Beyond the spherical Sun, ed. K. Fletcher & M. Thompson, 12
  * Dikpati & Gilman (2006) Dikpati, M. & Gilman, P. A. 2006, ApJ, 649, 498
  * Dikpati & Gilman (2009) Dikpati, M. & Gilman, P. A. 2009, Flux-Transport Solar Dynamos, Vol. 32, 67
  * Duvall (1979) Duvall, T. L., J. 1979, Sol. Phys., 63, 3
  * Fisher & Welsch (2008) Fisher, G. H. & Welsch, B. T. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 383, Subsurface and Atmospheric Influences on Solar Activity, ed. R. Howe, R. W. Komm, K. S. Balasubramaniam, & G. J. D. Petrie, 373
  * Fournier et al. (2018) Fournier, D., Hanson, C. S., Gizon, L., & Barucq, H. 2018, A&A, 616, A156
  * Getling et al. (2021) Getling, A. V., Kosovichev, A. G., & Zhao, J. 2021, ApJ, 908, L50
  * Giles et al. (1997) Giles, P. M., Duvall, T. L., Scherrer, P. H., & Bogart, R. S. 1997, Nature, 390, 52
  * Gizon (2003) Gizon, L. 2003, PhD thesis, STANFORD UNIVERSITY
  * Gizon (2004) Gizon, L. 2004, Sol. Phys., 224, 217
  * Gizon et al. (2017) Gizon, L., Barucq, H., Duruflé, M., et al. 2017, A&A, 600, A35
  * Gizon et al. (2010) Gizon, L., Birch, A. C., & Spruit, H. C. 2010, ARA&A, 48, 289
  * Gizon et al. (2020) Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020, Science, 368, 1469
  * Gizon et al. (2001) Gizon, L., Duvall, T. L., J., & Larsen, R. M. 2001, in IAU Symposium, Vol. 203, Recent Insights into the Physics of the Sun and Heliosphere: Highlights from SOHO and Other Space Missions, ed. P. Brekke, B. Fleck, & J. B. Gurman, 189
  * González Hernández et al. (2010) González Hernández, I., Howe, R., Komm, R., & Hill, F. 2010, ApJ, 713, L16
  * González Hernández et al. (2008) González Hernández, I., Kholikov, S., Hill, F., Howe, R., & Komm, R. 2008, Sol. Phys., 252, 235
  * Gottschling et al. (2021) Gottschling, N., Schunker, H., Birch, A. C., Löptien, B., & Gizon, L. 2021, A&A, 652, A148
  * Haber et al. (2004) Haber, D. A., Hindman, B. W., Toomre, J., & Thompson, M. J. 2004, Sol. Phys., 220, 371
  * Harvey et al. (1996) Harvey, J. W., Hill, F., Hubbard, R. P., et al. 1996, Science, 272, 1284
  * Hathaway & Rightmire (2010) Hathaway, D. H. & Rightmire, L. 2010, Science, 327, 1350
  * Hathaway & Rightmire (2011) Hathaway, D. H. & Rightmire, L. 2011, ApJ, 729, 80
  * Hindman et al. (2003) Hindman, B., Gizon, L., Haber, D., Duval, Thomas, J., & Toomre, J. 2003, in ESA Special Publication, Vol. 517, GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. Sawaya-Lacoste, 299–302
  * Hindman et al. (2009) Hindman, B. W., Haber, D. A., & Toomre, J. 2009, ApJ, 698, 1749
  * Howard & Labonte (1980) Howard, R. & Labonte, B. J. 1980, ApJ, 239, L33
  * Komm et al. (2015) Komm, R., González Hernández, I., Howe, R., & Hill, F. 2015, Sol. Phys., 290, 3113
  * Komm et al. (2020) Komm, R., Howe, R., & Hill, F. 2020, Sol. Phys., 295, 47
  * Komm et al. (1993) Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys., 147, 207
  * Liang et al. (2017) Liang, Z.-C., Birch, A. C., Duvall, Thomas L., J., Gizon, L., & Schou, J. 2017, A&A, 601, A46
  * Liang & Chou (2015) Liang, Z.-C. & Chou, D.-Y. 2015, ApJ, 805, 165
  * Liang et al. (2018) Liang, Z.-C., Gizon, L., Birch, A. C., Duvall, T. L., & Rajaguru, S. P. 2018, A&A, 619, A99
  * Lin & Chou (2018) Lin, C.-H. & Chou, D.-Y. 2018, ApJ, 860, 48
  * Lisle & Toomre (2004) Lisle, J. & Toomre, J. 2004, in ESA Special Publication, Vol. 559, SOHO 14 Helio- and Asteroseismology: Towards a Golden Future, ed. D. Danesy, 556
  * Liu et al. (2012) Liu, Y., Hoeksema, J. T., Scherrer, P. H., et al. 2012, Sol. Phys., 279, 295
  * Löptien et al. (2017) Löptien, B., Birch, A. C., Duvall, T. L., et al. 2017, A&A, 606, A28
  * Löptien et al. (2016) Löptien, B., Birch, A. C., Duvall, T. L., Gizon, L., & Schou, J. 2016, A&A, 590, A130
  * Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd edn. (USA: Cambridge University Press)
  * Scherrer et al. (1995) Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229
  * Spruit (2003) Spruit, H. C. 2003, Sol. Phys., 213, 1
  * Wang et al. (1991) Wang, Y. M., Sheeley, N. R., J., & Nash, A. G. 1991, ApJ, 383, 431
  * Welsch et al. (2004) Welsch, B. T., Fisher, G. H., Abbett, W. P., & Regnier, S. 2004, ApJ, 610, 1148
  * Woodard (1997) Woodard, M. F. 1997, ApJ, 485, 890

##  Appendix A Contour determination

Figure 10: Ensemble-averaged active-region flows. $\Delta\phi$ and
$\Delta\lambda$ are the longitude and latitude with respect to the center of
the active region. The flows close to the center of each polarity have been
masked out before averaging because LCT flows in highly-magnetic regions are
not reliable (Löptien et al. 2017). The background shows the ensemble-averaged
line-of-sight magnetic field. For a given value of $b$ that determines a
contour for an individual active region, we draw the line delimiting the areas
included in the contours for at least $90\%$ of the active regions used in the
ensemble average; we try the following values of $b$: $4$ G, $3.5$ G, and $3$
G (dotted, solid, and dashed green lines, respectively).

We first compute the ensemble average of active regions in a way similar to
Braun (2019) and Gottschling et al. (2021). A brief description of the
procedure is as follows. We track and remap the $B_{r}$ maps (obtained in
Section 2.2.1) in the same way as the flow maps in Section 2.3. We compute the
daily averages of unsigned magnetic flux density maps from $B_{r}$. We smooth
these maps by a 2D Gaussian with a FWHM of $10^{\circ}$ (heliographic
degrees). We identify the peaks (pixels of value higher than that of any of
the eight neighboring pixels), and discard the ones that are within $20\degr$
of other peaks with stronger flux density. This procedure ensures that all the
selected peaks are clearly separated from each other. We compute the total
unsigned magnetic flux over a box, spanning $20\degr$ in longitude and
$10\degr$ in latitude around each peak. Only if the total unsigned magnetic
flux is above $10^{21}$ Mx do we retain the active regions associated with the
peaks. We compute the center of mass of the pixels with positive $B_{r}$ and
that of the pixels with negative $B_{r}$, both weighted by the smoothed and
daily-averaged unsigned flux density; the average of the two centers of mass
is defined as the center of an active region. To limit the impact of noise, we
do not consider pixels with $|B_{r}|$ less than $50$ G in the computation of
the center of mass and of the magnetic flux. We align the selected active
regions with respect to their center. If an active region is in the southern
hemisphere, we flip it in the north-south direction and we reverse the sign of
$B_{r}$ and $v^{\text{AR}}_{\theta}$, so that it can be averaged together with
the active regions in the northern hemisphere.

Figure 10 shows the ensemble-averaged active-region flows using the
$v^{\text{AR}}_{\theta}$ and $v^{\text{AR}}_{\phi}$ from Section 2.3. The
inflow pattern is clearly seen, with the convergence center located
preferentially in the trailing polarity. The pattern extends on average up to
$10^{\circ}$ away from the center of each polarity, in both the latitudinal
and the longitudinal directions. This is consistent with the results of, e.g.,
Löptien et al. (2017), Braun (2019) and Gottschling et al. (2021).

For each active region, we define the contour as the line along which the
smoothed $|B_{r}|$ (see Section 2.2.1) is equal to a magnetic threshold $b$.
Three values of $b$, $3.5$ G, $4$ G, and $4.5$ G, are tested. For each value,
we draw in Fig. 10 the line delimiting the area included in the contours for
at least $90\%$ of the active regions. All three areas cover up to at least
$10\degr$ away from the active regions, that is, they include most of the
active-region inflows.

We want to choose the contour that extends as far away as possible from the
active regions, to make sure that all the active-region flows are included.
However, with $b=3$ G, the background flows contain no data for more than a
year at mid-latitudes around 2011 in the north and around 2014 in the south;
in that case, smoothing the background flow maps in time with a Gaussian of
FWHM of one year is not possible. Increasing the FWHM might also not be
adequate because the time scale of variation of the background flows is on the
order of years (see Appendix B). On the other hand, using $b=4$ G results in
weaker inflows, because the background flows may contain part of the outer
edge of the inflows. As a result, we choose $b=3.5$ G as the best compromise
to consider as much active-region flows as possible and still have a reliable
background estimation.

##  Appendix B Temporal variation of the background flows

We track, remap and average the background flows in the same way as the
active-region flows in order to study their temporal variation. The results
are presented in Fig. 11.

Figure 11: Longitudinally-averaged background flows $v_{\theta}^{\mathrm{QS}}$
(panel a; positive values are southward) and $v_{\phi}^{\mathrm{QS}}$ (panel
b; positive values are prograde). We smoothed in latitude with a Gaussian of
FWHM of $3.6^{\circ}$ and in time with a Gaussian of FWHM of $6$ months.
Panels c–d: Averages of the top panels over the latitudes
$\lambda=20^{\circ}-35^{\circ}$ in each hemisphere. The typical standard
errors of the mean, computed over the latitudinal interval, are about $0.06$
m/s for $v_{\theta}^{\mathrm{QS}}$ and about $0.03$ m/s for
$v_{\phi}^{\mathrm{QS}}$, and are not shown here. In all panels, the shaded
areas indicate the times when the edge effects become visible due to the one-
year temporal smoothing of the flows done in Section 2.2.2.

A large-scale modulation of the meridional component is present, with a peak-
to-peak amplitude of about $12$ m/s after we average over mid-latitudes (Fig.
11c). This trend is roughly antisymmetric with respect to the equator, and
resembles an equatorward flow in the first half of the time series, then a
poleward flow in the second half toward solar minimum (Fig. 11a). This is
consistent with Gottschling et al. (2021). These flows contain systematics,
but they might also contain true flows. If this is the case, this trend may
indicate that the global-scale meridional flows become stronger when
approaching the solar minimum, as was previously reported by, e.g., Hathaway &
Rightmire (2010), González Hernández et al. (2010), Komm et al. (2015), Lin &
Chou (2018), and Getling et al. (2021).

The toroidal component also shows a modulation that is roughly symmetric with
respect to the equator (Fig. 11b). There exists a decreasing trend with time,
with an amplitude that varies between $-3$ and $+5$ m/s (Fig. 11d). We remind
the reader that the mean value over the whole time period had been removed
from the LCT data (see Section 2.1), so the variation is with respect to the
time average. The prograde pattern during the first half of the time series
migrates toward the equator along with the activity belts, which resembles the
behavior of torsional oscillations (e.g., Komm et al. 2020). This suggests
that this component of the background flows might also contain true flows in
the quiet Sun.

We note that, as mentioned in Section 2.2.2 and Appendix A, the background
flow maps can contain little to no data at mid-latitudes during active times.
These cut-outs of active regions in the background flow maps are filled in
using neighboring data with Gaussian smoothing. This smoothing procedure may
result in a slight underestimate of the background flows at mid-latitudes.

##  Appendix C 3D Born sensitivity kernels

The computation of the 3D Born sensitivity kernels is based on the approaches
from Gizon et al. (2017) and Fournier et al. (2018). The wave field
$\psi(\mbox{\boldmath$r$},\omega)$ is solution of a scalar wave equation in
the frequency domain

$-(\omega^{2}+2{\rm i}\omega\gamma)\psi-2{\rm
i}\omega\mbox{\boldmath$u$}\cdot\nabla\psi-c\nabla\cdot\left(\frac{1}{\rho}\nabla(\rho
c\psi)\right)=s(\mbox{\boldmath$r$},\omega),$ (7)

where $\rho$ is the density, $c$ the sound speed, $\gamma$ the attenuation and
$s$ a stochastic source term describing wave excitation. The wave field is
related to the 3D wave displacement $\boldsymbol{\xi}$ through
$\psi(\mbox{\boldmath$r$},\omega)=c(\mbox{\boldmath$r$})\nabla\cdot\boldsymbol{\xi}(\mbox{\boldmath$r$},\omega)$.

The flow kernels $\mathbf{K}=(K_{r},K_{\theta},K_{\phi})$ can be computed from
the knowledge of the Green’s function solution of Eq. (7) with a Dirac on the
right hand side:

$\displaystyle\mathbf{K}(\mbox{\boldmath$r$};\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2})=2{\rm
i}\rho(r)\int_{-\infty}^{\infty}\omega$ $\displaystyle
W^{\ast}(\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2},\omega)\Bigl{[}G(\mbox{\boldmath$r$}_{2};\mbox{\boldmath$r$},\omega)\nabla
C(\mbox{\boldmath$r$}_{1};\mbox{\boldmath$r$},\omega)$
$\displaystyle-G^{\ast}(\mbox{\boldmath$r$}_{1};\mbox{\boldmath$r$},\omega)\nabla
C^{\ast}(\mbox{\boldmath$r$}_{2};\mbox{\boldmath$r$},\omega)\Bigr{]}\;\textrm{d}\omega,$
(8)

where $\mbox{\boldmath$r$}_{1}=(R,\theta_{1},\phi_{1})$ and
$\mbox{\boldmath$r$}_{2}=(R,\theta_{2},\phi_{2})$ are two observation points,
$R$ is the observation radius, $C$ is the cross-covariance, $W$ is a weighting
function in order to relate the travel-time perturbation to changes in the
cross-covariance, and $\nabla$ is the gradient operator with respect to $r$.
As in Gizon et al. (2017), we assume energy equipartition so that the cross-
covariance is related to the imaginary part of the Green’s function. For the
sake of simplicity, we now drop the $\omega$ in the notation of the cross-
covariance and of the Green’s function. In a spherically symmetric background,
the Green’s function depends only on the angular distance between source and
receiver and can be obtained from its Legendre coefficients (Fournier et al.
2018):

$G(r,\theta,\phi;\mbox{\boldmath$r$}_{i})=\frac{1}{\sqrt{2\pi}}\sum_{\ell}G_{\ell}(r;R)P_{\ell}(\cos\gamma_{i}),$
(9)

where $i\in\\{1,2\\}$, $P_{\ell}$ is the Legendre polynomial of order $\ell$,

$\cos\gamma_{i}=\cos\theta\cos\theta_{i}+\sin\theta\sin\theta_{i}\cos(\phi-\phi_{i}),$
(10)

and $G_{\ell}(r;R)$ is solution of

$-(\omega^{2}+2{\rm
i}\omega\gamma)G_{\ell}-c\frac{d}{dr}\left(\frac{1}{\rho}\frac{d}{dr}(\rho
cG_{\ell})\right)+\frac{\ell(\ell+1)c^{2}}{r^{2}}G_{\ell}=\delta(r-R).$ (11)

We keep the values of $\ell$ up to $300$.

Inserting Eq. (9) into Eq. (8), the kernel for the radial flow is given by

$\displaystyle
K_{r}(\mbox{\boldmath$r$};\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2})=\frac{2\rho(r)}{\pi}\sum_{\ell,\ell^{\prime}}\Bigl{[}$
$\displaystyle-
f_{\ell\ell^{\prime}}^{r}(r)P_{\ell^{\prime}}(\cos\gamma_{1})P_{\ell}(\cos\gamma_{2})$
$\displaystyle+g_{\ell\ell^{\prime}}^{r}(r)P_{\ell}(\cos\gamma_{1})P_{\ell^{\prime}}(\cos\gamma_{2})\Bigr{]},$
(12)

where

$\displaystyle f_{\ell\ell^{\prime}}^{r}(r)$
$\displaystyle=\int_{0}^{\infty}\omega\textrm{Im}\left[W^{\ast}(\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2},\omega)G_{\ell}(r;R)\partial_{r}C_{\ell^{\prime}}(r;R)\right]{\rm
d}\omega,$ (13) $\displaystyle g_{\ell\ell^{\prime}}^{r}(r)$
$\displaystyle=\int_{0}^{\infty}\omega\textrm{Im}\left[W^{\ast}(\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2},\omega)G_{\ell}^{\ast}(r;R)\partial_{r}C_{\ell^{\prime}}^{\ast}(r;R)\right]{\rm
d}\omega.$ (14)

Similarly, denoting

$\displaystyle f_{\ell\ell^{\prime}}^{j}(r)$
$\displaystyle=\int_{0}^{\infty}\omega\textrm{Im}\left[W^{\ast}(\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2},\omega)G_{\ell}(r;R)C_{\ell^{\prime}}(r;R)\right]{\rm
d}\omega,$ (15) $\displaystyle g_{\ell\ell^{\prime}}^{j}(r)$
$\displaystyle=\int_{0}^{\infty}\omega\textrm{Im}\left[W^{\ast}(\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2},\omega)G_{\ell}^{\ast}(r;R)C_{\ell^{\prime}}^{\ast}(r;R)\right]{\rm
d}\omega,$ (16)

where $j\in\\{\theta,\phi\\}$, the kernels for the horizontal flow components
are given by

$\displaystyle
K_{j}(\mbox{\boldmath$r$};\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2})=\frac{2\rho(r)}{\pi\;r}\sum_{\ell,\ell^{\prime}}\Bigl{[}$
$\displaystyle-
f_{\ell\ell^{\prime}}^{j}(r)\alpha_{1}^{j}P^{\prime}_{\ell^{\prime}}(\cos\gamma_{1})P_{\ell}(\cos\gamma_{2})$
$\displaystyle+g_{\ell\ell^{\prime}}^{j}(r)\alpha_{2}^{j}P_{\ell}(\cos\gamma_{1})P^{\prime}_{\ell^{\prime}}(\cos\gamma_{2})\Bigr{]},$
(17)

where $P^{\prime}$ is the derivative of the Legendre polynomials and

$\displaystyle\alpha_{i}^{\theta}$
$\displaystyle=-\cos\theta_{i}\sin\theta-\sin\theta_{i}\cos\theta\cos(\phi-\phi_{i}),$
(18) $\displaystyle\alpha_{i}^{\phi}$
$\displaystyle=-\sin\theta_{i}\sin(\phi-\phi_{i})$ (19)

are respectively the derivative of $\cos\gamma_{i}$ with respect to $\theta$,
and the derivative with respect to $\phi$ and divided by $\sin\theta$.

Let’s consider $\mbox{\boldmath$r$}_{1}$ and $\mbox{\boldmath$r$}_{2}$ two
foot points in the arc-to-arc geometry described by Liang et al. (2017) and
used in this paper. Let’s denote $(\theta_{0},\phi_{0})$ the colatitude and
longitude of the midpoint, $\Delta$ the separation distance between the foot
points, and $\psi$ between a meridian and the ray path connecting the paired
points on the arcs. Then we rewrite the flow kernels as

$\mathbf{K}(\mbox{\boldmath$r$};\mbox{\boldmath$r$}_{1},\mbox{\boldmath$r$}_{2})=\mathbf{K}(\mbox{\boldmath$r$};\theta_{0},\phi_{0},\Delta,\psi).$
(20)

We use this notation in the main text in order to make the averaging of the
travel-time perturbations over the arcs more explicit.

##  Appendix D Flow and travel-time modeling

The linear fitting procedure in Section 5 takes into account the errors in
both coordinates (Press et al. 2007, Section 15.3). We present the results of
the fits in Fig. 12. For Fig. 12a,c,d, because the data used were smoothed in
time and latitude, we use the points distant from each other by $6$ months and
$3.6\degr$ in latitude so that they are independent. In all panels, the value
of the intercept is smaller than the error in the vertical coordinate
estimated from the misfit, so we neglected it.

There is a clear linear correlation between
$\partial\langle|B_{r}|\rangle/\partial\theta$ and $\langle
v_{\theta}^{\mathrm{AR}}\rangle$ (Fig. 12a), between the HMI magnetic field
and the MDI magnetic field (Fig. 12b), and between
$\widetilde{v}_{\theta}^{\mathrm{\,AR}}$ and $\tau^{\mathrm{AR}}$ (Fig. 12c).
In contrast, the correlation between $\widetilde{v}_{\theta}^{\mathrm{\,AR}}$
and $\tau_{\mathrm{m}}^{\mathrm{AR}}$ (Fig. 12d) is not as high in other
panels, implying that a linear model might yield a poor fit. This is because
we cannot mask magnetic pixels in the modeled flows, which are a function of
time and latitude but not longitude. Therefore we have to use the same flows
for the fit of $\tau_{\mathrm{m}}^{\mathrm{AR}}$ as for the fit of
$\tau^{\mathrm{AR}}$. The larger relative error in the slope (shaded area in
red) in Fig. 12d reflects the poor linear fit.

Figure 12: Scatter plots of the data used for the determination of the
proportionality constant $c_{0}$ (panel a), of the scaling factor for the
magnetograms (panel b), of the conversion constant for
$\widetilde{\tau}^{\mathrm{\,AR}}$ (panel c), and the conversion constant for
$\widetilde{\tau}_{\mathrm{m}}^{\mathrm{\,AR}}$ (panel d). The red lines
indicate the best fits. The red shaded areas represent the errors on the
slopes, obtained from the fits. The numbers in the top left corners are the
values of the Pearson correlation coefficient, of the slope, and of the error
on the slope. In all cases, the $p$ values of the correlation coefficients are
too small to be given here ($\ll 0.01$).