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

Recent observations of Rossby waves and other more exotic forms of inertial oscillations in the Sun’s convection zone have kindled the hope that such waves might be used as a seismic probe of the Sun's interior. Here, we present a 3D numerical simulation in spherical geometry that models the Sun’s convection zone and upper radiative interior. This model features a wide variety of inertial oscillations, including both sectoral and tesseral equatorial Rossby waves, retrograde mixed inertial modes, prograde thermal Rossby waves, the recently observed high-frequency retrograde (HFR) vorticity modes, and what may be latitudinal overtones of these HFR modes. With this model, we demonstrate that sectoral and tesseral Rossby waves are ubiquitous within the radiative interior as well as within the convection zone. We suggest that there are two different Rossby-wave families in this simulation that live in different wave cavities: one in the radiative interior and one in the convection zone. Finally, we suggest that many of the retrograde inertial waves that appear in the convection zone, including the HFR modes, are in fact all related, being latitudinal overtones that are mixed modes with the prograde thermal Rossby waves.


INTRODUCTION
Although Rossby waves and other inertial oscillations have been extensively studied for the past eighty years, their recent observation on the Sun has revived an interest in the solar physics community, particularly in their potential uses for helioseismology.The Rossby waves themselves are likely to be a sensitive diagnostic of the convection zone's superadiabaticity (e.g., Gilman 1987), and the subsequent detection of high-latitude inertial oscillations (Gizon et al. 2021) has the potential to open up seismic exploration of the polar caps where many dynamo processes are likely to operate.
Rossby waves were first observed on the surface of the Sun by Löptien et al. (2018), who observed sectoral modes, i.e., modes with eigenfunctions that lack latitudinal nodes.This initial discovery has since been confirmed using a variety of techniques and data sets (e.g., Alshehhi et al. 2019;Hanasoge & Mandal 2019;Liang et al. 2019;Hanson et al. 2020;Proxauf et al. 2020;Gizon et al. 2021;Hathaway & Upton 2021;Waidele & Zhao 2023).More recently, Triana et al. (2022) noted the possible presence in the observational data of tesseral Rossby modes, which have one or more latitudinal nodes in their eigenfunctions.Other types of inertial modes have also been observed, such as high-latitude and criticallatitude inertial modes (Gizon et al. 2021) and the mysterious high-frequency retrograde (HFR) waves (Hanson et al. 2022), which currently lack a thorough theoretical explanation.
Conversely, there are theoretically expected inertial oscillations that have not yet been observed.For example, thermal Rossby waves are a class of prograde-propagating inertial waves that have long appeared in simulations (e.g., Busse 1970;Hindman et al. 2020;Bekki et al. 2022a;Hindman & Jain 2022, 2023;Jain & Hindman 2023) and in laboratory experiments (e.g., Mason 1975;Busse & Hood 1982;Azouni et al. 1985;Chamberlain & Carrigan 1986;Cordero & Busse 1992;Smith et al. 2014;Lonner et al. 2022), but have yet to be detected observationally in the Sun.Bekki et al. (2022b) additionally predicts the presence of a "mixed mode" comprised of a prograde branch of thermal Rossby waves with one latitudinal node and a retrograde branch of inertial waves with one radial node.This mode is mixed in the sense that the Yanai mode is mixed (e.g.Gill 1982), with differing behavior between the retrograde and prograde branches.The prograde branch of Yanai waves is essentially composed of internal gravity waves, while the retrograde branch is composed of planetary Rossby waves.Bekki et al. (2022b) refers to the retrograde branch of the mixed mode that they have identified as an equatorial Rossby wave with one radial node.However, this mode has significant vertical motion and thus does not behave like a typical equatorial Rossby wave.For this reason, we refer to them as retrograde mixed modes throughout.In this work, we additionally suggest that this mode is just a single member of a family of mixed modes that reside only in convection zones where the stratification is nearly neutrally stable.This family may include the observed HFR modes, which also evince strong vertical motions.
Although Rossby waves have been a mainstay in meteorology since their discovery in the late 1930s (Rossby 1939;Haurwitz 1940), the theoretical work on Rossby waves in a stellar context has been more sporadic.Papaloizou & Pringle (1978) introduced the concept of "r-modes" to the stellar physics community when they noticed Rossby waves amongst the solution set for low-frequency oscillations of rotating stars.Provost et al. (1981), Smeyers et al. (1981), and Saio (1982) studied the eigenfunctions and radial structure of these modes in the case of slow, spatially uniform rotation.Wolff & Blizard (1986) focused on the radial behavior of r-modes in the Sun's interior, concluding that the r-modes exist in two separate cavities: one cavity in the radiative interior and one in the convection zone.The splitting of the Rossby waves into two distinct families arises due to the sign of the Ledoux discriminant; waves in an unstable stratification propagate in a fundamentally different manner than those in a stable stratification (see also Albekioni et al. 2023).
The recent observations of inertial waves in the Sun have prompted a flurry of theoretical work to better understand which modes can exist in the convection zone and how different physical processes impact these modes.The bulk of this work has been in the form of linear eigenmode calculations (e.g., Gizon et al. 2020;Bekki et al. 2022b;Triana et al. 2022;Hindman & Jain 2022, 2023;Bhattacharya & Hanasoge 2023;Bekki 2023), but the work of Bekki et al. (2022a) examined the inertial modes appearing in a nonlinear numerical simulation with a spatial domain that spanned the solar convection zone.Here, we adopt a similar methodology and explore the inertial modes that manifest in a numerical simulation, but our spatial domain includes both a convection zone and a significant portion of the underlying stably stratified radiative interior.
We find a veritable menagerie of waves that have all been self-consistently excited in a solar-like environment.In particular, the model features equatorial Rossby waves in the radiative interior, tesseral equatorial Rossby waves in the convection zone, and high-frequency retrograde (HFR) waves of latitudinal orders that have not been noted previously.
The paper is organized as follows.Section 2 provides details of the simulation and code.Section 3 describes the spectra that make up the bulk of our analysis.Section 4 focuses on the equatorial Rossby waves found in the radiative interior, while Section 5 covers the many types of waves found in the convection zone.Section 6 summarizes our findings and discusses the implications of our work, including the possibility of two families of Rossby waves as well as the classification of HFR modes.

The Rayleigh Code
This numerical model was generated using the Rayleigh convection code (Featherstone et al. 2021;Featherstone & Hindman 2016a;Matsui et al. 2016), which solves the fluid equations in rotating spherical geometry using the pseudospectral 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 m l (θ, ϕ) and a Chebyshev polynomial T k (r).The three spherical coordinates r, θ, and ϕ, are the radius, colatitude, and longitude, respectively, and r, θ, and φ 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 timestep 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 where the primary variables are the velocity v, magnetic field B, and the fluctuations of the gas pressure P ′ and specific entropy density s ′ about the corresponding reference-state profiles.The angular velocity of the rotating reference frame, Ω 0 , is aligned with the axis of the spherical coordinate system.The gravitational acceleration is given by g(r); the viscous, thermal and magnetic diffusivities by ν(r), κ(r), and η(r), respectively; c P is the specific heat capacity at constant pressure; and e ij is the rate-of-strain tensor.The radial profiles ρ(r) and T (r) are the time-independent, spherically symmetric reference-state density and temperature, which satisfy the ideal gas law, and dŝ/dr 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).

The Numerical Experiment
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 ⊙ and r max = 0.947 R ⊙ , where R ⊙ = 6.957 × 10 10 cm is the solar radius.The transition from convective stability to instability is primarily controlled by the entropy gradient of the reference state dŝ/dr 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 de-aliasing of l max = 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 overshoot 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: where σ ≡ 10 −2 erg g −1 K −1 cm −1 , δ ≡ 0.05 R ⊙ , and r 0 = 0.719 R ⊙ .The gravitational acceleration is given by where G is the universal gravitational constant and M ⊙ = 1.989 × 10 33 is the solar mass.The rotation rate is given by Ω 0 /2π = 1370 nHz, and a solar luminosity L ⊙ = 3.85 × 10 33 erg s −1 is driven through the convection zone via the temporally steady internalheating 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 ν(r max ) = κ(r max ) = 5 × 10 12 cm 2 s −1 and η(r max ) = 1.25 × 10 12 cm 2 s −1 , respectively.All diffusivities increase with radius, varying like ∼ ρ−1/2 .At both boundaries, we use stress-free and impenetrable conditions on the velocity, i.e., We impose a fixed conductive flux at the inner and outer spherical boundaries, using 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 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 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.

Simulation Summary
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 Ω(r) = Ω 0 + ⟨v ϕ ⟩/r sin θ, where the angular brackets refer to a combined temporal and zonal average, and Ω 0 is the frame rate.Figure 1 shows nondimensionalized differential rotation profile (Ω − Ω 0 )/Ω 0 .Here and in Figure 2, the pink dash-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 Ω 0 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 Ω 0 .
Figures 2a and 2b 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 eight orders of magnitude, while the horizontal velocities only decrease by four 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 two orders of magnitude in the square of the field components).
Figure 2c, shows the radial behavior of the dynamic Elsasser number as defined by Soderlund et al. (2012): 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) has demonstrated that the magnetic field in this simulation is concentrated in the low-order sphericalharmonic components.The Elsasser number compares the strength of the Coriolis force to that of the Lorentz force.An Elsasser number greater than one indicates that the Lorentz force dominates the dynamics, whereas an Elsasser number of less than one indicates that the Coriolis force dominates.The Elsasser number in this simulation is always quite small (Λ d < 10 −2 ), 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.
3. SUMMARY OF WAVE SPECTRA We illustrate the wave fields in our model through spectra of the radial component of the vorticity.These spectra are obtained with Fourier transforms in time t and spherical harmonic transforms in the horizontal coordinates θ and ϕ.The resulting decomposition of the data is four-dimensional with each spectral component being a function of temporal frequency ω, radius r, harmonic degree l, and azimuthal order m.An individual spherical harmonic possesses a number of nodes in latitude equal to λ = l − |m|, which we call the latitudinal order.The waves with λ = 0, or l = ±m, are the sectoral modes, which have zero nodes in latitude.The tesseral modes have λ > 0 (or l > |m|) and have at least one latitudinal node.Without loss of generality, we often choose to illustrate only positive values of the azimuthal order (m > 0) and adopt the convention that negative frequencies correspond to modes propagating in the retrograde direction and positive frequencies in the prograde direction.
We generate spectra at fifteen different radii throughout the spatial domain, with five samples each in the radiative interior (0.502 R ⊙ -0.65 R ⊙ ), tachocline (0.678 R ⊙ -0.728 R ⊙ ), and convection zone (0.731 R ⊙ -0.935 R ⊙ ).The simulation was run with a maximum spherical harmonic degree of 255, but we only present spectra for m ≤ 80.The Fourier transform was taken over a single contiguous realization of about 200 years in duration (or 8770 rotational periods), resulting in spectra with a Nyquist frequency of 5000 nHz and a frequency resolution of 0.156 nHz.In terms of a dimensionless frequency ω/Ω 0 , this corresponds to a Nyquist frequency of 3.6 and frequency resolution of 1.1 × 10 −4 .
Figure 3 displays power spectra of the radial vorticity, averaged over the radiative interior (left), tachocline (center), and convection zone (right).To more clearly see each type of wave and its behavior, the top row displays the sum of wave modes with even λ = 0, 2, 4, and 6, which are symmetric across the equator, while the bottom row displays the sum of modes with odd λ = 1, 3, and 5, which are anti-symmetric across the equator.Summing the power only over large spatial scales, or low latitudinal orders (λ < 7), removes a significant amount of nonmodal power arising from convective noise.We additionally scale each m individually to make the mode relationships more visible across all length scales.
, where ⟨⟩ denotes a temporal and spherical average over the equilibrated state.The horizontal contributions are orders of magnitude larger than the radial contribution within the radiative interior.(b) Spherically-averaged magnetic energy density, where the fluctuating magnetic field components are defined in the same rms manner as the velocity fluctuations.The horizontal contributions of the magnetic field again dominate the vertical contribution, though this dominance is less than for the velocity field.(c) Dynamic Elsasser number, as defined by Equation ( 13).The Elsasser number is everywhere tiny, so the inertial waves are largely a hydrodynamic phenomenon with magnetism only acting as a small perturbation.
Figure 3. Waves in the radiative interior, tachocline, and convection zone-Power spectra of the radial vorticity as a function of azimuthal order m and dimensionless temporal frequency ω/Ω 0 , averaged over the radiative interior (left), tachocline (center), and convection zone (right).The top row displays the sum of power from modes that are symmetric about the equator with λ = l − m = 0, 2, 4, and 6, and the bottom row displays the sum of power from anti-symmetric modes with λ = l − m = 1, 3, and 5.The power is scaled separately for each m, so that the maximum power at any given m has a value of unity.The images are shown with a logarithmic color table.We note a rich wavefield in both the radiative interior and convection zone.In particular, sectoral and tesseral Rossby waves are present throughout the domain, along with ridges of high-frequency power in the convection zone corresponding to HFR modes.
Because the dynamic range of each plot was chosen to maximize the visibility of select features, we have elected not to display a color bar.Instead, corresponding line profiles for each wave type provide the relevant magnitudes.The resulting spectrum clearly illustrates in a single image the various types of inertial oscillations that are present in each portion of the spatial domain.
We discuss each of these wave types in the following sections.In the radiative interior, we find both sectoral and tesseral modes of equatorial Rossby waves (Section 4).In the tachocline and convection zone, sectoral and tesseral Rossby waves once again appear (Section 5.1), along with thermal Rossby waves (Section 5.2), retrograde mixed modes (Section 5.3), and HFR modes with what are potentially their latitudinal overtones (Section 5.4).

Doppler-Shift
If the modes live in a region rotating at a rate Ω different from our frame rate Ω 0 (i.e., Ω = Ω 0 + ∆Ω), then a wave of azimuthal order m and intrinsic frequency ω ′ will be observed in the reference frame of the simulation (rotating at rate Ω 0 ) as having the Doppler-shifted frequency (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.Since 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.

WAVES IN THE RADIATIVE INTERIOR
In the radiative interior, we observe a rich spectrum of equatorial Rossby waves.Both sectoral and tesseral modes are present, with power distributed relatively widely amongst the distinct modes.We explore the nature of these waves in the following subsections.

Rossby-Wave Spectra
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 4a 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, This dispersion relation is derived for purely hydrodynamic, equatorial Rossby waves in two horizontal dimensions with solid-body rotation at a frame rate of Ω 0 (Haurwitz 1940).
For such a spherically symmetric system, the eigenfunctions for the Rossby waves are pure spherical harmonics.Since the radiative interior of our model is rotating slightly slower than the frame rate (with a difference of ∆Ω ≈ −0.0015 Ω 0 , see Figure 1), we have applied a Doppler shift to the traditional dispersion relation.
In Figure 4a, the sectoral modes are the ridge labeled λ = 0.The tesseral modes have smaller frequencies in magnitude, 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 4b, 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 4b 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 evident in the latitudinally summed spectral images.
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 eight orders of magnitude above the power background, indicating that these equatorial Rossby waves are the most prominent contributor to the horizontal  15) for the labeled pairs (l, m).The green dotted line denotes power contributed by a mode with a different value of m (arising from a weak misalignment of the mean angular momentum vector and the spherical coordinate axis caused by numerical noise.)The modes blur together near zero frequency since their separation falls below the mode linewidths.Given that these peaks rise orders of magnitude above the noise floor, the vast majority of power in the radiative interior is concentrated in equatorial Rossby waves.b) m = 6.The power is taken at a radius of 0.595 R ⊙ and binned down in frequency by a factor of 32.The peaks from left to right correspond to increasing values of λ (or equivalently l).The higher λ modes are clearly present, but when their power is summed together (as in Figure 4b), their frequency separation falls below their linewidths, and they become an indistinguishable blur.The green dotted lines denote power contributed by modes with a different value of m, labeled by (l, m).
velocity components v θ and v ϕ .The rms velocity summed across all modes has a typical speed of 50 cm s −1 , 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 ω/Ω 0 ≈ 5 × 10 −3 ).This power feature is is likely due to the dynamo cycle which has a similar frequency (see Matilsky et al. 2022Matilsky et al. , 2023)).

Latitudinal Eigenfunctions
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 dis-cursions 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 0 (Y m l0 ), but will also contain weak contributions from spherical harmonics Y m l with the same azimuthal order m and nearby values of the harmonic degree l ≈ l 0 .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 m l0 , 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 0 = 3, 4, and 5.It is clear that the power peaks at the harmonic degree l 0 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.

WAVES IN THE CONVECTION ZONE
The right-hand column of Figure 3 reveals a wide variety of different inertial oscillations present in the convection zone.To begin the exploration of these modes, Figure 7 illustrates the power in radial vorticity at three different radii within this region.The top row is the sum over the power in equatorially symmetric modes with λ = 0 and λ = 2, while the bottom row is the sum over the power in the antisymmetric modes λ = 1 and λ = 3. Filled black circles mark equatorial Rossby waves (Section 5.1), numbers 0-2 mark thermal Rossby waves (Section 5.2), roman numeral I marks the retrograde mixed mode (Section 5.3), and roman numerals II-VI mark HFR modes and their latitudinal overtones (Section 5.4).

Equatorial Rossby Waves
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 since 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  (2022).The plus signs denote the observed HFR values from Hanson et al. (2022).A Doppler shift of 0.015 Ω 0 was applied to these four sets of values, as described in the text.This region features a wide variety of different inertial oscillations, some of which have been observed on the Sun and some of which have not (for a summary, see Table 2).
closest to the surface (right-hand panels), the only modes that are clearly visible are the sectoral modes for roughly m ≤ 5 (Figure 7c), λ = 1 tesseral modes for m < 7 (Figure 7f), and the λ = 2 tesseral modes for m < 3 (Figure 7c).Further discussion of mode properties is deferred to the Discussion (Section 6.1).

Thermal Rossby Waves
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 θ + v θ cos θ (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 ra-  3 and 4, respectively.The thermal Rossby waves in (a) appear to be stable with a Q factor of about 10.The prograde mixed modes in panel (b) and the retrograde mixed modes in panel (c) both have the m = 0 mode illustrated, and it is clear that they have identical power profiles with ω → −ω.
dial vorticity each convective column appears anti-symmetric 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 Ω 0 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 Ω 0 .
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 8a 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 8a with black lines.Widths ranging between 0.03 Ω 0 and 0.08 Ω 0 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 Appendix A.
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 ϕ .Since 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 Ω 0 .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), refer to this behavior where a prograde branch and a retrograde branch are connected as mixed modes.
Figure 8b 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 7f, there is a prograde, higher frequency ridge of power (labeled "2") that slopes downwards and asymptotes towards 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 8b, an additional symmetric mode appears as a high-frequency bump in power in m = 6 (near a frequency of ω = 1.3 Ω 0 ), 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 sphericalharmonic 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.

Retrograde Mixed Modes
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Ω 0 ) from Table 2 of Bekki et al. (2022b) are overplotted as empty black circles and match the ridges rather well.
Figure 8c 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 8c, and the values are listed in Table 4.The widths of these peaks are typically 0.05 Ω 0 , 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 linewidth, the fits to the wings and the power background are clearly inadequate.In particular, the wings are slightly asymmetric.Since 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 8c, and it occurs at the same absolute frequency as the corresponding m = 0 mode in Figure 8b.

HFR Modes
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 separate class from equatorial Rossby waves.Using a linear eigensolver for the Boussinesq equations, Triana et al. (2022) found a class of equatorially antisymmetric modes with similar frequencies 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 anti-symmetric modes.
We find evidence for both the symmetric and antisymmetric branches of the HFR modes.In Figure 7, there are three anti-symmetric bands (labeled II, IV, and VI) and two symmetric bands (III, 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 Ω 0 , are anti-symmetric in radial vorticity and found near the solar surface.The underlying anti-symmetric band of power in our spectra is most likely these HFR modes.However, since 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 anti-symmetric 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, VI) that to our knowledge have not been reported previously.These features are also visible in Figure 9, which displays the corresponding power profiles in radial vorticity for m = 9 and λ = 0 through 4.These additional enhancements in power are slight and extremely broad; hence, if modal, have an extremely low Q factor.
Figure 10 2023) (times signs), again all Doppler shifted by 0.015 Ω 0 .We note that Ridges II and III have the same equatorial symmetry and similar frequencies as previous observations and numerical calculations, with particular concordance with the results of Bhattacharya & Hanasoge (2023).
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.Panels (b) and (d) of Figure 10 show m = 9 power profiles for v θ and v r spectra at two different radii.The v θ spectra have been summed over equatorially symmetric (antisymmetric) power λ = 0, 2 (λ = 1, 3).v r exhibits the opposite symmetry, so the symmetric modes are the sum of λ = 1, 3 while the antisymmetric modes are the sum of λ = 0, 2. We note that while v θ is a couple of orders of magnitude greater than v r at the surface, the difference reduces to less than one order of magnitude at the deeper radial slice, which qualitatively matches previous findings (Triana et al. 2022;Bekki 2023) 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.

Impact of Differential Rotation on Convection
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 Ω 0 .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 Figure 9. Power profiles for the HFR modes-Power spectra in radial vorticity for m = 9 at a radius of 0.935R ⊙ for latitudinal orders that are (a) symmetric (λ = 0, 2) and (b) anti-symmetric (λ = 1, 3) across the equator.The power spectrum has been binned down by a factor of 32 in frequency.We clearly see two symmetric and three antisymmetric high-frequency enhancements in power, matching the labeled ridges in Figure 7.These five peaks are probably all in the same family, with higher-frequency power corresponding to modes of higher latitudinal order.convection zone, the equatorial rotation rate rises above the frame rate by as much as 0.036 Ω 0 , 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).
6. DISCUSSION This simulation simultaneously features a wide variety of self-consistently excited inertial oscillations, some of which have been seen previously in observations and models, and some of which have not.Table 2 lists the waves discussed in this work, where they exist in the simulation, if they have been observed, and whether they have been seen in recent models.Equatorial Rossby waves in the radiative interior, tesseral equatorial Rossby waves in the convection zone, and high-frequency retrograde vorticity waves have not been previously noted in a 3D convection simulation such as this one.Because this model has a solar-like rotation profile (fast equator and slow poles) complete with stable and unstable regions and a tachocline, it provides a useful tool to investigate how these waves manifest and where their wave cavities reside.

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), & 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. 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: that 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).
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 higher in magnitude than the purple dashed line (more negative), while the convection-zone profile has its mean at a lower-magnitude frequency (less negative).
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 referring to a powerspectrum-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Ω 0 /l(l + 1).Because ⟨∆ω⟩ is always positive, the convection zone modes always have a lower 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.
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  Critical-latitude modes -Yes 12 Yes 2,3,5 1 Löptien et al. (2018), 2 Bekki et al. (2022b,a), 3 Gizon et al. (2020), 4 Triana et al. (2022), 5 Fournier et al. (2022), 6 Waidele & Zhao (2023), 7 Hindman & Jain (2022), 8 Hindman & Jain (2023), 9 Bhattacharya & Hanasoge (2023), 10 Hanson et al. (2022), 11 Bekki (2023), 12 Gizon   et al. (2021)   Table 2.A summary of inertial waves, where they are found in this simulation-radiative interior (RI) and/or convection zone (CZ)-and the observations and models that have referred to them in the recent literature.).While there is not a clear relationship between latitudinal node number and frequency difference, the mean frequency for the convection zone is always less negative than those in the radiative interior.We can see that for the sectoral modes, this difference increases with m.
the order of 0.01 Ω 0 ) and would likely have gone unnoticed in observations.The observed linewidths are typically 20-40 nHz (e.g., Löptien et al. 2018), which corresponds to 0.3-0.6Ω 0 ; hence, an observed spectral peak is too broad to separate the two potential peaks.Furthermore, since 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.
In addition to splitting the power into two spectral peaks due to potentially distinct stable-zone and convection-zone mode cavities, the power profiles for the equatorial Rossby waves are clearly not Lorentzian and may actually be the superposition of many under-resolved peaks.The spectrum for the (l, m) = (4, 4) mode as seen in the radiative interior (Figure 12) evinces many narrow peaks.To our eyes, these do not look like realization noise spikes as many of the peaks smoothly span multiple nearby frequencies.However, we admit that the stochastic excitation of these modes has only begun to be explored (see Philidet & Gizon 2023), and frequency correlations in the realization noise are poorly understood.
If real (and not noise), the closely spaced, under-resolved peaks could be radial overtones of the equatorial Rossby waves in the radiative interior.Because they exist in a highly stratified, stable environment, these radial overtones would be very closely packed in frequency with a spacing likely to be inversely proportional to the square of the buoyancy frequency (see Vallis 2017;Albekioni et al. 2023).Our nondimensional frequency resolution is 1.1 × 10 −4 , but Wolff & Blizard (1986) estimate the radial splitting between peaks to be as small as 2 × 10 −4 .We do not have the frequency resolution, nor with the present modelling results can we average over sufficient realizations for noise reduction, to make a definitive statement or to distinguish between fine-scale modal structure and realization noise.
The issue of whether radial overtones are present is further complicated by the presence of differential rotation.It is well-known that shear flows permit the existence of additional families of inertial waves (e.g., Kuo 1949;Mack 1976), such as the critical-latitude modes (Gizon et al. 2021(Gizon et al. , 2020)).These additional modes also are expected to have a dense spectra with peaks that blend together (Philidet & Gizon 2023) and without detailed information about the eigenfunctions may be very difficult to disentangle from radial overtones.

HFR Modes May Be the Retrograde Branch of Thermal Rossby Waves
In Section 5.4, amongst the retrograde propagating waves, we find two equatorially symmetric bands of power and three anti-symmetric bands of power (labelled II-VI in Figures 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.Prograde-propagating waves have been shown with negative values of both the azimuthal order m and frequency ω/Ω 0 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 (Section 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).
In Figure 14, we note that the anti-symmetric 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 8b (see Section 5.2).This overtone was most visi-ble for high-m values, and 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 true, the picture becomes simple.There are only two unique families of inertial waves in the convection zone.We have equatorial Rossby waves which 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.

Absent Inertial Waves
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. 2020Gizon et al. , 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, as discussed in Sections 2 and 4, there is 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, and the resulting magnetic splittings are too small to discriminate (see Figure 2c).

CONCLUSION
In this work, we presented a 3D numerical simulation of the upper radiative interior and convection zone that features a wide variety of inertial oscillations.We summarize our primary findings as follows: 1. Rossby waves in the radiative interior: We presented a rich wave-field of equatorial Rossby waves in the radiative interior and demonstrated that they account for most of the horizontal motion in that region.Their presence raises questions about their impact on dynamo and transport properties in a region of the Sun that has historically been considered rather quiescent.
2. Rossby-wave cavities: We found both sectoral and tesseral equatorial Rossby waves throughout the convection zone that seem to occur at frequencies distinct from those in the radiative interior.This implies that there are potentially two unique families of equatorial Rossby waves living in separate radiative-interior and convection-zone cavities.The presence of both sectoral and tesseral equatorial Rossby waves of varying node number throughout the domain bodes well for their potential helioseismic utility, but it does raise some questions.We do not yet know how these modes are excited and whether this mechanism differs between the radiative interior and convection zone.We also need to develop a better understanding of the conditions under which these waves live and how they interact with each other.Stratification and superadiabaticity in particular may play a significant role in controlling the mode frequencies and determining the amplitude of any given radial overtones.
3. HFR modes and thermal Rossby waves: We noted the presence of both previously observed and modeled branches of HFR modes along with three additional branches that are probably their latitudinal overtones.These modes seem to be mixed with thermal Rossby waves.If true, the theoretical picture becomes simple and unified.The Busse modes, with no nodes in latitude, are the fundamental mode.The previously noted mixed mode, which relates Roberts modes to a 3D retrograde oscillation, is the first latitudinal overtone.The anti-symmetric HFR modes, mixed with prograde thermal Rossby waves with two nodes in latitude, are the second latitudinal overtone, and on.

Figure 1 .
Figure 1.Rotation profile achieved in our simulation-This model features a differentially rotating convection zone and rigidly rotating radiative interior.(a) The azimuthally and temporally averaged rotation profile measured relative to the frame rate shown in a meridional cross-section.The north pole of the coordinate system is vertical on the page.The dashed black semi-circle indicates the base of the convection zone and the dash-dotted magenta semi-circle shows the bottom of the tachocline.(b) Radial cuts at several latitudes (indicated by labels) through the azimuthally averaged profile.The inset shows a magnified view of the rotation rate in the radiative interior, which has about 0.0015 Ω 0 of latitudinal differential rotation.

Figure 2 .
Figure 2. Energy densities achieved in the simulation-The velocity and magnetic field are primarily horizontal within the radiative interior.(a) Spherically-averaged fluctuating kinetic energy density separated into the contributions from each velocity component.Fluctuating velocities are defined by v′ i ≡ ⟨(v i − ⟨v i ⟩) 2 ⟩ 1/2, where ⟨⟩ denotes a temporal and spherical average over the equilibrated state.The horizontal contributions are orders of magnitude larger than the radial contribution within the radiative interior.(b) Spherically-averaged magnetic energy density, where the fluctuating magnetic field components are defined in the same rms manner as the velocity fluctuations.The horizontal contributions of the magnetic field again dominate the vertical contribution, though this dominance is less than for the velocity field.(c) Dynamic Elsasser number, as defined by Equation (13).The Elsasser number is everywhere tiny, so the inertial waves are largely a hydrodynamic phenomenon with magnetism only acting as a small perturbation.

Figure 4 .
Figure 4. Equatorial Rossby waves in the radiative interior-(a) Power spectra in radial vorticity as a function of azimuthal order m and dimensionless temporal frequency ω/Ω 0 , averaged over the radiative interior and summed over all latitudinal orders.The λ = 0 and λ = 1 ridges are labeled.The power is scaled separately for each m, so that the maximum power at any given m has a value of unity, and shown with a logarithmic color table.The black dots indicate the theoretical values given by Equation (15).A large number of both sectoral and tesseral equatorial Rossby modes are present throughout the radiative interior.(b) m = 2 power profile at radius 0.595 R ⊙ , summed over all λ, binned down in frequency by a factor of 32.The dashed black lines mark the theoretical values given by Equation (15) for the labeled pairs (l, m).The green dotted line denotes power contributed by a mode with a different value of m (arising from a weak misalignment of the mean angular momentum vector and the spherical coordinate axis caused by numerical noise.)The modes blur together near zero frequency since their separation falls below the mode linewidths.Given that these peaks rise orders of magnitude above the noise floor, the vast majority of power in the radiative interior is concentrated in equatorial Rossby waves.

Figure 5 .
Figure 5. Rossby-wave power profiles-Power spectra of the radial vorticity as a function of dimensionless temporal frequency ω/Ω 0 for modes with 0 ≤ λ ≤ 10 and (a) m = 2 and (b) m = 6.The power is taken at a radius of 0.595 R ⊙ and binned down in frequency by a factor of 32.The peaks from left to right correspond to increasing values of λ (or equivalently l).The higher λ modes are clearly present, but when their power is summed together (as in Figure4b), their frequency separation falls below their linewidths, and they become an indistinguishable blur.The green dotted lines denote power contributed by modes with a different value of m, labeled by (l, m).

Figure 6 .
Figure 6.Eigenfunction perturbations from asphericity-Power for individual equatorial Rossby wave modes as a function of harmonic degree, demonstrating that the mode eigenfunctions are not a single spherical harmonic.(a) Radial vorticity power spectra for m = 2 in the radiative interior.(b) Power as a function of harmonic degree l, summed over a narrow frequency window.The central frequency of the band is found by computing the mean of the power distribution, and the width of the window is 10σ, where σ is the standard deviation.The power is clearly concentrated in a single spherical harmonic component with weak spread across nearby spherical harmonic components lying within ±2 of the principal component.

Figure 7 .
Figure 7. Waves in the convection zone-Power in radial vorticity as a function of azimuthal order m and dimensionless temporal frequency ω/Ω 0 at three different radii within the convection zone.The top row displays the sum of symmetric modes λ = 0 and 2, and the bottom row displays the sum of anti-symmetric modes λ = 1 and 3.The power is scaled separately for each m so that the maximum power at any given m has a value of unity.The images are shown on a logarithmic color table.The filled black circles are the expected values for equatorial Rossby waves from Equation (15) with no Doppler shift applied.Prograde thermal Rossby-wave ridges are numbered via their latitudinal node number in v ϕ , i.e., "0" indicates the Busse modes (no latitudinal nodes), "1" indicates the Roberts modes (one node), "2" indicates modes with three latitudinally stacked convective columns (two nodes), and so on.Roman numerals label retrograde inertial waves, with I denoting the retrograde mixed mode and II-VI denoting HFR modes.The empty black circles and the empty black squares show the numerical frequencies fromBekki et al.  (2022b).Filled black squares denote the semi-analytic dispersion relation for thermal Rossby waves fromHindman & Jain (2022).The plus signs denote the observed HFR values fromHanson et al. (2022).A Doppler shift of 0.015 Ω 0 was applied to these four sets of values, as described in the text.This region features a wide variety of different inertial oscillations, some of which have been observed on the Sun and some of which have not (for a summary, see Table2).

Figure 8 .
Figure 8. Power profiles for thermal Rossby waves and mixed modes-Selected power profiles from radial-vorticity power spectra, binned down in frequency by a factor of 32.(a) Busse modes, shown for λ = 1 at a radius of 0.889 R ⊙ .(b) Roberts modes, shown for λ = 2, at a radius of 0.935 R ⊙ .The power in the m = 0 mode is multiplied by a factor of two for clarity.The bump labeled "3" may indicate the presence of a thermal Rossby wave with three nodes in latitude in v ϕ .(c) The retrograde branch of the mixed mode of Bekki et al. (2022b) shown for λ = 0 for m = 2, 4, and 6 and in λ = 2 for m = 0, at a radius of 0.935 R ⊙ .The power in the m = 2 mode is multiplied by a constant factor of 3 for clarity.The (l, m) = (2, 0) mode, multiplied by a constant factor of 20, is plotted to show its relation to the m = 0 mode in (b).Lorentzian fits have been overplotted in black lines in Panels (a) and (c), and the corresponding values are listed in Tables3 and 4, respectively.The thermal Rossby waves in (a) appear to be stable with a Q factor of about 10.The prograde mixed modes in panel (b) and the retrograde mixed modes in panel (c) both have the m = 0 mode illustrated, and it is clear that they have identical power profiles with ω → −ω.
displays power spectra for the poloidal velocity |v θ | 2 + |v r | 2 , summed over two different radii, with the five HFR ridges again marked by roman numerals II-VI and I denoting the retrograde mixed mode.Because the equatorial symmetries of v θ and v r are opposite, the modes are denoted by the symmetry consistent with the radial vorticity; the "symmetric modes" are the sum of |v θ | 2 over λ = 0, 2 and |v r | 2 over λ = 1, 3, and the antisymmetric modes are the converse.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 (

Figure 10 .
Figure 10.HFR poloidal-velocity power-(a,c) Poloidal-velocity power spectra summed over two radii near the surface, presented for (a) equatorially symmetric modes and (c) equatorially antisymmetric modes consistent with the symmetries in radial vorticity.v θ shares the same symmetries as radial vorticity, but v r has the opposite.To make the "symmetric" spectra, we thus sum over λ = 0 and 2 in |v θ | 2 and λ = 1 and 3 in |v r | 2 .For the "antisymmetric" spectra we sum over the converse pairs.The retrograde mixed mode is again labeled I, with open circles denoting numerical values from Bekki et al. (2022b).Equatorial Rossby waves are marked with closed circles.HFR mode observations from Hanson et al. (2022) (plus signs) and numerical eigensolver results from Bhattacharya & Hanasoge (2023) (symmetric modes marked in stars and antisymmetric modes marked in times signs) and Triana et al. (2022) (triangles) are overplotted.Five identified HFR ridges are marked with roman numerals II-VI.(b, d) m = 9 power profiles for v θ and v r power spectra at the two different radii, with the corresponding HFR modes marked with roman numerals.(b) shows v θ summed over λ = 0 and 2 and v r summed over λ = 1 and 3. (d) shows the converse.While v θ is greater at both radii, the difference between v θ and v r decreases deeper in the interior.

Figure 11 .
Figure 11.Convective features in spectra-Radial-vorticity spectra for λ = 0 at three distinct radii in the convection zone: (a) 0.728 R ⊙ , (b) 0.835 R ⊙ , and (c) 0.935 R ⊙ .Each column has been normalized by the largest value for that m and shown on a logarithmic scale.The black dashed line shows the Doppler shift due to the radial differential rotation, given by m∆Ω.∆Ω is the average differential rotation in latitude at each depth and is valued at (a) 0.007 Ω 0 , (b) 0.015 Ω 0 , and (c) 0.036 Ω 0 .At large m, the convection is Doppler-shifted due to the differential rotation in the convection zone.

Figure 12 .
Figure12.Radial variation of the power distribution of a Rossby wave-Line profiles for power spectra of the radial vorticity for the single spherical harmonic component (l, m) = (4, 4).The blue spectrum is at a radius of 0.502 R ⊙ located at the base of our radiative interior, and the orange is at 0.776 R ⊙ within the lower portion of the convection zone.The purple line is the value Equation (15) gives when ∆Ω = 0, or what we would see if these waves lived in a region rotating at the frame rate of the simulation.We can clearly see that the equatorial Rossby waves in the radiative interior have a lower central frequency than those in the convection zone, potentially indicating two different wave families residing in two different cavities.

Figure 13 .
Figure13.Frequency differences-The difference in mean frequency ⟨ω⟩ = ⟨ω⟩ cz −⟨ω⟩ ri between modes averaged over the radiative interior and modes averaged over the convection zone.⟨∆ω⟩ has been non-dimensionalized by the dispersion relation value ω 2D = 2mΩ 0 /l(l + 1).While there is not a clear relationship between latitudinal node number and frequency difference, the mean frequency for the convection zone is always less negative than those in the radiative interior.We can see that for the sectoral modes, this difference increases with m.

Figure 14 .
Figure 14.Mixed modes-Power spectra for the poloidal velocity power, summed over λ = 0-3, averaged across two radial slices near the surface.Prograde-propagating waves have been shown with negative values of azimuthal order m and negative frequencies ω/Ω 0 .The Busse modes, labeled "0", are marked with values from Hindman & Jain (2022).Numerical values from Bekki et al. (2022b) mark the previously noted mixed modes, with the prograde Roberts modes (open squares), labeled "1", and the retrograde branch (open circles) labeled "I".Sectoral equatorial Rossby waves are marked with filled circles.Ridge II, which corresponds with the observed HFR modes, appears to cross m = 0 and merge with the prograde thermal Rossby wave with two latitudinal nodes in v ϕ , labeled "2".This potentially indicates that all of these waves are mixed modes in the same family, with thermal Rossby waves making up the prograde branch and HFR modes making up the retrograde branch.

Table 1
Non-dimensional Fluid Parameters

Table 1 .
The flux Rayleigh number RaF, Ekman number Ek, Prandtl number Pr, magnetic Prandtl number Prm, convective Rossby number Roc, 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.

Table 2
Summary of Observed and Modeled Solar Inertial Modes