© D. Fournier et al. 2022

1. Introduction

Gizon et al. (2021) have recently 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) discuss the importance of latitudinal differential rotation for 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 do not consider the equatorial β-plane approximation) in the presence of realistic latitudinal differential rotation and eddy viscosity. The intention is to study if critical and high-latitude modes can 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 of these studies showed the possibility of unstable modes when 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 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 2000, 2002; Gilman et al. 2007). The main difference with previous works is the introduction of a viscous term which is important in order 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 we 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 includes a short discussion on the stability of differential rotation for distant stars (Sect. 6).

2. Linear toroidal modes on a sphere

2.1. 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 colatitude θ measured from the rotation axis . We worked 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

(1)

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 was modeled by a horizontal Laplacian Δ. The force on the right-hand side is assumed to be the derivative of a potential Π. In the rotating frame, we decomposed the velocity into the mean axisymmetric flow U and the wave velocity u′,

(2)

where ϕ is the longitude, which increases in the prograde direction. For example, is the flow associated with latitudinal differential rotation. To first order in the wave amplitude, we have

(3)

The two horizontal components of this equation are

(4)

(5)

where

(6)

is the material derivative in the rotating frame. For purely toroidal modes, we can introduce the stream function Ψ(θ, ϕ, t), such that

(7)

Combining Eqs. (4) and (5), we obtained the equation of Watson (1981) which was modified on the right-hand side to include viscosity:

(8)

where

(9)

2.2. Modal decomposition

We looked for wave solutions of the form

(10)

where m is the longitudinal wavenumber and ω is the (complex) angular frequency measured in the rotating reference frame. Introducing the Ekman number

(11)

the function ψ satisfies

(12)

where we defined the relative differential rotation

(13)

and the operator L_m, such that

(14)

that is the horizontal Laplacian on the unit sphere (also called the associated Legendre operator). Equation (12) with E = 0 reduces to the equation from Watson (1981) 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 imposed the flow to vanish at the poles:

(15)

2.3. Numerical solutions

In the inviscid case (E = 0 in Eq. (12)), critical latitudes θ_c appear where mδ(θ_c)−ω/Ω_ref = 0, that is where

(16)

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 was replaced by a viscous layer whose width is proportional to E^1/3. The eigenfunctions are now regular. The solutions can be expended onto a series of normalized associated Legendre polynomials up to order ℓ = L,

(17)

where the are normalized, such that . We inserted this expansion into Eq. (12) and projected onto a particular polynomial . Using

(18)

we obtained the following matrix equation:

(19)

where b = [b_mb_m + 1… b_L]^T is a vector and matrix C has elements

(20)

Equation (19) defines an eigenvalue problem where ω is the eigenvalue and b is the associated eigenvector. The differential rotation (through δ and ζ) couples the different values of ℓ and ℓ′, so that matrix C is not diagonal (the eigenfunctions are not ). 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, the modes that have a north-south symmetric stream function are referred to as symmetric, that is to say they are symmetric in and antisymmetric in . 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 neither symmetric nor antisymmetric.

Thanks to the decomposition given by Eq. (17), the boundary conditions, as shown in Eq. (15), are automatically satisfied for m > 1. For m = 1, however, the boundary condition on the derivative is not satisfied because at the poles. The modifications for this case are discussed in Appendix A.

2.4. Input parameters: Viscosity and rotation profile

The only input quantities required to solve the eigenvalue problem are the viscosity and the rotation profile. For this paper, we chose an Ekman number at the surface E = 4 × 10⁻⁴, which means ν = 125 km² s⁻¹ and an associated Reynolds number of 300. This choice led to a good match with the observed eigenfunctions of the equatorial Rossby modes (Gizon et al. 2020). 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 wrote the profile as a truncated series of harmonic functions to smooth it:

(21)

For each r, the coefficients Ω_ℓ, 0 ≤ ℓ ≤ N were found by fitting the observations, taking the random errors into account in the minimization. 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 when employing the methods of the present paper. We 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) averaged over 1996–2018 is shown in the left panel of Fig. 1. Close to the surface, an expansion with N up to 30 is required to capture the sharp change in slope for 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 have an important effect on the spectrum as shown in the next section.

Fig. 1. Main quantities entering in the eigenvalue problem defined by Eq. (12). (a) Normalized solar differential rotation rate δ(θ) at different radii, inferred from helioseismology and averaged over 1996 to 2018 (the shaded area covers ±3σ), with respect to the Carrington rotation rate, Ωref = ΩCarr. The solid curves show the N = 30 fits to the solar rotation rates. The dashed curve shows the N = 2 fit to the solar surface rotation rate, i.e., Ω0 + Ω2cos2θ with Ω0/2π = 452 nHz and Ω2/2π = −114 nHz. (b) The function ζ defined by Eq. (9) is plotted for each of the three cases of panel (a). For r = R⊙ and N = 30, the function reaches a minimum at ζ(θ = 0)≈ − 22 (not shown on the plot). The oscillations at the surface are due to the m = 0 zonal flows (also known as torsional oscillations). The horizontal dashed-black line shows , where β = 2Ωref/r, which plays the role of ζ in the equatorial β plane (see Eq. (17) in Gizon et al. 2020).

3. Spectrum for simplified rotation laws

3.1. Uniform rotation

In the case of uniform rotation, Ω = Ω₀, the complex eigenvalues can be obtained directly from Eq. (19) (matrix C is diagonal because δ and ζ are constant):

(22)

The modes are classical Rossby modes with eigenfunctions (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.

3.2. Two-term differential rotation and E ≳ 10⁻⁴

Here we solve the eigenvalue problem given by Eq. (19) using Ω₀ + Ω₂cos²θ as an approximation for the solar differential rotation at the surface. Fitting the observed surface rotation profile, we used Ω₀/2π = 452 nHz and Ω₂/2π = −114 nHz. This simple profile allows one 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 used an Ekman number E = 4 × 10⁻⁴ which is close to the value for the eddy viscosity at the solar surface. As shown in Fig. 2 (top panels), the spectrum has 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 Fig. 5) the P family (after Pekeris 1948), the S family (after Schensted 1961), and the A family; they correspond to the “center modes”, the “damped modes”, and the “wall modes”, respectively. We refer the interested reader 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).

Fig. 2. Eigenfrequencies ω for m = 1–3 (from left to right). The Ekman number is E = 4 × 10−4. The top panels are for Ω(N = 2) = Ω0 + Ω2cos2θ, with Ω0/2π = 452 nHz and Ω2/2π = −114 nHz, and the bottom panels are for the solar rotation rate Ω(N = 30). Modes with north-south symmetric eigenfunctions are shown with filled symbols, and 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, R1, and R2 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) was obtained by tracking the modes from the N = 2 case (top panels). See the evolution between the two profiles in the online movies for m = 1 and 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. They are Rossby modes. In the case where m is 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 Fig. 4). We identified up to two additional Rossby modes in the frequency spectra (denoted R₁ and R₂). These can be traced back to the traditional Rossby modes in the case of uniform rotation. By progressively turning on the differential rotation, that is to say Ω(θ) = Ω₀ + ϵΩ₂cos²θ with ϵ increasing from 0 to 1, we found 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). The basic properties of the modes studied in this paper are summarized in Table 2.

3.3. Two-term differential rotation and E ≲ 10⁻⁴

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., Fig. 1 of Muñoz-Jaramillo et al. 2011). The top panels in Fig. C.1 show the spectrum of the modes m = 1, 2, and 3 when the Ekman number is set to E = 2 × 10⁻⁵ (instead of 4 × 10⁻⁴ 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 a nearly constant imaginary part. The two branches corresponding to the A and P families are still clearly visible.

4. Spectrum for solar rotation at the surface

When we used a very good approximation for the solar differential rotation at the surface (N = 30 terms in the expansion to describe Ω), we found 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 tracked their frequencies in the complex plane by slowly transitioning from the case Ω(N = 2) to the solar case Ω(N = 30), see, for example, the online movie for m = 2.

The real parts of the eigenfrequencies of the modes R, R₁, and R₂ do not change very much compared to the case of the two-term rotation profile. However, the imaginary parts of the eigenfrequencies do change. In particular, those of the (symmetric) R₂ mode with m = 1 and the (antisymmetric) R₁ mode with m = 2 are now positive, that is 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 they are thus more likely to be observable in the solar data.

Figure C.2 shows the eigenfunctions for selected modes with m = 2 shown in Fig. 2. The three Rossby modes have eigenfunctions that are close to , that is close to the case of uniform rotation. For m = 2, the critical latitudes are very close to the poles and have little effect on the eigenfunctions. For larger values of m, the eigenfunctions are confined between the critical layers and have a significant imaginary part (see Gizon et al. 2020, for m = 10).

The eigenfunctions of the high-latitude modes A₁ and A₂ have a modulus that peaks near the critical layers at ≈ ± 70° and the argument between the real and imaginary parts varies below this latitude (which implies a spiral pattern there). The eigenfunctions sensitively depend 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). In Sect. 5.2, we show that an improved match can be obtained when the solar rotation profile is taken deep in the convection zone.

The critical-latitude modes P₁ and P₂ have eigenfunctions that oscillate around and below the critical latitudes. Their real and imaginary parts are not in phase as functions of latitude and they 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).

4.1. 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 it increases drastically when the Ekman number becomes larger than 10⁻³. More surprisingly, the real part is also significantly affected for modes A₁ and P₁, by an amount of more than 100 nHz depending on the value for the viscosity. The R mode is hardly affected for this small value of m. For larger values of m, we observed a change in the R-mode frequency of, for example, 15 nHz for m = 6 and 30 nHz for m = 10 over the range of Ekman numbers covered by the plot. All of these variations would be measurable in the solar observations and thus the modes may be used as probes for the viscosity in the solar interior. The eigenfunctions are also affected (see Fig. 3) and, in particular, the imaginary part of mode A₁ changes sign with viscosity and mode P₁ becomes more and more confined between the viscous layers as viscosity decreases. The shape of the R-mode eigenfunction does not significantly depend on the viscosity for this small value of m, unlike for larger values of m as already discussed by Gizon et al. (2020).

Fig. 3. Influence of the viscosity on the eigenvalues and eigenfunctions of the R, A1, and P1 modes. (a): real (solid curves) and imaginary (dashed curves) parts of the frequencies of the modes R, A1, and P1 for m = 2 as a function of the Ekman number in the case of the observed solar surface rotation profile. Panels b–g: real (top panels) and imaginary (bottom panels) parts of the stream functions for the same modes at three selected values for the Ekman number. See also the top panels of Fig. C.3 for a 2D representation of the stream functions for these modes.

4.2. Effect of meridional flow

In order to study the importance of the meridional flow V(θ), we considered the background flow

(23)

At the surface, we used a poleward flow with maximum amplitude V_max = 15 m/s,

(24)

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₁ 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 in the imaginary part of the frequency of the R₁ 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 midlatitudes where the meridional flow is largest. The real parts of the eigenfrequencies of the m = 2 modes located on the P branch get closer to zero when the meridional flow is included with a shift ∼10 nHz. The meridional flow also stretches the eigenfunctions in latitude, see, for example, the right panel of Fig. C.5 for the P₁ mode for m = 2.

4.3. Effect of time-varying zonal flows

To assess the sensitivity of the modes to the details of the rotation profile, we computed the mode frequencies as a function of time over the last two solar cycles. We used the inferred rotation profile from helioseseismology (Larson & Schou 2018), obtained at a cadence of 72 days. The rotation profiles were averaged over five bins to obtain yearly averages from 1996–2018. The time variations of the frequencies of the least-damped modes with m = 2 are plotted in Fig. C.6. We found that the frequency of the R mode varies by less than 1 nHz, which is in agreement with the calculation of Goddard et al. (2020) using first-order perturbation theory. The critical-latitude mode P₁ also changes by a very small amount of ±2 nHz. The frequency of the high-latitude mode A₁ varies by up to ±7 nHz and shows a 22-year periodicity. Other modes have frequencies that strongly vary (such as P₅, not shown in 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.

5. 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.

5.1. Stability analysis

The hydrodynamical stability of solar latitudinal differential rotation was studied in two dimensions for purely toroidal disturbances 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 (Rayleigh 1879). Fjørtoft (1950) proved a more restrictive condition for instability: there exist a θ and a θ₀ (where θ₀ is a zero of ζ), such that [Ω(cos θ)−Ω(cos θ₀)]ζ(θ) < 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 Ω = Ω₀ + Ω₂cos²θ, Watson (1981) showed that a necessary condition for stability is −(2/7)Ω₀ < Ω₂ < 1.14Ω₀. This condition is met for two-term fits to the solar rotation profile. Charbonneau et al. (1999) added a third term in Ω₄cos⁴θ and found numerically that when Ω₄ ≈ Ω₂ ≈ −0.1Ω₀ (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.74 R_⊙, 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⁻⁴ (five modes with m = 1, five modes with m = 2, and two modes with m = 3). We find that for E < 10⁻⁴, the radius above which some modes become unstable (≈0.75 R_⊙ for m = 1, ≈0.77 R_⊙ for m = 2, and ≈0.78 R_⊙ for m = 3) does not sensitively depend on the viscosity and it is almost given by the E = 0 limit. Only values of E > 10⁻⁴ have a stabilising effect. For E > 5 × 10⁻³, all depths are stable. For 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 modes R₁ and R₂ for m = 1 and m = 2, and R and R₁ 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 the shape of the eigenfunctions much which resemble the classical 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₁ and R₂ for m = 1 and R₁ for m = 2) and no modes are unstable below 0.91 R_⊙ instead of 0.75 R_⊙ using E_η. The amplitudes of the inertial modes may thus have a diagnostic value to learn about eddy viscosity in the solar interior.

Fig. 4. Number of unstable modes as a function of the Ekman number at different radii for m = 1, 2, 3. For m ≥ 4, all modes are stable. The solar differential rotation is from Larson & Schou (2018). The dashed-black line gives the Ekman number EMLT as a function of the radius using the turbulent viscosity from mixing-length theory (Muñoz-Jaramillo et al. 2011). The solid black line shows the Ekman number Eη, using the quenched diffusivity model proposed by Muñoz-Jaramillo et al. (2011), see their Fig. 4b. The regions corresponding to the most unstable modes (R2, R1, and R) and to the two most unstable modes (R2 and R1, R1 and R2, and R and R1) are colored and labeled with these modes. The color bar indicates the total number of unstable modes at any point in the diagram.

5.2. Propagation diagram and comparison with observations

Figure 5 shows the propagation diagram for different modes when the radius is either r = R_⊙ or r = 0.75 R_⊙. 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 drew 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 they 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₁ and R₂ which are self-excited. 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₁, A₂, R₁, and R₂. 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 a 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). Figure C.8 shows a different representation of the propagation diagram where the frequency is divided by the wavenumber m. Modes with similar values of Re[ω]/m have similar critical latitudes for all values of m. This diagram highlights the separation in latitude between the different families of modes.

Fig. 5. Model dispersion relations for all the modes with −Im[ω]/2π < 100 nHz. Shown are the Rossby modes (black curves for R, R1, and R2), 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.75 R⊙ 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).

Figure C.3 is a comparison of the eigenfunctions of the modes R, P₁, and A₁ calculated at the surface and at the bottom of the convection zone for m = 2. 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₁ and – to a lesser extent – the critical-latitude mode P₁ have eigenfunctions that significantly change 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).

6. 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. Similar to the equatorial β plane, viscous critical layers appear and the spectrum contains different families of modes due to 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 for this problem, as in the β plane problem. Using the surface solar rotation profile, two additional Rossby modes, R₁ (for m ≤ 6) and R₂ (for m ≤ 2), are also present. Remarkably, the calculated spectrum closely resembles 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 they 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 Ω_2pcos^2pθ with p ≥ 2 are essential to destabilize the system. For the Sun, an oversimplified fit of the form Ω₀ + Ω₂cos²θ would lead to the wrong conclusion that the system is stable. Using a value for 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₁ and R₂ for m = 1 and m = 2, and R and R₁ for m = 3. If the modes R₁ and R₂ 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.78 R_⊙, 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 for only a few Kepler Sun-like stars (Benomar et al. 2018). For these stars, the Ω₂cos²θ term is very large compared to the equatorial value Ω₀, 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.

Acknowledgments

This work is supported by the ERC Synergy Grant WHOLE SUN #810218 and by the DFG Collaborative Research Center SFB 1456 (project C04). LG acknowledges NYUAD Institute Grant G1502. LH acknowledges an internship agreement between SUPAERO and the MPS as part of her Bachelor thesis. The source code is available at https://doi.org/10.17617/3.OM51HE.

References

Appendix A: Eigenvalue problem for m = 1

In order to ensure that ψ satisfies the boundary conditions (Eq. (15)) for m = 1, we expanded the solution on the basis of associated Legendre polynomials instead of :

(A.1)

Next, we needed to specify how the operators L₁ and (L₁)² act on :

(A.2)

(A.3)

The eigenvalue problem, as seen in Eq. (12), becomes

(A.4)

where I is the (L − 1)×(L − 1) identity matrix, C is defined by Eq. (20) with m = 1, and the elements of the matrices A, D, and G are given by the following:

(A.5)

(A.6)

(A.7)

with

(A.8)

Appendix B: Inclusion of meridional flow

Here we derive the equation for the stream function when the background flow includes both the differential rotation and a meridional flow V(θ), that is

(B.1)

The horizontal components of Eq. (3) become

(B.2)

(B.3)

Combining Eq. (B.2) and Eq. (B.3) and using the definition of the stream function (Eq. (7)), we obtain

(B.4)

where ζ(θ) is defined by Eq. (9). For each longitudinal mode, we have

(B.5)

where .

B.1. Case m ≥ 2

In this case, the function ψ(θ) is expanded according to Eq. (17). Using , we obtain

(B.6)

This leads to the linear system

(B.7)

where matrix C is given by Eq. (20) and matrix C^mer has elements

(B.8)

For the meridional flow defined by Eq. (24), the effects of matrix C^mer on the spectrum can be seen in Fig. C.4 for the case where m = 2. Example eigenfunctions are plotted in Fig. C.5.

B.2. Case m = 1

When m = 1, we used a decomposition of ψ on the basis of in order to enforce the boundary conditions (see Sect. A). The eigenvalue problem has thus become

(B.9)

where D and G are defined by Eqs. (A.6) and (A.7) and C^mer and D^mer are given by

(B.10)

(B.11)

Appendix C: Supplementary figures

Fig. C.1. Eigenfrequencies ω for m = 1 to 3 (from left to right) when the Ekman number is E = 2 × 10−5. The labels and the colors are the same as in Fig. 2. The top panels are for Ω(N = 2) = Ω0 + Ω2cos2θ, and the bottom panels are for the solar rotation rate Ω(N = 30) at a radius r = 0.75 R⊙.

Fig. C.2. Stream functions ψ of selected 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, and the dashed curves illustrate 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.

Fig. C.3. Stream functions Ψ(θ, ϕ, t = 0) for the m = 2 modes R, A1, and P1, using the surface solar differential rotation and the Ekman number E = 2 × 10−4 (top) and at the bottom of the convection zone (r = 0.75 R⊙ with E = 2 × 10−5). Latitudes and longitudes are highlighted every 30° with dotted curves and the equator is shown with a solid line. The red curves show the latitudes of the viscous critical layers. The eigenfrequencies are expressed in the Carrington frame of reference.

Fig. C.4. Effect of the meridional flow (Eq. 24) on the m = 2 mode frequencies, ω. The surface latitudinal differential rotation is from Larson & Schou (2018). The black squares show the mode frequencies when the meridional flow is included and the red circles are for when it is not. The full symbols correspond to the symmetric modes, and the open symbols correspond to the antisymmetric modes. The Ekman number is fixed at E = 4 × 10−4 and all frequencies are expressed in the Carrington frame of reference. The eigenfunctions of the modes R, A1, and P1 are shown in Fig. C.5.

Fig. C.5. Eigenfunctions of the R, A1, and P1 modes with (red) and without (blue) the meridional flow. The real parts are given by the solid curves and the imaginary parts are illustrated by the dashed curves. The computations are for m = 2 and E = 4 × 10−4.

Fig. C.6. Temporal changes in the frequencies ω of selected m = 2 modes caused by the Sun’s zonal flows (at the surface). The left panel shows the real parts of the frequencies and the right panel is for the imaginary parts. The Ekman number is E = 4 × 10−4. It is important to notice that mode R1 is unstable part of the time.

Fig. C.7. Stream functions Ψ(θ, ϕ, t = 0) for the m = 1, 2, and 3 Rossby modes R, R1, and R2, some of which may be unstable. The differential rotation is that of the Sun’s surface (temporal average) and the Ekman number is E = 2 × 10−4. The red curves show the latitudes of the viscous critical layers.

Fig. C.8. Dispersion relation diagrams for all the modes with |Im[ω]/2π|< 100 nHz. This figure is similar to Fig. 5, but the ordinate has been replaced by ω/m so that the critical latitudes occur on horizontal lines. A few selected critical (co)latitudes θc are highlighted (horizontal gray lines)

Movies

Movie 1 associated with Fig. 2 (movie_eigenvalues_m1) (Access here)

Movie 2 associated with Fig. 2 (movie_eigenvalues_m2) (Access here)

All Tables

All Figures

Fig. 1. Main quantities entering in the eigenvalue problem defined by Eq. (12). (a) Normalized solar differential rotation rate δ(θ) at different radii, inferred from helioseismology and averaged over 1996 to 2018 (the shaded area covers ±3σ), with respect to the Carrington rotation rate, Ωref = ΩCarr. The solid curves show the N = 30 fits to the solar rotation rates. The dashed curve shows the N = 2 fit to the solar surface rotation rate, i.e., Ω0 + Ω2cos2θ with Ω0/2π = 452 nHz and Ω2/2π = −114 nHz. (b) The function ζ defined by Eq. (9) is plotted for each of the three cases of panel (a). For r = R⊙ and N = 30, the function reaches a minimum at ζ(θ = 0)≈ − 22 (not shown on the plot). The oscillations at the surface are due to the m = 0 zonal flows (also known as torsional oscillations). The horizontal dashed-black line shows , where β = 2Ωref/r, which plays the role of ζ in the equatorial β plane (see Eq. (17) in Gizon et al. 2020).

In the text

Fig. 2. Eigenfrequencies ω for m = 1–3 (from left to right). The Ekman number is E = 4 × 10−4. The top panels are for Ω(N = 2) = Ω0 + Ω2cos2θ, with Ω0/2π = 452 nHz and Ω2/2π = −114 nHz, and the bottom panels are for the solar rotation rate Ω(N = 30). Modes with north-south symmetric eigenfunctions are shown with filled symbols, and 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, R1, and R2 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) was obtained by tracking the modes from the N = 2 case (top panels). See the evolution between the two profiles in the online movies for m = 1 and m = 2.

In the text

Fig. 3. Influence of the viscosity on the eigenvalues and eigenfunctions of the R, A1, and P1 modes. (a): real (solid curves) and imaginary (dashed curves) parts of the frequencies of the modes R, A1, and P1 for m = 2 as a function of the Ekman number in the case of the observed solar surface rotation profile. Panels b–g: real (top panels) and imaginary (bottom panels) parts of the stream functions for the same modes at three selected values for the Ekman number. See also the top panels of Fig. C.3 for a 2D representation of the stream functions for these modes.

In the text

Fig. 4. Number of unstable modes as a function of the Ekman number at different radii for m = 1, 2, 3. For m ≥ 4, all modes are stable. The solar differential rotation is from Larson & Schou (2018). The dashed-black line gives the Ekman number EMLT as a function of the radius using the turbulent viscosity from mixing-length theory (Muñoz-Jaramillo et al. 2011). The solid black line shows the Ekman number Eη, using the quenched diffusivity model proposed by Muñoz-Jaramillo et al. (2011), see their Fig. 4b. The regions corresponding to the most unstable modes (R2, R1, and R) and to the two most unstable modes (R2 and R1, R1 and R2, and R and R1) are colored and labeled with these modes. The color bar indicates the total number of unstable modes at any point in the diagram.

In the text

Fig. 5. Model dispersion relations for all the modes with −Im[ω]/2π < 100 nHz. Shown are the Rossby modes (black curves for R, R1, and R2), 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.75 R⊙ 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).

In the text

	Fig. C.1. Eigenfrequencies ω for m = 1 to 3 (from left to right) when the Ekman number is E = 2 × 10−5. The labels and the colors are the same as in Fig. 2. The top panels are for Ω(N = 2) = Ω0 + Ω2cos2θ, and the bottom panels are for the solar rotation rate Ω(N = 30) at a radius r = 0.75 R⊙.
In the text

Fig. C.2. Stream functions ψ of selected 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, and the dashed curves illustrate 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.

In the text

Fig. C.3. Stream functions Ψ(θ, ϕ, t = 0) for the m = 2 modes R, A1, and P1, using the surface solar differential rotation and the Ekman number E = 2 × 10−4 (top) and at the bottom of the convection zone (r = 0.75 R⊙ with E = 2 × 10−5). Latitudes and longitudes are highlighted every 30° with dotted curves and the equator is shown with a solid line. The red curves show the latitudes of the viscous critical layers. The eigenfrequencies are expressed in the Carrington frame of reference.

In the text

Fig. C.4. Effect of the meridional flow (Eq. 24) on the m = 2 mode frequencies, ω. The surface latitudinal differential rotation is from Larson & Schou (2018). The black squares show the mode frequencies when the meridional flow is included and the red circles are for when it is not. The full symbols correspond to the symmetric modes, and the open symbols correspond to the antisymmetric modes. The Ekman number is fixed at E = 4 × 10−4 and all frequencies are expressed in the Carrington frame of reference. The eigenfunctions of the modes R, A1, and P1 are shown in Fig. C.5.

In the text

	Fig. C.5. Eigenfunctions of the R, A1, and P1 modes with (red) and without (blue) the meridional flow. The real parts are given by the solid curves and the imaginary parts are illustrated by the dashed curves. The computations are for m = 2 and E = 4 × 10−4.
In the text

	Fig. C.6. Temporal changes in the frequencies ω of selected m = 2 modes caused by the Sun’s zonal flows (at the surface). The left panel shows the real parts of the frequencies and the right panel is for the imaginary parts. The Ekman number is E = 4 × 10−4. It is important to notice that mode R1 is unstable part of the time.
In the text

	Fig. C.7. Stream functions Ψ(θ, ϕ, t = 0) for the m = 1, 2, and 3 Rossby modes R, R1, and R2, some of which may be unstable. The differential rotation is that of the Sun’s surface (temporal average) and the Ekman number is E = 2 × 10−4. The red curves show the latitudes of the viscous critical layers.
In the text

	Fig. C.8. Dispersion relation diagrams for all the modes with |Im[ω]/2π|< 100 nHz. This figure is similar to Fig. 5, but the ordinate has been replaced by ω/m so that the critical latitudes occur on horizontal lines. A few selected critical (co)latitudes θc are highlighted (horizontal gray lines)
In the text