Issue 
A&A
Volume 664, August 2022



Article Number  A6  
Number of page(s)  16  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202243473  
Published online  04 August 2022 
Viscous inertial modes on a differentially rotating sphere: Comparison with solar observations^{⋆}
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: fournier@mps.mpg.de
^{2}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3}
Center for Space Science, NYUAD Institute, New York University Abu Dhabi, Abu Dhabi, UAE
^{4}
Institut Supérieur de l’Aéronautique et de l’Espace – SUPAERO, 10 avenue Édouard Belin, 31055 Toulouse, France
Received:
3
March
2022
Accepted:
18
April
2022
Context. In a previous paper, we studied the effect of latitudinal rotation on solar equatorial Rossby modes in the β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 highlatitude modes.
Aims. 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.
Methods. At a fixed radius, we solved the eigenvalue problem numerically using a spherical harmonics decomposition of the velocity stream function.
Results. Due to the presence of viscous critical layers, the spectrum consists of four different families: Rossby modes, highlatitude modes, criticallatitude modes, and strongly damped modes. For each longitudinal wavenumber m ≤ 3, up to three Rossbylike modes are present on the sphere, in contrast to the equatorial β 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.75 R_{⊙} and Ekman numbers E < 10^{−4}, at least one mode is unstable. For either m = 1 or m = 2, up to two Rossby modes (one symmetric and one antisymmetric) 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.
Conclusions. 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 selfexcited modes in the model have frequencies close to those of the observed modes with the largest amplitudes.
Key words: waves / hydrodynamics / instabilities / Sun: interior / Sun: rotation / methods: numerical
Movies associated to Fig. 2 are available at https://www.aanda.org
© D. Fournier et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model.
Open Access funding provided by Max Planck Society.
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), highlatitude modes (the m = 1 symmetric mode manifests itself as a spiral structure, see Hathaway et al. 2013; Bogart et al. 2015), and criticallatitude 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 criticallatitude modes, the highlatitude 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 highlatitude 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 shallowwater 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 fourterm 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 selfexcited (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 NavierStokes 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 was modeled by a horizontal Laplacian Δ. The force on the righthand 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′,
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
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 Eqs. (4) and (5), we obtained the equation of Watson (1981) which was modified on the righthand side to include viscosity:
where
2.2. Modal decomposition
We looked 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
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 fourthorder, which has profound implications for the spectrum. Four boundary conditions are required. We imposed the flow to vanish at the poles:
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
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,
where the are normalized, such that . We inserted this expansion into Eq. (12) and projected onto a particular polynomial . Using
we obtained the following matrix equation:
where b = [b_{m}b_{m + 1}… b_{L}]^{T} is a vector and matrix C has elements
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 northsouth symmetric stream function are referred to as symmetric, that is to say they are symmetric in and antisymmetric in . The modes with a northsouth antisymmetric stream function are called antisymmetric. For the sake of simplicity, we do not consider the case of a general rotation profile with northsouth 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^{−4}, which means ν = 125 km^{2} s^{−1} 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:
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 northsouth symmetric rotation, the odd coefficients are zero. Rotation from globalmode helioseismology is northsouth 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 solarcyclerelated zonal flows. Adding more terms would not change the surface profile significantly. Below the nearsurface 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} + Ω_{2}cos^{2}θ 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 dashedblack 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, Ω = Ω_{0}, the complex eigenvalues can be obtained directly from Eq. (19) (matrix C is diagonal because δ and ζ are constant):
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. Twoterm differential rotation and E ≳ 10^{−4}
Here we solve the eigenvalue problem given by Eq. (19) using Ω_{0} + Ω_{2}cos^{2}θ as an approximation for the solar differential rotation at the surface. Fitting the observed surface rotation profile, we used Ω_{0}/2π = 452 nHz and Ω_{2}/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 solarlike. We used an Ekman number E = 4 × 10^{−4} 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 Yshape in the complex plane. This spectrum is characteristic of the spectrum of the OrrSommerfeld 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} + Ω_{2}cos^{2}θ, with Ω_{0}/2π = 452 nHz and Ω_{2}/2π = −114 nHz, and the bottom panels are for the solar rotation rate Ω(N = 30). Modes with northsouth 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 OrrSommerfeld 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) 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. 
Correspondence of terminologies for viscous modes in shear flows.
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_{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, that is to say Ω(θ) = Ω_{0} + ϵΩ_{2}cos^{2}θ 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.
Description of the viscous modes discussed in this paper for the case where Ω = Ω_{0} + Ω_{2}cos^{2}θ.
3.3. Twoterm 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., Fig. 1 of MuñozJaramillo 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^{−5} (instead of 4 × 10^{−4} as in Fig. 2). For m = 2 and m = 3, the spectrum is not Yshaped 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 Yshape 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_{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, that is these modes are unstable (selfexcited). The second derivative of the rotation profile has a strong influence on the damping of the modes, especially the highlatitude and criticallatitude 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 solarlike Ω(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 highlatitude 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 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 criticallatitude 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 they have roughly the same amplitude. These modes are not very different from the case of the twoterm rotation profile. They resemble the observed m = 2 criticallatitude 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^{−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 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 Rmode 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_{1} changes sign with viscosity and mode P_{1} becomes more and more confined between the viscous layers as viscosity decreases. The shape of the Rmode 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, A_{1}, and P_{1} modes. (a): real (solid curves) and imaginary (dashed curves) parts of the frequencies of the modes R, A_{1}, and P_{1} 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
At the surface, we used a poleward flow with maximum amplitude V_{max} = 15 m/s,
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 in 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 criticallatitude 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_{1} mode for m = 2.
4.3. Effect of timevarying 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 leastdamped 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 firstorder perturbation theory. The criticallatitude mode P_{1} also changes by a very small amount of ±2 nHz. The frequency of the highlatitude mode A_{1} varies by up to ±7 nHz and shows a 22year periodicity. Other modes have frequencies that strongly vary (such as P_{5}, not shown in the plot, which varies by up to ±10 nHz). Since the frequency resolution corresponding to a 3year 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 θ_{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 twoterm 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 pmode 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^{−4} (five modes with m = 1, five modes with m = 2, and two modes with m = 3). We find that for E < 10^{−4}, 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^{−4} have a stabilising effect. For E > 5 × 10^{−3}, all depths are stable. For the same figure, we draw estimates of the Ekman number with depth, either under mixinglength theory or using a quenched diffusivity (MuñozJaramillo 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_{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 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_{1} and R_{2} for m = 1 and R_{1} 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 dashedblack line gives the Ekman number E_{MLT} as a function of the radius using the turbulent viscosity from mixinglength theory (MuñozJaramillo et al. 2011). The solid black line shows the Ekman number E_{η}, using the quenched diffusivity model proposed by MuñozJaramillo et al. (2011), see their Fig. 4b. The regions corresponding to the most unstable modes (R_{2}, R_{1}, and R) and to the two most unstable modes (R_{2} and R_{1}, R_{1} and R_{2}, and R and R_{1}) 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 highlatitude modes have frequencies that overlap with the region between the slowest and fastest highlatitude 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 highlatitude modes for m = 1, 2, and 3 have frequencies close to the modes R_{1} and R_{2} which are selfexcited. This may explain why the observed highlatitude 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 criticallatitude modes with frequencies between −150 nHz and 0 nHz may be associated with the dense spectrum of criticallatitude 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 criticallatitude 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, R_{1}, and R_{2}), the highlatitude modes (red area), the strongly damped modes (blue area), and the criticallatitude 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_{1}, and A_{1} calculated at the surface and at the bottom of the convection zone for m = 2. For this low m value, the Rmode 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 quasitoroidal 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 highlatitude mode A_{1} and – to a lesser extent – the criticallatitude mode P_{1} 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 highlatitude 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 highlatitude 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 highlatitude, criticallatitude, 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_{1} (for m ≤ 6) and R_{2} (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 selfexcited, 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 for the viscosity corresponding to the quenched diffusivity model from MuñozJaramillo et al. (2011), we find that six modes are selfexcited. 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 highlatitude 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 Sunlike 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 highorder terms that fully determine the rotation profile of these stars.
Rotational stability of selected Kepler stars for which Ω_{0} and Ω_{2} were measured using asteroseismology (Benomar et al. 2018).
Movies
Movie 1 associated with Fig. 2 (movie_eigenvalues_m1) Access here
Movie 2 associated with Fig. 2 (movie_eigenvalues_m2) Access here
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
 Bekki, Y., Cameron, R. H., & Gizon, L. 2022, A&A, 662, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benomar, O., Bazot, M., Nielsen, M. B., et al. 2018, Science, 361, 1231 [Google Scholar]
 Bogart, R. S., Baldner, C. S., & Basu, S. 2015, ApJ, 807, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P., Dikpati, M., & Gilman, P. A. 1999, ApJ, 526, 523 [NASA ADS] [CrossRef] [Google Scholar]
 Dikpati, M., & Gilman, P. A. 2001, ApJ, 551, 536 [NASA ADS] [CrossRef] [Google Scholar]
 Drazin, P. G., & Reid, W. H. 1981, Hydrodynamic Stability (Cambridge University Press) [Google Scholar]
 Dziembowski, W., & Kosovichev, A. 1987, Acta Astron., 37, 341 [Google Scholar]
 Fjørtoft, R. 1950, Geofys. Publ., 17, 1 [Google Scholar]
 Gilman, P. A., & Dikpati, M. 2000, ApJ, 528, 552 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A., & Dikpati, M. 2002, ApJ, 576, 1031 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A., Dikpati, M., & Miesch, M. S. 2007, ApJS, 170, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Fournier, D., & Albekioni, M. 2020, A&A, 642, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goddard, C. R., Birch, A. C., Fournier, D., & Gizon, L. 2020, A&A, 640, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grosch, C. E., & Salwen, H. 1968, J. Fluid Mech., 34, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 2009, A&A, 504, 303 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29 [Google Scholar]
 Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
 Mack, L. M. 1976, J. Fluid Mech., 73, 497 [NASA ADS] [CrossRef] [Google Scholar]
 Maslowe, S. 2003, Ann. Rev. Fluid Mech., 1, 1 [Google Scholar]
 MuñozJaramillo, A., Nandy, D., & Martens, P. C. H. 2011, ApJ, 727, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Pekeris, C. L. 1948, Phys. Rev., 74, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Rayleigh, L. 1879, Proc. London Math. Soc., s111, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., Valdettaro, L., & Georgeot, B. 2002, J. Fluid Mech., 463, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Rossby, C. G. 1939, J. Mar. Res., 2, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Saio, H. 1982, ApJ, 256, 717 [Google Scholar]
 Schensted, I. V. 1961, PhD Thesis, The University of Michigan, Ann Arbor [Google Scholar]
 Schmid, P. J., & Henningson, D. S. 2012, Stability and Transition in Shear Flows (Springer), 142 [Google Scholar]
 Watson, M. 1981, Geophys. Astrophys. Fluid Dyn., 16, 285 [Google Scholar]
 Zaqarashvili, T. V., Albekioni, M., Ballester, J. L., et al. 2021, Space Sci. Rev., 217, 15 [Google Scholar]
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 :
Next, we needed to specify how the operators L_{1} and (L_{1})^{2} act on :
The eigenvalue problem, as seen in Eq. (12), becomes
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:
with
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
The horizontal components of Eq. (3) become
Combining Eq. (B.2) and Eq. (B.3) and using the definition of the stream function (Eq. (7)), we obtain
where ζ(θ) is defined by Eq. (9). For each longitudinal mode, we have
where .
B.1. Case m ≥ 2
In this case, the function ψ(θ) is expanded according to Eq. (17). Using , we obtain
This leads to the linear system
where matrix C is given by Eq. (20) and matrix C^{mer} has elements
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
where D and G are defined by Eqs. (A.6) and (A.7) and C^{mer} and D^{mer} are given by
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} + Ω_{2}cos^{2}θ, 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, A_{1}, and P_{1}, 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, A_{1}, and P_{1} are shown in Fig. C.5. 
Fig. C.5. Eigenfunctions of the R, A_{1}, and P_{1} 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 R_{1} is unstable part of the time. 
Fig. C.7. Stream functions Ψ(θ, ϕ, t = 0) for the m = 1, 2, and 3 Rossby modes R, R_{1}, and R_{2}, 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) 
All Tables
Description of the viscous modes discussed in this paper for the case where Ω = Ω_{0} + Ω_{2}cos^{2}θ.
Rotational stability of selected Kepler stars for which Ω_{0} and Ω_{2} were measured using asteroseismology (Benomar et al. 2018).
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} + Ω_{2}cos^{2}θ 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 dashedblack 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} + Ω_{2}cos^{2}θ, with Ω_{0}/2π = 452 nHz and Ω_{2}/2π = −114 nHz, and the bottom panels are for the solar rotation rate Ω(N = 30). Modes with northsouth 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 OrrSommerfeld 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) 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, A_{1}, and P_{1} modes. (a): real (solid curves) and imaginary (dashed curves) parts of the frequencies of the modes R, A_{1}, and P_{1} 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 dashedblack line gives the Ekman number E_{MLT} as a function of the radius using the turbulent viscosity from mixinglength theory (MuñozJaramillo et al. 2011). The solid black line shows the Ekman number E_{η}, using the quenched diffusivity model proposed by MuñozJaramillo et al. (2011), see their Fig. 4b. The regions corresponding to the most unstable modes (R_{2}, R_{1}, and R) and to the two most unstable modes (R_{2} and R_{1}, R_{1} and R_{2}, and R and R_{1}) 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, R_{1}, and R_{2}), the highlatitude modes (red area), the strongly damped modes (blue area), and the criticallatitude 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} + Ω_{2}cos^{2}θ, 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, A_{1}, and P_{1}, 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, A_{1}, and P_{1} are shown in Fig. C.5. 

In the text 
Fig. C.5. Eigenfunctions of the R, A_{1}, and P_{1} 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 R_{1} 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, R_{1}, and R_{2}, 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.