Inertial Waves in a Nonlinear Simulation of the Sun's Convection Zone and Radiative Interior

This numerical model was generated using the Rayleigh convection code (Featherstone & Hindman 2016a; Matsui et al. 2016; Featherstone et al. 2021), which solves the fluid equations in rotating spherical geometry using the pseudo-spectral algorithm from Glatzmaier (1984) and Clune et al. (1999). Each physical variable is represented as a linear combination of spectral components, with each component consisting of the product of a spherical harmonic ${Y}_{l}^{m}(\theta ,\phi )$ and a Chebyshev polynomial T_k (r). The three spherical coordinates r, θ, and ϕ, are the radius, colatitude, and longitude, respectively, and $\hat{{\boldsymbol{r}}}$, $\hat{{\boldsymbol{\theta }}}$, and $\hat{{\boldsymbol{\phi }}}$ represent the local unit vectors pointing in these three directions. In the spectral representation, l is the harmonic degree, m is the azimuthal order, and k is the radial order. The fluid variables are updated at each time step using a second-order Crank–Nicolson scheme that handles most of the linear terms in the fluid equations directly in spectral space. The nonlinear terms and the Coriolis term are calculated using an Adams–Bashforth algorithm that is performed in physical space after transformation from the spectral representation.

Because our flows are low Mach number, we linearize the thermodynamic variables about a reference state and assume a solenoidal mass flux. Specifically, we use the Lantz–Braginsky–Roberts (LBR) formulation (Lantz 1992; Braginsky & Roberts 1995) of the anelastic MHD equations (e.g., Gough 1969; Gilman & Glatzmaier 1981) in the rotating frame, given by

$$\begin{eqnarray}&&\begin{array}{c}\hat{\rho }\left[\displaystyle \frac{{\rm{\partial }}{\boldsymbol{v}}}{{\rm{\partial }}t}+{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{v}}+2{{\boldsymbol{\Omega }}}_{0}\times {\boldsymbol{v}}\right]=\displaystyle \frac{g\hat{\rho }\,{s}^{{\rm{{\prime} }}}}{{c}_{P}}\hat{{\boldsymbol{r}}}\\ -\,\hat{\rho }{\boldsymbol{\nabla }}\left(\displaystyle \frac{{P}^{{\rm{{\prime} }}}}{\hat{\rho }}\right)+\displaystyle \frac{1}{4\pi }\left({\boldsymbol{\nabla }}\times {\boldsymbol{B}}\right)\times {\boldsymbol{B}}+{\boldsymbol{\nabla }}\cdot { \mathcal D },\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{c}\hat{\rho }\hat{T}\left[\frac{{\rm{\partial }}{s}^{{\rm{{\prime} }}}}{{\rm{\partial }}t}+{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{s}^{{\rm{{\prime} }}}+{v}_{r}\frac{d\hat{s}}{{dr}}\right]={\boldsymbol{\nabla }}\cdot \left[\hat{\rho }\hat{T}\kappa {\boldsymbol{\nabla }}{s}^{{\rm{{\prime} }}}\right]\\ +\,Q+{\rm{\Phi }}+\frac{\eta }{4\pi }{\left|{\boldsymbol{\nabla }}\times {\boldsymbol{B}}\right|}^{2},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times \left({\boldsymbol{v}}\times {\boldsymbol{B}}-\eta {\boldsymbol{\nabla }}\times {\boldsymbol{B}}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}_{{\rm{ij}}}=2\hat{\rho }\nu \left[{e}_{{\rm{ij}}}-\frac{1}{3}\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}\right){\delta }_{{\rm{ij}}}\right],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\rm{\Phi }}=2\hat{\rho }\nu \left[{e}_{{ij}}{e}_{{ij}}-\displaystyle \frac{1}{3}{\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}\right)}^{2}\right],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left(\hat{\rho }{\boldsymbol{v}}\right)=0,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 7 }$$

where the primary variables are the velocity v , magnetic field B , and the fluctuations of the gas pressure $P^{\prime}$ and specific entropy density $s^{\prime}$ about the corresponding reference-state profiles. The angular velocity of the rotating reference frame, Ω₀, is aligned with the axis of the spherical coordinate system. The gravitational acceleration is given by g(r); the viscous, thermal, and magnetic diffusivities are given by ν(r), κ(r), and η(r), respectively; c_P is the specific heat capacity at constant pressure; e_ij is the rate-of-strain tensor; and δ_ij is the Kronecker delta. The radial profiles $\hat{\rho }(r)$ and $\hat{T}(r)$ are the time-independent, spherically symmetric reference-state density and temperature, which satisfy the ideal gas law, and ${\rm{d}}\hat{s}/{\rm{d}}r$ is the reference-state entropy gradient. The equation of state for an ideal gas is linearized and expressed in terms of the fluctuations about this reference state. The momentum equation includes the Coriolis, buoyancy, pressure gradient, Lorentz, and viscous forces. The thermal energy equation includes conductive, ohmic, and viscous heating, as well as an internal-heating term Q(r) that represents radiative heating. The internal heating is normalized to deposit a solar luminosity of heat, which is distributed preferentially in roughly the bottom third of the convection zone (for details, see Featherstone & Hindman 2016a).

This MHD model evinces a self-consistently established tachocline that is in a statistical steady state. We previously explored the properties of this tachocline and the balances that enable its formation in Matilsky et al. (2022; see also Matilsky et al. 2023). The simulation spans a spherical shell covering the upper radiative interior and lower convection zone, extending between the radii ${r}_{\min }=0.491\,{R}_{\odot }$ and ${r}_{\max }=0.947\,{R}_{\odot }$, where R_⊙ = 6.957 × 10¹⁰ cm is the solar radius. The transition from convective stability to instability is primarily controlled by the entropy gradient of the reference state ${\rm{d}}\hat{s}/{\rm{d}}r$ and occurs at r_bcz = 0.729 R_⊙. Convective downflows overshoot into a thin layer within the stable region, the base of which is r_ov = 0.710 R_⊙.

The spherical shell spans roughly the upper two density scale heights of the radiative interior and the bottom three density scale heights of the Sun's convection zone. The horizontal grid resolution is N_θ = 384 and N_ϕ = 768, corresponding to a maximum spherical harmonic degree after dealiasing (using the two-thirds rule) of l = 255. In the radial direction, there are three stacked Chebyshev domains with 64 points each, with interior boundaries located at 0.669 R_⊙ and 0.719 R_⊙. The stacked grids provide substantially higher radial resolution in the tachocline and in the overshooting layer (Matilsky et al. 2022).

The background entropy gradient is a positive constant in the radiative interior and zero in the convection zone. The following smooth profile connects the two domains:

$$\begin{eqnarray}\displaystyle \frac{d\hat{s}}{{dr}}=\sigma \left\{\begin{array}{ll}1 & r\leqslant {r}_{0}-\delta ,\\ 1-{\left[1-{\left(\displaystyle \frac{r-{r}_{0}}{\delta }\right)}^{2}\right]}^{2} & {r}_{0}-\delta \lt r\lt {r}_{0},\\ 0 & r\geqslant {r}_{0},\end{array}\right.\end{eqnarray} \tag{ 8 }$$

where σ ≡ 10⁻² erg g⁻¹ K⁻¹ cm⁻¹, δ ≡ 0.05 R_⊙, and r₀ = 0.719 R_⊙. The gravitational acceleration is given by

$$\begin{eqnarray}&&g(r)=\displaystyle \frac{{{GM}}_{\odot }}{{r}^{2}},\end{eqnarray} \tag{ 9 }$$

where G is the universal gravitational constant and M_⊙ = 1.989 × 10³³ g is the solar mass.

The rotation rate is given by Ω₀/2π = 1370 nHz, and a solar luminosity L_⊙ = 3.85 × 10³³ erg s⁻¹ is driven through the convection zone via the temporally steady internal-heating profile Q(r) (for details, see Matilsky et al. 2022). The viscous, thermal, and magnetic diffusivities at the top of the domain are given by $\nu ({r}_{\max })=\kappa ({r}_{\max })=5\times {10}^{12}$ cm² s⁻¹ and $\eta ({r}_{\max })\,=1.25\times {10}^{12}$ cm² s⁻¹, respectively. All diffusivities increase with radius, varying like $\sim {\hat{\rho }}^{-1/2}$.

At both boundaries, we use stress-free and impenetrable conditions on the velocity, i.e.,

$$\begin{eqnarray}&&{v}_{r}=\displaystyle \frac{\partial \,}{\partial r}({v}_{\theta }/r)=\displaystyle \frac{\partial \,}{\partial r}({v}_{\phi }/r)=0.\end{eqnarray} \tag{ 10 }$$

We impose a fixed conductive flux at the inner and outer spherical boundaries, using

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial s^{\prime} }{\partial r}\right|}_{{r}_{\min }}=0,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial s^{\prime} }{\partial r}\right|}_{{r}_{\max }}={\left.-\displaystyle \frac{{L}_{\odot }}{4\pi {r}^{2}\kappa (r)\hat{\rho }(r)\hat{T}(r)}\right|}_{{r}_{\max }}.\end{eqnarray} \tag{ 12 }$$

No heat is allowed to enter at the bottom of the domain, and at the top of the domain, heat exits via thermal conduction. Potential field matching conditions are used on the magnetic fields at both boundaries.

Although we ran the simulation dimensionally, its parameter regime can be uniquely described by the nondimensional parameters given in Table 1, where the tildes indicate volume-averaged quantities and $H={r}_{\max }-{r}_{\mathrm{bcz}}$ is the depth of the convection zone. We adopt the flux Rayleigh number Ra_F, where F is the energy flux imposed by radiative heating, i.e., the flux not carried by the radiation field (see Featherstone & Hindman 2016b). Ek is the Ekman number, which controls the importance of viscosity relative to rotation. The convective Rossby number Ro_c characterizes the relative strength of the buoyancy and Coriolis forces, and it often determines the degree of rotational influence on the convection (Aurnou et al. 2020; Camisassa & Featherstone 2022). The Prandtl number Pr and the magnetic Prandtl number Pr_m specify the various ratios of the viscous, thermal, and magnetic diffusion coefficients. The buoyancy parameter Bu assesses the stiffness of the radiative interior.

Table 1. Nondimensional Fluid Parameters

Parameter	Definition	Value
RaF	$\tfrac{\tilde{g}\tilde{F}{H}^{4}}{{c}_{P}\tilde{\rho }\tilde{T}\tilde{\nu }{\tilde{\kappa }}^{2}}$	5.68 × 105
Ek	$\tfrac{\tilde{\nu }}{2{{\rm{\Omega }}}_{0}{H}^{2}}$	5.35 × 10−4
Pr	$\tfrac{\tilde{\nu }}{\tilde{\kappa }}$	1
Prm	$\tfrac{\tilde{\nu }}{\tilde{\eta }}$	4
Roc	${\left(\tfrac{{\mathrm{Ra}}_{{\rm{F}}}\,{\mathrm{Ek}}^{2}}{\Pr }\right)}^{\tfrac{1}{2}}$	0.403
Bu	$\tfrac{{\tilde{N}}^{2}}{4{{\rm{\Omega }}}_{0}^{2}}$	6.35 × 103

Notes. The flux Rayleigh number Ra_F , Ekman number Ek, Prandtl number Pr, magnetic Prandtl number Pr_m , convective Rossby number Ro_c , and buoyancy parameter Bu for the numerical simulation have been defined using volumetric averages of the reference-state profiles. Averages for the first five parameters are taken over the convection zone only, while the buoyancy frequency is averaged over the radiative interior only. These quantities are indicated by use of a tilde over the variable.

Download table as:  ASCIITypeset image

This model has a solar-like differential rotation profile, with the equator rotating faster than the polar regions, illustrated in Figure 1. Dynamo action creates a large-scale, cycling, non-axisymmetric magnetic field whose torque enforces solid-body rotation in the radiative interior and maintains a tachocline (see Matilsky et al. 2022). The rotation rate is given by ${\rm{\Omega }}(r)={{\rm{\Omega }}}_{0}+\langle {v}_{\phi }\rangle /r\sin \theta$, where the angular brackets refer to a combined temporal and zonal average, and Ω₀ is the frame rate. Figure 1 shows the nondimensionalized differential rotation profile (Ω − Ω₀)/Ω₀. Here and in Figure 2, the pink dashed–dotted lines indicate the base of the tachocline at r_tach = 0.641 R_⊙, and the black dashed lines indicate the base of the convection zone at 0.729 R_⊙. The model evinces 0.04 Ω₀ difference in the rotation rate from pole to equator in the convection zone, whereas in the radiative interior the pole-to-equator contrast is roughly 0.0015 Ω₀.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 2(a) and (b) display the kinetic and magnetic energy densities as functions of radius across the entire domain. The contributions to the kinetic energy from the horizontal and radial velocity components are similar in magnitude throughout the convection zone. Across the tachocline, the contribution from the radial velocity component drops by 8 orders of magnitude, while the horizontal velocities only decrease by 4 orders of magnitude. The fact that horizontal velocities tend to be significantly higher than vertical velocities in simulations has been noted multiple times in past work (e.g., Alvan et al. 2015; Lawson et al. 2015). In Section 4, we determine that this enhanced horizontal power is due to a rich spectrum of equatorial Rossby waves. Similarly, the magnetic field is primarily horizontal in the radiative interior compared to the convection zone, with the radial field component being smaller by roughly an order of magnitude (or 2 orders of magnitude in the square of the field components).

Figure 2(c) shows the radial behavior of the dynamic Elsässer number as defined by Soderlund et al. (2012):

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{d}(r)=\displaystyle \frac{1}{4\pi }\displaystyle \frac{{B}^{2}}{2\hat{\rho }{{\rm{\Omega }}}_{0}{UL}},\end{eqnarray} \tag{ 13 }$$

where B(r) is the spherically averaged rms field strength, L is the characteristic length scale for the magnetic field, and U(r) is the spherically averaged rms flow speed. We choose L = π R_⊙/2, or a length scale corresponding to a large-scale field component with a harmonic degree of l = 2. Matilsky et al. (2023) have demonstrated that the magnetic field in this simulation is concentrated in the low-order spherical harmonic components. The Elsässer number compares the strength of the Coriolis force to that of the Lorentz force. An Elsässer number greater than one indicates that the Lorentz force dominates the dynamics, whereas an Elsässer number of less than one indicates that the Coriolis force dominates. The Elsässer number in this simulation is always quite small (Λ_d < 10⁻²), hence the wave dynamics are rotationally rather than magnetically constrained. Thus, even though the magnetic field is crucial in the dynamics of the mean flows, the inertial waves can essentially be treated as nonmagnetic and hydrodynamic.

If the modes live in a region rotating at a rate Ω different from our frame rate Ω₀ (i.e., Ω = Ω₀ + ΔΩ), then a wave of azimuthal order m and intrinsic frequency $\omega ^{\prime}$ will be observed in the reference frame of the simulation (rotating at rate Ω₀) as having the Doppler-shifted frequency

$$\begin{eqnarray}&&\omega =\omega ^{\prime} +m{\rm{\Delta }}{\rm{\Omega }},\end{eqnarray} \tag{ 14 }$$

(e.g., Bretherton & Garrett 1968).

It is important to note that both the observed and intrinsic frequency of a normal mode do not change with spatial position; the same ΔΩ correction can be applied across the entire domain. The Doppler shift that relates these two frequencies is a spatial average of the rotation rate difference over the region where the mode's eigenfunction has a significant amplitude. Therefore, the Doppler shift depends on the radial and latitudinal eigenfunctions. For waves living in the radiative interior, which is rotating mostly like a solid body, we can easily apply a ΔΩ that is simply the difference between our frame rate and the rotation rate of the radiative interior. Because the radiative interior is rotating slightly slower than our frame rate, the correction is negative (ΔΩ < 0). For waves that reside in the convection zone, the situation is more complicated because the rotation profile of the convection zone varies significantly in both radius and latitude. Different modes that occupy different parts of the convection zone will sense a different mean rotation rate. Thus, without knowing the detailed shape of the eigenfunctions, we can only make informed guesses for the exact Doppler shift that should be applied.

The left-hand column of Figure 3 makes clear that both sectoral and tesseral equatorial Rossby waves are present in the radiative interior of our simulation. Figure 4(a) displays the power spectrum in radial vorticity, averaged over the radiative interior, summed over all latitudinal orders (or equivalently, over all harmonic degrees). The black dots denote a theoretical estimate of the frequencies,

$$\begin{eqnarray}&&\omega =m{\rm{\Delta }}{\rm{\Omega }}-\displaystyle \frac{2m{{\rm{\Omega }}}_{0}}{l(l+1)}.\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

This dispersion relation is derived for purely hydrodynamic, equatorial Rossby waves in two horizontal dimensions with solid-body rotation at a frame rate of Ω₀ (Haurwitz 1940). For such a spherically symmetric system, the eigenfunctions for the Rossby waves are pure spherical harmonics. Because the radiative interior of our model is rotating slightly slower than the frame rate (with a difference of ΔΩ ≈ −0.0015 Ω₀; see Figure 1), we have applied a Doppler shift to the traditional dispersion relation.

In Figure 4(a), the sectoral modes are the ridge labeled λ = 0. The tesseral modes have frequencies smaller in magnitude (i.e., less negative), with each ridge moving closer to zero frequency as the latitudinal order λ increases. The first set of tesseral modes with λ = 1 are labeled. As expected, all of the equatorial Rossby waves have a negative frequency and thus propagate in the retrograde direction. Ridges with λ = 2, 3, and 4 are also visible. Higher-order ridges are less clear because the spacing between adjacent ridges falls below their line widths, leading to blending of the individual peaks. This blurring is more obvious in Figure 4(b), which displays a cut through the radial-vorticity power spectra at m = 2 for 0.595 R_⊙, summed over all l. One can clearly identify peaks matching modes for m = 2 and l = 2, 3, 4, 5, and 6. In addition to these distinct peaks, at low frequencies, there is a collection of modes with l > 6 that have blended together. The green dotted line in Figure 4(b) denotes a contribution from a mode with a different m-value, specifically the m = 3 sectoral mode. This "leakage" results from an extremely weak misalignment of the mean angular momentum vector with the axis of the spherical coordinate system caused by numerical noise in the simulation.

Modes with λ as high as 10 are visible in the power profiles for individual spherical harmonic components, as shown in Figure 5 for m = 2 and m = 6. Each panel displays, from left to right, peaks with λ = 0 to λ = 10. The separation between the peaks decreases with increasing λ, but they are all clearly separated from the zero-frequency feature. With this in mind, the radiative interior of our model features far more equatorial Rossby-wave modes than are obvious when they are latitudinally summed (note the blurring in Figure 4(a) and the pileup of modes near zero frequency in Figure 4(b)).

Zoom In Zoom Out Reset image size

Almost all of the horizontal power in the radiative interior is in the inertial band of frequencies, and the power profiles in Figures 4 and 5 rise up to 8 orders of magnitude above the power background, indicating that these equatorial Rossby waves are the most prominent contributor to the horizontal velocity components v_θ and v_ϕ . The rms velocity summed across all modes has a typical speed of 50 cm s⁻¹, which explains the enhanced velocities in the radiative interior noted by Lawson et al. (2015) and Matilsky et al. (2022) and suggests that the radiative interior is not as quiescent as is often believed.

For completeness, we note that all of the temporal spectra for individual spherical harmonic components (see Figure 5) display a common power peak at a slightly negative frequency that is very near zero. In the spectra for m = 6 (panel (b)), this feature appears around −7 nHz (or ω/Ω₀ ≈ 5 × 10⁻³). This power feature is is likely due to the dynamo cycle, which has a similar frequency (see Matilsky et al. 2022, 2023).

An eigenmode for a Rossby wave living on a 2D spherical surface undergoing solid-body rotation has a single spherical harmonic for its eigenfunction and an eigenfrequency given by Equation (15). While the frequencies of the Rossby waves in our model are described well by this equation, small discursions from spherical symmetry in the radiative interior, such as differential rotation (see the inset in Figure 1) and magnetism, cause small perturbations to both the eigenfrequencies and the eigenfunctions. Thus, we anticipate that a given horizontal eigenfunction will be dominated by a single spherical harmonic with a given harmonic degree l₀ (${Y}_{{l}_{0}}^{m}$), but will also contain weak contributions from spherical harmonics Y_l ^m with the same azimuthal order m and nearby values of the harmonic degree l ≈ l₀. Figure 5 clearly evinces such behavior. A spectrum for a specific harmonic degree has not only a single large peak but also a sequence of small subsidiary peaks at the frequencies associated with other harmonic degrees. The upper panel of Figure 6 shows the power spectra in radial vorticity for modes with m = 2, plotted versus frequency and harmonic degree l. Though most power sits in the dominant spherical harmonic ${Y}_{{l}_{0}}^{m}$, there is some power that spreads across other l values for each mode. The lower panel displays power integrated around a narrow frequency band centered around the appropriate peak for three modes with m = 2 and l₀ = 3, 4, and 5. It is clear that the power peaks at the harmonic degree l₀ of the dominant spherical harmonic component, but the eigenfunctions have small perturbations (probably due to the differential rotation) that spread significant power over harmonic degrees within ±2 of the dominant degree.

Zoom In Zoom Out Reset image size

Within the convection zone, equatorial Rossby waves are still present, though they are far less prominent than in the radiative interior due to the presence of other convective and wave motions. From the right-hand column of Figure 3, we note that the sectoral modes for 2 ≤ m ≤ 8 remain, along with several tesseral modes for λ ≤ 4. This figure is averaged over the entire convection zone, so these modes are present somewhere in that region and not necessarily significant throughout.

In Figure 7, the black dots denote the theoretical values given by Equation (15). Because we do not know where these modes live and the shift appears to be small, we have opted not to apply a Doppler shift. At the base of the convection zone (left-hand panels), sectoral and tesseral modes are present and distinct for roughly m < 11. At the radius closest to the surface (right-hand panels), the only modes that are clearly visible are the sectoral modes for roughly m ≤ 5 (Figure 7(c), λ = 1 tesseral modes for m < 7 (Figure 7(f), and the λ = 2 tesseral modes for m < 3 (Figure 7(c). Further discussion of mode properties is deferred to the Discussion (Section 6.1).

Like equatorial Rossby waves, thermal Rossby waves result from the conservation of potential vorticity (e.g., Roberts 1968; Busse 1970). Rather than featuring motions confined to a spherical shell (i.e., v_θ and v_ϕ ), thermal Rossby waves manifest as prograde-propagating equatorial convective columns aligned with the rotation axis, hence possessing velocity components in all three directions, but with the dynamically active components being in longitude, v_ϕ , and in the cylindrical radius, ${v}_{r}\sin \theta +{v}_{\theta }\cos \theta$ (Hindman & Jain 2022; Jain & Hindman 2023). These waves can exist wherever the Coriolis force rivals or dominates buoyancy. Hence, they appear in the convection zone of our model and are excluded from the radiative interior.

The gravest frequency thermal Rossby wave, here called a Busse mode, is the one that lacks latitudinal nodes in v_ϕ (Busse 1970). Thus, the wave manifests as a parade of convective columns wrapped around the equator, with each column consisting of a single, axially aligned roll with unidirectional spin. Because of this latitudinal dependence, in radial vorticity each convective column appears antisymmetric across the equator and contributes to the spherical harmonic components with odd values of λ. In the lower panels of Figures 3 and 7, this prograde-propagating thermal Rossby wave appears as the positive frequency prong (labeled "0" in Figure 7) that rises in frequency and asymptotes toward the frame rotation rate Ω₀ as the azimuthal order increases. Hindman & Jain (2022) calculated the eigenfrequencies for a neutrally stable layer in a local, equatorial, f-plane model. We have adapted their calculation for the depth of the convection zone appropriate for the simulation; the resulting eigenfrequencies are overplotted in Figure 7 as filled black squares. Because the Busse mode lives in the convection zone near the equator, which is rotating faster than the frame rate, we apply a Doppler shift to these eigenfrequencies. We do not know the specifics of the radial eigenfunction, so we cannot easily pick a theoretically motivated ΔΩ. Instead, we take the volume-averaged equatorial rotation rate (between ±15° latitude) across the entire convection zone, resulting in ΔΩ = 0.015 Ω₀.

The narrow line widths of Busse's thermal Rossby waves suggest that the modes are stable and largely linear, with only weak nonlinear coupling to other convective modes. Figure 8(a) displays select power profiles with λ = 1 for m = 2, 4, and 6 at a radius of 0.889 R_⊙, binned down in frequency by a factor of 32. We fit these binned profiles to a Lorentzian with a linear background term, and these fits are listed in Table 3 and overplotted in Figure 8(a) with black lines. Widths ranging between 0.03 Ω₀ and 0.08 Ω₀ indicate lifetimes of 10–30 rotation periods. The quality factor, Q—i.e., the ratio of the line width to the mode frequency—is about ten across all eight modes. Further details regarding line fitting are found in the Appendix.

Zoom In Zoom Out Reset image size

In addition to Busse modes, i.e., the fundamental thermal Rossby wave consisting of a single, axially aligned roll in latitude (with a radial vorticity that is antisymmetric across the equator), there is a clear signature in our spectra for the first latitudinal overtone. Such thermal Rossby waves consist of two counter-rotating rolls stacked in latitude, with one node at the equator in v_ϕ . Because thermal Rossby waves of this type were first explored by Roberts (1968), we refer to these modes as Roberts modes, and they are labeled "1" in Figure 7. Such modes have a radial vorticity that is symmetric across the equator and thus appear in the even values of λ. We have overplotted the eigenfrequency calculations of Bekki et al. (2022b) with open squares, with an applied Doppler shift of 0.015 Ω₀. Roberts (1968) first examined these modes in a Boussinesq system, but Jain & Hindman (2023) have recently considered analytical solutions in a stratified atmosphere. Such modes have also been detected by Bekki et al. (2022b) in their linear-mode calculations, where they have noted that this family of waves smoothly transitions from the branch of prograde thermal Rossby waves to a branch of retrograde inertial waves with one radial node in radial vorticity as the azimuthal order m passes through zero from positive to negative values. We find the same (Section 5.3), and following Bekki et al. (2022b), we refer to this behavior where a prograde branch and a retrograde branch are connected as mixed modes.

Figure 8(b) shows select line profiles for the Roberts modes. The modes appear in both λ = 0 and λ = 2 and appear clearest in λ = 2 at a radius near the surface of 0.935R_⊙. Because the line profiles are asymmetric about the peak frequency and the background power has a strong nonlinear frequency dependence, we have elected not to fit them for this work.

Beyond these two easily identifiable types of thermal Rossby waves, we do find evidence of further latitudinal overtones, both symmetric and antisymmetric across the equator. In Figure 7(f), there is a prograde, higher-frequency ridge of power (labeled "2") that slopes downward and asymptotes toward the rotation rate. This is likely a thermal Rossby wave with two latitudinal nodes in v_ϕ , or three stacked columns, that is antisymmetric in radial vorticity (i.e., three latitudinal nodes in radial vorticity). In Figure 8(b), an additional symmetric mode appears as a high-frequency bump in power in m = 6 (near a frequency of ω = 1.3 Ω₀), labeled "3" here and in Figure 7. The eigenfunctions of these columnar waves are not spherical harmonics, so we expect latitudinal overtones to be spread across multiple spherical harmonic components. Because this high-frequency bump appears in an equatorially symmetric spherical harmonic, this feature could indicate the presence of the mode with four columns stacked in latitude and three latitudinal nodes in v_ϕ . (In radial vorticity, there would be four latitudinal nodes.) In Section 6.2, we speculate that these high-order modes correspond to the prograde branches of mixed HFR waves.

Throughout the tachocline and convection zone, there exists in the spectra for equatorially symmetric vorticity a strong ridge of power at low m with frequencies that lie between the ridges for the λ = 0 and λ = 1 equatorial Rossby modes (top rows of Figures 3 and 7, labeled with Roman numeral I). Bekki et al. (2022b) referred to these modes as sectoral equatorial Rossby modes with one radial node (n = 1) and describe them as having a dominant λ = 0 component in v_θ and a dominant λ = 1 component in v_ϕ , which would result in the equatorially symmetric signature that we see in the radial vorticity. The Doppler-shifted values (ΔΩ = 0.015 Ω₀) from Table 2 of Bekki et al. (2022b) are overplotted as empty black circles and match the ridges rather well.

Figure 8(c) displays line profiles for the m = 2, 4, and 6 modes for λ = 0 at a radius of 0.935 R_⊙. The equatorial Rossby wave (l, m) = (2, 2) is visible as the labeled tiny peak on the far left, emphasizing that the mixed modes are significantly wider and larger in amplitude. We again perform a Lorentzian fit with a linear background. These fits are overplotted as solid black lines in Figure 8(c), and the values are listed in Table 4. The widths of these peaks are typically 0.05 Ω₀, which corresponds to mode lifetimes of 20 rotations and a Q factor of between 2 and 10. While these fits provide a good fit to the core of the line profile, and hence a reasonable estimate for the line width, the fits to the wings and the power background are clearly inadequate. In particular, the wings are slightly asymmetric. Because we do not currently have a good physical argument for this asymmetry, we have elected to fit a symmetric Lorentzian with a linear background rather than fit an asymmetric profile.

As discussed in Bekki et al. (2022b), these modes (in this work, corresponding to Ridge I) are mixed modes. Specifically, these retrograde waves are smoothly joined through m = 0 to the Roberts modes (Ridge 1 in the upper panels of Figure 7), which have one node in latitude at the equator in the variable v_ϕ (Section 5.2). In our spectra, where the prograde waves are presented with positive frequencies and the retrograde waves with negative ones, the mixed nature appears as a common absolute frequency ∣ω∣ for the m = 0 mode of each branch. The (l, m) = (2, 0) mode is plotted in Figure 8(c), and it occurs at the same absolute frequency as the corresponding m = 0 mode in Figure 8(b).

Following the observational discovery of HFR modes on the Sun by Hanson et al. (2022), Triana et al. (2022) provided numerical evidence to support the identification of HFR waves as a class separate from equatorial Rossby waves. Using a linear eigensolver for the Boussinesq equations, Triana et al. (2022) found a class of equatorially antisymmetric modes with frequencies similar to those in Hanson et al. (2022) and featuring larger poloidal (i.e., vertical) kinetic energy than Rossby waves. Bhattacharya & Hanasoge (2023) used an anelastic solver and generated qualitatively similar high-frequency modes, along with an equatorially symmetric branch that has not yet been observed. Bekki (2023) corroborated the previous eigensolver results for the primary set of antisymmetric modes.

We find evidence for both the symmetric and antisymmetric branches of the HFR modes. In Figure 7, there are three antisymmetric bands (labeled II, IV, and VI) and two symmetric bands (III and V) of retrograde-propagating inertial waves occurring at frequencies of higher magnitude than the equatorial Rossby waves. The HFR waves observed by Hanson et al. (2022), marked as plus signs on the figure with a Doppler shift of ΔΩ = 0.015 Ω₀, are antisymmetric in radial vorticity and found near the solar surface. The underlying antisymmetric band of power in our spectra is most likely these HFR modes. However, because the nature of the HFR modes is still unclear from both an observational and theoretical standpoint, this identification is based purely on the equatorial symmetry and the similarity in frequencies.

In addition to this antisymmetric ridge, we see a strong ridge of power (III) in the spectra for modes that are symmetric across the equator, for λ = 0 and λ = 2, that is similar to the symmetric λ = 0 ridge seen by Bhattacharya & Hanasoge (2023). Beyond these two prominent ridges, we note additional bands of high-frequency power (IV, V, and VI) that to our knowledge have not been reported previously.

Figure 9 displays radial vorticity power spectra at 0.935R_⊙, effectively a zoomed-in version of Figures 7(c) and 7(f). The five HFR ridges are again marked by Roman numerals II–VI, and I denotes the retrograde mixed mode. Values from Bekki et al. (2022b) are marked with open circles, and equatorial Rossby waves are marked with closed circles. Panel (a) shows power with equatorial symmetry, with values from Bhattacharya & Hanasoge (2023) overplotted as stars, while Panel (c) shows the antisymmetric power, with observations from Hanson et al. (2022) overplotted as plus signs, along with numerical values from Triana et al. (2022; triangles) and Bhattacharya & Hanasoge (2023; multiplication signs), again all Doppler shifted by 0.015 Ω₀. We note that Ridges II and III have the same equatorial symmetry and frequencies similar to those of previous observations and numerical calculations, with particular concordance with the results of Bhattacharya & Hanasoge (2023).

Zoom In Zoom Out Reset image size

These features are also visible in Figures 9(b) and (d), which display the corresponding power profiles in radial vorticity for m = 9 and λ = 0 through 4. The black dashed line marks the Doppler-shifted m = 9 observed frequency from Hanson et al. (2022). These additional enhancements in power are slight and extremely broad; hence, if modal, they have an extremely low Q factor. We note that ridges II and IV are visible in λ = 1 but are barely present in λ = 3, while ridge VI is more visible in λ = 3 than λ = 1. Ridges III and V are visible in both λ = 0 and 2. The eigenfunctions of these modes are unlikely to be spherical harmonics, so we expect power to be spread across multiple spherical harmonic components. However, the manner in which this occurs is still unknown. These features also appear to be heavily dependent on the general background profile of the spectrum. Because we do not yet have a good enough physical understanding of the background to model it, and due to the low and broad nature of these features, we have elected not to fit these line profiles.

Figure 10 shows the m = 9 power spectral line profiles for the individual velocity components at two different radii near the surface. Because the equatorial symmetry of v_θ is opposite that of v_r and v_ϕ , the modes are denoted by the symmetry consistent with the radial vorticity (which also shares the same symmetry as v_θ ); the "symmetric modes" are the sum of ∣v_θ ∣² over λ = 0, 2, while ∣v_r ∣² and ∣v_ϕ ∣² are summed over λ = 1, 3, and the antisymmetric modes are the converse. A key feature of all of these modes is a significant vertical velocity component, both deep within the convection zone and at the surface. This property has been recognized in previous studies (Triana et al. 2022; Bekki 2023) and is quite evident from the spectra in Figure 10. In Figure 10, we note that, while v_θ and v_ϕ are a couple of orders of magnitude greater than v_r at the surface, the difference reduces to less than 1 order of magnitude at the deeper radial slice, which qualitatively matches previous findings (Triana et al. 2022; Bekki 2023).

Zoom In Zoom Out Reset image size

It is highly probable that these five ridges (II–VI) and the mixed retrograde wave (I) are all in the same family, with higher-frequency ridges corresponding to modes of higher latitudinal order. Further analysis is deferred to Section 6.2.

The final feature that we identify in the spectra is the convection, which exhibits a Doppler shift due to the radial variation of the mean rotation rate. In Figure 3, this feature has been averaged over several different radii and manifests as the smear of power around and above zero frequency. Figure 11 shows the spectra in radial vorticity, for λ = 0, out to m = 80 at three different radii in the convection zone. The overplotted dashed black line shows the expected frequency resulting from the Doppler shift due to the difference between the local mean rotation rate and the frame rate Ω₀. At the base of the convection zone, the equatorial rotation rate is almost the same as the frame rate, so the convective power remains near zero frequency for all m. Higher in the convection zone, the equatorial rotation rate rises above the frame rate by as much as 0.036 Ω₀, resulting in a significant Doppler shift and the sloped distribution of power present in the center and right panels of Figure 11. The keen-eyed reader will note the sectoral equatorial Rossby waves and the retrograde mixed mode in the lower left of each panel, along with a squished thermal Rossby-wave branch (prograde frequencies) and an equatorially symmetric HFR branch (retrograde frequencies) in panel (c).

Zoom In Zoom Out Reset image size

6.1.1. Are There Two Wave Cavities?

In their study of the radial behavior of r-modes, Wolff & Blizard (1986) numerically calculated the first-order frequency correction to Equation (15) for a solar interior model, with a stably stratified radiative interior beneath an unstable convection zone. They found that there are two separate wave cavities, one in the radiative interior and another in the convection zone, with the frequency correction being of opposite sign in the two regions. Specifically, if we add the frequency correction that arises from stratification, δ ω_strat, to Equation (15):

$$\begin{eqnarray}&&\omega =m{\rm{\Delta }}{\rm{\Omega }}+\delta {\omega }_{\mathrm{strat}}-\displaystyle \frac{2m{{\rm{\Omega }}}_{0}}{l(l+1)}.\end{eqnarray} \tag{ 16 }$$

Wolff & Blizard (1986) found that the the frequency correction is positive for the convection zone cavity, δ ω_strat > 0, and negative for the cavity in the stably stratified radiative interior, δ ω_strat < 0. These corrections are constants that depend on the particular eigenmode and cavity within which the waves reside. Further, Wolff & Blizard (1986) found that such cavities can support radial overtones with differing numbers of nodes in radius. In the following subsections, we explore both of these possibilities: if the Rossby modes that we have observed in our numerical simulation live in two distinct cavities, and whether radial overtones might be present.

A resonant normal mode has a single frequency that is independent of spatial position. Hence, we can attempt to distinguish between modes in the radiative interior and modes in the convection zone by their frequencies. In addition to the frequency shift mentioned previously that arises from the nature of the stratification (negative in the radiative zone and positive in the convection zone), we expect a Doppler shift based on the difference in the rotation rate between the two zones. Because the radiative interior is rotating slightly slower than the frame rate of our model, we expect a negative ΔΩ correction in Equation (15). Conversely, the convection zone is rotating faster than the frame rate, necessitating a positive ΔΩ correction. The expected Doppler shifts thus have the same sign as the frequency shift identified by Wolff & Blizard (1986). Hence, we can identify modes living in a radiative interior cavity as those with a total negative frequency shift from mΔΩ + δ ω_strat and those in the convective zone cavity as those with a positive frequency shift.

Figure 12 displays power profiles for the (l, m) = (4, 4) sectoral mode at full frequency resolution at the base of the radiative interior and in the lower convection zone. The purple dashed line marks the expected frequency for modes in a region rotating at the frame rate (where ΔΩ = 0) and with no stratification correction, δ ω_strat = 0. We can clearly see that the power profile within the radiative interior has a mean frequency that is greater in magnitude than the purple dashed line (more negative), while the convection-zone profile has its mean at a lower-magnitude frequency (less negative).

Zoom In Zoom Out Reset image size

With this expected frequency difference in mind, we average the power in each mode over the radiative interior and over the convection zone and calculate the first frequency moment 〈ω〉 of the power distribution separately in each region (the angular brackets here refer to a power-spectrum-weighted average). Figure 13 shows the difference, 〈Δω〉 = 〈ω〉_cz − 〈ω〉_ri , in the mean frequencies for the convection zone and radiative interior, scaled by the 2D dispersion relation value ω_2D = 2mΩ₀/l(l + 1). Because 〈Δω〉 is always positive, the convection zone modes always have a lower-magnitude frequency (less negative) than those in the radiative interior, as predicted by Wolff & Blizard (1986). There is not a clear relationship between frequency difference and latitudinal node number λ, but for the sectoral modes, 〈Δω〉 does increase with m.

Zoom In Zoom Out Reset image size

The distinct shift in frequencies between equatorial Rossby waves in the radiative interior and those in the convection zone suggests that the equatorial Rossby waves are behaving as Wolff & Blizard (1986) predicted. The waves are split into two different families. One family has frequencies that are more negative than the classic dispersion relation would indicate, and these live in the radiative interior, i.e., they propagate in the radiative interior and are evanescent in the convection zone. Conversely, the other family has frequencies that are less negative than the classic 2D dispersion relation and live in the convection zone (i.e., they propagate in the convection zone and are evanescent in the radiative interior). The frequency shift between these two families is tiny (on the order of 0.01 Ω₀) and would likely have gone unnoticed in observations. The observed line widths are typically 20–40 nHz (e.g., Löptien et al. 2018), which corresponds to 0.3–0.6 Ω₀; hence, an observed spectral peak is too broad to separate the two potential peaks. Furthermore, because the deeply seated modes that live in the radiative interior are likely to have lower amplitude at the solar surface due to their evanescence in the convection zone, current photospheric observations are inherently less sensitive to them.

6.1.2. Do We See Radial Overtones?

In Section 5.4, among the retrograde propagating waves, we find two equatorially symmetric bands of power and three antisymmetric bands of power (labeled II–VI in Figure 7 and 10). The lowest-frequency symmetric and antisymmetric bands are consistent with previous HFR mode observations and simulations, but the additional ridges are new. These lesser features are likely latitudinal overtones of these HFR modes. Given that we are seeing the same features across different latitudinal orders (Figure 9), the eigenfunctions of these waves are linear combinations of the spherical harmonics.

While we suspect that these five ridges are all HFR modes, there are still open questions about the physical mechanism that produces them. Jain & Hindman (2023) used an analytic model to study the radial and latitudinal propagation of thermal Rossby waves in an isentropically stratified atmosphere. For modes with latitudinal propagation, they find a new set of retrograde-propagating inertial waves whose eigenfrequencies and latitudinal behavior bear a qualitative resemblance to the HFR modes—see Figure 5 in Jain & Hindman (2023). These retrograde inertial waves are the retrograde branch of the prograde thermal Rossby waves. They only exist for waves that possess latitudinal nodes, and for m = 0, the motions consist of rolls with axes aligned with lines of constant latitude. For m ≠ 0, the modes have fully 3D motions. Notably, these modes have similar latitudinal behavior to the potential HFR modes that we have identified in our simulation, with latitudinal overtones occurring at higher and higher (negative) frequencies without limit.

If the HFR modes are the same type of mode as identified by Jain & Hindman (2023), then the HFR modes are in the thermal Rossby-wave family and are mixed modes. By this we mean that the prograde thermal Rossby-wave branch smoothly transitions to a retrograde HFR mode as the azimuthal order m passes through zero from positive to negative values.

Figure 14 plots the power spectra for the poloidal velocity power, summed over λ = 0–3, averaged across two radial slices near the surface. Each individual m has been individually scaled by the mean power in the displayed range and shown with a logarithmic color table to help boost visibility of the lower-power features. Prograde-propagating waves have been shown with negative values of both the azimuthal order m and frequency ω/Ω₀ to better make mixed-mode relationships clear. The prograde-propagating thermal Rossby waves are again labeled 0, 1, and 2, while the retrograde branches are labeled with Roman numerals. For example, the previously discussed mixed mode (Sections 5.3) is denoted with open circles for the retrograde modes (labeled I) and open squares for the prograde modes, i.e., Roberts modes (labeled 1), with numerical values obtained from Bekki et al. (2022a). The Busse mode (labeled 0) is again marked by closed squares with numerical values from Hindman & Jain (2022).

Zoom In Zoom Out Reset image size

In Figure 14, we note that the antisymmetric HFR branch (II) seems to flow into the prograde thermal Rossby wave with two latitudinal nodes (2), which directly relates the primary HFR observations to the previously noted latitudinal overtone. Additionally, while it is not as visible in this figure, we remind the reader that we identified a potential thermal Rossby wave with high frequency and three latitudinal nodes in Figure 8(b) (see Section 5.2). This overtone was most visible for high-m values, and it has a frequency of opposite sign but suspiciously similar magnitude to numerical HFR mode results. Ridge III in Figure 14 could potentially transition into this higher-order thermal Rossby wave. Ridges IV–VI may exhibit similar patterns, though it is difficult to make a definitive conclusion due to their low amplitudes.

Furthermore, the retrograde inertial wave identified by Bekki et al. (2022b) has already been shown to be a mixed mode that transitions (when prograde propagating) to a thermal Rossby wave, specifically to one with one node in v_ϕ in latitude. The equatorial symmetry of this mixed mode and its frequency suggests that it might be a member of the HFR family. If this is true, the picture becomes simple. There are only two unique families of inertial waves in the convection zone. We have equatorial Rossby waves that are primarily horizontal in motion. Separately, there exists a sequence of latitudinal overtones of a truly 3D mode that is a prograde thermal Rossby wave "mixed" with a retrograde HFR mode, with the Busse mode as the gravest member. The mode found by Bekki et al. (2022b) is just one member of this series.

While this simulation features many types of inertial waves, several varieties are notably missing. We find no evidence of critical-latitude modes (e.g., Gizon et al. 2020, 2021), most likely due to the relatively weak differential rotation in this model. Additionally, we find no evidence of high-latitude modes (e.g., Gizon et al. 2021), also potentially due to the weak differential rotation. Through thermal wind balance, weak differential rotation leads to weak latitudinal entropy gradients, and the high-latitude modes have been demonstrated to reach large amplitude due to baroclinic instability (Bekki et al. 2022b). However, we further acknowledge that the spectral decomposition into spherical harmonics that we perform is not the most robust way to identify high-latitude features.

Finally, we find no evidence of the splitting of Rossby waves into fast and slow MHD Rossby waves by the magnetic field. This absence is most certainly due to the rather weak magnetic fields generated by this particular dynamo, resulting in magnetic splittings that are too small to discriminate here. To confirm that we are in the weak-field limit, we first compare the Alfvén velocity v_A with the rotational velocity ΩR_⊙. Taking a magnetic field magnitude of about 10³ G (see Figure (2b)), we calculate the Alfvén velocity to be on the order of 10³ cm s⁻¹. Because ΩR ∼ 10⁵ cm s⁻¹, the ratio of the two, v_A/ΩR_⊙ ∼ 10⁻², places our simulation squarely within the weak-field limit. Zaqarashvili et al. (2021) calculate the frequency shift between the hydrodynamic mode and the "fast" magnetic mode for a purely toroidal magnetic field to be

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\omega }{{\rm{\Omega }}}=-\displaystyle \frac{{v}_{{\rm{A}}}^{2}m}{2{{\rm{\Omega }}}^{2}{R}^{2}}\left[l(l+1)-2\right].\end{eqnarray} \tag{ 17 }$$

For example, using the prior estimate for the Alfvén speed, v_A ∼ 10³ cm s⁻¹, the fractional shift for (l, m) = (4, 4) is about 10⁻⁵, making it far too small to be visible in our spectra or on the line profiles shown in Figure 12.

In addition to lacking some of the modes that have been observed on the Sun, we see inertial waves in this simulation that have not yet been observed. As seen in Figures 7(c) and (f), retrograde inertial waves and thermal Rossby waves feature as prominently as, or moreso than, equatorial Rossby waves at the surface. Observationally, the opposite occurs. The sectoral equatorial Rossby waves have by for the most power. There are several potential reasons for this discrepancy. While the influence of radial stratification in a convection zone is not well understood, preliminary results suggest that the inertial waves are highly sensitive to that stratification (e.g., Gizon et al. 2021; Bekki et al. 2022b; Bekki 2023). This simulation is far less turbulent than the Sun and therefore has a larger superadiabatic gradient. Thus, the radial cavity of all of our modes is probably spatially shifted from that in the Sun. This can influence mode excitation, damping, and surface amplitude.

Additionally, we do not yet understand the excitation mechanisms for these inertial waves. Though recent work has shown that turbulent convection can excite equatorial Rossby waves (Philidet & Gizon 2023), we would need to ensure that our convection spectrum matches observations to achieve the correct distribution of power among the various modes. However, like all global simulations of convection that we are aware of, our spectrum has far too much large-scale power (for a synopsis of this issue, see Featherstone & Hindman 2016b).