Effect of latitudinal differential rotation on solar Rossby waves: Critical layers, eigenfunctions, and momentum fluxes in the equatorial $\beta$ plane

Retrograde-propagating waves of vertical vorticity with longitudinal wavenumbers between 3 and 15 have been observed on the Sun with a dispersion relation close to that of classical sectoral Rossby waves. The observed vorticity eigenfunctions are symmetric in latitude, peak at the equator, switch sign near $20^\circ$-$30^\circ$, and decrease at higher latitudes. We search for an explanation that takes into account solar latitudinal differential rotation. In the equatorial $\beta$ plane, we study the propagation of linear Rossby waves (phase speed $c<0$) in a parabolic zonal shear flow, $U = - \overline{U}\ \xi^2<0$, where $\overline{U} = 244$ m/s and $\xi$ is the sine of latitude. In the inviscid case, the eigenvalue spectrum is real and continuous and the velocity stream functions are singular at the critical latitudes where $U = c$. We add eddy viscosity in the problem to account for wave attenuation. In the viscous case, the stream functions are solution of a fourth-order modified Orr-Sommerfeld equation. Eigenvalues are complex and discrete. For reasonable values of the eddy viscosity corresponding to supergranular scales and above (Reynolds number $100 \le Re \le 700$), all modes are stable. At fixed longitudinal wavenumber, the least damped mode is a symmetric mode with a real frequency close to that of the classical Rossby mode, which we call the R mode. For $Re \approx 300$, the attenuation and the real part of the eigenfunction is in qualitative agreement with the observations (unlike the imaginary part of the eigenfunction, which has a larger amplitude in the model. Conclusion: Each longitudinal wavenumber is associated with a latitudinally symmetric R mode trapped at low latitudes by solar differential rotation. In the viscous model, R modes transport significant angular momentum from the dissipation layers towards the equator.


Introduction
In the earth atmosphere, Rossby (1939) waves are global-scale waves of radial vorticity that propagate in the direction opposite to rotation (retrograde). They find their origin in the conservation of vertical absolute vorticity, i.e. 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 P m m (sin λ) ∝ (cos λ) m where λ is 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 Sun's 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 will show that there exists a critical latitude 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 choose to work in two dimensions and in the equatorial β plane. This simplification is acceptable for wavenumbers that are large enough (say ≥ 6).
The stability and dynamics of parabolic (Poiseuille) shear flows in the presence of critical layers was summarized by, e.g., Drazin & Reid (2004). 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. Since we also wish to explain the lifetime of the modes (Löptien et al. 2018), we introduce a viscous term in the Article number, page 1 of 10 arXiv:2008.02185v1 [astro-ph.SR] 5 Aug 2020 A&A proofs: manuscript no. aanda 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 with coefficients 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). Note that in the present paper we focus on the eigenvalue problem. 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 earth 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 via 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 Orszag (1971)'s method, 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, we discuss in Sect. 8 implications of the model for angular momentum transport and the dynamics of solar differential rotation.

Waves in a sheared zonal flow: Equations of motion in the equatorial β plane
In the equatorial β plane, we study the propagation of twodimensional 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 choose the sine transformation: where φ is longitude, λ is 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 respectively equal to their φ and λ components on the sphere (Ripa 1997). The total velocity is the sum of the background zonal flow U(y)x and the horizontal wave velocity u(x, y, t) with By choice, U(0) = 0 at the equator. The latitudinal shear is specified by the parabolic (Poiseuille) profile where the amplitude U is chosen to approximate the Sun's surface differential rotation at low and mid latitudes. From the 'standard' solar angular velocity profile given by Beck (2000), ∆Ω = −0.35[(y/R) 2 + (y/R) 4 ] µrad/s, we find that the value U = 244 m/s is a good choice, see Fig. 1. The two-dimensional Navier-Stokes equations in the equatorial β plane are 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 Assuming a barotropic fluid, we can eliminate the pressure term by combining the two components of the momentum equation, to obtain

Linear inviscid case: critical points
In the linear inviscid case, we search for wave solutions of the form Dropping the nonlinear and viscous terms, Eq. (8) becomes the Rayleigh-Kuo equation (Kuo 1949): where c = ω/k is the phase speed. The above equation differs from Rayleigh's (1879) equation only through the β term. It can be rewritten as a Helmholtz equation: with 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 For equatorial propagation, K(0) = 0 is a natural choice, and we may consider the eigenvalue as an example. In our case, β − U = 1.12β, so waves propagate faster than in the no-flow case. The critical points y = ±y c where U(±y) = c 0 are given by Thus, for kR > 4.31, there are critical latitudes at λ = ±λ c , such that λ c = arcsin(4.31/kR).
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), however its first and second derivatives diverge at the critical layer (see, e.g., Haynes 2003).

Viscous critical layer
Bulk viscosity removes singularities. The linear viscous equation for ψ is For β = 0 we recognize the Orr-Sommerfeld equation (Orr 1907;Sommerfeld 1909). Equation (17) is a fourth-order differential equation, which requires four boundary conditions, e.g. ψ(±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 U − c ≈ U (y c )(y − y c ) ∼ (Uy c /R 2 )δ, so that the width of the viscous layer is approximately where 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 2 /s (Simon & Weiss 1997;Duvall, Jr. & 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.

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 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 Re * = (u max /U) −3/2 (ky c ) 1/2 , the ratio between the widths of the viscous and nonlinear critical layers is 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 larger than the Reynolds number Re ≈ 700 estimated in the previous section for solar supergranulation. Hence, it is not unreasonable 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.

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 0 (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≥0 a p ξ 2p with ξ = y/R for |y| < y c and ψ(y) = p≥1 b p (ξ 2 − 1) p for |y| > y c . The coefficients a p and b p are computed by Article number, page 3 of 10 A&A proofs: manuscript no. aanda recurrence. Setting κ = kR and ξ c = y c /R, the inner solution is given by For the outer solution, The coefficient b 1 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, thus 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 2 ψ − ψ diverges at the critical latitude and thus the comparison with the observations is difficult.
We note that for each eigenvalue c = c 0 there exists also an antisymmetric eigenfunction, ψ a . This solution can be obtained like above, but with different boundary conditions at the equator. We set ψ a (0) = 0 and ψ a (0) = ψ 0 , where ψ 0 can be chosen to control the maximum value of ψ a , for example unity. The other conditions are the same as before, i.e. ψ a (±R) = 0 and ψ a continuous at the critical point. The antisymmetric solution can be expanded as ψ a (y) = p≥0 A p ξ 2p+1 for |y| < y c and ψ a (y) = p≥1 B p ξ(ξ 2 − 1) p for |y| > y c . The coefficients A p an B p are again obtained by recurrence. Figure 3 shows example antisymmetric eigenfunctions.

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 β = βR 2 /U. The dimensionless eigenvaluesĉ = c/U and eigenfunctionsψ(ξ) = ψ(y)/(RU) solve the equation whereD = −κ 2 + d 2 /dξ 2 is the Laplacian and the prime now denotes derivation with respect to ξ. For U = 244 m/s, we havê β = 16.4. 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).
In both cases, the numerical solutions (eigenvalues and eigenfunctions) are complex.  Fig. 5. Eigenfunctions for the symmetric modes R, P 4 , S 2 , and A 6 (top row, ψ s ) and for the antisymmetric modes P 3 , S 3 , and A 5 (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.

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 | > 0. All modes are stable since 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 have been called A, P and S by Mack (1976). Example eigenfunctions are displayed in Fig. 5 and Fig. 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 labelled 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.
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 0 . The R mode is the mode with the longest lifetime in the spectrum. It is also the mode for which 1 −1 |ψ | 2 dξ is the smallest. As seen 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).
Notice that there are modes in the P branch with c r close to c 0 (P 3 and P 4 ), however these modes have a much shorter lifetime than the R mode (by a factor of two to three at Re = 300) and eigenfunctions that differ significantly from the observations. 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 small effect on the dispersion relation. At fixed wavenumber, the real part of the eigenfrequency changes with Re by less than ten percent 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 linewidths (Γ/2π = −kc i /π in nHz) are in the range 70 -100 nHz, i.e. in reasonable agreement with the ob-served mode linewidths (Γ/2π from Liang et al. 2019). Note that a mode's e-folding lifetime is given by τ = 2/Γ.

R modes
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 gets smaller in amplitude. The bottom row of Fig. 8 shows the vertical vorticity eigenfunctions. Notice the fast variations near the viscous layer for Re = 700, where ψ is largest.

Influence of the meridional flow
In addition to the rotational shear U, here we include the effect of the meridional flow V on the R mode. The total velocity is where the meridional flow is approximated by The value of V is chosen such that the maximum value of V is 15 m/s (near latitude λ = 35 • ).
Including the meridional flow, the two-dimensional linearized Navier-Stokes equations in the equatorial β plane become Combining these equations, the fourth-order differential equation for the stream function (linear problem) is with D = −k 2 + d 2 /dy 2 . 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 appear. 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 | = −0.986 − 0.446i when U and V are included, compared to c/|c 0 | = −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 influence on the ψ eigenfunction: it is stretched towards higher latitudes and its real part has a larger amplitude near the viscous layers. Proxauf et al. (2020) measured and characterized the vorticity eigenfunctions of the solar Rossby modes using ring-diagram helioseismic analysis. To enable a direct comparison between observations and theory, some smoothing must be applied to Table 1. Parameters W, λ 0 , λ min and ζ min characterizing the functional form of the real part of the R-mode vorticity eigenfunctions. The smoothed eigenfunctions from theory (Re = 300) are given for the cases when (a) only the zonal flow U is included and (b) both U and the meridional flow V are included. The theoretical values are compared with the observations from Proxauf et al. (2020 −0.28 −0.11 −0.14 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 seen in Fig. 10, the smoothing has a significant effect on the theoretical vorticity near the viscous layer.

Comparison with observations
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 λ 0 where the sign changes, the latitude λ min where 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 λ 0 = 17 • for kR = 12 to λ 0 = 21 • for kR = 8. For (b) the values of λ 0 are slightly higher by ≈ 2 • , while the observed λ 0 ≈ 28 • does not vary much with k.
The theoretical values of λ min vary with kR a little faster than λ 0 , and theory is in better agreement with observations for kR = 8. The imaginary part of the vorticity eigenfunctions is a lot more difficult to measure (Proxauf et al. 2020). It is significantly different from zero for some values of m, however 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 smallest 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 much higher in the model than in the observations.

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.
Let us consider a time-dependent zonal flow U(y, t), which evolves slowly according to the x-component of the momentum equation averaged over x (see, e.g., Geisler & Dickinson 1974): where the prime denotes a y-derivative and angle brackets · denote the average over x. We used ∂ x u x + ∂ y u y = 0 to obtain the term involving u x u 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. In our problem the horizontal Reynolds stress has two components, the first one due to rotating turbulent convection (the Λ effect) and the other one due to the Rossby waves, u x u y = u x u y R + u x u y Λ .
Let us estimate u x u y R in the model for a single R mode. It is related to the stream function through 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, 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 2 /s 2 at latitude 20 • . This means that R modes transport angular momentum from the dissipation layer to the equator, i.e. reinforce latitudinal differential rotation. This is the expected result for idealized Rossby waves incident on a critical layer (see, e.g., Vallis 2006).
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 figure 10). The acceleration −∂ y Q xy is plotted in Fig. 12 and compared to the above mentioned measurements. We see that Rossby waves in our model contribute to the equatorial acceleration at a significant level.  (Hanasoge et al. 2016, their figure 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 (see, e.g., Lekshmi et al. 2018).

Conclusion
Using a simple two-dimensional 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 to prior results in the fluids literature. For example, we used a wellestablished 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 the identification of the viscous R mode, due to the β term in the equation. We find that the combination of the shear flow and the viscosity lead to chevronshaped eigenfunctions, and thus to non-zero 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 the study of the influence 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 the 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 helio-seismology. Also, a better understanding of Rossby waves will 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 due to the solar-cycle variations in the zonal flows (Goddard et al. 2020).