© The Authors 2026

1. Introduction

Internal differential rotation is crucial for understanding the Sun’s dynamo and angular momentum transport mechanisms (e.g., Käpylä et al. 2023). Our current knowledge of the Sun’s differential rotation comes primarily from p-mode helioseismology, which exploits the frequency splittings of modes with different azimuthal orders (e.g., Thompson et al. 2003). The Sun’s convection zone (CZ), which constitutes the outer 30% of the Sun, rotates differentially, varying with both radius and latitude, while the radiative interior rotates almost rigidly. The tachocline exists at the transition between these two zones (e.g., Spiegel & Zahn 1992), and there is also a near-surface shear layer in the outermost 5% of the solar radius (e.g., Thompson et al. 1996).

Despite remarkable achievements, p-mode helioseismology still faces significant limitations. It provides limited information about the Sun’s rotation below ∼0.2 R_⊙ (e.g., Gizon et al. 1997) and in regions within ∼20° of the rotation axis (e.g., Schou et al. 1998). This limitation arises because very few p-modes have appreciable kinetic energy density in these regions, as their sensitivity kernels peak near the solar surface. As a consequence, inversions of frequency splittings obtained using various methods can differ significantly in the polar regions (Thompson et al. 2003, their Figure 7).

Solar inertial modes are low-frequency global modes of oscillation that are restored by the Coriolis force. These are retrograde modes with a strong toroidal component and periods on timescales comparable to the solar rotation period (∼25 days at the equator). Maps of horizontal flows near the solar surface have been measured for more than a decade using p-mode local helioseismology and local correlation tracking applied to data from the Helioseismic and Magnetic Imager aboard the Solar Dynamics Observatory (HMI/SDO) and Global Oscillations Network Group (GONG). These long-duration observations allow for the unambiguous detection and characterization of solar inertial modes (Löptien et al. 2018; Gizon et al. 2021; Hanson et al. 2022).

In this work we focus on the high-latitude mode with azimuthal wavenumber m = 1 and north-south symmetric radial vorticity (hereafter the HL1 mode). This mode has been observed to have the largest amplitude on the surface of the Sun, in the range 5–25 m/s at latitudes above 60° (Gizon et al. 2021). It has been consistently detected in multiple datasets and is also present in the Mount Wilson Observatory Dopplergrams spanning the last five sunspot cycles (Ulrich 2001; Liang & Gizon 2025). The frequency of the HL1 mode is known to a high degree of precision. Over the period of HMI observations, its mean frequency in the Carrington frame is ω_obs/2π = −87.9 ± 1.9 nHz (Fig. 1).

Fig. 1. Observed HL1 frequencies in the Carrington reference frame extracted from HMI Dopplergrams at a cadence of one year with a three-year sliding window (blue points with error bars, Liang & Gizon 2025). The mean frequency ωobs/2π = −87.9 nHz is indicated by the black horizontal line, and the gray shaded region shows the ±1σ range (±1.9 nHz). The frequency uncertainty was estimated from measurements over independent time intervals centered on 2012, 2015, 2018, and 2021, thereby accounting for correlations introduced by the sliding window. The solid green line marks the real part of the eigenfrequency computed using eigenvalue Solver 1 (−95.6 nHz), and the dashed green line marks the same for Solver 2 (−96.1 nHz). The azimuthal velocity of the observed mode at the surface is shown as an inset (adapted from Gizon et al. 2021).

Because their physics is intrinsically linked to rotation, inertial modes are particularly well suited for diagnosing differential rotation (see Goddard et al. 2020, for the case of Rossby modes). Linear eigenmode computations have indicated that the kinetic energy density of the HL1 mode peaks inside the convection zone (Gizon et al. 2021; Bekki et al. 2022), making it a very promising candidate for directly probing rotation in the solar interior.

2. Mismatch between model and observed frequencies of HL1 inertial mode

In order to place constraints on the solar internal rotation using the observed frequency of the HL1 mode, we needed to solve the forward problem, that is, to compute the mode frequency for a given prescription of the solar differential rotation (see Appendix A). We denoted the solar angular velocity measured in an inertial frame by Ω(r, θ), where r is the radius and θ is the colatitude (θ = 0 at the north pole). We adopted the reference differential rotation rate Ω_ref(r, θ) derived from p-mode helioseismology (Larson & Schou 2018). In the Carrington frame rotating at angular frequency Ω_Carr/2π = 456.0 nHz, we linearized the equations of momentum, entropy, and continuity about a standard solar model (Model S, Christensen-Dalsgaard et al. 1996). The computational domain extends radially over the range 0.65 R_⊙ < r < 0.985 R_⊙. We assumed perturbations proportional to exp(imϕ − iωt), where m = 1 and ϕ is the longitude. Under horizontal stress-free boundary conditions, the system yields a discrete spectrum of complex eigenvalues ω and corresponding eigenvectors. We used an eigenvalue solver based on Bekki et al. (2022) to compute the eigenmodes in the inertial frequency range (Solver 1; see Appendix A). This solver was validated through a code-to-code benchmark comparison with the independent implementation of Mukhopadhyay et al. (2025) (Solver 2; see Appendix A). The m = 1 high-latitude eigenmode is identified in the spectrum as the only self-excited low-frequency mode (ℑω > 0) that exhibits north-south symmetric radial vorticity. Using Solver 1, we obtained a model eigenfrequency of ω_ref/2π = −95.6 + 3.9 i nHz in the Carrington frame. The real part of ω_ref thus differs from the observed ω_obs by four standard deviations (see Fig. 1). The velocity eigenfunction at the surface, however, is similar to the observation, with power confined to a similar latitudinal range above 60° and a qualitatively similar spiral pattern.

We performed a set of eigenvalue computations that indicate that this discrepancy is likely due to an incorrect reference differential rotation profile. This discrepancy cannot be fixed by varying the other parameters of the background solar model, such as the viscous or thermal diffusivities or the superadiabaticity. In particular, we computed eigenmodes around the commonly used value of turbulent diffusivity in the convection zone ν_CZ = 10¹² cm² s⁻¹. Within the range 10¹¹ cm² s⁻¹ ≤ ν_CZ ≤ 2 × 10¹² cm² s⁻¹, we find that the eigenfrequency remains within the bounds 5.0 nHz < (ω_obs − ℜω)/2π < 10.0 nHz, indicating that the diffusion cannot be used as a tunable parameter to bring the eigenfrequency within observational errorbars. Further, we find that the superadiabaticity (δ = ∇ − ∇_ad) affects the eigenfrequency only weakly, by at most 3 nHz when the superadiabatic gradient is varied in the range |δ|< 5 × 10⁻⁷. This motivates focusing on the background rotation rate as the primary physical parameter of the model that can be adjusted to bring the model into agreement with the observations, while keeping δ = 0 and ν_CZ = 10¹² cm² s⁻¹ fixed.

3. Sensitivity of HL1 mode frequency to the local rotation rate

The linear sensitivity of the mode frequency to a spatial change in rotation is specified through the integral equation

$$\begin{array}{c}\hfill \delta \omega ={\displaystyle {\int }_0^{\pi }{\int }_0^{R_{\odot }}}K(r,\theta )\delta \mathrm{\Omega }(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,\end{array}$$(1)

where δω = ω − ω_ref is the perturbation to the (complex) eigenfrequency that results from a small perturbation in the rotation rate, δΩ = Ω − Ω_ref, and K is the (complex) sensitivity kernel. The discretized version of this equation is

$$\begin{array}{c}\hfill \delta \omega ={\displaystyle \sum _{i,j}}K(r_i,{\theta }_j)\delta \mathrm{\Omega }(r_i,{\theta }_j)\mathrm{d}\mathrm{\Sigma },\end{array}$$(2)

where dΣ = r dr dθ is the surface element in the meridional plane (expressed in terms of the radial and latitudinal grid spacings, dr = 2.3 Mm and dθ = 2.0°) used in our calculations. The most straightforward method for computing the sensitivity kernel is numerical. For a given spatial pixel (r_i, θ_j) where the kernel is to be evaluated, we considered a rotation profile that is perturbed only at that pixel,

$$\begin{array}{c}\hfill \delta \mathrm{\Omega }(r_p,{\theta }_q)=ϵ{\delta }_{\mathit{pi}}{\delta }_{\mathit{qj}},\end{array}$$(3)

where the grid (r_p, θ_q) spans the full numerical domain. Here, δ_pi and δ_qj are Kronecker delta functions. The perturbation amplitude ϵ must be sufficiently small to ensure that the resulting change in oscillation frequency remains in the linear regime. We chose ϵ/2π = 1.0 nHz. The latitudinal entropy gradient was also perturbed numerically to maintain thermal wind balance. Using Solver 1, we computed the eigenfrequency of the high-latitude mode, ω_ϵ(r_i, θ_j), for this perturbed system. The sensitivity kernel at pixel (i, j) was then given by

$$\begin{array}{c}\hfill K(r_i,{\theta }_j)\mathrm{d}\mathrm{\Sigma }=[{\omega }_{ϵ}(r_i,{\theta }_j)-{\omega }_{\mathrm{ref}}]/ϵ.\end{array}$$(4)

The real part of the kernel, presented in Fig. 2, shows a prominent peak centered at (r₀, θ₀) = (0.82 R_⊙, 15° ) and two smaller sidelobes, one positive and one negative, at 0.72 R_⊙. The main peak of the kernel has full widths at half maximum of 7° in latitude and 0.13 R_⊙ in radius, as shown in the cuts in Figs. B.1 and B.2. We found almost no sensitivity in the radiative interior.

Fig. 2. Spatial sensitivity of the frequency of the HL1 mode to localized changes in the background rotation rate. The kernel K(r, θ) is symmetric across the equator. Grid sizes dr = 2.3 Mm and dθ = 2° were used for this calculation. For reference, a patch of area 10dr × 2rdθ = 20dΣ is outlined with a solid black contour. The dashed black contour corresponds to the value of the kernel at half maximum.

Because the kernel is reasonably well localized in both radius and latitude, the inferred frequency shift can be interpreted as providing a direct estimate of the local rotation perturbation at location (r₀, θ₀). Specifically, the observed frequency difference can be related to the local change in rotation rate through the integral of the kernel over the meridional plane. The total integral of the kernel over the meridional plane, ∑_ijK(r_i, θ_j)dΣ = 0.95 + 0.03 i, has a real part that is close to unity¹. We did not attempt a quantitative interpretation of the imaginary part since there is presently no established method to relate the imaginary part of the mode frequency to an observable (such as the mode amplitude, see Mushtaq et al. 2026, for a preliminary weakly nonlinear approach in 2D). Using the constraint (ω_obs − ℜω_ref)/2π = 7.7 ± 1.9 nHz, we estimated the perturbation to the rotation rate near (r₀, θ₀) as

$$\begin{array}{c}\hfill \delta \mathrm{\Omega }(r_0,{\theta }_0)/2\pi \approx (7.7\pm 1.9)/0.95\mathrm{nHz}=8.1\pm 2.0\mathrm{nHz}.\end{array}$$(5)

This small but significant enhancement of the rotation rate relative to the p-mode reference value (Ω_ref(r₀, θ₀)/2π = 357.2 nHz) shows that the HL1 mode provides a sensitive probe of differential rotation at high latitudes in the deep convection zone.

The above constraint is compared in Fig. 4 with a p-mode rotation estimate at (r₀, θ₀) obtained using optimally localized averaging (OLA; Howe 2009). This method combines p-mode frequency splittings to produce a resolution kernel localized at the target position (Pijpers & Thompson 1992, 1994). A significant discrepancy can be observed at a colatitude where OLA inversions are most uncertain. Remarkably, a single inertial mode yields a kernel with a comparable latitudinal extent.

4. Parametric model for rotation correction

To study the effect of a smooth modification to the rotation rate at high latitudes, we defined a simple parametric correction that depends only on the colatitude:

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Omega }(\theta ;A)=\{\begin{array}{cc}A\hfill & \theta <{10}^{°},\hfill \\ A[1+sin(4.5({30}^{°}-\theta ))]/2\hfill & {10}^{°}\le \theta <{50}^{°},\hfill \\ 0\hfill & {50}^{°}\le \theta <{90}^{°}.\hfill \end{array}\end{array}$$(6)

This rotation perturbation is symmetric across the equator. In the above equation, the profile reaches a maximum amplitude (A) at the pole and has a width of 30° (see Fig. 3a). This perturbation is radially constant, and it therefore also affects the rotation rate in the radiative zone. However, we showed in the previous section that the mode is insensitive to rotation in the radiative zone, and we chose not to include radial dependence of the profile for simplicity.

Fig. 3. Smooth high-latitude modification to the background rotation rate, ΔΩ(θ; A), and its effect on the HL1 mode frequency. (a) Functional form of ΔΩ(θ; A) as defined by Eq. (6), symmetric across the equator. (b) Variation of the real part of the mode frequency with A when using Solver 1 (green curve). By construction ω = ωref at A = 0. Agreement with the observed eigenfrequency (gray shaded region) is obtained for 10.3 nHz < A/2π < 16.0 nHz (range between the dotted green arrows). The value of A/2π = 13.6 nHz corresponding to the best fit is marked by the dashed green arrow. The dotted black line shows the linear approximation Eq. (7).

For infinitesimal perturbations in the angular velocity parameterized as above (i.e., as A → 0), the response of the mode frequency is

$$\begin{array}{c}\hfill \delta \omega (A)\approx {\displaystyle {\int }_{\mathrm{\Sigma }}}K(r,\theta )\mathrm{\Delta }\mathrm{\Omega }(\theta ;A)\mathrm{d}\mathrm{\Sigma }=(0.65+0.20\mathrm{i})A.\end{array}$$(7)

In the linear regime, A ≈ (ω_obs − ℜω_ref)/0.65. However, the range of validity of this linear approximation is not yet established, motivating the direct numerical computations of ℜδω(A) for finite-amplitude values of A.

We recomputed the eigenfrequency of the mode with the modified rotation rate Ω(r, θ) = Ω_ref(r, θ)+ΔΩ(θ; A) for values of A/2π ranging from 0 to 20 nHz. In all cases, the latitudinal entropy gradient was recomputed such that the thermal wind balance is maintained. The variation of the real part of the eigenfrequency with A is presented in Fig. 3b. The relationship between the real part of the eigenfrequency and A is roughly linear for A/2π ≲ 16 nHz. For A/2π between 10.3 nHz and 16.0 nHz, the real part of the eigenfrequency is within ±1σ of the observed frequency and the imaginary part of the frequency is positive (the mode is linearly unstable). A profile of the modified rotation rate at r = 0.8 R_⊙ is shown in Fig. 4 for $A/2\pi =13.6_{-3.3}^{+2.4}$ nHz, giving a polar rotation rate of 364.0 nHz.

Fig. 4. Constraints on high-latitude solar rotation rate at r = 0.8 R⊙. The green cross shows the constraint from the HL1 mode frequency (2010–2024), according to Eq. (5), with the horizontal bar indicating the width of the kernel. The green curve shows the rotation rate Ωref(r, θ)+ΔΩ(θ; A) for A/2π = 13.6 nHz, and the green shaded region indicates the associated uncertainty (see Fig. 3). The p-mode reference rotation rate Ωref(0.8 R⊙, θ) is given by the red dashed curve. The red cross shows a p-mode OLA estimate at the target location (0.8 R⊙, θ = 15° ) obtained from yearly HMI inversions (Howe 2009) averaged over 2010–2023. The horizontal bar represents the width of the averaging kernel, while the vertical bar represents the random error.

5. Conclusion

We find that the rotation kernel of the HL1 mode frequency is highly localized in the convection zone near r = 0.8 R_⊙ and θ = 15°. Without performing any inversion, the integrated HL1 kernel implies a local rotation rate of approximately 365.3 ± 2.0 nHz at this location, that is, 8.1 ± 2.0 nHz above the reference estimate from p-mode helioseismology (Larson & Schou 2018) and above a p-mode OLA inversion (Howe 2009). Using the smooth parametric profile prescribed in Sect. 4, the inferred enhancement reaches 13.6 ± 3 nHz at the pole relative to the reference model. These estimates are absolute measurements in the sense that we used the frequency of a single retrograde inertial mode. In contrast, acoustic modes infer the rotation rate using a differential measurement between prograde and retrograde modes, which reduces dependencies on the background model.

The frequencies of many other solar inertial modes have been measured (Gizon et al. 2021), as well as their variations with the solar cycle (Waidele & Zhao 2023; Lekshmi et al. 2026). In future work, we will aim to constrain the solar differential rotation at multiple positions in the convection zone using all modes simultaneously by solving an inverse problem. In this respect, the current situation is reminiscent of the early days of helioseismology, when Christensen-Dalsgaard & Gough (1976) foresaw the potential of p modes as probes of the solar interior. The validated eigenvalue solvers described in this paper will be valuable tools in this respect, particularly for interpreting the full spectrum of inertial modes and for assessing the contribution of time-varying zonal flows to their temporal variations (see Goddard et al. 2020).

Acknowledgments

PD and LG designed and performed research, PD and YB updated an earlier version of Solver 1, PD wrote the first draft of the paper and all authors contributed to the final manuscript. The authors acknowledge the essential contribution of Suprabha Mukhopadhyay, which enabled the cross-validation of Solver 1 with Solver 2. We thank Zhi-Chao Liang for making available the HL1 mode frequencies shown in Fig. 1 and Rachel Howe for providing the p-mode OLA rotation inversion shown in Fig. 4. We acknowledge useful discussions with Damien Fournier and Robert Cameron. We acknowledge partial financial support from ERC Synergy Grant WHOLESUN (grant agreement No. 810218) to LG. PD is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen.

References

Appendix A: Computation of eigenmodes of a rotating solar model

A.1. Linearized equations of motion

We work in the Carrington frame of reference, rotating at frequency Ω_Carr/2π = 456.0 nHz with respect to an inertial frame. We use the linearized forms of the momentum, entropy, and continuity equations, where perturbations in velocity (v′), pressure (p′), density (ρ′), and entropy (s′) are proportional to exp[i(mϕ − ωt)], where ϕ is the longitude and ω is the angular frequency of the mode of oscillation, to obtain

$$\begin{array}{cc}\hfill \mathrm{i}\omega {\mathit{v}}^{\prime }=& \frac{\mathrm{\nabla }{p}^{\prime }}{\rho }+\frac{{\rho }^{\prime }}{\rho }g\widehat{\mathit{r}}-2{\mathit{v}}^{\prime }\times (\mathrm{\Omega }\widehat{\mathit{z}})+rsin\theta ({\mathit{v}}^{\prime }·\mathrm{\nabla }\mathrm{\Omega })\widehat{\mathit{\varphi }}\hfill \\ \hfill & -\frac{\mathrm{\nabla }·\mathcal{D}}{\rho }+\mathrm{i}m(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{Carr}}){\mathit{v}}^{\prime },\hfill \end{array}$$(A.1)

$$\begin{array}{c}\hfill \mathrm{i}\omega s^{\prime }=\&{\mathit{v}}^{\prime }·\mathrm{\nabla }s-\frac{\mathrm{\nabla }·(\kappa \rho T\mathrm{\nabla }{s}^{\prime })}{\rho T}+\mathrm{i}m(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{Carr}})s^{\prime },\end{array}$$(A.2)

$$\begin{array}{c}\hfill \mathrm{i}\omega {\rho }^{\prime }=\&\mathrm{\nabla }·(\rho {\mathit{v}}^{\prime })+\mathrm{i}m(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{Carr}}){\rho }^{\prime }.\end{array}$$(A.3)

In the above equations, the unprimed variables (ρ, p, g, s, T) refer to the background reference state, and κ is the thermal diffusivity. The unit vectors $\widehat{\mathit{r}}$, $\widehat{\mathit{\varphi }}$, and $\widehat{\mathit{z}}$ point in the radial, longitudinal, and rotation-axis directions, respectively. An ideal gas equation of state is used to close the equations. The perturbations are complex functions of radius and colatitude. The components of the viscous stress tensor 𝒟 are given by

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mathit{ij}}=\rho {\nu }_{\mathrm{turb}}(S_{\mathit{ij}}-\frac{2}{3}({\delta }_{\mathit{ij}}\mathrm{\nabla }·{\mathit{v}}^{\prime })),\end{array}$$(A.4)

where the S_ij are the components of the deformation tensor in spherical coordinates (Fan & Fang 2014) and ν_turb is the turbulent viscosity.

A.2. Boundary conditions

Complemented by boundary conditions, the above equations constitute an eigenvalue problem. We choose the computational domain to span 0.65 R_⊙ to 0.985 R_⊙ in radius. We use impenetrable horizontal stress-free boundary conditions for the velocity at the top and bottom boundaries, and assume there is no entropy flux at either boundary:

$$\begin{array}{c}\hfill v_r^{\prime }=0,\frac{\partial }{\partial r}\left(\frac{v_{\theta }^{\prime }}{r}\right)=0,\frac{\partial }{\partial r}\left(\frac{v_{\varphi }^{\prime }}{r}\right)=0,\frac{\partial s^{\prime }}{\partial r}=0.\end{array}$$(A.5)

Further, we assume the perturbations v_r′, v_ϕ′, ρ′, and s′ to vanish at the poles (θ = 0 and π). For the specific case of m = 1 (our case), we require ∂v_θ′/∂θ = 0. This allows for a polar crossing flow. With these boundary conditions, Eqs. (A.1)–(A.3) form an eigenvalue problem, which is non-singular because of the viscous stresses Eq. (A.4). The spectrum of eigenvalues ω is discrete.

A.3. Reference solar model and reference differential rotation

We use the background density, pressure, and gravitational acceleration from the standard solar model S (Christensen-Dalsgaard et al. 1996). The reference internal rotation is provided by p-mode helioseismology (Larson & Schou 2018) applied to data from HMI/SDO. The time series of the angular velocity Ω(r, θ; t) is derived from successive data segments of length 72 days and 36 a-coefficients, a standard product available for download at the Joint Science Operations Center (JSOC) website under the filename hmi.v_sht_2drls. We adopt as the background rotation rate the mean profile from 30 April 2010 to 16 February 2024 and denote it by Ω_ref(r, θ), see Fig. A.1.

Fig. A.1. The reference solar angular velocity in the Carrington frame, Ωref − ΩCarr, from p-mode helioseismology (Larson & Schou 2018). This meridional cut shows the computational domain for the eigenvalue solvers, 0.65 < r/R⊙ < 0.985.

We need to specify ∇s in Eq. (A.2). Invoking the thermal wind balance for the latitudinal derivative of the entropy (Rempel 2005), we obtain

$$\begin{array}{c}\hfill {\mathit{v}}^{\prime }·\mathrm{\nabla }s=\frac{c_p\delta }{H_p}{v}_r^{\prime }+\frac{c_{p}rsin\theta }{\mathit{rg}}\frac{\partial ({\mathrm{\Omega }}_{\mathrm{ref}}^2)}{\partial z}{v}_{\theta }^{\prime },\end{array}$$(A.6)

where c_p is the specific heat at constant pressure, H_p is the pressure scale height, and δ = ∇ − ∇_ad is the superadiabatic temperature gradient. Because the convection zone is nearly adiabatically stratified, we set δ = 0 there and smoothly transition in the radiative interior (r/R_⊙ < 0.70) to the value of δ(r) from Model S.

We also need to specify ν_turb in Eq. (A.4). We adopt a turbulent viscosity profile in the form of a step function:

$$\begin{array}{c}\hfill {\nu }_{\mathrm{turb}}(r)={\nu }_{\mathrm{RZ}}+\frac{{\nu }_{\mathrm{CZ}}-{\nu }_{\mathrm{RZ}}}{2}[1+tanh(\frac{r/R_{\odot }-0.7}{0.015})],\end{array}$$(A.7)

where ν_RZ = 10⁹ cm² s⁻¹ and ν_CZ = 10¹² cm² s⁻¹ are representative values of the turbulent viscosities in the radiative and convection zones. We use the same profile for the thermal diffusivity.

A.4. Two eigenvalue solvers

Numerical eigenmodes may depend on the computational methods and codes. Therefore, we considered two independent codes, which use different discretization methods to solve the same problem. The first eigenvalue solver, referred to as ‘Solver 1’, is based on the finite-difference solver of Bekki et al. (2022). The second solver, ‘Solver 2’, is the Dedalus-based spectral solver described by Mukhopadhyay et al. (2025). Resolution-based convergence tests were performed to ensure numerical accuracy.

A.5. Solver 1 based on Bekki et al. (2022)

Solver 1 uses a uniform (staggered) grid in radius and colatitude, and a second-order finite difference scheme to compute derivatives. We use 100 radial and 90 latitudinal gridpoints, corresponding to resolutions of 2.3 Mm and 2°. The details of the grid and the computation methods was presented originally in Bekki et al. (2022). During validation, we found a number of bugs in the code, related to the treatment of the viscous and thermal diffusivities. We use a revised version of the code, where these bugs are corrected for. In addition, we upgraded the Fortran eigenvalue solver from the LAPACK package to ARPACK (a sparse solver, see Lehoucq et al. 1998) to speed up the computations.

At fixed m = 1, we compute the 100 eigenfrequencies nearest to the (real) target frequency of −90 nHz in the Carrington frame. We identify the mode with the largest growth rate (largest imaginary part of the eigenvalue) as the mode of interest. This mode has north-south symmetric radial vorticity, as observed on the Sun. Solver 1 gives a mode frequency of ℜω_ref/2π = −95.6 nHz in the Carrington frame (see solid green line in Fig. 1). This is also the only self-excited mode with symmetric radial vorticity in the computed eigenspectrum.

A.6. Solver 2 from Mukhopadhyay et al. (2025)

Dedalus (Burns et al. 2020) is an open-source, parallel code that solves partial differential equations in spectral space. It uses Chebyshev polynomials and spherical harmonics as radial and horizontal basis functions respectively to compute the eigenvalues of a system. Previously, Mukhopadhyay et al. (2025) used Dedalus to develop an eigenvalue solver for solar inertial modes (hereafter Solver 2). To validate Solver 1, we used Solver 2 to compute the eigenmodes within the same computational domain. Solver 2 used 84 and 32 gridpoints in radius and colatitude respectively. We specify the same target frequency and compute the same number of eigenmodes as for Solver 1. The eigenmode with the largest growth rate has a real frequency of −96.1 nHz (see dashed line in Fig. 1) and a similar eigenfunction to that from Solver 1.

A.7. Importance of the position of the lower boundary

As part of our numerical investigations, we found that the inclusion of the upper radiative zone r = 0.65 R_⊙ to 0.70 R_⊙ in the computational domain plays an important role in unambiguously identifying the HL1 mode. When considering only the convection zone, as in Mukhopadhyay et al. (2025) or Bekki et al. (2022), the system produces two pairs of self-excited eigenmodes within ∼10 nHz; each pair consisting of two modes with opposite latitudinal symmetries in their eigenfunction. This is at odds with the observation of a single mode dominating the power spectrum. When the upper layers of the radiative zone are included in the model, however, we find only one pair of self-excited eigenmodes, with the fastest growing mode (eigenfrequency with the largest imaginary part) having the same north-south symmetry as that dominant in the observations. The surface profile of this mode is very similar to that extracted from surface observations.

Appendix B: Additional figures

Fig. B.1. Cuts at different radii though the kernel presented in Fig. 2.

Fig. B.2. Cuts at different latitudes of the kernel presented in Fig. 2.

All Figures

Fig. 1. Observed HL1 frequencies in the Carrington reference frame extracted from HMI Dopplergrams at a cadence of one year with a three-year sliding window (blue points with error bars, Liang & Gizon 2025). The mean frequency ωobs/2π = −87.9 nHz is indicated by the black horizontal line, and the gray shaded region shows the ±1σ range (±1.9 nHz). The frequency uncertainty was estimated from measurements over independent time intervals centered on 2012, 2015, 2018, and 2021, thereby accounting for correlations introduced by the sliding window. The solid green line marks the real part of the eigenfrequency computed using eigenvalue Solver 1 (−95.6 nHz), and the dashed green line marks the same for Solver 2 (−96.1 nHz). The azimuthal velocity of the observed mode at the surface is shown as an inset (adapted from Gizon et al. 2021).

In the text

Fig. 2. Spatial sensitivity of the frequency of the HL1 mode to localized changes in the background rotation rate. The kernel K(r, θ) is symmetric across the equator. Grid sizes dr = 2.3 Mm and dθ = 2° were used for this calculation. For reference, a patch of area 10dr × 2rdθ = 20dΣ is outlined with a solid black contour. The dashed black contour corresponds to the value of the kernel at half maximum.

In the text

Fig. 3. Smooth high-latitude modification to the background rotation rate, ΔΩ(θ; A), and its effect on the HL1 mode frequency. (a) Functional form of ΔΩ(θ; A) as defined by Eq. (6), symmetric across the equator. (b) Variation of the real part of the mode frequency with A when using Solver 1 (green curve). By construction ω = ωref at A = 0. Agreement with the observed eigenfrequency (gray shaded region) is obtained for 10.3 nHz < A/2π < 16.0 nHz (range between the dotted green arrows). The value of A/2π = 13.6 nHz corresponding to the best fit is marked by the dashed green arrow. The dotted black line shows the linear approximation Eq. (7).

In the text

Fig. 4. Constraints on high-latitude solar rotation rate at r = 0.8 R⊙. The green cross shows the constraint from the HL1 mode frequency (2010–2024), according to Eq. (5), with the horizontal bar indicating the width of the kernel. The green curve shows the rotation rate Ωref(r, θ)+ΔΩ(θ; A) for A/2π = 13.6 nHz, and the green shaded region indicates the associated uncertainty (see Fig. 3). The p-mode reference rotation rate Ωref(0.8 R⊙, θ) is given by the red dashed curve. The red cross shows a p-mode OLA estimate at the target location (0.8 R⊙, θ = 15° ) obtained from yearly HMI inversions (Howe 2009) averaged over 2010–2023. The horizontal bar represents the width of the averaging kernel, while the vertical bar represents the random error.

In the text

	Fig. A.1. The reference solar angular velocity in the Carrington frame, Ωref − ΩCarr, from p-mode helioseismology (Larson & Schou 2018). This meridional cut shows the computational domain for the eigenvalue solvers, 0.65 < r/R⊙ < 0.985.
In the text

	Fig. B.1. Cuts at different radii though the kernel presented in Fig. 2.
In the text

	Fig. B.2. Cuts at different latitudes of the kernel presented in Fig. 2.
In the text