© The Authors 2026

1. Introduction

Inertial modes are global oscillations in a rotating fluid whose restoring force is the Coriolis force. Different classes of inertial modes have been identified in the solar near-surface flows using various observational datasets and techniques (see Gizon et al. 2024 for a review). These include equatorial Rossby modes (Löptien et al. 2018), critical-layer and high-latitude modes (Gizon et al. 2021), and so-called high-frequency retrograde (HFR) modes (Hanson et al. 2022). All observed modes propagate retrograde with respect to the fastest local rotation rate, i.e. the surface equatorial rotation rate (or the Carrington rotation rate). The solar inertial modes are found to be very sensitive to the solar differential rotation in the convection zone (CZ; Fig. 1a), as well as to several unknown quantities in the deep interior, such as turbulent viscosity and super-adiabaticity (e.g. Gizon et al. 2021; Bekki et al. 2022b; Fournier et al. 2022; Hanson & Hanasoge 2024). Diagnosing these physical quantities using inertial modes may improve our understanding of the global-scale flow dynamics. Apart from their diagnostic role, the solar inertial modes may also play a significant dynamical role by providing non-linear feedback to regulate the Sun’s differential rotation (Bekki et al. 2024).

Fig. 1. Profiles of differential rotation used in this study. (a) Observed 2D profile in a meridional plane, taken from Larson & Schou (2018). (b) Simplified differential rotation profile, given by Eq. (5), used to model the observed differential rotation. The dotted black lines indicate the base of the CZ at r = 0.71 R⊙. The Carrington rotation rate, Ω0/2π = 456 nHz, is marked on the colour bar. (c) Cuts of observed (solid) and simplified (dashed) differential rotation profiles at fixed latitudes: 0° (red), 30° (blue), 45° (brown), and 75° (dark green). The green-shaded area denotes the CZ, the grey-shaded region the OS, and the orange-shaded region the RZ.

Linear models in 2D and 3D have helped identify the inertial modes and understand their physics (e.g. Gizon et al. 2020; Bekki et al. 2022b; Fournier et al. 2022; Bhattacharya & Hanasoge 2023; Mukhopadhyay et al. 2025). So far, most theoretical studies of inertial modes have been conducted only in the CZ and with impenetrable radial boundary conditions. It remains unclear how the radiative interior affects the properties of solar inertial modes in the convective envelope.

Blume et al. (2024) recently presented a non-linear simulation including the radiative interior and showed that both sectoral and tesseral Rossby modes are ubiquitously found in the radiative zone (RZ). It is implied that these toroidal Rossby modes are likely present in a wave cavity different from the one in the CZ. Matilsky et al. (2022) further proposed that these Rossby modes might play a role in the confinement of the tachocline or in generating the local dynamo in the RZ. Neither study, however, reported the spatial eigenfunctions of these Rossby modes, and the mode coupling between the radiative and convective zones is still not well understood.

The presence of a convectively stable stratification allows for the existence of (gravito-)inertial modes whose restoring forces are both buoyancy and Coriolis forces. Gravito-inertial modes are important for inferring rotation profiles in γ Doradus stars, which have convective cores and radiative envelopes (e.g. Ouazzani et al. 2020; Tokuno & Takata 2022). These modes can be classified into two types: inertial modes modified by buoyancy and gravity modes modified by rotation (e.g. Dintrans & Rieutord 2000). In the Sun’s radiative interior, the Brunt-Väisälä frequency is about three orders of magnitude higher than the rotation rate. Thus, most attention has so far been focused on perturbations of gravity modes by rotation (e.g. Alvan et al. 2014), rather than on inertial modes perturbed by gravity or buoyancy.

In this study, we conducted a linear eigenmode analysis for the solar case, including not only quasi-toroidal inertial modes but also non-toroidal inertial modes (HFR and prograde columnar modes), and examined how they couple to the radiative interior. Given the high precision and accuracy of solar inertial mode observations, accurate computations of both frequencies and eigenfunctions are essential for interpreting these observations. Coupling to the RZ may potentially influence the eigenmodes restricted to the CZ shell alone. We also studied Rossby modes in the radiative interior, which are perturbed by buoyancy due to the sub-adiabatic stratification. The paper is organized as follows. In Sect. 2, we describe our model of small-amplitude inertial modes, which covers both the CZ and the upper RZ. The eigenfunctions and spectra of the selected inertial modes computed in the setup are analysed in Sect. 3. Finally, we discuss the implications of our results in Sect. 4.

2. Inertial modes in a model that combines the CZ and RZ

We used the flexible spectral code Dedalus (Burns et al. 2020) to solve the linear eigenvalue problem for low-frequency modes of oscillation in the Sun. We employed the same approach as in Mukhopadhyay et al. (2025). The setup used in this paper encompasses the upper RZ and CZ of the Sun, with free-surface boundary conditions applied at the top and bottom of the domain. We used a fully compressible model, as gravity modes in the radiative interior have been found to be affected by the anelastic approximation (e.g. Brown et al. 2012; Vasil et al. 2013; Hindman & Julien 2024).

2.1. Equations for linear oscillations

We solved the linearized compressible fluid dynamical equations,

$$\begin{array}{cc}& \frac{\partial \mathit{u}}{\partial t}+{\mathit{v}}_0·\mathbf{\nabla }\mathit{u}+\mathit{u}·\mathbf{\nabla }{\mathit{v}}_0\hfill \\ \hfill & +2{\mathrm{\Omega }}_0{\mathit{e}}_z\times \mathit{u}+\frac{\mathbf{\nabla }{p}_1}{{\rho }_0}+\frac{{\rho }_1}{{\rho }_0}g{\mathit{e}}_r-\frac{1}{{\rho }_0}\mathbf{\nabla }·\left(\nu {\rho }_0\mathsf{S}\right)=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial s_1}{\partial t}+{\mathit{v}}_0·\mathbf{\nabla }{s}_1+u_r\frac{\partial s_0}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial s_0}{\partial \theta }-\frac{1}{{\rho }_0{T}_0}\mathbf{\nabla }·(\kappa {\rho }_0{T}_0\mathbf{\nabla }{s}_1)=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial {\rho }_1}{\partial t}+\mathbf{\nabla }·({\rho }_0\mathit{u})+\mathbf{\nabla }·({\rho }_1{\mathit{v}}_0)=0,\hfill \end{array}$$(3)

in combination with the linearized quasi-adiabatic equation of state (e.g. Unno et al. 1979),

$$\begin{array}{c}\hfill \frac{p_1}{p_0}={\mathrm{\Gamma }}_1\frac{{\rho }_1}{{\rho }_0}+{\mathrm{\Gamma }}_1{\mathrm{\nabla }}_{\mathrm{ad}}\frac{{\rho }_0{T}_0}{p_0}{s}_1.\end{array}$$(4)

In the above equations, u, ρ₁, p₁, and s₁ denote the perturbations of velocity, density, pressure, and entropy with respect to the background. We took the background density (ρ₀), pressure (p₀), temperature (T₀), gravitational acceleration (g), first adiabatic exponent (Γ₁), and adiabatic gradient (∇_ad) from the standard solar model, Model S (Christensen-Dalsgaard et al. 1996). The equations were formulated in the Carrington frame of reference rotating at Ω₀/2π = 456 nHz. In this frame, the base flow v₀ = (Ω(r, θ)−Ω₀) r sin θ e_ϕ denotes the Sun’s differential rotation. The background latitudinal entropy gradient (∂s₀/∂θ) is estimated under the assumption of thermal wind balance for solar differential rotation (e.g. Pedlosky 1982; Thompson et al. 2003). The remaining quantities in the model equations are described in the following subsections.

2.2. Differential rotation model

In most of our calculations, we used a simplified analytic profile for the differential rotation given by

$$\begin{array}{c}\hfill \mathrm{\Omega }(r,\theta )={\mathrm{\Omega }}_{\mathrm{CZ}}(\theta )+\frac{({\mathrm{\Omega }}_{\mathrm{RZ}}-{\mathrm{\Omega }}_{\mathrm{CZ}}(\theta ))}{2}[1+\mathrm{erf}\left(\frac{r_{\mathrm{tc}}-r}{d_{\mathrm{tc}}}\right)],\end{array}$$(5)

where r_tc = 0.71 R_⊙ and d_tc = 0.045 R_⊙ are the position and thickness of the tachocline, respectively. We assumed that the RZ rotates rigidly at Ω_RZ/2π = 433.5 nHz, while the CZ has a rotational shear in latitude as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{CZ}}(\theta )={\mathrm{\Omega }}_1+{\mathrm{\Omega }}_2{cos}^2\theta +{\mathrm{\Omega }}_3{cos}^4\theta ,\end{array}$$(6)

with Ω₁/2π = 464 nHz, Ω₂/2π = −76 nHz, and Ω₃/2π = −50 nHz. This simplified profile allowed the radial and latitudinal dependences to be separated and reduced memory usage in the Dedalus solver, which enabled us to use a higher grid resolution to better resolve modes with many radial nodes.

However, when modelling the high-latitude modes, we used the actual data of the Sun’s differential rotation measured by helioseismology (from Larson & Schou 2018). This is because the baroclinically unstable high-latitude modes are known to sensitively depend on the differential rotation profile and the associated latitudinal entropy gradient (Gizon et al. 2021; Bekki et al. 2022a). The observed and simplified profiles of the differential rotation are compared in Fig. 1. The shaded regions in Fig. 1c are described in the following subsection.

2.3. Turbulent diffusivities and super-adiabaticity

The strain-rate tensor is given by

$$\begin{array}{c}\hfill \mathsf{S}=\mathbf{\nabla }\mathit{u}+\mathbf{\nabla }{\mathit{u}}^T-\frac{2}{3}(\mathbf{\nabla }·\mathit{u}){\mathsf{I}}_3.\end{array}$$(7)

We used radially varying profiles of turbulent viscosity ν and turbulent thermal diffusivity κ, similar to Nandy & Choudhuri (2002), given by

$$\begin{array}{c}\hfill \nu (r)={\nu }_{\mathrm{RZ}}+\frac{({\nu }_{\mathrm{CZ}}-{\nu }_{\mathrm{RZ}})}{2}[1+\mathrm{erf}\left(\frac{r-r_{\mathrm{OS}}}{d_{\mathrm{OS}}}\right)],\end{array}$$(8)

$$\begin{array}{c}\hfill \kappa (r)={\kappa }_{\mathrm{RZ}}+\frac{({\kappa }_{\mathrm{CZ}}-{\kappa }_{\mathrm{RZ}})}{2}[1+\mathrm{erf}\left(\frac{r-r_{\mathrm{OS}}}{d_{\mathrm{OS}}}\right)].\end{array}$$(9)

Here, the turbulent viscosity and thermal diffusivity in the CZ (ν_CZ, κ_CZ) were both set to be 10¹² cm² s⁻¹. In the RZ below, we adopted the values ν_RZ = 10¹⁰ cm² s⁻¹ and κ_RZ = 4 × 10¹⁰ cm² s⁻¹. We assumed a slightly higher thermal diffusivity than turbulent viscosity in the RZ to account for radiative effects (e.g. Brun et al. 2011; Brown et al. 2012). We used the parameters for the transition depth r_OS = 0.7 R_⊙ and the transition thickness d_OS = 0.04 R_⊙. The left panel of Fig. 2 shows the radial profiles of ν and κ used in this study.

Fig. 2. Left: Profiles of turbulent viscosity, ν(r), and turbulent thermal diffusivity, κ(r), as functions of radius. They are expressed by Eqs. (8) and (9). Right: Profile of super-adiabaticity, δ(r), as a function of radius (Eq. 10). The super-adiabaticity profile from the standard Model S is represented by the dashed red curve (Christensen-Dalsgaard et al. 1996). The y-axis is linear for |δ|< 10−7 and logarithmic beyond that. In both panels, the green-shaded area denotes the CZ, the grey-shaded region the OS, and the orange-shaded region the RZ.

The background radial entropy gradient, ∂s₀/∂r, denotes the degree of convective instability in the background stratification. Although the RZ is strongly sub-adiabatic (convectively stable), uncertainty remains regarding super-adiabaticity in the bulk CZ. While Model S predicts a super-adiabatic stratification throughout the whole CZ, this prediction has been questioned by recent observations and numerical models of solar large-scale convection (e.g. Käpylä et al. 2017; Gizon et al. 2021; Hotta et al. 2022; Bekki 2024). In this study, we employed the following step function for ∂s₀/∂r, which smoothly connects the strongly sub-adiabatic RZ and the adiabatic bulk CZ,

$$\begin{array}{c}\hfill \frac{\partial s_0}{\partial r}=\frac{1}{4}{\left(\frac{\partial s_0}{\partial r}\right)}_{\mathrm{RZ}}[1-\mathrm{erf}\left(\frac{r-r_{\delta 1}}{d_{\delta }}\right)][1-\mathrm{erf}\left(\frac{r-r_{\delta 2}}{d_{\delta }}\right)],\end{array}$$(10)

where r_δ1 = 0.61 R_⊙, r_δ2 = 0.68 R_⊙, d_δ = 0.04 R_⊙, and (∂s₀/∂r)_RZ = 1.32 × 10⁻² erg g⁻¹ K⁻¹ cm⁻¹. The right panel of Fig. 2 shows the corresponding radial profile of super-adiabaticity δ = −H_pc_p⁻¹(∂s₀/∂r) used in our computations. We note that our δ(r) profile has a smoother transition between the sub-adiabatic RZ and the adiabatic CZ than that of Model S. This is required in our model to avoid oscillations in the projection of the profile onto the Chebyshev basis.

Based on the profiles of turbulent diffusivities and super-adiabaticity shown in Fig. 2, our numerical domain can be divided into three regions:

• the RZ (r < 0.64 R_⊙), where turbulent motions are suppressed (ν and κ are substantially small), and the stratification is strongly sub-adiabatic (δ ∼ −0.1),

• the overshooting (OS) layer (0.64 R_⊙ ≤ r ≤ 0.71 R_⊙), where turbulent motions penetrate but the stratification is significantly sub-adiabatic (e.g. Hotta 2017; Käpylä et al. 2017),

• the CZ (r > 0.71 R_⊙) where the vigorous convective motions lead to significant values of ν and κ and the stratification is close to adiabatic (|δ|≲10⁻⁶).

These three regions are highlighted by different colours in Fig. 2. Figure 1c demonstrates that the tachocline with significant radial differential rotation overlaps with the OS defined in our setup.

2.4. Free-surface boundary conditions

We set the lower boundary in the upper RZ at r_i = 0.5 R_⊙ and the upper boundary slightly below the photosphere at r_o = 0.985 R_⊙. At both radial boundaries, in contrast to the conventionally used impenetrable boundary condition, we imposed the free-surface boundary condition for the radial velocity u_r. This is done in the presence of base flow v₀ by setting the Lagrangian derivative of total pressure to 0 (e.g. Mei et al. 2005), as

$$\begin{array}{c}\hfill \frac{D}{\mathit{Dt}}(p_0+p_1)\approx \frac{\partial p_1}{\partial t}+\mathit{u}·\mathbf{\nabla }{p}_0+{\mathit{v}}_0·\mathbf{\nabla }{p}_1=0.\end{array}$$(11)

Using the wave Ansatz where all the perturbations are assumed to be proportional to exp(imϕ − iωt), this translates to

$$\begin{array}{c}\hfill -\mathrm{i}\omega p_1-u_r{\rho }_{0}g+\mathrm{i}m(\mathrm{\Omega }(r,\theta )-{\mathrm{\Omega }}_0)p_1=0.\end{array}$$(12)

Equation (12) was forced at both radial boundaries. We validated this free-surface boundary condition by reproducing the dispersion relation of the f modes on a spherical surface (e.g. Christensen-Dalsgaard 2002; Antia 1998). For low-frequency inertial modes whose oscillation periods are much longer than those of f modes, the free-surface boundary condition is found to almost reduce to the impenetrability condition (e.g. Hindman & Jain 2022, 2023). Nevertheless, the free-surface boundary condition is expected to work better for g modes and gravito-inertial modes in the RZ. In addition, we applied the conventional horizontal stress-free condition, i.e. S_rθ = S_rϕ = 0, at both radial boundaries. We assumed that the flux of entropy perturbations vanishes at the bottom boundary, and entropy perturbations vanish at the top boundary. We note that this treatment is consistent with the free-surface boundary condition, given that the radial entropy gradient ∂s₀/∂r vanishes at the top boundary in our setup.

2.5. Eigenvalue solver

We followed the same approach as in Mukhopadhyay et al. (2025) to solve the eigenvalue problem of the governing equations (Eqs. 1–4) with the boundary conditions (Sect. 2.4), considering the wave Ansatz exp(imϕ − iωt) for each azimuthal order m. We solved the sparse problem on the spherical shell basis of Dedalus and selected the modes of interest, in the same way as in Mukhopadhyay et al. (2025). We used a typical resolution of 180 points in radius and 60 points in latitude in calculations with the simplified differential rotation. When the observed solar differential rotation was adopted, we used a grid with 90 radial points and 36 latitudinal points. All eigenmodes presented in the paper are numerically converged to errors well below the observational uncertainties (Gizon et al. 2021).

3. Inertial modes of the extended model

3.1. Modes confined primarily to the CZ

First, we examined how the inertial modes of the CZ are affected by the presence of the RZ below. To this end, we performed calculations with and without the RZ and compared the resulting eigenmodes. The model without the RZ has the lower boundary at r_i = 0.71 R_⊙ with an impenetrable boundary condition, thereby excluding both the OS and the RZ. This CZ-only model has been traditionally employed in previous studies (Bekki et al. 2022b; Mukhopadhyay et al. 2025). Except for this lower boundary treatment, the two models are identical.

Figures 3a and 3b compare the eigenfunctions of the m = 1 high-latitude mode, m = 3 equatorial Rossby modes with no radial node (n_CZ = 0) and with one radial node in the CZ (n_CZ = 1), the m = 10 HFR mode, and the m = 3 prograde columnar mode. All of these modes have already been discussed in Mukhopadhyay et al. (2025) within the CZ-only framework. We find that the inclusion of the RZ has a limited impact on the eigenfunctions of these inertial modes within the CZ. In fact, these modes are found to decay in the RZ (Fig. 3b). Figure 3c compares the surface velocity eigenfunctions between the two models, where the eigenfunctions are normalized such that the modes have the same total horizontal velocity power at the surface in the two models. The differences in the velocity eigenfunctions at the surface (r = r_o) can be measured by the correlation coefficient

$$\begin{array}{c}\hfill \mathcal{C}=\frac{{\int }_S(u_{\theta ,\mathrm{RZ}}{u}_{\theta ,\mathrm{CZ}}+u_{\varphi ,\mathrm{RZ}}{u}_{\varphi ,\mathrm{CZ}})dS}{\sqrt{{\int }_S(u_{\theta ,\mathrm{RZ}}^2+u_{\varphi ,\mathrm{RZ}}^2)dS}\sqrt{{\int }_S(u_{\theta ,\mathrm{CZ}}^2+u_{\varphi ,\mathrm{CZ}}^2)dS}},\end{array}$$(13)

Fig. 3. Comparison of selected inertial modes computed in the setups that include or exclude the RZ. (a) Velocity eigenfunctions in meridional cross-sections from the CZ-only model. The real part of the eigenfunction corresponds to a longitude ϕ0 (where uθ is maximum), while the imaginary part corresponds to the longitude ϕ0 − π/2m. The m = 1 high-latitude mode is normalized to have the maximum surface velocity of 10 m s−1, while the other modes are normalized to have a maximum surface velocity of 1 m s−1. Solid black curves denote the critical layers of the mode where ℜ[ω]=m(Ω − Ω0). The frequencies are measured in the Carrington frame. (b) Same as panel a but from the extended model that includes the RZ. Dotted black lines denote the base of the CZ. (c) Horizontal velocity eigenfunctions as functions of latitude at the surface. Dashed blue and solid red curves represent the results with and without the RZ, respectively. The correlation coefficient between the eigenfunctions, 𝒞, defined in Eq. (13), is noted above each subplot. (d) Radial profiles of the RMS velocity.

where S represents the spherical surface and u_RZ and u_CZ are from the models with and without the RZ, respectively. We find that 𝒞 ≳ 0.9 for the inertial modes discussed here, demonstrating that the way these inertial modes are observed at the top of the CZ is almost unaffected by the inclusion of the RZ. Figure 3d compares the radial variations of mode power measured by root-mean-square (RMS) velocity u_RMS. It is shown that the power distribution in the bulk CZ is only marginally affected by the inclusion of the RZ, except near the base of the CZ, where the deviation can become significant (e.g. m = 1 high-latitude mode). Although the modes decay in the RZ, they can exhibit finite motion in the OS. This can affect mode frequencies and growth rates by altering mode dissipation (see Appendix A). Figure 3 suggests that CZ inertial modes are generally decoupled from the RZ, as modes with finite radial motions are excluded from the strongly sub-adiabatic RZ. Although purely toroidal modes can form global eigenfunctions under uniform rotation (Appendix B), the inclusion of differential rotation inevitably makes the CZ modes decoupled from the RZ. Thus, we note that the reported decoupling of CZ inertial modes from the RZ is a generic result, largely independent of model details such as viscosity profiles and lower boundary conditions.

Figure 4 compares the frequencies and growth rates of the above-discussed inertial modes in the two setups that include and exclude the RZ. The inclusion of RZ affects the mode frequencies by at most a few nanohertz (< 10 nHz) for most inertial modes except some non-toroidal modes (HFR modes at m ≥ 12 and prograde columnar modes at m ≤ 6). The differences in frequencies are remarkably large (> 20 nHz) only for HFR modes at m > 12. We note that these HFR modes are special in several aspects: they are strongly non-toroidal, and their radial motions tend to be more localized near the base of the CZ as m increases (Bekki 2024). Furthermore, their critical layers penetrate into the OS (see Fig. 3b). These various factors change the effective rotation rate associated with the mode and cause an apparent shift in the frequency measured in the Carrington frame (refer to Appendix C for details). This Doppler shift of the frequency in the Carrington frame is further amplified for high m. In contrast to Hindman & Jain (2023), the dispersion relation of the prograde columnar modes in our setup is not significantly affected by the inclusion of the RZ because the strongly sub-adiabatic RZ itself acts like a lower boundary to these convective modes.

Fig. 4. Top: Dispersion relations of the inertial modes in the CZ computed in the setups that include the RZ (solid lines with points) and exclude the RZ (dashed lines with open circles). The colours denote the various types of inertial modes. Bottom: Growth rates of the same modes, with the same notations.

Although the changes in frequencies are rather insignificant for most inertial modes, the inclusion of the RZ has a stronger impact on the mode growth rates. As shown in the bottom panel of Fig. 4, the modes tend to be more damped when the RZ is included. This results from the dissipation in the OS caused by the inclusion of the RZ (see Fig. A.1). For high m, the modes tend to be confined to the outer part of the CZ, so their growth rates are not significantly affected by the inclusion of the RZ. This is not the case for the HFR modes, which still exhibit significant motions throughout the CZ at m = 16, as their critical layers are located near the CZ base.

3.2. Modes confined primarily to the RZ

Next, we discuss the modes that have dominant motions in the RZ. We find that there exist quasi-toroidal inertial modes trapped inside the RZ whose eigenfunctions are well represented by the spherical harmonics and follow the dispersion relations of the classical Rossby modes

$$\begin{array}{c}\hfill {\omega }_{\ell ,m}^{\mathrm{R}}=-\frac{2m{\mathrm{\Omega }}_{\mathrm{RZ}}}{\ell (\ell +1)}+m({\mathrm{\Omega }}_{\mathrm{RZ}}-{\mathrm{\Omega }}_0),\end{array}$$(14)

where ℓ is the harmonic degree. Here, the second term represents the Doppler frequency shift in the Carrington frame. For each combination of (ℓ, m), these RZ Rossby modes can have any number of radial nodes n_RZ while sharing the same frequencies ${\omega }_{\ell ,m}^{\mathrm{R}}$, where n_RZ is counted for u_θ in the region 0.5 R_⊙ ≤ r ≤ 0.71 R_⊙. Figure 5 displays their dispersion relations for ℓ − m = 0, 1, and 2, with two different n_RZ, clearly demonstrating that the computed frequencies agree strikingly well with the theoretical prediction (Eq. 14) and that the frequencies are independent of n_RZ.

Fig. 5. Dispersion relations of the Rossby modes inside the RZ with ℓ = m (navy blue), ℓ = m + 1 (red), and ℓ = m + 2 (dark green) for azimuthal orders 1 ≤ m ≤ 16. The frequencies are measured in the Carrington frame. Plus symbols (nRZ = 0) and open circles (nRZ = 4) denote modes with different numbers of radial nodes, nRZ, in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. Solid grey curves represent the theoretical dispersion relations of the classical Rossby modes ${\omega }_{\ell ,m}^{\mathrm{R}}$ given by Eq. (14). The dashed black line denotes the maximum possible value of ${\omega }_{\ell ,m}^{\mathrm{R}}$.

Figure 6a shows the complex eigenfrequency spectrum at ℓ = m = 3. The equatorial Rossby mode inside the CZ, discussed in Sect. 3.1, can be identified at ℜ[ω]/2π = −263.3 nHz (denoted by the cross symbol). Apart from that, near the classical Rossby mode frequency ${\omega }_{\ell =m=3}^{\mathrm{R}}$, there exist many RZ Rossby modes with different n_RZ (denoted by star symbols). Their growth rates decrease with the increase in n_RZ as dissipation increases with n_RZ (see Appendix A). Figures 6b and 6c further show the meridional eigenfunctions of the RZ Rossby modes with n_RZ = 0 and 3, respectively. The eigenfunctions indicate that the modes are strongly confined within the RZ and decay in the CZ. Since the RZ is uniformly rotating and no critical layers exist, the horizontal eigenfunctions in the RZ show no major deviations from the corresponding spherical harmonics. We also note that these RZ Rossby modes have substantial entropy perturbations arising from the strongly sub-adiabatic background stratification.

Fig. 6. Position of the RZ Rossby modes with ℓ = m = 3 in the complex frequency space and their eigenfunctions. (a) Complex frequency spectrum of equatorial Rossby modes in the RZ with azimuthal order m = 3. Star symbols denote the RZ Rossby modes, while the cross symbol denotes the CZ Rossby mode. The red (blue) star marks the RZ Rossby mode with the number of radial nodes, nRZ = 0 (nRZ = 3), in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. The vertical solid blue line denotes the observed Rossby mode frequency, while the vertical dashed green line denotes the frequency given by Eq. (14). (b) Meridional eigenfunction of the RZ Rossby mode with nRZ = 0 (red star). The longitudes corresponding to the real and imaginary components are designated in the same way as in Fig. 3. Dotted lines denote the position of the base of the CZ. Eigenfunctions are normalized such that the maximum of uθ is 1 m s−1. The unit of s1 is erg g−1 K−1. (c) Same as b but for the RZ Rossby mode with nRZ = 3 (blue star).

As has been reported in Blume et al. (2024), the Rossby modes in the sub-adiabatic RZ can have any allowed combination of (ℓ, m), including the non-sectoral (ℓ ≠ m) components (see also Fig. 5). This is in striking contrast to the Rossby modes in the adiabatically stratified CZ, where only sectoral (ℓ = m) modes can exist, as only sectoral modes can sustain the required radial force balance between the Coriolis and pressure gradient forces (Provost et al. 1981; Damiani et al. 2020). In the sub-adiabatic RZ, on the other hand, the radial component of the Coriolis force can be balanced by the combined effect of buoyancy and pressure to maintain the toroidal nature of these Rossby modes. To confirm this, we analysed the force balance in the RZ Rossby modes. The various forces are given by

$$\begin{array}{cc}& {\mathit{F}}_{\mathrm{Coriolis}}=-2{\mathrm{\Omega }}_0({\mathit{e}}_z\times \mathit{u}),\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\mathit{F}}_{\mathrm{Pressure}}=-\mathbf{\nabla }\left(\frac{p_1}{{\rho }_0}\right),\hfill \end{array}$$(16)

$$\begin{array}{cc}& {\mathit{F}}_{\mathrm{Buoyancy}}=\frac{s_1}{c_p}g{\mathit{e}}_r,\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\mathit{F}}_{\mathrm{Viscous}}=\frac{1}{{\rho }_0}\mathbf{\nabla }·\left(\nu {\rho }_0\mathsf{S}\right).\hfill \end{array}$$(18)

Figure 7 presents the force balance analysis for a non-sectoral RZ Rossby mode with (ℓ,m) = (4, 3). Figure 7a shows that the mode is quasi-toroidal, where the radial motions are several orders of magnitude smaller than the horizontal motions. However, the small radial motions induce a substantial buoyancy force in the sub-adiabatic RZ to balance the radial pressure gradient and Coriolis forces (Fig. 7b). Nevertheless, the geostrophic nature of the traditional Rossby modes persists in the horizontal force balance (Fig. 7c). This causes the Rossby modes in the RZ to have frequencies determined entirely by the rotation rate, given by Eq. (14).

Fig. 7. Analyses of the non-sectoral Rossby mode in the RZ with ℓRZ = 4, m = 3, nRZ = 3. (a) Radial profiles of the RMS velocity in the horizontal (red) and radial (blue) directions. The profiles are normalized to have a maximum total RMS velocity of unity. (b) RMS of the radial component of various forces, defined in Eqs. (15)–(18). The pressure gradient force (solid blue), Coriolis force (dashed black), buoyancy force (dashed red), and viscous force (solid orange) are normalized so that the maximum of the pressure gradient force is unity. (c) Same as (b) but for the horizontal component of the forces. Colours of the shaded regions follow the same notation as in Fig. 2.

3.3. CZ–RZ mixed modes

We now explore modes with significant motion in both the CZ and RZ. Although we call them CZ–RZ mixed modes in this paper, they correspond to low-frequency modes arising from the mixing between Rossby modes in the RZ and inertial modes in the CZ, and should not be confused with the p − g mixed modes commonly discussed in asteroseismology. In practice, we find that significant mode mixing only occurs when |Δω|≲10 nHz, where |Δω| is the frequency separation between the RZ Rossby mode and the nearest CZ inertial mode, consistent with avoided-crossing behaviour. Such mode-mixing can occur only when the RZ Rossby mode and its nearest CZ inertial mode have the same north–south symmetry. In such special cases, we encounter CZ–RZ mixed modes, whose eigenfunctions have substantial amplitudes in both the CZ and RZ. A caveat should be noted: Although we refer to them as mixed modes, they are still essentially RZ Rossby modes that have non-negligible amplitudes in the CZ, and only a subset of RZ Rossby modes are allowed to behave in this way. Figure 8 highlights 11 examples illustrating how mode mixing between the RZ and CZ can occur in the dispersion diagrams for both equatorially symmetric and antisymmetric modes. We note that the possible options of mode mixing are quite limited, as the frequencies of the RZ Rossby modes are restricted to a finite range by Eq. (14). For instance, the mixing between the ℓ_RZ = m Rossby modes in the RZ and the n_CZ = 0 equatorial Rossby modes in the CZ can only occur at m = 1 and 2. Here, ℓ_RZ denotes the harmonic degree of the RZ Rossby component of these mixed modes. The n_CZ = 0 Rossby modes in the CZ with m = 5 and 6 can also couple with the RZ Rossby modes but with ℓ_RZ − m = 2 and 6, respectively. A similar coupling is allowed for the n_CZ = 1 Rossby modes in the CZ with m = 4 and 5. As for the high-latitude modes, only the m = 1 mode with north–south symmetric vorticity and the m = 2 mode with the opposite symmetry can couple with the RZ Rossby modes. Finally, the HFR modes with m = 7, 8, and 9 are allowed to couple with the RZ Rossby modes with north–south antisymmetric vorticity. The mode masses of all these mixed modes are, however, dominated by the RZ contribution (as quantified in Appendix D). We note that there is no mixed mode between the retrograde-propagating RZ Rossby modes and the prograde columnar modes.

Fig. 8. Frequencies of the retrograde-propagating inertial modes in the CZ (dashed coloured; the same as in Fig. 4) overplotted with the dispersion relations of the RZ Rossby modes (solid grey; given by Eq. (14)). The left and right panels correspond to the modes with north–south symmetric and antisymmetric vorticity, respectively. Open circles mark the subset of RZ Rossby modes that can form mixed modes with the CZ inertial modes.

In principle, CZ–RZ mixed modes are allowed to have any number of radial nodes n_RZ in the region 0.5 R_⊙ ≤ r ≤ 0.71 R_⊙, since the frequencies of the RZ Rossby modes are independent of n_RZ. However, n_RZ has a significant impact on how the mode power is distributed between the CZ and RZ: As n_RZ increases, the relative velocity amplitude in the CZ increases. A possible explanation for this is given in Appendix A. Figure 9 displays the normalized kinetic energy within the CZ of the RZ Rossby modes as a function of the absolute frequency difference |Δω| between each RZ Rossby mode and its nearest CZ inertial mode. It clearly shows that the mode tends to be more strongly coupled and has more power inside the CZ when |Δω| is smaller. It is also demonstrated that the fraction of kinetic energy in the CZ increases with n_RZ. These results collectively indicate that, for an RZ Rossby mode to be a CZ–RZ mixed mode, it must have high n_RZ as well as a frequency sufficiently close (|Δω|≲10 nHz) to that of the corresponding CZ inertial mode. Note that the threshold |Δω|≲10 nHz for coupling is a rough estimate, phenomenologically obtained for the setup using Fig. 9.

Fig. 9. Fraction of kinetic energy density, E = ⟨ρ0u2⟩, of the RZ Rossby modes contained in the CZ as a function of the absolute difference between the eigenfrequency of each mode and that of the nearest CZ inertial mode (see Fig. 8). For clarity, we show only the RZ Rossby modes with ℓRZ = m (circles), ℓRZ = m + 1 (diamonds), and ℓRZ = m + 2 (squares). The colour of the symbols denotes the number of radial nodes (nRZ) in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. Only modes with nRZ ≤ 5 are shown.

The velocity eigenfunctions of selected CZ–RZ mixed modes are compared in Fig. 10 with those of the CZ inertial modes (discussed in Sect. 3.1), showing overall remarkable similarity in the CZ. Notably, the surface eigenfunctions exhibit only minor differences, suggesting that it is challenging to determine from surface observations alone whether the mode is mixed with RZ Rossby modes or not. As a robustness check, we verified numerical convergence for two representative CZ–RZ mixed modes that exhibit strong spatial variations in radial and latitudinal directions (see Appendix E).

Fig. 10. Selected inertial modes computed in the extended setup, which includes the RZ. (a) Modes with motions predominantly in the CZ and decaying in the RZ. (b) Modes with comparable motions in the RZ and CZ, with flows in the CZ similar to those of the corresponding modes in panel a. (c) and (d): Surface and radial variations of the modes, respectively, as in Fig. 3. The spherical harmonic degree and the number of radial nodes of the RZ Rossby modes in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙ are denoted by ℓRZ and nRZ, respectively. The m = 1 high-latitude modes are normalized to have a maximum surface velocity of 10 m s−1, the m = 2 high-latitude modes to 4 m s−1, and the other modes to 1 m s−1. The notation and layout are the same as in Fig. 3.

4. Discussion

4.1. Most CZ modes are weakly affected by the RZ

In this work, we extended the linear model of inertial modes in the Sun’s differentially rotating CZ of Mukhopadhyay et al. (2025) to include the sub-adiabatic RZ. To study how the RZ affects the solar inertial modes in the CZ, we first compared the eigenmodes computed in models that include and exclude the RZ. We find that the mode frequencies are only marginally affected (a few nanohertz) by the inclusion of the RZ, besides some non-toroidal modes (HFR modes at m ≥ 12 and prograde columnar modes at m ≤ 6). In the model with the RZ, the eigenfunctions are found to slowly decay in the OS before vanishing in the RZ. The associated viscous dissipation in the OS leads to larger damping rates of these inertial modes. Within the CZ, the eigenfunctions are primarily altered near the base of the CZ. On the other hand, the surface eigenfunctions are shown to be almost insensitive to the inclusion of the RZ. Therefore, while the RZ introduces only minor modifications, the CZ-only model remains a valid and practical framework for studying the dispersion relations and surface eigenfunctions of the solar inertial modes and for interpreting the solar observations. If very accurate frequencies and eigenfunctions are desired, extending the computational domain down to 0.55 R_⊙ will suffice as the eigenfunctions vanish by 0.6 R_⊙ (see Fig. 3). Finally, it should be noted that we used a more gradually varying super-adiabaticity profile than that predicted by Model S, corresponding to a thicker OS layer. A more realistic OS model will be needed for better modelling.

4.2. RZ Rossby modes and their coupling with CZ modes

As found by Blume et al. (2024), the sub-adiabatic and uniformly rotating RZ can host a wide variety of quasi-toroidal Rossby modes, including the non-sectoral ones (which do not exist in the quasi-adiabatic CZ). We find that they follow the classical Rossby mode dispersion relation (Eq. 14) and have horizontal eigenfunctions well represented by spherical harmonics. They can have all possible combinations of (ℓ_RZ, m, n_RZ), because the radial force balance can be maintained with the help of buoyancy. However, since their frequencies are entirely dictated by their rotation rate, they should be interpreted as inertial modes perturbed by buoyancy. The perturbations in structure and frequency of the modes due to buoyancy are rather minor because the Brunt-Väisälä frequency is significantly larger than the rotation rate (Dintrans & Rieutord 2000).

As the number of radial nodes (n_RZ) increases, the RZ Rossby modes are more damped and, in turn, have more power in the CZ (see Appendix A). We find that, when the frequencies of the RZ Rossby modes are within ∼10 nHz of the nearest CZ inertial mode with the same north–south symmetry, these RZ Rossby modes couple with the CZ and develop substantial flows in the CZ remarkably similar to those of the CZ inertial modes. We refer to such modes as CZ–RZ mixed modes. We remark that these couplings would have been very different if one used uniform rotation throughout the domain, as the frequency of the RZ Rossby modes is entirely determined by the rotation rate in the RZ (see Appendix B). Blume et al. (2024) discusses two separate Rossby wave cavities in the stably stratified RZ and adiabatic CZ. Our studies suggest that the separation of Rossby-mode frequencies in the CZ and RZ is primarily due to differential rotation. Furthermore, CZ Rossby modes cannot mix with sectoral Rossby modes in the RZ for m > 2 due to their different frequencies. Nevertheless, the RZ Rossby modes can mix with other inertial modes in the CZ under the concurrent conditions identified in our work.

4.3. Detectability of CZ–RZ mixed modes

One may ask if any of the CZ–RZ mixed inertial modes can be observed, as many of these modes have comparable growth rates to the CZ inertial modes (see Fig. 10). The surface eigenfunctions and frequencies of many CZ inertial modes and their nearest CZ–RZ mixed modes are likely too similar to be distinguished using surface observations alone. Figure 11 shows the frequencies of the observed solar inertial modes. Many of the observed modes have frequencies close to those of the RZ Rossby modes, which share the same north–south symmetry. This allows them to couple with the RZ Rossby modes. It may be that the marked modes with close frequencies are weakly mixed with the RZ Rossby modes. Many such modes include critical-layer modes, which could plausibly be extensions of some RZ Rossby modes with high n_RZ.

Fig. 11. Frequencies of the observed solar inertial modes (open circles; from Löptien et al. 2018; Gizon et al. 2021; Hanson et al. 2022) overplotted with the dispersion relations of the RZ Rossby modes (solid grey; given by Eq. (14)). The left and right panels correspond to modes with north–south symmetric and antisymmetric vorticity, respectively. The observed inertial modes are high-latitude (blue), Rossby (red), HFR (green), and critical-layer (magenta) modes. Observed modes whose frequencies are within 10 nHz of an RZ Rossby mode frequency with the same symmetry (|Δω|≲10 nHz) are shown with a translucent fill.

A model for the excitation and damping of these modes would certainly help determine whether some of the mixed modes can reach observable amplitudes. For the linearly stable modes, stochastic excitation by convection is likely the dominant mechanism (Philidet & Gizon 2023). However, as reported in Table 1 (see also Appendix D), the mixed modes are shown to possess a mode mass (see e.g. Christensen-Dalsgaard 2002; Baudin et al. 2005) much higher than that of the CZ inertial modes. This makes it energetically difficult to excite such mixed modes by turbulent convection in the CZ. For linearly unstable modes, dedicated non-linear simulations will be required to reliably estimate saturation amplitudes (Bekki et al. 2024). Latitudinal shear instability in the tachocline may excite some RZ Rossby modes (Garaud 2001). However, we do not find such linear instabilities in our study.

4.4. Implications of horizontal motions below the CZ

We find that the CZ inertial modes have notably significant motions in the OS (see Fig. 3). The m = 1 high-latitude mode has RMS velocity ∼3 m s⁻¹ in the OS, when normalized to the observed amplitude at the surface. The m = 3 equatorial Rossby mode also has RMS velocity ∼0.7 m s⁻¹ in the OS, while the other modes shown in Fig. 3 have RMS velocities ≲0.2 m s⁻¹ in the OS. Furthermore, CZ–RZ mixed modes have higher velocities in the RZ as compared to the CZ because of less dissipation in the RZ (refer to Appendix A). The presence of these (mostly horizontal) motions in the sub-adiabatic OS may contribute significantly to horizontal turbulent diffusion in the solar tachocline (see discussion in Garaud et al. 2025). They could also aid in the confinement of the tachocline through local dynamo effects (e.g. Matilsky et al. 2022).

Acknowledgments

Author contributions: LG and XZ initiated this project; SM implemented all the equations in Dedalus and performed the computations and analysis; YB provided supervision to validate the results; all authors discussed the results; SM wrote the initial draft, and all authors contributed to the final manuscript. We thank P. Dey, D. Fournier, and R. H. Cameron for helpful discussions. SM is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen. YB and LG acknowledge support from ERC Synergy Grant WHOLE SUN 810218. XZ acknowledges the financial support from the German Research Foundation (DFG) through grants 521319293, 540422505, and 550262949.

References

Appendix A: Effects of dissipation in the RZ

In this appendix, we examine the effects of dissipation in the RZ. For many of the CZ–RZ mixed modes presented in the paper, the velocity amplitude tends to be larger in the RZ than in the CZ. We provide a physical explanation for this behaviour. We also illustrate why the CZ inertial modes are more damped when the RZ is included.

The equation for the kinetic energy density e_kin = (ρ₀/2)|u|² is obtained by taking a dot product of Eq. (1) with u. Neglecting the advection terms, this yields

$$\begin{array}{c}\hfill \frac{\partial e_{\mathrm{kin}}}{\partial t}\approx -\mathbf{\nabla }·(p_1\mathit{u}+\nu {\rho }_0\mathit{u}·\mathsf{S})-{\rho }_{1}g{u}_r-Q_{\mathrm{vis}},\end{array}$$(A.1)

where Q_vis denotes the viscous dissipation of kinetic energy:

$$\begin{array}{c}\hfill Q_{\mathrm{vis}}=\frac{{\rho }_0}{2}\nu \mathsf{S}:\mathsf{S}.\end{array}$$(A.2)

By taking a volume integral under appropriate radial boundary conditions, we obtain the equation for the total kinetic energy of the toroidal Rossby modes (with u_r ≈ 0) as

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\displaystyle {\int }_V}{e}_{\mathrm{kin}}dV\approx -{\displaystyle {\int }_V}{Q}_{\mathrm{vis}}dV.\end{array}$$(A.3)

Therefore, the mode damping rate γ = −ℑ[ω] can be expressed as

$$\begin{array}{c}\hfill 2\gamma \approx \frac{{\displaystyle {\int }_V}{Q}_{\mathrm{vis}}dV}{{\displaystyle {\int }_V}{e}_{\mathrm{kin}}dV}.\end{array}$$(A.4)

The order-of-magnitude estimate S : S ≈ k²|u|², where k is the local wavenumber, implies that the damping γ increases with the number of radial nodes.

The mode damping rate is given as the radial average of νk² weighted by the kinetic energy density e_kin. In general, there is a tendency for eigenmodes to distribute more (less) kinetic energy into a region where the dissipation νk² is small (large) in order to minimize the mode damping rate γ. The RZ, having a lower ν, tends to have higher kinetic energy for CZ–RZ mixed modes with a low number of radial nodes. This further helps explain the increase in motions in the CZ for CZ–RZ mixed modes as the number of radial nodes in the RZ increases (because k is increased). Thereby, the RZ Rossby modes also become more damped with the increase in n_RZ (see Fig. 6).

Figures 3 and 4 demonstrate that the CZ inertial modes are significantly damped by the inclusion of the RZ. We demonstrate that this is caused by dissipation arising from the decay of modes in OS. Figure A.1 shows the spherically averaged radial profile of dissipation Q_vis for the m = 3 prograde columnar mode, which is the most affected mode among the CZ inertial modes whose eigenfunctions are presented in Fig. 3. The radial profile of Q_vis demonstrates that the dissipation is much higher in OS as compared to the whole of the CZ. This is because of the sharp decay of the modes in OS, where the viscosity is finite. The generated dissipation significantly damps the CZ inertial modes when the RZ is included (see Fig. 4).

Fig. A.1. Radial profile of dissipation, Qvis, defined in Eq. (A.2), for the m = 3 prograde columnar mode in the models with (red) and without the RZ (blue). The eigenmodes are shown in Fig. 3. They are normalized to have the same total energy density in the CZ.

Appendix B: Case of uniform rotation

In Sect. 3.2, we show that the Rossby modes in the CZ and RZ are separated in frequency because of differential rotation. This appendix discusses the case of uniform rotation throughout the Sun, i.e. Ω(r, θ) = Ω₀, while maintaining the radial profiles of all other background quantities. We focus here on the sectoral (ℓ = m) Rossby modes, because the non-sectoral modes cannot exist in the CZ. Figure B.1 compares the eigenfunctions of the sectoral modes with 1 ≤ m ≤ 4 under uniform rotation and those under differential rotation. In the case of uniform rotation, the mode frequencies are given by ω_m = −2Ω₀/(m + 1) in both the CZ and RZ, so the eigenmodes extend globally. In other words, there is essentially no cavity separating the Rossby modes in the RZ from those in the CZ. We note, however, that the underlying force balance maintaining their toroidal nature is different in these two regions: The radial pressure gradient force is balanced by the Coriolis force in the CZ but by buoyancy in the RZ (see Fig. 7). Consequently, the radial eigenfunctions do not show the well-known r^m dependence (predicted assuming the geostrophic balance) inside the RZ. Figure B.1 also shows that the mode tends to have a larger velocity amplitude in the RZ than in the CZ as m increases. This is likely because, when the local wavenumber k increases with m, the eigenmode tends to concentrate more of its kinetic energy in the less viscous RZ in order to minimize the damping rate (see Appendix A).

Fig. B.1. Analysis of the effects of differential rotation on sectoral Rossby modes with no nodes in radius (nRZ = nCZ = 0, ℓRZ = m). (a) Meridional cross-sections of uθ and uϕ of the Rossby modes computed under uniform rotation for different azimuthal orders, m. The longitudes corresponding to the real and imaginary components are designated in the same way as in Fig. 3. The frequencies and the growth rates measured in the Carrington frame are noted below the eigenfunctions. All eigenfunctions are normalized to have a maximum uθ of 1 m s−1 at the surface. The dotted lines indicate the position of the CZ base. (b) Same as (a) but under the simplified differential rotation given by Eq. (5), shown in Fig. 1. (c) Radial profiles of the RMS velocity of the modes computed with and without differential rotation. The colours of the shaded regions have the same meaning as in Fig. 2.

In the case of differential rotation in the CZ, on the other hand, the situation is markedly different. The sectoral Rossby modes in the RZ and the equatorial Rossby modes in the CZ (distorted from the original sectoral spherical harmonics by critical layers) have different frequencies, and thus they cannot form a single global eigenmode but instead exist only as isolated, distinct modes. A similar result has been found by Blume et al. (2024), who reported that the CZ and RZ host separate wave cavities for Rossby modes. We note that the only exception is the ℓ = m = 1 mode, in which the CZ and RZ are strongly coupled. This particular mode, known as the spin-over mode (e.g. Kerswell 1993), has a fixed frequency of −Ω₀ and exhibits global motions with quasi-uniform vorticity in a direction perpendicular to the rotational axis (Greenspan 1968). Owing to its special flow structure, this mode experiences negligible viscous dissipation. Consequently, unlike the other modes discussed above, the spin-over mode does not need to redistribute its kinetic energy in response to local viscosity, allowing it to have nearly the same eigenfunction in the CZ and the RZ. This mode, however, is difficult to detect because its frequency is nearly zero in an inertial frame (−31.7 nHz in the Earth’s frame), placing it in a spectral region often dominated by systematics and signals associated with Earth’s revolution around the Sun.

Appendix C: Analysis of HFR modes

In this appendix, we investigate the origin of the significant difference in the dispersion relations of HFR modes at m ≥ 12 between the models that include and exclude the RZ (see Fig. 4). We analyse the effective rotation rate of HFR modes with 10 ≤ m ≤ 16 because the dispersion relations begin to diverge at m = 10. Although the differences are small for m < 12, they become large, exceeding 20 nHz, only for m > 12. The effective rotation rate of the mode can be measured by

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{eff}}=\frac{\int \mathrm{\Omega }(r,\theta )e_{\mathrm{kin}}dV}{\int e_{\mathrm{kin}}dV}.\end{array}$$(C.1)

In the Carrington frame, a change in Ω_eff causes a Doppler shift in the frequency given by

$$\begin{array}{c}\hfill \omega =\omega ({\mathrm{\Omega }}_{\mathrm{eff}})+m({\mathrm{\Omega }}_{\mathrm{eff}}-{\mathrm{\Omega }}_0),\end{array}$$(C.2)

where ω is the frequency in the Carrington frame and ω(Ω_eff) is that in the frame rotating at Ω_eff. At high m, small differences in Ω_eff lead to significant differences in ω.

The top panel of Fig. C.1 shows that the effective rotation rates Ω_eff of the HFR modes diverge between the models with and without the RZ as m increases. The bottom panel of Fig. C.1 compares the frequencies of the HFR modes in the frame rotating at Ω_eff, as obtained from Eq. (C.2). It indicates that ω(Ω_eff) is indeed very similar for the HFR modes in the two models. The differences in ω(Ω_eff) between the two models are of the same order as those in Ω_eff. Therefore, it is demonstrated that the large deviation in the HFR mode frequencies seen in Fig. 4 is simply caused by the difference in the effective rotation rates.

Fig. C.1. Top: Effective rotation rate (Ωeff), given by Eq. (C.1), for the HFR modes in the models that include the RZ (solid line with points) and exclude the RZ (dashed line with open circles). Bottom: Frequency, ω(Ωeff), of the HFR modes in the frame rotating at Ωeff, calculated using Eq. (C.2), with the same notation.

Appendix D: Mode mass of the mixed modes

In this appendix, we present the mode masses of the computed eigenmodes (see e.g. Baudin et al. 2005). The mode mass of an eigenmode is defined as

$$\begin{array}{c}\hfill \mathcal{M}=\frac{1}{u_{\mathrm{RMS}}^2(r_s)}{\displaystyle {\int }_V}{\rho }_0{|\mathit{u}|}^{2}dV,\end{array}$$(D.1)

where u_RMS²(r_s) is the observed mode amplitude at the surface r_s. Note that the mode mass is independent of the choice of mode normalization. Modes with higher mode masses are energetically more difficult to excite (Samadi et al. 2001).

Figure D.1 presents the mode masses of the CZ–RZ mixed modes and their corresponding CZ inertial modes. It is shown that the mixed modes possess higher mode masses compared to the CZ inertial modes, and the mode mass tends to decrease as n_RZ increases. We note that, owing to very large velocity amplitudes in the RZ, mode masses of some mixed modes can exceed the solar mass M_⊙ by several orders of magnitude. These modes are very unlikely to be excited stochastically by turbulent convection in the CZ. Only mixed modes with sufficiently large n_RZ can have mode masses close to those of the CZ inertial modes. On the other hand, our study also reveals that the growth rate of the mixed modes decreases as the number of radial nodes, n_RZ, increases (see Fig. 6). These competing factors may determine the allowed range of n_RZ for the mixed modes that can be excited to observable amplitudes at the surface.

Fig. D.1. Mode mass (Eq. D.1) of the mixed modes (red stars) and of their nearest inertial modes with motions predominantly in the CZ (blue dots). The 11 mixed modes marked in Fig. 8 are shown. The colour of the star symbols denotes the number of radial nodes (nRZ) in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. For simplicity, we show only the mixed modes with nRZ ≤ 5.

Appendix E: Convergence in resolution

In this appendix, we present supplementary computations to demonstrate that both the complex eigenfrequencies and the eigenfunctions of the modes presented in this study are well resolved. To this end, we examine the convergence with resolution for two selected modes: the CZ–RZ mixed HFR mode (corresponding to the mode 5. (b) in Fig. 10), which has the largest radial variations (with n_RZ = 13), and the CZ–RZ mixed n_CZ = 1 Rossby mode (corresponding to the mode 4. (b) in Fig. 10), which has the largest latitudinal variations (with ℓ_RZ − m = 8). Figure E.1 shows that the eigenfunctions of these modes are almost unaltered for N_r ≥ 120 and N_θ ≥ 50, exhibiting radial and latitudinal resolution convergence for both modes. Tables E.1 and E.2 further show that the eigenfrequencies are also well converged within a difference of 0.1 nHz for N_r ≥ 120 and N_θ ≥ 50. This ensures that all other modes presented in this study, which have lower radial and latitudinal variations, are also converged.

Fig. E.1. Convergence with respect to grid resolution for two representative CZ–RZ mixed modes with strong spatial variations. The top and bottom rows show the results for the modes labelled 5. (b) and 4. (b) in Fig. 10, respectively. Panel 1: Mode structures for various radial resolutions (Nr) with fixed latitudinal resolution Nθ = 60. Panel 2: Mode structures for various Nθ with fixed Nr = 180. Surface and radial variations of the modes are represented in panels (a) and (b), respectively, in the same way as in Fig. 10.

All Tables

All Figures

Fig. 1. Profiles of differential rotation used in this study. (a) Observed 2D profile in a meridional plane, taken from Larson & Schou (2018). (b) Simplified differential rotation profile, given by Eq. (5), used to model the observed differential rotation. The dotted black lines indicate the base of the CZ at r = 0.71 R⊙. The Carrington rotation rate, Ω0/2π = 456 nHz, is marked on the colour bar. (c) Cuts of observed (solid) and simplified (dashed) differential rotation profiles at fixed latitudes: 0° (red), 30° (blue), 45° (brown), and 75° (dark green). The green-shaded area denotes the CZ, the grey-shaded region the OS, and the orange-shaded region the RZ.

In the text

Fig. 2. Left: Profiles of turbulent viscosity, ν(r), and turbulent thermal diffusivity, κ(r), as functions of radius. They are expressed by Eqs. (8) and (9). Right: Profile of super-adiabaticity, δ(r), as a function of radius (Eq. 10). The super-adiabaticity profile from the standard Model S is represented by the dashed red curve (Christensen-Dalsgaard et al. 1996). The y-axis is linear for |δ|< 10−7 and logarithmic beyond that. In both panels, the green-shaded area denotes the CZ, the grey-shaded region the OS, and the orange-shaded region the RZ.

In the text

Fig. 3. Comparison of selected inertial modes computed in the setups that include or exclude the RZ. (a) Velocity eigenfunctions in meridional cross-sections from the CZ-only model. The real part of the eigenfunction corresponds to a longitude ϕ0 (where uθ is maximum), while the imaginary part corresponds to the longitude ϕ0 − π/2m. The m = 1 high-latitude mode is normalized to have the maximum surface velocity of 10 m s−1, while the other modes are normalized to have a maximum surface velocity of 1 m s−1. Solid black curves denote the critical layers of the mode where ℜ[ω]=m(Ω − Ω0). The frequencies are measured in the Carrington frame. (b) Same as panel a but from the extended model that includes the RZ. Dotted black lines denote the base of the CZ. (c) Horizontal velocity eigenfunctions as functions of latitude at the surface. Dashed blue and solid red curves represent the results with and without the RZ, respectively. The correlation coefficient between the eigenfunctions, 𝒞, defined in Eq. (13), is noted above each subplot. (d) Radial profiles of the RMS velocity.

In the text

	Fig. 4. Top: Dispersion relations of the inertial modes in the CZ computed in the setups that include the RZ (solid lines with points) and exclude the RZ (dashed lines with open circles). The colours denote the various types of inertial modes. Bottom: Growth rates of the same modes, with the same notations.
In the text

Fig. 5. Dispersion relations of the Rossby modes inside the RZ with ℓ = m (navy blue), ℓ = m + 1 (red), and ℓ = m + 2 (dark green) for azimuthal orders 1 ≤ m ≤ 16. The frequencies are measured in the Carrington frame. Plus symbols (nRZ = 0) and open circles (nRZ = 4) denote modes with different numbers of radial nodes, nRZ, in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. Solid grey curves represent the theoretical dispersion relations of the classical Rossby modes ${\omega }_{\ell ,m}^{\mathrm{R}}$ given by Eq. (14). The dashed black line denotes the maximum possible value of ${\omega }_{\ell ,m}^{\mathrm{R}}$.

In the text

Fig. 6. Position of the RZ Rossby modes with ℓ = m = 3 in the complex frequency space and their eigenfunctions. (a) Complex frequency spectrum of equatorial Rossby modes in the RZ with azimuthal order m = 3. Star symbols denote the RZ Rossby modes, while the cross symbol denotes the CZ Rossby mode. The red (blue) star marks the RZ Rossby mode with the number of radial nodes, nRZ = 0 (nRZ = 3), in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. The vertical solid blue line denotes the observed Rossby mode frequency, while the vertical dashed green line denotes the frequency given by Eq. (14). (b) Meridional eigenfunction of the RZ Rossby mode with nRZ = 0 (red star). The longitudes corresponding to the real and imaginary components are designated in the same way as in Fig. 3. Dotted lines denote the position of the base of the CZ. Eigenfunctions are normalized such that the maximum of uθ is 1 m s−1. The unit of s1 is erg g−1 K−1. (c) Same as b but for the RZ Rossby mode with nRZ = 3 (blue star).

In the text

Fig. 7. Analyses of the non-sectoral Rossby mode in the RZ with ℓRZ = 4, m = 3, nRZ = 3. (a) Radial profiles of the RMS velocity in the horizontal (red) and radial (blue) directions. The profiles are normalized to have a maximum total RMS velocity of unity. (b) RMS of the radial component of various forces, defined in Eqs. (15)–(18). The pressure gradient force (solid blue), Coriolis force (dashed black), buoyancy force (dashed red), and viscous force (solid orange) are normalized so that the maximum of the pressure gradient force is unity. (c) Same as (b) but for the horizontal component of the forces. Colours of the shaded regions follow the same notation as in Fig. 2.

In the text

Fig. 8. Frequencies of the retrograde-propagating inertial modes in the CZ (dashed coloured; the same as in Fig. 4) overplotted with the dispersion relations of the RZ Rossby modes (solid grey; given by Eq. (14)). The left and right panels correspond to the modes with north–south symmetric and antisymmetric vorticity, respectively. Open circles mark the subset of RZ Rossby modes that can form mixed modes with the CZ inertial modes.

In the text

Fig. 9. Fraction of kinetic energy density, E = ⟨ρ0u2⟩, of the RZ Rossby modes contained in the CZ as a function of the absolute difference between the eigenfrequency of each mode and that of the nearest CZ inertial mode (see Fig. 8). For clarity, we show only the RZ Rossby modes with ℓRZ = m (circles), ℓRZ = m + 1 (diamonds), and ℓRZ = m + 2 (squares). The colour of the symbols denotes the number of radial nodes (nRZ) in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. Only modes with nRZ ≤ 5 are shown.

In the text

Fig. 10. Selected inertial modes computed in the extended setup, which includes the RZ. (a) Modes with motions predominantly in the CZ and decaying in the RZ. (b) Modes with comparable motions in the RZ and CZ, with flows in the CZ similar to those of the corresponding modes in panel a. (c) and (d): Surface and radial variations of the modes, respectively, as in Fig. 3. The spherical harmonic degree and the number of radial nodes of the RZ Rossby modes in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙ are denoted by ℓRZ and nRZ, respectively. The m = 1 high-latitude modes are normalized to have a maximum surface velocity of 10 m s−1, the m = 2 high-latitude modes to 4 m s−1, and the other modes to 1 m s−1. The notation and layout are the same as in Fig. 3.

In the text

Fig. 11. Frequencies of the observed solar inertial modes (open circles; from Löptien et al. 2018; Gizon et al. 2021; Hanson et al. 2022) overplotted with the dispersion relations of the RZ Rossby modes (solid grey; given by Eq. (14)). The left and right panels correspond to modes with north–south symmetric and antisymmetric vorticity, respectively. The observed inertial modes are high-latitude (blue), Rossby (red), HFR (green), and critical-layer (magenta) modes. Observed modes whose frequencies are within 10 nHz of an RZ Rossby mode frequency with the same symmetry (|Δω|≲10 nHz) are shown with a translucent fill.

In the text

	Fig. A.1. Radial profile of dissipation, Qvis, defined in Eq. (A.2), for the m = 3 prograde columnar mode in the models with (red) and without the RZ (blue). The eigenmodes are shown in Fig. 3. They are normalized to have the same total energy density in the CZ.
In the text

Fig. B.1. Analysis of the effects of differential rotation on sectoral Rossby modes with no nodes in radius (nRZ = nCZ = 0, ℓRZ = m). (a) Meridional cross-sections of uθ and uϕ of the Rossby modes computed under uniform rotation for different azimuthal orders, m. The longitudes corresponding to the real and imaginary components are designated in the same way as in Fig. 3. The frequencies and the growth rates measured in the Carrington frame are noted below the eigenfunctions. All eigenfunctions are normalized to have a maximum uθ of 1 m s−1 at the surface. The dotted lines indicate the position of the CZ base. (b) Same as (a) but under the simplified differential rotation given by Eq. (5), shown in Fig. 1. (c) Radial profiles of the RMS velocity of the modes computed with and without differential rotation. The colours of the shaded regions have the same meaning as in Fig. 2.

In the text

	Fig. C.1. Top: Effective rotation rate (Ωeff), given by Eq. (C.1), for the HFR modes in the models that include the RZ (solid line with points) and exclude the RZ (dashed line with open circles). Bottom: Frequency, ω(Ωeff), of the HFR modes in the frame rotating at Ωeff, calculated using Eq. (C.2), with the same notation.
In the text

	Fig. D.1. Mode mass (Eq. D.1) of the mixed modes (red stars) and of their nearest inertial modes with motions predominantly in the CZ (blue dots). The 11 mixed modes marked in Fig. 8 are shown. The colour of the star symbols denotes the number of radial nodes (nRZ) in the region 0.5 R⊙ ≤ r ≤ 0.71 R⊙. For simplicity, we show only the mixed modes with nRZ ≤ 5.
In the text

Fig. E.1. Convergence with respect to grid resolution for two representative CZ–RZ mixed modes with strong spatial variations. The top and bottom rows show the results for the modes labelled 5. (b) and 4. (b) in Fig. 10, respectively. Panel 1: Mode structures for various radial resolutions (Nr) with fixed latitudinal resolution Nθ = 60. Panel 2: Mode structures for various Nθ with fixed Nr = 180. Surface and radial variations of the modes are represented in panels (a) and (b), respectively, in the same way as in Fig. 10.

In the text