# Improved Measurements of the Sun’s Meridional Flow and Torsional Oscillation
from Correlation tracking on MDI & HMI magnetograms

Sushant S. Mahajan  Department of Physics & Astronomy, Georgia State
University, Atlanta, GA, USA.  David H. Hathaway  W.W. Hansen Experimental
Physics Laboratory, Stanford University, Stanford, CA, USA.  Andrés Muñoz-
Jaramillo  Southwest Research Institute, Boulder, CO, USA.  High Altitude
Observatory, Boulder, CO, USA  National Solar Observatory, Boulder, CO, USA
Petrus C. Martens  Department of Physics & Astronomy, Georgia State
University, Atlanta, GA, USA.

###### Abstract

The Sun’s axisymmetric flows, differential rotation and meridional flow,
govern the dynamics of the solar magnetic cycle and variety of methods are
used to measure these flows, each with its own strengths and weaknesses. Flow
measurements based on cross-correlating images of the surface magnetic field
have been made since the 1970s which require advanced numerical techniques
that are capable of detecting movements of less than the pixel size in images
of the Sun. We have identified several systematic errors in addition to the
center-to-limb effect that influence previous measurements of these flows and
propose numerical techniques that can minimize these errors by utilizing
measurements of displacements at several time-lags. Our analysis of line-of-
sight magnetograms from the Michelson Doppler Imager (MDI) on the ESA/NASA
Solar and Heliospheric Observatory (SOHO) and Helioseismic and Magnetic Imager
(HMI) on the NASA Solar Dynamics Observatory (SDO) shows long-term variations
in the meridional flow and differential rotation over two sunspot cycles from
1996 to 2020. These improved measurements can serve as vital inputs for solar
dynamo and surface flux transport simulations.

Correlation Tracking, Meridional Flow, Torsional Oscillation, Center-to-limb
effect

##  1 Introduction

The Sun’s axisymmetric flows, differential rotation (DR) and meridional flow
(MF) play key roles in virtually all models of the Sun’s magnetic dynamo
(Charbonneau, 2010). The differential rotation (variation in rotation rate
with both latitude and depth) stretches the radial and latitudinal components
of the magnetic field in the longitudinal direction - thereby increasing the
field strength and changing the field direction. The meridional flow
transports the radial and longitudinal components of the magnetic field in the
latitudinal direction - thereby building up the polar fields on the surface
(Sheeley et al., 1989; Jiang et al., 2010) and annihilating oppositely
directed longitudinal fields across the equator in the solar interior. Thus,
measurements of the meridional flow are important for the prediction of the
build up of polar fields which in turn are good predictors for the next solar
cycle (Muñoz-Jaramillo et al., 2013).

Measurements of these axisymmetic flows have been made with a wide variety of
methods. Carrington (1859) measured the positions of sunspots to determine the
latitudinal variation in rotation rate through the limited range of sunspot
latitudes. Since those earliest measurements using the only known tracer at
the time (sunspots), many more measurements have been made using a variety of
tracers and/or methods. Tracers include: sunspots, magnetic features (network
magnetic elements), Doppler features (supergranules), and intensity features
(granules). Measurement methods (beside feature tracking) include direct
Doppler (spectrographic) and helioseismic measurements. Beck (2000) has
compared and contrasted these various methods. Each has its own set of
advantages and disadvantages. This is particularly true for measurements of
the meridional flow, primarily due to the intrinsic low amplitude of this
global flow.

Sunspot tracking is limited in both time and space. Sunspots only appear in
the lower latitudes and at times don’t appear at all. Most other methods have
much better coverage in space (latitude) and time but present other problems.

Direct Doppler measurements of the meridional flow (Hathaway et al., 1996;
Ulrich, 2010) require an accurate determination of the convective blue shift
signal and its dependence on local magnetic fields. This signal is stronger
than the meridional flow signal and is similar in spatial characteristics -
making the meridional flow measurement via its direct Doppler signal very
difficult.

Feature tracking the Doppler features (supergranules) has an advantage in
having a strong signal near the limb - and hence high latitudes, but suffers
from systematic changes in the features due to line-of-sight effects. Hathaway
et al. (2006) showed that these projection effects make the supergranule
pattern appear to rotate faster than it actually does. The impact of line-of-
sight effects on meridional flow measurements has, to our knowledge, not been
investigated.

Helioseismology also has difficulty near the limb, but more importantly,
exhibits a systematic center-to-limb shift in the apparent acoustic field.
This effect manifests itself as a pseudo flow away from the disk center which
increases as a function of center to limb distance and interferes with
measurements of the meridional flow. Zhao et al. (2012) accounted for it by
measuring the rotation rate as a function of longitude from the central
meridian at the equator and subtracting it from meridional flow measurements
assuming that the effect is symmetric on the solar disk.

Feature tracking the photometric intensity features (granules) requires high
spatial and temporal resolution which have only been available recently. These
features are difficult to track near the limb and their short lifetimes make
accurate measurements of slow flows, such as the meridional flow, quite
difficult. This kind of tracking shows a “shrinking” Sun effect (an apparent
flow toward disk center) with an amplitude of $\sim 1000\ m\ s^{-1}$ at a
heliocentric angle of $60^{\circ}$. Löptien et al. (2016) investigated this
effect and showed that in large parts it originates from the apparent
asymmetry of granulation due to radiative transfer effects when observing at a
viewing angle. However, their model of the “shrinking Sun” effect was not
accurate enough to measure the much weaker meridional flow. They also pointed
out that the shrinking Sun effect was not exactly symmetric across the central
meridian implying an east-west dependence.

Feature tracking the weak (network) magnetic features has the advantages of
good coverage over the visible disk (but with the same difficulty near the
limb as with granules and acoustic waves) coupled with long lifetimes that are
advantageous for meridional flow measurements. Komm et al. (1993) used local
correlation tracking (LCT) on magnetograms separated by a day in time. By
averaging their daily measurements over two year intervals they acquired
enough accuracy to show a systematic variation in the meridional flow with
solar cycle phase. Hathaway & Rightmire (2010) used a variant of the Komm et
al. (1993) method on magnetograms obtained with SOHO/MDI (Scherrer et al.,
1995) at eight-hour intervals. The sheer number of observations allowed for
accurate measurements over individual 27-day solar rotations. Dikpati et al.
(2010) suggested that these measurements were compromised by a systematic
effect due to supergranule diffusion. They argued that this diffusion would
produce an apparent outflow from magnetic flux concentrations in active
regions and tested their argument using a purely diffusive model with a
diffuse magnetic field. Hathaway & Rightmire (2011) investigated this
possibility with a more realistic model - magnetic elements transported by
evolving supergranule flows - and found that this systematic effect was too
small to be measured. They further noted that the nature of the meridional
flow variations with solar cycle phase produced inflows toward active
latitudes, not outflows.

Until now, a systematic center-to-limb effect was not reported in tracking
network magnetic elements making it the most promising method for measuring
the meridional flow. Here we describe our measurements of the meridional flow
and torsional oscillation using a feature tracking method (LCT) on weak
magnetic features seen in solar magnetograms. This method involves mapping
data from full-disk magnetograms onto heliographic longitude-latitude maps and
then cross-correlation of blocks of data between two maps separated in time by
a few hours. The displacement at which two blocks separated in time show the
best correlation then provides a measure of the displacement vector. Our
original intent was to incorporate the daily magnetograms from NSO/Kitt Peak
to extend the meridional flow measurements to cover four solar cycles but we
discovered early on that several improvements to the existing LCT methods were
needed (as described in Section 3).

We, in this paper describe the discovery of a center-to-limb effect similar to
the systematic flow away from disk center seen with helioseismic methods or
the opposite of the “shrinking Sun” effect seen in ganule tracking which
compromises feature tracking measurements from magnetograms reported by Lamb
(2017) and Imada & Fujiyama (2018) and to some extent Hathaway & Rightmire
(2010).

Testing our LCT algorithm with synthetic magnetograms also revealed the need
to improve the interpolation algorithm used for projecting magnetograms onto a
heliographic grid (uniformly spaced in both latitude and longitude).
Histograms of the meridional flow velocity measurements indicated: 1) the need
to broaden the search for the correlation maximum to larger distances at high
latitudes so as not to truncate the normal distribution (and thereby influence
the value of the average), 2) the need to obtain multiple measurements in
longitude with larger search windows at higher latitudes so as to cover
similar areas on the surface of the Sun, and 3) the need to more accurately
determine the position of the correlation maximum with fractional pixel
accuracy.

In the following sections we describe the data, our displacement measurement
method, discovery of systematic errors and the improvements we made. We end
with a synopsis of our improved measurements of the meridional flow and
torsional oscillation profile as a function of time through solar cycles 23
and 24.

##  2 Data

###  2.1 MDI/HMI Magnetograms

We obtained full disk (1024$\times$1024 pixel) line-of-sight magnetograms from
the Michelson Doppler Imager (MDI) (Scherrer et al., 1995) with a cadence of
96 minutes from May 4, 1996 to Jan 11, 2011 and full disk (4096$\times$4096
pixel) line-of-sight magnetograms from the Helioseismic and Magnetic Imager
(HMI) (Scherrer et al., 2012) with a cadence of one hour from May 20, 2010 to
June 27, 2020. The MDI (HMI) magnetograms are obtained at a 60-second
(45-second) cadence and averaged (by the instrument teams) over 5 minutes (12
minutes) to reduce the 5-minute oscillation signal.

These magnetograms are created from polarization measurements along the
profiles of magnetically sensitive photospheric absorption lines (the Ni I
$676.8$ nm line for MDI and the Fe I $617.3$ nm line for HMI). We use the
line-of-sight magnetic measurements from both instruments and project them
onto grids uniformly spaced in longitude and latitude. The methods used in
projecting the data onto these grids were tested for accuracy as described in
section 3.1.

Figure 1: A synthetic solar magnetogram with 4000 positive polarity magnetic
features produced by the Analytical Magnetogram Creator (AMC) and sampled by
the CCD emulator. The magnetic field shown here is in $Gauss$.

###  2.2 Synthetic Magnetograms for Testing

In order to calibrate the performance of any tracking algorithm, a “ground-
truth” dataset with known movements of features is needed. For this purpose,
we developed an algorithm which allows us to create synthetic magnetograms of
the Sun by placing magnetic features at predetermined locations with
analytical precision.

Given the central latitude ($\lambda_{1}$) and longitude ($\phi_{1}$) of a
magnetic feature, its 2D Gaussian magnetic field distribution is modelled as
in Eqs. 1 & 2 below:

$B(\lambda,\phi)=Ae^{-\frac{d^{2}}{w^{2}}}\cos(\rho+\eta)$ (1)

$d=2R_{\odot}\sin^{-1}\sqrt{\sin^{2}\bigg{(}\frac{\lambda-\lambda_{1}}{2}\bigg{)}+\cos(\lambda_{1})\cos(\lambda)\sin^{2}\bigg{(}\frac{\phi-\phi_{1}}{2}\bigg{)}}$
(2)

where A is the amplitude i.e. the maximum magnetic field flux density within
the feature, $\phi$ is the longitude, $\rho$ is the heliocentric angle from
disk center, $\eta$ is the angle between the line-of-sight and the line
joining the observer to the center of the Sun, $w$ is the width of the
Gaussian and $d$ is the great circle distance between the location of the peak
of the feature ($\lambda_{1},\phi_{1}$) and an evaluation point
($\lambda,\phi$).

We call this algorithm the “Analytical Magnetogram Creator (AMC)” and its
output is an analytical function which describes the radial component of
magnetic flux density on the surface of the synthetic Sun. Note here that
being analytical functions, magnetograms created by AMC have infinite
resolution as they can be sampled at any desired resolution.

To create MDI/HMI-like magnetograms from the output of AMC, we created a
Charge-Coupled Device (CCD) emulator which emulates a camera placed at a
prescribed location with respect to the Sun, with a prescribed plate scale
(resolution), P angle (position angle of the solar rotation axis relative to
the “vertical”), B angle (tilt of the solar rotation axis toward or away from
the observer) and coordinates of the center of the solar disk. It numerically
integrates the magnetic flux density within the boundaries of every image
pixel by dividing each pixel into 10$\times$10 sub-pixels and evaluating the
line-of-sight magnetic flux density in the total area of the pixel which is
equivalent to the measured magnetic field.

Features placed on synthetic magnetograms of the Sun can be moved precisely by
a prescribed distance with insignificant numerical errors (verified
empirically) because the magnetic field values come from an exact analytical
formula. Such synthetic magnetograms of the Sun enable the accurate
calibration of our correlation tracking algorithm. One such synthetic
magnetogram generated by placing 4000 magnetic features randomly on the solar
disk is shown in Fig. 1.

##  3 Improvements to the Correlation Tracking Algorithm

###  3.1 Latitude-Longitude interpolation used for Heliographic projection

Figure 2: The performance of various interpolation algorithms used in
projection routines. The first two panels in the top row show the heliographic
projection of a well resolved 2D Gaussian magnetic feature analytically placed
on a synthetic magnetogram and sampled at different resolutions. The third
panel in the first row and all panels in the second row show the same magnetic
feature when different interpolation algorithms are used to reconstruct the
higher resolution projection of the feature by oversampling its low resolution
image. In the third and fourth rows, corresponding plots for a feature which
is barely resolved in the high-resolution image and unresolved in the low
resolution image are shown.

Most correlation tracking studies on solar data project the solar images onto
a two dimensional plane before the movements of features are tracked. We
choose the equirectangular Heliographic projection (uniformly spaced in
latitude and longitude) for our task following the prescription of Hathaway &
Rightmire (2010). It is extremely important to minimize uncertainty in
projected locations of features and distortion in their shape to maximize the
accuracy of correlation tracking. For this purpose, we have comprehensively
tested our Heliographic projection algorithm using the AMC and the CCD
emulator.

In Heliographic projection algorithms, it is common to interpolate solar
images using a bicubic or bicubic spline interpolant (for example in Liang et
al. (2018)) which operates on the uniformly spaced image pixel grid. However,
these algorithms incorrectly assume that the physical information in each
pixel is equally spaced in physical distance as well. To evaluate the accuracy
of these interpolation techniques, we performed the following experiments. In
two separate experiments two magnetic features, one well resolved (FWHM = 5.33
Mm) and one barely resolved (FWHM = 1.78 Mm), were placed at $60^{\circ}$
latitude and $60^{\circ}$ longitude using the AMC and each of them were then
sampled by the CCD emulator at two different resolutions: $4096\times
4096$(high) and $1024\times 1024$(low). Both bicubic and bicubic spline
interpolation were then used to oversample the low resolution magnetograms in
order to reconstruct their high resolution counterparts. We found that
interpolation over the image pixel grid distorts the oversampled magnetograms
as shown in Fig. 2 with the smaller feature suffering more distortion.
Correlation between the oversampled and the high-resolution magnetogram listed
in Table 1 was used as a measure of accuracy.

To minimize distortion in the projection algorithm, we adopted the biharmonic
spline interpolation algorithm of Sandwell (1987) and employed it on the non-
uniformly spaced ($\lambda$, $\phi\cos(\lambda)$) coordinates of the centers
of image pixels. This interpolation algorithm relies on data values in
neighboring pixels and their physical distance from the target location. Two
variants of this algorithm were used, one which utilizes data from a $3\times
3$ pixel area and another which uses data from a $4\times 4$ pixel area
surrounding a target location that is in between the second and the third
pixel in both dimensions. All interpolation methods magnified the barely
resolved feature (Fig. 2) more than the well resolved feature. Bicubic and
bicubic spline interpolation completely changed alignment of the barely
resolved feature while they only distorted the shape of the well resolved
feature. Biharmonic spline interpolation utilizing data from $4\times 4$
neighboring pixels outperformed all other methods as is suggested by the
higher correlation coefficient between the oversampled and the high-resolution
images listed in Table 1 and minimal visual distortion in Fig. 2. This
experiment was repeated at various locations on the solar disk with similar
results.

In our correlation tracking algorithm, each magnetogram is projected onto a
uniformly spaced latitude-longitude heliographic grid ($-90^{\circ}$ to
$+90^{\circ}$ in both directions) with number of grid points equivalent to the
dimension of the magnetogram itself. So MDI magnetograms are projected onto a
$1024\times 1024$ and HMI magnetograms on a $4096\times 4096$ heliographic
grid.

Table 1:  Correlation coefficient between the oversampled and high resolution image when different interpolation algorithms are used to oversample a low resolution image into a high resolution image of a feature on a synthetic magnetogram.  Interpolation method | Correlation (well-resolved feature) | Correlation (barely resolved feature)
---|---|---
Bicubic (pixels) | 0.9604 | 0.7881
Bicubic spline (pixels) | 0.9664 | 0.7802
Biharmonic spline 3 $\times$ 3 (lat-long) | 0.9763 | 0.8380
Biharmonic spline 4 $\times$ 4 (lat-long) | 0.9833 | 0.8463

###  3.2 Choosing block sizes and changing search area with latitude

Figure 3: Blocks selected for tracking from an HMI magnetogram (in $Gauss$) projected onto a uniform grid in latitude-longitude are marked in red. The horizontal extent of the projected magnetogram stretches from $-90^{\circ}$ on the left to $+90^{\circ}$ longitude on the right and the vertical extent from the south pole to north pole. The yellow line indicates the solar equator. Table 2: Search area at solar disk center as a function of time-lag is shown here chosen so that velocities upto $\pm 400$ $m/s$ can be detected. The time unit for MDI is $96$ minutes and for HMI is one hour. The longitudinal search area at the equator is equal to the values in this table and varies with a factor of the secant of the latitude.  No. of time units | HMI search area | MDI search area
---|---|---
1 | $\pm 0.132^{\circ}$ | $\pm 0.352^{\circ}$
2 | $\pm 0.308^{\circ}$ | $\pm 0.527^{\circ}$
4 | $\pm 0.659^{\circ}$ | $\pm 1.055^{\circ}$
8 | $\pm 1.362^{\circ}$ | $\pm 2.285^{\circ}$

Hathaway & Rightmire (2010) used thin longitudinal strips to measure
meridional flow in line-of-sight magnetograms. This approach, however, does
not provide us any information about longitudinal variation in the
measurements. We have chosen to split the projected magnetogram into
rectangular blocks in latitude-longitude such that each block is large enough
to contain at least one supergranule ($\approx 30Mm$). Each of these blocks
has a latitudinal extent of $4.44^{\circ}$ so that there are 40 blocks
covering all latitudes. To avoid tracking a block all the way to the limb of
the Sun, all blocks were chosen within $60^{\circ}$ longitude from the central
meridian and their longitudinal extent was $8.93^{\circ}$ at the equator. To
keep the physical area inside each block constant, their longitudinal extent
was varied as a function of latitude as shown in Fig. 3. Tracking such blocks
provides us less latitudinal resolution than thin longitudinal strips, but we
gain longitudinal information in flow patterns.

Each of the blocks from one magnetogram ($mag1$) was cross-correlated within a
search area around their expected location (estimated by the average
differential rotation) in another magnetogram ($mag2$) a certain time apart.
This process was done both forward (select a block in $mag1$ and search around
its expected location in $mag2$) and backwards (select a block in $mag2$ and
search around its expected location in $mag1$) in time to double the number of
measurements and minimize any systematic errors. This search area was varied
with time-lag between successive magnetograms and by a factor of secant of the
heliocentric latitude of block center (see Table 2) in order to keep the
search area constant in physical units. The search area was chosen so that it
will encompass the distribution of velocities beyond three standard
deviations. In post-processing, we drop velocity measurements beyond two
standard deviations for each block while calculating the average velocity
(averaged over all magnetogram pairs in each Carrington rotation). Averaging
without excluding these measurements results in uneven coverage of the tails
on either side of the distribution and leads to underestimation of the
velocities by biasing the measurements towards the center of the search area.

To track the network magnetic elements, any 9$\times$9 pixel area in projected
HMI or 3$\times$3 pixel area in projected MDI magnetograms with an average
magnetic field in excess of 150 $Gauss$ was masked.

###  3.3 Peak-locking

Figure 4: Left: A histogram of all measurements of meridional flow velocity
from Carrington rotation 2098 at 8 hr time-lag shows peak-locking at $18$
$m/s$ intervals. Moreover, the measurements in the original algorithm were
clipped due to a small search area. Right: Corresponding histogram of
meridional flow measurements from our improved algorithm mitigates peak-
locking and captures the full range of the velocities making both the
individual and average measurements more accurate.

Hathaway & Rightmire (2010) found the location of the maximum correlation by
finding the parabola through the peak of the correlation coefficient to
measure the displacement of each block within any time-lag. Such a parabolic
fit has an accuracy of $~{}0.1$ pixel shift. A histogram of all measurements
of meridional flow velocities at a time-lag of eight hours obtained using this
method is shown in Fig. 4. We found that a simple parabolic (or 2D Gaussian)
fit to the spatial distribution of correlation coefficient results in a
phenomenon known as peak-locking in the field of particle image velocimetry
(Huang et al., 1997; Westerweel, 1997; Chen & Katz, 2005). We get more
measurements at integer pixel shifts (equivalent to every $18$ $m/s$ at HMI
resolution for 8 hr time-lag) than at fractional pixel shifts. Half-pixel
shifts tend to be far less likely to be found. This phenomenon appears to be
known in the solar-physics community (Langfellner et al., 2015) but has not
been specifically addressed by name. While Langfellner et al. (2015) and B.T.
Welsch (private communication) address this by an iterative procedure which
translates the cross-correlation, we resolve this problem by taking the result
from the parabolic fit as an estimate of the displacement, and then
interpolating the reference block to shift it by this fractional pixel
estimate and iterating the procedure until we find a near zero shift. The
total displacement is then a sum of shifts through all iterations.
Implementing such an iterative procedure improved our histogram of meridional
flow measurements (see Fig. 4) and increased the accuracy of tracking five
times to $0.02$ pixel shift.

###  3.4 Discovery and Correction of the Center-to-limb effect

Tracking On Signed Magnetograms Tracking On Unsigned Magnetograms Tracking On
Synthetic Magnetograms

Figure 5: A comparison of tracking results from signed magnetograms (top row)
and unsigned magnetograms (middle row) for Carrington rotation 2098. The
bottom row shows tracking results from synthetic data created with a noise
level of 10 Gauss (FWHM). The first column in each row shows the residual
rotation rate measured at one hour time-lag obtained after removing the
average rotation rate measured at each latitude, while the second column shows
the residual rotation rate measured at the equator at several time-lags. The
third column in each row shows the residual meridional flow speed measured at
one hour time-lag obtained after removing the average meridional flow speed
measured at each latitude, while the fourth column shows the latitudinal
variation of meridional flow speed at several time-lags.

Flow measurements from correlation tracking performed on signed magnetic field
at several time-lags averaged over Carrington rotation 2098 are shown in the
top row of Fig. 5. These measurements indicate that the east limb of the Sun
is rotating $\sim 100$ m/s slower than the west limb of the Sun at one hour
time-lag. This difference in rotation rate between the east and west limb
decreases as a function of time-lag. The meridional flow measurements after
removal of the longitudinally averaged flow speed at each latitude do not show
any limb to limb variation in the longitudinal direction. The meridional flow
profiles themselves, however, vary significantly as a function of time-lag.
This pattern of variation with time-lag was also seen when the algorithm of
Hathaway & Rightmire (2010) was used for tracking. This combination of
variation in meridional flow and the east-west asymmetry in rotation rate with
time-lag is reminiscent of the center-to-limb variation reported by Zhao et
al. (2012). It appears like a nearly instantaneous displacement of all
features away from disk center which contributes less to velocity measurements
at longer time-lags.

A very similar variation in meridional flow (measured from one month of data)
was reported by Lamb (2017) by tracking individual magnetic features with
varying lifetimes: 2 hr, 10 hr and 20 hr. Lamb (2017) found faster flow speeds
for the shorter lived features. This variation in meridional flow speed as a
function of feature lifetimes was attributed to short-lived flows and apparent
motion due to feature-feature interactions. However, the longitudinal
variation of the flow velocities was not analyzed in this study. The fact that
we see unrealistic longitudinal variation in rotation rate along with
meridional flow speed variation rules out short-lived flows as the culprit and
indicates that the source of this discrepancy must be a systematic error
because the Sun cannot know which meridian is centered in HMI magnetograms and
cannot adjust its rotation rate accordingly.

Interestingly, this systematic center-to-limb effect reduces significantly (by
$\sim 50\%$) when unsigned magnetograms are used for tracking as is shown in
the second row of Fig. 5. The third row of Fig. 5 indicates the absence of
these patterns in flow speeds measured from synthetic magnetograms suggesting
that these effects do not originate from our projection or tracking algorithm.
In order to understand why tracking the unsigned magnetic field improves the
measurements, we must take a closer look at the dynamics of the solar
photosphere.

Small scale magnetic features on the photosphere of the Sun tend to accumulate
near converging flows in granules and supergranules and tend to be absent at
cell centers where the plasma flow is diverging (Simon & Leighton, 1964). When
we use the signed value of magnetic field for correlation tracking, we track
the positive and negative magnetic features themselves without giving any
consideration to places where they are not found. When we use the absolute
value of the magnetic field we also track cell centers as negative features
when the average is removed. The reduction in the apparent shift away from
disk center suggests that this shift is directly associated with the magnetic
features themselves.

Even though correlation tracking on unsigned magnetograms shows a reduced
systematic shift, some residual variation of meridional flow persists and the
rotation rate still shows some asymmetry in the east-west direction. So, we
develop a robust methodology for removing systematic errors from flow
measurements which would be generally applicable regardless of whether the
magnetograms are signed or unsigned and even independent of the features
(magnetic, doppler, photometric, etc.) being used for tracking.

We begin by making certain simplifying assumptions about the composition of
measured displacements of blocks. The simplest way would be to assume that
there exists a fixed systematic shift, $\Delta$, away from disk center at all
time-lags. With no intrinsic variation in baseline (true) flow speed ($v$)
with time-lag, this constant shift affects the shorter time-lag measurements
more than the longer time-lag measurements. This, however, proved to be too
simple. When we decomposed the measured displacements for each block at
several time-lags into a displacement due to a baseline velocity ($v\Delta t$)
and a constant shift ($\Delta$) we found different $\Delta$ for different
time-lags which violated our initial assumption. We then added a possible
intrinsic variation in flow speed with time-lag by adding another term,
$\delta$, representing a linear variation in flow speed with time-lag. We thus
included:

  * •

$\Delta$ component: A constant shift independent of time-lag defined as in Eq.
3.

  * •

$\delta$ component: A quantity with the dimensions of acceleration which
accounts for changes in flow speed as a linear function of time-lag.

Using the definitions above, we can write the displacement, $D_{n}$ of a
particular block, detected by correlation tracking at time-lag n$\Delta t$,
where $\Delta t$ is 1-hour for HMI and 96-minutes for MDI, as

$D_{n}=(v+\delta\ n\Delta t)n\Delta t+\Delta$ (3)

so that $\delta$ and $\Delta$ can be determined via the expressions

$\delta_{1-2-4}=(D_{4}-3D_{2}+2D_{1})/(6\Delta t^{2})$ (4)
$\centering\delta_{2-4-8}=(D_{8}-3D_{4}+2D_{2})/(24\Delta
t^{2})\@add@centering$ (5)

$\Delta_{1-2}=2D_{1}-D_{2}+2\delta\ \Delta t^{2}$ (6)
$\Delta_{2-4}=2D_{2}-D_{4}+8\delta\ \Delta t^{2}$ (7)
$\centering\Delta_{4-8}=2D_{4}-D_{8}+32\delta\ \Delta t^{2}.\@add@centering$
(8)

Finally, we can derive “baseline” velocity profiles (a velocity independent of
the systematic shift and time-lag variations)

$v_{1-2}=(D_{2}-D_{1})/(1\Delta t)-3\delta\ \Delta t\\\ $ (9)
$v_{2-4}=(D_{4}-D_{2})/(2\Delta t)-6\delta\ \Delta t\\\ $ (10)
$v_{4-8}=(D_{8}-D_{4})/(4\Delta t)-12\delta\ \Delta t.\\\ $ (11)

Here the subscripts, e.g. 1-2-4, indicate which time-lag displacements were
used to derive the quantity, e.g. 1$\Delta$t, 2$\Delta$t, 4$\Delta$t.

Figure 6:  From left to right: $v$, the baseline meridional flow velocity,
constant shift $\Delta_{MF}$ and the time-lag dependency parameter
$\delta_{MF}$ derived from combinations of measurements at different time-lags
as in Eqs. 4 to 11 and averaged over one year from May 2010 to April 2011 for
HMI. Their agreement supports our assumption that $\Delta$, $\delta$ and the
baseline flow are independent of time-lag. Figure 7: Baseline flow velocity
$v$, $\Delta$, and $\delta$ averaged over one year (from May 2010 to April
2011) and over different combinations of time-lags. The top row shows the
baseline flow velocity (in $m/s$), middle row shows $\Delta$ (in $km$) and
$\delta$ (in $m/s/hr$) is shown in the bottom row. The first and third columns
show the residuals after removing the average quantity at each latitude. The
second column shows the longitudinal dependence of the quantity at the solar
equator and the fourth column shows the latitudinal dependence (averaged over
all longitudes) of the quantity.

The baseline velocity profiles, $v$, time-lag dependency parameter, $\delta$,
and constant shift, $\Delta$, for the meridional flow and rotation rate are
derived from measurements of displacements at several time-lags and averaged
over each Carrington rotation. One year temporal averages of longitudinally
averaged $v,\delta$ & $\Delta$ shown in Fig. 6 do not violate our initial
assumption that these parameters are independent of the combination of time-
lags used to derive them up to time-lags of 8 hr for meridional flow and 4 hr
for rotation rate. The simplified assumption of $\delta$ component
representing a linear variation in baseline velocity with time-lag may not
hold beyond time-lags of 4 hr for rotation rate. This is thought to be due to
the fact that at longer time-lags, our measurements would be biased towards
tracking larger, stronger and longer-lived magnetic features that may be
rooted deeper inside the Sun and hence move with the rotation rate of deeper
layers which does not strictly vary linearly with depth (Antia et al., 1998;
Schou et al., 1998; Howe, 2009).

The spatial distribution of one year averages of the derived baseline flow
profiles ($v$), $\delta$ and $\Delta$ averaged over different combinations of
time-lags is shown in Fig. 7. Amazingly, the systematic center-to-limb error
is automatically captured by $\Delta$ in both meridional flow and rotation
rate while the baseline flow profiles and the time-lag dependency parameter do
not show a significant, unrealistic East-West dependence.

##  4 Results

We now discuss the measurements obtained after making aforementioned
improvements to our LCT algorithm and speculate on the physical meaning of
baseline flow profiles ($v$) and that of $\Delta$ and $\delta$. The constant
shift $\Delta$ shows an unrealistic East-West variation which is a clear
indication that this component is a systematic error that varies as a function
of position on the solar disk just like the center-to-limb effect. The time
histories of the derived components ($v$, $\Delta$ & $\delta$) are shown as
functions of latitude in Fig. 8. Since the latitude at the center of the solar
disk varies annually, we see annual variations in $\Delta$ in Fig. 8. If this
$\Delta$ component were not removed, these annual variations would contaminate
the flow velocity measurements. The time-lag dependency parameter $\delta$
however does not show any unrealistic dependence on longitude in Fig. 7 and
thus is probably not a systematic error but rather an indication of variation
in meridional flow and rotation rate with feature lifetime. The interpretation
and properties of $\Delta$ and $\delta$ are further discussed in section 4.3.

###  4.1 Meridional Flow

The baseline meridional flow velocity averaged during the MDI and HMI overlap
period from May 20, 2010 to Jan 11, 2011 (shown in Fig. 9) was nearly
identical (within error bars) for latitudes $<45^{\circ}$ and no cross-
calibration was required to create the top left panel of Fig. 8. It is
expected that the meridional flow measurements at higher latitudes are more
trustworthy from HMI due to its higher resolution and lower noise level
magnetograms. The footprint of solar active regions is clearly visible in the
meridional flow panel of Fig. 8 with meridional flow being slower (by around 5
$m/s$) at mid-latitudes during the sunspot maximum. This variation is
represented in the fitting coefficient of the associated Legendre polynomial
$P_{2}^{1}(\lambda)=2\sin\lambda\cos\lambda$ to the meridional flow profile
shown in the left panel of the last row of Fig. 8 along with error bars that
span two standard deviations. Note that the temporal trend in this fitting
coefficient for the baseline velocity $v$ is remarkably similar to the fitting
coefficient for 8-hour time-lags reported by Hathaway & Rightmire (2010) but
the amplitude in our measurements is systematically higher. This is due to the
small effect of $\Delta$ and $\delta$ at 8-hour time-lags which reduces the
apparent meridional flow. The fit coefficient shows that the amplitude of the
meridional flow was slowest near the sunspot maximum of solar cycle 23. The
top left panel of Fig. 8 captures the spatiotemporal variations in meridional
flow in greater detail.

###  4.2 Torsional Oscillation

The baseline rotation rate measured in MDI magnetograms during the overlap
period (shown in Fig. 9) also was nearly identical to the one measured from
HMI magnetograms except within $10^{\circ}$ of the equator. Above $50^{\circ}$
latitude, rotation rates obtained from HMI were much less noisy (by a factor
of three) compared to MDI.

We subtracted the average rotation rate during the HMI era from HMI
measurements and the average rotation rate during MDI era from MDI
measurements to derive the torsional oscillation profile (deviations from the
average rotation rate) in Fig. 8. This measured torsional oscillation pattern
is remarkably similar to the one obtained from helioseismic techniques (Howe
et al., 2000; Antia & Basu, 2000; Vorontsov et al., 2002; Howe, 2009) and
shows the equatorward propagating branches of faster and slower than average
rotation rate in great detail while the signal to noise ratio is not high
enough to see the poleward propagation of high-latitude branches clearly. The
equatorward propagating branch with faster than average rotation rate (in red)
which the active regions of solar cycle 25 (expected to start from year 2020)
are expected to follow is very clearly visible before the end of solar cycle
24 centered around $40^{\circ}$ latitude from late 2017 onwards and it appears
disconnected from what looks like a weak patch of faster than average rotation
rate at latitudes above $50^{\circ}$ preceding it.

Figure 8:  Latitude-time plots of corrected flow velocities and systematic
error along with temporal variation of Legendre coefficients fitted to the
flow velocities. The top row shows latitude-time plot of the baseline
meridional flow and torsional oscillation obtained from over 23 years of
measurements from MDI and HMI data in units of $m/s$. The boundaries of active
latitudes are marked in black using the Royal Greenwich Observatory sunspot
records. The corresponding plots for $\Delta$ (in $km$) and $\delta$ (in
$m/s/hr$) are shown in the second and the third row respectively. The bottom
row shows the scaled sunspot number (in black) and our fitting coefficients of
associated Legendre polynomial of order 1, degree 2 to the meridional flow
profile and of order 1, degree 1 to the torsional oscillation (translated by
55 $m/s$) with $2\sigma$ errors in blue whereas the fitting coefficients from
Hathaway & Rightmire (2010) are shown with red dots and standard error.
Figure 9:  Measured quantities temporally averaged during the overlap period
between MDI (blue) and HMI (red). Top row: baseline meridional flow and
rotation rate (in Carrington frame). Middle row: center-to-limb error in MF
and DR. Bottom row: time-lag dependency parameter $\delta$ in MF and DR. The
horizontal axis for $\Delta_{DR}$ is longitude near the equator.

###  4.3 Properties of constant shift error $\Delta$ and time-lag dependency
$\delta$

Interestingly, most of the differences between the displacement measurements
($D_{\Delta t}$) from MDI and HMI (see Fig. 9) were captured by the constant
shift ($\Delta$) i.e. the systematic center-to-limb error which improved the
consistency of measurements across instruments. These differences may be due
to dependence on certain properties of the different spectral lines which are
observed to create the magnetograms e.g.: their height of formation in the
solar atmosphere (Löptien et al., 2016). The fact that we are using different
sets of timelags (multiples of 96 minutes for MDI, multiples of 1 hr for HMI)
may also be responsible for some part of the difference.

Most of the annual variations in uncorrected displacements were also captured
by $\Delta$, leaving behind purified baseline meridional flow and torsional
oscillation profiles shown in Fig. 8. The center-to-limb effect ($\Delta$) is
defined in Eq. 3 and can be re-written as

$lim_{\Delta t\to 0}D_{\Delta t}=\Delta$ (12)

which means it is a constant shift away from disk center that would be there
even if we correlate magnetograms with near zero time-lag. When this constant
shift is divided by time-lag, it dominates the measured velocities at small
time-lags. For HMI at $60^{\circ}$ from disk center, we find a constant shift
of $\approx 100$ $km$ in both the latitudinal and longitudinal displacements
(see Fig. 7). This implies an apparent velocity of $27.8$ $m/s$ at 1 hr and
$3.5$ $m/s$ at 8 hr time-lags. While the contribution from the constant shift
error decreases as a function of time-lag, the contribution from time-lag
dependency parameter $\delta$ increases due to the $\Delta t^{2}$ factor.

The footprint of active regions is clearly visible in the long-term
measurements of $\Delta_{MF}$ and $\Delta_{DR}$ in Fig. 8 suggesting that this
systematic error is higher in the quiet Sun regions than in the vicinity of
active regions. This explains why Lamb (2017) measured a very high meridional
flow for weak features with life times up to two hours which primarily are
found in quiet Sun regions. Löptien et al. (2016) also suggested that the
center-to-limb effect may depend on the magnetic field which agrees with our
observation here. Contamination from this constant shift/center-to-limb error
may also be the cause of deviation between the actual flow measurements and
the fitted functional form of the meridional flow profile shown in Fig. 3 of
Imada & Fujiyama (2018).

The time-lag dependency parameter, on the other hand accounts for variation in
flow velocity with time-lag so that the effective flow velocity at any time-
lag ($v^{\prime}_{\Delta t}$) is slightly different from the baseline flow
velocity as

$v^{\prime}_{\Delta t}=v+\delta\Delta t.$ (13)

The patterns of time-lag dependency parameter in meridional flow $\delta_{MF}$
and rotation rate $\delta_{DR}$ shown in Fig. 8 show some abrupt change
between the two instruments (see Fig. 9). $\delta_{MF}$ hovers around $\sim
0.6$ $m/s/hr$ with a sign opposite to that of the mean meridional flow in each
hemisphere indicating that the meridional flow slows down as a function of
time-lag. Whereas $\delta_{DR}$ is in the range $0.0$ to $2.0$ $m/s/hr$ over
most of the solar disk indicating that the rotation rate increases with time-
lag. The one year average of $\delta_{MF}$ and $\delta_{DR}$ plotted in Fig. 7
does not show any unrealistic longitudinal variation and appears to be
axisymmetric. Therefore, there is no convincing reason to suspect that
$\delta_{MF}$ and $\delta_{DR}$ are some kind of systematic errors.

##  5 Discussion & Conclusion

While it is clear that the constant shift ($\Delta$) comprises of the center-
to-limb error, there isn’t a convincing, straightforward physical
interpretation of time-lag dependency parameter ($\delta$). We speculate that
$\delta_{DR}$ accounts for a physical increase in rotation rate with increase
in time-lag. At longer time-lags the measurements are biased towards tracking
larger, stronger magnetic features with longer lifetimes that are anchored
deeper inside the Sun and hence move with the velocities of the deeper layers
as suggested by Hathaway (2012). The behavior of $\delta_{DR}$ also agrees
with helioseismic inversions of rotation rate which show a near surface shear
layer where the rotation rate increases as we look deeper below the
photosphere (Antia et al., 1998; Schou et al., 1998; Antia et al., 2008; Howe,
2009).

On the other hand, for $\delta_{MF}$, our measurements suggest that meridional
flow decreases with depth below the photosphere. This was first noted by
Hathaway (2012) in the measured meridional motions of supergranules which
indicated that the poleward meridional flow seen at the photosphere decreases
with depth across the surface shear layer and becomes equatorward at the base
of this layer and this was confirmed by helioseismic studies (Zhao et al.,
2013; Jackiewicz et al., 2015; Chen & Zhao, 2017) and Imada & Fujiyama (2018)
as well. However, Basu & Antia (2010) found that the meridional flow speed
actually increased with depth in the outer $2\%$ of the Sun. Even though the
flow velocities derived in the above mentioned studies may be affected by the
center-to-limb effect, their qualitative conclusion (rotation rate increases
with field strength and depth; meridional flow decreases with field strength
and depth) agrees with ours.

In the light of our speculation, the visibility of the footprint of active
regions in the pattern of $\delta_{MF}$ and $\delta_{DR}$ then suggests that
the difference between the velocities of the deeper layers and the photosphere
is less in the presence of high magnetic field strength features. This, again,
is explained by the fact that stronger magnetic features which are rooted
deeper below the photosphere but appear at the photosphere to move with the
velocity of the deeper layers in which they are anchored, reduce the apparent
difference in the velocity between the photosphere and deeper layers in their
neighborhood.

If our speculation on the meaning of $\delta_{DR}$ and $\delta_{MF}$ is
correct then the implied meaning of the baseline flow velocities is that: this
is the velocity with which the weakest magnetic features move and the stronger
magnetic features (predominantly tracked at longer time-lags) move with
$v^{\prime}_{\Delta t}$ (Eq. 13).

Regardless of whether our speculations are correct, the improvements described
in sec 3 are generally applicable to many other tracking algorithms and they
have increased the accuracy of our LCT algorithm five fold so that it can
detect a $0.02$ pixel shift reliably. Our proposed methodology for the removal
of the center-to-limb effect is computationally more expensive than the
traditional approach (removing some constant profile) as it utilizes
measurements of displacements at several time-lags but it works with very few
assumptions. To use it effectively, one does not have to assume that the
center-to-limb effect is symmetric across the entire solar disk, that it is
independent of the distribution of magnetic flux on solar disk, or that it is
constant in time (Zhao et al., 2012; Komm et al., 2015). The variation of
center-to-limb error with the phase of the solar cycle seen in Fig. 8 clearly
indicates that the center-to-limb effect cannot be assumed to be constant in
time and it has to be calculated synchronously with the flow profiles. If it
is assumed to be a constant, its temporal variation may contaminate measured
flow profiles.

All our final baseline velocity, $\Delta$ and $\delta$ measurements as well as
the Legendre fitting coefficients to baseline flow profiles shown in Fig. 8
are available in an easy to use form at
https://dataverse.harvard.edu/dataverse/lct-on-solar-magnetograms (Mahajan,
2020a, b).

The authors were supported by grant NNX14AO83G from the NASA Heliophysics
Grand Challenges Program. The MDI data used are courtesy of the ESA-NASA/SOHO
and the MDI science teams. The HMI data used are courtesy of the NASA/SDO and
the HMI science teams. S.S.M. is grateful to NASA Ames Research Center for
hosting him during the summer of 2016. The authors are also grateful to Lisa
Upton for her insights and discussions that helped this paper.

## References

  * Antia & Basu (2000) Antia, H. M., & Basu, S. 2000, ApJ, 541, 442
  * Antia et al. (1998) Antia, H. M., Basu, S., & Chitre, S. M. 1998, MNRAS, 298, 543
  * Antia et al. (2008) —. 2008, ApJ, 681, 680
  * Basu & Antia (2010) Basu, S., & Antia, H. M. 2010, ApJ, 717, 488
  * Beck (2000) Beck, J. G. 2000, Sol. Phys., 191, 47
  * Carrington (1859) Carrington, R. C. 1859, Monthly Notices of the Royal Astronomical Society, 19, 81
  * Charbonneau (2010) Charbonneau, P. 2010, Living Reviews in Solar Physics, 7, 3
  * Chen & Katz (2005) Chen, J., & Katz, J. 2005, Measurement Science and Technology, 16, 1605
  * Chen & Zhao (2017) Chen, R., & Zhao, J. 2017, ApJ, 849, 144
  * Dikpati et al. (2010) Dikpati, M., Gilman, P. A., & Ulrich, R. K. 2010, ApJ, 722, 774
  * Hathaway (2012) Hathaway, D. H. 2012, ApJ, 760, 84
  * Hathaway & Rightmire (2010) Hathaway, D. H., & Rightmire, L. 2010, Science, 327, 1350
  * Hathaway & Rightmire (2011) —. 2011, ApJ, 729, 80
  * Hathaway et al. (2006) Hathaway, D. H., Williams, P. E., & Cuntz, M. 2006, ApJ, 644, 598
  * Hathaway et al. (1996) Hathaway, D. H., Gilman, P. A., Harvey, J. W., et al. 1996, Science, 272, 1306
  * Howe (2009) Howe, R. 2009, Living Reviews in Solar Physics, 6, 1
  * Howe et al. (2000) Howe, R., Christensen-Dalsgaard, J., Hill, F., et al. 2000, ApJ, 533, L163
  * Huang et al. (1997) Huang, H., Dabiri, D., & Gharib, M. 1997, Measurement Science and Technology, 8, 1427
  * Imada & Fujiyama (2018) Imada, S., & Fujiyama, M. 2018, ApJ, 864, L5
  * Jackiewicz et al. (2015) Jackiewicz, J., Serebryanskiy, A., & Kholikov, S. 2015, ApJ, 805, 133
  * Jiang et al. (2010) Jiang, J., Isik, E., Cameron, R. H., Schmitt, D., & Schüssler, M. 2010, ApJ, 717, 597
  * Komm et al. (2015) Komm, R., González Hernández, I., Howe, R., & Hill, F. 2015, Sol. Phys., 290, 1081
  * Komm et al. (1993) Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys., 147, 207
  * Lamb (2017) Lamb, D. A. 2017, ApJ, 836, 10
  * Langfellner et al. (2015) Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67
  * Liang et al. (2018) Liang, Z.-C., Gizon, L., Birch, A. C., Duvall, T. L., & Rajaguru, S. P. 2018, A&A, 619, A99
  * Löptien et al. (2016) Löptien, B., Birch, A. C., Duvall, T. L., Gizon, L., & Schou, J. 2016, A&A, 590, A130
  * Mahajan (2020a) Mahajan, S. 2020a, Level 3 Flow Measurements: Longitudinally averaged f(latitude, time), vDRAFT VERSION, Harvard Dataverse, doi:10.7910/DVN/I0OWHG.  https://doi.org/10.7910/DVN/I0OWHG
  * Mahajan (2020b) —. 2020b, Level 4: Fitted Photospheric Flow Profiles for modelers, vDRAFT VERSION, Harvard Dataverse, doi:10.7910/DVN/0IZBHY.  https://doi.org/10.7910/DVN/0IZBHY
  * Muñoz-Jaramillo et al. (2013) Muñoz-Jaramillo, A., Dasi-Espuig, M., Balmaceda, L. A., & DeLuca, E. E. 2013, ApJ, 767, L25
  * Sandwell (1987) Sandwell, D. T. 1987, Geophys. Res. Lett., 14, 139
  * Scherrer et al. (1995) Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * Schou et al. (1998) Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390
  * Sheeley et al. (1989) Sheeley, N. R., J., Wang, Y. M., & DeVore, C. R. 1989, Sol. Phys., 124, 1
  * Simon & Leighton (1964) Simon, G. W., & Leighton, R. B. 1964, ApJ, 140, 1120
  * Ulrich (2010) Ulrich, R. K. 2010, ApJ, 725, 658
  * Vorontsov et al. (2002) Vorontsov, S. V., Christensen-Dalsgaard, J., Schou, J., Strakhov, V. N., & Thompson, M. J. 2002, Science, 296, 101
  * Westerweel (1997) Westerweel, J. 1997, Measurement Science and Technology, 8, 1379
  * Zhao et al. (2013) Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, T. L., J., & Hartlep, T. 2013, ApJ, 774, L29
  * Zhao et al. (2012) Zhao, J., Nagashima, K., Bogart, R. S., Kosovichev, A. G., & Duvall, Jr., T. L. 2012, ApJ, 749, L5