A Linear Model for Inertial Modes in a Differentially Rotating Sun

Inertial wave modes in the Sun are of interest owing to their potential to reveal new insight into the solar interior. These predominantly retrograde-propagating modes in the solar subsurface appear to deviate from the thin-shell Rossby–Haurwitz model at high azimuthal orders. We present new measurements of sectoral inertial modes at m > 15 where the modes appear to become progressively less retrograde compared to the canonical Rossby–Haurwitz dispersion relation in a corotating frame. We use a spectral eigenvalue solver to compute the spectrum of solar inertial modes in the presence of differential rotation. Focussing specifically on equatorial Rossby modes, we find that the numerically obtained mode frequencies lie along distinct ridges, one of which lies strikingly close to the observed mode frequencies in the Sun. We also find that the n = 0 ridge is deflected strongly in the retrograde direction. This suggests that the solar measurements may not correspond to the fundamental n = 0 Rossby–Haurwitz solutions as was initially suspected, but to those for a higher n. The numerically obtained eigenfunctions also appear to sit deep within the convection zone—unlike those for the n = 0 modes—which differs substantially from solar measurements and complicates inference.


INTRODUCTION
Global-scale oscillations such as Rossby waves (Rossby 1945) are characteristic of rotating atmospheres, and they have been studied in terrestrial settings (Pedlosky 1987(Pedlosky , 2003)), as well as in astrophysical settings (Lou 2000;Lanza et al. 2009;Zaqarashvili et al. 2021).Waves with similar temporal and spatial characteristics have also recently been detected on the Sun (Löptien et al. 2018;Liang et al. 2019;Hanasoge & Mandal 2019;Mandal & Hanasoge 2020;Hanson et al. 2020;Proxauf et al. 2020;Mandal et al. 2021;Waidele & Zhao 2023).A key characteristic of these largely toroidal oscillation modes is that the primary restoring contribution to sustain them is provided by the Coriolis force.Several authors have expounded on the linearized theory of inertial modes in slowly rotating stars (Provost et al. 1981;Saio 1982).Typically, the angular profiles of modes on the Sun are interpreted in a basis of spherical harmonics, Y ℓ,m (θ, ϕ), where ℓ is the angular degree and m is the azimuthal order, and θ and ϕ represent the co-latitude and azimuth respectively.The dominant modes that are observed on the Sun are sectoral, that is, most of the power in the mode is concentrated in the channel ℓ = m, and these modes are symmetric about the equator.More recently, observations of equatorially antisymmetric modes of oscillation by Hanson et al. (2022) have been interpreted as the fundamental full-sphere oscillations in the solar convection zone corresponding to ℓ = m + 1 (Triana et al. 2022).Several early studies had chosen to analyze linear oscillations about a uniformly rotating background, but the differential rotation of the Sun may be expected to contribute significantly to the spectrum of modes.This is because, firstly, the Doppler shift encountered by a mode with an azimuthal order m -given by m∆Ω for a constant rate of differential rotation ∆Ω about a reference rotation rate of Ω 0 -may dominate the canonical dispersion relation of sectoral Rossby modes −2Ω 0 /(m + 1) at high m, and secondly, non-uniform rotation couples modes across the harmonic orders to a much larger degree than uniform rotation, which makes mode identification and labelling challenging.The former leads to critical latitudes, where the inviscid eigenvalue problem becomes formally singular (Watts et al. 2004;Gizon et al. 2020).In certain conditions, differential rotation might induce instabilities in the modes of oscillation (Zaqarashvili et al. 2021;Fournier et al. 2022).The use of these modes in helioseismology, therefore, requires a comprehensive understanding of the dispersion relation and stability of inertial modes in a solar-like background rotation.
Several authors (Gizon et al. 2020;Fournier et al. 2022;Bekki et al. 2022a,b) have carried out numerical analysis of inertial modes in the presence of a solar-like differential rotation profile.Each set of authors reported r-mode dispersion relations where the real parts of the frequencies were close to the central frequencies of observed solar inertial modes.Gizon et al. (2020) suggest that equatorial modes are confined between critical latitudes, and that other categories of modes -such as critical-latitude modes and high-latitude modes -might exist in a differentially rotating medium.These modes have putatively been detected on the Sun as well (Gizon et al. 2021).Hathaway et al. (2013), and later Bogart et al. (2015) and Hathaway & Upton (2021), found evidence for polar vortices, which might correspond to a self-excited, high-latitude m = 1 mode driven by baroclinicity (Bekki et al. 2022b).Bekki et al. (2022a) have detected several classes of inertial modes in nonlinear simulations of rotating convection, including equatorial modes, the high-latitude m = 1 feature, as well as columnar convective modes.Of these, equatorial inertial modes are the easiest to identify, as their frequencies are close to the canonical thin-shell Rossby-Haurwitz dispersion relation for sectoral modes (ω = −2Ω 0 /(m + 1) for a uniform rotation rate Ω 0 ).
In this work, we compute the spectrum of the equatorial inertial modes numerically in the anelastic approximation by including a solar-like differential rotation profile.This approach differs from previous studies, where differential rotation had been studied either as a radius-dependent latitudinal profile (Fournier et al. 2022), or as a radially and latitudinally varying profile that is confined to a truncated radial domain (Bekki et al. 2022b).The latter of these studies comes the closest to our analysis, although the specifics of the analysis differ in the way compressibility is addressed, the profile of rotation that is used in the analysis, the tracking rate used to compute the mode frequencies, as well as in the numerical approach (finite difference in a meridional plane in (Bekki et al. 2022b), whereas spectral in ours).We use a Petrov-Galerkin discretization of the equations of mass, momentum and energy, following the approach presented by Bhattacharya & Hanasoge (2023).A spectral approach often leads to exponential convergence for smooth functions (Gottlieb et al. 1978), which may reduce the resolutions required to obtain convergent solutions, and consequently make the eigenvalue problem in a non-uniform background more tractable.The nominally higher accuracy expected from spectral methods also affords more trust in the solutions.One of the main questions that we seek to answer in this work is why the dispersion relation of the Rossby waves that are observed on the Sun is consistent with that seen in a thin-shell such as the earth's atmosphere, when the radial extent of the solar convection zone spans many scale heights, and is therefore far from a thin shell.Furthermore, the measured dispersion relation appears to progressively deviate from the Rossby-Haurwitz model (becoming less retrograde than expected at m > 10), and we explore whether this may be explained in terms of the solar differential rotation.

WAVE EQUATION
In this analysis, we use a spherical coordinate system, with r representing the radius, and θ and ϕ standing for the co-latitude and azimuth respectively.We model waves as small perturbations in a spherical shell with a lower radial bound denoted by r in and an outer bound r out .We choose the background to correspond to Model S (Christensen-Dalsgaard et al. 1996) that is augmented by a super-adiabaticity profile inspired by Rempel (2005) (see (Bhattacharya & Hanasoge 2023) for details).The background profile may be represented in terms of its radially varying profiles of density ρ(r), pressure p(r), temperature T (r), gravitational acceleration g(r) that is directed inwards, as well as a radially and latitudinally varying profile of entropy S(r, θ).We denote perturbations to the background parameters by using primes as superscripts.We denote the wave velocity by u(r, θ, ϕ, t), and its vorticity ∇ × u(r, θ, ϕ, t) by ω(r, θ, ϕ, t).Wherever unambiguous, we drop the spatial and temporal coordinates in the interest of brevity.
We assume that the medium is close to the anelastic limit (Gough 1969;Gilman & Glatzmaier 1981;Braginsky & Roberts 1995), where the equation of mass continuity takes the form (1) In this scenario, we assume that fluid velocities are highly subsonic, and therefore filter sound waves out of the analysis.This decouples pressure and density from the fluid velocity.Pressure perturbations do not evolve in time and may be solved for at each instant given the velocity profile, while the density perturbation may be obtained from an equation of state.This reduces the number of unknowns that we solve for to the velocity profile and entropy perturbation S ′ .Imposition of the anelastic approximation further restricts the number of independent components of the velocity field to two.We represent the rotation velocity of the medium by Ω = Ω(r, θ) e z , and assume that this is being tracked in a frame rotating at a constant angular velocity Ω 0 = Ω 0 e z .In a uniformly rotating medium where Ω is a constant in space, one may choose the tracking rate Ω 0 to be equal to Ω, and analyze wave propagation in a co-rotating frame.In the Sun, however, differential rotation contributes additional terms to the wave equation.We denote the e z -component of the differential rotation velocity ∆Ω = Ω − Ω 0 by ∆Ω(r, θ).
We represent the momentum equation in terms of the velocity u, the vorticity ω = ∇ × u, and the entropy perturbation S ′ , in the form presented by Lantz (1992) and Braginsky & Roberts (1995).The Navier-Stokes equation in a rotating frame in the anelastic approximation may be expressed as where the viscous force term F ν is given by and we assume that ν is a constant.We have applied the vector identity u • ∇u = 1 2 ∇ |u| 2 − u × ω to obtain terms on the right-hand side.This form of the equation emphasizes that only the u × ω term of the two contributes to the time-evolution of vorticity, as the curl of a gradient is zero.Bhattacharya & Hanasoge (2023) had analyzed the terms that correspond to a uniformly rotating background, and in this work, we focus on the additional terms arising from differential rotation.The u × ω term is second-order in u in a uniformly rotating medium.In a differentially rotating medium, however, there arises a first-order contribution due to a coupling between the background rotation and the fluid velocity.We may expand the velocity field in terms of the intrinsic fluid velocity u f and the velocity u Ω that is associated with differential rotation as u = u f + u Ω , where u Ω = ∆Ω × r.Defining ω Ω = ∇ × u Ω , we may expand the advection term as where we have assumed that the centrifugal component is negligible in comparison with gravity, and we have ignored the second-order terms in u f as well.In addition to the advection term, the Coriolis force due to differential rotation may be expressed as −2∆Ω × u f .Collectively, the net contribution of differential rotation becomes The second term in Equation (5) contributes −im∆Ω ω f r to the time-evolution of the radial component of vorticity, which we may recognize as a Doppler shift.Since this term is linear in m whereas the sectoral Rossby-Haurwitz dispersion relation 2Ω/(m + 1) varies inversely with m, we may expect the Doppler-shift to contribute significantly for m ≫ 1.To analyze the impact of the first term in Equation ( 5), we note that if ∆Ω is a constant in space (which we may denote by ∆Ω const ), the corresponding vorticity may be expressed as ω Ω = 2∆Ω const .The first term in Equation (5) consequently becomes −4∆Ω const × u f , which acts as an augmented Coriolis force.The impact of this term on the spectrum would be a scaling of the rotation velocity in the Rossby-Haurwitz dispersion relation.The rotation profile in the Sun is far from constant, but we may interpret the result in terms of an appropriately averaged rotation rate over the domain.
We do not consider the contribution of differential rotation to the viscous force.Fournier et al. (2022) have demonstrated that the real and the imaginary parts of the eigenfrequencies are related, and increasing the Ekman number induces a change in both (especially for Ekman numbers above 10 −3 ), although, out of all the solutions, the r mode is the least affected by an increase in the Ekman number.Bekki et al. (2022b) also find a mild dependence of the equatorial Rossby modes on the viscosity coefficient.In this analysis, we are in a weak-viscosity regime (ν = 5 × 10 11 cm 2 /s, corresponding to an Ekman number of around 10 −5 ), and therefore this impact may be expected to be insignificant.In any case, our assumption of the model of how turbulent viscosity impacts inertial waves is simplistic, and more realistic models that incorporate the depth-dependence of viscosity might be necessary to conclusively comment on the impact of differential rotation on the line widths on modes (see e.g., Muñoz-Jaramillo et al. (2011) for a discussion on various models of depth-dependent diffusivity).However, it remains to be determined whether these sorts of subtleties induce observable changes to the modes.
The energy equation may be expressed in terms of the entropy perturbation S ′ and the wave velocity u f as We have disregarded other terms such as radiative and viscous heating in this analysis, but such terms may be included in a straightforward manner.The radiative heating term contributes to thermal transport, so an improved treatment of the energy equation would include such a term (Featherstone & Hindman 2016).We may simplify Equation ( 6) further by noting that the radial derivative of the entropy profile may be expressed in terms of the super-adiabatic gradient δ as ∂ r S = −γδ/H ρ , where γ is the adiabatic exponent, and H ρ is the density scale-height.An estimate of the latitudinal gradient of entropy is found to be ∂ θ S = 2(c p /g) r 2 sin θ Ω 0 • ∇∆Ω by requiring thermal-wind balance.Bekki et al. (2022b) find that this term contributes significantly to the instability of the m = 1 high-latitude mode, and, as a consequence, the degree of latitudinal non-adiabaticity may be probed through solar inertial modes.The influence of this term on equatorial modes may be limited (Bhattacharya & Hanasoge 2023), since the term is proportional to ∂ θ ∆Ω, which is small near the equator compared to the reference rotation rate Ω 0 , and increases with latitude (see Fig 1).Owing to computational limitations, we do not include the latitudinal entropy gradient term in our analysis.
We choose a poloidal-toroidal decomposition of the wave velocity u f that automatically satisfies the continuity equation.The azimuthal symmetry of the rotational angular velocity Ω implies that ρu Ω is divergence-free, so the continuity equation translates to that for the wave velocity ∇ • (ρu f ) = 0. We express the wave velocity u f in terms of stream functions as where V (x) corresponds to the radial component of the wave vorticity, and W (x) corresponds to the radial component of the wave velocity.We may express the velocity components in terms of the stream functions as where ∇ h represents the angular component of the Laplacian operator.This suggests that in the special case where the field is purely toroidal, sectoral, and featuring a single azimuthal order m, the stream function V would have an angular profile proportional to sin m θ, and v θ would have a latitudinal dependence given by sin m−1 θ.In this case, v θ would be symmetric about the equator, whereas v ϕ would be antisymmetric.Substituting Equation ( 7) in ( 2) and ( 6), we obtain coupled equations for V (x), W (x) and S ′ (x).Transforming to temporal frequency, the equations reduce to an eigenvalue problem where the eigenvalues are mode frequencies ω, and the eigenvectors are [V (x), W (x), S ′ (x)].

Discretizing operators
Bhattacharya & Hanasoge (2023) had presented an approach to discretize the eigenvalue equation for inertial waves in a spectral basis spanned by Chebyshev polynomials in the radial direction, and spherical harmonics in the lateral coordinates.They applied the approach to a background medium that was either rotating uniformly or differentially in radius, but latitudinally uniform.In this work, we extend their approach to a latitudinally differentially rotating background, bearing in mind that the solar rotation profile is azimuthally symmetric.We choose to work in spherical polar coordinates, with the North Pole directed along the Sun's axis of rotation.In this convention, θ represents the co-latitude, and ϕ the azimuth.We limit our radial domain to r in = 0.6R ⊙ to r out = 0.985R ⊙ owing to concerns of numerical expense.We define the non-dimensional radius r as which maps the radial domain The Doppler shift induced by differential rotation on the inertial-wave spectrum depends critically on the rotation profile, and we take extra care to incorporate the profile in its entirety in our domain.We choose a second scaled radial coordinate In terms of this radius, we define our modified rotation profile Ω (r, θ) as which squeezes the outer envelope of the rotation profile in the Sun (Larson & Schou 2018) to the radial domain used in this analysis.This procedure ensures that the impact of the near-surface shear layer is incorporated in the analysis while avoiding the steep density gradient associated with r > 0.985R ⊙ .From this point onward, we treat the modified profile Ω (r, θ) as the one that features in Equation ( 2).We may expand this profile in a Chebyshev-Legendre tensor-product basis as where T n (r) is the Chebyshev polynomial of degree n, and P ℓ (cos θ) is the normalized Legendre polynomial of degree ℓ.The ℓ = 0 term in the sum corresponds to the latitudinally uniform component of rotation, that had been studied by Bhattacharya & Hanasoge (2023), of which the n = 0 term corresponds to uniform rotation.We compute the expansion coefficients Ω nℓ by firstly fitting a smoothing spline to the rotation angular velocity, and secondly performing an adaptive orthogonal polynomial fit to the spline.The smoothing procedure reduces the number of coefficients that contribute in Equation ( 13).We set our reference tracking rate to Ω 0 = 2π × 453.1 nHz to match the commonly used mean equatorial rotation rate at the solar surface that has been determined using f-mode seismology (Schou 1999).This matches the rate used by Löptien et al. (2018) and subsequent papers.We plot the rotation profile, as well as the equatorial profile of the differential component of the profile, in Figure 1.We expand the stream functions V and W and the entropy S ′ in a Chebyshev-spherical-harmonic basis as where T n (r) are Chebyshev polynomials and Pℓm (cos θ) are normalized associated Legendre polynomials.Substituting Equation ( 14) into ( 2) and ( 6), we recast the eigenvalue problem in terms of the expansion coefficients V nℓm , W nℓm and S ′ nℓm .Azimuthal symmetry of the background leads to the equations decoupling in the azimuthal order m, and we may solve a set of smaller eigenvalue problems in each m, where the eigenvectors are the collated coefficients [V nℓm , W nℓm , S ′ nℓm ] for that m.We describe the approach to computing the matrix elements of the angular terms in Appendix A. The discretized representation of a separable operator is given by the tensor product of the corresponding radial and angular matrices, whereas for a non-separable operator such as the rotation profile, we may express the result as a sum over separable terms, as in Equation ( 13).We use the Julia package ApproxFun.jl(version 0.13, Olver & Townsend 2014) to obtain this discretized representation, and subsequently use the LAPACK wrappers presented by the Julia standard library LinearAlgebra to perform the eigen-decomposition (using Julia version 1.9).The code for this is published freely on GitHub under the MIT license1 .

Boundary conditions
We impose impenetrable, stress-free boundary conditions at the radial extremities.We also require that there be no entropy flux across the boundaries.This may be expressed in terms of the velocity components as We enforce the boundary conditions by following a basis-recombination approach as presented by Bhattacharya & Hanasoge (2023).
We do not impose additional boundary conditions at the poles aside from the ones that naturally emerge from the choice of spherical harmonics as the basis, which ensures smoothness of the stream functions at the poles.This may be seen as a behavioral boundary condition rather than a numerical one (Boyd 2000).

DATA ANALYSIS
The mode properties (frequency, line width etc.) of the sectoral Rossby waves have been studied extensively up to m = 15 (Löptien et al. 2018, and subsequent papers).Here, we aim to characterize higher-order modes, i.e., m > 15.
There are a number of factors that have made it difficult to analyze these higher-order modes.Firstly, the amplitude of the modes decreases with increasing m, leading to poor signal-to-noise ratios for the fitting of the modes.This makes it difficult to characterize these high degree modes with techniques that have poor signal-to-noise across the spectrum (e.g.time-distance measurements).Furthermore, some data products have a low spatial Nyquist frequency.For instance, the GONG ring diagram pipeline (Corbard et al. 2003) has tiles separated by 15 • , leading to a spatial Nyquist of m = 12.Beyond m = 15 aliasing pollutes the spectrum and makes mode characterization impossible for the Global Oscillation Network Group (GONG) data.However, the flow maps from the Helioseismic and Magnetic Imager (HMI) ring pipeline have a Nyquist frequency of m = 24, and here we average all depths down to ∼ 20 Mm below the surface to improve the signal-to-noise ratio (e.g.Hanson et al. 2022).Global mode coupling and local correlation tracking (LCT) are not as limited in the spatial frequencies, though as reported by Löptien et al. (2018) the latter has weak power at high m.Here, we utilize toroidal flows computed from the global mode coupling analysis (MCA) for both odd and even m (Mandal et al. 2021), as well as the HMI ring diagram analysis (RDA) flows averaged for all depths down to ∼ 20 Mm.Note-Fits for each m are performed in a frequency window 400 nHz wide, centered on the dispersion relation −2Ω0/(m+ 1).Negative frequencies indicate retrograde motion.We were unsuccessful in performing fits to the ℓ = m = 20 mode.
Figure 2 shows the radial vorticity (toroidal flow) power spectrum of the HMI MCA and RDA.In both cases, the mode ridge begins to shift away from the thin-shell dispersion relation, becoming more prograde with increasing m.At these relatively high azimuthal order, any small Doppler shift will become significant due to the factor m in the frequency shift.In both the MCA and RDA, there is an apparent mode power beyond m = 15 up to m = 21, where the modes appear to become slightly prograde (relative to the equatorial rotation frame of reference).Beyond the azimuthal order m = 21, there is no discernible signal in any of the measurements.
Using a Markov Chain Monte-Carlo method for fitting (eg.Hanson et al. 2022), we fit the mode ridges for both HMI MCA and RDA and compute the associated errors.Figure 3 shows slices of the power spectra, and the fitted Lorenztians for 15 < m < 22.For ℓ = m = 20 the fitting routine failed to provide reliable results that were not dependent on the initial guess.This is due to low signal-to-noise ratio, and we do not provide the fit information.Table 1 shows the fit parameters and their respective 68% confidence intervals.Figure 4 visualizes the mode fit parameters from Table 1, as well as the modes m > 2 which have been previously characterized.The mode frequencies from the RDA and MCA appear to have a trending difference beyond m = 10.This results in the MCA modes appearing to travel prograde relative to the tracking rate.Systematic differences in the mode frequencies between different methods are well known, and we show the travel time measurements (Liang et al. 2019) to further demonstrate this (see also Fig. 3 of Hanson et al. 2020).Small differences in the tracking rate, or low frequency noise, could lead to strong biases at large m.Further investigation is required to identify the differences between different data sets and methods, that lead to these discrepancies at high m.While there are differences in the mode frequencies, the line width measurements agree within errors.Finally, in the bottom panel of Fig. 4 we show how the signal-to-noise ratio of the modes rapidly declines with increasing m, reaching values close to 1 beyond m = 19.(Mandal et al. 2021), on HMI data (2010-2020).Right panel: Radial vorticity spectrum computed through ring diagram analysis and averaging all depths, e.g.Hanson et al. (2022).The red line represents the sectoral Rossby dispersion in the thin shell approximation (−2Ω0/(m + 1)), while the blue dashed line shows the maximum azimuthal order to which previous studies have characterized the modes.

Tracking at the equatorial rate
We have computed the spectrum of inertial modes by choosing 60 Chebyshev polynomials along the radial direction, and 40 harmonic degrees in the latitudinal direction.Since we can separate out symmetric and antisymmetric modes, the effective angular resolution doubles, and we only retain harmonic degrees corresponding to even ℓ+m while solving for the symmetric modes.We plot the spectrum of latitudinally symmetric (about the equator) modes in Figure 5.We only select sufficiently smooth solutions that have 90% of their power concentrated below the radial Chebyshev order n = 10 and the harmonic degree ℓ = m + 15 for a specific azimuthal order m.We have also restricted ourselves to modes that have two radial nodes at most.Alongside the spectrum that we obtain numerically, we have indicated the observed mode frequencies as measured by Liang et al. (2019) (time-distance (TD) analysis applied to combined HMI and MDI data, blue error bars), Hanson et al. (2020) (ring-diagram (RD) analysis applied to GONG data, red error bars), Proxauf et al. (2020) (ring-diagram analysis applied to HMI data, green error bars) and the new mode frequencies obtained by applying mode-coupling analysis to HMI data ((HMI, MC), denoted by black error bars).On the right panel, we plot the spectrum with the Rossby-Haurwitz dispersion relation −2Ω 0 /(m + 1) subtracted from the frequencies.Distinct ridges in the spectrum stand out, along which the eigenfunctions bear qualitative resemblance.We have indicated four of these ridges to which we refer as "ridge-1", "ridge-2", "ridge-3" and "ridge-4".Ridge-2 lies strikingly close to the observed mode frequencies.
We plot the doubled imaginary parts of the Rossby-ridge eigenfrequencies in Figure 6, along with the measured mode line-widths on the Sun.In this analysis, we have set the kinematic viscosity coefficient ν to 50 km 2 /s.We find that there is a correspondence between the observed line widths and the numerically obtained values, even though the level of turbulent damping chosen is less than the 100 km 2 /s estimate by Gizon et al. (2021) for surface modes.
A subtlety that arises in interpreting the components of the eigenfunctions is that they are only determined to within an overall phase factor exp(iψ), where ψ is an arbitrary constant.To account for this, we ensure that the normalization process sets the (n = 0, ℓ = m) component for the toroidal component V to be purely real.The eigenfunctions corresponding to ridge-2 appear to change their nature with increasing m, starting off as predominantly localized at the surface for m ≤ 4 to peaking near the base of the convection zone for m ≥ 8, with a transition around m = 6.We plot the real part of the eigenfunction V for various values of m in Figure 8.We show eigenfunctions for ridge 1 in Fig. 7, those for ridge-3 in Fig. 9 and for ridge-4 in Figure 10.Ridge-3 eigenfunctions appear to change their nature with m in a manner that is similar to those for ridge-2, with the functions being peaked at the surface at m < 4, and progressively shifting to the bottom of the convection zone at higher m.However, this effect is much less pronounced than that for ridge-2, and at m > 4, there appear to be numerical artifacts that develop at the surface, possibly originating from the lower boundary condition.The functions corresponding to ridge-1 as well as ridge-4 are of interest, as these are localized at the surface.In particular, the functions corresponding to ridge-4 display latitudinal nodes in a manner similar to the fitting function used by Proxauf et al. (2020).However, the mismatch between the dispersion relation for ridge-4 and the data is considerable, so we refrain from establishing a correspondence between the two.At low m, eigenfunctions for ridge-3 appear to correspond to n = 0, whereas those for ridge-2 correspond to the radial order n = 1, although such an identification becomes difficult at higher m.
To trace the spectral ridges back to those for a uniformly rotating medium, we compute the inertial-mode spectrum by scaling the differential rotation ∆ Ω(r, θ) by a constant f for various values of f lying in (0, 1).We plot the spectrum for four different values of f in Figure 11, after adding 2Ω 0 /(m + 1) to the mode frequencies.We have highlighted three distinct ridges that progressively drift away from the Rossby-Haurwitz dispersion relation (indicated by the horizontal dotted line at an ordinate of zero).The largest deviation from the linear behavior is seen at low m, where the expectation that the 1/(m + 1) dependence of the dispersion relation dominates the Doppler shift (which  1.We show the fits for the HMI RDA (red), HMI MCA (black) and for reference the measurements from time-distance power spectra (Liang et al. 2019, blue).Top panel: Difference between the measured mode frequency and the thin-shell dispersion relation.Middle panel: Measured line widths.Bottom panel: signal-tonoise ratio of the modes measured from the spectra.The shaded area shows the 68% confidence level.
is proportional to m), and more so if the differential rotation rate is small compared to the reference rotation rate Ω 0 .The filled patches line in each sub-plot represents the lowest m for which there exist critical latitudes, where the phase speed of sectoral Rossby modes equals the local differential rotation speed.Modes to the right of the filled patch line are influenced strongly by differential rotation, whereas the ones to the left are impacted weakly.We find that this roughly matches the point where ridge-3 starts to depart from the zero ordinate.For the solar rotation profile, critical latitudes exist for all m > 2 (Bekki et al. 2022b), which explains the shape of ridge-3 in Figure 5.Despite this temptingly simple interpretation in terms of Doppler shifts, we note that not all the ridges appear to scale linearly with the factor f , and the eigenfunctions at high m are significantly altered from those in a uniformly rotating medium.
We note that there is no reason to expect distinct ridges of eigenvalues that vary linearly with m, given that the eigenfunctions span a large spatial domain and sense a wide range of rotation rates.An analysis of the dispersion relations starting from the sensitivity kernels for these inertial modes may shed light on why this happens to be the case.

Tracking at the Carrington rate
We repeat the analysis in a truncated radial domain spanning r = [0.71R⊙ , 0.985R ⊙ ] to compare our results with Bekki et al. (2022a).In this analysis, we truncate the solar rotation profile to the domain without scaling the radius, so the near-surface shear layer as well as the tachocline are both left out of the profile (although we find that including the near-surface shear layer does not impact the spectrum significantly).We further choose the tracking rate Ω 0 to be the Carrington rotation rate 2π × 456 nHz.We plot the spectrum in Figure 12.There is a correspondence between the ridges for n = 0 and n = 1 that are reported by Bekki et al. (2022b) and that in our result, although the match is not exact.Further work is necessary to improve the quantitative match between the two approaches by identifying the sources of systematic differences.

DISCUSSION AND CONCLUSION
In this work, we have computed the spectrum of equatorial inertial modes in the Sun in the anelastic approximation, assuming a solar-like profile of differential rotation.We find that the mode frequencies often lie along distinct ridges.By tuning the strength of differential rotation, we may trace them back to those obtained in a uniformly rotating system, for which the dispersion relations are easier to interpret.We find that the spectral ridges vary nearly linearly with m, which suggests that the shifts in dispersion relations may be interpreted as those due to a constant effective rotation rate.This result is surprising because the spatial extent of the mode eigenfunctions would indicate that each mode ought to sample the rotation rate differently.Further study into the sensitivities of the modes to the rotation profile might shed light on this.
An interesting result that we obtain is that the sectoral ridge with the radial order n = 0 -which coincides with the canonical Rossby-Haurwitz dispersion relation ω = −2Ω 0 /(m + 1) for a uniformly rotating medium -is strongly deflected in the retrograde direction.It is likely, therefore, that the observed mode frequencies correspond to one of the other ridges that we identify in the spectrum.A similar conclusion was reached by Bekki et al. (2022a), where the authors suggested that the observed frequencies might correspond to modes with n = 1.Based on the proximity of the mode frequencies to the measured frequencies, we speculate this to be the case as well.However, we note that the eigenfunctions for these modes are localized near the base of the convection zone for m > 8, unlike the near-surface measurements on the Sun.Moreover, the mode frequencies at m < 7 for ridge-2 appear to differ considerably from solar measurements, although this might be because these modes extend deep into the Sun and are therefore more sensitive to the lower boundary condition.Ridges 1 and 4, on the other hand, have eigenfunctions that are localized near the surface, and, in particular, the latter has eigenfunctions that seem to resemble the latitudinal form suggested by Proxauf et al. (2020).However, the dispersion relation for modes along this ridge lies considerably away from the observed frequencies, so these are unlikely to be candidates for the observed modes.
To get a complete picture of inertial modes under the effect of differential rotation, we have also presented the measurements of solar sectoral Rossby-Haurwitz waves with azimuthal orders greater than 15, using ring-diagram and mode-coupling analysis.These modes have eluded previous studies due to signal-to-noise or aliasing effects, and the first measurements to our knowledge were reported recently by Waidele & Zhao (2023) using time-distance seismology.The most interesting finding from these measurements is the increasing discrepancy between the thin-shell dispersion relation and the observations.While MCA and RDA have a small disagreement -which warrants further investigation into data sets and methodology -they both suggest a strong divergence of the solar modes from the theory towards less-retrograde frequencies.This trend with increasing m is qualitatively similar to the behavior of the ridge-2 from our models.
A few words are in order on the questions left untouched in this work.We have focussed strictly on the equatorial modes, ignoring the wider spectrum of critical-latitude and high-latitude modes that have purportedly been detected.We have also chosen to focus on modes with m > 1, which excludes the high-latitude spiral mode that was detected by Hathaway & Upton (2021).Moreover, we have limited our analysis strictly to equatorially symmetric modes, thereby omitting the high-frequency Rossby modes that were detected by Hanson et al. (2022).Further analysis using this approach may choose to focus on these aspects.Additionally, our analysis is purely hydrodynamic, but magnetic fields may play an important role governing the dispersion relation (Dikpati et al. 2020).In particular, total vorticity is no longer conserved in the presence of Lorentz forces (Zaqarashvili et al. 2021), although the impact of this may not be significant in the solar near-surface layers, where the dispersion relation of the observed modes appear to closely resemble that of classical Rossby waves.The approach presented here may be extended to include magnetic fields, which might unveil a richer spectrum of modes, such as atmospheric modes (McIntosh et al. 2017), as well as those that may exist in the solar tachocline (Zaqarashvili et al. 2010;Zaqarashvili 2018;Dikpati et al. 2018).
We start by recognizing that normalized associated Legendre polynomials Pℓm (x) are identical to the weighted and normalized ultraspherical (Gegenbauer) polynomials (1 − x 2 ) m/2 Ĉm+1/2 ℓ−m (x).One way to understand this correspondence is that the Legendre polynomial P ℓ (x) is equal to the ultraspherical polynomial C 1/2 ℓ (x).The m-th derivative of C 1/2 ℓ (x) is equal to C m+1/2 ℓ−m (x) within a normalization factor.Since the associated Legendre polynomial P ℓm (x) is proportional to the m-th derivative of the Legendre polynomials P ℓ (x) weighted by (1 − x 2 ) m/2 , the correspondence between the former and the weighted ultraspherical polynomials is naturally established.The matrix element O ℓℓ ′ ,m may therefore be obtained by evaluating the corresponding matrix element in a basis of the weighted ultraspherical polynomials instead.The latter permits fast orthogonal polynomial transforms (Olver & Townsend 2013), and an algorithm to obtain the coefficients in the weighted ultraspherical basis is already implemented in the Julia package ApproxFun.jl.We borrow the results directly in our application.

Figure 1 .
Figure 1.Left: Profile of the smoothed solar rotation velocity Ω(r, θ) used in the analysis.Top right: Equatorial cross-section of the differential rotation (Ω(r, θ = π/2) − Ω0)/2π.The smoothed profile (black) is used in the analysis.For reference, the dashed gray line depicts the rescaled solar rotation profile.Bottom right: the rotation profile used by (Bekki et al. 2022b), and the corresponding solar rotation profile truncated to the radial domain [0.71R⊙, 0.985R⊙] that is used by the authors.The reference rate Ω0 in the bottom right panel is chosen to be the Carrington rate 2π × 456 nHz to match that in the work.

Figure 2 .
Figure 2. Sectoral power spectrum (ℓ = m) of the radial vorticity, two different helioseismic techniques.Left panel: Radial vorticity spectrum computed using the global Mode Coupling Analysis (MCA), with even azimuthal orders (Mandal et al. 2021), on HMI data (2010-2020).Right panel: Radial vorticity spectrum computed through ring diagram analysis and averaging all depths, e.g.Hanson et al. (2022).The red line represents the sectoral Rossby dispersion in the thin shell approximation (−2Ω0/(m + 1)), while the blue dashed line shows the maximum azimuthal order to which previous studies have characterized the modes.

Figure 3 .
Figure 3.Comparison of the HMI MCA (black) and HMI RDA (red) radial vorticity power spectra for the same 10-year observation period (2010-2020).The harmonic degree ℓ and azimuthal order m are shown in the top right of each panel.Lorentzian-like ridges are characterized as sectoral Rossby-Haurwitz waves.The expected mode frequency in the thin-shell approximation (ω = −2Ω0/(m + 1)) is shown by the green dashed lines.Lorentzian fits to the HMI RDA are shown by the cyan line.

Figure 4 .
Figure 4. Plot of the mode fit parameters shown in Table1.We show the fits for the HMI RDA (red), HMI MCA (black) and for reference the measurements from time-distance power spectra(Liang et al. 2019, blue).Top panel: Difference between the measured mode frequency and the thin-shell dispersion relation.Middle panel: Measured line widths.Bottom panel: signal-tonoise ratio of the modes measured from the spectra.The shaded area shows the 68% confidence level.

Figure 5 .
Figure 5. Top left panel: Spectrum of latitudinally symmetric solar inertial modes, obtained numerically (this work, white circles) compared with those from solar measurements.The error bars represent measured mode frequencies on the Sun reported by Liang et al. (2019); Proxauf et al. (2020); Hanson et al. (2020), as well as new measurements in this work (black).Top right panel: Same as the left panel, with the frequencies corresponding to the Rossby-Haurwitz dispersion relation for sectoral modes (ω = −2Ω0/(m + 1)) subtracted from the eigenvalues.Bottom left panel: Zoom into the top right panel, focussing on the observed modes.In each panel, only a small section of modes that are the least-damped and have smooth eigenfunctions have been plotted.

Figure 6 .
Figure 6.Full-width at half maximum of observed Rossby modes on the Sun (Proxauf et al. (2020), error bars), and the doubled imaginary components of the numerically obtained eigenfrequencies for various ridges in Figure 5.

Figure 7 .Figure 8 .Figure 9 .Figure 10 .
Figure 7. Real parts of the normalized toroidal stream function V from "ridge 1" for various values of m.

Figure 11 .Figure 12 .
Figure 11.Spectra of symmetric equatorial Rossby modes obtained by scaling the differential rotation profile of the Sun by a constant factor.The panels correspond to different ridges, and the dots joined by lines in each subplot correspond to different rotation rates.The filled patches indicate the range in m beyond which the sectoral Rossby modes are strongly impacted by differential rotation.

Table 1 .
Mode fit parameters of the Rossby-Haurwitz Sectoral modes measured from the HMI MCA and RDA radial vorticity power spectrum.