EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 662 (June 2022) A&A, 662 (2022) A16 Full HTML
Open Access
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/0004-6361/202243164
Published online 02 June 2022

© Y. Bekki et al. 2022

Licence Creative CommonsOpen 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 quasi-toroidal 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 well-defined 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 follow-up studies that confirm these observations (e.g., Liang et al. 2019; Hanasoge & Mandal 2019; Proxauf et al. 2020; Hanson et al. 2020). Using a one-dimensional β-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 mid-latitudes that are localized near their critical latitudes. Several tens of critical-latitude modes have been identified in the range m ≤ 10. Another family of inertial modes introduced by the Sun’s differential rotation are the high-latitude modes (Gizon et al. 2021). The highest amplitude mode (∼10−20 m s−1 above 50° latitude) is the m = 1 mode with north-south antisymmetric longitudinal velocity vϕ with respect to the equator. This m = 1 mode was identified by Gizon et al. (2021) using linear calculations in two-dimensional model, which are further discussed in the present paper and Fournier et al. (2022). It corresponds to the spiral-like velocity feature reported at high latitudes by Hathaway et al. (2013), although it was reported as giant-cell convection in that work.

The equatorial-Rossby and high-latitude modes involve mostly toroidal motions with a radial velocity that is small compared to the horizontal velocity components. Non-toroidal inertial modes have also been theoretically studied, mainly in the context of incompressible fluids. These modes tend to be localized onto so-called 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 solar-like 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 one-dimensional 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 high-latitude 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 non-adiabatic stratification on these modes. We note that the critical-latitude 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 1012 cm2 s−1, the equatorial Rossby modes with no radial node (n = 0) strongly depart from the expected rm 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 high-latitude 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 low-frequency modes are discussed in Sect. 3 for the inviscid, adiabatically stratified, and uniformly rotating case. Then, the effects of turbulent diffusion and a non-adiabatically 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 fully-compressible hydrodynamic equations in a spherical coordinate (r, θ, ϕ).

2.1. Linearized equations

The linearized equations of motion, continuity, and energy conservation are:

(1)

(2)

(3)

where, v = (vr, vθ, vϕ) is the 1st-order 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 p0, ρ0, T0, g, and Hp 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 adiabatically-stratified 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 (Christensen-Dalsgaard et al. 1996). The variables with subscript 1, p1, ρ1, and s1, represent the first-order 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:

(4)

where γ = 5/3 is the specific heat ratio and cv 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 double-logarithmic 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 ∂s0/∂θ 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

(5)

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, p1, and s1 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

(6)

(7)

(8)

(9)

(10)

where Cs = (γp0/ρ0)1/2 is the sound speed and cp = γcv is the constant specific heat at constant pressure. Here, the longitudinal velocity vϕ, density perturbation ρ1, and entropy perturbation s1 are out of phase with the meridional components of velocity (vr and vθ) in the inviscid limit (v = κ = 0).

Equations (6)–(10) can be combined into an eigenvalue problem:

(11)

where

(12)

and M is the linear differential operator represented by the right-hand 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 rmin = 0.71 R⊙ to rmax = 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 second-order (in both the radial and latitudinal directions) PDEs and one first-order 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 stress-free conditions for the velocity and assume there is no entropy flux (∝κ∂s1/∂r) across the boundary:

(13)

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 non-axisymmetric cases (m ≠ 0), at the poles we impose

(14)

to make the quantities single valued. For the axisymmetric case (m = 0), at both poles we assume instead

(15)

2.4. Numerical scheme

We numerically solve the above eigenvalue problem using a finite differencing method in the meridional plane. We use a spatially-uniform grids. The grids for vϕ, ρ1, and s1 are staggered grids by half a grid point in radius for vr 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 second-order accurate scheme. By converting the two dimensional grid (Nr, Nθ) into one dimensional array with the size Nr, Nθ for all variables, V is defined as a one dimensional vector with size ∼5NrNθ. Once the boundary conditions are properly set, M can be constructed as a two-dimensional complex matrix with the size approximately (5NrNθ × 5NrNθ). 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 right-hand side of the Eqs. (6)–(10). In most of the calculations, we use the grid resolution of (Nr, Nθ)=(16, 72). We have also carried out higher-resolution calculations with (Nr, 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 large-scale modes which have been observed on the Sun, the results are converged with (Nr, Nθ)=(16, 72).

thumbnail Fig. 1.

Layout of the staggered grid used to solve the eigenvalue equation. The grid locations where vϕ, ρ1, and s1 are defined are denoted by red circles. The blue and green circles represent the grid locations of vr 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 (vr, vθ, vϕ, ρ1, s1) 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 retrograde-propagating modes and ℑ[ω] > 0 corresponds to exponentially growing modes.

2.5. Example spectrum for uniform rotation

For each m, there are 5NrNθ 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 × 1012 cm2 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.

thumbnail 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 (ν = κ = 1011 cm2 s−1). Shown are the results at m = 8. Upper panels: complex eigenfrequencies ω in the co-rotating frame. Panel a: real frequencies in the range ±10 mHz showing the acoustic modes (p modes). Panel b: zoom-in focusing on the inertial range |ℜ[ω]| < 2Ω0. Panel c: zoom-in focusing on the convectively-unstable modes (ℑ[ω]> 0). Lower panels: example eigenfunctions of pressure (acoustic), non-toroidal 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 low-frequency 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.

thumbnail 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 ∂s0/∂θ = 0). We present the dispersion relations and eigenfunctions of various types of global-scale 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, non-adiabatic stratification, and differential rotation are compared to these reference results.

In the inviscid case with uniform rotation, M is self adjoint, thus the physically-meaningful 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 vr and vθ have the same complex phase on each meridional plane, and those of vϕ, ρ1 are 90° out of phase with respect to vr and vθ. In presenting the results in this section, we choose a meridional plane where vr and vθ are real.

In the following sections, we conduct a mode-by-mode 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 north-south symmetries, and the high-latitude modes with both north-south symmetries. Fundamental properties of these modes are summarized in Table 1. Their dispersion relations are presented in Table 2.

Table 1.

Summary of the properties of the modes of the models discussed in this paper.

Table 2.

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:

(16)

(17)

where lmax = 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 (|a0|> |aj| for all j > 0); the l = m + 1 component of vϕ is dominant (|b1|> |bj| 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, uniformly-rotating 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 theoretically-expected dispersion relation, ω = −2Ω0/(m + 1), respectively. We find that the differences in the normalized frequencies are less than 10−2 at all m.

thumbnail 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 vr and vθ are purely real and vϕ, p1, and ζr are purely imaginary. The black dashed line denotes the meridional plane at ϕ = −π/2m where vϕ, p1 and ζr are real. (c) Meridional cuts of the m = 2 eigenfunctions for the velocity v(r, θ)exp[i(mϕ − ωt)], the pressure p1(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, 105 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 vr is about 103 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 vr. The pressure perturbation, p1, 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θ ∝ rm. It is seen that computed eigenfunctions exhibit the rm dependence that agree with the analytical solutions. We also confirm the same rm 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 stress-free boundary condition, ∂(vθ/r)/∂r = 0, at the top and bottom boundaries, which conflicts with the rm 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θ ∝ sinm − 1θ.

thumbnail 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θ ∝ rm. The eigenfunctions are normalized to unity at the surface r = rmax. (b) Latitudinal structure of the eigenfunction of vθ at the surface. Dashed lines are the theoretical solution in the form of Legendre polynomials, vθ ∝ sinm − 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 equatorially-trapped axisymmetric inertial mode. We show in Sect. 3.2.2 that this mode is connected to a prograde-propagating 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.

thumbnail 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 vr 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 three-dimensional.

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.

thumbnail 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 = rmax. (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 mode-by-mode analysis for the columnar convective modes (thermal Rossby waves) with both hemispheric symmetries. Here, we define the north-south symmetry based on the eigenfunction of z-vorticity ζz. The “banana cell” convection pattern can be essentially regarded as the north-south symmetric part of these convective modes. We will also show that the north-south ζz-antisymmetric modes are essentially mixed with the n = 1 equatorial modes.

3.2.1. North-south ζz-symmetric modes

North-south ζ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 vr 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 north-south ζz-symmetric columnar convective modes. For comparison, we overplot in black dashed line the dispersion relation derived from the one-dimensional 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 non-dispersive 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 one-dimensional 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).

thumbnail Fig. 8.

Dispersion relation and eigenfunctions of the north-south ζz-symmetric columnar convective modes in the case of uniform rotation, no viscosity, and adiabatic stratification. (a) Dispersion relation of the north-south ζz-symmetric columnar convective modes in red points. For comparison, dispersion relation analytically derived using one-dimensional 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 p1(r, θ)exp[i(mϕ − ωt)], and the z-vorticity ζr(r, θ)exp[i(mϕ − ωt)]. The solutions are shown in the meridional plane at ϕ = 0 and t = 0 where vr and vθ are purely real and vϕ, p1 and ζz are purely imaginary. The units of velocity, pressure, and vorticity are m s−1, 105 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 north-south ζz-symmetric columnar convective modes. The flow structure is dominantly characterized by the longitudinal velocity shear outside the tangential cylinder, leading to a strong z-vorticity (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 z-vortex 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 p1 is generally positive (negative) where z-vorticity ζ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. North-south ζz-antisymmetric modes

North-south ζ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.

thumbnail Fig. 9.

Dispersion relation and eigenfunctions of the north-south ζ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 north-south ζz-symmetric columnar convective modes is shown in black dashed line for comparison.

Figures 9c and d show the example eigenfunctions of the north-south ζz-antisymmetric columnar convective modes. The flow structure is dominantly characterized by z-vortex 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 non-propagating 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 north-south ζ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 so-called Yanai waves which are mixed modes between retrograde-propagating Rossby modes and prograde-propagating inertial-gravity modes (Matsuno 1966; Vallis 2006).

thumbnail Fig. 10.

Dispersion relation of the “mixed modes” between the n = 1 equatorial Rossby modes (red) and the north-south ζ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 dot-dashed lines represent the dispersion relation of the n = 0 equatorial Rossby modes and north-south ζ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 high-latitude inertial modes with both hemispheric symmetries.

3.3.1. North-south ζ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 ϖ = rmin, namely, it is the cylinder aligned with the rotation axis which touches the radiative interior at the equator. North-south ζz-symmetric high-latitude modes can be selected by applying the following criteria: The kinetic energy is predominantly inside the tangential cylinder, namely, Ein/ECZ > 0.5, where Ein and ECZ are the volume-integrated 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 z-nodes of vθ is zero at ϖ = 0.5 R⊙.

Figure 11a shows the dispersion relation of the north-south ζz-symmetric high-latitude modes. We find the high-latitude modes are much more dispersive than the columnar convective modes at low m. The dispersion relation is found to be roughly approximated by the non-sectoral 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.

thumbnail Fig. 11.

Dispersion relation and eigenfunctions of the high-latitude modes with north-south 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 high-latitude modes. The fluid motion is predominantly characterized by z-vortices inside the tangential cylinder in both hemispheres, as schematically illustrated in Fig. 11b. The power of ζz peaks at the tangential cylinder ϖ = rmin. We note that the longitudinal velocity vϕ extends slightly outside the cylinder. Again, ℑ[p1]ℑ[ζz]< 0 follows from the mode being in geostrophic balance.

3.3.2. North-south ζz-antisymmetric modes

North-south ζz-antisymmetric high-latitude 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 z-nodes is zero for vθ at ϖ = 0.5 R⊙.

The example modes are presented in Fig. 12. The eigenfunctions show very similar properties of the high-latitude modes discussed in Fig. 11 except for the north-south 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 high-latitude modes are found to be similar to that of l = m + 2 Rossby modes, as shown in Fig. 12a.

thumbnail Fig. 12.

Dispersion relation and eigenfunctions of north-south ζz-antisymmetric high-latitude 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

(18)

On the other hand, typical diffusive timescale can be estimated as:

(19)

where lm denotes the typical length scale of the Rossby mode. Figure 13 compares PRo and τdiff as functions of m. Two representative values of turbulent diffusitivies in the solar convection zone ν = 1012 and 1013 cm2 s−1 are shown (e.g., Ossendrijver 2003). When PRo ≪ τdiff, viscous diffusion is almost negligible. However, if PRo ≳ τdiff, diffusion can have a dominant effect on the Rossby modes. For a given turbulent diffusivity ν, the critical azimuthal order, mcrit, can be defined as:

(20)

thumbnail Fig. 13.

Comparison between the oscillation periods of Rossby modes PRo and the diffusive timescale, τdiff, for two representative values of turbulent diffusivities of ν = 1012 and 1013 cm2 s−1. The horizontal black dashed line represents the length of the SDO/HMI observational record of Tobs ≈ 12 years.

The Rossby modes are dominated by diffusive effects for m > mcrit. 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 uniformly-rotating adiabatic fluid with varying diffusivities; ν = 109, 1010, 1011, 1012, and 1013 cm2 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 high-latitude 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 e-folding times of these modes |ℑ[ω]|−1 are shown in Fig. 15 for the two representative values of turbulent viscosity ν = 1011 and 1012 cm2 s−1.

thumbnail Fig. 14.

Eigenfrequency spectra of the low-frequency 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.

thumbnail Fig. 15.

e-folding lifetimes of various low-frequency modes for a viscous diffusivity (a) ν = 1011 cm2 s−1 and (b) ν = 1012 cm2 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 (Tobs ≈ 12 yr as of today). In both cases, rotation is uniform and the stratification is adiabatic. The lifetimes of the convective modes and high-latitude 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 well-known rm 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.).

thumbnail 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 non-adiabatic stratification

In this section, the effects of non-adiabatic 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 (ν = κ = 1012 cm2 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 (north-south ζ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.

thumbnail Fig. 17.

(a) Dispersion relations of the north-south ζ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 ν = κ = 1012 cm2 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 vr and s1 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 s1 is always ahead of the phase with positive vr in longitude, leading to a negative correlation between ℜ[vr] and ℑ[s1]. In a physical sense, this means that, in this case, the buoyancy force acts as an additional restoring force for prograde-propagating columnar convective modes. In other words, these modes share a property of prograde-propagating 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 vr and s1 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 s1 is behind the phase with positive vr in longitude, leading to a positive correlation between ℜ[vr] and ℑ[s1]. 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 vr and s1 now coincide and they both have the same sign at each phase, leading to ⟨vrs1⟩> 0.

thumbnail Fig. 18.

Radial velocity vr (upper plots) and entropy perturbation s1 (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 ν = κ = 1012 cm2 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 ⟨vrT1⟩ in Fig. 19a manifests that the enthalpy flux is transported upward. The Reynolds stress components ⟨vrvϕ⟩ 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.).

thumbnail Fig. 19.

Transport properties of thermal energy and angular momentum by the north-south ζz-symmetric columnar convective modes for m = 16. (a) Correlation between radial velocity velocity and temperature perturbation ⟨vrT1⟩, (b) Reynolds stress between radial and longitudinal velocities ⟨vrvϕ⟩, and (c) Reynolds stress between latitudinal and longitudinal velocities ⟨vθvϕ⟩. The background is weakly superadiabatic (δ = 2 × 1016), rotation is uniform, and moderate diffusivities are used (ν = κ = 1012 cm2 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 = rmin and rmax, 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

(21)

thumbnail 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 theoretically-expected 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 global-scale 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 retrograde-propagating inertial modes are affected by critical latitudes at higher m (see Fig. 21b).

thumbnail 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 ν = 1012 cm2 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 ν = 1011 and 1012 cm2 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 ∂s0/∂θ 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.

thumbnail 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 (ν = 1011 cm2 s−1) and strong diffusion (ν = 1012 cm2 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)−iEkm(m + 1), where Ek = 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 ν = 1012 cm2 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 uniformly-rotating case, strong radial and longitudinal flows are driven around the critical latitudes, which leads to strong concentrations of z-vorticity 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 chevron-like inclination towards the equator. However, ζr has prominent power peaks around the critical layers.

thumbnail 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.

thumbnail 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.

thumbnail 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 (ν = 1011 cm2 s−1) and strong diffusion (ν = 1012 cm2 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 ⟨vrvθ⟩ and ⟨vθvθ⟩ for n = 0 modes at m = 8. The Reynolds stresses become substantially non-zero near the viscous critical layers. It is striking that even n = 0 mode, which in the case of uniform rotation is toroidal and non-convective, 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.

thumbnail Fig. 26.

Reynolds stress components (a), (d) ρ0⟨vrvϕ⟩ 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 so-called gyroscopic pumping equation (e.g., Elliott et al. 2000; Miesch & Hindman 2011)

(22)

where FRS, FMC, and FVD are the angular momentum fluxes transported by the Reynolds stress, meridional circulation, and turbulent viscous diffusion, respectively. They are defined by

(23)

(24)

(25)

where vm 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 FMC, θ (black dot-dashed line), we use the observational meridional circulation data obtained by Gizon et al. (2020b). For FVD, θ (black dashed line), we assume the spatially-uniform viscosity of ν = 1012 cm2 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 FRS, θ 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 FMC, θ + FVD, θ. 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.

thumbnail Fig. 27.

Radially averaged latitudinal angular momentum fluxes. Solid, dot-dashed, and dashed lines represent the angular momentum fluxes associated with the Reynolds stress of the equatorial Rossby modes FRS, θ, advection by meridional circulation FMC, θ, and diffusion by turbulent viscosity FVD, θ, 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 high-latitude inertial modes

We assume the solar differential rotation is in thermal wind balance; where the deviation from the Taylor-Proudman’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 high-latitude 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 high-latitude inertial modes by varying the amplitude of the imposed latitudinal entropy gradient. Of particular interest is this effect on the m = 1 mode with north-south antisymmetric ζz. Here, we assume the latitudinal dependence of the background entropy profile is

(26)

where Δθs = s0,eq − s0,pole(< 0) represents the entropy difference between the cooler equator and the hotter poles. For simplicity, the radial dependence is ignored (s0 is uniform in radius and thus convectively neutral). We use moderately viscous and thermal diffusivities ν = κ = 1012 cm2 s−1.

Figure 28a shows the eigenfrequencies of the m = 1 north-south ζz-antisymmetric high-latitude 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 high-latitude 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.

thumbnail Fig. 28.

(a) Growth rate versus frequency of the ζz-antisymmetric high-latitude mode with m = 1, for different values of the latitudinal entropy difference |Δθs|=s0, pole − s0, 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 ν = κ = 1012 cm2 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 helioseismically-constrained differential rotation profile using:

(27)

With this realistic baroclinicity included in our model (Eq. (10)), we find that the m = 1 high-latitude mode is self-excited: the growth rate is ℑ[ω]/2π = 14.1 nHz, which translates into the growing timescale of 4.3 months. This may explain why the high-latitude 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 high-latitude 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 z-vortical motion and is quasi-toroidal (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 ring-diagram measurements presented by Gizon et al. (2021). This means that the observations may miss a fraction of the mode power.

thumbnail Fig. 29.

Eigenfunctions of vr, vθ, vϕ, and s1 of the m = 1 north-south ζz-antisymmetric high-latitude 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 ν = κ = 1012 cm2 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 uniformly-rotating spherical shell, Gilman (1975) found convectively-unstable modes near the poles with a spiral structure (see his Fig. 17). As with the high-latitude 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 high-latitude modes, both with different north-south 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 well-known rm 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 non-zero 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 high-latitude 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 non-adiabaticity 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 (rmin, rmax)=(0.71 R⊙, 0.985 R⊙) and thus both the tachocline and the near-surface 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 three-dimensional numerical simulations of rotating convection in the strongly non-linear regime. The aim is to have a physical understanding of all the modes in the low-frequency spectrum and thus a reliable identification of the observed modes, including the critical-latitude 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 Max-Planck 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 long-term scholarship program for degree-seeking graduate students abroad from the Japan Student Services Organization (JASSO).

References

  1. 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]
  2. Balbus, S. A., Bonart, J., Latter, H. N., & Weiss, N. O. 2009, MNRAS, 400, 176 [NASA ADS] [CrossRef] [Google Scholar]
  3. Baruteau, C., & Rieutord, M. 2013, J. Fluid Mech., 719, 47 [Google Scholar]
  4. Bekki, Y., & Yokoyama, T. 2017, ApJ, 835, 9 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bessolaz, N., & Brun, A. S. 2011, ApJ, 728, 115 [NASA ADS] [CrossRef] [Google Scholar]
  7. Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
  8. Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
  9. Busse, F. H. 1970, J. Fluid Mech., 44, 441 [NASA ADS] [CrossRef] [Google Scholar]
  10. Busse, F. H. 2002, Phys. Fluids, 14, 1301 [NASA ADS] [CrossRef] [Google Scholar]
  11. Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [Google Scholar]
  12. Elliott, J. R., Miesch, M. S., & Toomre, J. 2000, ApJ, 533, 546 [NASA ADS] [CrossRef] [Google Scholar]
  13. Evonuk, M. 2008, ApJ, 673, 1154 [NASA ADS] [CrossRef] [Google Scholar]
  14. Evonuk, M., & Samuel, H. 2012, Earth Planet. Sci. Lett., 317, 1 [CrossRef] [Google Scholar]
  15. Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
  16. Featherstone, N. A., & Hindman, B. W. 2016, ApJ, 830, L15 [Google Scholar]
  17. Fournier, D., Gizon, L., & Hyest, L. 2022, A&A, in press http://dx.doi.org/10.1051/0004-6361/202243473 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156 [Google Scholar]
  19. Gilman, P. A. 1975, J. Atm. Sci., 32, 1331 [NASA ADS] [CrossRef] [Google Scholar]
  20. 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]
  21. Gilman, P. A. 1987, ApJ, 318, 904 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gilman, P., & Dikpati, M. 2014, ApJ, 787, 60 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gizon, L., Fournier, D., & Albekioni, M. 2020a, A&A, 642, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020b, Science, 368, 1469 [Google Scholar]
  25. Gizon, L., Cameron, R. H., Bekki, Y., et al. 2021, A&A, 652, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Glatzmaier, G. A., & Gilman, P. A. 1981, ApJS, 45, 381 [NASA ADS] [CrossRef] [Google Scholar]
  27. Glatzmaier, G., Evonuk, M., & Rogers, T. 2009, Geophys. Astrophys. Fluid Dyn., 103, 31 [NASA ADS] [CrossRef] [Google Scholar]
  28. Greenspan, H., Batchelor, C., Ablowitz, M., et al. 1968, The Theory of Rotating Fluids, Cambridge Monographs on Mechanics (Cambridge: Cambridge University Press) [Google Scholar]
  29. Guenel, M., Baruteau, C., Mathis, S., & Rieutord, M. 2016, A&A, 589, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Guenther, D. B., & Gilman, P. A. 1985, ApJ, 295, 195 [Google Scholar]
  31. Hanasoge, S., & Mandal, K. 2019, ApJ, 871, L32 [Google Scholar]
  32. Hanson, C. S., Gizon, L., & Liang, Z.-C. 2020, A&A, 635, A109 [EDP Sciences] [Google Scholar]
  33. Hathaway, D. H., & Upton, L. A. 2021, ApJ, 908, 160 [NASA ADS] [CrossRef] [Google Scholar]
  34. Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  35. Hindman, B. W., Featherstone, N. A., & Julien, K. 2020, ApJ, 898, 120 [NASA ADS] [CrossRef] [Google Scholar]
  36. Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
  37. Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 803, 42 [NASA ADS] [CrossRef] [Google Scholar]
  38. Ingersoll, A. P., & Pollard, D. 1982, Icarus, 52, 62 [NASA ADS] [CrossRef] [Google Scholar]
  39. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011, Astron. Nachr., 332, 883 [CrossRef] [Google Scholar]
  40. Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23 [CrossRef] [Google Scholar]
  41. Käpylä, P. J., Viviani, M., Käpylä, M. J., Brandenburg, A., & Spada, F. 2019, Geophys. Astrophys. Fluid Dyn., 113, 149 [CrossRef] [Google Scholar]
  42. Karak, B. B., Miesch, M., & Bekki, Y. 2018, Phys. Fluids, 30, 046602 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kitchatinov, L. L. 2013, Astron. Lett., 39, 561 [NASA ADS] [CrossRef] [Google Scholar]
  44. Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
  45. Larson, T. P., & Schou, J. 2018, Sol. Phys., 293, 29 [Google Scholar]
  46. Liang, Z.-C., Gizon, L., Birch, A. C., & Duvall, T. L. 2019, A&A, 626, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat. Astron., 2, 568 [Google Scholar]
  48. Lorenzani, S., & Tilgner, A. 2001, J. Fluid Mech., 447, 111 [NASA ADS] [CrossRef] [Google Scholar]
  49. Maas, L. R. M., & Lam, F.-P. A. 1995, J. Fluid Mech., 300, 1 [NASA ADS] [CrossRef] [Google Scholar]
  50. Matilsky, L. I., Hindman, B. W., & Toomre, J. 2020, ApJ, 898, 111 [NASA ADS] [CrossRef] [Google Scholar]
  51. Matsuno, T. 1966, J. Meteorol. Soc. Japan. Ser. II, 44, 25 [NASA ADS] [CrossRef] [Google Scholar]
  52. Miesch, M. S., & Hindman, B. W. 2011, ApJ, 743, 79 [NASA ADS] [CrossRef] [Google Scholar]
  53. Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [CrossRef] [Google Scholar]
  54. Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618 [Google Scholar]
  55. Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
  56. Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
  57. Papaloizou, J., & Pringle, J. E. 1978, MNRAS, 182, 423 [NASA ADS] [CrossRef] [Google Scholar]
  58. Proxauf, B., Gizon, L., Löptien, B., et al. 2020, A&A, 634, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
  60. Rieutord, M., & Valdettaro, L. 1997, J. Fluid Mech., 341, 77 [Google Scholar]
  61. Rieutord, M., & Valdettaro, L. 2018, J. Fluid Mech., 844, 597 [Google Scholar]
  62. Rieutord, M., Georgeot, B., & Valdettaro, L. 2001, J. Fluid Mech., 435, 103 [Google Scholar]
  63. Rüdiger, G. 1989, Differential Rotation and Stellar Convection (Berlin: Akademie Verlag) [Google Scholar]
  64. Saio, H. 1982, ApJ, 256, 717 [Google Scholar]
  65. Sibgatullin, I. N., & Ermanyuk, E. V. 2019, J. Appl. Mech. Tech. Phys., 60, 284 [NASA ADS] [CrossRef] [Google Scholar]
  66. Smeyers, P., Craeynest, D., & Martens, L. 1981, Ap&SS, 78, 483 [NASA ADS] [CrossRef] [Google Scholar]
  67. Spruit, H. C., & Knobloch, E. 1984, A&A, 132, 89 [NASA ADS] [Google Scholar]
  68. Tilgner, A. 2007, in Treatise on Geophysics, ed. G. Schubert (Amsterdam: Elsevier), 8.07, 207 [NASA ADS] [Google Scholar]
  69. Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  70. Verhoeven, J., & Stellmach, S. 2014, Icarus, 237, 143 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Summary of the properties of the modes of the models discussed in this paper.

In the text
Table 2.

Dispersion relations of the modes of the model with uniform rotation (Ω = Ω0), ν = κ = 0, and δ = 0.

In the text

All Figures

thumbnail Fig. 1.

Layout of the staggered grid used to solve the eigenvalue equation. The grid locations where vϕ, ρ1, and s1 are defined are denoted by red circles. The blue and green circles represent the grid locations of vr and vθ, respectively. The grid resolution is reduced for a visualization purpose.
In the text
thumbnail 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 (ν = κ = 1011 cm2 s−1). Shown are the results at m = 8. Upper panels: complex eigenfrequencies ω in the co-rotating frame. Panel a: real frequencies in the range ±10 mHz showing the acoustic modes (p modes). Panel b: zoom-in focusing on the inertial range |ℜ[ω]| < 2Ω0. Panel c: zoom-in focusing on the convectively-unstable modes (ℑ[ω]> 0). Lower panels: example eigenfunctions of pressure (acoustic), non-toroidal 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
thumbnail 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
thumbnail 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 vr and vθ are purely real and vϕ, p1, and ζr are purely imaginary. The black dashed line denotes the meridional plane at ϕ = −π/2m where vϕ, p1 and ζr are real. (c) Meridional cuts of the m = 2 eigenfunctions for the velocity v(r, θ)exp[i(mϕ − ωt)], the pressure p1(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, 105 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
thumbnail 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θ ∝ rm. The eigenfunctions are normalized to unity at the surface r = rmax. (b) Latitudinal structure of the eigenfunction of vθ at the surface. Dashed lines are the theoretical solution in the form of Legendre polynomials, vθ ∝ sinm − 1θ. All the eigenfunctions are normalized at the equator.
In the text
thumbnail 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
thumbnail 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 = rmax. (b) Latitudinal structure of the eigenfunction of vθ at the surface normalized at the equator.
In the text
thumbnail Fig. 8.

Dispersion relation and eigenfunctions of the north-south ζz-symmetric columnar convective modes in the case of uniform rotation, no viscosity, and adiabatic stratification. (a) Dispersion relation of the north-south ζz-symmetric columnar convective modes in red points. For comparison, dispersion relation analytically derived using one-dimensional 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 p1(r, θ)exp[i(mϕ − ωt)], and the z-vorticity ζr(r, θ)exp[i(mϕ − ωt)]. The solutions are shown in the meridional plane at ϕ = 0 and t = 0 where vr and vθ are purely real and vϕ, p1 and ζz are purely imaginary. The units of velocity, pressure, and vorticity are m s−1, 105 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
thumbnail Fig. 9.

Dispersion relation and eigenfunctions of the north-south ζ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 north-south ζz-symmetric columnar convective modes is shown in black dashed line for comparison.
In the text
thumbnail Fig. 10.

Dispersion relation of the “mixed modes” between the n = 1 equatorial Rossby modes (red) and the north-south ζ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 dot-dashed lines represent the dispersion relation of the n = 0 equatorial Rossby modes and north-south ζz-symmetric columnar convective modes.
In the text
thumbnail Fig. 11.

Dispersion relation and eigenfunctions of the high-latitude modes with north-south 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
thumbnail Fig. 12.

Dispersion relation and eigenfunctions of north-south ζz-antisymmetric high-latitude 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
thumbnail Fig. 13.

Comparison between the oscillation periods of Rossby modes PRo and the diffusive timescale, τdiff, for two representative values of turbulent diffusivities of ν = 1012 and 1013 cm2 s−1. The horizontal black dashed line represents the length of the SDO/HMI observational record of Tobs ≈ 12 years.
In the text
thumbnail Fig. 14.

Eigenfrequency spectra of the low-frequency 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
thumbnail Fig. 15.

e-folding lifetimes of various low-frequency modes for a viscous diffusivity (a) ν = 1011 cm2 s−1 and (b) ν = 1012 cm2 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 (Tobs ≈ 12 yr as of today). In both cases, rotation is uniform and the stratification is adiabatic. The lifetimes of the convective modes and high-latitude 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
thumbnail 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
thumbnail Fig. 17.

(a) Dispersion relations of the north-south ζ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 ν = κ = 1012 cm2 s−1.
In the text
thumbnail Fig. 18.

Radial velocity vr (upper plots) and entropy perturbation s1 (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 ν = κ = 1012 cm2 s−1.
In the text
thumbnail Fig. 19.

Transport properties of thermal energy and angular momentum by the north-south ζz-symmetric columnar convective modes for m = 16. (a) Correlation between radial velocity velocity and temperature perturbation ⟨vrT1⟩, (b) Reynolds stress between radial and longitudinal velocities ⟨vrvϕ⟩, and (c) Reynolds stress between latitudinal and longitudinal velocities ⟨vθvϕ⟩. The background is weakly superadiabatic (δ = 2 × 1016), rotation is uniform, and moderate diffusivities are used (ν = κ = 1012 cm2 s−1). The eigenfunctions are normalized such that the maximum radial velocity is 10 m s−1 at the equator.
In the text
thumbnail 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 theoretically-expected 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
thumbnail 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 ν = 1012 cm2 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
thumbnail 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 (ν = 1011 cm2 s−1) and strong diffusion (ν = 1012 cm2 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)−iEkm(m + 1), where Ek = 4 × 10−4 is the Ekman number at the solar surface (Fournier et al. 2022).
In the text
thumbnail 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
thumbnail 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
thumbnail 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 (ν = 1011 cm2 s−1) and strong diffusion (ν = 1012 cm2 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
thumbnail Fig. 26.

Reynolds stress components (a), (d) ρ0⟨vrvϕ⟩ 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
thumbnail Fig. 27.

Radially averaged latitudinal angular momentum fluxes. Solid, dot-dashed, and dashed lines represent the angular momentum fluxes associated with the Reynolds stress of the equatorial Rossby modes FRS, θ, advection by meridional circulation FMC, θ, and diffusion by turbulent viscosity FVD, θ, 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
thumbnail Fig. 28.

(a) Growth rate versus frequency of the ζz-antisymmetric high-latitude mode with m = 1, for different values of the latitudinal entropy difference |Δθs|=s0, pole − s0, 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 ν = κ = 1012 cm2 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
thumbnail Fig. 29.

Eigenfunctions of vr, vθ, vϕ, and s1 of the m = 1 north-south ζz-antisymmetric high-latitude 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 ν = κ = 1012 cm2 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.