© The Authors 2025

1. Introduction

The Sun undergoes various oscillations through which we can probe deep into its interior. The pressure modes (p modes) are acoustic modes of oscillation restored by pressure and have a timescale of minutes (e.g. Christensen-Dalsgaard 2002). They have been successfully used to infer the Sun’s large-scale flows, internal differential rotation (e.g. Schou et al. 1998; Larson & Schou 2018) and meridional circulation (e.g. Giles et al. 1997; Gizon et al. 2020a). Modes in the inertial frequency range – ‘inertial modes’ for short – are low-frequency modes of oscillation restored by the Coriolis force (e.g. Greenspan 1968). Their periods are comparable to the Sun’s rotation period.

All observed inertial modes were first detected in the horizontal flow field on the Sun’s surface. They are most easily described as waves of radial vorticity that propagate retrograde in a frame that co-rotates with the Sun’s Carrington rotation rate (approximately the equatorial rotation rate). By convention, in this paper, we use the Carrington rotation rate as a reference when using the terms retrograde and prograde. The observed inertial modes belong to several classes (see Gizon et al. 2024, for a review). The equatorial Rossby modes follow the well-known classical dispersion relation for the sectoral Rossby modes of a uniformly rotating fluid for azimuthal wavenumbers in the range 3 ≤ m ≲ 30 (Löptien et al. 2018; Liang et al. 2019; Proxauf et al. 2020; Hanson & Hanasoge 2024). Their retrograde propagation results from the planetary β effect, which arises due to the latitudinal variation of the tangential component of the Coriolis force (e.g. Vallis 2017). Gizon et al. (2021) found many additional quasi-toroidal global inertial modes of oscillation in the solar data. Among these modes are high-latitude modes with m = 1, 2, 3. These global modes can have either symmetric or antisymmetric north-south vorticity. Their velocity eigenfunctions display a characteristic spiral pattern in the polar regions. The mode with the largest velocity amplitude (∼15 m/s) is the m = 1 high-latitude mode with north-south symmetric radial vorticity; it is present in daily-cadence solar Dopplergrams over the last five solar cycles (Liang & Gizon 2025). We note that the velocity pattern associated with this m = 1 mode had been reported before but misidentified as giant convection cells (Hathaway et al. 2013). In addition to high-latitude modes, many other modes, with m up to at least 10, have been observed with amplitudes that peak at mid-latitudes, near their critical latitudes. The critical latitude of a mode is the latitude where its phase speed equals the local rotation velocity. Modes that are retrograde (with respect to the equatorial rotation rate) can have critical latitudes since the Sun’s rotation rate decreases with latitude (with the sharpest decrease above 60°). All the modes described above, including the equatorial Rossby and high-latitude modes, have critical latitudes.

Furthermore, Hanson et al. (2022) detected retrograde modes with equatorially-antisymmetric radial vorticity. These modes do not follow the dispersion relation of any classical Rossby modes. Since they propagate with phase speeds about three times that of the equatorial Rossby modes, Hanson et al. (2022) referred to them as high-frequency retrograde (HFR) modes. Many HFR modes also have critical latitudes, as they propagate in the retrograde direction.

Under the assumption that the modes have small enough amplitudes, linear modelling is key to identifying the modes and revealing their physical nature. It was recognized early on that solar latitudinal differential rotation and turbulent viscosity are crucial ingredients to describe the purely toroidal modes on the sphere (Gizon et al. 2020b; Fournier et al. 2022). The frequency spectrum under differential rotation differs significantly from that under uniform rotation. Bekki et al. (2022b) also solved the eigenvalue problem in a spherical shell representative of the differentially rotating convection zone. The linear modes in 1.5D of Fournier et al. (2022) and 2.5D of Bekki et al. (2022b) were both essential in identifying the observed modes and providing a unified explanation for the equatorial Rossby modes and the mid- and high-latitude modes, including the m = 1 high-latitude mode (see Gizon et al. 2021, 2024). It was also found that the high-latitude modes are baroclinically unstable and sensitive to the latitudinal entropy gradient in the convection zone. Linear modelling has also been instrumental in identifying the HFR modes as non-toroidal retrograde mixed inertial modes with anti-symmetric radial vorticity (Triana et al. 2022; Bhattacharya & Hanasoge 2023; Bekki 2024; Jain et al. 2024).

Furthermore, some linear studies have described retrograde modes that are similar to the equatorial Rossby modes but with n = 1, where n is the number of nodes in radius (Bekki et al. 2022b; Bhattacharya & Hanasoge 2023; Jain et al. 2024). These modes have been referred to as n = 1 equatorial Rossby modes (Bekki et al. 2022b). They also belong to the class of retrograde mixed modes, similar to the HFR modes, but with north-south symmetric radial vorticity (Bekki et al. 2022b; Jain et al. 2024). In order to label this particular subset of mixed modes, we stick to the name n = 1 equatorial Rossby mode. These modes have frequencies comparable to the sectoral power of radial vorticity in the solar surface flows for m ≳ 8 (Hanson & Hanasoge 2024). However, the observed power cannot be conclusively attributed to these modes due to the lack of comparison between their eigenfunctions and the solar flow structures. Adding to that, the frequencies and growth rates of these linear modes are sensitive to the superadiabaticity, which is not well-constrained for the Sun (Hanson & Hanasoge 2024; Bekki et al. 2022b).

Linear calculations previously predicted modes that are excited by thermal instability in a rotating fluid (Roberts 1968; Busse 1970). These are non-toroidal modes consisting of convective columns that propagate in the prograde direction under the effect of compressional and topographical β effects (e.g. Glatzmaier & Gilman 1981; Gastine et al. 2014; Bekki et al. 2022b; Hindman & Jain 2022. These modes are also referred to as columnar convective modes, Busse columns, banana cells, and thermal Rossby modes. In our paper, we refer to them as prograde columnar modes. In non-linear numerical simulations, prograde columnar modes play a crucial role in the dynamics of the convection zone as they transport angular momentum equatorward and contribute to determining the differential rotation profile (e.g. Miesch et al. 2006). They also transport heat poleward, which helps establish the thermal wind balance (e.g. Matilsky et al. 2020). Despite their important function in simulations, they have not yet been identified in the solar surface flows. There are several possible reasons why they have not been observed. They may be present below the near-surface layers but shielded by small-scale convection close to the solar surface, or their spatial scale may be too small to be identified in the solar surface flows (Hindman & Jain 2022; Featherstone & Hindman 2016). Also, the modes may simply not be present at measurable amplitudes in the solar convection zone, which could be related to the unresolved convective conundrum (e.g. Hotta et al. 2023).

In our paper, we study the above-discussed modes as summarized in Table 1. Understanding the physics of inertial modes can help diagnose different properties of the solar interior to complement p mode helioseismology, as they are sensitive to several parameters of the deep convection zone, including the turbulent viscosity, superadiabaticity, and latitudinal entropy gradient (Gizon et al. 2021; Bekki et al. 2022b; Hanson & Hanasoge 2024). More importantly, the baroclinically unstable high-latitude modes likely play a key role in the dynamics of the Sun by controlling the solar differential rotation (Bekki et al. 2024).

2. Previous theoretical studies: Various assorted assumptions

Rieutord & Valdettaro (1997) numerically analysed the spectrum of inertial waves trapped in a spherical shell containing a uniformly rotating incompressible fluid. Baruteau & Rieutord (2013) extended this setup without stratification to a differentially rotating fluid. Following the discovery of solar quasi-toroidal inertial modes (Löptien et al. 2018; Gizon et al. 2021), several attempts have been made to model the various classes of solar inertial modes employing varying simplifying assumptions. Bekki et al. (2022b) modelled the observed solar inertial modes using a differentially rotating compressible fluid in the solar convection zone. Furthermore, Bhattacharya & Hanasoge (2023), Bhattacharya et al. (2024) modelled the solar inertial modes in the spherical shell geometry considering the anelastic approximation (Glatzmaier & Gilman 1981). This approximation assumes that the sound speed is much faster than the other characteristic flow speeds in the system, thereby filtering out the acoustic modes. Although these sound-proof setups often do not conserve energy for gravity waves in the sub-adiabatic radiative interior, they generally conserve energy in the nearly adiabatic convection zone (Brown et al. 2012). Triana et al. (2022) also used a simpler uniformly rotating unstratified incompressible fluid in a spherical shell to model the solar inertial modes. The study aimed to identify the HFR modes previously discovered by Hanson et al. (2022). Such an incompressible model is widely used to study planetary eigenmodes employing the Boussinesq approximation to include the buoyancy effects (e.g. Marti et al. 2016; Barik et al. 2023). Furthermore, several studies have been conducted in simplified geometries, such as on spherical surfaces (Fournier et al. 2022; Gizon et al. 2020b) and in cylindrical geometry (Jain & Hindman 2023; Jain et al. 2024). We refer to the model of Jain & Hindman (2023) as sound-proof because it assumes the limit of low frequency and high sound speed to filter out the sound waves.

The setups in these linear models for the solar inertial modes differ in various aspects, such as the numerical domain, the degree of superadiabaticity, the choice of background rotation profile, and the latitudinal entropy gradient in the convection zone. Many studies, including Triana et al. (2022), Bhattacharya & Hanasoge (2023), Jain & Hindman (2023), and Jain et al. (2024) assumed solid-body rotation for the Sun. However, as discussed earlier, latitudinal differential rotation implies the existence of critical latitudes and additional families of modes (e.g. the high-latitude modes). Table 2 summarizes the assorted fundamental model assumptions in the various linear simulations of the solar inertial modes. Although the solar inertial modes have also been studied using nonlinear simulations (e.g. Bekki et al. 2022a, 2024; Blume et al. 2024), we limit our studies to the linear modes in the solar convection zone.

Despite employing several differing simplifications, the effects of these simplifying assumptions on the linear inertial modes have not been tested. Bekki (2024) compared only the HFR modes obtained from the different models used in the literature employing spherical shell geometry. The effects of various sound-proof approximations were tested earlier for gravity waves in stellar interiors (Brown et al. 2012; Vasil et al. 2013), but not for inertial modes.

3. Description of our numerical computations: Compressible, anelastic, and Boussinesq cases

In this work, we study the effects of the anelastic and Boussinesq approximations on the solar inertial modes by comparison to the fully compressible case used as the reference. We implement the compressible, anelastic, and Boussinesq model setups in Dedalus (Burns et al. 2020) with the same boundary conditions (stress-free impenetrable). The background stratification is the same for the compressible and the anelastic cases. The background temperature, superadiabaticity, and latitudinal entropy variations are identical in all computations, including the Boussinesq case. This helps us understand the main physical ingredients needed to describe the inertial modes correctly. We note that we always represent the effects of small-scale turbulence in terms of turbulent eddy viscosity. The effects of magnetic fields on the inertial modes are ignored in our study. We discuss the five types of inertial modes described in Table 1.

3.1. Linearized equations of fluid dynamics

The linearized Navier-Stokes equation in the Carrington frame, which rotates at Ω₀/2π = 456 nHz, can be expressed as

$$\begin{array}{cc}& \frac{\partial \mathbf{u}}{\partial t}+{\mathbf{v}}_0·\mathrm{\nabla }\mathbf{u}+\mathbf{u}·\mathrm{\nabla }{\mathbf{v}}_0+2{\mathrm{\Omega }}_0{\mathbf{e}}_z\times \mathbf{u}\hfill \\ \hfill & +\mathrm{\nabla }\left(\frac{p_1}{{\rho }_0}\right)-\frac{s_1}{c_p}g{\mathbf{e}}_r-\frac{1}{{\rho }_0}\mathrm{\nabla }·\overleftrightarrow{D}=0.\hfill \end{array}$$(1)

Here, u, ρ₁, p₁, and s₁ denote the perturbations of velocity, density, pressure, and entropy with respect to the background, while c_p is the specific heat at constant pressure. The viscous stress tensor is given by $\overleftrightarrow{D}={\rho }_0\nu (\mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }{\mathbf{u}}^T-\frac{2}{3}(\mathrm{\nabla }·\mathbf{u})\overleftrightarrow{I}),$ with ν = 10¹² cm²/s being the spatially constant isotropic turbulent viscosity and ⃡I the identity tensor. The background flow v₀ = (Ω(r, θ)−Ω₀) r sin θ e_ϕ represents the solar differential rotation Ω(r, θ) in the convection zone taken from Larson & Schou (2018). In this study, we omit the effects of meridional circulation, which has been shown to have a small impact (Gizon et al. 2020b; Fournier et al. 2022). For compressible and anelastic setups, the background stratifications of density ρ₀, pressure p₀, and temperature T₀ are taken from the standard solar model S (Christensen-Dalsgaard et al. 1996). The acceleration due to gravity g is determined using the hydrostatic balance in model S. For the anelastic setup, Eq. (1) is the form of the Navier-Stokes equation that is recommended by Brown et al. (2012), as it conserves energy both in the convection zone and the radiative zone. For the incompressible setup, we apply the Boussinesq approximation following Spiegel & Veronis (1960) to account for buoyancy effects, ensuring similar conditions to the other models and allowing baroclinic instability (e.g. Molemaker et al. 2005). This is in contrast to the model used in Triana et al. (2022), which ignored the buoyancy effects. Hence, we refer to it as the Boussinesq model instead of the incompressible model in our study. Under this approximation, we use the mean density, ρ_m = 0.0589 g cm⁻³, in the convection zone as the background density ρ₀. Figure 1 depicts the difference between the mean density and the background density from model S. However, the temperature T₀ is used from the Solar model S. The background pressure p₀ is given by the hydrostatic equilibrium, similar to Spiegel & Veronis (1960). See Appendix A for more details.

Fig. 1. Variation of the background density in the solar convection zone from the standard model S (Christensen-Dalsgaard et al. 1996), denoted by the blue curve. The black dashed line shows the mean density ρm in the solar convection zone. The Boussinesq setup uses the mean density, whereas the other stratified setups use the density from model S.

The following linearized equation of entropy (s) is used

$$\begin{array}{c}\hfill \frac{\partial s_1}{\partial t}+{\mathbf{v}}_0·\mathrm{\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}\mathrm{\nabla }·(\kappa {\rho }_0{T}_0\mathrm{\nabla }{s}_1)=0,\end{array}$$(2)

where s₀ denotes the background entropy. In this study, we assume the adiabatic stratification in radius, ∂s₀/∂r = 0, for simplicity. The latitudinal entropy gradient ∂s₀/∂θ represents the thermal wind balance of the solar differential rotation and is estimated as (e.g. Pedlosky 1982; Thompson et al. 2003)

$$\begin{array}{c}\hfill \frac{\partial s_0}{\partial \theta }=r^{2}sin\theta \frac{g}{c_p}\frac{\partial ({\mathrm{\Omega }}^2)}{\partial z},\end{array}$$(3)

where z is the coordinate along the rotational axis. The turbulent thermal diffusivity κ is assumed to be spatially constant, with a value equal to ν.

For a compressible fluid, the linearized continuity equation

$$\begin{array}{c}\hfill \frac{\partial {\rho }_1}{\partial t}+\mathrm{\nabla }·({\rho }_0\mathbf{u})+\mathrm{\nabla }·({\rho }_1{\mathbf{v}}_0)=0,(\mathrm{compressible})\end{array}$$(4)

is used in combination with the linearized equation of state

$$\begin{array}{c}\hfill \frac{p_1}{p_0}=\gamma \frac{{\rho }_1}{{\rho }_0}+\frac{s_1}{c_v},\end{array}$$(5)

where γ = 5/3 is the specific heat ratio, and c_v is the specific heat at constant volume. Under the anelastic approximation, we assume that the sound speed is much faster than the rotational or advective flow speed in the Sun. Thus, the continuity equation is reduced to the anelastic equation,

$$\begin{array}{c}\hfill \mathrm{\nabla }·({\rho }_0\mathbf{u})=0.(\mathrm{anelastic})\end{array}$$(6)

Under the Boussinesq approximation, this is further reduced to

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{u}=0.(\mathrm{Boussinesq})\end{array}$$(7)

The above linearized fluid dynamical equations are used in this work to model and understand the physics of solar inertial modes.

3.2. Boundary conditions

The Sun does not have any hard boundaries, but to execute the computations, one needs to implement some boundary conditions. As in some previous linear studies, we perform our computations in the spherical shell of the solar convection zone. The numerical domain in our study extends from the base of the convection zone (r_i = 0.71R_⊙) to slightly below the photosphere (r_o = 0.985R_⊙). There are various possible options for the boundary conditions, such as free-surface boundary conditions, stress-free boundary conditions, and no-slip boundary conditions. It is not very clear which boundary condition would be the best for modelling the solar inertial modes. The no-slip boundary condition is not expected to be good in modelling the inertial modes as it does not allow any motions at the boundaries. The free-surface boundary condition might be adequate to model the inertial modes as it allows a flexible boundary height. However, to avoid the computational complexity, we use the stress-free impenetrable boundary conditions, which have been used in previous studies of the solar inertial modes (e.g. Bekki et al. 2022b; Bhattacharya et al. 2024; Triana et al. 2022). This works better than the no-slip boundary conditions as it allows horizontal motions at the boundaries, although it assumes no vertical motions at the boundaries. Also, we assume no flux of entropy across each boundary. The boundary conditions can be expressed as

$$\begin{array}{cc}& u_r(r_i)=u_r(r_o)=0,\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\overleftrightarrow{D}}_{r\theta }(r_i)={\overleftrightarrow{D}}_{r\theta }(r_o)=0,\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\overleftrightarrow{D}}_{r\varphi }(r_i)={\overleftrightarrow{D}}_{r\varphi }(r_o)=0,\hfill \end{array}$$(10)

$$\begin{array}{cc}& {\frac{\partial s_1}{\partial r}|}_{r_i}={\frac{\partial s_1}{\partial r}|}_{r_o}=0.\hfill \end{array}$$(11)

Since we use spherical harmonics as basis functions, the problem has no singularities at the poles and does not require any boundary conditions there (Boyd 2001).

3.3. Formulation of eigenvalue problem in Dedalus

We use Dedalus (Burns et al. 2020), a flexible open-source spectral code, to solve the linear eigenvalue problem of the solar inertial modes. The above Eqs. (1), (2), and (4) – (7) are solved as an eigenvalue problem using the wave Ansatz where each perturbed physical quantity is proportional to exp(imϕ − iωt), with m the azimuthal order and ω the frequency. We implement Eqs. (8) – (11) as the boundary conditions. The calculations are performed on the spherical shell basis of Dedalus using the sparse eigenvalue solver (Burns et al. 2020). We solve the sparse problem spanning the range of the inertial frequencies, |ℜ[ω]| ≤ 2Ω₀. This helps us filter out the low-frequency modes of our interest. The obtained eigenfrequencies ω and eigenfunctions of the inertial modes are analysed for different model assumptions. All the eigenmodes we present in this paper are obtained using a grid on the spherical shell basis with 24 radial and 124 latitudinal points. We detect them based on their north-south symmetries, growth rates, and other known properties, such as the number of radial or latitudinal nodes. The convergence errors in the eigenfrequencies obtained by halving or doubling the resolutions are orders of magnitude smaller than the frequency resolution of observations (which is on the order of a few nHz, see Gizon et al. 2021). The reported eigenfunctions are also well-converged.

4. Effects of simplifying the continuity equation

4.1. Effects of anelastic and Boussinesq approximations under uniform rotation

We start by comparing the eigenmodes obtained from the different models under uniform rotation. Figure 2 compares the eigenfunctions of the n = 0 equatorial Rossby mode with m = 3, the HFR mode with m = 10, the n = 1 equatorial Rossby mode with m = 3, and the prograde columnar mode with m = 3 for the three different models. First, we note that the eigenmodes in the compressible and anelastic models are nearly identical. They have negligible differences in the eigenfunctions of less than 1%, and their frequencies are the same correct up to 0.1 nHz. However, the eigenmodes in the Boussinesq model differ significantly from those obtained using the other models, except for the n = 0 equatorial Rossby modes. In all the models, the n = 0 equatorial Rossby modes are almost perfectly toroidal and are driven by the planetary β effect, which is unaffected by the background density stratification. However, the other modes, having a substantial non-toroidal nature, are quite different in the Boussinesq model as compared to the other models. Significant differences observed in the spatial eigenfunctions of the HFR modes, the n = 1 equatorial Rossby modes, and the prograde columnar modes can be understood as follows: In contrast to the n = 0 equatorial Rossby modes, these non-toroidal modes exhibit a radial node in the horizontal velocity eigenfunctions within the convection zone (CZ). Under the background density stratification, the constraint of local mass conservation (Eq. 6) requires the mass fluxes in the lower and upper parts of CZ (across the nodal plane) to be balanced with each other, leading to slower velocities in the lower CZ. In the Boussinesq model, on the other hand, the velocity eigenfunctions can have comparable amplitudes throughout the CZ.

Fig. 2. Comparison of the eigenmodes of the different classes of inertial modes (see Table 1) computed using compressible, anelastic, and Boussinesq models under uniform rotation. Note that the unstable high-latitude inertial modes are absent under solid body rotation. Here, we plot the real part of uθ and the imaginary part of uϕ of the computed eigenmodes. The longitudes corresponding to the real and imaginary phases of the eigenfunctions are ϕ = ϕ0 and ϕ = ϕ0 − π/2m, respectively, where ϕ0 is a longitude, where uθ attains its maximum. The corresponding frequencies measured in the Carrington frame are stated below each eigenmode. The imaginary parts of the frequencies indicate the growth rates of the modes. All eigenfunctions are normalized such that the maximum of uθ is 1 m/s at the surface.

Figure 3 presents the dispersion relations of the modes displayed in Fig. 2 for the three models under uniform rotation, covering azimuthal orders m from 0 to 15. Similar to the eigenmodes, the dispersion relations of the various modes are identical between the anelastic and compressible models across all values of m. The deviations in the dispersion relations of the Boussinesq model relative to the other models are evident for the non-toroidal modes discussed above. The deviations are enormous for the prograde columnar modes and the HFR modes. In both cases, the frequencies are significantly shifted to the negative (more retrograde) direction. Due to their less toroidal nature, the deviations are smaller for the n = 1 equatorial Rossby modes. The purely toroidal n = 0 equatorial Rossby modes have the same dispersion relations for all the models.

Fig. 3. Dispersion relations of the studied classes of inertial modes under uniform rotation for the range of azimuthal orders 0 ≤ m ≤ 15. The blue, green, orange, and red colours represent the n = 0 equatorial Rossby, the HFR, the n = 1 equatorial Rossby, and the prograde columnar modes, respectively. Open circles, cross symbols, and dotted points show the results computed from the fully compressible, anelastic, and Boussinesq models, respectively.

To investigate the effects of using the Boussinesq approximation on mode frequencies, we examine the linearized vorticity equation. By taking a curl of Eq. (1) under uniform rotation, we express the radial and z components of the vorticity equation as

$$\begin{array}{cc}& \frac{\partial {\zeta }_r}{\partial t}\approx \underset{\mathrm{planetary}\beta \mathrm{effect}}{\underbrace{\frac{2{\mathrm{\Omega }}_{0}sin\theta }{r}{u}_{\theta }}}-\frac{2{\mathrm{\Omega }}_{0}cos\theta }{H_{\rho }}{u}_r\hfill \\ \hfill & +2{\mathrm{\Omega }}_{0}cos\theta \frac{\partial u_r}{\partial r}-\frac{2{\mathrm{\Omega }}_{0}sin\theta }{r}\frac{\partial u_r}{\partial \theta },\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill & \frac{\partial {\zeta }_z}{\partial t}\approx \underset{\mathrm{compressional}\beta \mathrm{effect}}{\underbrace{-\frac{2{\mathrm{\Omega }}_0}{H_{\rho }}{u}_r}}+\underset{\mathrm{topographic}\beta \mathrm{effect}}{\underbrace{2{\mathrm{\Omega }}_0\frac{\partial u_z}{\partial z}}}-\frac{g}{c_{p}r}\frac{\partial s_1}{\partial \varphi },\hfill \end{array}$$(13)

where ζ = ∇ × u and H_ρ is the density scale height. For simplicity, we omit the viscous diffusive terms. Here, the first term on the right-hand side of Eq. (12) represents the planetary β effect, while the first term on the right-hand side of the Eq. (13) represents the compressional β effect. The second term on the right-hand side of the Eq. (13) corresponds to the topographic β effect (when integrated over z). To assess the relative importance of these β effects on the mode frequencies, we further transform the above equations into the following form:

$$\begin{array}{cc}\hfill \omega |{\zeta }_r{|}^2& \approx {\mathcal{G}}_{\mathrm{planetary}}+{\mathcal{G}}_{\mathrm{comp}}+{\mathcal{G}}_{\mathrm{other}},\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill \omega |{\zeta }_z{|}^2& \approx {\mathcal{W}}_{\mathrm{comp}}+{\mathcal{W}}_{\mathrm{topographic}}+{\mathcal{W}}_{\mathrm{other}},\hfill \end{array}$$(15)

with

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\mathrm{planetary}}& =i\frac{2{\mathrm{\Omega }}_{0}sin\theta }{r}{u}_{\theta }{\zeta }_r^{\ast },\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\mathrm{comp}}& =-i\frac{2{\mathrm{\Omega }}_{0}cos\theta }{H_{\rho }}{u}_r{\zeta }_r^{\ast },\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\mathrm{other}}& =2i{\mathrm{\Omega }}_{0}cos\theta \frac{\partial u_r}{\partial r}{\zeta }_r^{\ast }-i\frac{2{\mathrm{\Omega }}_{0}sin\theta }{r}\frac{\partial u_r}{\partial \theta }{\zeta }_r^{\ast },\hfill \end{array}$$(18)

and

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\mathrm{comp}}& =-i\frac{2{\mathrm{\Omega }}_0}{H_{\rho }}{u}_r{\zeta }_z^{\ast },\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\mathrm{topographic}}& =2i{\mathrm{\Omega }}_0\frac{\partial u_z}{\partial z}{\zeta }_z^{\ast },\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\mathrm{other}}& =\frac{mg}{c_{p}r}{s}_1{\zeta }_z^{\ast }.\hfill \end{array}$$(21)

Here, ^* denotes the complex conjugate of the quantities. A negative real part of 𝒢 or 𝒲 indicates that the associated physical effect promotes retrograde propagation, whereas a positive real part implies that it promotes prograde propagation. Figures B.1–B.4 present all terms in Eqs. (14) and (15) for the n = 0 equatorial Rossby mode (m = 3), the HFR mode (m = 10), the n = 1 equatorial Rossby mode (m = 3), and the prograde columnar mode (m = 3) from compressible and Boussinesq setups. Note that 𝒢_comp = 𝒲_comp = 0 in the Boussinesq setup because the assumption of constant density leads to an infinite density scale height H_ρ.

Our analysis reveals that the primary cause of frequency changes from the anelastic (or compressible) to Boussinesq models is the absence of 𝒲_comp in the z-vorticity equation. Figure 4 presents ℜ[𝒲_comp] in the compressible setup for the HFR mode with m = 10, n = 1 equatorial Rossby mode with m = 3, and the prograde columnar mode with m = 3. As for the n = 1 equatorial Rossby modes, ℜ[𝒲_comp] is negative, suggesting that the compressional β effect promotes their retrograde propagation. This occurs because, in these modes, the radial vortical motions are dominant and z-vorticity ζ_z is primarily generated by the strong radial shear of longitudinal flows and is not associated with their radial motions (ζ_z and −∂u_r/∂ϕ have the opposite sign). The absence of 𝒲_comp in the Boussinesq model slightly shifts their mode frequencies towards the positive (more prograde) direction, as seen in Fig. 3 for m ≤ 4. On the other hand, both HFR modes and prograde columnar modes have strongly positive 𝒲_comp, indicating that the compressional β effect enforces their prograde propagation. This is because their radial motions are strongly associated with their z-vortices (ζ_z and −∂u_r/∂ϕ have the same sign). Without the compressional β effect, the frequencies of the HFR and prograde columnar modes shift towards the negative (more retrograde) direction, as seen in Fig. 3. We also note that, in prograde columnar modes, the planetary β effect plays an additional role in decreasing their prograde frequencies in the Boussinesq model (see Fig. B.4). In the absence of density stratification, upflows converge towards the equator due to spherical curvature, enhancing the planetary β effect. This does not occur in compressible or anelastic models because the upflows tend to expand and drive horizontally diverging motions due to the background density stratification.

Fig. 4. Estimates of the importance of the compressional β effect under uniform rotation using ℜ[𝒲comp], as defined in Eq. (19), for the compressible setup. A negative value of ℜ[𝒲comp] implies that the compressional β effect promotes retrograde propagation, while a positive value implies that it promotes prograde propagation. Column (a) denotes the HFR mode with m = 10, column (b) denotes the n = 1 equatorial Rossby mode with m = 3, and column (c) denotes the prograde columnar mode with m = 3. For all modes, we normalize ℜ[𝒲comp] such that its maximum absolute value is 1.

4.2. Effects of anelastic and Boussinesq approximations under solar differential rotation and latitudinal entropy gradient

Now, we include solar differential rotation and its associated latitudinal entropy gradient in our computations. Figure 5 presents the dispersion relations of the different modes with 0 ≤ m ≤ 15 under this setup for all models. Here, we also consider the high-latitude modes, which are baroclinically unstable and owe their existence to the Sun’s differential rotation and the latitudinal entropy gradient (Bekki et al. 2022b). For simplicity, we only study the fastest-growing high-latitude mode for each azimuthal order. As in the uniform rotation case, the compressible and anelastic models yield the same dispersion relations for the different modes. In the Boussinesq model, the frequencies of the prograde columnar modes and the HFR modes shift in the negative (more retrograde) direction. In contrast, the frequencies of the n = 0 equatorial Rossby and n = 1 equatorial Rossby modes experience only marginal changes in the Boussinesq model. These results are consistent with the uniform rotation case (Sect. 4.1). The retrograde-propagating high-latitude modes exhibit a trend of negative frequency shift in the Boussinesq model, similar to that of the HFR modes. This arises from the absence of the compressional β effect in the Boussinesq model, as high-latitude modes are essentially non-toroidal.

Fig. 5. Dispersion relations of the various classes of inertial modes listed in Table 1, computed using different models with solar differential rotation and the associated latitudinal entropy gradient. The same notation applies to colours and symbols as in Fig. 3. The only addendum is the high-latitude mode with north-south symmetric radial vorticity, represented in purple.

Figure 6 compares the velocity eigenfunctions of the n = 0 equatorial Rossby mode with m = 3, the high-latitude mode with m = 1, the HFR mode with m = 10, the n = 1 equatorial Rossby mode with m = 3, and the prograde columnar mode with m = 3 in the differentially rotating case across the three models. As before, the eigenmodes obtained from the compressible and anelastic models are identical. In contrast to the uniform rotation case, the eigenfunctions in the Boussinesq model differ even for the n = 0 equatorial Rossby modes under solar differential rotation. We note that under differential rotation ΔΩ(r, θ) = Ω(r, θ)−Ω₀, the modes have critical latitudes θ_c, where their phase speeds match the local differential rotation speed (i.e. ΔΩ(r, θ_c) = ℜ[ω]/m). Black solid curves in Fig. 6 indicate the locations of the critical latitudes. Previous studies (Gizon et al. 2020b; Fournier et al. 2022) have shown that the mode eigenfunctions undergo strong distortion near the critical latitudes. n = 0 equatorial Rossby, HFR, and n = 1 equatorial Rossby modes are confined to the equatorial region by the critical latitudes. In contrast, the high-latitude modes exist at latitudes above the critical latitudes. The columnar modes propagate prograde and lack critical latitudes.

Fig. 6. Comparison of the eigenmodes for different classes of inertial modes computed using compressible, anelastic, and Boussinesq models under solar differential rotation and the associated latitudinal entropy gradient. The different rows denote: a) n = 0 equatorial Rossby mode (m = 3), b) high-latitude mode with north-south symmetric radial vorticity (m = 1), c) HFR mode (m = 10), d) n = 1 equatorial Rossby mode (m = 3), e) prograde columnar mode (m = 3). Here, we plot the real part of uθ and the imaginary part of uϕ of the computed eigenmodes in the meridional plane. The corresponding longitudes are chosen in the same way as in Fig. 2. The black solid curves on the meridional cross-sections denote the critical latitudes where ℜ[ω]=m(Ω − Ω0). We also show the surface velocity: uϕ for the high-latitude mode and the prograde columnar mode, and uθ for the other modes. The corresponding frequencies in the Carrington frame are stated below each eigenmode. The imaginary parts of the frequencies are the growth rates of the modes. All eigenfunctions, except for the high-latitude mode, are normalized such that the maximum of uθ is 1 m/s at the surface. For the high-latitude mode, the maximum velocity at the surface is set to 10 m/s.

Under differential rotation, substantial radial motions develop near the critical latitudes (Bekki et al. 2022b). To estimate the relative impact of the radial motions on the eigenmodes, we define the non-toroidicity Γ based on the ratio of the kinetic energy in the radial direction to the total kinetic energy of the mode:

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\sqrt{\frac{{\int }_{\mathit{CZ}}{\rho }_0{u}_r^{2}dV}{{\int }_{\mathit{CZ}}{\rho }_0(u_r^2+u_{\theta }^2+u_{\varphi }^2)dV}}.\end{array}$$(22)

Figure 7 presents the non-toroidicity Γ of the various modes for different azimuthal orders m. Compared to uniform rotation, the non-toroidicity Γ increases significantly due to critical latitudes under differential rotation. This is most significant for the n = 0 equatorial Rossby modes at m ≥ 4 but is also seen in the HFR modes (m ≥ 4) and the n = 1 equatorial Rossby modes (m ≥ 3). For prograde columnar and high-latitude modes, the non-toroidicity Γ increases almost linearly with m. Unlike the prograde columnar modes, the retrograde high-latitude modes are strongly influenced by the critical latitudes. In the Boussinesq model, these modes exhibit more retrograde frequencies due to the absence of the compressional β effect. Consequently, their critical latitudes shift polewards, confining the mode power to higher latitudes than in the compressible model (see Fig. 6).

Fig. 7. Estimates of the significance of the radial motions associated with the different inertial modes for the compressible setup, quantified using non-toroidicity Γ defined in Eq. (22). The values of non-toroidicity for different modes are plotted against azimuthal orders ranging from m = 1 to m = 10. Different colours represent different modes, as in Fig. 5. Solid lines denote modes under differential rotation, while dotted lines denote modes under uniform rotation. Open circles indicate modes affected by critical latitudes under differential rotation, while filled circles indicate other modes under differential rotation.

5. Summary and outlook

In this work, we use a fixed framework in Dedalus to assess the validity of the anelastic and Boussinesq approximations for modelling the inertial modes in the solar convection zone. Firstly, we find that the eigenmodes obtained by the compressible and anelastic models are almost identical, regardless of the presence of the background solar differential rotation. The anelastic assumption does not affect the properties of the solar inertial modes since a scale separation between low-frequency inertial modes and very high-frequency sound waves occurs. Hence, one can safely use the anelastic approximation to simplify the calculations and reduce the numerical cost. Further, our results show that most of the inertial modes computed using the Boussinesq model are appreciably different from those calculated using the other models. The Boussinesq model reproduces identical inertial modes as the other models only for the n = 0 equatorial Rossby modes under uniform rotation, which are purely toroidal. Otherwise, the absence of density stratification in the Boussinesq model leads to significant differences in non-toroidal modes, as it lacks the compressional β effect. Therefore, a permissible setup for the studied solar inertial modes is the anelastic model with solar differential rotation because of the substantial effects of differential rotation on the eigenmodes.

There are still some issues remaining in the current modelling of the inertial modes. One prominent issue is the absence of the radiative interior below and the near-surface layer above the computational domain. Whether the anelastic formulation is valid for modelling the inertial modes in the radiative zone remains an open problem. The anelastic models have been found to fail to conserve energy for the gravity waves in the sub-adiabatic radiative interior, and pseudo-incompressible models have been shown to work better in such cases (Brown et al. 2012; Vasil et al. 2013). Thus, the anelastic approximation can likely affect the non-toroidal inertial modes in the radiative interior if they couple with gravity modes (gravito-inertial modes, see, e.g. Mathis 2009; Dintrans & Rieutord 2000). Nonetheless, the anelastic formulation we use in our paper is expected to be valid for studying quasi-toroidal inertial modes in the radiative interior, as done through nonlinear numerical simulations (Blume et al. 2024). On the other hand, it is well known that the anelastic approximation breaks down when modelling the solar convection in the near-surface layer where the stratification is strongly superadiabatic and the convective speed becomes as high as the sound speed (e.g. Nordlund et al. 2009). The extent to which the anelastic formulation can be used to model the low-frequency inertial modes in this near-surface layer requires future verification.

Acknowledgments

Author contributions: LG and XZ initiated this project, SM implemented all the equations in Dedalus and performed the computations, YB provided close 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 R. Cameron, P. Dey, V. Kannan, and J. Schou for helpful discussions. We also acknowledge the hospitality of Nordita, Stockholm, during the 2024 program on stellar convection. 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. The codes and data used in the manuscript are available in the Edmond database at https://doi.org/10.17617/3.NP26AI.

References

Appendix A: Sets of equations for the different models

Since there are subtle differences in the equations for the different models, we explicitly state the equations for the different models. The terms in the equations below are explained in Sect. 3.

A.1. Compressible model

The complete set of linearized equations for the compressible model is:

$$\begin{array}{cc}& \frac{\partial \mathbf{u}}{\partial t}+{\mathbf{v}}_0·\mathrm{\nabla }\mathbf{u}+\mathbf{u}·\mathrm{\nabla }{\mathbf{v}}_0+2{\mathrm{\Omega }}_0{\mathbf{e}}_z\times \mathbf{u}+\mathbf{\nabla }\left(\frac{p_1}{{\rho }_0}\right)\hfill \\ \hfill & -\frac{s_1}{c_p}g{\mathbf{e}}_r-\frac{1}{{\rho }_0}\mathrm{\nabla }·\left({\rho }_0\nu (\mathrm{\nabla }(\mathbf{u})+\mathrm{\nabla }({\mathbf{u}}^T)-\frac{2}{3}(\mathrm{\nabla }·\mathbf{u})\overleftrightarrow{I})\right)=0,\hfill \end{array}$$(A.1)

$$\begin{array}{c}\hfill \frac{\partial s_1}{\partial t}+{\mathbf{v}}_0·\mathrm{\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}\mathrm{\nabla }·(\kappa {\rho }_0{T}_0\mathrm{\nabla }{s}_1)=0,\end{array}$$(A.2)

$$\begin{array}{c}\hfill \frac{\partial {\rho }_1}{\partial t}+\mathrm{\nabla }·({\rho }_0\mathbf{u})+\mathrm{\nabla }·({\rho }_1{\mathbf{v}}_0)=0,\end{array}$$(A.3)

$$\begin{array}{c}\hfill \frac{p_1}{p_0}=\gamma \frac{{\rho }_1}{{\rho }_0}+\frac{s_1}{c_v}.\end{array}$$(A.4)

For the compressible model, we solve Eqs. (A.1) – (A.3) with the constraint Eq. (A.4).

A.2. Anelastic approximation

We use the anelastic formulation, which conserves energy in both the convection and radiative zones and is recommended by Brown et al. (2012). The linearized equations for the anelastic model are:

$$\begin{array}{cc}& \frac{\partial \mathbf{u}}{\partial t}+{\mathbf{v}}_0·\mathrm{\nabla }\mathbf{u}+\mathbf{u}·\mathrm{\nabla }{\mathbf{v}}_0+2{\mathrm{\Omega }}_0{\mathbf{e}}_z\times \mathbf{u}+\mathrm{\nabla }\left(\frac{p_1}{{\rho }_0}\right)\hfill \\ \hfill & -\frac{s_1}{c_p}g{\mathbf{e}}_r-\frac{1}{{\rho }_0}\mathrm{\nabla }·\left({\rho }_0\nu (\mathrm{\nabla }(\mathbf{u})+\mathrm{\nabla }({\mathbf{u}}^T)-\frac{2}{3}(\mathrm{\nabla }·\mathbf{u})\overleftrightarrow{I})\right)=0,\hfill \end{array}$$(A.5)

$$\begin{array}{c}\hfill \frac{\partial s_1}{\partial t}+{\mathbf{v}}_0·\mathrm{\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}\mathrm{\nabla }·(\kappa {\rho }_0{T}_0\mathrm{\nabla }{s}_1)=0,\end{array}$$(A.6)

$$\begin{array}{c}\hfill \mathrm{\nabla }·({\rho }_0\mathbf{u})=0.\end{array}$$(A.7)

In this case, we solve Eqs. (A.5) – (A.6), with the anelastic constraint Eq. (A.7). In the Navier-Stokes equation, we implicitly assume that it follows the equation of state (A.4), but we do not use it as an additional constraint when solving the eigenvalue problem.

A.3. Boussinesq approximation

For the Boussinesq model, the linearized model equations are:

$$\begin{array}{cc}& \frac{\partial \mathbf{u}}{\partial t}+{\mathbf{v}}_0·\mathrm{\nabla }\mathbf{u}+\mathbf{u}·\mathrm{\nabla }{\mathbf{v}}_0+2{\mathrm{\Omega }}_0{\mathbf{e}}_z\times \mathbf{u}+\frac{1}{{\rho }_m}\mathrm{\nabla }{p}_1\hfill \\ \hfill & -\frac{s_1}{c_p}g{\mathbf{e}}_r-\mathrm{\nabla }·\left(\nu (\mathrm{\nabla }(\mathbf{u})+\mathrm{\nabla }({\mathbf{u}}^T))\right)=0,\hfill \end{array}$$(A.8)

$$\begin{array}{c}\hfill \frac{\partial s_1}{\partial t}+{\mathbf{v}}_0·\mathrm{\nabla }{s}_1+\frac{u_{\theta }}{r}\frac{\partial s_0}{\partial \theta }-\frac{1}{T_0}\mathrm{\nabla }·(\kappa T_0\mathrm{\nabla }{s}_1)=0,\end{array}$$(A.9)

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{u}=0.\end{array}$$(A.10)

Here, we again solve Eqs. (A.8) – (A.9) with the incompressibility constraint Eq. (A.10). Here, ρ_m is the mean density of the solar convection zone from the Solar model S (see Fig. 1). Unlike Triana et al. (2022), we include the Boussinesq approximation, which considers the effect of density perturbations only in the buoyancy force term. The equations for the Boussinesq approximation are derived for a compressible fluid following procedures similar to those in Spiegel & Veronis (1960).

We use the same background temperature T₀, which is almost adiabatic in the convection zone and taken from model S, for all three models. We consider no deviation from the adiabatic temperature gradient for simplicity and to match the conditions across all models. We also use the same latitudinal entropy gradient from thermal wind balance in all three models, given by

$$\begin{array}{c}\hfill \frac{\partial s_0}{\partial \theta }=r^{2}sin\theta \frac{g}{c_p}\frac{\partial ({\mathrm{\Omega }}^2)}{\partial z}.\end{array}$$(A.11)

Appendix B: Analysis of linearized vorticity equation for the inertial modes

In this appendix, we present detailed analyses of the terms in the radial and z vorticity equations (Eqs. (14) and (15)). Figure B.1 compares the significance of the left and right-hand side terms of the linearized vorticity equations for the n = 0 equatorial Rossby modes with m = 3 from the compressible and Boussinesq models under uniform rotation. The planetary β effect represented by ℜ[𝒢_planetary] drives the retrograde propagation of these modes. Since the n = 0 equatorial Rossby modes are toroidal under uniform rotation, no difference is observed between the compressible and Boussinesq results. Figure B.2 shows the results for the m = 10 HFR mode. It indicates that the retrograde propagation of the HFR modes is driven by combined effects of 𝒢_planetary, 𝒢_comp, 𝒢_other, and 𝒲_topographic. The radial motions associated with the z-vorticity induce the compressional β effect ℜ[𝒲_comp], which counteracts the other terms and promotes prograde propagation. In the Boussinesq model, the mode frequencies become more retrograde due to the absence of 𝒲_comp. Figure B.3 presents the results for the n = 1 equatorial Rossby mode with m = 3. In contrast to the n = 0 equatorial Rossby mode, the n = 1 equatorial Rossby mode exhibits both radial and z vortical motions. For the radial component of the vorticity ζ_r, the planetary β effect is the primary driver of their retrograde propagation. We also show that the topographic β effect represented by ℜ[𝒲_topographic] enhances the retrograde propagation of their z-vorticity ζ_z. In the Boussinesq model, the mode frequencies become slightly more prograde due to the absence of the retrograde-promoting 𝒲_comp. Figure B.4 shows the results for the m = 3 prograde columnar mode. The z-vortical motions dominate over the radial-vortical motions in these modes. In the compressible model, their prograde propagation is primarily driven by the compressional β effect represented by ℜ[𝒲_comp] and weakly influenced by the topographic β effect represented by ℜ[𝒲_topographic]. In the Boussinesq model, ℜ[𝒲_topographic] is the sole driver of their prograde propagation. Furthermore, the lack of diverging (converging) tendency of upflows (downflows) in the Boussinesq model amplifies the radial vortical motions at the surface, thereby enhancing the planetary β effect ℜ[𝒢_planetary], which reinforces the retrograde mode propagation. Thus, in the Boussinesq model, the prograde columnar modes exhibit significantly lower frequencies, as shown in Fig. 3.

Fig. B.1. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the n = 0 equatorial Rossby mode (m = 3) in the compressible and Boussinesq models, under uniform rotation. Refer to Eqs. (16) - (21) for the definition of the various quantities. A negative quantity implies that the associated physical effect promotes retrograde propagation, while a positive quantity implies that its physical effect promotes prograde propagation. The respective frequencies of the modes are specified on the left. We normalize the eigenfunctions from the different models to have the same integrated kinetic energy density such that the maximum absolute value of ℜ[ω|ζr|2] is 100 for the compressible model.

Fig. B.2. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the HFR mode (m = 10) in the compressible and Boussinesq models, under uniform rotation. We use the same definition and interpretation for the different quantities as in Fig. B.1. The respective frequencies of the modes are specified on the left. The eigenfunctions have the same normalization as in Fig. B.1.

Fig. B.3. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the n = 1 equatorial Rossby mode (m = 3) in the compressible and Boussinesq models, under uniform rotation. We use the same definition and interpretation for the different quantities as in Fig. B.1. The respective frequencies of the modes are specified on the left. The eigenfunctions have the same normalization as in Fig. B.1.

Fig. B.4. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the prograde columnar mode (m = 3) in the compressible and Boussinesq models, under uniform rotation. We use the same definition and interpretation for the different quantities as in Fig. B.1. The respective frequencies of the modes are specified on the left. We normalize the eigenfunctions from the different models to have the same integrated kinetic energy density such that the maximum of ℜ[ω|ζz|2] is 100 for the compressible model.

Appendix C: Quantitative comparison of anelastic and Boussinesq models with respect to compressible model

Fig. C.1. Effects of the non-toroidicity of the eigenmodes on the absolute difference between the real part of the eigenfrequencies for the compressible and Boussinesq models. The non-toroidicity Γ of the eigenmodes is calculated using Eq. (22) for the compressible model. The left column corresponds to uniform rotation, while the right corresponds to solar differential rotation. The different colours denote the different modes, similar to Fig. 5. Points of the same colour correspond to modes with different azimuthal orders m.

The following Tables C.1 – C.4 present the frequencies ω for the three models and non-toroidicity Γ of different modes under uniform rotation. Tables C.5 – C.9 report the same for the modes under solar differential rotation and its corresponding entropy gradient. The non-toroidicity, represented as Γ, is calculated using the eigenfunction for the compressible model with Eq. (22). The tables indicate that the compressible and anelastic models yield the same frequencies accurate to 0.1 nHz. In contrast, the Boussinesq model produces significantly different frequencies for almost all the modes except the n = 0 equatorial Rossby modes under uniform rotation. The computation of the eigenmodes by the compressible model requires about 1.33 times the computational time needed for the anelastic model. Thus, the anelastic model reduces computational cost without altering the solar inertial modes.

Figure C.1 shows the non-toroidicity of the modes computed using the compressible model plotted against the absolute difference in frequencies between the Boussinesq and compressible models under uniform and solar differential rotation. In general, the frequency difference increases with the non-toroidicity of the different inertial modes. The trend is most pronounced for the prograde columnar modes under both uniform and differential rotation and the high-latitude modes under differential rotation. It is also evident for the n = 1 equatorial Rossby modes under uniform rotation. However, the frequency difference decreases with increasing azimuthal order for these modes (see Fig. 3). Under differential rotation, the non-toroidal effects of critical latitudes become more significant with increasing azimuthal order. These competing effects complicate the relationship between non-toroidicity and frequency difference for the n = 1 equatorial Rossby modes under differential rotation. The variations in frequency difference and non-toroidicity remain relatively small for the HFR modes under uniform rotation. However, under differential rotation, critical latitudes introduce variations in the HFR modes. The n = 0 equatorial Rossby modes are nearly toroidal and exhibit no frequency difference under uniform rotation. Under differential rotation, they acquire a slight non-toroidal character and exhibit small frequency differences due to the radial motions near critical latitudes at higher values of m. Although non-toroidicity provides a rough estimate of the differences in eigenmodes caused by the Boussinesq approximation, it is not a precise measure, as additional physical effects contribute to the characteristics of the inertial modes.

All Tables

All Figures

	Fig. 1. Variation of the background density in the solar convection zone from the standard model S (Christensen-Dalsgaard et al. 1996), denoted by the blue curve. The black dashed line shows the mean density ρm in the solar convection zone. The Boussinesq setup uses the mean density, whereas the other stratified setups use the density from model S.
In the text

Fig. 2. Comparison of the eigenmodes of the different classes of inertial modes (see Table 1) computed using compressible, anelastic, and Boussinesq models under uniform rotation. Note that the unstable high-latitude inertial modes are absent under solid body rotation. Here, we plot the real part of uθ and the imaginary part of uϕ of the computed eigenmodes. The longitudes corresponding to the real and imaginary phases of the eigenfunctions are ϕ = ϕ0 and ϕ = ϕ0 − π/2m, respectively, where ϕ0 is a longitude, where uθ attains its maximum. The corresponding frequencies measured in the Carrington frame are stated below each eigenmode. The imaginary parts of the frequencies indicate the growth rates of the modes. All eigenfunctions are normalized such that the maximum of uθ is 1 m/s at the surface.

In the text

Fig. 3. Dispersion relations of the studied classes of inertial modes under uniform rotation for the range of azimuthal orders 0 ≤ m ≤ 15. The blue, green, orange, and red colours represent the n = 0 equatorial Rossby, the HFR, the n = 1 equatorial Rossby, and the prograde columnar modes, respectively. Open circles, cross symbols, and dotted points show the results computed from the fully compressible, anelastic, and Boussinesq models, respectively.

In the text

Fig. 4. Estimates of the importance of the compressional β effect under uniform rotation using ℜ[𝒲comp], as defined in Eq. (19), for the compressible setup. A negative value of ℜ[𝒲comp] implies that the compressional β effect promotes retrograde propagation, while a positive value implies that it promotes prograde propagation. Column (a) denotes the HFR mode with m = 10, column (b) denotes the n = 1 equatorial Rossby mode with m = 3, and column (c) denotes the prograde columnar mode with m = 3. For all modes, we normalize ℜ[𝒲comp] such that its maximum absolute value is 1.

In the text

	Fig. 5. Dispersion relations of the various classes of inertial modes listed in Table 1, computed using different models with solar differential rotation and the associated latitudinal entropy gradient. The same notation applies to colours and symbols as in Fig. 3. The only addendum is the high-latitude mode with north-south symmetric radial vorticity, represented in purple.
In the text

Fig. 6. Comparison of the eigenmodes for different classes of inertial modes computed using compressible, anelastic, and Boussinesq models under solar differential rotation and the associated latitudinal entropy gradient. The different rows denote: a) n = 0 equatorial Rossby mode (m = 3), b) high-latitude mode with north-south symmetric radial vorticity (m = 1), c) HFR mode (m = 10), d) n = 1 equatorial Rossby mode (m = 3), e) prograde columnar mode (m = 3). Here, we plot the real part of uθ and the imaginary part of uϕ of the computed eigenmodes in the meridional plane. The corresponding longitudes are chosen in the same way as in Fig. 2. The black solid curves on the meridional cross-sections denote the critical latitudes where ℜ[ω]=m(Ω − Ω0). We also show the surface velocity: uϕ for the high-latitude mode and the prograde columnar mode, and uθ for the other modes. The corresponding frequencies in the Carrington frame are stated below each eigenmode. The imaginary parts of the frequencies are the growth rates of the modes. All eigenfunctions, except for the high-latitude mode, are normalized such that the maximum of uθ is 1 m/s at the surface. For the high-latitude mode, the maximum velocity at the surface is set to 10 m/s.

In the text

Fig. 7. Estimates of the significance of the radial motions associated with the different inertial modes for the compressible setup, quantified using non-toroidicity Γ defined in Eq. (22). The values of non-toroidicity for different modes are plotted against azimuthal orders ranging from m = 1 to m = 10. Different colours represent different modes, as in Fig. 5. Solid lines denote modes under differential rotation, while dotted lines denote modes under uniform rotation. Open circles indicate modes affected by critical latitudes under differential rotation, while filled circles indicate other modes under differential rotation.

In the text

Fig. B.1. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the n = 0 equatorial Rossby mode (m = 3) in the compressible and Boussinesq models, under uniform rotation. Refer to Eqs. (16) - (21) for the definition of the various quantities. A negative quantity implies that the associated physical effect promotes retrograde propagation, while a positive quantity implies that its physical effect promotes prograde propagation. The respective frequencies of the modes are specified on the left. We normalize the eigenfunctions from the different models to have the same integrated kinetic energy density such that the maximum absolute value of ℜ[ω|ζr|2] is 100 for the compressible model.

In the text

Fig. B.2. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the HFR mode (m = 10) in the compressible and Boussinesq models, under uniform rotation. We use the same definition and interpretation for the different quantities as in Fig. B.1. The respective frequencies of the modes are specified on the left. The eigenfunctions have the same normalization as in Fig. B.1.

In the text

Fig. B.3. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the n = 1 equatorial Rossby mode (m = 3) in the compressible and Boussinesq models, under uniform rotation. We use the same definition and interpretation for the different quantities as in Fig. B.1. The respective frequencies of the modes are specified on the left. The eigenfunctions have the same normalization as in Fig. B.1.

In the text

Fig. B.4. Relative importance of planetary β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{planetary}}]$), compressional β effect ($\mathfrak{R}[{\mathcal{G}}_{\mathrm{comp}}],\mathfrak{R}[{\mathcal{W}}_{\mathrm{comp}}]$), topographical β effect ($\mathfrak{R}[{\mathcal{W}}_{\mathrm{topographic}}]$) and other terms in the vorticity equation to determine the propagation and frequency of the prograde columnar mode (m = 3) in the compressible and Boussinesq models, under uniform rotation. We use the same definition and interpretation for the different quantities as in Fig. B.1. The respective frequencies of the modes are specified on the left. We normalize the eigenfunctions from the different models to have the same integrated kinetic energy density such that the maximum of ℜ[ω|ζz|2] is 100 for the compressible model.

In the text

Fig. C.1. Effects of the non-toroidicity of the eigenmodes on the absolute difference between the real part of the eigenfrequencies for the compressible and Boussinesq models. The non-toroidicity Γ of the eigenmodes is calculated using Eq. (22) for the compressible model. The left column corresponds to uniform rotation, while the right corresponds to solar differential rotation. The different colours denote the different modes, similar to Fig. 5. Points of the same colour correspond to modes with different azimuthal orders m.

In the text