# Minuscule corrections to near-surface solar internal rotation using mode-
coupling

∗Srijan Bharati Das Department of Geosciences
Princeton University
Princeton, New Jersey, USA  ∗Samarth G. Kashyap Department of Astronomy and
Astrophysics
Tata Institute of Fundamental Research
Mumbai, India  Deniz Oktay Department of Computer Science
Princeton University
Princeton, New Jersey, USA  Shravan M. Hanasoge Department of Astronomy and
Astrophysics
Tata Institute of Fundamental Research
Mumbai, India  Jeroen Tromp Department of Geosciences
and Program in Applied & Computational Mathematics
Princeton University
Princeton, New Jersey, USA

(Received 2022 July 18; Revised 2023 March 18; Accepted 2023 March 23)

###### Abstract

The observed solar oscillation spectrum is influenced by internal
perturbations such as flows and structural asphericities. These features
induce splitting of characteristic frequencies and distort the resonant-mode
eigenfunctions. Global axisymmertric flow — differential rotation — is a very
prominent perturbation. Tightly constrained rotation profiles as a function of
latitude and radius are products of established helioseismic pipelines that
use observed Dopplergrams to generate frequency-splitting measurements at high
precision. However, the inference of rotation using frequency-splittings do
not consider the effect of mode-coupling. This approximation worsens for high-
angular-degree modes, as they become increasingly proximal in frequency. Since
modes with high angular degrees probe the near-surface layers of the Sun,
inversions considering coupled modes could potentially lead to more accurate
estimates of rotation very close to the surface. In order to investigate if
this is indeed the case, we perform inversions for solar differential
rotation, considering coupling of modes for angular degrees $160\leq\ell\leq
300$ in the surface gravity $f$-branch and first-overtone $p$ modes. In
keeping with the character of mode coupling, we carry out a non-linear
inversion using an eigenvalue solver. Differences in inverted profiles for
frequency splitting measurements from MDI and HMI are compared and discussed.
We find that corrections to the near-surface differential rotation profile,
when accounting for mode-coupling effects, are smaller than 0.003 nHz and
hence are insignificant. These minuscule corrections are found to be
correlated with the solar cycle. We also present corrections to even-order
splitting coefficients, which could consequently impact inversions for
structure and magnetic fields.

Sun: helioseismology — Sun: oscillations — Sun: interior — differential
rotation — QDPT

††journal: ApJS††thanks: Both authors have contributed equally to this study.

##  1 Introduction

Helioseismology has enabled high-precision measurements of solar internal
rotation. Turbulence in the convection zone excites modes of oscillation which
propagate through the solar interior and are sensitive to the prevalent
structure and flows. The extent to which these modes of oscillation, called
the solar normal modes, are coupled depends on their proximity in frequency.
Rotation breaks spherical symmetry, resulting in prograde modes with higher
frequencies and retrograde modes with lower frequencies than the corresponding
eigenfrequencies predicted in static, non-rotating standard models (such as
Model-S as defined by Christensen-Dalsgaard et al., 1996). This is called
rotational frequency splitting.

Degenerate perturbation theory (DPT; Lavely & Ritzwoller, 1992) has
traditionally been employed for estimating rotation, where coupling between
distinctly different modes is ignored. DPT has been employed for inversions
using helioseismic data from BBSO (Libbrecht, 1989; Brown et al., 1989), GONG
(e.g., Thompson et al., 1996), MDI (e.g., Kosovichev et al., 1997; Schou et
al., 1998) and HMI (e.g., Larson & Schou, 2018). Differential rotation (DR) is
an extensively studied feature in the Sun. It is beyond the scope of this
paper to discuss the many important studies on DR, and the reader is referred
to Howe (2008) for a comprehensive review of all notable contributions. Of
particular importance, from a methodological perspective, is Schou et al.
(1998), which compared inversions via seven different methods and identified
various robust properties of differential rotation. These studies culminated
in the consensus that solar internal rotation is zonal with a solidly rotating
core, a differentially rotating convection zone with two radial shear layers,
one at the bottom of the convection zone, termed the “tachocline” (Spiegel &
Zahn, 1992) and one near the surface, termed the near-surface shear layer
(NSSL; Thompson et al., 1996). These radial shear layers have drawn attention
given their potential importance for driving the solar dynamo and their
prominent role in maintaining the global angular momentum budget.

Quasi-degenerate perturbation theory (QDPT; Lavely & Ritzwoller, 1992, ,
henceforth LR92), on the other hand, accounts for cross-coupling between modes
when computing frequency splittings. Although it represents a more accurate
model, very few studies have employed QDPT in inferring DR because of its
computational complexity as compared to DPT. Schad & Roth (2020) formulated
measurements in terms of mode-amplitude ratios to infer DR. Woodard et al.
(2013) and Kashyap et al. (2021) fit mode-amplitude spectra using the
prescription of Vorontsov (2011) (hereafter V11) to infer $a$-coefficients,
which are polynomial expansion coefficients of frequency splittings. Higher-
angular-degree modes are closely spaced in frequency, implying the worsening
of the underlying assumption of DPT. Kashyap et al. (2021) showed that the
difference in splittings estimated by DPT and QDPT are statistically
significant for larger angular degrees in the $f$ and $p_{1}$ branches. Here,
we attempt to infer the rotation profile through the application of the more
general QDPT formalism, while still using the $a$-coefficients as our primary
measurement. Given that only high-angular degrees are coupled, we propose a
“hybrid” inversion method, where DPT is used for low-angular-degree modes and
QDPT is used for high-angular-degree modes. Low-angular-degree modes are
sensitive to greater depths, while the high angular degrees are trapped very
close to the surface. Consequently, we expect only near-surface corrections
from our study.

In this study, we use almost 22-years of $a$-coefficient measurements to infer
DR. The period between May 1996 and April 2010 was constrained by measurements
from MDI (Scherrer et al., 1995), onboard SOHO, and the period between April
2010 and January 2018 was constrained by measurements from HMI (Schou et al.,
2012), onboard the SDO. While the odd $a$-coefficients are almost completely
governed by differential rotation, the even $a$-coefficients contain
contributions from structure perturbations such as solar oblateness (Woodard,
2016), magnetic fields (Antia et al., 2000; Baldner et al., 2009), sound-speed
asphericity (Antia et al., 2001; Baldner & Basu, 2008) and second-order
effects of rotation (Gough & Thompson, 1990). We carry out forward
calculations via an exact eigenvalue solver to estimate first-order
contributions to even $a$-coefficients due to solar rotation and compare our
estimates with V11, which adopted a semi-analytic approach.

The outline of this paper is as follows. In Section 2, we introduce basic
notations for mode coupling in the context of DR as well as the
$a$-coefficient formalism. Section 3 furnishes details pertaining to (a)
categorization of QDPT and DPT modes in Section 3.1, (b) setting up the cost-
function to be minimized during inversion in Section 3.2, and, (c) determining
optimal truncation in the angular degree of perturbation and the spectral
window of coupled modes in Section 3.3. We present the results of corrections
for time variation in DR in Section 4.1 and compare even-splitting
coefficients with V11 results in Section 4.2. We summarize the key points in
Section 5.

##  2 Theoretical formulation: Isolated and coupled multiplets

Linear perturbation analysis of the hydrodynamic equations of mass continuity,
conservation of momentum, and energy results in an eigenvalue problem which
enables calculation of eigenfrequencies and eigenfunctions of standard solar
models, such as Model-S (Christensen-Dalsgaard et al., 1996). Such standard
models do not account for effects of asphericity, flows, anisotropy, non-
adiabaticity, and magnetic fields. In the absence of such perturbations, these
standard models have theoretically predicted “degenerate” eigenfrequencies
${}_{n}\omega_{\ell}$ and eigenfunctions ${}_{n}\xi_{\ell m}$. Here $n$ is the
radial order and $\ell$ is the spherical harmonic degree. For a given
multiplet ${}_{n}\mathrm{S}_{\ell}\,$, the constituent $2\ell+1$ modes are
labelled by a third quantum number $m\in[-\ell,\ell]$. The eigenfrequencies
are degenerate in $m$ because of spherical symmetry.

The influence of the above-mentioned effects needs to be accounted for as
additional perturbations to the standard models. This results in splitting of
eigenfrequencies and distortion of the eigenfunctions as follows

${}_{n}\omega_{\ell}\to\omega_{\mathrm{ref}}+\delta{}_{n}\omega_{\ell
m}\,,\qquad{}_{n}\mbox{\boldmath$\bf\xi$}_{\ell
m}\to{}_{n}\mbox{\boldmath$\bf\xi$}_{\ell
m}+\delta{}_{n}\mbox{\boldmath$\bf\xi$}_{\ell m}\,,$ (1)

where $\delta{}_{n}\omega_{\ell m}$ and
$\delta{}_{n}\mbox{\boldmath$\bf\xi$}_{\ell m}$ are the corrections to
eigenfrequencies and eigenfunctions and $\omega_{\mathrm{ref}}$ is a reference
frequency which maybe chosen to be close to the unperturbed frequency
${}_{n}\omega_{\ell}$ of the multiplet ${}_{n}\mathrm{S}_{\ell}\,$ whose
perturbed eigenstate we are interested in. This results a new eigenvalue
problem for the “supermatrix” ${\sf Z}$ (see Appendix B in Das et al. (2020))

$\sum_{k\in\mathcal{K}}Z_{k^{\prime}k}\,c_{k}=2\omega_{\mathrm{ref}}\,\delta\omega_{k^{\prime}}\,c_{k^{\prime}}\,.$
(2)

Elements of ${\sf Z}$ encode the coupling of modes in the presence of the
aforementioned perturbations. Here, $k$ is conveniently used as a combined
index to denote the mode $(n,\ell,m)$. For a temporal bandwidth $\Delta\omega$
and spectral bandwidth $\Delta\ell$, the supermatrix ${\sf Z}$ is built around
a “central multiplet” ${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$ by considering a set
of modes $\mathcal{K}$ that obey $k\in\mathcal{K}$,
$|{}_{n}\omega_{\ell}-{}_{n_{0}}\omega_{\ell_{0}}|<\Delta\omega$ and
$|\ell-\ell_{0}|<\Delta\ell$. A sufficiently large $\Delta\omega$ and
$\Delta\ell$ ensures that all significant couplings with
${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$ are accounted for (see Section 3.3). The
perturbed eigenfrequencies of the $2\ell_{0}+1$ modes estimated from ${\sf Z}$
thereby accurately account for cross-coupling of multiplet $k_{0}$ with its
proximal neighbours. The unperturbed mass density is denoted by $\rho_{0}$
(Model S). Without loss of generality, we have chosen
$\omega_{\mathrm{ref}}=\omega_{k_{0}}$, for our inversions.

Following the convention in LR92 and V11, we represent the 3D rotational
velocity as

$\mathbf{v}_{\mathrm{rot}}(r,\theta,\phi)=-\sum_{s}w_{s}(r)\,\hat{\mathbf{r}}\times\mbox{\boldmath$\bf\nabla$}_{1}Y_{s,0}(\theta,\phi)\,,$
(3)

where $w_{s}(r)$ are the respective odd-degree toroidal coefficients and
$Y_{\ell,m}(\theta,\phi)$ are spherical harmonics labelled by angular degree
$\ell$ and azimuthal order $m$. The elements of the supermatrix due to the
rotation field ${\mathbf{v}}_{\mathrm{rot}}$ may be expressed as (see Eqns.
[135-136] in LR92 and Eqn. [A1] in V11)

$\begin{split}Z^{(n_{0},\ell_{0})}_{k^{\prime}k,m}=2\omega_{\mathrm{ref}}\,\kappa_{\ell^{\prime}}\kappa_{\ell}\sum_{s}\gamma_{s}\,\bigg{(}\begin{smallmatrix}\ell^{\prime}&s&\ell\\\
1&0&-1\end{smallmatrix}\bigg{)}\,\bigg{(}\begin{smallmatrix}\ell^{\prime}&s&\ell\\\
m&0&-m\end{smallmatrix}\bigg{)}\,\int_{r=0}^{R_{\odot}}\rho_{0}\,w_{s}(r)\,T^{k^{\prime}k}_{s}(r)\,r\,{\mathrm{d}}r+\left(\omega_{k}^{2}-\omega_{\rm{ref}}^{2}\right)\,\delta_{k^{\prime}k}.\end{split}$
(4)

Note that the superscript $(n_{0},\ell_{0})$ denotes the particular central
multiplet for which the supermatrix is constructed and the subscript $m$ is
used to explictly imply that only modes with the same azimuthal order are
coupled in the presence of an axisymmetric flow field. We also have
$\kappa_{\ell}=\left[\ell(\ell+1)(2\ell+1)\right]^{1/2}$,
$\gamma_{s}=\sqrt{(2s+1)/4\pi}$ and $T^{k^{\prime}k}_{s}(r)$, the sensitivity
kernel for $w_{s}(r)$, given by

$T^{k^{\prime}k}_{s}(r)=\left[U_{k^{\prime}}U_{k}-U_{k^{\prime}}V_{k}-V_{k^{\prime}}U_{k}+\frac{\ell^{\prime}(\ell^{\prime}+1)+\ell(\ell+1)-s(s+1)}{2}V_{k^{\prime}}V_{k}\right].$
(5)

$U_{k}(r)$ and $V_{k}(r)$ are the radial and horizontal eigenfunctions
corresponding to mode $k$, respectively. $T_{s}^{k^{\prime}k}(r_{0})$
therefore encodes the degree to which the multiplets $k$ and $k^{\prime}$
couple in the presence of an axisymmetric, degree-$s$ rotation field of unit
strength at $r=r_{0}$.

The eigenvalues of acoustic modes in Model-S are degenerate in $m$, i.e.,
${}_{n_{0}}\omega_{\ell_{0}m}={}_{n_{0}}\omega_{\ell_{0}}.$ Differential
rotation is an axisymmetric perturbation that breaks the spherical symmetry of
the system, thus lifting the degeneracy in $m$ and “splitting” the
frequencies. As shown in Eqn. (2), the frequency splitting associated with
mode $k_{0}$ may be estimated from eigenvalues of the supermatrix
$Z^{(n_{0},\ell_{0})}_{kk^{\prime},m}$. The supermatrix for
${}_{0}\mathrm{S}_{200}\,$ due to differential rotation, considering couplings
for $196\leq\ell\leq 204$, is shown in Fig. 1.

It is standard practice in helioseismology to project frequency splittings
$\delta{}_{n}\omega{}_{\ell m}={}_{n}\omega{}_{\ell m}-{}_{n}\omega{}_{\ell}$
on to a basis of orthogonal polynomials $\mathcal{P}^{(\ell)}_{j}(m)$. The
resultant fitting coefficients $a_{j}^{n\ell}$ are the so-called
“$a$-coefficients”:

${}_{n}\omega{}_{\ell
m}={}_{n}\omega{}_{\ell}+\sum_{j=0}^{j_{\text{max}}}a^{n\ell}_{j}\,\mathcal{P}^{(\ell)}_{j}(m).$
(6)

In practice, $a$-coefficients are recorded for $j_{\text{max}}=36$ (e.g.,
Schou, 1999). A recipe for obtaining these may be found in Appendix A of Schou
et al. (1994). Harnessing the orthonormality of the basis polynomials,
$\sum_{m=-\ell}^{\ell}\mathcal{P}^{(\ell)}_{j}(m)\mathcal{P}^{(\ell)}_{k}(m)=\delta_{jk}$,
we write the $a$-coefficients as

$a^{n\ell}_{j}=\sum_{m=-\ell}^{\ell}\delta{}_{n}\omega{}_{\ell
m}\,\mathcal{P}^{(\ell)}_{j}(m).$ (7)

The isolated multiplet approximation for a zonal perturbation implies that the
frequency splittings $\delta_{n}\omega_{\ell m}$ in Eqn. (7) are equal to the
diagonal elements of $Z^{(n_{0},\ell_{0})}_{k_{0}k_{0},m}$ scaled by
$1/(2\,{}_{n_{0}}\omega_{\ell_{0}})$. While this is a good approximation for
low angular-degree multiplets, Kashyap et al. (2021) showed that there is an
error (in an L2-norm sense) that is significant enough to be detectable (about
2$\sigma$), when using the isolated-multiplet approximation for high-$\ell$
multiplets. This error arises when the off-diagonal components of the
supermatrix (as shown in Fig. 1 for the central multiplet
${}_{0}\mathrm{S}_{200}\,$) are ignored. The rest of the paper will explore
the effects of considering these off-diagonal elements on the $a$-coefficients
predicted from forward calculations and the systematic changes in inverted
profiles from 22 years of MDI and HMI measurements.

Figure 1: Supermatrix corresponding to the central multiplet
${}_{0}\mathrm{S}_{200}\,$. The colors are spaced out in logarithmic scale for
better visibility of the weak off-diagonal coupling components. Grids in the
matrix denote the submatrices capturing coupling of multiplets, which can be
read off from the top of the column and left of the row.
${}_{0}\mathrm{S}_{200}\,-{}_{0}\mathrm{S}_{|200\pm 2|}\,$ couplings are 100
times weaker, whereas ${}_{0}\mathrm{S}_{200}\,-{}_{0}\mathrm{S}_{|200\pm
4|}\,$ couplings are 1000 times weaker than self-coupling. We have used the
symlog functionality of the matplotlib Python package to represent negative
values in logarithmic scale.

##  3 Inversion

###  3.1 Selection of coupled multiplets

We carry out inversions for differential rotation using $a$-coefficients
obtained from 72-day data-sets over the period between 1996-05-01 and
2018-01-06 (all dates in this study use the yyyy-mm-dd format). We do away
with the isolated multiplet approximation for modes with $\ell\geq 160$ in the
$n=0,1$ radial branches. The L2-norm error decreases rapidly with decreasing
$\ell$ (Kashyap et al., 2021). High-precision observations from HMI possess
uncertainties as small as 0.03% of the observed $a$-coefficients. Hence, a
safe cutoff for the L2-norm error is taken to be $10^{-4}$ from Fig. 9 of
Kashyap et al. (2021), and all modes with larger errors are considered for the
full-coupling problem. Henceforth, we refer to these as “coupled” multiplets,
while the others are “isolated” multiplets. The inversions are “hybrid” in
nature since frequency splittings are modeled in two ways:

  1. (A)

For isolated multiplets, the main diagonal of the submatrix corresponding to
the central multiplet $k_{0}$ is used, i.e., $\delta_{n}\omega_{\ell m}=$
diag$\left(Z_{k_{0},k_{0}}^{(n_{0},\ell_{0})}\right)/(2\,{}_{n_{0}}\omega_{\ell_{0}})$.
We drop the superscript label $(n_{0},\ell_{0})$ on the supermatrix from now
on. It is implied that any supermatrix corresponds to a central multiplet
${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$.

  2. (B)

For coupled multiplets, the eigenvalue problem is solved for the entire
supermatrix ${\sf Z}$ and frequency splittings corresponding to the central
multiplet ${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$ are read off. This one-to-one
mapping between eigenfrequencies and the corresponding modes is possible since
the perturbed eigenfunctions are very close to the unperturbed analogs —
arising from the diagonally dominant nature of the supermatrix. We use float32
for all computation, since any difference with higher precision is
significantly lower than observational noise.

In hybrid inversions, based on our defined cutoff, we see that the typical
number of observed multiplets is $\sim 1800$, and $\sim 200$ of those are
found to be coupled, i.e., roughly 10% of all the observed modes are treated
as coupled, and these are responsible for corrections to the DPT-rotation
profiles. To understand the depth up to which coupled multiplets can result in
such corrections, we first estimate the lower turning points of coupled modes.
Invoking the Cowling approximation, it may be shown that the lower turning
point $r_{\mathrm{LTP}}$ of modes ($n>0$) may be expressed as

$\frac{c^{2}(r_{\mathrm{LTP}})}{r_{\mathrm{LTP}}^{2}}=\frac{\omega^{2}}{\ell\,(\ell+1)},$
(8)

where $c(r)$ is the sound speed at radial distance $r$ from the center of the
Sun and $\omega$ the mode frequency. The lowest angular degree for the two
radial branches $n=0,1$ where mode-coupling is considered is $\ell=160$. The
above equation is valid for $n>0$ radial orders. So, using Eqn. (8) to
estimate the $r_{\mathrm{LTP}}$ for ${}_{1}\mathrm{S}_{160}\,$, we get
$r^{n=1}_{\mathrm{LTP}}=0.971\,R_{\odot}$. Since eigenfunctions corresponding
to the $n=0$ branch, has no nodes in radius, it peaks close to the surface and
dies off in an evanescent manner. At $r=0.971\,R_{\odot}$, the
${}_{0}\mathrm{S}_{160}\,$ mode eigenfunctions are approximately two orders of
magnitude weaker than the maxima near the surface. Going deeper to
$r=0.9\,R_{\odot}$, we see that these eigenfunctions are six orders of
magnitude smaller than the near-surface maximum. It is therefore expected that
differential rotation below $0.9\,R_{\odot}$ is weakly sensitive to the effect
of coupled modes. In our inversions, we fit for differential rotation in the
spherical shell bounded by $r\in[0.9\,R_{\odot},1.0\,R_{\odot}]$.

###  3.2 Inversion methodology

We parameterize $w_{s}(r)$ on a radial grid on a basis comprising cubic
B-splines $\beta_{p}(r)$, where index $p$ refers to the knot location
corresponding to a specific B-spline polynomial with coefficient $c_{s}^{p}$,
i.e.,

$w_{s}(r)=\begin{cases}\sum_{p}c_{s}^{p,\,\mathrm{fixed}}\,\beta_{p}(r)\qquad\mathrm{for}\qquad
r<0.9\\\ \sum_{p}c_{s}^{p}\,\beta_{p}(r)\qquad\qquad\mathrm{for}\qquad r\geq
0.9\end{cases},$ (9)

where $c_{s}^{p,\,\mathrm{fixed}}$ cooresponds to the spline coefficients of
the profiles of 2D RLS inversions. We use the same radial grid as the 2D RLS
inversions hosted on Stanford University’s Joint Science Operations Center
(JSOC) database 111Radial grid on JSOC.. Therefore, the inverse problem of
finding $\Omega(r,\theta,\phi)$ reduces to inferring the coefficients
$c_{s}^{p}$. Using Eqns. (4) and (9), we express the supermatrix more
explicitly in terms of the spline coefficients as

$Z_{k^{\prime}k}=\sum_{p}\sum_{s=1,3,5}c_{s}^{p}\,\Lambda_{p,k^{\prime}k}^{s}+\left(\omega_{k}^{2}-\omega_{\rm{ref}}^{2}\right)\,\delta_{k^{\prime}k},$
(10)

where,

$\Lambda_{p,k^{\prime}k}^{s}=2\,\omega_{\mathrm{ref}}\,\kappa_{\ell^{\prime}}\,\kappa_{\ell}\,\gamma_{s}\,\bigg{(}\begin{smallmatrix}\ell^{\prime}&s&\ell\\\
1&0&-1\end{smallmatrix}\bigg{)}\,\bigg{(}\begin{smallmatrix}\ell^{\prime}&s&\ell\\\
m&0&-m\end{smallmatrix}\bigg{)}\,\int_{r=0}^{R_{\odot}}\rho_{0}(r)\,\beta_{p}(r)\,T^{k^{\prime}k}_{s}(r)\,r\,{\mathrm{d}}r.$
(11)

The separation of terms that do not depend on $c_{s}^{p}$ into
$\Lambda_{p,k^{\prime}k}^{s}$ is essential. This allows for a one-time
precomputation of $\Lambda_{p,k^{\prime}k}^{s}$, saving both time and
computational expense during the non-linear inversion when considering coupled
modes. Since we only fit for differential rotation for $r\geq 0.9$, we
classify the spline coefficients $c_{s}^{p}$ into a fixed component, denoted
by $c_{s}^{p,\,\mathrm{fixed}}$, corresponding to the basis functions with
local support in $r<0.9$ and the inverted component, denoted by
$c_{s}^{p,\,\mathrm{fit}}$, corresponding to the basis functions with local
support in $r\geq 0.9$. The $c_{s}^{p,\,\mathrm{fixed}}$ is computed using
profiles available on JSOC. Accordingly, $\Lambda_{p,k^{\prime}k}^{s}$ would
also split up into $\Lambda_{p,k^{\prime}k}^{s,\,\mathrm{fixed}}$ and
$\Lambda_{p,k^{\prime}k}^{s,\,\mathrm{fit}}$. The supermatrix in Eqn. (10) may
then be expressed in terms of a fixed part $Z_{k^{\prime}k}^{\mathrm{fixed}}$
and the $c_{s}^{p,\,\mathrm{fit}}$ dependent part,

$Z_{k^{\prime}k}=Z_{k^{\prime}k}^{\mathrm{fixed}}+\sum_{p}\sum_{s=1,3,5}c_{s}^{p,\,\mathrm{fit}}\,\Lambda_{p,k^{\prime}k}^{s,\,\mathrm{fit}},$
(12)

where
$Z_{k^{\prime}k}^{\mathrm{fixed}}=\sum_{p}\sum_{s=1,3,5}c_{s}^{p,\,\mathrm{fixed}}\,\Lambda_{p,k^{\prime}k}^{s,\,\mathrm{fixed}}+\left(\omega_{k}^{2}-\omega_{\rm{ref}}^{2}\right)\,\delta_{k^{\prime}k}$
and is precomputed along with $\Lambda_{p,k^{\prime}k}^{s,\,\mathrm{fit}}$.
Finally, the modeled data $d_{\mathrm{mod}}$ is expressed as

$d_{\mathrm{mod}}=\begin{cases}\mathcal{A}\left[\texttt{diag}(Z)\right]/(2\,\omega_{k_{0}})&\text{for
isolated modes,}\\\
\mathcal{A}\left[\texttt{eig}(Z)\right]/(2\,\omega_{k_{0}})&\text{for coupled
modes},\end{cases}$ (13)

where $\texttt{diag}(Z)$ returns elements on the main diagonal of the
supermatrix $Z$ corresponding to the self-coupling of the central multiplet
${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,-\,{}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$, and
$\texttt{eig}(Z)$ returns eigenvalues corresponding to the central multiplet
${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$ after solving the eigenvalue problem for
the supermatrix $Z$. The operator $\mathcal{A}$ indicates the selection of
modes belonging to the central multiplet ${}_{n_{0}}\mathrm{S}_{\ell_{0}}\,$.

For the inversion, we define a regularized misfit function $\chi^{2}$ in an
L2-norm sense

$\chi^{2}=\overbrace{\sum_{i\,\in\,\mathrm{all\,modes}}\left(\frac{d^{i}-d^{i}_{\mathrm{mod}}}{\sigma_{d}^{i}}\right)^{2}}^{\chi^{2}_{d}=\text{data
misfit}}+\sum_{s=1,3,5,7,9}\mu_{s}\overbrace{\left[\frac{{\mathrm{d}}^{2}w_{s}(r)}{{\mathrm{d}}r^{2}}\right]^{2}}^{\chi^{2}_{m}=\text{model
smoothness}},$ (14)

where $d_{i}$ and $\sigma_{d}^{i}$ are the $a$-coefficients and their
corresponding uncertainties, measured from HMI or MDI, and $\mu_{s}$ are
angular-degree-dependent regularization parameters. The second term dictates
the smoothness of the inverted profile. When using splines, it may be shown
that this term can be recast in terms of the spline coefficients $c_{s}^{p}$
and an operator $D$ which captures the second derivative of the basis
functions $\beta_{p}(r)$ with respect to $r$:

$\left[\frac{{\mathrm{d}}^{2}w_{s}(r)}{{\mathrm{d}}r^{2}}\right]^{2}=\sum_{i,j\in\mathrm{all\,knots}}c_{s}^{i}\,D^{ij}\,c_{s}^{j}.$
(15)

We further use a non-dimensional regularization parameter
$\overline{\mu}_{s}$, which is related to $\mu_{s}$ as below

$\mu_{s}=\overline{\mu}_{s}\,N\,\frac{\mathrm{tr}(H_{d})}{\mathrm{tr}(H_{m})},$
(16)

where $N$ is the total number of data $a$-coefficients $d^{i}$, $H_{d}$ the
Hessian associated with the data misfit $\chi^{2}_{d}$, and $H_{m}$ the
Hessian associated with the model misfit $\chi^{2}_{m}$. These Hessians are
computed as second derivatives of the respective $\chi^{2}$ with respect to
the vector of the fitted spline coefficients $c_{s}^{p,\,\mathrm{fit}}$. The
inverse problem is non-linear due to the eigenvalue operation as shown in Eqn.
(13) and the model parameters $c_{s}^{p,\,\mathrm{fit}}$ undergo an iterative
march towards the final solution via the standard Newton’s method to minimize
the misfit, $\chi^{2}$ (Tarantola, 1987). The marching in parameter-space
requires the computation of gradient and hessian. These quantities are
computed numerically using the autograd functionality of jax (Bradbury et al.,
2018). Further discussions on the accuracy and robustness on the eigenvalue
operator and inversion tools can be found in Appendix C.

The non-dimensional regularization parameter $\overline{\mu}_{s}$ used for the
inversions presented in this paper is chosen by the traditional L-curve
method. For this, we chose 50 logarithmically-spaced values of
$\overline{\mu}_{s}$ in a sufficiently large window of $[10^{-6},10^{6}]$ for
each $s$. For each of these values of $\overline{\mu}_{s}$, we perform the
inversion and calculate the data misfit $\chi_{d}^{2}$ and model misfit
$\chi^{2}_{m}$ from the inverted model parameters. Subsequently, we made the
L-curve by plotting $\chi^{2}_{d}$ vs. $\chi^{2}_{m}$ and interpolated the
curve with smooth cubic splines using interp1d function of Python’s scipy
module. The $\overline{\mu}_{s}$ for the final inversion results presented in
this paper was chosen from the knee of this L-curve. A representative plot for
the L-curve and a demonstration of the convergence for our hybrid inversions
is provided in Appendix D.

###  3.3 Determining $s_{\mathrm{max}}$ of $w_{s}(r)$ and supermatrix
dimensionality

Figure 2: Left panel: Consecutive differences of $a$-coefficients when
considering increasing $s_{\mathrm{max}}$ for non-zero $w_{s}$. To illustrate
further, in (B), the red block on the lower left shows the quantity
$\log_{10}\left\\{[a_{2}(s_{\mathrm{max}}=3)-a_{2}(s_{\mathrm{max}}=1)]/\sigma(a_{2})\right\\}$.
These consecutive differences for the odd and even $a$-coefficients are shown
in the top and bottom rows, respectively. All blocks in the area to the right
of the yellow line demarcate statistically insignificant consecutive
differences. Right panel: Same as the left panel, but for an increasing number
of neighbouring multiplets $\Delta\ell_{\mathrm{max}}$ instead of increasing
maximum-angular-degree $s_{\mathrm{max}}$. Forward calculations were done for
the multiplet ${}_{0}\mathrm{S}_{280}\,$ using inverted rotation profiles
available on JSOC, which uses HMI data measured during the 72-day solar minima
period between 2019-03-14 and 2019-05-25. Since this figure demonstrates the
saturation of $a$-coefficients and odd $a$-coefficients $a_{j}$ are
significant only when $j\geq s_{\mathrm{max}}$, we hatch the upper-triangle in
(A).

Coupling of modes is permissible within a certain frequency window
$\Delta\omega$. Further, the Wigner $3-j$ symbols in Eqn. (4) impose the
selection rule $|\ell-\ell_{0}|\leq s$, which disallows a central multiplet
with angular degree $\ell_{0}$ from coupling with modes outside the window
$\ell^{\prime}\in[\ell_{0}-s,\,\ell_{0}+s]$. This means that submatrices
corresponding to coupling of angular degrees $(\ell_{0},\ell^{\prime})$ would
be non-zero. Consequently, with an increase in the maximum angular degree of
rotation $s_{\mathrm{max}}$, farther and farther bands in the supermatrix
${\sf Z}$ are filled. The total number of submatrices that constitutes a
supermatrix, depends on the angular degree of the multiplet farthest from the
central multiplet. We define this maximum offset in angular degree between the
central multiplet and its farthest neighbour as
$\Delta\ell_{\mathrm{max}}=\mathrm{max}(\ell-\ell_{0}$). Although at first it
might seem that $s_{\mathrm{max}}=\Delta\ell_{\mathrm{max}}$, that is not
true. $\Delta\ell_{\mathrm{max}}$, infact, controls the frequency window
$\Delta\omega$ and in order for the eigenfrequencies to converge to a stable
value, it is essential to choose a large enough $\Delta\omega$ and hence a
large enough $\Delta\ell_{\mathrm{max}}$. The supermatrix would therefore have
non-zero bands upto the submatrices where $\ell\leq\ell_{0}\pm s$ and have
zeros in all submatrices thereafter. Since each multiplet contains $2\ell+1$
modes, large values of $s$ or $\Delta\ell_{\mathrm{max}}$ correspond to large
sizes of supermatrix ${\sf Z}$. The non-linear inversion involves solving
these eigenvalue problems for all coupled multiplets at every iteration.
Consequently, with increasing size of ${\sf Z}$, the computational cost of
inversions grows significantly. In principle, choosing $s_{\mathrm{max}}$ and
$\Delta\ell_{\mathrm{max}}$ ensures all components of differential rotation
has been accounted for as well as coupling with all modes have been
considered. Therefore, the “true” estimate of the perturbed eigenstates
require $s_{\mathrm{max}},\,\Delta\ell_{\mathrm{max}}\to\infty$. However, this
is (A) not computationally tractable, and (B) not practically necessary since
the eigenfrequencies and eigenfunctions converge to a stable solution with
increasing values of both $s_{\mathrm{max}}$ and $\Delta\ell_{\mathrm{max}}$.
This necessitates a judicious choice of truncation of the maximum angular
degree $s_{\mathrm{max}}$ of $w_{s}(r)$ and the farthest neighbours in $\ell$
that need to be considered for the eigenvalue problem to converge
sufficiently. Since the maximum coupling is expected where
$(\partial\omega/\partial\ell)_{n}$ is the smallest along a radial branch $n$,
we choose $\ell_{0}=281$ on the $f$-branch (see Fig. 1 in Gizon & Birch,
2005). Since we use modes upto $\ell=300$ on the radial branch with
$\Delta\ell_{\mathrm{max}}=19$, ${}_{0}\mathrm{S}_{281}\,$ is the maximally-
coupled central multiplet that could be analyzed.

We test the optimal $s_{\mathrm{max}}$ for differential rotation via forward
calculations using 2D-RLS profiles from JSOC during the solar minimum
corresponding to the 72-day period between 2019-03-14 and 2019-05-25. To do
this, we first vary $s_{\mathrm{max}}$ from 1 to 19 holding
$\Delta\ell_{\mathrm{max}}=19$ fixed. Although modes beyond $\ell_{0}\pm
s_{\mathrm{max}}$ do not couple with the central multiplet, this test allows
us to check for convergence in the $a$-coefficients as a function of
$s_{\mathrm{max}}$. The result of this test is summarized in the left panel of
Fig. 2. Each tile in the colour-map represents the quantity
$\log_{10}\left\\{[a_{j}(s_{\mathrm{max}}=n+2)-a_{j}(s_{\mathrm{max}}=n)]/\sigma(a_{j})\right\\}$.
This serves as a measure of the amount of change in the $a$-coefficients when
$s_{\mathrm{max}}$ is increased by 2. The difference is scaled by the observed
uncertainty of the $a$-coefficients to indicate whether the changes are
statistically significant. We have used a diverging colourbar with red patches
indicating non-negligible differences between successive $s_{\mathrm{max}}$
cases, while blue patches indicate negligible differences. Fig. 2(A) shows
that, for odd $a$-coefficients, the consecutive differences are negligible for
all $s_{\mathrm{max}}$. For completeness, we also present the consecutive
differences of the even $a$-coefficients in Fig. 2(B). It may be safely
inferred that from the case “11-9” onwards, the consecutive differences are
negligible. Therefore, considering $w_{s}(r)=0$ for $s>9$ should give us
$a$-coefficients that are not statistically different from solving the problem
with $s_{\mathrm{max}}=19$.

Next, we vary $\Delta\ell_{\mathrm{max}}$ from 1 to 19 holding
$s_{\mathrm{max}}=19$ fixed. This serves as an independent test for how many
neighbouring multiplets we need to consider in order to ensure saturation of
$a$-coefficients. The results are summarized in the right panel of Fig. 2.
Similar to the previous test, the coloured tiles represent the quantity
$\log_{10}\left\\{[a_{j}(\Delta\ell_{\mathrm{max}}=n+2)-a_{j}(\Delta\ell_{\mathrm{max}}=n)]/\sigma(a_{j})\right\\}$.
From Fig. 2(C) & (D), we see that all red patches occur for
$\Delta\ell_{\mathrm{max}}\leq 8$.

##  4 Results

###  4.1 Inverse problem: DPT vs. hybrid

Figure 3: Comparison between differential-rotation profiles $w_{s}(r)$ up to
angular degree $s=9$ obtained from JSOC in black and our DPT inversions in
dashed red lines. The JSOC profiles are labelled “2D RLS” since they were
obtained from projecting the 2D inversions onto a 1.5D basis, as discussed in
Appendix B. Since the modes that we use in hybrid inversions are almost
insensitive below $r_{\mathrm{th}}=0.9R_{\odot}$, we fix $w_{s}(r)$ to the
JSOC values below this depth, as indicated by the vertical blue lines.

Figure 4: Offset in surface differential rotation inferred using hybrid
inversions as compared to DPT inversions spanning 22 years. We have used MDI
measurements from May 1996 to April 2010 and thereafter HMI measurements up to
January 2018. The vertical dashed-white line indicates the separation between
MDI and HMI regions.

Figure 5: Offset in 22 years of differential rotation as a function of depth
for the chosen latitude $15^{\circ}$. As in Fig. 5, the vertical dashed-white
line marks the separation between MDI and HMI regions.

As demonstrated in Section. 3.3, we may choose $s_{\mathrm{max}}=9$ and
$\Delta\ell_{\mathrm{max}}=8$ without incurring any statistically significant
errors due to truncation. The resultant supermatrix $\mbox{${\sf
Z}$}^{n_{0},\ell_{0}}$ contains contributions from perturbations due to
$w_{1}(r)$ through $w_{9}(r)$ which couples neighbouring multiplets in the
spectral window $\ell_{0}-9<\ell<\ell_{0}+9$. As illustrated in Eqn. (13), we
model the data as the diagonal elements of ${\sf Z}$ for the isolated modes.
For coupled modes, the data is modeled via the full eigenvalue solver. As is
customary, we have validated our inversion methodology with synthetic data
(with and without noise) in Appendix D. Having presented this verification of
our hybrid inversion methodology using artificially generated data for a known
rotation profile, we proceed to perform hybrid inversions on real data.

We carry out the inversion using this “hybrid” modeling method and compare the
results against a purely DPT inversion, which is linear and considers all
multiplets as isolated. The former represents an exact model and consequently
an “exact” inversion, while the latter is an approximation that has been
traditionally used in inferring differential-rotation profiles.

In Fig. 3, we show the consistency of our 1.5D DPT inversions as compared to
the 2D RLS inversions available on JSOC, both of which ignore cross-coupling
of modes. The black line is a projection of the rotation profile
$\Omega(r,\theta)$ inferred via 2D RLS onto the 1.5D profile $w_{s}(r)$, as
outlined in Appendix (B). The red line is from our 1.5D inversions. We note
that the two profiles are not expected to match exactly because of differences
such as inversion methodologies, model parameterization and exact choice of
regularization parameters. Despite this, the two profiles agree well both
qualitatively and quantitatively. The plotted $w_{s}(r)$ profiles are from the
72-day MDI $a$-coefficient measurement starting on $12^{\mathrm{th}}$ July,
1996.

Having established the validity of our DPT inversions against traditionally
accepted results, we present the differences between differential rotation
profiles inferred from our hybrid and DPT inversions in Figs. 5 and 5. Fig. 5
shows the difference at the solar surface as a function of latitude over 22
years, in chunks of 72-day periods, between May 1996 and January 2018. For
this, we ignore temporal variations of differential rotation within each
72-day period. We perform separate inversions of $\mathbf{v}_{\mathrm{rot}}$
(see Eqn. [3]) from $a$-coefficient measurements estimated from the
corresponding 72-day time-series. Fig. 5 is the same, but plotted as a
function of depth from the solar surface down to $0.9R_{\odot}$ at latitude
$\lambda=15^{\circ}$. We used MDI measurements up to April 2010 and HMI
measurements thereafter. This is indicated by the vertical dashed-white line
in both figures. Comparing our Figs. 5 and 5 with Figs. 3(C) and 3(D) in
Vorontsov et al. (2002), respectively, we note that the correction in
differential rotation due to mode coupling is approximately three orders of
magnitude smaller than the torsional-oscillation signal.

Despite this small difference, we note a systematic pattern as a function of
time. In Fig. 5, for the MDI years, we see bands of fast rotation around
$20^{\circ}$ latitudes and slow rotation around $40^{\circ}$ latitudes. The
poles show an alternating pattern of 5 years of slow followed by fast
rotation. Similar features can be discerned in the HMI years. However, the
bands of slower rotation, which were at $40^{\circ}$ for MDI, are less
pronounced and seem to have moved to around $60^{\circ}$ for HMI. Moreover,
the sub-polar branches of fast rotation in MDI that start around 1998 and
merge with the polar fast-rotation band around 2002 are missing for the HMI
years (which were expected between 2011 and 2014). This discrepancy could
possibly be attributed to the change in instruments and the lower SNR for
$a_{7}^{n\ell}$ in HMI as compared to MDI measurements. The depth profile in
Fig. 5 broadly shows a layer of faster rotation around
$0.98R_{\odot}<r<1.0R_{\odot}$ and a layer of slower rotation between
$0.93R_{\odot}<r<0.97R_{\odot}$. Once again, there are differences between the
MDI and HMI years, which may be (at least in part) attributed to systematic
differences in measurements. Although these demonstrate that there are
systematic differences that may be spotted upon careful scrutiny, these are
minuscule and may be disregarded for practical purposes.

###  4.2 Effect on even $a$-coefficients

Figure 6: Same as Figures 1(a) and (b) in Vorontsov (2011) but using a
numerical eigenvalue solver to estimate mode coupling. For each subplot, the
black-solid line shows the contribution from mode coupling (QDPT) and the
black-dashed line marks the contribution from centrifugal effects. The red
line denotes the combined contribution from mode coupling and centrifugal
effects. $a_{2}$ and $a_{4}$ HMI measurements (from the 360-day period between
2010-04-30 and 2011-04-25, which coincides with the solar minimum) are marked
as black dots. The pink lines denote second order effects as presented in
Antia et al. (2000). For the pink lines, we have used the same limits of
angular degree $\ell$ in our plots as used by Antia et al. (2000).

Figure 7: Same as Fig. 7 but for the 360-day period of solar maxima between
2014-04-09 and 2015-04-04.

V11 developed an asymptotic description of mode coupling for high angular
degrees and presented its effects on even $a$-coefficients. In our study,
solving the eigenvalue problem not only gives us the odd $a$-coefficients, but
also numerical estimates of even $a$-coefficients. Therefore, for the sake of
completeness, we also tally our results, which were computed using a numerical
eigenvalue solver, with those of V11. To elaborate further on how we obtained
our estimates of $a_{2}^{n\ell}$ and $a_{4}^{n\ell}$, the following steps may
be considered: (A) using Eqn. (10), we construct the supermatrix $\mbox{${\sf
Z}$}_{k^{\prime}k}$ using the differential rotation profile for the chosen
period (solar minima or maxima as discussed below), (B) since we are
interested in considering mode-coupling, we carry out an eigenvalue problem to
estimate the frequency splittings $\delta{}_{n}\omega_{\ell m}$, and (C) using
Eqn. (7), we estimate the even $a$-coefficients for $j=2,4$.

Figs. 1(a) & (b) in V11 compare contributions of mode coupling (using a semi-
analytic treatment) and centrifugal effects on even coefficients $a_{2}$ and
$a_{4}$ for the $f$-mode. Reproducing the calculations of V11, the distortion
of solar surface due to centrifugal effects and the resulting change in
gravitational potential are written as

$R(\theta)=R_{\odot}\left[1+\epsilon_{2}P_{2}(\cos\theta)+\epsilon_{4}P_{4}(\cos\theta)\right]$
(17)
$\Psi=-\frac{GM_{\odot}}{r}\left[1-\left(\frac{R_{\odot}}{r}\right)^{3}J_{2}P_{2}(\cos\theta)-\left(\frac{R_{\odot}}{r}\right)^{5}J_{4}P_{4}(\cos\theta)\right],$
(18)

where $\epsilon_{2},\,\epsilon_{4}$ are the oblateness coefficients,
$J_{2},\,J_{4}$ are gravitational moments and $P_{2},\,P_{4}$ are Legendre
polynomials. V11 showed that, the even-ordered $a-$coefficients can be written
as

$\delta
a_{2}=\left[\frac{5}{2}\epsilon_{2}+3J_{2}-\frac{10}{3}(\epsilon_{4}+J_{4})\right]\left(\frac{\partial\omega}{\partial\ell}\right)_{n}$
(19) $\delta
a_{4}=\frac{3}{8}\left(7\epsilon_{4}+10J_{4}\right)\left(\frac{\partial\omega}{\partial\ell}\right)_{n}$
(20)

Using the measurements of $J_{2},\,J_{4}$ from Roxburgh (2001), V11 showed
that the corrections to the even-ordered $a-$coefficients can be written as

$\delta a_{2}=-1.20\times
10^{-5}\left(\frac{\partial\omega}{\partial\ell}\right)_{n},\qquad\delta
a_{4}=-1.47\times
10^{-6}\left(\frac{\partial\omega}{\partial\ell}\right)_{n}.$ (21)

From observed data, $\partial\omega/\partial\ell$ can be computed to be
${}_{n}\omega_{\ell+1}-{}_{n}\omega_{\ell}$. These calculations are shown in
Fig. 7. V11 used one year of SOHO/MDI measurements over the solar minimum to
diminish contributions from magnetism. We carry out the same analysis using
results from an eigenvalue solver, but using measurements from SDO/HMI during
a 360-day period over a solar minimum (2010-04-30 to 2011-04-25). From Fig. 7
we conclude that, in agreement with V11, mode coupling and centrifugal effects
together provide an adequate fit to the observed $a_{2}$ coefficients. Mode
coupling also overwhelmingly dominates $a_{4}$ measurements, providing a good
fit to the data. However, for $a_{2}$, V11 observed mode coupling to dominate
over centrifugal effects beyond $\ell=200$, whereas we find this to happen
beyond $\ell\sim 170$. We also plot the second-order effect of rotation as in
Antia et al. (2000) in the dashed magenta line and the combined contribution
from mode-coupling and these second order effects in the solid magneta line.
The second-order correction in Antia et al. (2000) is a more complete
treatment that the asymptotic calculation of centrifugal effects in V11. For
both $a_{2}^{n\ell}$ and $a_{4}^{n\ell}$, just mode-coupling added with V11’s
asymptotic centrifugal effects seems to explain the HMI data at solar minima,
accounting for a more rigorous second-order correction due to rotation clearly
grazes the outer envelop of the black dots. This suggests that there are other
structure perturbations that need to be accounted for such as sound-speed
anomaly or a weak background solar minima magnetic field.

Fig. 7 is plotted in the same spirit as Fig. 7 but for the 360-day solar
minima period between 2014-04-09 and 2015-04-04. Neither $a_{2}^{n\ell}$ nor
$a_{4}^{n\ell}$ are explain by the effects from rotation alone. Since solar
maxima has significantly more magnetic activity, this departure in the plotted
even $a$-coefficients could be attributed to strengthened magnetic fields.
According to selection rules of $a$-coefficient kernels for Lorentz-stresses
(see Das et al., 2020), these departures in even $a$-coefficients may be using
to infer the following magnetic quantities:
$B_{r}^{2},(B_{\theta}^{2}+B_{\phi}^{2}),(B_{\theta}^{2}-B_{\phi}^{2})$ and
$B_{r}\,B_{\theta}$.

Figure 8: Departure of QDPT-estimated $a_{2}^{0\ell}$ and $a_{4}^{0\ell}$
coefficients from the DPT estimates. Differences are scaled with respect to
the standard deviations of the respective measurements. Statistically
significant differences of magnitudes $>1\sigma$ lie outside the black hatched
area. Forward calculations were performed using inverted rotation profiles
available on JSOC, which are associated with HMI data measured during the
72-day period between 2014-04-09 and 2014-06-20 which corresponds to a solar
maxima.

Measurements during solar maxima are expected to possess stronger signatures
of magnetic fields as compared to the solar minima. Therefore, we also provide
an estimate of the even $a$-coefficients (scaled by the measured 360-day
uncertainties) using JSOC rotation profiles from the solar maximum between
2014-04-09 and 2015-04-09. We carried out full eigenvalue solutions using
$s_{\mathrm{max}}=19$ for all the modes observed by HMI. Only the $f$-modes
were found to have predicted $a_{2}$ and $a_{4}$ which were larger than the
uncertainties from a one-year period. The result for the modes resolved by HMI
for $n=0$ is shown in Fig. 8. While all measured modes have
$a_{4}/\sigma(a_{4})>1$, only $\ell>160$ modes have $a_{2}/\sigma(a_{2})>1$.
The SNR for $a_{4}$ is about twice as large as $a_{2}$. If the same
calculation is repeated for 360-day measurements instead of the currently used
72-day measurements, the SNR would be scaled up by $\sqrt{5}$.

##  5 Discussion

Measurement of differential rotation is one of the triumphs of
helioseismology. Traditional methods of inference have used frequency
splittings $\delta\omega_{n\ell}^{m}$ for global modes of oscillations under
the assumption that modes are self-coupled. While this is acceptable for most
modes in observed power spectra, this approximation worsens for high angular
degrees in the $f$ and $p_{1}$ branches. Kashyap et al. (2021) investigated
the applicability of the isolated-multiplet approximation across all the modes
resolved by HMI up to $\ell=300$. They found deviations between
$\delta\omega^{m,\mathrm{QDPT}}_{n\ell}$ and
$\delta\omega^{m,\mathrm{DPT}}_{n\ell}$ beyond $\ell=160$ of the order of
$\approx\pm 2\,\sigma(\delta\omega_{n\ell})$, where
$\sigma(\delta\omega_{n\ell})$ are the uncertainties in the measured frequency
splittings. The frequency splittings may be expressed in terms of odd and even
$a$-coefficients, which constrain the equatorially symmetric part of
differential rotation and axisymmetric structure perturbations, respectively.
This motivates the current study which is aimed at: (i) investigating the
changes in predicted odd and even $a$-coefficients via an eigenvalue treatment
of the operator ${\sf Z}$, (ii) finding corrections to differential rotation
when considering mode coupling as compared to isolated multiplets, and (iii)
comparing the expected second-order changes in the even $a$-coefficients
computed from the exact eigenvalue solver as compared to asymptotic analysis
from previous studies, such as V11.

In this study, we adopt the QDPT approach to modeling the frequency splittings
due to differential rotation. To the best of our knowledge, all studies that
use $a$-coefficients to infer differential rotation assume that multiplets are
isolated. We carry out forward modeling as well as an inverse problem by
constructing the full supermatrix ${\sf Z}$ that accounts for coupling, with
all neighbouring multiplets adhering to selections rules imposed by Wigner
$3-j$ symbols. Section 8(c) in LR92 discusses the vanishingly small effect of
QDPT for low angular-degree modes, quantified by the coupling-strength
coefficient. Using an asymptotic analytical treatment, V11 proposed the
maximum correction in frequency splittings as a function of $\ell$ and
mentioned that the effect of mode coupling on oscillation frequencies might be
observable towards high $\ell$. Our study demonstrates that the mode-coupling-
induced corrections in odd $a$-coefficients introduce a relatively tiny change
in differential rotation or torsional oscillations. We present this in Section
4.1 by performing 1.5D inversions for 22 years of MDI and HMI data. The
corrections show a systematic variation over solar cycles, although the
correction amplitudes are minuscule. This implies that accounting for mode
coupling when inferring differential rotation solely from odd $a$-coefficients
is insignificant. In Section 4.2, we also compare the results of our predicted
even $a$-coefficients arising from mode coupling with those presented in the
semi-analytic treatment of V11 and find them to be consistent with only very
minor differences. Therefore, we believe that future studies accounting for
even $a$-coefficients due to mode coupling can reliably use the corrections
suggested by V11 instead of the more computationally expensive method
performed here.

Apart from the correction in torsional oscillation (which is weak but has a
systematic variation over solar cycles), we confirm that considering the
contribution of differential rotation to even splitting coefficients due to
mode coupling is statistically significant. Fig. (8) shows that this effect in
$a_{2}^{n\ell}$ and $a_{4}^{n\ell}$ is as large as $1.5\,\sigma$ and
$2.5\,\sigma$, respectively. The effect of mode coupling on $a_{2}^{n\ell}$
coefficients was accounted for in a comprehensive study by Chatterjee & Antia
(2009) to put limits on flow velocities in giant cells. However, most studies
that use even $a$-coefficients to infer global magnetism have not accounted
for this correction. For instance, Antia et al. (2000) calculated estimates of
magnetic field using forward modelling, after considering distortion due to
second-order corrections of rotation on even $a$-coefficients. They predicted
a magnetic field of strength 20 kG (or an equivalent acoustic perturbation)
located 30 Mm below the solar surface using even coefficients $a_{2}^{n\ell}$
and $a_{4}^{n\ell}$. Baldner et al. (2009, 2010) too carried out inversions
for magnetic field strengths without accounting for mode coupling. Therefore,
improving such prior estimates when the effect of QDPT is considered remains
an important area in which to make progress.

S.B.D and S.G.K thank Jesper Schou (Max Planck Institute for Solar System
Research) for numerous insights and Pritwiraj Moulik (Princeton University)
for kindly providing the B-splines package as well as multiple discussions on
inversion methodology.

## References

  * Antia et al. (2001) Antia, H. M., Basu, S., Hill, F., et al. 2001, MNRAS, 327, 1029, doi: 10.1046/j.1365-8711.2001.04819.x
  * Antia et al. (2000) Antia, H. M., Chitre, S. M., & Thompson, M. J. 2000, A&A, 360, 335
  * Baldner et al. (2009) Baldner, C. S., Antia, H. M., Basu, S., & Larson, T. P. 2009, ApJ, 705, 1704, doi: 10.1088/0004-637X/705/2/1704
  * Baldner et al. (2010) —. 2010, Astronomische Nachrichten, 331, 879, doi: 10.1002/asna.201011418
  * Baldner & Basu (2008) Baldner, C. S., & Basu, S. 2008, ApJ, 686, 1349, doi: 10.1086/591514
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs, 0.3.13.  http://github.com/google/jax
  * Brown et al. (1989) Brown, T. M., Christensen-Dalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526, doi: 10.1086/167727
  * Chatterjee & Antia (2009) Chatterjee, P., & Antia, H. M. 2009, ApJ, 707, 208, doi: 10.1088/0004-637X/707/1/208
  * Christensen-Dalsgaard et al. (1996) Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286
  * Das et al. (2020) Das, S. B., Chakraborty, T., Hanasoge, S. M., & Tromp, J. 2020, ApJ, 897, 38, doi: 10.3847/1538-4357/ab8e3a
  * Gizon & Birch (2005) Gizon, L., & Birch, A. C. 2005, Living Reviews in Solar Physics, 2, 6
  * Gough & Thompson (1990) Gough, D. O., & Thompson, M. J. 1990, MNRAS, 242, 25, doi: 10.1093/mnras/242.1.25
  * Howe (2008) Howe, R. 2008, Advances in Space Research, 41, 846, doi: 10.1016/j.asr.2006.12.033
  * Kashyap et al. (2021) Kashyap, S. G., Das, S. B., Hanasoge, S. M., Woodard, M. F., & Tromp, J. 2021, ApJS, 253, 47, doi: 10.3847/1538-4365/abdf5e
  * Kosovichev et al. (1997) Kosovichev, A. G., Schou, J., Scherrer, P. H., et al. 1997, Sol. Phys., 170, 43, doi: 10.1023/A:1004949311268
  * Larson & Schou (2018) Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29, doi: 10.1007/s11207-017-1201-5
  * Lavely & Ritzwoller (1992) Lavely, E. M., & Ritzwoller, M. H. 1992, Philosophical Transactions of the Royal Society of London Series A, 339, 431, doi: 10.1098/rsta.1992.0048
  * Libbrecht (1989) Libbrecht, K. G. 1989, ApJ, 336, 1092, doi: 10.1086/167079
  * Ritzwoller & Lavely (1991) Ritzwoller, M. H., & Lavely, E. M. 1991, ApJ, 369, 557, doi: 10.1086/169785
  * Roxburgh (2001) Roxburgh, I. W. 2001, A&A, 377, 688, doi: 10.1051/0004-6361:20011104
  * Schad & Roth (2020) Schad, A., & Roth, M. 2020, ApJ, 890, 32, doi: 10.3847/1538-4357/ab65ec
  * Scherrer et al. (1995) Scherrer, P. H., Bogart, R. S., Bush, R. I., et al. 1995, Sol. Phys., 162, 129, doi: 10.1007/BF00733429
  * Schou (1999) Schou, J. 1999, Index of / schou/anavw72z, http://quake.stanford.edu/~schou/anavw72z/?C=M;O=A
  * Schou et al. (1994) Schou, J., Christensen-Dalsgaard, J., & Thompson, M. J. 1994, ApJ, 433, 389, doi: 10.1086/174653
  * Schou et al. (1998) Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390, doi: 10.1086/306146
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Spiegel & Zahn (1992) Spiegel, E. A., & Zahn, J. P. 1992, A&A, 265, 106
  * Tarantola (1987) Tarantola, A. 1987, Inverse problem theory. Methods for data fitting and model parameter estimation
  * Thompson et al. (1996) Thompson, M. J., Toomre, J., Anderson, E. R., et al. 1996, Science, 272, 1300, doi: 10.1126/science.272.5266.1300
  * Vorontsov (2011) Vorontsov, S. V. 2011, Monthly Notices of the Royal Astronomical Society, 418, 1146, doi: 10.1111/j.1365-2966.2011.19564.x
  * Vorontsov et al. (2002) Vorontsov, S. V., Christensen-Dalsgaard, J., Schou, J., Strakhov, V. N., & Thompson, M. J. 2002, Science, 296, 101, doi: 10.1126/science.1069190
  * Woodard et al. (2013) Woodard, M., Schou, J., Birch, A. C., & Larson, T. P. 2013, Sol. Phys., 287, 129, doi: 10.1007/s11207-012-0075-9
  * Woodard (2016) Woodard, M. F. 2016, MNRAS, 460, 3292, doi: 10.1093/mnras/stw1223

##  Appendix A Noisy $s=7$ HMI measurements

Figure 9: Normalized histograms for signal-to-noise of odd $a$-coefficients
$a_{1}-a_{9}$. The gray, blue and red colors indicate modes with low,
intermediate and high angular degree $\ell$. Histograms in the top panel were
constructed from MDI data measured during the 72-day period between 1996-07-12
and 1996-09-22. Histograms in the bottom panel were constructed from HMI data
measured during the 72-day period between 2010-04-30 and 2010-07-11. Both of
these periods coincides with solar minima.

Figure 10: 1.5D counterparts of the 2D-RLS rotation profiles obtain from JSOC
during two different solar minima phases. The profiles in the top panel
correspond to MDI measurements between 1996-07-12 and 1996-09-22 and those in
the bottom panel correspond to HMI measurements between 2010-04-30 and
2010-07-11.

Inference of rotation as a function of depth depends critically on the signal-
to-noise (SNR) of the measured $a$-coefficients of modes sensitive to
perturbations at that depth. High (low) angular-degree modes have shallow
(deep) lower turning points and are sensitive to shallow (deep) perturbations.
Therefore, it is necessary to measure the SNR of the coupled modes. In doing
so, we found that some of the 72-day HMI measurements have poor SNR for
$a_{7}^{n\ell}$. This is reflected in the considerably wiggly $w_{7}(r)$ found
from decomposing the 2D-RLS results on JSOC into 1.5D rotation profiles. Fig.
9 shows normalized histograms for $a$-coefficients measured during phases of
solar minima by MDI and HMI in three ranges of angular degrees: (a) low $\ell$
modes in gray where $0\leq\ell<100$, (b) intermediate $\ell$ modes in blue
where $100\leq\ell<200$, and, (c) high $\ell$ modes in red where
$200\leq\ell<300$. For $a_{1}^{n\ell},a_{3}^{n\ell},a_{5}^{n\ell}$ and
$a_{9}^{n\ell}$, the SNR for the shallow-sensitive modes is significantly
larger than unity for both MDI and HMI. The SNR for only the HMI-measured
$a_{7}^{n\ell}$, however, is largely contained within $\pm 1$. Consequently,
as seen in Fig. 10, the inferred $w_{s}(r)$ profiles from MDI have the desired
smoothness imposed by the regularization term. This is also true for
$w_{1}(r),w_{3}(r),w_{5}(r)$ and $w_{9}(r)$ profiles from HMI. However, the
$w_{7}(r)$ profile from HMI is unusually wiggly.

##  Appendix B Converting 2D rotation profiles to 1.5D

In the 1.5D inversions presented in Section (4.1), we fix $w_{s}(r)$ below
$r_{\mathrm{th}}=0.9R_{\odot}$ to the corresponding JSOC profiles (which were
obtained via 2D RLS). This requires converting $\Omega(r,\theta)$ to its 1.5D
equivalent $w_{s}(r)$. We follow the prescription of Ritzwoller & Lavely
(1991) which is outlined here for completeness and ease of reference of the
reader. The rotation profile $\Omega(r,\theta)$ may be written in terms of
Legendre polynomials $P_{k}$ as

$\Omega(r,\theta)=\sum_{k=0,2,4...}\Omega_{k}(r)\,P_{k}(\cos\theta)$ (B1)

Further, we know that rotational velocity
${\mathbf{v}}_{\mathrm{rot}}(r,\theta)=\hat{\mathbf{r}}\times\hat{\mathbf{z}}\,\Omega(r,\theta)$.
Using this and Eqn. (3), we have

$\mathrm{v}_{\mathrm{rot}}=r\sin\theta\sum_{k=0,2,4...}\Omega_{k}(r)\,P_{k}(\cos\theta)=-\sum_{s=1,3,5,...}\,w_{s}(r)\,\partial_{\theta}Y_{s0}$
(B2)

Since our equations are written in terms of $w_{s}(r)$, we need to convert
$\Omega_{k}(r)\to w_{s}(r)$. To do this, we may project
$\mathrm{v}_{\mathrm{rot}}$ onto the basis of $P_{k}$, i.e.,

$r\,\Omega_{k}(r)=-\frac{2}{2k+1}\int\frac{1}{\sin\theta}\left(\sum_{s=1,3,5,...}w_{s}(r)\,\partial_{\theta}Y_{s0}(\theta)\right)P_{k}(\cos\theta)\sin\theta\,{\mathrm{d}}\theta\,,$
(B3)

which gives us the following matrix equation

$r\,\Omega_{k}(r)=\sum_{s}\alpha_{ks}\,w_{s}(r)$ (B4)

Finally, $w_{s}(r)$ can be obtained by
$w_{s}(r)=\bm{\alpha}^{-1}\,\Omega_{k}(r)$. To obtain upto $w_{9}(r)$, we use
the following matrix inverse

$\bm{\alpha}^{-1}=2\sqrt{\pi}\left(\begin{array}[]{cccccc}1/\sqrt{3}&-1/5\sqrt{3}&0&0&0&0\vspace{0.2cm}\\\
0&1/5\sqrt{7}&-1/9\sqrt{7}&0&0&0\vspace{0.2cm}\\\
0&0&1/9\sqrt{11}&-1/13\sqrt{11}&0&0\vspace{0.2cm}\\\
0&0&0&1/13\sqrt{15}&-1/17\sqrt{15}&0\vspace{0.2cm}\\\
0&0&0&0&1/17\sqrt{19}&-1/21\sqrt{19}\vspace{0.2cm}\\\
0&0&0&0&0&1/21\sqrt{23}\\\ \end{array}\right)$ (B5)

##  Appendix C Accuracy and robustness of the numerical eigenvalue solver

Figure 11: Ratio of the L2-norm of $E$, the difference between the hybrid and
DPT supermatrices, and the L2-norm of the DPT supermatrix $A$. The blue, red
and black dashed lines correspond to angular degrees $\ell=200,240,280$,
respectively.

We use the linalg.eigh module in Python’s numpy package for computing
eigenvalues. The lialg.eigh module is a wrapper for its LAPACK implementation
of evaluating eigenvalues of real symmetric or complex Hermitian matrices. As
mentioned in the documentation webpage of numpy.linalg.eigh, for a symmetric
and real matrix (as is the case for our supermatrix), the _syevd routine is
used for solving a real symmetric matrix using divide and conquer algorithm.
Under “Application Notes” in the _syevd webpage, the developers mention that
“The computed eigenvalues and eigenvectors are exact for a matrix A+E such
that $||E||_{2}=\mathcal{O}(\epsilon)\,||A||_{2}$, where $\epsilon$ is the
machine precision”. In our case, $A$ may be regarded as the DPT supermatrix
and $E$ would then be the difference between the hybrid and DPT supermatrices.
Note that the hybrid supermatrices are constructed from the final converged
$c_{s}^{p,\mathrm{fit}}$ arrays (see Eqn. [12]) obtained from our inversion
algorithm — the coefficients we use to construct the final plots Figs. 5 & 5.

For our problem, we compute the DPT supermatrix (which ignores coupling across
multiplets) and the hybrid supermatrix (which accounts for coupling across
multiplets). Since the coupling across multiplets induced due to differential
rotation is weak — as evidenced by the smallness of our corrections, it is
reasonable to ask if the matrices themselves are different enough for LAPACK’s
_syevd algorithm to yield distinctly different eigenvalues. We can frame this
question, in light of the Application Note provided by the LAPACK developers,
as: “Is the ratio of L2-norm of our $E$ matrix (difference between hybrid and
DPT supermatrix) and the L2-norm of our DPT matrix $A$ significantly larger
than the machine precision?” We carry out our calculations in 64-bits
precision, meaning $\mathcal{O}(\epsilon)\approx 10^{-16}$. So, to investigate
the above question, we have calculated $||E||_{2}\big{/}||A||_{2}$ across
multiple years for three large angular degrees $(\ell=200,240,280)$ where
hybrid fitting is applicable. The results are presented in Fig. 11. We see
that the ratio is consistently larger than $\mathcal{O}(10^{-6})$ which is
atleast $10^{10}$ times larger than the order of machine precision
$\mathcal{O}(\epsilon)$. Therefore, our DPT supermatrix $A$ and the
corresponding hybrid supermatrix $A+E$ are different enough to yield
distinctly different eigenvalues. This shows that the eigenvalue routine used
in Eqn. (14) is accurate enough to yield eigenvalues which are not garbled by
machine errors for the level of difference between DPT and hybrid
supermatrices.

Eqn. (15) shows the cost function and the inversion involves fitting the model
parameters $c_{s}^{p,\mathrm{fit}}$ while trying to minimize this cost
function. Since our problem involves an eigenvalue operation, this is non-
linear and so we adopt the standard Newton’s method (Tarantola, 1987). This
involves starting from a guess solution (which is close enough to the true
solution than minimizes the cost function) and iteratively stepping towards
the minima by computing a gradient vector and Hessian tensor at each updated
model parameter. These are commonplace in machine learning community and
Google’s jax package has emerged as an efficient tool for the same (Bradbury
et al., 2018). In our study, we have used the automatic differentiation
routines jax.grad to compute gradients and the routines jax.jacfwd and
jax.jacrev to compute the Hessian. Note that jax uses 32-bits precision by
default but we have used the additional switch config.update(’jax_enable_x64’,
True) to use 64-bit machine precision in our calculations.

##  Appendix D Validation of inversion

This section demonstrates the validity of the non-linear hybrid inversions
using synthetic, yet realistic rotation profiles. To do so, we choose a
reference profile $\mathbf{v}_{\mathrm{rot}}^{\mathrm{ref}}(r,\theta)$
available through the JSOC pipeline corresponding to the 72-day MDI
$a$-coefficient measurements starting on $12^{\mathrm{th}}$ July, 1996. We
then use Eqns. (9), (12) & (13) to generate $a$-coefficients using the hybrid
method, i.e., isolated multiplet treatment for DPT modes and full QDPT
treatment for coupled modes. We first carry out an inversion using these
noise-free synthetically generated $a$-coefficients. Fig. 12 compares the
spherical harmonic components of the inverted profile
$\mathbf{v}_{\mathrm{rot}}^{\mathrm{inv}}$ with those of the synthetic profile
$\mathbf{v}_{\mathrm{rot}}^{\mathrm{ref}}$ (see Eqn. [3] for the spherical
harmonic decomposition). The two profiles are seen to be exactly the same
indicating the validity of noise-free inversion methodology. Similarly,
inversions using data corrupted with synthetic noise (at 0.1 $\sigma$, where
$\sigma$ represents the level of noise from observed $a$-coefficients) is
shown in Fig. 14. In this case, the $\mathbf{v}_{\mathrm{rot}}^{\mathrm{inv}}$
and $\mathbf{v}_{\mathrm{rot}}^{\mathrm{ref}}$ are within the expected errors
with neither large fluctuations nor overly smoothened profiles (both of which
are usually seen for improperly regularized inversions in the presence of
noisy data). Hence, we deem the inversions to be successfully benchmarked
using realistic rotation profiles. We have performed this 2-pronged test
(first clean and then noisy inversion) over a variety of different profiles to
validate the robustness of our inversion methodology.

Figure 12: Validation of noise-free non-linear inversion. We construct
synthetic dataset with $a-$coefficients generated using Eqn. (13), shown by
the solid black curve. We carry out non-linear inversions for the $w_{s}(r)$
profiles using these synthetically generated $a-$coefficients. The inverted
profiles are shown by the dashed red lines.

Figure 13: Validation of inversion with noise. As in Fig. 12, the solid black
line indicates the synthetic $w_{s}(r)$ profiles to which we add noise before
carrying out a non-linear inversion to infer the dashed red line.

Figure 14: The first five panels show L-curves for the different angular
degrees $s=\\{1,3,5,7,9\\}$. These plots are representative of L-curves
corresponding to a typical 72-day MDI or HMI inversion. The red marker in each
plot shows our chosen knee of the curve and the corresponding value of
$\mu_{s}$ is mentioned in the title of each subplot. The last panel in the
lower-right, shows a histogram (normalized to yield unit area) of
$\chi^{2}_{d}/N$ for all MDI inversions.

The first five panels in Fig. 14 show representative L-curves for determining
regularization parameters $\mu_{s}$ corresponding to the different angular
degrees $s=\\{1,3,5,7,9\\}$. The red marker indicates the chosen optimal value
of $\mu_{s}$ located at the knee of the curve. This choice of regularization
parameter represents an optimal balance between the data misfit (the degree to
which theoretical predictions from our inferred profile matches the observed
data) and the model smoothness (the degree of smoothness we impose in order to
generate physically meaningful solutions and avoid unrealistic oscillatory
behaviour of inferred profiles). The lower-right panel in Fig. 14 shows a
distribution of the data misfit (scaled by the total number of data points
$N$) where the histogram is normalized to yield unit area. The distribution
peaks around a scaled chi-square value of 1.4 which validates that our
inversion is not over-fitting while producing inferred profiles which generate
predictions having high-fidelity to the observed $a$-coefficients. The
histogram is created using chi-square values from all MDI inversions.

Figure 15: Histograms showing the change in $\chi^{2}$ value between
successive iterations across all 72-day MDI inversions. Histograms are color-
coded according to their iteration number, as shown in the figure legend. Each
histogram is normalized to unit area.

Fig. 15 is presented to demonstrate the convergence of our iterative non-
linear Newton inversion. Each histogram shows the distribution of the change
in total chi-square value between successive iterations. The histogram
corresponding to the first iteration shows that, on an average, the total
misfit drops by $\mathcal{O}(10^{3})$ in the first iteration. Thereafter, the
following iterations are seen to cluster around a much smaller change in
misfit — around 8 orders of magnitude smaller than the change in the first
iteration. For the hybrid inversions, we start from the DPT profile which is
obtained from carrying out a linear inversion under the isolated multiplet
approximation. Therefore, Fig. 15 shows that during the iterative inversions
the rotation profiles undergo an initial non-trivial change from the DPT
profiles, followed by negligible or insignificant changes in the following
iterations. We run each of our inversions for five iterations to ensure the
above convergence in misfit is achieved for all the cases.