Viscous inertial modes on a differentially rotating sphere: Comparison with solar observations

In a previous paper we studied the effect of latitudinal rotation on solar equatorial Rossby modes in the beta-plane approximation. Since then, a rich spectrum of inertial modes has been observed on the Sun, which is not limited to the equatorial Rossby modes and includes high-latitude modes. Here we extend the computation of toroidal modes in 2D to spherical geometry, using realistic solar differential rotation and including viscous damping. The aim is to compare the computed mode spectra with the observations and to study mode stability. At fixed radius, we solve the eigenvalue problem numerically using a spherical harmonics decomposition of the velocity stream function. Due to the presence of viscous critical layers, the spectrum consists of four different families: Rossby modes, high-latitude modes, critical-latitude modes, and strongly damped modes. For each longitudinal wavenumber m<4, up to three Rossby-like modes are present on the sphere, in contrast to the equatorial beta plane where only the equatorial Rossby mode is present. The least damped modes in the model have eigenfrequencies and eigenfunctions that resemble the observed modes; the comparison improves when the radius is taken in the lower half of the convection zone. For radii above 0.75R and Ekman numbers E<10^{-4}, at least one mode is unstable. For either m=1 or m=2, up to two Rossby modes are unstable when the radial dependence of the Ekman number follows a quenched diffusivity model (E=2. 10^{-5} at the base of the convection zone). For m=3, up to two Rossby modes can be unstable, including the equatorial Rossby mode. Although the 2D model discussed here is highly simplified, the spectrum of toroidal modes appears to include many of the observed solar inertial modes. The self-excited modes in the model have frequencies close to those of the observed modes with the largest amplitudes.


Introduction
reported the observation and identification of a large set of global modes of solar oscillations in the inertial frequency range. These include the equatorial Rossby modes (Löptien et al. 2018), high-latitude modes (the m = 1 symmetric mode manifests itself as a spiral structure, see Hathaway et al. 2013;Bogart et al. 2015), and critical-latitude modes. The observed modes were identified by comparison with (a) inertial modes computed for an axisymmetric model of the convection zone (using a 2D solver) and (b) purely toroidal modes computed on the solar surface (using a 1D solver). In this paper, we provide additional results based on the 1D solver. Additional results based on the 2D solver are provided in a companion paper (Bekki et al. 2022).
Under a simplified setup in the equatorial β plane, Gizon et al. (2020) discussed the importance of latitudinal differential rotation in this problem. The eigenvalue problem for purely toroidal modes is singular in the absence of viscosity. In the presence of eddy viscosity, the eigenmodes can be grouped into several families of modes: the critical-latitude modes, the high-latitude modes, the strongly damped modes, and the equatorial Rossby modes. These modes arise from the existence of viscous critical layers where the wave speed is equal to the zonal shear flow. In particular, the equatorial Rossby modes (the R modes) are trapped below the critical layer. In the present paper, we wish to extend the work of Gizon et al. (2020) to the study of toroidal modes on the sphere (i.e. we drop the equatorial β-plane approximation) in the presence of realistic latitudinal differential rotation and eddy viscosity. The intention is to include critical and high-latitude modes in the discussion: can these modes also be captured by simple physics on the sphere?
Instabilities and waves on a differentially rotating sphere have been studied before. The basic equation for the oscillations (Sect. 2) was established (in an inertial frame) by Watson (1981) and later extended to more realistic solar rotation profiles by Dziembowski & Kosovichev (1987) and Charbonneau et al. (1999). All these studies showed the possibility of unstable modes when the differential rotation is strong enough, which is the case for the Sun in the upper convection zone. Unstable modes also appear in the overshooting part of the tachocline when considering a shallow-water model allowing for radial motions instead of purely toroidal modes (Dikpati & Gilman 2001).
Here, we restrict ourselves to the 2D approximation which is valid for strongly stratified media when the rotation rate is much smaller than the buoyancy frequency (Watson 1981). For the Sun, it is valid in the radiative zone, the A&A proofs: manuscript no. manuscript_arxiv atmosphere, and the outer part of the convective envelope (Dziembowski & Kosovichev 1987). Moreover, the 2D approximation is a valuable approximation to the 3D problem for the modes that vary slowly with depth (Kitchatinov & Rüdiger 2009) and are away from the axis of symmetry (Rieutord et al. 2002).
In this work we study only the hydrodynamical modes and do not consider the influence of a magnetic field, even if it was shown to be a significant ingredient (Gilman & Dikpati 2000Gilman et al. 2007). The main difference with previous works is the introduction of a viscous term which is important to understand the shape of the eigenfunctions and the lifetime of the subcritical modes. Also we do not restrict latitudinal differential rotation to a two-or four-term profile, but consider the solar rotation profile as measured by helioseismology at the solar surface (Sect. 4) and in the interior (Sect. 5). We study the spectrum and the eigenfunctions of the normal modes of the system. Depending on the parameters of the problem, we find that some modes may be self-excited (Sect. 5). The conclusion, Sect. 6, includes a short discussion on the stability of differential rotation for distant stars.

Equation for the velocity stream function
We study the propagation of purely toroidal modes on a sphere of radius r under the influence of latitudinal differential rotation. In an inertial frame, the rotation rate Ω(θ) depends on the co-latitude θ measured from the rotation axisẑ. We work in a frame rotating at the reference angular velocity Ω ref (chosen to be the Carrington rate later in the paper). In the rotating frame, the Navier-Stokes equation is whereẑ is the unit vector along the rotation axis, u the horizontal velocity, and ν the eddy viscosity. For the sake of simplicity, wave damping by turbulence is modelled by a horizontal Laplacian ∆. The force on the right-hand side is assumed to derive from a potential Π. In the rotating frame, we decompose the velocity into the mean axisymmetric flow U and the wave velocity where φ is the longitude (increases in the prograde direction). For example, U = (Ω − Ω ref )ẑ × r is the flow associated with latitudinal differential rotation. To first order in the wave amplitude, we have The two horizontal components of this equation are where is the material derivative in the rotating frame. For purely toroidal modes, we can introduce the stream function Ψ(θ, φ, t) such that Combining Eq. (4) and Eq. (5), we obtain the equation of Watson (1981) modified on the right-hand side to include viscosity: where (Ω sin 2 θ) .

Modal decomposition
We look for wave solutions of the form where m is the longitudinal wavenumber and ω is the (complex) angular frequency measured in the rotating reference frame. Introducing the Ekman number the function ψ satisfies where we defined the relative differential rotation and the operator L m such that i.e. the horizontal Laplacian on the unit sphere (also called the associated Legendre operator). Equation (12) with E = 0 reduces to Watson (1981)'s equation when written in a rotating frame. However, when E 0, Eq. (12) is fourth-order, which has profound implications for the spectrum. Four boundary conditions are required. We impose that the flow vanishes at the poles:

Numerical solutions
In the inviscid case (E = 0 in Eq. 12), critical latitudes θ c appear where mδ( where The stream eigenfunctions are continuous but not regular at θ c , thus u ′ is singular there (see Gizon et al. 2020, for the β-plane problem). When including viscosity, the critical latitude is replaced by a viscous layer whose width is proportional to E 1/3 . The eigenfunctions are now regular. The solutions can be expended onto 0 30 60 90 120 150 180 where the P m ℓ are normalized such that 1 −1 [P m ℓ (µ)] 2 dµ = 1. We insert this expansion into Eq. (12) and project onto a particular polynomial P m ℓ . Using L m P m ℓ (cos θ) = −ℓ(ℓ + 1)P m ℓ (cos θ), we obtain the following matrix equation: where b = [b m b m+1 · · · b L ] T is a vector and the matrix C has elements Equation (19) defines an eigenvalue problem where ω is the eigenvalue and b the associated eigenvector. The differential rotation (through δ and ζ) couples the different values of ℓ and ℓ ′ , so that the matrix C is not diagonal (the eigenfunctions are not P m ℓ ). Under the assumption that the rotation profile is symmetric about the equator, the problem decouples into odd and even values of ℓ and has to be solved separately for the symmetric and antisymmetric eigenfunctions. Throughout this paper, we call symmetric the modes that have a north-south symmetric stream function, i.e. that are symmetric in u ′ θ and antisymmetric in u ′ φ . The modes with a north-south antisymmetric stream function are called antisymmetric. For the sake of simplicity, we do not consider the case of a general rotation profile with North-South asymmetries. This would result in eigenfunctions that are not symmetric nor antisymmetric.
Thanks to the decomposition given by Eq. (17) the boundary conditions, Eqs. (15), are automatically satisfied for m > 1. For m = 1, however, the boundary condition on the derivative is not satisfied because dP 1 ℓ (cos θ)/dθ 0 at the poles. The modifications for this case are discussed in Appendix A.

Input parameters: viscosity and rotation profile
The only input quantities required to solve the eigenvalue problem are the viscosity and the rotation profile. In this paper, we choose an Ekman number at the surface E = 4 × 10 −4 , which means ν = 125 km 2 /s and an associated Reynolds number of 300. This choice leads to a good match with the observed eigenfunctions of the equatorial Rossby modes . The influence of the viscosity on the spectrum and the eigenfunctions is studied in Sect. 4.1.
Regarding rotation, we use the profile inferred by helioseismology from Larson & Schou (2018). Since the first and second derivatives of Ω(r, θ) with respect to θ are required to compute ζ, we first write the profile as a truncated series of harmonic functions to smooth it: For each r, the coefficients Ω ℓ , 0 ≤ ℓ ≤ N are found by fitting the observations taking the random errors into account. These coefficients change with the solar cycle (torsional oscillations). For north-south symmetric rotation, the odd coefficients are zero. Rotation from global-mode helioseismology is north-south symmetric by construction, however general rotation profiles could be used using the methods of the present paper. Note that the coefficients in this decomposition depend on N because the powers of cos θ are not orthogonal. The (symmetric) rotation profile from Larson & Schou (2018)
High-latitude Modes whose eigenfunctions have largest amplitudes in the polar regions. Their frequencies are the most negative in the rotating frame. The least damped mode at fixed m has frequency ω ≈ mΩ 2 .
Critical-latitude Modes whose eigenfunctions have the largest amplitudes at mid-or low-latitudes, between their critical latitudes. Their frequencies are the smallest in absolute value, with Re[ω] ≈ Im[ω].

Strongly damped
Modes with very large attenuation (|Im[ω]| ≥ |Re[ω]|), whose eigenfunctions are highly oscillatory around their critical latitudes and frequencies satisfy Re[ω] ≈ mΩ 2 /2 (for the β-plane analogy, see Grosch & Salwen 1968). Table 2: Description of the viscous modes discussed in this paper, for the case Ω = Ω 0 +Ω 2 cos 2 θ. The high-latitude, critical-latitude, and strongly damped modes owe their existence to the presence of viscous critical layers. Rossby modes would exist even in the case of vanishing differential rotation (i.e. without critical layers).
of Fig. 1. Close to the surface, an expansion with N up to 30 is required to capture the sharp change of slope of the profile around 60 • latitude, as well as the solar-cycle related zonal flows. Adding more terms would not change the surface profile significantly. Below the near-surface shear layer, an expansion with N = 4 is generally sufficient to capture the rotation profile. The quantity ζ, which involves the first and second derivatives of Ω(θ; r, N) with respect to θ shows strong variations near the surface (right panel of Fig. 1). These will have an important effect on the spectrum as shown in the next section.

Uniform rotation
In the case of uniform rotation, Ω = Ω 0 , the complex eigenvalues can be obtained directly from Eq. (19) (the matrix C is diagonal because δ and ζ are constant): The modes are classical Rossby modes with eigenfunctions P m ℓ (cos θ). The modes are stable as their eigenvalues all have negative imaginary parts. For a given m, the sectoral Rossby mode (ℓ = m) is the least damped mode.

Two-term differential rotation and E 10 −4
Here we solve the eigenvalue problem given by Eq. (19) using Ω 0 + Ω 2 cos 2 θ as an approximation to the solar differential rotation at the surface. Fitting the observed surface rotation profile, we use Ω 0 /2π = 452 nHz and Ω 2 /2π = −114 nHz. This simple profile allows to make the connection with the study of Gizon et al. (2020), where a parabolic flow was used in the β−plane, and serves as a reference to track and identify the modes when the rotation profile is steeper and more solar-like. We use an Ekman number E = 4 × 10 −4 which is close to the value of the eddy viscosity at the solar surface. As shown in Fig. 2 (top panels), the spectrum takes a Y-shape in the complex plane. This spectrum is characteristic of the spectrum of the Orr-Sommerfeld equation for a parabolic (Poiseuille) flow. The three branches of the spectrum have been called by Mack (1976, his figure 5) the P family (after Pekeris 1948), the S family (after Schensted 1961), and the A family; they correspond respectively to the 'center modes', the 'damped modes', and the 'wall modes'. We refer to the monographs by Drazin & Reid (1981) and Schmid & Henningson (2012) for more details about this problem in the context of hydrodynamics, as well as to the review paper by Maslowe (2003) for discussions about critical layers in shear flows. The correspondence between the terminologies used in hydrodynamics and in solar physics is given in Table 1. We find that the three branches are still present for low values of m, even though the β-plane approximation is no longer justified. We label the symmetric modes with an even index and the antisymmetric modes with an odd index. They are ordered on the branches such that the index increases with the attenuation (see Fig. 2, top panel for m = 2).
In addition to the above modes, up to three other modes are present in the spectrum at low m values. These modes are outside the three branches of the spectrum and are strongly affected by the Coriolis force. These are Rossby modes. For m large enough, only the equatorial Rossby mode (denoted by the letter R) is easily identifiable outside the branches (as in the β−plane approximation, see Gizon et al. 2020, their figure 4). We identify up to two additional Rossby modes in the frequency spectra (denoted R 1 and R 2 ). These can be traced back to the traditional Rossby modes in the case of uniform rotation. By progressively turning on the differential rotation, i.e. Ω(θ) = Ω 0 + ǫΩ 2 cos 2 θ with ǫ increasing from 0 to 1, we find that these correspond to the modes with ℓ = m + 1 and ℓ = m + 2 in the case of uniform rotation (ǫ = 0), see the online movies for m = 1 and m = 2. The study of Rossby waves in stars and more generally in astrophysics is a broad topic. We refer the reader to the recent review by Zaqarashvili et al. (2021). Basic properties of the modes studied in this paper are summarized in Table 2. 3.3. Two-term differential rotation and E 10 −4 While the value of the viscosity is well constrained at the solar surface, its value inside the Sun remains largely unknown. Several models presume that the viscosity should decrease with depth (see e.g. Eigenfrequencies ω for m = 1 to 3 (from left to right). The Ekman number is E = 4×10 −4 . The top panels are for Ω(N = 2) = Ω 0 + Ω 2 cos 2 θ, with Ω 0 /2π = 452 nHz and Ω 2 /2π = −114 nHz and the bottom panels for the solar rotation rate Ω(N = 30). Modes with north-south symmetric eigenfunctions are shown with filled symbols, those with antisymmetric eigenfunctions are shown with open symbols. The frequencies are expressed in the Carrington frame of reference. We note that the imaginary part of the frequency is connected to the mode linewidth Γ through Γ = −2ω i . The top axis of each plot marks the frequencies corresponding to the critical latitudes 30 • , 50 • , and 70 • (this is not a statement about the latitude of maximum power). The three branches for the A, P, and S families are highlighted with different colors (see Mack 1976, for a comparison with the spectrum of the Orr-Sommerfeld equation). The modes denoted R, R 1 , and R 2 are closest to the traditional r-modes with ℓ = m, ℓ = m + 1, and ℓ = m + 2 respectively (see, e.g., Saio 1982). The denomination of the modes in the solar case (bottom panels) is obtained by tracking the modes from the N = 2 case (top panels). See the evolution between the two profiles on the online movies for m = 1 and m = 2. of 4 × 10 −4 as in Fig. 2). For m = 2 and m = 3, the spectrum is not Y-shaped but more complex as many S modes are now situated on a plateau with nearly constant imaginary part. The two branches corresponding to the A and P families are still clearly visible.

Spectrum for solar rotation at the surface
When we use a very good approximation to the solar differential rotation at the surface (N = 30 terms in the expansion to describe Ω), we find that the Y-shape of the spectrum is much harder to identify and that many modes move away from the branches, see bottom panels in Fig. 2 for 1 ≤ m ≤ 3. To label the modes, we track their frequencies in the complex plane by slowly transition-ing from the case Ω(N = 2) to the solar case Ω(N = 30), see, e.g., the online movie for m = 2.
The real parts of the eigenfrequencies of the modes R, R 1 and R 2 do not change very much compared to the case of the twoterm rotation profile. However the imaginary parts of the eigenfrequencies do change. In particular, those of the (symmetric) R 2 mode with m = 1 and the (antisymmetric) R 1 mode with m = 2 are now positive, i.e. these modes are unstable (self-excited). The second derivative of the rotation profile has a strong influence on the damping of the modes, especially the high-latitude and critical-latitude modes (see, e.g, online movies). While the damping rates of the A and P modes are large for Ω(N = 2), some of these modes become less damped for the solar-like Ω(N = 30) and are thus more likely to be observable in the solar data.  The eigenfunctions of the high-latitude modes A 1 and A 2 have a modulus that peaks near the critical layers at ≈ ±70 • and the argument between the real and imaginary parts varies below this latitudes (which implies a spiral pattern there). The eigenfunctions depend sensitively on the rotation profile. The computed frequencies of these modes (−224 nHz) are not very far from the observed frequencies at −171 nHz and −151 nHz (Gizon et al. 2021, their table 1). We will show in Sect. 5.2, that an improved match can be obtained when the solar rotation profile is taken deep in the convection zone.
The critical-latitude modes P 1 and P 2 have eigenfunctions that oscillate around and below the critical latitudes. Their real and imaginary parts are not in phase as functions of latitude and have roughly the same amplitude. These modes are not very different from the case of the two-term rotation profile. They resemble the observed m = 2 critical-latitude modes at frequencies −12 nHz and −24 nHz reported by Gizon et al. (2021).

Effect of turbulent viscosity
In the uniform rotation case, the viscosity only influences the imaginary part of the eigenfrequencies as shown by Eq. (22). However, the picture changes when differential rotation is included. The left panel of Fig. 3 shows the eigenfrequencies of several modes with m = 2 as a function of the Ekman number when the surface solar differential rotation profile is considered. The imaginary part goes to zero as the viscosity tends to 0 and increases drastically when the Ekman number becomes larger than 10 −3 . More surprisingly, the real part is also significantly affected for modes A 1 and P 1 , by an amount of more than 100 nHz depending on the value of the viscosity. The R mode is little affected for this small value of m. For larger values of m, we observe a change in the R-mode frequency of, e.g., 15 nHz for m = 6 and 30 nHz for m = 10 over the range of Ekman numbers covered by the plot. All these variations would be measurable in the solar observations and thus the modes may be used as probes of the viscosity in the solar interior. The eigenfunctions are also affected (see Fig. 3) in particular the imaginary part of the mode A 1 changes sign with viscosity and the mode P 1 becomes more and more confined between the viscous layers as viscosity decreases. The shape of the R-mode eigenfunction does not depend significantly on the viscosity for this small value of m, unlike for larger values of m as already discussed by Gizon et al. (2020).

Effect of meridional flow
In order to study the importance of the meridional flow V(θ), we consider the background flow At the surface, we take a poleward flow with maximum amplitude V max = 15 m/s, V(θ) = −1.54 × V max sin 2 θ sin 2θ.
The derivation and the discretization of this problem are presented in Appendix B. The stream function ψ satisfies Eq. (B.5). For m 5, the effect of the meridional flow on the R mode was already discussed by Gizon et al. (2020) in the β plane. Figure C.4 shows the entire frequency spectrum with and without the meridional flow in spherical geometry. For m = 2, the eigenfrequencies of the modes R and A 1 are shifted by only a few nanohertz and their eigenfunctions are not significantly affected (see Fig. C.5, left and middle panels). An interesting effect caused by the meridional flow is the change of the imaginary part of the frequency of the R 1 mode. This mode becomes even more unstable as a consequence of the meridional flow.
We find that the critical-latitude modes are the most affected by the meridional flow. This is not surprising as their eigenfunctions vary fast at mid-latitudes where the meridional flow is largest. The real parts of the eigenfrequencies of the m = 2 modes located on the P branch are getting closer to zero when the meridional flow is included, with a shift of ∼ 10 nHz. The The regions corresponding to the most unstable modes (R 2 , R 1 , R) and to the two most unstable modes (R 2 and R 1 , R 1 and R 2 , R and R 1 ) are colored and labelled with these modes. The color bar indicates the total number of unstable modes at any point in the diagram. meridional flow also stretches the eigenfunctions in latitude, see, e.g., the right panel of Fig. C.5 for the P 1 mode for m = 2.

Effect of time-varying zonal flows
To assess the sensitivity of the modes to the details of the rotation profile, we compute the mode frequencies as a function of time over the last two solar cycles. We use the inferred rotation profile from helioseseismology (Larson & Schou 2018), obtained at a cadence of 72 days. The rotation profiles are averaged over five bins, to obtain yearly averages during 1996-2018. The time variations of the frequencies of the least-damped modes with m = 2 are plotted in Fig. C.6. We find that the frequency of the R mode varies by less than 1 nHz, in agreement with the calculation of Goddard et al. (2020) using first-order perturbation theory. The critial-latitude mode P 1 also changes by a very small amount of ±2 nHz. The frequency of the high-latitude mode A 1 varies by up to ±7 nHz and shows a 22-year periodicity. Other modes have frequencies that strongly vary (such as P 5 , not shown on the plot, which varies by up to ±10 nHz). Since the frequency resolution corresponding to a 3 year time series is ≈ 10 nHz and the typical linewidth of a mode is also ≈ 10 nHz, the frequency shifts of some of the modes may be detectable in observed time series.

Spectrum for solar rotation at different radii
In the previous section we only considered the rotation profile on the solar surface. Here, we consider latitudinal differential rotation at different depths in the solar interior. We study mode stability as a function of depth, and how the dispersion relations of the different modes are affected.

Stability analysis
The hydrodynamical stability of solar latitudinal differential rotation was studied in two dimensions for purely toroidal distur-bances in the inviscid case. A necessary condition for instability is that the latitudinal gradient of vorticity (ζ, see Eq. (9)) must change sign at least once with latitude (Lord Rayleigh 1879). Fjørtoft (1950) proved a more restrictive condition for instability: there exist a θ and a θ 0 (where θ 0 is a zero of ζ) such that [Ω(cos θ) − Ω(cos θ 0 )]ζ(θ) < 0. As seen in Fig. 1b, the function ζ(θ) switches sign multiple times at the surface, and thus unstable modes could exist. In the special case of a rotation law Ω = Ω 0 + Ω 2 cos 2 θ, Watson (1981) showed that a necessary condition for stability is −(2/7)Ω 0 < Ω 2 < 1.14Ω 0 . This condition is met for two-term fits to the solar rotation profile. Charbonneau et al. (1999) added a third term in Ω 4 cos 4 θ and found numerically that when Ω 4 ≈ Ω 2 ≈ −0.1Ω 0 (as is the case for the Sun) then two modes (one symmetric and one antisymmetric) are unstable for each of the cases m = 1 and m = 2. They considered different solar rotation profiles inferred from the LOWL, GONG, and MDI p-mode splittings and performed a stability analysis on spheres at different depths. They found that all modes are stable below ≈ 0.74R ⊙ while some modes with m = 1 and m = 2 become unstable in the upper convection zone.
We extended the stability analysis of Charbonneau et al. (1999) by including viscosity and using the latest rotation profile from helioseismology (Larson & Schou 2018). Figure 4 shows the number of unstable modes as a function of the Ekman number and depth for m = 1, 2 and 3. For m ≥ 4, all modes are stable. Several modes can be unstable, with up to 12 modes in a small layer below the solar surface for E < 10 −4 (5 modes with m = 1, 5 modes with m = 2, and 2 modes with m = 3). We find that for E < 10 −4 , the radius above which some modes become unstable (≈ 0.75R ⊙ for m = 1, ≈ 0.77R ⊙ for m = 2, and ≈ 0.78R ⊙ for m = 3) does not depend sensitively on viscosity and is almost given by the E = 0 limit. Only values E > 10 −4 have a stabilising effect. For E > 5 × 10 −3 all depths are stable. On the same figure we draw estimates of the Ekman number with depth, either under mixing-length theory or using a quenched diffusivity (Muñoz-Jaramillo et al. 2011). Using the value E η at each depth, six modes are unstable, two for each value of m. These are the (b) r = 0.75R ⊙ , E = 2 × 10 −5 Fig. 5: Model dispersion relations for all the modes with −Im[ω]/2π < 100 nHz. Shown are the equatorial Rossby modes (black curves for R, R 1 and R 2 ), the high-latitude modes (red area), the strongly damped modes (blue area), and the critical-latitude modes (orange area). The differential rotation is solar (N = 30), both at the surface (r = R ⊙ and E = 4 × 10 −4 , left panel) and at the bottom of the convection zone (r = 0.75R ⊙ and E = 2 × 10 −5 , right panel). Frequencies are expressed in the Carrington frame of reference. The symbols show the observed modes reported by Gizon et al. (2021).
modes R 1 and R 2 for m = 1 and m = 2, and R and R 1 for m = 3. Their eigenfunctions at the surface are shown in Fig. C.7. Since the critical latitude is very close to the pole, it does not affect much the shape of the eigenfunctions which resemble the classical P m ℓ with ℓ = m, m + 1, m + 2 obtained in the case of uniform rotation. Using the values of E MLT , only three modes are unstable (R 1 and R 2 for m = 1, and R 1 for m = 2) and no modes are unstable below 0.91R ⊙ instead of 0.75R ⊙ using E η .
The amplitudes of the inertial modes may thus have diagnostic value to learn about eddy viscosity in the solar interior. Figure 5 shows the propagation diagram for different modes when the radius is either r = R ⊙ or r = 0.75R ⊙ . We consider all the modes in the model whose imaginary frequencies are less than 100 nHz (in absolute value) and which would thus be easier to detect on the Sun. Since several modes of each family are present for each value of m, we draw areas in the propagation diagram where these modes are present. The observed frequencies reported by Gizon et al. (2021) are overplotted. We see that observed high-latitude modes have frequencies that overlap with the region between the slowest and fastest high-latitude modes of the 2D model and are best approximated when the rotation profile is taken at the base of the convection zone. This is not surprising as these modes have most of their kinetic energy there according to 3D modeling (see Gizon et al. 2021;Bekki et al. 2022). The observed high-latitude modes for m = 1, 2, and 3 have frequencies close to the modes R 1 and R 2 that are selfexcited. This may explain why the observed high-latitude modes have the largest amplitudes. However, the eigenfunctions from the model differ from the observations. It is thus difficult to clearly identify these modes among the modes A 1 , A 2 , R 1 , and R 2 . The observed critical-latitude modes with frequencies between −150 nHz and 0 nHz may be associated with the dense spectrum of critical-latitude modes. However, the model does not have any mode with positive frequency (in the Carrington frame). The model is also unable to explain some of the observed critical-latitude modes with Re[ω]/2π < −200 nHz (for example the modes with m = 9 and 10 and frequencies around −280 nHz).  For this low m value, the R-mode eigenfunction changes little with depth. This suggests that the 2D problem captures the essential physics of R modes and that the assumption that these modes are quasi-toroidal is well justified. This conclusion is confirmed by solving the 3D eigenvalue problem in a spherical shell: Gizon et al. (2021) and Bekki et al. (2022) find that the radial velocity of R modes is about two orders of magnitude smaller than their horizontal velocity. On the other hand, the high-latitude mode A 1 and (to a lesser extent) the critical-latitude mode P 1 have eigenfunctions that change significantly between the surface and the base of the convection zone (middle panels in Fig. C.3). In the case of the high-latitude modes, this is not too surprising (the 2D approximation is poor close to the axis of symmetry, see, e.g., Rieutord et al. 2002). A better understanding of the high-latitude modes requires not only a 3D geometrical setup, but also the inclusion of a latitudinal entropy gradient (Gizon et al. 2021;Bekki et al. 2022).

Discussion
We extended the work of Gizon et al. (2020) in the equatorial β plane to a spherical geometry, in order to study the effects of differential rotation and viscosity on the modes with the lowest m values. Like in the equatorial β plane, viscous critical layers appear and the spectrum contains different families of modes due to the latitudinal differential rotation (Fig. 2). The high-latitude, critical-latitude, and strongly damped modes (Mack 1976) are still present when a realistic solar rotation profile is used. Due to the Coriolis force, the Rossby equatorial modes (R modes) are also present in this problem, like in the β plane problem. Using the surface solar rotation profile, two additional Rossby modes, R 1 (for m ≤ 6) and R 2 (for m ≤ 2), are also present. Remarkably, the calculated spectrum resembles closely that of the inertial modes observed on the Sun by Gizon et al. (2021), especially when the differential rotation is taken deep in the convection zone, see Fig. 5b. The frequencies of some modes in the model are sensitive to the value of the viscosity (E MLT versus E η ) and could be a good probe of this parameter (Fig. 3, left panel). We also find that the zonal flows (Fig. C.6) and the meridional flow (Fig. C.4) can affect the modes. The 2D model presented here is useful to understand the basic physics of toroidal modes on the Sun, however a 2D model cannot replace a more sophisticated 3D model that includes radial motions, realistic stratification, superadiabaticity (i.e. the strength of the convective driving), and entropy gradients (for a treatment of these effects, see Bekki et al. 2022).
We confirm the result by Charbonneau et al. (1999) that some modes are self-excited, even when viscosity and realistic solar differential rotation are included. The angular velocity terms in Ω 2p cos 2p θ with p ≥ 2 are essential to destabilize the system. For the Sun, an oversimplified fit of the form Ω 0 + Ω 2 cos 2 θ would lead to the wrong conclusion that the system is stable. Using a value of the viscosity corresponding to the quenched diffusivity model from Muñoz-Jaramillo et al. (2011), we find that six modes are self-excited. These are the Rossby modes R 1 and R 2 for m = 1 and m = 2, and R and R 1 for m = 3. If the modes R 1 and R 2 correspond to the high-latitude modes observed by Gizon et al. (2021), it is understandable why they should have the largest amplitudes. Above r = 0.78R ⊙ , the m = 3 equatorial R mode is unstable; this mode is the R mode with the lowest m value that is observed on the Sun. The excitation and damping by turbulent convection of the subcritical modes will be studied in upcoming papers.
In the case of distant stars, latitudinal differential rotation is detectable with asteroseismology only for a few Kepler Sun-like stars (Benomar et al. 2018). For these stars, the Ω 2 cos 2 θ term is very large compared to the equatorial value Ω 0 , and they are unstable according to Watson's criterion (Table 3). Unfortunately we do not have any information about the high-order terms that fully determine the rotation profile of these stars.
A&A proofs: manuscript no. manuscript_arxiv Appendix A: Eigenvalue problem for m = 1 In order to ensure that ψ satisfies the boundary conditions (Eq. (15)) for m = 1, we expand the solution on a basis of associated Legendre polynomials P 2 ℓ instead of P 1 ℓ : Next, we need to specify how the operators L 1 and (L 1 ) 2 act on P 2 ℓ (cos θ): L 2 1 P 2 ℓ (cos θ) = − 12 cos θ sin 3 θ ∂ θ P 2 ℓ (cos θ) The eigenvalue problem, Eq. (12), becomes (20) with m = 1, and the elements of the matrices A, D and G are given by: modes for m = 2, using the solar surface differential rotation profile and the Ekman number E = 4 × 10 −4 . The solid curves show the real parts, the dashed curves the imaginary parts. In each case, the normalization is such that ψ = 1 where |ψ| is maximum (this happens at different latitudes depending on the mode). In each panel, the name of the mode and the real part of its frequency are given (see also Fig. 2e). The vertical red arrows give the latitudes of the viscous layers. r = R ⊙ R (-336.5-0.6i nHz) A 1 (-233.0-9.2i nHz) P 1 (-26.1-23.6i nHz) r = 0. 75R ⊙ R (-327.9-0.1i nHz) A 1 (-182.8-2.0i nHz) P 1 (6.5-8.0i nHz)