Issue 
A&A
Volume 662, June 2022



Article Number  A16  
Number of page(s)  23  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202243164  
Published online  02 June 2022 
Theory of solar oscillations in the inertial frequency range: Linear modes of the convection zone
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: bekki@mps.mpg.de
^{2}
Institut für Astrophysik, GeorgAugustUniverstät Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3}
Center for Space Science, NYUAD Institute, New York University Abu Dhabi, Abu Dhabi, UAE
Received:
20
January
2022
Accepted:
9
March
2022
Context. Several types of globalscale inertial modes of oscillation have been observed on the Sun. These include the equatorial Rossby modes, criticallatitude modes, and highlatitude modes. However, the columnar convective modes (predicted by simulations and also known as banana cells or thermal Rossby waves) remain elusive.
Aims. We aim to investigate the influence of turbulent diffusivities, nonadiabatic stratification, differential rotation, and a latitudinal entropy gradient on the linear global modes of the rotating solar convection zone.
Methods. We numerically solved for the eigenmodes of a rotating compressible fluid inside a spherical shell. The model takes into account the solar stratification, turbulent diffusivities, differential rotation (determined by helioseismology), and the latitudinal entropy gradient. As a starting point, we restricted ourselves to a superadiabaticity and turbulent diffusivities that are uniform in space. We identified modes in the inertial frequency range, including the columnar convective modes as well as modes of a mixed character. The corresponding mode dispersion relations and eigenfunctions are computed for azimuthal orders of m ≤ 16.
Results. The three main results are as follows. Firstly, we find that, for m ≳ 5, the radial dependence of the equatorial Rossby modes with no radial node (n = 0) is radically changed from the traditional expectation (r^{m}) for turbulent diffusivities ≳10^{12} cm^{2} s^{−1}. Secondly, we find mixed modes, namely, modes that share properties of the equatorial Rossby modes with one radial node (n = 1) and the columnar convective modes, which are not substantially affected by turbulent diffusion. Thirdly, we show that the m = 1 highlatitude mode in the model is consistent with the solar observations when the latitudinal entropy gradient corresponding to a thermal wind balance is included (baroclinically unstable mode).
Conclusions. To our knowledge, this work is the first realistic eigenvalue calculation of the global modes of the rotating solar convection zone. This calculation reveals a rich spectrum of modes in the inertial frequency range, which can be directly compared to the observations. In turn, the observed modes can inform us about the solar convection zone.
Key words: convection / Sun: interior / Sun: rotation
© Y. Bekki et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
Based on ten years of observations from the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO), Gizon et al. (2021) discovered that the Sun supports a large number of global modes of inertial oscillations. The restoring force for these inertial modes is the Coriolis force and thus the modes have periods comparable to the solar rotation period (∼27 days). The inertial modes can potentially be used as a tool to probe the interior of the Sun because they are sensitive to properties of the deep convection zone that the p modes are insensitive to. In order to achieve this goal, we need a better understanding of the mode physics.
The low frequency modes of solar oscillation have been described in a rotating frame (angular velocity Ω_{ref}). Because the Sun is essentially symmetric about its rotation axis, the velocity of each mode in the rotating frame has the form of v(r, θ)exp[i(mϕ − ωt)], where r is the radius, θ is the colatitude, ϕ is the longitude, m is the azimuthal order, and ω is the mode eigenfrequency. Gizon et al. (2021) provide all observed eigenfrequencies ω for each m, along with the eigenfunctions (v_{θ} and v_{ϕ} at the surface) for a few selected modes.
The first family of inertial modes observed on the Sun consists of the quasitoroidal equatorial Rossby modes (Löptien et al. 2018). They are analogous to the sectoral r modes described by, for instance, Papaloizou & Pringle (1978), Smeyers et al. (1981), and Saio (1982). On the Sun, these modes have 3 ≤ m ≤ 15 with a welldefined dispersion relation close to ω = −2Ω_{ref}/(m + 1), where ω is the mode angular frequency and Ω_{ref}/2π = 453.1 nHz is the equatorial rotation rate at the surface. For positive values of m, a negative ω indicates retrograde propagation. There have been several followup studies that confirm these observations (e.g., Liang et al. 2019; Hanasoge & Mandal 2019; Proxauf et al. 2020; Hanson et al. 2020). Using a onedimensional βplane model with a parabolic shear flow and viscosity, Gizon et al. (2020a) showed that these modes, among others, are affected by differential rotation and are trapped between the critical latitudes where the phase speed of a mode is equal to the local rotational velocity. Fournier et al. (2022) extended this model to a spherical geometry using a realistic differential rotation model and found that some Rossby modes can be unstable for m ≤ 3.
Gizon et al. (2021) also report a family of modes at midlatitudes that are localized near their critical latitudes. Several tens of criticallatitude modes have been identified in the range m ≤ 10. Another family of inertial modes introduced by the Sun’s differential rotation are the highlatitude modes (Gizon et al. 2021). The highest amplitude mode (∼10−20 m s^{−1} above 50° latitude) is the m = 1 mode with northsouth antisymmetric longitudinal velocity v_{ϕ} with respect to the equator. This m = 1 mode was identified by Gizon et al. (2021) using linear calculations in twodimensional model, which are further discussed in the present paper and Fournier et al. (2022). It corresponds to the spirallike velocity feature reported at high latitudes by Hathaway et al. (2013), although it was reported as giantcell convection in that work.
The equatorialRossby and highlatitude modes involve mostly toroidal motions with a radial velocity that is small compared to the horizontal velocity components. Nontoroidal inertial modes have also been theoretically studied, mainly in the context of incompressible fluids. These modes tend to be localized onto socalled attractors, closed periodic orbits of rays reflecting off the spherical boundaries (Maas & Lam 1995; Rieutord & Valdettaro 1997, 2018; Rieutord et al. 2001; Sibgatullin & Ermanyuk 2019). They are also strongly affected by critical latitudes when differential rotation is included (e.g., Baruteau & Rieutord 2013; Guenel et al. 2016).
In numerical simulations of solarlike rotating convection, equatorial convective columns aligned with the rotation axis are prominent (e.g., Miesch et al. 2008; Bessolaz & Brun 2011; Matilsky et al. 2020). They are known as “Busse columns” (after Busse 1970), or “thermal Rossby waves”, or “banana cells” in the literature. We refer to them as “columnar convective modes” in the remainder of this paper. These convective columns propagate in the prograde direction owing either to the “topographic βeffect” originating from the geometrical curvature (e.g., Busse 2002) or to the “compressional βeffect” originating from the strong density stratification (Ingersoll & Pollard 1982; Evonuk 2008; Glatzmaier et al. 2009; Evonuk & Samuel 2012; Verhoeven & Stellmach 2014). Glatzmaier & Gilman (1981) numerically derived the dispersion relation and the radial eigenfunctions of these convective modes using a onedimensional cylinder model. They showed that the fundamental (n = 0) mode is the fastest of these prograde propagating modes with an eigenfunction that is localized near the surface, where the compressional βeffect is strongest.
In the parameter regime of the various numerical simulations, the columnar convective modes are the structures that are the most efficient to transport thermal energy upward under the rotational constraint (e.g., Gilman 1986; Miesch et al. 2000, 2008; Brun et al. 2004; Käpylä et al. 2011; Gastine et al. 2013; Hotta et al. 2015; Featherstone & Hindman 2016; Matilsky et al. 2020; Hindman et al. 2020). Furthermore, it is often argued that these convective modes play a critical role in transporting the angular momentum equatorward to maintain the differential rotation of the Sun (e.g., Gilman 1986; Miesch et al. 2000; Balbus et al. 2009). The dominant columnar convective modes seen in simulations have not been detected in the velocity field at the surface of the Sun. However, we show in this paper that some retrograde inertial modes have a mixed character and share some properties with columnar convection.
In this paper, we study the properties of the equatorial Rossby modes, the highlatitude inertial modes, and the columnar convective modes in the linear regime. We are mainly interested in the effects of turbulent diffusion, solar differential rotation, and nonadiabatic stratification on these modes. We note that the criticallatitude modes, discussed in Fournier et al. (2022), are not dealt with in depth in this paper.
Firstly, here we show that when the turbulent viscosity is above approximately 10^{12} cm^{2} s^{−1}, the equatorial Rossby modes with no radial node (n = 0) strongly depart from the expected r^{m} dependence and the radial vorticity at the surface is no longer maximum at the equator at azimuthal wavenumbers m ≳ 5. Secondly, we report a new class of modes with frequencies close to that of the classical Rossby modes. They share properties of both equatorial Rossby modes and convective modes. Thirdly, we provide a physical explanation for the properties of the m = 1 highlatitude modes in terms of the baroclinic instability due to the latitudinal entropy gradient in the convection zone.
The organization of the paper is as follows. In Sect. 2 we specify the linearized equations and solve the eigenvalue problem. The lowfrequency modes are discussed in Sect. 3 for the inviscid, adiabatically stratified, and uniformly rotating case. Then, the effects of turbulent diffusion and a nonadiabatically stratified background are discussed in Sects. 4 and 5. We discuss how the solar differential rotation and the associated baroclinicity affect the mode properties in Sect. 6. The results are summarized in Sect. 7.
2. Eigenvalue problem
In order to investigate the properties of various inertial modes in the Sun, a new numerical code has been developed. We consider the linearized fullycompressible hydrodynamic equations in a spherical coordinate (r, θ, ϕ).
2.1. Linearized equations
The linearized equations of motion, continuity, and energy conservation are:
where, v = (v_{r}, v_{θ}, v_{ϕ}) is the 1storder velocity perturbation. In this paper, we only consider the differential rotation for the mean flow and ignore meridional circulation. Thus, the background velocity is U = r sin θ(Ω − Ω_{0})e_{ϕ}. Here, Ω is a function of r and θ, denoting the rotation rate in the Sun’s convection zone, and Ω_{0} is the rotation rate of the observer’s frame. We note that, in this work, we start by analyzing the case with uniform rotation for simplicity. In this case, Ω_{0} represents the rotation rate of the unperturbed background state. For the case with the solar differential rotation, we chose to use the Carrington rotation rate Ω_{0}/2π = 456.0 nHz.
The unperturbed model is given by p_{0}, ρ_{0}, T_{0}, g, and H_{p} which are the pressure, density, temperature, gravitational acceleration, and pressure scale height of the background state. The background is assumed to be spherically symmetric and in an adiabaticallystratified hydrostatic balance. All of these variables are functions of r alone. We use the same analytical model as Rempel (2005) and Bekki & Yokoyama (2017) for the background stratification which nicely mimics the solar model S (ChristensenDalsgaard et al. 1996). The variables with subscript 1, p_{1}, ρ_{1}, and s_{1}, represent the firstorder perturbations of pressure, density, and entropy that are associated with velocity perturbation, v. Here, to close the equations, the linearized equation of state is used as follows:
where γ = 5/3 is the specific heat ratio and c_{v} denotes the specific heat at constant volume.
Although the background is approximated to be adiabatic, we can still introduce a small deviation from the adiabatic stratification in terms of the superadiabaticity δ = ∇ − ∇_{ad}, where ∇ = dlnT/dlnp is the doublelogarithmic temperature gradient. In the solar convection zone, superadiabaticity is estimated as δ ≈ 10^{−6} (e.g., Ossendrijver 2003). Also, when the solar differential rotation is included, we may add a latitudinal entropy variation ∂s_{0}/∂θ that is associated with the thermal wind balance of the differential rotation (e.g., Rempel 2005; Miesch et al. 2006; Brun et al. 2011).
We assume that the viscous stress tensor, 𝒟, is given by
where 𝒮 is the deformation tensor. We refer to Fan & Fang (2014) (their Eqs. (8)–(13)) for more detailed expressions of 𝒮_{ij} in spherical coordinates. The viscous and thermal diffusivities are denoted by v and κ respectively.
2.2. Eigenvalue problem
We assume that the ϕ and t dependence of all the perturbations v, ρ_{1}, p_{1}, and s_{1} is given by the waveform exp[i(mϕ − ωt)], where m is the azimuthal order (an integer) and ω is the complex angular frequency. With this representation, Eqs. (1)–(3) give
where C_{s} = (γp_{0}/ρ_{0})^{1/2} is the sound speed and c_{p} = γc_{v} is the constant specific heat at constant pressure. Here, the longitudinal velocity v_{ϕ}, density perturbation ρ_{1}, and entropy perturbation s_{1} are out of phase with the meridional components of velocity (v_{r} and v_{θ}) in the inviscid limit (v = κ = 0).
Equations (6)–(10) can be combined into an eigenvalue problem:
where
and M is the linear differential operator represented by the righthand side of the Eqs. (6)–(10). The operator M depends on azimuthal order m and the model parameters such as differential rotation, Ω(r, θ), superadiabaticity, δ, and diffusivities, ν and κ.
2.3. Boundary conditions
In this study, we confine our numerical domain from r_{min} = 0.71 R_{⊙} to r_{max} = 0.985 R_{⊙} in the radial direction to avoid the strong density stratification near the solar surface and gravity modes in the radiative interior. Because of viscosity, in this problem we have four secondorder (in both the radial and latitudinal directions) PDEs and one firstorder PDE. Equation (9) does not increase the order of the system as ρ_{1} can be eliminated from the system without increasing the order of the other equations. Thus, eight boundary conditions are required in the radial direction (four at the top, four at the bottom). At the top and bottom, we use impenetrable horizontal stressfree conditions for the velocity and assume there is no entropy flux (∝κ∂s_{1}/∂r) across the boundary:
All latitudes are covered in the numerical scheme, from the north pole (θ = 0) to the south pole (θ = π). We need another eight boundary conditions in the θ direction. For nonaxisymmetric cases (m ≠ 0), at the poles we impose
to make the quantities single valued. For the axisymmetric case (m = 0), at both poles we assume instead
2.4. Numerical scheme
We numerically solve the above eigenvalue problem using a finite differencing method in the meridional plane. We use a spatiallyuniform grids. The grids for v_{ϕ}, ρ_{1}, and s_{1} are staggered grids by half a grid point in radius for v_{r} and half a grid point in colatitude for v_{θ} (following Gilman 1975), as is illustrated in Fig. 1. Spatial derivatives are evaluated with a centered secondorder accurate scheme. By converting the two dimensional grid (N_{r}, N_{θ}) into one dimensional array with the size N_{r}, N_{θ} for all variables, V is defined as a one dimensional vector with size ∼5N_{r}N_{θ}. Once the boundary conditions are properly set, M can be constructed as a twodimensional complex matrix with the size approximately (5N_{r}N_{θ} × 5N_{r}N_{θ}). This method is similar to that of Guenther & Gilman (1985). In practice, each element of M can be computed by substituting a corresponding unit vector V into the righthand side of the Eqs. (6)–(10). In most of the calculations, we use the grid resolution of (N_{r}, N_{θ})=(16, 72). We have also carried out higherresolution calculations with (N_{r}, N_{θ})=(24, 180) for a uniform rotation case to check the grid convergence of the results. When the grid resolution is increased, the total number of eigenmodes increases accordingly. The additional modes have higher radial and latitudinal wavenumbers and are more finely structured. For the interpretation of the largescale modes which have been observed on the Sun, the results are converged with (N_{r}, N_{θ})=(16, 72).
Fig. 1. Layout of the staggered grid used to solve the eigenvalue equation. The grid locations where v_{ϕ}, ρ_{1}, and s_{1} are defined are denoted by red circles. The blue and green circles represent the grid locations of v_{r} and v_{θ}, respectively. The grid resolution is reduced for a visualization purpose. 
We use the LAPACK routines (Anderson et al. 1999) to numerically compute the eigenvalues and eigenvectors of M(m, ν, κ, δ, Ω), corresponding to the mode frequencies, ω, and the eigenfunctions (v_{r}, v_{θ}, v_{ϕ}, ρ_{1}, s_{1}) of linear modes in the Sun. In this study, we limit the range of azimuthal orders to m ≥ 0 and allow the real frequency to take a negative value. This means that ℜ[ω] < 0 corresponds to retrogradepropagating modes and ℑ[ω] > 0 corresponds to exponentially growing modes.
2.5. Example spectrum for uniform rotation
For each m, there are 5N_{r}N_{θ} eigensolutions with frequencies ω and eigenfunctions V. As an example, we show the typical distribution of the output eigenfrequencies in a complex plane for the case with m = 1, δ = 10^{−6} (weakly superadiabatic) and ν = κ = 2 × 10^{12} cm^{2} s^{−1} in Fig. 2. We note that the differential rotation is not included here for simplicity; the uniform rotation rate Ω is equal to the Carrington rotation rate Ω_{0}.
Fig. 2. Overview of our linear eigenmode calculation in the case of uniform rotation (Ω = Ω_{0}), a weakly superadiabatic stratification (δ = 10^{−6}), and moderate turbulent viscous and thermal diffusivities (ν = κ = 10^{11} cm^{2} s^{−1}). Shown are the results at m = 8. Upper panels: complex eigenfrequencies ω in the corotating frame. Panel a: real frequencies in the range ±10 mHz showing the acoustic modes (p modes). Panel b: zoomin focusing on the inertial range ℜ[ω] < 2Ω_{0}. Panel c: zoomin focusing on the convectivelyunstable modes (ℑ[ω]> 0). Lower panels: example eigenfunctions of pressure (acoustic), nontoroidal inertial, toroidal inertial (equatorial Rossby), and columnar convective modes, from left to right. The eigenfrequencies of these modes are highlighted by orange, green, blue, and red dots in the upper panels. 
The modes belong to one of several regions in the complex eigenfrequency spectrum. The modes seen in Fig. 2a are acoustic modes (p modes) that are slightly damped due to the viscous and thermal diffusion. On this plot, the effect of rotation is not visible to the eye. In the rest of this paper, we focus on the lowfrequency modes in the inertial frequency range. Inertial oscillations are confined within the range ℜ[ω] < 2Ω_{0} (e.g., Greenspan et al. 1968). Figure 2b shows the spectrum of inertial modes in the complex plane. The sectoral Rossby mode with no radial node (n = 0) is easy to identify by comparison with the analytical frequency ω = −2Ω_{0}/(m + 1). Owing to the slightly superadiabatic background (δ > 0), we can see that some modes have positive imaginary frequencies (ℑ[ω]> 0) at very low frequencies and thus are unstable. These convective modes are shown in Fig. 2c.
When the background is weakly subadiabatic (e.g., δ = −10^{−6}), all the modes become stable (ℑ[ω]< 0) and some inertial modes are partially mixed with gravity modes (g modes). When Ω_{0} = 0, the modes are either purely convective modes or purely g modes depending on the sign of δ as shown in Fig. 3. The frequency of the g modes depends on δ and, depending on Ω_{0}, can lie in the inertial range.
Fig. 3. Eigenfrequency spectrum in the complex plane at m = 8 for (a) δ = 10^{−6}, Ω = Ω_{0}, (b) δ = −10^{−6}, Ω = Ω_{0}, (c) δ = 10^{−6}, Ω = 0, and (d) δ = −10^{−6}, Ω = 0, respectively. Here, Ω_{0} is the Carrington rotation rate. Only inertial frequency range is shown. Upper and lower panels: cases with and without uniform rotation. Left and right panels: cases with superadiabatic and subadiabatic background. Panel a: is the same as Fig. 2b. 
3. Reference case: No diffusion, adiabatic stratification, uniform rotation
In this section, we report the results of an ideal case where turbulent viscous and thermal diffusivities are set to zero (ν = κ = 0), the background is convectively neutral (δ = 0), and no differential rotation is included (Ω(r, θ)=Ω_{0} and ∂s_{0}/∂θ = 0). We present the dispersion relations and eigenfunctions of various types of globalscale vorticity modes that might be relevant to the Sun. We go on to use the results of this ideal setup as references and, then, the effects of turbulent diffusion, nonadiabatic stratification, and differential rotation are compared to these reference results.
In the inviscid case with uniform rotation, M is self adjoint, thus the physicallymeaningful solutions must have real eigenfrequencies. We find about 10% of the eigenfrequencies to have a nonzero imaginary part; these correspond to numerical artifacts due to truncation errors, and the corresponding eigenfunctions have most of their power at high spatial frequencies. For the solutions with purely real eigenfrequencies, the eigenfunctions of v_{r} and v_{θ} have the same complex phase on each meridional plane, and those of v_{ϕ}, ρ_{1} are 90° out of phase with respect to v_{r} and v_{θ}. In presenting the results in this section, we choose a meridional plane where v_{r} and v_{θ} are real.
In the following sections, we conduct a modebymode analysis for the equatorial Rossby modes with no radial nodes (n = 0) and one radial node (n = 1), columnar convective modes (thermal Rossby waves) with both northsouth symmetries, and the highlatitude modes with both northsouth symmetries. Fundamental properties of these modes are summarized in Table 1. Their dispersion relations are presented in Table 2.
Summary of the properties of the modes of the models discussed in this paper.
Dispersion relations of the modes of the model with uniform rotation (Ω = Ω_{0}), ν = κ = 0, and δ = 0.
3.1. Equatorial Rossby modes
In this section, we discuss the equatorial Rossby modes (r modes). The modes with no radial nodes (n = 0) and one radial node (n = 1) are reported.
3.1.1. n = 0 modes
In order to extract the n = 0 equatorial Rossby mode at each m, we apply the following procedure to the computed eigenfunctions V. The latitudinal and longitudinal velocities at the surface are projected onto a basis of associated Legendre polynomials:
where l_{max} = 2N_{θ}/3 − 1 = 47. We also compute the number of radial nodes, n, of v_{θ} at the equator. We select the modes that satisfy all of the following three criteria: The l = m component of v_{θ} is dominant (a_{0}> a_{j} for all j > 0); the l = m + 1 component of v_{ϕ} is dominant (b_{1}> b_{j} for all j ≠ 1); and the number of radial nodes of v_{θ} is zero at the equator, n = 0.
Figure 4a shows the dispersion relation of the selected n = 0 equatorial Rossby modes for this ideal setup for m = 1 − 16. It should be noted that these modes are the only type of inertial modes where a simple analytical solution can be found in the inviscid, uniformlyrotating limit (e.g., Saio 1982). Therefore, we use this analytical solution to verify our code. The red points and black dashed lines represent the computed eigenfrequencies in our model and the theoreticallyexpected dispersion relation, ω = −2Ω_{0}/(m + 1), respectively. We find that the differences in the normalized frequencies are less than 10^{−2} at all m.
Fig. 4. Dispersion relation and eigenfunctions of the equatorial Rossby modes without radial nodes in the case of uniform rotation, no viscosity, and adiabatic stratification. (a) Dispersion relation from the calculated modes (red). Overplotted black dashed line represents the theoretical dispersion relation of the sectoral (l = m) Rossby modes ω = 2Ω_{0}/(m + 1). (b) Schematic illustration of flow structure of the mode with m = 6. The red and blue volume rendering shows the structure of ℜ[ζ_{r}(r, θ)exp(imϕ − iωt)]. The black solid curve shows the meridional plane at ϕ = 0 and at t = 0 where v_{r} and v_{θ} are purely real and v_{ϕ}, p_{1}, and ζ_{r} are purely imaginary. The black dashed line denotes the meridional plane at ϕ = −π/2m where v_{ϕ}, p_{1} and ζ_{r} are real. (c) Meridional cuts of the m = 2 eigenfunctions for the velocity v(r, θ)exp[i(mϕ − ωt)], the pressure p_{1}(r, θ)exp[i(mϕ − ωt)], and the radial vorticity ζ_{r}(r, θ)exp[i(mϕ − ωt)]. The solutions are shown in the meridional plane at ϕ = 0 and t = 0. The units of the color bars are m s^{−1} for the three velocity components, 10^{5} dyn cm^{−2} for the pressure, and 10^{−8} s^{−1} for the vorticity. The eigenfunctions are normalized such that the maximum of v_{θ} is 2 m s^{−1}. (d) The same as panel c but for m = 8. 
The typical flow structure of this mode is schematically illustrated in Fig. 4b where the volume rendering of the radial vorticity ζ_{r} is shown by red and blue. Figures 4c and d show the real eigenfunctions for m = 2 and 8, respectively. The eigenfunctions are normalized such that the maximum of v_{θ} is 2 m s^{−1} at the surface. The amplitude of radial velocity v_{r} is about 10^{3} times smaller than those of horizontal velocities, v_{θ} and v_{ϕ}, implying that the fluid motion is essentially toroidal. We find that using a higher resolution leads to even smaller values of v_{r}. The pressure perturbation, p_{1}, is positive (negative) where the radial vorticity, ζ_{r}, is negative (positive) in the northern (southern) hemisphere, which is consistent with the modes being in geostrophic balance. As m increases, the n = 0 equatorial Rossby modes are shifted to the surface and to the equator. The horizontal eigenfunction of ζ_{r} becomes more elongated in latitude, which means that v_{θ} becomes much stronger than v_{ϕ} to keep the mass conservation horizontally.
Figure 5a shows the radial structure of the eigenfunctions of v_{θ} at the equator for selected azimuthal orders m. Solid and dashed lines compare our results with the analytical solution v_{θ} ∝ r^{m}. It is seen that computed eigenfunctions exhibit the r^{m} dependence that agree with the analytical solutions. We also confirm the same r^{m} dependence for the eigenfunctions of v_{θ} in the middle latitudes (not shown). For higher m, the radial eigenfunction shows a slight deviation (within a few percent error) from the analytical solution. This is possibly due to the stressfree boundary condition, ∂(v_{θ}/r)/∂r = 0, at the top and bottom boundaries, which conflicts with the r^{m} dependence. Figure 5b shows the latitudinal eigenfunctions of v_{θ} at the surface. Again, an agreement can be seen between our results and the analytical solutions v_{θ} ∝ sin^{m − 1}θ.
Fig. 5. (a) Radial structure of the eigenfunction of v_{θ} at the equator for the n = 0 equatorial Rossby modes in the inviscid, uniformly rotating, and adiabatically stratified case. Overplotted dashed lines represent theoretically predicted radial dependence v_{θ} ∝ r^{m}. The eigenfunctions are normalized to unity at the surface r = r_{max}. (b) Latitudinal structure of the eigenfunction of v_{θ} at the surface. Dashed lines are the theoretical solution in the form of Legendre polynomials, v_{θ} ∝ sin^{m − 1}θ. All the eigenfunctions are normalized at the equator. 
3.1.2. n = 1 modes
The equatorial Rossby modes with one radial node (n = 1) can be selected by applying the following filters for latitudinal and longitudinal velocity eigenfunctions: we select the modes where the l = m component of v_{θ} is dominant at the surface; the l = m + 1 component of v_{ϕ} is dominant at the surface; and the number of radial nodes of v_{θ} is one at the equator.
Figure 6a shows the dispersion relation of the selected n = 1 equatorial Rossby modes for 0 ≤ m ≤ 16. It should be noted that we successfully identify the axisymmetric mode (m = 0) at ℜ[ω]= − 0.63Ω_{0}. This m = 0 mode is an equatoriallytrapped axisymmetric inertial mode. We show in Sect. 3.2.2 that this mode is connected to a progradepropagating columnar convective mode. The n = 1 Rossby modes propagate in a retrograde direction with slower phase speed than that of n = 0 Rossby modes at low m. However, for m ≥ 8, the mode frequencies become so close to those of n = 0 modes that they are almost indistinguishable.
Fig. 6. Dispersion relation and eigenfunctions of the equatorial Rossby modes with one radial node (n = 1) in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 4 is used here. 
Figure 6b shows a schematic sketch of typical flow motion of the n = 1 equatorial Rossby mode. Figures 6c and d further shows the obtained eigenfunctions of n = 1 equatorial Rossby modes plotted in the same way as in Fig. 4. It is clearly shown that v_{θ} has a nodal plane in the middle convection zone at the equator which extends in the direction of the rotation axis. One of the most striking consequences of the existence of the radial node is that substantial v_{r} is involved, owing to the radial shear of v_{θ}. Therefore, unlike the n = 0 modes, the associated fluid motions are no longer purely toroidal and become essentially threedimensional.
Figure 7a shows the radial structure of the eigenfunctions of v_{θ} at the equator for selected m. It is clearly seen that the location of the radial node shifts towards the surface as m increases. Figure 7b shows the latitudinal structure of the eigenfunctions of v_{θ} at the surface. The eigenfunctions peak at the equator and change their sign in the middle latitudes (25° −50°) and decay at higher latitudes.
Fig. 7. (a) Radial structure of the eigenfunction of v_{θ} at the equator of the n = 1 equatorial Rossby modes in the inviscid, uniformly rotating, adiabatically stratified case. The eigenfunctions are normalized to unity at the surface r = r_{max}. (b) Latitudinal structure of the eigenfunction of v_{θ} at the surface normalized at the equator. 
3.2. Columnar convective modes
In this section, we carry out a similar modebymode analysis for the columnar convective modes (thermal Rossby waves) with both hemispheric symmetries. Here, we define the northsouth symmetry based on the eigenfunction of zvorticity ζ_{z}. The “banana cell” convection pattern can be essentially regarded as the northsouth symmetric part of these convective modes. We will also show that the northsouth ζ_{z}antisymmetric modes are essentially mixed with the n = 1 equatorial modes.
3.2.1. Northsouth ζ_{z}symmetric modes
Northsouth ζ_{z}symmetric columnar convective modes can be selected by applying the following filters on the velocity eigenfunctions: the l = m component of v_{ϕ} is dominant at the surface; the l = m + 1 component of v_{θ} is dominant at the surface; the number of radial nodes of v_{r} is zero at the equator; and the number of radial nodes of v_{ϕ} is one at the equator.
Figure 8a shows the dispersion relation of the selected northsouth ζ_{z}symmetric columnar convective modes. For comparison, we overplot in black dashed line the dispersion relation derived from the onedimensional cylinder model of Glatzmaier & Gilman (1981) (their Fig. 2). Qualitatively, they both show similar features: columnar convective modes propagate in a prograde direction at all m. The modes are almost nondispersive at low m (≤7), but at higher m, the mode frequencies become almost constant at ℜ[ω] ≈ 0.85Ω_{0}. Quantitatively, our model produces the mode frequencies slightly higher (less than 10%) than that of the onedimensional cylinder model. This difference likely comes from the spherical geometry of our model: Our model takes into account both compressional and topographic βeffects that both lead to a prograde phase propagation, whereas only compressional βeffect is included in the cylinder model of Glatzmaier & Gilman (1981).
Fig. 8. Dispersion relation and eigenfunctions of the northsouth ζ_{z}symmetric columnar convective modes in the case of uniform rotation, no viscosity, and adiabatic stratification. (a) Dispersion relation of the northsouth ζ_{z}symmetric columnar convective modes in red points. For comparison, dispersion relation analytically derived using onedimensional cylinder model by Glatzmaier & Gilman (1981) is overplotted in black dashed line. (b) Schematic illustration of flow structure of the mode. Red and blue volume rendering shows the structure of ℜ[ζ_{z}(r, θ)exp(imϕ − iωt)] for m = 6 at t = 0. (c) Meridional cuts of the m = 2 eigenfunctions for the velocity v(r, θ)exp[i(mϕ − ωt)], the pressure p_{1}(r, θ)exp[i(mϕ − ωt)], and the zvorticity ζ_{r}(r, θ)exp[i(mϕ − ωt)]. The solutions are shown in the meridional plane at ϕ = 0 and t = 0 where v_{r} and v_{θ} are purely real and v_{ϕ}, p_{1} and ζ_{z} are purely imaginary. The units of velocity, pressure, and vorticity are m s^{−1}, 10^{5} dyn cm^{−2}, and 10^{−8} s^{−1}, respectively. The eigenfunctions are normalized such that maximum of v_{ϕ} is 2 m s^{−1}. (d) The same as panel c but for m = 8. 
Figures 8c and d show example eigenfunctions of the northsouth ζ_{z}symmetric columnar convective modes. The flow structure is dominantly characterized by the longitudinal velocity shear outside the tangential cylinder, leading to a strong zvorticity (where z is a coordinate in the direction of the rotation axis). Substantial radial motions are involved where v_{ϕ} converges or diverges in longitudes, as schematically illustrated in Fig. 8b. Owing to the spherical curvature of the top boundary, equatorward (poleward) latitudinal flows are involved where radial flows are outward (inward). The zvortex tubes outside the tangential cylinder are often referred to as Taylor columns or Busse columns in the geophysical context (Busse 1970, 2002) or banana cells in the solar context (Miesch et al. 2000). The pressure perturbation p_{1} is generally positive (negative) where zvorticity ζ_{z} is negative (positive), as the modes are in geostrophic balance. As m increases, the modes are more concentrated towards the surface and towards the equator.
3.2.2. Northsouth ζ_{z}antisymmetric modes
Northsouth ζ_{z}antisymmetric columnar convective modes can be selected by filtering out the eigenfunctions that satisfy the following: the l = m component of v_{θ} is dominant at the surface; the l = m + 1 component of v_{ϕ} is dominant at the surface; and the number of radial nodes of v_{θ} is one at the equator.
The dispersion relation of the ζ_{z}antisymmetric columnar convective modes is shown in Fig. 9a. For comparison, we also show the dispersion relation of the ζ_{z}symmetric modes in black dashed line. The modes propagate in a prograde direction with faster phase speed than that of the ζ_{z}symmetric modes. At high m, the dispersion relation asymptotically approaches that of the ζ_{z}symmetric modes.
Fig. 9. Dispersion relation and eigenfunctions of the northsouth ζ_{z}antisymmetric columnar convective modes in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 8 is used. Panel a: the dispersion relation of the northsouth ζ_{z}symmetric columnar convective modes is shown in black dashed line for comparison. 
Figures 9c and d show the example eigenfunctions of the northsouth ζ_{z}antisymmetric columnar convective modes. The flow structure is dominantly characterized by zvortex tubes that are antisymmetric across the equator. It should be noted that strong latitudinal motions are involved at the equator at the surface.
We find that the eigenfunctions of the m = 0 mode are the complex conjugate of the n = 1 equatorial Rossby mode, which means that these two modes are identical at m = 0 (we note the phase speed no longer matters in the case of the nonpropagating axisymmetric mode). To better illustrate this point, we show in Fig. 10 the dispersion relations of these two modes in the full (m, ℜ[ω]) domain extended to negative azimuthal orders. It is seen that the dispersion relations of these two modes connects across m = 0 and form a single continuous curve. This implies that these two modes are essentially mixed with each other: the n = 1 equatorial Rossby modes and the northsouth ζ_{z}antisymmetric columnar convective modes should be regarded as retrograde and prograde branches of the “mixed” (Rossby) modes. It is instructive to note that this mode mixing can be understood as analogous to the socalled Yanai waves which are mixed modes between retrogradepropagating Rossby modes and progradepropagating inertialgravity modes (Matsuno 1966; Vallis 2006).
Fig. 10. Dispersion relation of the “mixed modes” between the n = 1 equatorial Rossby modes (red) and the northsouth ζ_{z}antisymmetric columnar convective modes (blue) in the inviscid, uniformly rotating, adiabatically stratified case. The black points denote the axisymmetric mode at m = 0. Black solid dashed and dotdashed lines represent the dispersion relation of the n = 0 equatorial Rossby modes and northsouth ζ_{z}symmetric columnar convective modes. 
The flow structure itself of the ζ_{z}antisymmetric columnar convective mode has been recognized to be convectively unstable in the previous literature (Lorenzani & Tilgner 2001; Tilgner 2007). However, its relation to the n = 1 equatorial Rossby modes has never been reported.
3.3. High latitude modes
In this subsection, we present the eigenmodes of the highlatitude inertial modes with both hemispheric symmetries.
3.3.1. Northsouth ζ_{z}symmetric modes
To discuss these modes, it is useful to introduce a cylindrical coordinate system (ϖ,ϕ, z). In this coordinate system, the tangent cylinder is located at ϖ = r_{min}, namely, it is the cylinder aligned with the rotation axis which touches the radiative interior at the equator. Northsouth ζ_{z}symmetric highlatitude modes can be selected by applying the following criteria: The kinetic energy is predominantly inside the tangential cylinder, namely, E_{in}/E_{CZ} > 0.5, where E_{in} and E_{CZ} are the volumeintegrated kinetic energies inside the tangent cylinder and in the entire convection zone, respectively; the l = m + 1 component of v_{θ} is dominant at the bottom of the convection zone; and the number of znodes of v_{θ} is zero at ϖ = 0.5 R_{⊙}.
Figure 11a shows the dispersion relation of the northsouth ζ_{z}symmetric highlatitude modes. We find the highlatitude modes are much more dispersive than the columnar convective modes at low m. The dispersion relation is found to be roughly approximated by the nonsectoral Rossby modes’ dispersion relation with one latitudinal node (l = m + 1), as shown in the black dashed line in Fig. 11a. This is because the horizontal flows at the bottom boundary behave like the l = m + 1 (classical) Rossby modes. We note, however, that this is not regarded as the mode mixing as discussed in Sect. 3.2.2.
Fig. 11. Dispersion relation and eigenfunctions of the highlatitude modes with northsouth symmetric ζ_{z} in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 8 is used. Panel a: the dispersion relation of the l = m + 1 Rossby modes is shown in black dashed line. 
Figures 11c and d show example eigenfunctions of the ζ_{z}symmetric highlatitude modes. The fluid motion is predominantly characterized by zvortices inside the tangential cylinder in both hemispheres, as schematically illustrated in Fig. 11b. The power of ζ_{z} peaks at the tangential cylinder ϖ = r_{min}. We note that the longitudinal velocity v_{ϕ} extends slightly outside the cylinder. Again, ℑ[p_{1}]ℑ[ζ_{z}]< 0 follows from the mode being in geostrophic balance.
3.3.2. Northsouth ζ_{z}antisymmetric modes
Northsouth ζ_{z}antisymmetric highlatitude modes are selected using the following filters: The kinetic energy is predominantly inside the tangential cylinder; the l − m = 1 (or 3) component of v_{ϕ} is dominant at the bottom of the convection zone; and the number of znodes is zero for v_{θ} at ϖ = 0.5 R_{⊙}.
The example modes are presented in Fig. 12. The eigenfunctions show very similar properties of the highlatitude modes discussed in Fig. 11 except for the northsouth symmetry. It should be pointed out that there exists a latitudinal flow along the tangential cylinder and across the equator. The dispersion relation of the ζ_{z}antisymmetric highlatitude modes are found to be similar to that of l = m + 2 Rossby modes, as shown in Fig. 12a.
Fig. 12. Dispersion relation and eigenfunctions of northsouth ζ_{z}antisymmetric highlatitude modes in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 11 is used. Panel a: the dispersion relation of the l = m + 2 Rossby mode is shown in black dashed line. 
4. Effect of turbulent diffusion
So far, we have discussed the results for an inviscid case. In this section, we examine the effects of viscous and thermal diffusion arising from turbulent mixing of momentum and entropy in the Sun (e.g., Rüdiger 1989). We start our discussion by estimating the impact of the turbulent diffusion on (classical) Rossby modes. The oscillation period of the equatorial Rossby mode at the azimuthal order m is given by
On the other hand, typical diffusive timescale can be estimated as:
where l_{m} denotes the typical length scale of the Rossby mode. Figure 13 compares P_{Ro} and τ_{diff} as functions of m. Two representative values of turbulent diffusitivies in the solar convection zone ν = 10^{12} and 10^{13} cm^{2} s^{−1} are shown (e.g., Ossendrijver 2003). When P_{Ro} ≪ τ_{diff}, viscous diffusion is almost negligible. However, if P_{Ro} ≳ τ_{diff}, diffusion can have a dominant effect on the Rossby modes. For a given turbulent diffusivity ν, the critical azimuthal order, m_{crit}, can be defined as:
Fig. 13. Comparison between the oscillation periods of Rossby modes P_{Ro} and the diffusive timescale, τ_{diff}, for two representative values of turbulent diffusivities of ν = 10^{12} and 10^{13} cm^{2} s^{−1}. The horizontal black dashed line represents the length of the SDO/HMI observational record of T_{obs} ≈ 12 years. 
The Rossby modes are dominated by diffusive effects for m > m_{crit}. Figure 13 implies that the Rossby modes in the Sun are substantially affected by the turbulent diffusion, especially for m ≥ 5 − 6.
In this paper, we carry out a set of calculations of uniformlyrotating adiabatic fluid with varying diffusivities; ν = 10^{9}, 10^{10}, 10^{11}, 10^{12}, and 10^{13} cm^{2} s^{−1}. For simplicity, we fixed the Prandtl number to unity so that κ = ν. In this case, both the eigenfrequencies and eigenfunctions are complex. Figure 14 shows the eigenfrequencies of the six types of Rossby modes discussed in Sect. 3 for different viscous diffusivities in a complex plane. Figures 14a and b show the cases for m = 2 and 16, respectively. In general, the modes are damped by diffusion so that the imaginary frequencies are shifted towards more negative values. At small m (e.g., m = 2), diffusion tends to act predominantly on the columnar convective modes with both symmetries and n = 1 equatorial Rossby modes, whereas the n = 0 Rossby modes and the highlatitude modes remain almost unaffected. At large m (e.g., m = 16), however, all the modes are damped to a similar degree. We note that a strong diffusion modifies not only the imaginary part but also the real part of the mode frequencies. The computed efolding times of these modes ℑ[ω]^{−1} are shown in Fig. 15 for the two representative values of turbulent viscosity ν = 10^{11} and 10^{12} cm^{2} s^{−1}.
Fig. 14. Eigenfrequency spectra of the lowfrequency vorticity modes in a complex plane with different values of diffusivities for (a) m = 2 and (b) m = 16. Different colors represent different classes of inertial modes. Different symbols represent different values of the viscous and thermal diffusivities. In all cases, rotation is uniform and the stratification is adiabatic. 
Fig. 15. efolding lifetimes of various lowfrequency modes for a viscous diffusivity (a) ν = 10^{11} cm^{2} s^{−1} and (b) ν = 10^{12} cm^{2} s^{−1}. We note that all the modes selected here are stable modes (ℑ[ω]< 0). Different colors represent different types of inertial modes. The horizontal black dashed line shows the length of the SDO/HMI observational record (T_{obs} ≈ 12 yr as of today). In both cases, rotation is uniform and the stratification is adiabatic. The lifetimes of the convective modes and highlatitude modes are very sensitive to the radial and latitudinal entropy gradients, a point that is discussed in Sects. 5 and 6.2. 
We go on to focus on the n = 0 equatorial Rossby modes to see how eigenfunctions are affected by the viscous diffusion. Figure 16a shows the real (top row) and imaginary (bottom row) eigenfunctions of radial vorticity ζ_{r} at m = 16 for different values of viscous diffusivities ν. As ν increases, the n = 0 equatorial Rossby modes are shifted towards the base of the convection zone. This is clearly illustrated in Fig. 16b, where the absolute amplitudes of radial vorticity at the equator are shown as functions of radius. When ν becomes sufficiently large, the radial eigenfunction substantially deviates from the wellknown r^{m} dependence. This can be explained as follows: with the moderate diffusion included, the radial force balance between Coriolis force and pressure gradient force is no longer maintained. Consequently, radial flows are driven and the diffusive momentum flux becomes directed radially inward. In fact, the confinement of the n = 0 equatorial Rossby modes near the base is also seen in rotating convection simulations where the diffusion can be significantly enhanced by turbulent convection (Bekki et al., in prep.).
Fig. 16. (a) Meridional eigenfunctions of radial vorticity ζ_{r} of the n = 0 equatorial Rossby mode at m = 16 for different values of viscous diffusivities ν. Upper and lower panels: normalized real and imaginary eigenfunctions, respectively. (b) Radial eigenfunctions of ζ_{r} at the equator (normalized by their maximum amplitudes). Different colors represent different values of diffusivities. In all cases, rotation is uniform and the stratification is adiabatic. 
5. Effect of nonadiabatic stratification
In this section, the effects of nonadiabatic stratification are investigated. While theoretical models of the Sun conventionally assume a slightly positive superadiabaticity value 0 < δ ≲ 10^{−6} (e.g., Ossendrijver 2003), recent numerical simulations of solar convection imply that the lower half of the convection zone might be slightly subadiabatic (Hotta 2017; Käpylä et al. 2017, 2019; Bekki et al. 2017; Karak et al. 2018). To this end, we vary the superadiabaticity from weakly subadiabatic to weakly superadiabatic, δ = −2 × 10^{−6}, −10^{−6}, 0, 10^{−6}, 2 × 10^{−6}, while keeping the diffusivities fixed (ν = κ = 10^{12} cm^{2} s^{−1}). The solar differential rotation and latitudinal entropy gradient are not included. Since the entropy perturbation is generated by the radial flow, in this section, we focus on the (northsouth ζ_{z}symmetric) columnar convective modes where strong radial motions are involved.
Figure 17a shows the dispersion relations of the ζ_{z}symmetric columnar convective modes for different δ. As the background becomes more subadiabatic (superadiabatic), the mode frequencies become higher (lower), that is, the modes propagate in a prograde direction with faster (slower) phase speed. When δ is sufficiently large, the imaginary mode frequencies become positive, namely, the modes become convectively unstable. This is clearly manifested in Fig. 17b where the mode frequencies are plotted in a complex plane. Each point denotes each mode with the associated azimuthal order labeled nearby. The stable and unstable modes are distinguished by circles and diamonds, respectively. For δ > 0 (blue and purple), a sudden transition occurs from stable to unstable branches (at m = 5 and 4). The critical azimuthal order for this transition depends on the superadiabaticity δ via the Rayleigh number criterion for the convective instability.
Fig. 17. (a) Dispersion relations of the northsouth ζ_{z}symmetric columnar convective modes with different background superadiabaticity values δ. Different colors represent different values of superadiabaticity. Circles and diamonds denote the stable (ℑ[ω]< 0) and unstable (ℑ[ω]> 0) modes, respectively. (b) Eigenfrequencies in the complex plane. Each circle (diamond) represent a mode with azimuthal order m, which is labeled with small integers from m = 1 to 16. In all cases, rotation is uniform and the diffusivities are set ν = κ = 10^{12} cm^{2} s^{−1}. 
The changes in the dispersion relation can be understood by considering whether the buoyancy force acts as a restoring force or the opposite. Figures 18a and b present the snapshots of v_{r} and s_{1} in an equatorial plane seen from the north pole for a weakly subadiabatic background (δ = −2 × 10^{−6}) for m = 3 and 8, respectively. It is seen that the phase with positive s_{1} is always ahead of the phase with positive v_{r} in longitude, leading to a negative correlation between ℜ[v_{r}] and ℑ[s_{1}]. In a physical sense, this means that, in this case, the buoyancy force acts as an additional restoring force for progradepropagating columnar convective modes. In other words, these modes share a property of progradepropagating g modes. Consequently, the mode frequencies become higher for δ < 0. The opposite situation happens for δ > 0. Figures 18c and d show the same equatorial cuts of v_{r} and s_{1} for a weakly superadiabatic background. When m is not large enough for the convective instability to occur, it is seen that the phase with positive s_{1} is behind the phase with positive v_{r} in longitude, leading to a positive correlation between ℜ[v_{r}] and ℑ[s_{1}]. Therefore, the buoyancy force acts against the original restoring force of the compressional βeffect, which weakens the prograde propagation of columnar convective modes. As a consequence, the mode frequencies become lower for δ > 0. This effect was first studied in Gilman (1987) using a simplified cylindrical model. Figure 18d shows the case where m is sufficiently large and the mode becomes convectively unstable. It is obviously seen that the phases of v_{r} and s_{1} now coincide and they both have the same sign at each phase, leading to ⟨v_{r}s_{1}⟩> 0.
Fig. 18. Radial velocity v_{r} (upper plots) and entropy perturbation s_{1} (lower plots) of the columnar convective modes along the rotational axis displayed in the equatorial plane for subadiabatic and superadiabatic background. Panels a and b: cases with subadiabatic background δ = −2 × 10^{−6} for m = 3 and m = 8, respectively. Panels c and d: are the same plots for a superadiabatic background δ = 2 × 10^{−6}. The eigenfunctions are normalized such that the maximum radial velocity is 10 m s^{−1} at the equator. In all cases, rotation is uniform and the diffusivities are set ν = κ = 10^{12} cm^{2} s^{−1}. 
Figure 19 further shows the transport properties of thermal energy and angular momentum by convectively unstable columnar convective modes. We show the case with δ = 2 × 10^{−6} and for m = 16. Positive ⟨v_{r}T_{1}⟩ in Fig. 19a manifests that the enthalpy flux is transported upward. The Reynolds stress components ⟨v_{r}v_{ϕ}⟩ and ⟨v_{θ}v_{ϕ}⟩ are representative of the radial and latitudinal angular momentum fluxes, respectively. It is shown that the columnar convective modes can transport the angular momentum radially upward in the bulk of the convection zone and eqautorward near the surface. This agrees with the results found in the rotating convection simulation (Bekki et al., in prep.).
Fig. 19. Transport properties of thermal energy and angular momentum by the northsouth ζ_{z}symmetric columnar convective modes for m = 16. (a) Correlation between radial velocity velocity and temperature perturbation ⟨v_{r}T_{1}⟩, (b) Reynolds stress between radial and longitudinal velocities ⟨v_{r}v_{ϕ}⟩, and (c) Reynolds stress between latitudinal and longitudinal velocities ⟨v_{θ}v_{ϕ}⟩. The background is weakly superadiabatic (δ = 2 × 10^{16}), rotation is uniform, and moderate diffusivities are used (ν = κ = 10^{12} cm^{2} s^{−1}). The eigenfunctions are normalized such that the maximum radial velocity is 10 m s^{−1} at the equator. 
6. Effect of solar differential rotation
Finally, in this section, we take into account the effects of solar differential rotation. For prescribing Ω(r, θ), we use the data obtained from global helioseismology inversions from MDI and HMI (Larson & Schou 2018) as shown in Fig. 20a. We note that the observational data is truncated at r = r_{min} and r_{max}, and therefore the effects of strong radial shear layers such as tachocline and the near surface shear layer of the Sun are not included. The observing frame is chosen to be Carrington frame with the rotation rate Ω_{0}/2π = 456.0 nHz. Figure 20b shows the latitudinal profiles of the differential rotation at different depths. Horizontal dashed lines indicate the estimated phase speed of the n = 0 equatorial Rossby modes, −2Ω_{0}/[m(m+1)], for selected m values. For m > 2, there emerge critical latitudes where the phase speed of the Rossby mode matches with the differential rotation speed. As discussed in Gizon et al. (2020a) and Fournier et al. (2022), turbulent viscous diffusion is required to get rid of the singularities at the critical latitudes, leading to a formation of viscous critical layers with the typical latitudinal extent δ_{crit} given by
Fig. 20. Solar differential rotation profile used in this study. (a) Differential rotation Ω(r, θ) in a meridional plane, deduced from the global helioseismology (Larson & Schou 2018). (b) Latitudinal profiles of differential rotation at different depths. Horizontal dashed lines indicate the theoreticallyexpected phase speed of the sectoral (l = m) classical Rossby modes for selected azimuthal orders m = 2, 3, 4, 8, 16. The observing frame is chosen to be the Carrington frame rotating at Ω_{0}/2π = 456.0 nHz. 
Figure 21 shows the distribution of eigenfrequencies of the globalscale inertial modes in a complex plane for m = 2 and 16. Shown in shaded area represent the range of mode frequencies where differential rotation can have a strong impact by producing the critical layers. At higher values of m, the number of eigenmodes that are affected by differential rotation increases: In fact, most of the retrogradepropagating inertial modes are affected by critical latitudes at higher m (see Fig. 21b).
Fig. 21. Eigenfrequencies of inertial modes under the solar differential rotation in the complex plane for (a) m = 2 and (b) m = 16. The diffusivity is set to ν = 10^{12} cm^{2} s^{−1} and the background is assumed to be adiabatic, δ = 0. The shaded areas indicate the frequency range associated with the surface differential rotation, namely, m(Ω_{pole} − Ω_{0}) < ℜ[ω]< m(Ω_{eq} − Ω_{0}). 
6.1. Rossby modes with viscous critical layers
In this section, we carry out a set of calculations for ν = 10^{11} and 10^{12} cm^{2} s^{−1} with the differential rotation included to study how the equatorial Rossby modes are affected by the viscous critical layers. For the sake of simplicity, the background is set to be perfectly adiabatic and the latitudinal entropy variation ∂s_{0}/∂θ is switched off.
Figure 22a shows the dispersion relation of the equatorial Rossby modes with n = 0 (red) and n = 1 (blue) for weak (dashed) and strong (solid) viscous diffusivities, respectively. Shown in white circles, squares, and diamonds are the frequencies of the Rossby modes observed on the Sun (Löptien et al. 2018; Liang et al. 2019; Proxauf et al. 2020). The viscous diffusivity value is found to have a rather small effect on the real part of their eigenfrequencies. At m = 3, the observed frequency agrees almost perfectly with the n = 0 equatorial Rossby mode’s frequency. However, for m > 3, the observed frequencies lie in between the frequencies of n = 0 and n = 1 modes. Figure 22b shows the computed eigenfrequencies in a complex plane. Unlike the n = 1 modes, the n = 0 modes are substantially damped only for m ≥ 4, which is likely owing to the emergence of the critical latitudes that significantly modify the n = 0 modes’ eigenfunctions. The linewidths of the equatorial Rossby modes in our model are within the same order of magnitude as the observations as shown in Fig. 22b.
Fig. 22. (a) Dispersion relations of the equatorial Rossby modes for the cases with solar differential rotation. Red and blue curves represent the modes with no radial nodes (n = 0) and one radial node (n = 1), respectively. Solid and dashed lines denote the cases with weak diffusion (ν = 10^{11} cm^{2} s^{−1}) and strong diffusion (ν = 10^{12} cm^{2} s^{−1}), respectively. For comparison, the observed Rossby mode frequencies reported in Löptien et al. (2018), Liang et al. (2019), and Proxauf et al. (2020) are plotted by white hexagons and squares. All the presented frequencies are measured in the Carrington frame rotating at Ω_{0}/2π = 456.0 nHz. (b) Mode linewidths versus mode frequencies. Each point represents a mode with azimuthal order m, which is labeled with small integers from m = 1 to 16. Overplotted (open circles connected by black lines) are the mode linewidths and frequencies measured by Proxauf et al. (2020) for m = 3 to 15. In both panels, the green dots connected by line segments refer to a theoretical model for the simplified case of uniform rotation, ω/Ω_{0} = −2/(m + 1)−iE_{k}m(m + 1), where E_{k} = 4 × 10^{−4} is the Ekman number at the solar surface (Fournier et al. 2022). 
Figure 23 shows the velocity eigenfunctions of the n = 0 modes for the case with ν = 10^{12} cm^{2} s^{−1}. Figures 23a and c show meridional cuts through the eigenfunctions for m = 5 and 12, respectively. As already discussed in Sect. 4, the latitudinal velocity is confined close to the base of the convection zone. With differential rotation included, they are further trapped in the equatorial region bounded by the viscous critical layers. Unlike the uniformlyrotating case, strong radial and longitudinal flows are driven around the critical latitudes, which leads to strong concentrations of zvorticity there. Figures 23b and d show the latitudinal velocity v_{θ} (top rows) and radial vorticity ζ_{r} (bottom rows) at the top of the domain. They both have a similar chevronlike inclination towards the equator. However, ζ_{r} has prominent power peaks around the critical layers.
Fig. 23. Eigenfunctions of the equatorial Rossby modes with no radial nodes (n = 0). (a) Real (upper) and imaginary (lower) eigenfunctions of three components of velocity shown in a meridional plane for m = 5. The eigenfunctions are normalized such that the maximum latitudinal velocity is 2 m s^{−1} at the surface. The solid black line indicates the location of the critical latitudes where the phase speed of a Rossby mode matches to the differential rotation sped. (b) Horizontal eigenfunctions of latitudinal velocity v_{θ} (upper) and radial vorticity ζ_{r} (lower) at the surface r = 0.985 R_{⊙} for m = 5. The horizontal black dashed lines indicate the location of the critical latitudes at the surface. (c) and (d) are the counterparts of the panels a and b for m = 12. 
Figure 24 is the same figure as Fig. 23 but for the n = 1 equatorial Rossby modes. Unlike the n = 0 modes, the eigenfunctions of v_{θ} (and ζ_{r}) peak at the surface and at the equator. Although the critical layers exist similarly to the n = 0 modes, they are found to have a rather limited impact on the n = 1 Rossby modes.
Fig. 24. Eigenfunctions of the equatorial Rossby modes with one radial node (n = 1). Same figure notation as in Fig. 23 is used. 
To see the diffusivity dependence, we show ζ_{r} at the surface for weak (top rows) and strong viscous diffusivities (bottom rows) in Fig. 25. The left and right panels are for the n = 0 and n = 1 equatorial Rossby modes, respectively. The solid and dashed lines denote the real and imaginary parts, and the vertical red line indicates the location of the critical latitudes. The phase is defined such that ℜ[ζ_{r}] = 0 at the equator and the maximum amplitudes are normalized to unity. Substantial structure is observed associated with the viscous critical layers. In general, this structure becomes broader and weaker as the viscosity ν is increased. It is also seen that amplitudes of the imaginary parts of the eigenfunctions are larger for n = 0 modes than for n = 1 modes.
Fig. 25. Radial vorticity ζ_{r} eigenfunctions at the surface (r = 0.985 R_{⊙}) for m = 8. Left and right panels: cases for the equatorial Rossby modes with no radial nodes (n = 0) and with one radial node (n = 1). Upper and lower panels: cases with weak diffusion (ν = 10^{11} cm^{2} s^{−1}) and strong diffusion (ν = 10^{12} cm^{2} s^{−1}). Black solid and dashed lines represent real and imaginary eigenfunctions, respectively. The real part of the eigenfunctions are defined to be zero at the equator. The vertical red lines denote the location of critical latitudes where the phase speed of a Rossby mode is equal to the differential rotation velocity. 
Next, let us examine the impact of the net angular momentum transport by the equatorial Rossby modes under the influences of solar differential rotation. Figures 26a and b show the Reynolds stress components ⟨v_{r}v_{θ}⟩ and ⟨v_{θ}v_{θ}⟩ for n = 0 modes at m = 8. The Reynolds stresses become substantially nonzero near the viscous critical layers. It is striking that even n = 0 mode, which in the case of uniform rotation is toroidal and nonconvective, can transport the angular momentum radially upward around the viscous critical layers. Latitudinally, the angular momentum is transported equatorward in both hemispheres. Figure 26c shows the ⟨v_{θ}v_{θ}⟩ at the surface for all m. It is seen that the correlations become small as m increases because the n = 0 modes are more and more confined closer to the base of the convection zone. The counterparts for n = 1 modes are shown in Figs. 26d–f. It is clear that the n = 1 modes also transport angular momentum radially upward and equatorward at higher m. However, unlike the n = 0 modes, the Reynolds stress ⟨v_{θ}v_{θ}⟩ peaks slightly below the surface. Therefore, the correlation at the surface becomes more prominent as m increases, as shown in Fig. 26f.
Fig. 26. Reynolds stress components (a), (d) ρ_{0}⟨v_{r}v_{ϕ}⟩ and (b), (e) ρ_{0}⟨v_{θ}v_{ϕ}⟩ associated with the equatorial Rossby modes for m = 8. The units are g cm^{−1} s^{−2}. Black solid lines denote the location of critical latitudes at each height. The eigenfunctions are normalized such that the maximum horizontal velocity at the top boundary is 2 m s^{−1}. Panels c and f: horizontal Reynolds stress averaged over radius where the overbar denotes the radial average. Different colors represent different azimuthal orders. Upper and lower panels: cases for n = 0 modes and n = 1 Rossby modes, respectively. 
It is instructive to examine how significant the angular momentum transport by these equatorial Rossby modes can be in the Sun. To this end, we consider the socalled gyroscopic pumping equation (e.g., Elliott et al. 2000; Miesch & Hindman 2011)
where F_{RS}, F_{MC}, and F_{VD} are the angular momentum fluxes transported by the Reynolds stress, meridional circulation, and turbulent viscous diffusion, respectively. They are defined by
where v_{m} is the meridional flow. Figure 27 shows the each term of the latitudinal component of the Eq. (22) averaged over radius. The eigenfunctions are normalized such that the maximum horizontal velocity amplitude at the surface is 2 m s^{−1}, as inferred from observations (Löptien et al. 2018). To estimate F_{MC, θ} (black dotdashed line), we use the observational meridional circulation data obtained by Gizon et al. (2020b). For F_{VD, θ} (black dashed line), we assume the spatiallyuniform viscosity of ν = 10^{12} cm^{2} s^{−1}. It is shown that the equatorward angular momentum transport by the Reynolds stress is balanced by the poleward transport by meridional flow and by turbulent diffusion. The amplitude of F_{RS, θ} associated with n = 1 modes are found to be almost negligible, whereas that of n = 0 modes is substantial and accounts for about 30 − 40% of the other two contributions F_{MC, θ} + F_{VD, θ}. The difference between the n = 0 and n = 1 modes comes from that fact that the velocity eigenfunctions of the n = 1 modes peak at the surface, whereas those of the n = 0 modes peak near the base. Therefore, when the eigenfunctions are normalized by the surface velocity speed, only n = 0 modes become important for the convection zone dynamics. Some caution must be given here as the eigenfunctions can also be highly sensitive to various model parameters (such as ν and δ), and thus, a different set of parameters might lead to a different angular momentum balance. Furthermore, the model assumes that the diffusivities are uniform and isotropic, which will also affect the eigenfunctions. Nonetheless, it is suggested that the equatorial Rossby modes might potentially play a role in transporting the angular momentum equatorward in the Sun.
Fig. 27. Radially averaged latitudinal angular momentum fluxes. Solid, dotdashed, and dashed lines represent the angular momentum fluxes associated with the Reynolds stress of the equatorial Rossby modes F_{RS, θ}, advection by meridional circulation F_{MC, θ}, and diffusion by turbulent viscosity F_{VD, θ}, defined by the Eqs. (23)–(25). Red and blue lines represent the modes with no radial nodes (n = 0) and with one radial node (n = 1), respectively. Here, eigenfunctions are normalized such that the maximum horizontal velocity amplitude at the surface is 2 m s^{−1}. 
6.2. Effect of baroclinicity on highlatitude inertial modes
We assume the solar differential rotation is in thermal wind balance; where the deviation from the TaylorProudman’s state is balanced by the latitudinal entropy variation (e.g., Rempel 2005; Miesch et al. 2006; Brun et al. 2011). In other words, the solar convection zone is essentially baroclinic. Since the highlatitude modes are located at high latitudes, they are subject to the imposed baroclinicity in the convection zone and potentially become unstable (Knobloch & Spruit 1982; Spruit & Knobloch 1984; Kitchatinov 2013; Gilman & Dikpati 2014).
In this section, we study the effect of baroclinicity in the convection zone on the highlatitude inertial modes by varying the amplitude of the imposed latitudinal entropy gradient. Of particular interest is this effect on the m = 1 mode with northsouth antisymmetric ζ_{z}. Here, we assume the latitudinal dependence of the background entropy profile is
where Δ_{θ}s = s_{0,eq} − s_{0,pole}(< 0) represents the entropy difference between the cooler equator and the hotter poles. For simplicity, the radial dependence is ignored (s_{0} is uniform in radius and thus convectively neutral). We use moderately viscous and thermal diffusivities ν = κ = 10^{12} cm^{2} s^{−1}.
Figure 28a shows the eigenfrequencies of the m = 1 northsouth ζ_{z}antisymmetric highlatitude modes in a complex plane with varying Δ_{θ}s from 0 to 2000 erg g^{−1} K^{−1}. It is shown that as the baroclinicity is increased, the modes become unstable (ℑ[ω]> 0). In this sense, these modes can also be called baroclinic (Rossby) modes. Figures 28b–f show the eigenfunctions of v_{ϕ} both at the surface and at the central meridian for different values of Δ_{θ}s. It is clearly seen that, as Δ_{θ}s increases and the highlatitude mode becomes more and more baroclinically unstable, it begins to exhibit a spiralling flow structure around the poles. The spatial extent and the tilt of this spiral agree strikingly well with the observations. We refer to Hathaway & Upton (2021) and Gizon et al. (2021) for more details.
Fig. 28. (a) Growth rate versus frequency of the ζ_{z}antisymmetric highlatitude mode with m = 1, for different values of the latitudinal entropy difference Δ_{θ}s=s_{0, pole} − s_{0, eq}. The red star is for the case of a realistic latitudinal entropy gradient that depends on position, estimated according to Eq. (27). Realistic solar differential rotation is included. The background stratification is adiabatic (δ = 0) and the diffusivities are set to ν = κ = 10^{12} cm^{2} s^{−1}. (b)–(f) Eigenfunctions of the longitudinal velocity v_{ϕ} at the surface (r = 0.985 R_{⊙}) and at the central meridian for some selected Δs_{θ}. The eigenfunctions are normalized so that the maximum flow amplitudes at the surface is v_{ϕ} = 10 m s^{−1}. 
In order to assess if the baroclinicity in the Sun is large enough for the baroclinic instability to occur, we estimate the latitudinal entropy variation using the helioseismicallyconstrained differential rotation profile using:
With this realistic baroclinicity included in our model (Eq. (10)), we find that the m = 1 highlatitude mode is selfexcited: the growth rate is ℑ[ω]/2π = 14.1 nHz, which translates into the growing timescale of 4.3 months. This may explain why the highlatitude flow feature on the Sun has a much larger flow amplitude than the equatorial Rossby modes. Its mode frequency is ℜ[ω]/2π = −90.9 nHz (measured in the Carrington frame), which is close to the observed propagation frequency of the highlatitude flow feature of −86.3 nHz (Gizon et al. 2021). The eigenfunctions of this m = 1 mode are shown in Fig. 29. The mode is characterized by its dominant zvortical motion and is quasitoroidal (the vertical flow is about ten times weaker than the horizontal ones). It is clearly seen that, unlike the case without baroclinicity, a strong entropy perturbation is associated with this mode. The mode is strongly localized near the poles. Although this m = 1 mode has successfully been detected on the Sun, we find that about 30% of the total kinetic energy exists above 80° latitude, which is the observational upper limit of the ringdiagram measurements presented by Gizon et al. (2021). This means that the observations may miss a fraction of the mode power.
Fig. 29. Eigenfunctions of v_{r}, v_{θ}, v_{ϕ}, and s_{1} of the m = 1 northsouth ζ_{z}antisymmetric highlatitude inertial mode with the solar differential rotation and the corresponding latitudinal entropy gradient (Eq. (27)) included. The background stratification is adiabatic (δ = 0) and the diffusivities are set to ν = κ = 10^{12} cm^{2} s^{−1}. The eigenfunctions are normalized such that the maximum v_{ϕ} is 10 m s^{−1} at the surface. 
Using a linear model of Boussinesq convection in a uniformlyrotating spherical shell, Gilman (1975) found convectivelyunstable modes near the poles with a spiral structure (see his Fig. 17). As with the highlatitude modes we have found, the modes in his study lie mostly inside the tangent cylinder. However, those modes are convectively driven, whereas ours are baroclinically driven due to the solar differential rotation and the latitudinal entropy variation.
7. Summary
In this paper, we present a linear modal analysis of the oscillations of the solar convection zone at low frequencies. We report the dispersion relations and eigenfunctions of the equatorial Rossby modes without a radial node (n = 0) and with one radial node (n = 1), along with the columnar convective modes and the highlatitude modes, both with different northsouth symmetries. We find “mixed Rossby modes” that share properties of the n = 1 equatorial Rossby modes and the ζ_{z}antisymmetric columnar convective modes.
We studied the effects of the turbulent diffusion and the solar differential rotation on the equatorial Rossby modes. Our main findings are summarized as follows. One effect of turbulent diffusion is to radically change the radial force balance of the n = 0 equatorial Rossby modes. The modes are confined closer to the base of the convection zone and their eigenfunctions strongly deviate from the wellknown r^{m} radial dependence. When the solar differential rotation is taken into account, viscous critical layers are formed in latitudes where the phase speed of the equatorial Rossby mode is equal to the differential rotation speed. Strong radial and longitudinal flows are present in the viscous critical layers and the eigenfunctions are complex, implying nonzero Reynolds stresses. We also find that, unlike the n = 0 equatorial Rossby modes, the “mixed modes” are almost unaffected by the presence of solar differential rotation and strong viscous diffusivity. The retrograde frequencies of the observed Rossby modes of the Sun have values in between the model eigenfrequencies of the n = 0 and n = 1 modes for m ≥ 5 (see Fig. 22a).
We further demonstrate that the dispersion relations of the columnar convective modes are very sensitive to the background superadiabaticity δ. These modes are convectively unstable and transport thermal energy and angular momentum when δ > 0. We have also shown that the m = 1 highlatitude mode is substantially modified by a latitudinal entropy gradient. When the latitudinal entropy gradient exists, the mode is baroclinically unstable and the eigenfunction at the surface matches the observations with the correct geometry, including the correct sense for the spiral seen in v_{ϕ}. These results imply that the observations of the solar inertial modes can be used to measure the degree of nonadiabaticity in the Sun, as proposed by Gilman (1987).
Several simplifying assumptions were made in this study. For instance, the viscous and thermal diffusivities, ν and κ, and the superadiabaticity δ were all assumed to be spatially uniform, which is not realistic. Moreover, we set the bottom and top boundaries at (r_{min}, r_{max})=(0.71 R_{⊙}, 0.985 R_{⊙}) and thus both the tachocline and the nearsurface shear layer of the Sun were excluded from our model. Future work will be to include the radiative interior and the photosphere and to allow for a radial dependence of ν, κ and δ. In addition, it will be important to compare the present results to modes extracted from threedimensional numerical simulations of rotating convection in the strongly nonlinear regime. The aim is to have a physical understanding of all the modes in the lowfrequency spectrum and thus a reliable identification of the observed modes, including the criticallatitude modes.
Acknowledgments
We thank A.C. Birch and B. Proxauf for helpful discussions. Y.B. did most of the work, R.H.C. and L.G. provided supervision. Y.B. is a member of the International MaxPlanck Research School for Solar System Science at the University of Göttingen. We acknowledge support from ERC Synergy Grant WHOLE SUN 810218 and the hospitality of the Institut Pascal in March 2022. Y.B. is the beneficiary of a longterm scholarship program for degreeseeking graduate students abroad from the Japan Student Services Organization (JASSO).
References
 Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users’ Guide, 3rd edn. (Philadelphia, PA: Society for Industrial and Applied Mathematics) [CrossRef] [Google Scholar]
 Balbus, S. A., Bonart, J., Latter, H. N., & Weiss, N. O. 2009, MNRAS, 400, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Rieutord, M. 2013, J. Fluid Mech., 719, 47 [Google Scholar]
 Bekki, Y., & Yokoyama, T. 2017, ApJ, 835, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Bessolaz, N., & Brun, A. S. 2011, ApJ, 728, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1970, J. Fluid Mech., 44, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 2002, Phys. Fluids, 14, 1301 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [Google Scholar]
 Elliott, J. R., Miesch, M. S., & Toomre, J. 2000, ApJ, 533, 546 [NASA ADS] [CrossRef] [Google Scholar]
 Evonuk, M. 2008, ApJ, 673, 1154 [NASA ADS] [CrossRef] [Google Scholar]
 Evonuk, M., & Samuel, H. 2012, Earth Planet. Sci. Lett., 317, 1 [CrossRef] [Google Scholar]
 Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Featherstone, N. A., & Hindman, B. W. 2016, ApJ, 830, L15 [Google Scholar]
 Fournier, D., Gizon, L., & Hyest, L. 2022, A&A, in press http://dx.doi.org/10.1051/00046361/202243473 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156 [Google Scholar]
 Gilman, P. A. 1975, J. Atm. Sci., 32, 1331 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1986, in Physics of the Sun, eds. P. A. Sturrock, T. E. Holzer, D. M. Mihalas, & R. K. Ulrich, 1, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1987, ApJ, 318, 904 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P., & Dikpati, M. 2014, ApJ, 787, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Fournier, D., & Albekioni, M. 2020a, A&A, 642, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020b, Science, 368, 1469 [Google Scholar]
 Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Glatzmaier, G. A., & Gilman, P. A. 1981, ApJS, 45, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Glatzmaier, G., Evonuk, M., & Rogers, T. 2009, Geophys. Astrophys. Fluid Dyn., 103, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Greenspan, H., Batchelor, C., Ablowitz, M., et al. 1968, The Theory of Rotating Fluids, Cambridge Monographs on Mechanics (Cambridge: Cambridge University Press) [Google Scholar]
 Guenel, M., Baruteau, C., Mathis, S., & Rieutord, M. 2016, A&A, 589, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guenther, D. B., & Gilman, P. A. 1985, ApJ, 295, 195 [Google Scholar]
 Hanasoge, S., & Mandal, K. 2019, ApJ, 871, L32 [Google Scholar]
 Hanson, C. S., Gizon, L., & Liang, Z.C. 2020, A&A, 635, A109 [EDP Sciences] [Google Scholar]
 Hathaway, D. H., & Upton, L. A. 2021, ApJ, 908, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Hindman, B. W., Featherstone, N. A., & Julien, K. 2020, ApJ, 898, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 803, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Ingersoll, A. P., & Pollard, D. 1982, Icarus, 52, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011, Astron. Nachr., 332, 883 [CrossRef] [Google Scholar]
 Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23 [CrossRef] [Google Scholar]
 Käpylä, P. J., Viviani, M., Käpylä, M. J., Brandenburg, A., & Spada, F. 2019, Geophys. Astrophys. Fluid Dyn., 113, 149 [CrossRef] [Google Scholar]
 Karak, B. B., Miesch, M., & Bekki, Y. 2018, Phys. Fluids, 30, 046602 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L. 2013, Astron. Lett., 39, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
 Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29 [Google Scholar]
 Liang, Z.C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
 Lorenzani, S., & Tilgner, A. 2001, J. Fluid Mech., 447, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Maas, L. R. M., & Lam, F.P. A. 1995, J. Fluid Mech., 300, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Matilsky, L. I., Hindman, B. W., & Toomre, J. 2020, ApJ, 898, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuno, T. 1966, J. Meteorol. Soc. Japan. Ser. II, 44, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Hindman, B. W. 2011, ApJ, 743, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [CrossRef] [Google Scholar]
 Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618 [Google Scholar]
 Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
 Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 182, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Proxauf, B., Gizon, L., Löptien, B., et al. 2020, A&A, 634, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., & Valdettaro, L. 1997, J. Fluid Mech., 341, 77 [Google Scholar]
 Rieutord, M., & Valdettaro, L. 2018, J. Fluid Mech., 844, 597 [Google Scholar]
 Rieutord, M., Georgeot, B., & Valdettaro, L. 2001, J. Fluid Mech., 435, 103 [Google Scholar]
 Rüdiger, G. 1989, Differential Rotation and Stellar Convection (Berlin: Akademie Verlag) [Google Scholar]
 Saio, H. 1982, ApJ, 256, 717 [Google Scholar]
 Sibgatullin, I. N., & Ermanyuk, E. V. 2019, J. Appl. Mech. Tech. Phys., 60, 284 [NASA ADS] [CrossRef] [Google Scholar]
 Smeyers, P., Craeynest, D., & Martens, L. 1981, Ap&SS, 78, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C., & Knobloch, E. 1984, A&A, 132, 89 [NASA ADS] [Google Scholar]
 Tilgner, A. 2007, in Treatise on Geophysics, ed. G. Schubert (Amsterdam: Elsevier), 8.07, 207 [NASA ADS] [Google Scholar]
 Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Verhoeven, J., & Stellmach, S. 2014, Icarus, 237, 143 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Dispersion relations of the modes of the model with uniform rotation (Ω = Ω_{0}), ν = κ = 0, and δ = 0.
All Figures
Fig. 1. Layout of the staggered grid used to solve the eigenvalue equation. The grid locations where v_{ϕ}, ρ_{1}, and s_{1} are defined are denoted by red circles. The blue and green circles represent the grid locations of v_{r} and v_{θ}, respectively. The grid resolution is reduced for a visualization purpose. 

In the text 
Fig. 2. Overview of our linear eigenmode calculation in the case of uniform rotation (Ω = Ω_{0}), a weakly superadiabatic stratification (δ = 10^{−6}), and moderate turbulent viscous and thermal diffusivities (ν = κ = 10^{11} cm^{2} s^{−1}). Shown are the results at m = 8. Upper panels: complex eigenfrequencies ω in the corotating frame. Panel a: real frequencies in the range ±10 mHz showing the acoustic modes (p modes). Panel b: zoomin focusing on the inertial range ℜ[ω] < 2Ω_{0}. Panel c: zoomin focusing on the convectivelyunstable modes (ℑ[ω]> 0). Lower panels: example eigenfunctions of pressure (acoustic), nontoroidal inertial, toroidal inertial (equatorial Rossby), and columnar convective modes, from left to right. The eigenfrequencies of these modes are highlighted by orange, green, blue, and red dots in the upper panels. 

In the text 
Fig. 3. Eigenfrequency spectrum in the complex plane at m = 8 for (a) δ = 10^{−6}, Ω = Ω_{0}, (b) δ = −10^{−6}, Ω = Ω_{0}, (c) δ = 10^{−6}, Ω = 0, and (d) δ = −10^{−6}, Ω = 0, respectively. Here, Ω_{0} is the Carrington rotation rate. Only inertial frequency range is shown. Upper and lower panels: cases with and without uniform rotation. Left and right panels: cases with superadiabatic and subadiabatic background. Panel a: is the same as Fig. 2b. 

In the text 
Fig. 4. Dispersion relation and eigenfunctions of the equatorial Rossby modes without radial nodes in the case of uniform rotation, no viscosity, and adiabatic stratification. (a) Dispersion relation from the calculated modes (red). Overplotted black dashed line represents the theoretical dispersion relation of the sectoral (l = m) Rossby modes ω = 2Ω_{0}/(m + 1). (b) Schematic illustration of flow structure of the mode with m = 6. The red and blue volume rendering shows the structure of ℜ[ζ_{r}(r, θ)exp(imϕ − iωt)]. The black solid curve shows the meridional plane at ϕ = 0 and at t = 0 where v_{r} and v_{θ} are purely real and v_{ϕ}, p_{1}, and ζ_{r} are purely imaginary. The black dashed line denotes the meridional plane at ϕ = −π/2m where v_{ϕ}, p_{1} and ζ_{r} are real. (c) Meridional cuts of the m = 2 eigenfunctions for the velocity v(r, θ)exp[i(mϕ − ωt)], the pressure p_{1}(r, θ)exp[i(mϕ − ωt)], and the radial vorticity ζ_{r}(r, θ)exp[i(mϕ − ωt)]. The solutions are shown in the meridional plane at ϕ = 0 and t = 0. The units of the color bars are m s^{−1} for the three velocity components, 10^{5} dyn cm^{−2} for the pressure, and 10^{−8} s^{−1} for the vorticity. The eigenfunctions are normalized such that the maximum of v_{θ} is 2 m s^{−1}. (d) The same as panel c but for m = 8. 

In the text 
Fig. 5. (a) Radial structure of the eigenfunction of v_{θ} at the equator for the n = 0 equatorial Rossby modes in the inviscid, uniformly rotating, and adiabatically stratified case. Overplotted dashed lines represent theoretically predicted radial dependence v_{θ} ∝ r^{m}. The eigenfunctions are normalized to unity at the surface r = r_{max}. (b) Latitudinal structure of the eigenfunction of v_{θ} at the surface. Dashed lines are the theoretical solution in the form of Legendre polynomials, v_{θ} ∝ sin^{m − 1}θ. All the eigenfunctions are normalized at the equator. 

In the text 
Fig. 6. Dispersion relation and eigenfunctions of the equatorial Rossby modes with one radial node (n = 1) in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 4 is used here. 

In the text 
Fig. 7. (a) Radial structure of the eigenfunction of v_{θ} at the equator of the n = 1 equatorial Rossby modes in the inviscid, uniformly rotating, adiabatically stratified case. The eigenfunctions are normalized to unity at the surface r = r_{max}. (b) Latitudinal structure of the eigenfunction of v_{θ} at the surface normalized at the equator. 

In the text 
Fig. 8. Dispersion relation and eigenfunctions of the northsouth ζ_{z}symmetric columnar convective modes in the case of uniform rotation, no viscosity, and adiabatic stratification. (a) Dispersion relation of the northsouth ζ_{z}symmetric columnar convective modes in red points. For comparison, dispersion relation analytically derived using onedimensional cylinder model by Glatzmaier & Gilman (1981) is overplotted in black dashed line. (b) Schematic illustration of flow structure of the mode. Red and blue volume rendering shows the structure of ℜ[ζ_{z}(r, θ)exp(imϕ − iωt)] for m = 6 at t = 0. (c) Meridional cuts of the m = 2 eigenfunctions for the velocity v(r, θ)exp[i(mϕ − ωt)], the pressure p_{1}(r, θ)exp[i(mϕ − ωt)], and the zvorticity ζ_{r}(r, θ)exp[i(mϕ − ωt)]. The solutions are shown in the meridional plane at ϕ = 0 and t = 0 where v_{r} and v_{θ} are purely real and v_{ϕ}, p_{1} and ζ_{z} are purely imaginary. The units of velocity, pressure, and vorticity are m s^{−1}, 10^{5} dyn cm^{−2}, and 10^{−8} s^{−1}, respectively. The eigenfunctions are normalized such that maximum of v_{ϕ} is 2 m s^{−1}. (d) The same as panel c but for m = 8. 

In the text 
Fig. 9. Dispersion relation and eigenfunctions of the northsouth ζ_{z}antisymmetric columnar convective modes in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 8 is used. Panel a: the dispersion relation of the northsouth ζ_{z}symmetric columnar convective modes is shown in black dashed line for comparison. 

In the text 
Fig. 10. Dispersion relation of the “mixed modes” between the n = 1 equatorial Rossby modes (red) and the northsouth ζ_{z}antisymmetric columnar convective modes (blue) in the inviscid, uniformly rotating, adiabatically stratified case. The black points denote the axisymmetric mode at m = 0. Black solid dashed and dotdashed lines represent the dispersion relation of the n = 0 equatorial Rossby modes and northsouth ζ_{z}symmetric columnar convective modes. 

In the text 
Fig. 11. Dispersion relation and eigenfunctions of the highlatitude modes with northsouth symmetric ζ_{z} in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 8 is used. Panel a: the dispersion relation of the l = m + 1 Rossby modes is shown in black dashed line. 

In the text 
Fig. 12. Dispersion relation and eigenfunctions of northsouth ζ_{z}antisymmetric highlatitude modes in the inviscid, uniformly rotating, adiabatically stratified case. The same notation as Fig. 11 is used. Panel a: the dispersion relation of the l = m + 2 Rossby mode is shown in black dashed line. 

In the text 
Fig. 13. Comparison between the oscillation periods of Rossby modes P_{Ro} and the diffusive timescale, τ_{diff}, for two representative values of turbulent diffusivities of ν = 10^{12} and 10^{13} cm^{2} s^{−1}. The horizontal black dashed line represents the length of the SDO/HMI observational record of T_{obs} ≈ 12 years. 

In the text 
Fig. 14. Eigenfrequency spectra of the lowfrequency vorticity modes in a complex plane with different values of diffusivities for (a) m = 2 and (b) m = 16. Different colors represent different classes of inertial modes. Different symbols represent different values of the viscous and thermal diffusivities. In all cases, rotation is uniform and the stratification is adiabatic. 

In the text 
Fig. 15. efolding lifetimes of various lowfrequency modes for a viscous diffusivity (a) ν = 10^{11} cm^{2} s^{−1} and (b) ν = 10^{12} cm^{2} s^{−1}. We note that all the modes selected here are stable modes (ℑ[ω]< 0). Different colors represent different types of inertial modes. The horizontal black dashed line shows the length of the SDO/HMI observational record (T_{obs} ≈ 12 yr as of today). In both cases, rotation is uniform and the stratification is adiabatic. The lifetimes of the convective modes and highlatitude modes are very sensitive to the radial and latitudinal entropy gradients, a point that is discussed in Sects. 5 and 6.2. 

In the text 
Fig. 16. (a) Meridional eigenfunctions of radial vorticity ζ_{r} of the n = 0 equatorial Rossby mode at m = 16 for different values of viscous diffusivities ν. Upper and lower panels: normalized real and imaginary eigenfunctions, respectively. (b) Radial eigenfunctions of ζ_{r} at the equator (normalized by their maximum amplitudes). Different colors represent different values of diffusivities. In all cases, rotation is uniform and the stratification is adiabatic. 

In the text 
Fig. 17. (a) Dispersion relations of the northsouth ζ_{z}symmetric columnar convective modes with different background superadiabaticity values δ. Different colors represent different values of superadiabaticity. Circles and diamonds denote the stable (ℑ[ω]< 0) and unstable (ℑ[ω]> 0) modes, respectively. (b) Eigenfrequencies in the complex plane. Each circle (diamond) represent a mode with azimuthal order m, which is labeled with small integers from m = 1 to 16. In all cases, rotation is uniform and the diffusivities are set ν = κ = 10^{12} cm^{2} s^{−1}. 

In the text 
Fig. 18. Radial velocity v_{r} (upper plots) and entropy perturbation s_{1} (lower plots) of the columnar convective modes along the rotational axis displayed in the equatorial plane for subadiabatic and superadiabatic background. Panels a and b: cases with subadiabatic background δ = −2 × 10^{−6} for m = 3 and m = 8, respectively. Panels c and d: are the same plots for a superadiabatic background δ = 2 × 10^{−6}. The eigenfunctions are normalized such that the maximum radial velocity is 10 m s^{−1} at the equator. In all cases, rotation is uniform and the diffusivities are set ν = κ = 10^{12} cm^{2} s^{−1}. 

In the text 
Fig. 19. Transport properties of thermal energy and angular momentum by the northsouth ζ_{z}symmetric columnar convective modes for m = 16. (a) Correlation between radial velocity velocity and temperature perturbation ⟨v_{r}T_{1}⟩, (b) Reynolds stress between radial and longitudinal velocities ⟨v_{r}v_{ϕ}⟩, and (c) Reynolds stress between latitudinal and longitudinal velocities ⟨v_{θ}v_{ϕ}⟩. The background is weakly superadiabatic (δ = 2 × 10^{16}), rotation is uniform, and moderate diffusivities are used (ν = κ = 10^{12} cm^{2} s^{−1}). The eigenfunctions are normalized such that the maximum radial velocity is 10 m s^{−1} at the equator. 

In the text 
Fig. 20. Solar differential rotation profile used in this study. (a) Differential rotation Ω(r, θ) in a meridional plane, deduced from the global helioseismology (Larson & Schou 2018). (b) Latitudinal profiles of differential rotation at different depths. Horizontal dashed lines indicate the theoreticallyexpected phase speed of the sectoral (l = m) classical Rossby modes for selected azimuthal orders m = 2, 3, 4, 8, 16. The observing frame is chosen to be the Carrington frame rotating at Ω_{0}/2π = 456.0 nHz. 

In the text 
Fig. 21. Eigenfrequencies of inertial modes under the solar differential rotation in the complex plane for (a) m = 2 and (b) m = 16. The diffusivity is set to ν = 10^{12} cm^{2} s^{−1} and the background is assumed to be adiabatic, δ = 0. The shaded areas indicate the frequency range associated with the surface differential rotation, namely, m(Ω_{pole} − Ω_{0}) < ℜ[ω]< m(Ω_{eq} − Ω_{0}). 

In the text 
Fig. 22. (a) Dispersion relations of the equatorial Rossby modes for the cases with solar differential rotation. Red and blue curves represent the modes with no radial nodes (n = 0) and one radial node (n = 1), respectively. Solid and dashed lines denote the cases with weak diffusion (ν = 10^{11} cm^{2} s^{−1}) and strong diffusion (ν = 10^{12} cm^{2} s^{−1}), respectively. For comparison, the observed Rossby mode frequencies reported in Löptien et al. (2018), Liang et al. (2019), and Proxauf et al. (2020) are plotted by white hexagons and squares. All the presented frequencies are measured in the Carrington frame rotating at Ω_{0}/2π = 456.0 nHz. (b) Mode linewidths versus mode frequencies. Each point represents a mode with azimuthal order m, which is labeled with small integers from m = 1 to 16. Overplotted (open circles connected by black lines) are the mode linewidths and frequencies measured by Proxauf et al. (2020) for m = 3 to 15. In both panels, the green dots connected by line segments refer to a theoretical model for the simplified case of uniform rotation, ω/Ω_{0} = −2/(m + 1)−iE_{k}m(m + 1), where E_{k} = 4 × 10^{−4} is the Ekman number at the solar surface (Fournier et al. 2022). 

In the text 
Fig. 23. Eigenfunctions of the equatorial Rossby modes with no radial nodes (n = 0). (a) Real (upper) and imaginary (lower) eigenfunctions of three components of velocity shown in a meridional plane for m = 5. The eigenfunctions are normalized such that the maximum latitudinal velocity is 2 m s^{−1} at the surface. The solid black line indicates the location of the critical latitudes where the phase speed of a Rossby mode matches to the differential rotation sped. (b) Horizontal eigenfunctions of latitudinal velocity v_{θ} (upper) and radial vorticity ζ_{r} (lower) at the surface r = 0.985 R_{⊙} for m = 5. The horizontal black dashed lines indicate the location of the critical latitudes at the surface. (c) and (d) are the counterparts of the panels a and b for m = 12. 

In the text 
Fig. 24. Eigenfunctions of the equatorial Rossby modes with one radial node (n = 1). Same figure notation as in Fig. 23 is used. 

In the text 
Fig. 25. Radial vorticity ζ_{r} eigenfunctions at the surface (r = 0.985 R_{⊙}) for m = 8. Left and right panels: cases for the equatorial Rossby modes with no radial nodes (n = 0) and with one radial node (n = 1). Upper and lower panels: cases with weak diffusion (ν = 10^{11} cm^{2} s^{−1}) and strong diffusion (ν = 10^{12} cm^{2} s^{−1}). Black solid and dashed lines represent real and imaginary eigenfunctions, respectively. The real part of the eigenfunctions are defined to be zero at the equator. The vertical red lines denote the location of critical latitudes where the phase speed of a Rossby mode is equal to the differential rotation velocity. 

In the text 
Fig. 26. Reynolds stress components (a), (d) ρ_{0}⟨v_{r}v_{ϕ}⟩ and (b), (e) ρ_{0}⟨v_{θ}v_{ϕ}⟩ associated with the equatorial Rossby modes for m = 8. The units are g cm^{−1} s^{−2}. Black solid lines denote the location of critical latitudes at each height. The eigenfunctions are normalized such that the maximum horizontal velocity at the top boundary is 2 m s^{−1}. Panels c and f: horizontal Reynolds stress averaged over radius where the overbar denotes the radial average. Different colors represent different azimuthal orders. Upper and lower panels: cases for n = 0 modes and n = 1 Rossby modes, respectively. 

In the text 
Fig. 27. Radially averaged latitudinal angular momentum fluxes. Solid, dotdashed, and dashed lines represent the angular momentum fluxes associated with the Reynolds stress of the equatorial Rossby modes F_{RS, θ}, advection by meridional circulation F_{MC, θ}, and diffusion by turbulent viscosity F_{VD, θ}, defined by the Eqs. (23)–(25). Red and blue lines represent the modes with no radial nodes (n = 0) and with one radial node (n = 1), respectively. Here, eigenfunctions are normalized such that the maximum horizontal velocity amplitude at the surface is 2 m s^{−1}. 

In the text 
Fig. 28. (a) Growth rate versus frequency of the ζ_{z}antisymmetric highlatitude mode with m = 1, for different values of the latitudinal entropy difference Δ_{θ}s=s_{0, pole} − s_{0, eq}. The red star is for the case of a realistic latitudinal entropy gradient that depends on position, estimated according to Eq. (27). Realistic solar differential rotation is included. The background stratification is adiabatic (δ = 0) and the diffusivities are set to ν = κ = 10^{12} cm^{2} s^{−1}. (b)–(f) Eigenfunctions of the longitudinal velocity v_{ϕ} at the surface (r = 0.985 R_{⊙}) and at the central meridian for some selected Δs_{θ}. The eigenfunctions are normalized so that the maximum flow amplitudes at the surface is v_{ϕ} = 10 m s^{−1}. 

In the text 
Fig. 29. Eigenfunctions of v_{r}, v_{θ}, v_{ϕ}, and s_{1} of the m = 1 northsouth ζ_{z}antisymmetric highlatitude inertial mode with the solar differential rotation and the corresponding latitudinal entropy gradient (Eq. (27)) included. The background stratification is adiabatic (δ = 0) and the diffusivities are set to ν = κ = 10^{12} cm^{2} s^{−1}. The eigenfunctions are normalized such that the maximum v_{ϕ} is 10 m s^{−1} at the surface. 

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.