Issue 
A&A
Volume 642, October 2020



Article Number  A178  
Number of page(s)  10  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202038525  
Published online  19 October 2020 
Effect of latitudinal differential rotation on solar Rossby waves: Critical layers, eigenfunctions, and momentum fluxes in the equatorial β plane
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: gizon@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}
Ilia State University, Kakutsa Cholokashvili Ave. 3/5, Tbilisi, 0162, Georgia
^{5}
Evgeni Kharadze Georgian National Astrophysical Observatory, Abastumani, Adigeni, 0301, Georgia
Received:
28
May
2020
Accepted:
11
August
2020
Context. Retrogradepropagating 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°–30°, and decrease at higher latitudes.
Aims. We search for an explanation that takes solar latitudinal differential rotation into account.
Methods. In the equatorial β plane, we studied the propagation of linear Rossby waves (phase speed c < 0) in a parabolic zonal shear flow, U = − U̅ ξ^{2} < 0, where U̅ = 244 m s^{−1}, and ξ is the sine of latitude.
Results. 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 to the problem to account for wave attenuation. In the viscous case, the stream functions solve a fourthorder modified OrrSommerfeld equation. Eigenvalues are complex and discrete. For reasonable values of the eddy viscosity corresponding to supergranular scales and above (Reynolds number 100 ≤ Re ≤ 700), all modes are stable. At fixed longitudinal wavenumber, the least damped mode is a symmetric mode whose real frequency is close to that of the classical Rossby mode, which we call the R mode. For Re ≈ 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).
Conclusions. 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 toward the equator.
Key words: hydrodynamics / waves / turbulence / Sun: rotation / Sun: interior / Sun: oscillations
© L. Gizon et al. 2020
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.
Open Access funding provided by Max Planck Society.
1. Introduction
In the atmosphere of Earth, Rossby (1939) waves are globalscale 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 midlatitudes 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 NavierStokes equations to model damping by turbulent convection. As shown in Sect. 2, this leads to a new equation for the stream function: an OrrSommerfeld 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,
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
By choice, U(0) = 0 at the equator. The latitudinal shear is specified by the parabolic (Poiseuille) profile,
where the amplitude is chosen such that U(y) approximates the solar surface differential rotation at low and midlatitudes. From the standard solar angular velocity profile given by Beck (2000), ΔΩ = −0.35[(y/R)^{2} + (y/R)^{4}] μrad s^{−1}, we find that the value m s^{−1} 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, c_{0} = −(β − U″)/k^{2}, for longitudinal wavenumbers kR = 6 and kR = 10. At critical latitudes, U = c_{0}. 
The 2D NavierStokes 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
3. Critical latitudes
3.1. Linear inviscid case: critical points
In the linear inviscid case, we search for wave solutions of the form
When the nonlinear and viscous terms are dropped, Eq. (8) becomes the RayleighKuo equation (Kuo 1949),
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,
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 lowlatitude region where the solution is locally oscillatory (K > 0 for y< y_{c}) from the highlatitude 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 RayleighKuo 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β, therefore waves propagate faster than in the noflow 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
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
For β = 0, this is the OrrSommerfeld equation (Orr 1907; Sommerfeld 1909). Equation (17) is a fourthorder differential equation, which requires four boundary conditions, such as ψ(±R) = 0 and ψ′(± R) = 0 for a noslip 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 lefthand 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
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^{−1} (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
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^{−1} 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
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 RayleighKuo 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 recurrence. When we set κ = 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, 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 zerocrossing 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.
Fig. 2. Symmetric inviscid eigenfunctions for the eigenvalue c_{0} = −(β − U″)/k^{2} and kR = 6, 8, 10, 12 and 14. The stars mark the critical points. 
We note that for each eigenvalue c = c_{0} 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 ≥ 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.
Fig. 3. Antisymmetric inviscid eigenfunctions for c = c_{0} 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 OrrSommerfeld equation, Eq. (17), it is convenient to introduce dimensionless quantities. We define ξ = y/R, κ = kR, and . The dimensionless eigenvalues and eigenfunctions solve the equation
where is the Laplacian, and the prime now denotes derivation with respect to ξ. For m s^{−1}, 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 OrrSommerfeld 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
The antisymmetric solutions are obtained by setting
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}> 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 = c_{r} + ic_{i} in the complex plane for kR = 10 and Re = 300, normalized by the reference eigenvalue c_{0}=(β − U″)/k^{2}. 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 c_{r} < 0 and c_{i} < 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, c_{r} = c_{0}. 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_{0}. 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 Rmode 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 Rmode eigenfunctions look like chevrons in real space (Fig. 6).
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. 
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_{0} (P_{3} and P_{4}), 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 Rmode 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 efolding lifetime of a mode is given by τ = 2/Γ.
Fig. 7. Left panel: Rmode dispersion relations ω = kc_{r} for Re = 100, 300, and 700 (red, blue, and green). The dashed black curve is the reference dispersion relation ω = kc_{0}. 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 Γ = −2kc_{i} 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 Rmode 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
where the meridional flow is approximated by
with
The value of is chosen such that the maximum value of V is 15 m s^{−1} (near latitude λ = 35°).
When the meridional flow is included, the 2D linearized NavierStokes equations in the equatorial β plane become
When these equations are combined, the fourthorder 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 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}= − 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 Rmode 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 Rmode 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 ringdiagram 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 longitudelatitude 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 Rmode 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 ringdiagram 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 λ_{0} 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 zerocrossing 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 slightly faster with kR than λ_{0}, 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, λ_{0}, and λ_{min} in the table are not very sensitive to the Reynolds number for 100 ≤ Re ≤ 700.
Parameters W, λ_{0}, λ_{min}, and ζ_{min} that characterize the functional form of the real part of the Rmode vorticity eigenfunctions.
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. Rmode 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 timedependent zonal flow U(y, t) that evolves slowly according to the xcomponent of the momentum equation averaged over x (see, e.g., Geisler & Dickinson 1974),
where the prime denotes a yderivative 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 righthand 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 steadystate 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,
We 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 Rmode 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, 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 Q_{xy} 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 −∂_{y}Q_{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 −∂_{y}Q_{xy} obtained by superposition of nine viscous R modes for kR = 7, 8, …15 (solid blue curve). The dotdashed and dashed curves show −∂_{y}⟨u_{x}u_{y}⟩ for supergranulation (Hanasoge et al. 2016, their Fig. 10) and largerscale 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 globalscale 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 wellestablished method to solve the OrrSommerfeld 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 chevronshaped Rmode 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 superrotation 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 solarcycle 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
 Balmforth, N. J., & Morrison, P. J. 1995, Ann. N.Y. Acad. Sci., 773, 80 [CrossRef] [Google Scholar]
 Beck, J. G. 2000, Sol. Phys., 191, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Bekki, Y., Cameron, R., & Gizon, L. 2019, Poster at conference “Physics at the equator: from the lab to the stars”, ENS Lyon, France, 16–18 October, https://equatorialphys.sciencesconf.org/data/Bekki_poster.pdf [Google Scholar]
 Bennett, J. R., & Young, J. A. 1971, Mon. Weather Rev., 99, 202 [CrossRef] [Google Scholar]
 Boyd, J. P. 2018, Dynamics of the Equatorial Ocean (Berlin: Springer) [CrossRef] [Google Scholar]
 Damiani, C., Cameron, R. H., Birch, A. C., & Gizon, L. 2020, A&A, 637, A65 [EDP Sciences] [Google Scholar]
 Dellar, P. J. 2011, J. Fluid Mech., 674, 174 [CrossRef] [MathSciNet] [Google Scholar]
 Drazin, P., & Howard, L. 1966, Adv. Appl. Mech., 9, 1 [CrossRef] [Google Scholar]
 Drazin, P. G., & Reid, W. H. 2004, Hydrodynamic Stability, 2nd ed. (Cambridge: Cambridge Univ. Press) [CrossRef] [Google Scholar]
 Drazin, P. G., Beaumont, D. N., & Coaker, S. A. 1982, J. Fluid Mech., 124, 439 [CrossRef] [Google Scholar]
 Driscoll, T. A., Hale, N., & Trefethen, L. N. 2014, Chebfun Guide (Oxford: Pafnuty Publications), https://www.chebfun.org/ [Google Scholar]
 Duvall, T. L., Jr, & Gizon, L. 2000, Sol. Phys., 192, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Frederiksen, J. S., & Webster, P. J. 1988, Rev. Geophys., 26, 459 [CrossRef] [Google Scholar]
 Geisler, J., & Dickinson, R. 1974, J. Atmos. Sci., 31, 946 [CrossRef] [Google Scholar]
 Gill, A. E. 1982, AtmosphereOcean Dynamics (New York: Academic Press) [Google Scholar]
 Goddard, C. R., Birch, A. C., Fournier, D., & Gizon, L. 2020, A&A, 640, L10 [CrossRef] [EDP Sciences] [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Haynes, P. H. 2003, in Encyclopedia of Atmospheric Sciences, eds. J. R. Holton, J. A. Pyle, & J. A. Curry (London: Elsevier) [Google Scholar]
 Keller, H. B. 1968, Numerical Methods for Twopoint Boundaryvalue Problems (Waltham: Blaisdell) [Google Scholar]
 Kuo, H.L. 1949, J. Atmos. Sci., 6, 105 [Google Scholar]
 Lekshmi, B., Nandy, D., & Antia, H. M. 2018, ApJ, 861, 121 [CrossRef] [Google Scholar]
 Liang, Z.C., Gizon, L., Birch, A. C., & Duvall, T. L., Jr. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liu, J., & Schneider, T. 2011, J. Atmos. Sci., 68, 2742 [CrossRef] [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 [CrossRef] [Google Scholar]
 Orr, W. M’F. 1907, Proc. R. Irish Acad., 69, A27 [Google Scholar]
 Orszag, S. A. 1971, J. Fluid Mech., 50, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Platzman, G. W. 1968, Q. J. R. Meteorol. Soc., 94, 225 [CrossRef] [Google Scholar]
 Proxauf, B., Gizon, L., Löptien, B., et al. 2020, A&A, 634, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rayleigh, L. 1879, Proc. London Math. Soc., 57, s1 [Google Scholar]
 Read, P. L., & Lebonnois, S. 2018, Ann. Rev. Earth Planet. Sci., 46, 175 [CrossRef] [Google Scholar]
 Ripa, P. 1997, J. Phys. Oceanogr., 27, 633 [CrossRef] [Google Scholar]
 Rossby, C. G. 1939, J. Mar. Res., 2, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G. 1989, Differential Rotation and Stellar Convection (Berlin: AkademieVerlag) [Google Scholar]
 Saio, H. 1982, ApJ, 256, 717 [Google Scholar]
 Schensted, I.V. 1961, PhD Thesis, The University of Michigan, Ann Arbor, USA [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, G. W., & Weiss, N. O. 1997, ApJ, 489, 960 [NASA ADS] [CrossRef] [Google Scholar]
 Sommerfeld, A. 1909, Atti del IV Congresso Internazionale dei Matematici (Roma, 6–11 Apr 1908), 116 [Google Scholar]
 Stewartson, K. 1977, Geophys. Astrophys. Fluid Dyn., 9, 185 [CrossRef] [Google Scholar]
 Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge Univ. Press) [CrossRef] [Google Scholar]
 Watts, A. L., Andersson, N., & Williams, R. L. 2004, MNRAS, 350, 927 [CrossRef] [Google Scholar]
 Webster, P. J. 1973, Mon. Weather Rev., 101, 58 [CrossRef] [Google Scholar]
All Tables
Parameters W, λ_{0}, λ_{min}, and ζ_{min} that characterize the functional form of the real part of the Rmode vorticity eigenfunctions.
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, c_{0} = −(β − U″)/k^{2}, for longitudinal wavenumbers kR = 6 and kR = 10. At critical latitudes, U = c_{0}. 

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

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

In the text 
Fig. 4. Eigenvalues c = c_{r} + ic_{i} in the complex plane for kR = 10 and Re = 300, normalized by the reference eigenvalue c_{0}=(β − U″)/k^{2}. 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 c_{r} < 0 and c_{i} < 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, c_{r} = c_{0}. The modes in each family are labeled with integers that increase with attenuation. 

In the text 
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. 

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: Rmode dispersion relations ω = kc_{r} for Re = 100, 300, and 700 (red, blue, and green). The dashed black curve is the reference dispersion relation ω = kc_{0}. 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 Γ = −2kc_{i} 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 Rmode 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 Rmode 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 ringdiagram helioseismic observations for m = 10 near the surface (Proxauf et al. 2020). 

In the text 
Fig. 11. Horizontal Reynolds stress Q_{xy} 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 −∂_{y}Q_{xy} obtained by superposition of nine viscous R modes for kR = 7, 8, …15 (solid blue curve). The dotdashed and dashed curves show −∂_{y}⟨u_{x}u_{y}⟩ for supergranulation (Hanasoge et al. 2016, their Fig. 10) and largerscale 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 
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.