© L. Gizon et al. 2020

1. Introduction

In the atmosphere of Earth, Rossby (1939) waves are global-scale waves of radial vorticity that propagate in the direction opposite to rotation (retrograde). They originate in the conservation of vertical absolute vorticity, that is, the sum of planetary and wave vorticity (see, e.g., Platzman 1968; Gill 1982).

Equatorial Rossby waves have recently been observed on the Sun with longitudinal wavenumbers in the range 3 ≤ m ≤ 15 (Löptien et al. 2018; Liang et al. 2019). In the corotating frame, their dispersion relation is close to that of classical sectoral (l = m) Rossby waves, ω = −2Ω/(m + 1), where Ω/2π = 453.1 nHz is the equatorial rotation rate.

The observed variation of the eigenfunctions with latitude, however, differs noticeably from , where λ is the latitude, which is the expected answer for sectoral modes in a uniformly rotating sphere (e.g., Saio 1982; Damiani et al. 2020). Instead, the observed eigenfunctions have real parts that peak at the equator, switch sign near 20°–30°, and decay at higher latitudes (Löptien et al. 2018). Their imaginary parts are small (Proxauf et al. 2020).

An ingredient that is obviously missing in models of solar Rossby waves is latitudinal differential rotation. The solar rotation rate decreases fast with latitude: the difference between the rotation rate at mid-latitudes and at the equator is not small compared to the frequencies of the Rossby waves of interest. For m larger than 5, we show that a critical latitude exists at which the (negative) wave frequency equals the (negative) differential rotation rate counted from the equator.

To capture the essential physics while keeping the problem simple, we chose to work in two dimensions and in the equatorial β plane. This simplification is acceptable for wavenumbers that are large enough (about ≥6).

The stability and dynamics of parabolic (Poiseuille) shear flows in the presence of critical layers was summarized by Drazin & Reid (2004), for example. Kuo (1949) included the β effect in the problem. In the inviscid case, critical layers lead to a singular eigenvalue problem (see, e.g., Balmforth & Morrison 1995). The stream functions are continuous but not differentiable (Sect. 4), thus they cannot be compared directly to actual observations of the vorticity. Because we also wish to explain the lifetime of the modes (Löptien et al. 2018), we introduce a viscous term in the Navier-Stokes equations to model damping by turbulent convection. As shown in Sect. 2, this leads to a new equation for the stream function: an Orr-Sommerfeld equation whose coefficients are modified by the β effect. The viscosity removes singularities, and the eigenfunctions are regular across the viscous critical layer (see Sect. 3). To solve the eigenvalue problem accurately, we use a numerical method based on the Chebyshev decomposition proposed by Orszag (1971). As shown in Sect. 5, the eigenvalue spectrum includes a symmetric Rossby mode in addition to the other modes that are known to exist in the β = 0 case (Mack 1976). We focus on the eigenvalue problem here. We do not discuss the nonlinear dynamics in the critical layers, which would require solving the nonlinear evolution equation (e.g., Stewartson 1977).

In addition to the practical advantages of studying 2D Rossby waves in the β plane, the physics of this problem has been extensively discussed in Earth and planetary sciences. In the terrestrial atmosphere and oceans, Rossby waves encounter critical layers (see Frederiksen & Webster 1988; Vallis 2006; Boyd 2018). They play a role in the global dynamics by transporting angular momentum through Reynolds stresses, and they modify the mean flow (e.g., Bennett & Young 1971; Webster 1973; Geisler & Dickinson 1974).

Our model can be further extended to include the effect of the solar meridional flow using the method of Orszag (1971), see Sect. 6. In Sect. 7 we compare the predictions of the model to solar observations of the mode frequencies and damping rates and to observations of the vorticity eigenfunctions. Finally, in Sect. 8 we discuss the implications of the model for angular momentum transport and the dynamics of solar differential rotation.

2. Waves in a sheared zonal flow: Equations of motion in the equatorial β plane

In the equatorial β plane, we study the propagation of 2D Rossby waves through a steady zonal flow representative of solar differential rotation. Several β-plane approximations have been proposed (see, e.g., Dellar 2011). Here we chose the sine transformation,

(1)

(2)

where ϕ is the longitude, λ is the latitude, and R = 696 Mm is the solar radius. The x coordinate increases in the prograde direction, and the y coordinate increases northward. To first order in y/R, the x and y components of the velocity vector in the β plane are equal to their ϕ and λ components on the sphere (Ripa 1997), respectively.

The total velocity is the sum of the background zonal flow and the horizontal wave velocity u(x, y, t) with

(3)

By choice, U(0) = 0 at the equator. The latitudinal shear is specified by the parabolic (Poiseuille) profile,

(4)

where the amplitude is chosen such that U(y) approximates the solar surface differential rotation at low and mid-latitudes. From the standard solar angular velocity profile given by Beck (2000), ΔΩ = −0.35[(y/R)² + (y/R)⁴] μrad s⁻¹, we find that the value m s⁻¹ is a good choice, see Fig. 1.

Fig. 1. Parabolic zonal flow U (red curve, Eq. (4)) in the frame rotating at the equatorial rotation rate, which approximates the solar rotational velocity at the photosphere (blue curve). The horizontal black lines indicate the phase speed of Rossby waves, c0 = −(β − U″)/k2, for longitudinal wavenumbers kR = 6 and kR = 10. At critical latitudes, U = c0.

The 2D Navier-Stokes equations in the equatorial β plane are

(5)

(6)

where ν is the viscosity, and the equatorial Coriolis parameter is f = βy with β = 2Ω/R. The prime denotes a derivative, for example, U′ = dU/dy.

To enforce mass conservation, we introduce the stream function Ψ such that

(7)

Assuming a barotropic fluid, we can eliminate the pressure term by combining the two components of the momentum equation to obtain

(8)

3. Critical latitudes

3.1. Linear inviscid case: critical points

In the linear inviscid case, we search for wave solutions of the form

(9)

When the nonlinear and viscous terms are dropped, Eq. (8) becomes the Rayleigh-Kuo equation (Kuo 1949),

(10)

where c = ω/k is the phase speed. This equation differs from the Rayleigh (1879) equation only through the β term. It can be rewritten as a Helmholtz equation,

(11)

with

(12)

The squared latitudinal wavenumber, K(y), is singular at the critical points y = ±y_c such that U(y_c) = c. The critical points divide the low-latitude region where the solution is locally oscillatory (K >  0 for |y|< y_c) from the high-latitude regions, where it is locally evanescent (K <  0 for |y|> y_c).

Equation (10), supplemented by boundary conditions ψ(±R) = 0, is an eigenvalue problem that can be solved in the complex plane. It admits a continuum of neutral modes with real eigenfrequencies. According to Rayleigh’s theorem (adapted for the Rayleigh-Kuo equation, see Kuo 1949), there is no discrete mode because β − U″≠0 everywhere. We can thus solve the initial value problem given by Eq. (10) for any particular real value of c to obtain the associated eigenfunction (e.g., Drazin & Howard 1966; Drazin et al. 1982; Balmforth & Morrison 1995). These eigenfunctions are singular at the critical points.

Because U″ is constant in our problem and U(0) = 0, Eq. (12) implies that each mode can be associated with a value of K(0) such that

(13)

For equatorial propagation, K(0) = 0 is a natural choice, and we may consider the eigenvalue

(14)

as an example. In our case, β − U″ = 1.12β, therefore waves propagate faster than in the no-flow case. The critical points y = ±y_c where U(±y) = c₀ are given by

(15)

Thus, for kR >  4.31, there are critical latitudes at λ = ±λ_c, such that

(16)

To obtain the eigenfunctions, Eq. (10) should be solved separately for |y|< y_c and |y|> y_c. The analytical and numerical solutions are discussed in Sect. 4. The stream function is continuous (u_y = −ikψ is continuous), but its first and second derivatives diverge at the critical layer (see, e.g., Haynes et al. 2003).

3.2. Viscous critical layer

Bulk viscosity removes singularities. The linear viscous equation for ψ is

(17)

For β = 0, this is the Orr-Sommerfeld equation (Orr 1907; Sommerfeld 1909). Equation (17) is a fourth-order differential equation, which requires four boundary conditions, such as ψ(±R) = 0 and ψ′(± R) = 0 for a no-slip boundary condition. The critical layer of the inviscid case is replaced by a viscous critical layer around y = ±y_c. The width of this viscous layer, δ, is obtained by balancing the dominant viscous term with the dominant term on the left-hand side, (U − c)ψ″ ∼ νψ″″/k. Close to the viscous layer, we write d/dy ∼ 1/δ and , so that the width of the viscous layer is approximately

(18)

where

(19)

is the Reynolds number.

The width of the layer is controlled by Re. As discussed by Rüdiger (1989), ν should be understood as an eddy viscosity due to turbulent convection. As an example, we may estimate the Reynolds number associated with solar supergranulation. For supergranulation, the turbulent viscosity ν ≈ 250 km² s⁻¹ (Simon & Weiss 1997; Duvall & Gizon 2000) implies Re ≈ 700 and δ ≈ 0.07R for kR = 10. Not surprisingly, the width of the viscous layer is comparable in this case to the spatial scale of supergranulation.

3.3. Nonlinear critical layer

In order to assess whether it is legitimate to drop the nonlinear term (u⋅∇)ΔΨ in Eq. (8), we estimate the width of the nonlinear critical layer δ_NL. It is obtained by balancing the advection term k(U − c)ψ″ and the nonlinear term u_yψ‴. We find

(20)

where u_max is a characteristic velocity amplitude for the Rossby waves. On the Sun, Liang et al. (2019) measured u_max ≈ 2 m s⁻¹ for the maximum latitudinal velocity of a mode at the equator. For kR = 10, the width of the nonlinear critical layer is δ_NL ≈ 0.04R.

Introducing the threshold , the ratio between the widths of the viscous and nonlinear critical layers is

(21)

For Re <  Re_* the critical layer is linear and dominated by dissipation over the width δ. For kR = 10, we find Re_* ≈ 3000, which is much higher than the Reynolds number Re ≈ 700 estimated in the previous section for solar supergranulation. It is therefore reasonable to study the linear problem, as we do in the remainder of this paper. However, we caution that there is some uncertainty about the appropriate value for the viscosity.

4. Inviscid modes of oscillation

In the inviscid case, the spectrum for the Rayleigh-Kuo equation, Eq. (10), is real and continuous. We fix the value of c to c₀ (Eq. (14)) and compute the corresponding eigenfunctions. There are two singular solutions, both real: a solution that is symmmetric in latitude, and an antisymmetric solution. These solutions can be obtained by solving the equation analytically or numerically in two distinct intervals, 0 ≤ y <  y_c and y_c ≤ y ≤ R.

In the inner region, we solve the ODE with the conditions ψ(0) = 1 and ψ′(0) = 0 to obtain the symmetric solution. For the region y >  y_c, we impose continuity with the inner solution and use the boundary condition ψ(R) = 0. The symmetric solution can be expressed as a series in each regions. We write ψ(y) = ∑_p ≥ 0a_pξ^2p with ξ = y/R for |y|< y_c and ψ(y) = ∑_p ≥ 1b_p(ξ² − 1)^p for |y|> y_c. The coefficients a_p and b_p are computed by recurrence. When we set κ = kR and ξ_c = y_c/R, the inner solution is given by

(22)

(23)

(24)

For the outer solution,

(25)

(26)

(27)

The coefficient b₁ is chosen such that the solution is continuous at the critical point.

To validate the series solution, we also solved the problem numerically. In the inner region, we have an initial value problem that can be solved using classical libraries, for example, +odeint+ from +SciPy+. In the outer region, the problem is a boundary value problem that can be converted into an initial value problem using the shooting method (Keller 1968).

We find an excellent agreement between the analytical and the numerical solutions, therefore we only plot the analytical solution in Fig. 2 for kR = 10. The symmetric eigenfunction ψ_s switches sign before the critical latitude and evanesces above it. The location of the critical layer is close to the zero-crossing of the observed vorticity eigenfunctions from Löptien et al. (2018) and Proxauf et al. (2020). Unfortunately, in the inviscid case, the vorticity ζ = k²ψ − ψ″ diverges at the critical latitude and thus the comparison with the observations is difficult.

Fig. 2. Symmetric inviscid eigenfunctions for the eigenvalue c0 = −(β − U″)/k2 and kR = 6, 8, 10, 12 and 14. The stars mark the critical points.

We note that for each eigenvalue c = c₀ an antisymmetric eigenfunction, ψ_a, also exists. This solution can be obtained like above, but with different boundary conditions at the equator. We set ψ_a(0) = 0 and , where can be chosen to control the maximum value of ψ_a, for example, unity. The other conditions are the same as before, that is, ψ_a(±R) = 0 and ψ_a continuous at the critical point. The antisymmetric solution can be expanded as ψ_a(y) = ∑_p ≥ 0A_pξ^2p + 1 for |y|< y_c and ψ_a(y) = ∑_p ≥ 1B_pξ(ξ² − 1)^p for |y|> y_c. The coefficients A_p an B_p are again obtained by recurrence. Figure 3 shows example antisymmetric eigenfunctions.

Fig. 3. Antisymmetric inviscid eigenfunctions for c = c0 and kR = 6, 8, 10, 12 and 14. The stars mark the critical points.

5. Viscous modes of oscillation

5.1. Numerical method

In order to remove the singularities at the critical latitudes, we now include viscosity. The viscosity is specified through the Reynolds number, Re, which is a free parameter in our problem. For example, Re = 700 for supergranular turbulent viscosity.

To facilitate the numerical resolution of the modified Orr-Sommerfeld equation, Eq. (17), it is convenient to introduce dimensionless quantities. We define ξ = y/R, κ = kR, and . The dimensionless eigenvalues and eigenfunctions solve the equation

(28)

where is the Laplacian, and the prime now denotes derivation with respect to ξ. For m s⁻¹, we have . We consider values of the dimensionless longitudinal wavenumber in the range 8 ≤ κ ≤ 15.

We follow the numerical approach of Orszag (1971), originally developed for the Orr-Sommerfeld eigenvalue problem. We use the Matlab package +Chebfun+ to project functions onto Chebyshev polynomials and to compute spatial derivatives analytically (Driscoll et al. 2014). This package also provides practical tools to solve differential equations and eigenvalue problems (only a few lines of codes are needed).

To obtain the symmetric solutions, we solve the above eigenvalue problem on [0, 1] with the boundary conditions

(29)

The antisymmetric solutions are obtained by setting

(30)

In both cases, the numerical solutions (eigenvalues and eigenfunctions) are complex.

5.2. Spectrum

In Fig. 4 the eigenvalues c = c_r + ic_i are shown in the complex plane for κ = 10 and Re = 300. In the figure, these eigenvalues are normalized by |c₀|> 0. All modes are stable because none of the imaginary parts of the eigenfrequencies are positive. In the complex plane, the eigenvalues are distributed along three main branches that correspond to different types of eigenfunctions. The same branches appear in the standard plane Poiseuille problem; they were called A, P, and S by Mack (1976). Example eigenfunctions are displayed in Figs. 5 and 6. The A branch, for which the eigenfunctions have large amplitudes at high latitudes, refers to the “wall modes”. The P branch refers to the “center modes”, which oscillate near the viscous layers. The S branch corresponds to the “damped modes” (Schensted 1961). Schensted (1961) showed that the A and P branches have a finite number of eigenvalues and the S branch has an infinite number of eigenvalues. She obtained approximate equations for the three branches. We labeled the modes with integers in Fig. 4, such that even integers refer to the symmetric eigenfunctions and odd integers to the antisymmetric eigenfunctions. As noted by Drazin & Reid (2004), the even and odd modes in the A branch have nearly the same eigenfrequencies. As seen in Fig. 6, the A modes have significant amplitudes only at latitudes above the viscous layers.

Fig. 4. Eigenvalues c = cr + ici in the complex plane for kR = 10 and Re = 300, normalized by the reference eigenvalue |c0|=(β − U″)/k2. The R mode is clearly identifiable (the least damped mode), as are the three families of modes: the wall modes (A family), the center modes (P family), and the damped modes (S family). The modes with symmetric eigenfunctions are shown with full circles and those with antisymmetric eigenfunctions with open circles. All modes have cr <  0 and ci <  0, i.e., they are stable and propagate in the retrograde direction. The vertical line corresponds to the phase speed of the standard Rossby wave, cr = c0. The modes in each family are labeled with integers that increase with attenuation.

Our problem differs from the standard plane Poiseuille problem through the β term. As a result, one additional mode appears in the eigenvalue spectrum (Fig. 4). This mode, which we call the R mode for an obvious reason, is symmetric and has an eigenfrequency whose real part is close to that of the classical equatorial Rossby mode, c_r ≈ c₀. The R mode has the longest lifetime in the spectrum. It is also the mode for which is smallest. As shown in Fig. 5, the real part of the R-mode eigenfunction resembles the eigenfunction of the symmetric mode found in the inviscid case (Fig. 2), except that it is smooth everywhere (no infinite derivative at the critical points). In the viscous case, the complex R-mode eigenfunctions look like chevrons in real space (Fig. 6).

Fig. 5. Eigenfunctions for the symmetric modes R, P4, S2, and A6 (top row, ψs) and for the antisymmetric modes P3, S3, and A5 (bottom row, ψa). See Fig. 4 for the position of the corresponding eigenvalues in the complex plane. The real and imaginary parts are plotted with solid and dashed lines, respectively. The modulus of the eigenfunctions is normalized to one. By choice, all imaginary parts are zero at the equator.

Fig. 6. Stream functions in real space for all the modes shown in Fig. 5. The horizontal black lines show the central latitudes of the viscous layers, λ = ±25° for kR = 10. The R mode (top left panel) is confined to the equatorial region between the viscous layers.

For some modes in the P branch, c_r is close to c₀ (P₃ and P₄), but these modes have a far shorter lifetime than the R mode (by a factor of two to three at Re = 300) and eigenfunctions that differ significantly from the observations.

5.3. R modes

Figure 7 shows the R-mode eigenfrequencies as a function of wavenumber kR for different values of the viscosity. The value of the viscosity has a rather weak effect on the dispersion relation. At fixed wavenumber, the real part of the eigenfrequency changes with Re by less than 10% over the range 100 ≥ Re ≥ 700. On the other hand, the imaginary part c_i changes significantly with the value of Re. For Re <  700, the attenuation Γ = −2kc_i increases with k. For Re = 300, we find that the theoretical mode line widths (Γ/2π = −kc_i/π in nHz) are in the range 70–100 nHz, that is, they agree reasonably well with the observed mode line widths (Γ/2π from Liang et al. 2019). The e-folding lifetime of a mode is given by τ = 2/Γ.

Fig. 7. Left panel: R-mode dispersion relations ω = kcr for Re = 100, 300, and 700 (red, blue, and green). The dashed black curve is the reference dispersion relation ω = kc0. For comparison, the frequencies of each mode m observed by Löptien et al. (2018) and Liang et al. (2019) are multiplied by the factor (m + 1)/kR and plotted at abscissa kR (this simple conversion factor is derived from the dispersion relations for classical Rossby waves in spherical and local Cartesian geometries). Right panel: plot of Γ = −2kci for Re = 100, 300, and 700. The observed full widths at half maximum for each mode m are plotted at abscissa kR for comparison.

The top row of Fig. 8 shows the R-mode stream functions for different values of the Reynolds number. The normalization is such that at ξ = 0, the real part is one and the imaginary part is zero. As Re decreases, the stream function varies more slowly around the viscous layer, and its imaginary part decreases in amplitude. The bottom row of Fig. 8 shows the vertical vorticity eigenfunctions. For Re = 700, the vorticity varies fast near the viscous layer; here ψ″ is largest.

Fig. 8. Stream function (top row) and vertical vorticity (bottom row) for R modes with kR = 10. The Reynolds number is Re = 100, 300, and 700 from left to right. The solid and dashed curves correspond to the real and imaginary parts. The shaded areas indicate the locations of the viscous critical layers for each value of Re.

6. Effect of the meridional flow

In addition to the rotational shear U, we include the effect of the meridional flow V on the R mode. The total velocity is

(31)

where the meridional flow is approximated by

(32)

with

(33)

The value of is chosen such that the maximum value of V is 15 m s⁻¹ (near latitude λ = 35°).

When the meridional flow is included, the 2D linearized Navier-Stokes equations in the equatorial β plane become

(34)

(35)

When these equations are combined, the fourth-order differential equation for the stream function (linear problem) is

(36)

with D = −k² + d²/dy². Compared to the previous problem with U only, the term in front of Dψ is now complex, and an additional term involving the first and third derivatives of the stream function appears.

We consider only the symmetric solutions and focus on the R mode. The boundary conditions are ψ′(0) = ψ‴(0) = 0 and ψ(R) = ψ′(R) = 0. Like before, we follow the procedure by Orszag (1971) to solve the eigenvalue problem. The eigenfrequencies are not affected significantly by the meridional flow. For Re = 300 and kR = 10, we find c/|c₀|= − 0.986 − 0.446i when U and V are included, compared to c/|c₀|= − 0.982 − 0.398i when only U is included. Figure 9 shows the real and imaginary parts of the R-mode eigenfunction at fixed kR = 10 for Re = 300. The meridional flow V has a small but measurable effect on the ψ eigenfunction: it is stretched toward higher latitudes, and its real part has a larger amplitude near the viscous layers.

Fig. 9. Effect of the meridional flow V on the R-mode stream function for Re = 300 and kR = 10 (red curves). For comparison, the blue curves show the stream function when only the zonal flow U is included. The real and imaginary parts correspond to the solid and dashed curves.

7. Comparison with observations

Proxauf et al. (2020) measured and characterized the vorticity eigenfunctions of the solar Rossby modes using a ring-diagram helioseismic analysis. To enable a direct comparison between observations and theory, some smoothing must be applied to the model because the observed flows have a resolution of only σ = 6° in ϕ and λ. After remapping the theoretical flows on a longitude-latitude grid, we convolve u_ϕ and u_λ with a 2D Gaussian kernel with standard deviation σ in both coordinates. From the smoothed velocities, we compute the vertical vorticity. As shown in Fig. 10, the smoothing has a significant effect on the theoretical vorticity near the viscous layer.

Fig. 10. Real (blue solid lines) and imaginary (blue dashed lines) parts of the R-mode vertical vorticity at kR = 10 after smoothing the maps of horizontal velocities with a 2D Gaussian kernel with σ = 6°. Three different values of Re are shown. For comparison, the red curves with points show the ring-diagram helioseismic observations for m = 10 near the surface (Proxauf et al. 2020).

The functional form of the real part of the observed eigenfunction was described by several parameters by Proxauf et al. (2020): a full width at half maximum W, the latitude λ₀ at which the sign changes, the latitude λ_min at which the vorticity is most negative, and the value ζ_min of the vorticity at λ_min. These parameters are provided in Table 1 for three different values of kR in three different cases: (a) smoothed theoretical vorticity when the zonal shear flow U is included, (b) smoothed theoretical vorticity when U and V are both included, and (c) observations from Proxauf et al. (2020). For Re = 300, we find that the widths for cases (a) and (b) are within ≈3° of the observed values. For (a), the zero-crossing latitude varies from λ₀ = 17° for kR = 12 to λ₀ = 21° for kR = 8. For (b), the values of λ₀ are slightly higher by ≈2°, while the observed λ₀ ≈ 28° does not vary much with k. The theoretical values of λ_min vary slightly faster with kR than λ₀, and the theory agrees better with observations for kR = 8. The observed values of ζ_min range from −0.3 to −0.1, about half of the theoretical values. However, the values of ζ_min depend on Re, with more negative values for higher values of Re. On the other hand, the values of W, λ₀, and λ_min in the table are not very sensitive to the Reynolds number for 100 ≤ Re ≤ 700.

The imaginary part of the vorticity eigenfunctions is far more difficult to measure (Proxauf et al. 2020). It is significantly different from zero for some values of m, but it is noisy, and its functional form cannot be described by a simple parametric function. According to Proxauf et al. (2020), the sign of the observed imaginary part appears to be positive for the lowest values of m and negative for m >  5. The comparison provided in Fig. 10 shows that at latitudes below the viscous layer, the amplitude of the imaginary part is far higher in the model than in the observations.

8. R-mode momentum fluxes

The complex eigenfunctions in our model imply that R modes transport a net momentum flux in latitude. It is interesting to obtain an estimate of this momentum flux (even though the model eigenfunctions have imaginary parts that differ from the observations). The purpose of this exercise is to establish whether R modes play a significant role in the balance of forces that shape differential rotation.

We consider a time-dependent zonal flow U(y, t) that evolves slowly according to the x-component of the momentum equation averaged over x (see, e.g., Geisler & Dickinson 1974),

(37)

where the prime denotes a y-derivative and angle brackets ⟨ ⋅ ⟩ denote the average over x. We used ∂_xu_x + ∂_yu_y = 0 to obtain the term involving ⟨u_xu_y⟩ on the right-hand side of the equation. The various terms in Eq. (37) have been discussed in detail by Rüdiger (1989) in spherical geometry. These terms must balance exactly to explain the steady-state differential rotation.

The horizontal Reynolds stress has two components, the first due to rotating turbulent convection (the Λ effect), and the second due to the Rossby waves,

(38)

We estimate ⟨u_xu_y⟩_R in the model for a single R mode. It is related to the stream function through

(39)

where τ is the mode lifetime mentioned in Sect. 5.3. We normalize the R-mode stream function such that the velocity amplitude at the equator is equal to its observed value,

(40)

The Reynolds stress Q_xy is plotted in Fig. 11 for R modes with Re = 300 and kR = 8, 10, and 12. We find that Q_xy <  0 in the north, below the viscous critical layer. For example, for kR = 10, Q_xy reaches the minimum value of −2 m² s⁻² at latitude 20°. This means that R modes transport angular momentum from the dissipation layer to the equator, or in other words, they reinforce latitudinal differential rotation. This is the expected result for idealized Rossby waves incident on a critical layer (see, e.g., Vallis 2006).

Fig. 11. Horizontal Reynolds stress Qxy vs. latitude for R modes with kR = 8, 10, and 12. The Reynolds number is Re = 300, and mode amplitudes were normalized according to Eq. (40). Both the zonal flow U and the meridional flow V were included to compute the mode eigenfunctions. The stars show the locations of the viscous critical layer for the different values of kR.

Summing Q_xy over several R modes would lead to a horizontal Reynolds stress that is comparable in amplitude and sign to the value reported at large spatial scales by Hathaway et al. (2013). On the other hand, Q_xy has the opposite sign and is much smaller in amplitude than the (viscous) Reynolds stress associated with convective flows at supergranulation scales (Hanasoge et al. 2016, their Fig. 10). The acceleration −∂_yQ_xy is plotted in Fig. 12 and compared to the measurements described above. The Rossby waves in our model contribute to the equatorial acceleration at a significant level.

Fig. 12. Equatorial acceleration −∂yQxy obtained by superposition of nine viscous R modes for kR = 7, 8, …15 (solid blue curve). The dot-dashed and dashed curves show −∂y⟨uxuy⟩ for supergranulation (Hanasoge et al. 2016, their Fig. 10) and larger-scale convection (Hathaway et al. 2013), respectively. For reference, the solar “torsional oscillation” has a typical amplitude of about ±5 m s−1 (see, e.g., Lekshmi et al. 2018).

9. Conclusion

Using a simple 2D setup in the β plane, we have shown that latitudinal differential rotation and viscosity must play an important role in shaping the horizontal eigenfunctions of global-scale Rossby modes. Viscous critical layers form around latitudes where c = U. We find that only one symmetric mode, which we called the R mode, has an eigenvalue whose real part is close to that of the classical sectoral Rossby mode and whose imaginary part is close to the observed value when Re ≈ 300. The real part of the vorticity eigenfunctions can be made to agree qualitatively with solar observations (unlike the imaginary part).

Treating this problem as a stability problem for a viscous Poiseuille flow in the β plane enabled us to connect our results to prior results in the fluids literature. For example, we used a well-established method to solve the Orr-Sommerfeld equation numerically, and we easily identified the known families of modes A, P, and S in the complex plane (eigenvalues and eigenfunctions). A new aspect of our work is that we identified the viscous R mode, which owes its existence to the β term in the equation. We find that the combination of the shear flow and the viscosity lead to chevron-shaped R-mode eigenfunctions, and thus to nonzero angular momentum transport by R modes. In our model, horizontal Reynolds stresses due to R modes lead to significant equatorial acceleration. Another original aspect of our work is that we studied the effect of the solar meridional flow on R modes. We found that the meridional flow affects the eigenfunctions to measurable levels. Reynolds stresses have significantly larger amplitudes when the meridional flow is included, although the meridional flow plays a much smaller role than differential rotation in shaping the eigenfunctions.

Sophisticated 3D models (but without background shear flow) also include Rossby waves as a possible mechanism to produce equatorial super-rotation in the atmospheres of planets (e.g., Liu & Schneider 2011; Read & Lebonnois 2018) and exoplanets (e.g., Showman & Polvani 2011). Clearly, the present work will have to be extended to three dimensions (see, e.g., Watts et al. 2004, for the eigenvalues in the inviscid case) to account for the radial gradients of solar rotation measured by helioseismology. A better understanding of Rossby waves will also benefit from more realistic numerical experiments (see Bekki et al. 2019). Finally, we note that the model developed here can be used to estimate the temporal changes in the Rossby wave frequencies caused by solar-cycle variations in the zonal flows (Goddard et al. 2020).

Acknowledgments

We thank Aaron Birch, Leonie Frantzen, Shravan Hanasoge (CSS), and Bastian Proxauf for useful discussions. Author contributions: L. G. proposed and designed research. L. G. and M. A. solved the inviscid problem analytically. D. F. solved the viscous problem numerically (code available at https://edmond.mpdl.mpg.de/imeji/collection/50wepJqIODg7Jg13). L. G. wrote the draft paper. All authors reviewed the final manuscript. Funding: L. G. acknowledges partial support from ERC Synergy Grant WHOLE SUN 810218 and NYUAD Institute Grant G1502. M. A. acknowledges funding from the Volkswagen Foundation and the Shota Rustaveli National Science Foundation of Georgia (SRNSFG Grant N04/46). The computational resources were provided by the German Data Center for SDO through grant 50OL1701 from the German Aerospace Center (DLR).

References

All Tables

All Figures

	Fig. 1. Parabolic zonal flow U (red curve, Eq. (4)) in the frame rotating at the equatorial rotation rate, which approximates the solar rotational velocity at the photosphere (blue curve). The horizontal black lines indicate the phase speed of Rossby waves, c0 = −(β − U″)/k2, for longitudinal wavenumbers kR = 6 and kR = 10. At critical latitudes, U = c0.
In the text

	Fig. 2. Symmetric inviscid eigenfunctions for the eigenvalue c0 = −(β − U″)/k2 and kR = 6, 8, 10, 12 and 14. The stars mark the critical points.
In the text

	Fig. 3. Antisymmetric inviscid eigenfunctions for c = c0 and kR = 6, 8, 10, 12 and 14. The stars mark the critical points.
In the text

Fig. 4. Eigenvalues c = cr + ici in the complex plane for kR = 10 and Re = 300, normalized by the reference eigenvalue |c0|=(β − U″)/k2. The R mode is clearly identifiable (the least damped mode), as are the three families of modes: the wall modes (A family), the center modes (P family), and the damped modes (S family). The modes with symmetric eigenfunctions are shown with full circles and those with antisymmetric eigenfunctions with open circles. All modes have cr <  0 and ci <  0, i.e., they are stable and propagate in the retrograde direction. The vertical line corresponds to the phase speed of the standard Rossby wave, cr = c0. The modes in each family are labeled with integers that increase with attenuation.

In the text

Fig. 5. Eigenfunctions for the symmetric modes R, P4, S2, and A6 (top row, ψs) and for the antisymmetric modes P3, S3, and A5 (bottom row, ψa). See Fig. 4 for the position of the corresponding eigenvalues in the complex plane. The real and imaginary parts are plotted with solid and dashed lines, respectively. The modulus of the eigenfunctions is normalized to one. By choice, all imaginary parts are zero at the equator.

In the text

	Fig. 6. Stream functions in real space for all the modes shown in Fig. 5. The horizontal black lines show the central latitudes of the viscous layers, λ = ±25° for kR = 10. The R mode (top left panel) is confined to the equatorial region between the viscous layers.
In the text

Fig. 7. Left panel: R-mode dispersion relations ω = kcr for Re = 100, 300, and 700 (red, blue, and green). The dashed black curve is the reference dispersion relation ω = kc0. For comparison, the frequencies of each mode m observed by Löptien et al. (2018) and Liang et al. (2019) are multiplied by the factor (m + 1)/kR and plotted at abscissa kR (this simple conversion factor is derived from the dispersion relations for classical Rossby waves in spherical and local Cartesian geometries). Right panel: plot of Γ = −2kci for Re = 100, 300, and 700. The observed full widths at half maximum for each mode m are plotted at abscissa kR for comparison.

In the text

	Fig. 8. Stream function (top row) and vertical vorticity (bottom row) for R modes with kR = 10. The Reynolds number is Re = 100, 300, and 700 from left to right. The solid and dashed curves correspond to the real and imaginary parts. The shaded areas indicate the locations of the viscous critical layers for each value of Re.
In the text

	Fig. 9. Effect of the meridional flow V on the R-mode stream function for Re = 300 and kR = 10 (red curves). For comparison, the blue curves show the stream function when only the zonal flow U is included. The real and imaginary parts correspond to the solid and dashed curves.
In the text

	Fig. 10. Real (blue solid lines) and imaginary (blue dashed lines) parts of the R-mode vertical vorticity at kR = 10 after smoothing the maps of horizontal velocities with a 2D Gaussian kernel with σ = 6°. Three different values of Re are shown. For comparison, the red curves with points show the ring-diagram helioseismic observations for m = 10 near the surface (Proxauf et al. 2020).
In the text

	Fig. 11. Horizontal Reynolds stress Qxy vs. latitude for R modes with kR = 8, 10, and 12. The Reynolds number is Re = 300, and mode amplitudes were normalized according to Eq. (40). Both the zonal flow U and the meridional flow V were included to compute the mode eigenfunctions. The stars show the locations of the viscous critical layer for the different values of kR.
In the text

Fig. 12. Equatorial acceleration −∂yQxy obtained by superposition of nine viscous R modes for kR = 7, 8, …15 (solid blue curve). The dot-dashed and dashed curves show −∂y⟨uxuy⟩ for supergranulation (Hanasoge et al. 2016, their Fig. 10) and larger-scale convection (Hathaway et al. 2013), respectively. For reference, the solar “torsional oscillation” has a typical amplitude of about ±5 m s−1 (see, e.g., Lekshmi et al. 2018).

In the text