Overstable Convective Modes in a Polytropic Stellar Atmosphere

Within the convection zone of a rotating star, the presence of the Coriolis force stabilizes long-wavelength convective modes. These modes, which would have been unstable if the star lacked rotation, are called overstable convective modes or thermal Rossby waves. We demonstrate that the Sun's rotation rate is sufficiently rapid that the lower half of its convection zone could possess overstable modes. Further, we present an analytic solution for atmospheric waves that reside within a polytropic stratification. We explore in detail the properties of the overstable and unstable wave modes that exist when the polytrope is weakly unstable to convective overturning. Finally, we discuss how the thermal Rossby waves that reside within the convection zone of a star might couple with the prograde branch of the $g$ modes that are trapped within the star's radiative zone. We suggest that such coupling might enhance the photospheric visibility of a subset of the Sun's $g$ modes.


Introduction
The Coriolis force acts to inhibit convection, in both obvious and subtle ways. In addition to suppressing the efficiency of turbulent heat transport (e.g., Julien et al. 1996;Brandenburg et al. 2009), reducing the dominant spatial scale of the convective motions (e.g., Featherstone & Hindman 2016;Vasil et al. 2021), and inducing anisotropy and latitudinal variation in the convective heat flux (e.g., Cowling 1951;Tayler 1973), rotation has long been known to delay the onset of the convective instability, increasing the thermodynamic gradients required for spontaneous convective overturning. In a Boussinesq system this appears as an increase in the critical Rayleigh number that marks the onset of the instability (see Chandrasekhar 1961). In a gravitationally stratified fluid, the local Ledoux criterion (Ledoux 1947) is modified by rotation such that a steeper specific entropy or compositional gradient is required for instability. For example, Cowling (1951) argued that non-axisymmetric motions are locally stable if the square of the buoyancy frequency N 2 exceeds a rotationally dependent threshold (which is negative), (1) In Cowling's criteria, Ω is the star's rotation rate and k Ω and k h are two components of the wavenumber vector of the convective waveform. The component parallel to the rotation axis is k Ω and the horizontal component, perpendicular to gravity, is k h . This criterion makes it clear that certain modes of convection (i.e, certain wavenumbers) can be stabilized by rotation even if the atmosphere is globally unstable by the Ledoux criteria, namely N 2 < 0.
Cowling's formula further suggests that rotation only influences the stability of those convective motions that have variation parallel to the rotation axis, i.e., k Ω = 0. Without such variation, Cowling's criterion reduces to that of Ledoux, N 2 > 0. We now know that Cowling's analysis was incomplete. An entire class of gravito-inertial waves was excluded from the dispersion relation that was used to derive the stability criterion. These waves, called thermal Rossby waves in geophysics and planetary science and referred to as overstable convective modes or low-frequency prograde waves in astrophysics (see Ando 1989;Unno et al. 1989), are essentially convective modes that have been stabilized and made oscillatory by rotation. Stable forms of these waves exist that disobey Cowling's criterion. The existence of such waves requires either curvature of the boundaries to provide a topological β-effect (e.g., Roberts 1968;Busse 1970Busse , 1986 or radial stratification of the density producing a compressional β-effect (e.g., Hide 1966;Glatzmaier & Gilman 1981;Unno et al. 1989;Hindman & Jain 2022). Both cases enable a prograde-propagating vorticity wave because a spinning convective column that is pushed toward the axis of rotation must spin faster either from vortex stretching induced by the shape of the boundaries or by vortex narrowing caused by density stratification.
Such thermal Rossby waves are absent from Cowling's calculation because his analysis was a purely local one that ignored the presence of physical boundaries and did not consider spatial variation in the atmosphere beyond the term involving the buoyancy frequency. Ando (1989) expanded on Cowling's derivation by deriving a local dispersion relation that includes density stratification. From this dispersion relation one can derive the following local stability criterion for overstable convective modes, where k φ is the azimuthal component of the wavenumber, k is the total wavenumber, θ is the colatitude, and H is a length scale that depends on the density scale height H, the buoyancy frequency N , and the gravitational acceleration g, In a stellar convection zone, the buoyancy frequency is small and this scale length is nearly equal to the density scale height H ≈ H. If the wavelength is much shorter that the density scale height (kH 1), the criterion described by Equation (2) reduces to that of Cowling,Equation (1), unless the wavenumber is purely perpendicular to the rotation axis (k Ω = 0). For waves with such an alignment, Ando's correction to Cowling's criterion becomes the dominant term. If one follows Hindman & Jain (2022) and includes the effect of the acoustic cut-off frequency, ω c , the stability condition for thermal Rossby waves near the equator (sin θ ≈ 1) becomes where k c is a cut-off wavenumber that depends solely on the density scale height and its radial derivative, where c is the sound speed. The cut-off wavenumber encapsulates the reflection that can occur when a wave travels vertically into a region where its vertical wavelength becomes long compared to the density height. As such, the inclusion of the cut-off wavenumber is particularly important for the long wavelength waves for which kH 1. In this limit, the stability criteria becomes exceedingly simple. For example, deep within a stellar convection zone where H ≈ H and k 2 c ≈ 1/4H 2 , we find stability for N 2 > −4Ω 2 .
How likely are stable thermal Rossby waves to exist in a star like the Sun? From an examination of Equation (4) it becomes clear that stable waves can exist in a nominally unstable stratification (N 2 < 0) as long as the rotation rate is of the same order as (or larger than) the modulus of the buoyancy frequency. Figure 1 illustrates the square of the buoyancy frequency N 2 as a function of radius within model S, a standard model of the Sun's interior structure (Christensen-Dalsgaard et al. 1996). The square of the Sun's Carrington rotation rate, Ω 2 , is indicated by the horizontal red lines. Within the Sun's convection zone (shaded gray), the square of the buoyancy frequency is negative and thus the atmosphere is globally unstable to convection. However, within the lower half of the Sun's convection zone, the rotation rate becomes larger than the modulus of the buoyancy frequency and stable thermal Rossby waves become possible for long wavelengths. Hindman & Jain (2022) estimated that waves with an azimuthal order m = k φ R less than 30 might be stabilized by the rotation.
Since thermal Rossby waves are low frequency they will have frequencies commensurate with a subset of the g modes that reside in the Sun's radiative interior. Thus, there exists the intriguing possibility that the thermal Rossby waves will be coupled with the prograde branch of the g modes. In β Cephei stars, the coupling of thermal Rossby waves that reside within the convective core with the g modes of the overlying stable envelope has been proposed as a mechanism for the excitation of pulsations (Osaki 1974;Lee & Saio 1986, 1987. But, such coupling has not been explored previously in a solar context. The nonlinear waveforms, often called banana cells, that result when thermal Rossby waves are unstable have been studied extensively in the solar and stellar context through numerical simulations (e.g., Miesch et al. 2000;Hotta et al. 2015;Nelson et al. 2018;Hindman et al. 2020a,b). While a few studies of waves in atmospheres with adiabatic stratification have been performed in the past (Glatzmaier & Gilman 1981;Bekki et al. 2022;Hindman & Jain 2022), far less attention has been paid to the stable form of these waves when the atmosphere is convectively unstable (N 2 < 0). This will be our goal here, to derive the eigenfrequencies and eigenfunctions for the thermal Rossby waves that exist in a weakly unstable stratification. In particular, we explore the analytic solution presented by Hindman & Jain (2022) for thermal Rossby waves within a polytropic atmosphere. That solution is valid irrespective of whether the polytropic stratification is stable, unstable, or neutrally stable. Here we examine the potentially complex eigenfrequencies and eigenfunctions that apply for globally unstable stratifications. Section 2 of this paper provides the governing equation for the thermal Rossby waves within a stratified atmosphere. Section 3 describes the polytropic atmosphere and presents the analytic solutions. Section 4 explores how the thermal Rossby waves that reside within the solar convection zone might couple with the g modes that are trapped within the radiative interior. Finally, in Section 5 we discuss the implications of our results on mode stability and the visibility of g modes.

Atmospheric Waves in a Rotating Star
We seek solutions for sectoral modes that propagate longitudinally and potentially radially. We prohibit propagation in the latitudinal direction and assume that any atmospheric variation in that direction can be ignored. Such a 2D approximation is consistent with waves that are rotationally dominated and satisfy the Taylor-Proudman constraint. Further, we ignore the curvature of the star's isopycnals and assume constant gravity within a local plane-parallel model. Since we are interested in sectoral modes, we place the local Cartesian coordinate system at the star's equator. We orient the axes of this coordinate system such that x points in the longitudinal direction,ẑ points radially and antiparallel to the gravitational acceleration g = −gẑ, andŷ points in the invariant latitudinal direction, parallel to the star's rotation vector, Ω = Ωŷ. Under such conditions, Hindman & Jain (2022) have shown that linear atmospheric waves in the rotating reference frame satisfy the following governing equation, where k z is a height-dependent vertical wavenumber, This local vertical wavenumber depends on the atmospheric profiles of the sound speed c, buoyancy frequency N , acoustic cutoff frequency ω c , and density scale height H. The variable Ψ(z) expresses the radial behavior of the wave's Lagrangian pressure fluctuation δP , scaled by the square root of the atmosphere's mass density ρ 0 , We seek plane-wave solutions in the longitudinal direction where ω is the temporal frequency and k x is the zonal wavenumber. We have adopted the sign convention that for a positive wavenumber, k x > 0, the waves propagate in the prograde direction if the frequency is positive, ω > 0, and in the retrograde direction for negative frequencies, ω < 0. For comparison to waves in spherical geometry, the zonal wavenumber is related to the azimuthal order of the concomitant spherical harmonic m = k x R, where R is the stellar radius. Equation (6) describes the propagation of acoustic waves and gravito-inertial waves within an atmosphere with a general stratification. In the next section we specialize to a polytropic stratification which is particularly relevant to stellar convection zones.

Polytropic Stratification
A polytrope is an atmosphere for which the atmospheric pressure P 0 is related to the mass density ρ 0 through a power-law relation, where the exponent is usually expressed in terms of a polytropic index α, -5 - In such an atmosphere all of the thermodynamic profiles become power-law functions of the height coordinate, z. Hence, the pressure, density, and temperature all vanish at the height z = 0 and the atmosphere exists only within the half-space z < 0. For a polytropic stratification, the governing equation can be rewritten in the form of Whittaker's Equation (for details see Hindman & Jain 2022), where ζ = −2k x z is a dimensionless depth, ν ≡ (α + 2) /2 depends on the stratification, and κ is a constant eigenvalue, In the expression for κ above, γ is the fluid's adiabatic exponent andα is the value of the polytropic index that corresponds to an atmosphere that is neutrally stable to convective overturning,α ≡ (γ − 1) −1 . If α >α the atmosphere is stably stratified and if α <α the atmosphere is unstable to convective motions. The constant κ is the radial eigenvalue of the ODE. Since κ depends on the frequency and zonal wavenumber, once any given eigenvalue is obtained by applying boundary conditions, Equation (11) serves as a global dispersion relation for that radial mode. The term that is colored red in Equation (11) engenders acoustic oscillations and arises from the compressibility of the fluid with a weak correction by the Coriolis force. The blue term is due to buoyancy and is responsible for internal gravity waves. Finally, the orange term results from the Coriolis force and produces inertial waves.
Whittaker's Equation (Abramowitz & Stegun 1964) has two solutions that can be expressed in terms of Kummer's confluent hypergeometric functions of the first and second kind, M and U . The general solution for the Lagrangian pressure fluctuation is therefore a linear combination of these two solutions, with arbitrary constants C M and C U whose ratio is determined by the boundary conditions. We will examine two different sets of boundary condition, one that is appropriate for a semi-infinite domain, z ∈ (−∞, 0 ] and the other for a finite domain with a depth of D, z ∈ [−D, 0 ].

Eigenvalues for a Semi-infinite Domain
For illustrative purposes, we will consider a semi-infinite domain where we impose regularity conditions at the two singular points of Whittaker's equation (i.e., at z = 0 and z → −∞). Of course, a real star does not possess an infinite domain. But, as long as the turning points that correspond to the edges of the wave cavity are far from the physical boundaries, the solution is insensitive to whether we apply the boundary conditions at a finite depth or not. From Whittaker's Equation (10) we can deduce that deep in the polytropic atmosphere (z → −∞ or ζ → ∞) the solutions are evanescent and decay exponentially with the leading behavior, Ψ ∼ exp (k x z). Further, by setting the coefficient in square brackets in Equation (10) to zero, we find that the wave cavity has a lower boundary (or turning point) that predominantly scales inversely with the zonal wavenumber, Hence, for short zonal wavelengths, the wave cavity is confined within the upper portion of the polytrope and the eigenvalues are insensitive to a boundary condition that is applied at a deep but finite depth.
When the previously discussed conditions of regularity are imposed, the eigenvalue κ is restricted to a discrete spectrum of values labelled by their radial order n, which can be any non-negative integer, n ∈ 0, 1, 2, 3, · · ·. Note, the eigenvalue is insensitive to the zonal wavenumber k x because a semi-infinite polytropic atmosphere is self-similar and lacks an imposed scale length (see Hindman & Jain 2022). In addition to the simple analytic expression for the eigenvalue, the eigenfunctions also become more tractable. Regularity at both singular points collapses both Kummer functions into Associated Laguerre Polynomials (Abramowitz & Stegun 1964). The solution for the eigenfunction of the nth radial order is therefore given by where C n is an arbitrary constant and L (a) n is the nth-order Associated Laguerre Polynomial. In sections 3.3 and 3.4 we will discuss the resulting eigenfrequencies and eigenfunctions in conjunction with those appropriate for the finite domain that we will discuss in the next subsection.

Eigenvalues for a Finite Domain
In this subsection we model the convection zone of a low-mass star as a polytropic layer of finite depth D, spanning z ∈ [−D, 0 ]. For boundary conditions, we impose regularity at the origin (a singular point of Whittaker's Equation) and vanishing of the Lagrangian pressure fluctuation at the bottom of the layerwhich for low frequency is consistent with a condition of impenetrability (see Hindman & Jain 2022). Since, the U hypergeometric function is singular at the origin, the eigenfunction only depends on the M confluent hypergeometric function, and the eigenvalue κ n is determined by a transcendental dispersion relation that enforces the boundary condition at the bottom of the layer, We have numerically solved for the roots of this equation for a convection zone depth of D = 200 Mm and a stellar radius equal to that of the Sun, R = 696 Mm. Figure 2 illustrates the resulting eigenvalues κ n as a function of azimuthal order m = k x R for modes with a radial order n less than 9. The solid blue curves show κ n for the finite domain, whereas the black dotted lines illustrate the eigenvalue for the semi-infinite domain, κ n = n + α/2 + 1. The eigenvalues for the finite domain begin with a large value at low zonal wavenumber and decrease with increasing wavenumber, eventually asymptoting to a constant value that corresponds to the eigenvalue that holds for the semi-infinite domain. This asymptotic behavior occurs because the lower turning point moves higher in the atmosphere as the zonal wavenumber increases and once it has passed into the domain (i.e., when z turn > −D), the boundary condition at the bottom of the convection zone becomes increasing unimportant.

Eigenfrequencies for Slow Rotation
Given a numerical value for the eigenvalue κ, Equation (11) can be solved to obtain the corresponding eigenfrequencies. In order to understand the relative importance of each of the terms that appear in the definition of κ, it is useful to define several nondimensional parameters.
The nondimensional wave frequencyω is based on the star's rotational period and A is a coefficient for the term that leads to acoustic waves. The quantity S acts as a convective stability criterion for the stratification. When S > 0 the atmosphere is stable to convective overturning and when S < 0 the atmosphere is unstably stratified. The magnitude of S indicates the importance of buoyancy to rotation. Since the buoyancy frequency, N , diverges at the origin of a polytrope (where the density vanishes), there is always a region in the upper portion of the polytrope where |N | > Ω. Conversely, deep in the atmosphere the buoyancy frequency vanishes as the depth increases and is therefore small compared to the rotation rate, |N | < Ω. The stability parameter S determines the depth of the transition between these two regimes, z tran = −2SR.
The parameter is a nondimensional measure of the rotational speed that characterizes the centrifugal deformation of a star. For a star rotating much slower than its break up rotation rate, is a small parameter. For example, for the Sun = 2 × 10 −5 (Gizon et al. 2016). The smallness of will prove useful, enabling the expansion of the dispersion relations into distinct branches for the acoustic waves and the gravito-inertial waves.
By using these nondimensional numbers, the definition of κ can be rewritten more succinctly Once again, the red, blue, and orange terms are due to compressibility, buoyancy, and the Coriolis force, respectively. This expression, when multiplied byω 2 , is a quartic polynomial in the frequency. Hence, for any given value for the eigenvalue κ n , there are four solutions for the frequency. Furthermore, since the equation has the form of a depressed quartic and the coefficients are all real-valued, two of the roots always have a real frequency and the other two are either real or form a complex conjugate pair. The two roots that are always real correspond to high-frequency acoustic waves and these arise from the red term in Equation (22). The other two solutions (with real or complex frequencies) are low-frequency gravito-inertial waves and they arise jointly from the blue and orange terms in Equation (22).
When is small (slow rotation), the high-frequency and low-frequency branches can be easily separated by expanding the eigenfrequency in powers of the small parameter . The high-frequency solutions are found by expanding in terms of 1/2 ,ω By solving the global dispersion relation (22) order by order, we find the first two terms in the expansion, In dimensional variables, this translates to As expected, the eigenfrequencies for the acoustic waves are always real and hence the acoustic modes are always stable.
The form of Equation (25) is agnostic of the boundary conditions. The eigenfrequencies depend on the boundaries only through the eigenvalue κ. If we employ the solution for κ n that applies for a semi-infinite atmosphere, with n = n + 1, we obtain the well-known eigenfrequencies for p modes in a polytropic atmosphere (Lamb 1945) but with a rotational correction, The mode with radial order n produces the p mode called p n , i.e., n = 0 corresponds to p 1 , n = 1 to p 2 , and so on.
The low-frequency gravito-inertial waves are obtained with an expansion in powers of , which to leading order produces We immediately deduce that the eigenfrequencies are complex if the atmosphere possesses a sufficiently strong unstable stratification, The condition for stability is a global one that depends not only on the stratification, rotation rate, and wavelength, but also on the boundary conditions through the eigenvalue κ. Even for an atmosphere with an unstable stratification (S < 0 or α <α), a mode can be stabilized by rotation if the mode has sufficiently low zonal wavenumber k x and low radial order n (i.e., low m and low κ n ). Equivalently, for a given radial order n, stability parameter S, and rotation rate Ω, the mode is stable for low zonal wavenumbers and unstable for azimuthal orders exceeding a threshold, This behavior is exhibited in Figures 3 and 4. The polytropic atmosphere for each is weakly unstable with S = −5 × 10 −3 or, equivalently, α −α = −5 × 10 −7 . Figure 3 illustrates the complex eigenfrequencies for the low-frequency gravito-inertial waves in a semi-infinite domain (see subsection 3.1). The left-hand panel shows the real part of the eigenfrequencies, while the right-hand panel presents the imaginary part or growth rate. The first four radial orders (n = 0, 1, 2, 3) are shown as colored curves, respectively, black, green, blue, and red. Modes of higher radial order (up through n = 9) are displayed in light blue. There are two gravito-inertial wave solutions and both are prograde propagating with positive frequencies. The higher frequency solution, or the fast gravito-inertial wave, is shown with the solid curves. The slow gravito-inertial wave is displayed with dashed curves. As the wavenumber increases, the slow and fast gravito-inertial waves have real frequencies that slowly converge. At the threshold of convergence (at marginal stability) and beyond, the two modes become a complex conjugate pair with the same real part of their frequencies and oppositely signed imaginary parts (only the positive root is illustrated). When the two modes have purely real frequencies the modes correspond to overstable convective modes. Both the fast and slow waves would have been convectively unstable and non-oscillatory in the absence of rotation. When the mode frequencies are complex, the two solutions correspond to unstable, oscillatory convective modes that travel prograde at the same speed. For moderate values of the azimuthal order m, these unstable modes have a growth rate that is commensurate with their oscillation frequency. Hence, the modes grow in amplitude over a time scale of roughly the rotation period. At marginal stability and beyond, the two gravito-inertial waves have been called thermal Rossby waves in the convection literature (e.g., Busse 1986). Here, we use the term thermal Rossby wave to describe both the unstable and stable solutions. Figure 4 displays the complex eigenfrequencies for the finite domain. The color of the curves and the line styles have the same meaning as in Figure 3. We have included the eigenfrequencies for the semi-infinite domain as the dotted curves. Only the lowest four radial orders possess low-wavenumber overstable modes. For n ≥ 4 the gravito-inertial waves are unstable for all azimuthal orders. The lowest order modes look very similar to those of the semi-infinite domain, except for the fast gravito-inertial wave at the lowest wavenumbers. This occurs because the lower turning point is below the bottom of the radial domain for low wavenumbers. Hence, these modes are sensitive to the boundary condition at the bottom of the layer. Similar behavior for the thermal Rossby wave where the eigenfrequency approaches zero as the zonal wavenumber approaches zero has been seen previously (Glatzmaier & Gilman 1981;Bekki et al. 2022;Hindman & Jain 2022). As was noted for the semi-infinite domain, the growth rate of the unstable modes is comparable to the rotation rate of the star for moderate azimuthal and radial orders.

Eigenfunctions
Oddly, if the boundary conditions lack explicit dependence on the wave frequency, ω, the two acoustic modes and the two gravito-inertial wave modes have identical eigenfunction for the Lagrangian pressure fluctuation for the same zonal wavenumber k x . This arises because the ODE for the Lagrangian pressure fluctuation depends on the eigenfrequency only through the eigenvalue κ, and all four wave modes have the same degenerate eigenvalue. Where the four wave modes differ is the eigenfunctions for the other physical variables. For example, the zonal and radial velocity components, u and w respectively, can be obtained from the Lagrangian pressure fluctuation through the following differential operators (Hindman & Jain 2022), where σ 2 ≡ gk x −2Ωω. Since, these differential operators are frequency dependent, the velocity eigenfunctions for the acoustic waves and the gravito-inertial waves of the same radial order n will differ even though they possess the same eigenvalue κ n .
Further, in the low-frequency limit that is appropriate for the two gravito-inertial modes, the radial velocity, w, is proportional to the reduced Lagrangian Pressure fluctuation w ∝ δP/ρ 0 and the zonal velocity, u, is identical to within a constant amplitude for the fast and slow gravito-inertial waves. This becomes apparent when the low-frequency limit, 1, of Equations (32) and (33) is taken, Finally, we note that even for the unstable modes the Lagrangian pressure fluctuation has a real eigenfunction and from Equations (34) and (35) we can see that the eigenfunctions for the two velocity components are complex, but possess a complex phase that is independent of height. Furthermore, since we have adopted a weakly unstable polytrope with S = −5 × 10 −3 and α −α = −5 × 10 −7 , the eigenfunctions are nearly indistinguishable from the eigenfunctions for a neutrally stable atmosphere (S = 0 and α =α). Hence, instead of providing redundant illustrations here, we refer the reader to Figures 5 and 7 of Hindman & Jain (2022). Figure 5 from Hindman & Jain (2022) shows eigenfunctions for the semi-infinite atmosphere developed in 3.1 and Figure 7 presents eigenfunctions for the finite layer discussed in subsection 3.2

Coupling of Thermal Rossby Waves with Prograde g Modes
In the previous section we applied a perfectly-reflecting boundary condition (i.e., δP = 0) at the bottom of the convection zone. Hence, the thermal Rossby modes were completely confined to the convection zone and any potential coupling to the g modes of the star's stably-stratified radiative interior was neglected. Since, the square of the buoyancy frequency increases dramatically over the thin boundary between the convection zone and radiative interior, we expect that the reflection is indeed almost total and the coupling between g modes of the interior and thermal Rossby modes of the convection zone is extremely weak, except when the thermal Rossby wave and the gravity wave happen to have a common frequency. Hence, in a dispersion diagram we should see two distinct families of dispersion curves and where those curves cross (i.e., have common frequencies) we expect to see avoided crossings.
In order to demonstrate how the avoided crossings might appear, we have developed a simple, illustrative model consisting of a weakly-unstable polytropic layer of finite depth D that overlays an isothermal layer of thickness L. The polytrope represents the Sun's convection zone and the isothermal atmosphere its radiative interior. As before, we adopt D = 200 Mm and S = −5 × 10 −3 . We choose the buoyancy frequency of the isothermal atmosphere such that it is comparable to the buoyancy frequency of the Sun's interior, N = 10 −3 s −1 . For illustrative purposes we set the depth of the isothermal layer L to an unrealistically thin value, L = 500 km. We do this to control the frequency spacing between g modes of nearby radial order. If we were to use a more realistic value, say L = 500 Mm, the separation between g modes would be so tiny that the density of g-mode ridges in the dispersion diagram would obscure the effect that we are seeking to illustrate. The reader should note that the fine spacing between g modes is not the result of our choice of an isothermal atmosphere. Instead, it is a consequence of the smallness of inertial wave frequencies in comparison to the radiative interior's buoyancy frequency, ω N .
For boundary conditions, we require that the solution is regular at the origin, z = 0, and that the vertical derivative of the Lagrangian pressure fluctuation vanishes at the bottom of the isothermal atmosphere, z = −(D +L). This later condition is appropriate if the bottom boundary were to correspond to the radial center of the star. At the interface between the polytrope and isothermal atmosphere (z = −D), we require that the Lagrangian pressure fluctuation and its radial derivative are continuous. These conditions are consistent with the continuity of the pressure fluctuation and the normal velocity component. In the isothermal atmosphere the solutions for Ψ are sinusoidal (see Hindman & Jain 2022) and the boundary condition at the bottom of that layer fixes the phase of the sinusoid. Within the polytrope, the upper boundary condition fixes the solution to be the one proportional to the M Kummer function. The eigenfunction can therefore be written in the form where C iso and C poly are constant amplitudes and K iso is the radial wavenumber within the isothermal atmosphere, The global dispersion relation arises from the two continuity conditions and has the form, where M = dM/dζ is the derivative of the M Kummer function with respect to its third argument (the spatially varying argument). Figure 5a illustrates the eigenfrequencies for this "toy" model. As predicted there are two families of solutions. Superimposed over the family of thermal Rossby waves that we found in the previous section for the finite domain (see Figure 4), one observes a set of g modes that have frequencies that monotonically increase with the zonal wavenumber k x . The thermal Rossby waves that we found before by ignoring the presence of the g mode cavity are indicated by the black dotted curves. Wherever the dispersion curve for a thermal Rossby mode trapped in the convection zone crosses a dispersion curve for an interior g mode an avoided crossing occurs. Since, the jump in the buoyancy frequency between the isothermal atmosphere and the polytrope is strong with a ratio of 7 × 10 8 , the reflection coefficient between the two layers is nearly unity and the avoided crossings are extremely tight. Therefore, they are not easily visible in panel a. Figure 5b provides a zoom-in view of the crossing indicated in panel a by the small square box. As can be seen, the coupling between a g mode and a thermal Rossby waves occurs over a very narrow range of frequencies where the two distinct cavities have a common resonance.

Discussion
We have developed an analytic solution for atmospheric waves of all types that is valid in a compressible, polytropically-stratified atmosphere. In total there are four atmospheric wave solutions: two high-frequency acoustic waves and two low-frequency gravito-inertial waves. Our analytic solution has the form of a longitudinal plane wave multiplied by a radial eigenfunction comprised of Kummer functions. This solution is valid for any value of the polytropic index, independent of whether the polytrope is stable or unstable to convective overturning. We have explicitly illustrated solutions only for a weakly unstable stratification appropriate for a stellar convection zone. In such a stratification, both gravito-inertial waves propagate in the prograde direction. In a neutrally stable atmosphere (N 2 → 0), the slow wave becomes a degenerate, zero-frequency wave that is stationary in the rotating frame and the fast wave persists as a prograde-propagating wave (Hindman & Jain 2022).
Previous studies have explicitly considered low-frequency waves in a neutrally stable atmosphere. Glatzmaier & Gilman (1981) solved for the eigenfrequencies in the anelastic limit for a polytropic layer of finite radial extent using Frobenius expansions to describe the eigenfunctions. Bekki et al. (2022) solved for the linear eigenmodes by numerically solving the spatially discretized fluid equations in spherical geometry. Hindman & Jain (2022) derived the Kummer-function solutions in the limit of neutral stability (N 2 = 0). For all of these previous efforts, the wave solutions correspond to the fast gravito-inertial wave that we have explored here, which becomes a pure inertial wave when the buoyancy frequency is identically zero throughout the radial domain.

Wave Nomenclature
The two gravito-inertial waves that we find here, and which we have called the fast and slow thermal Rossby wave, are completely analogous to the two wave solutions found by Busse (1986) for a Boussinesq fluid in a rotating cylindrical shell. Just as we have assumed here, Busse (1986) ignored variations and motions parallel to the rotation axis. While we have a compressional β effect due to the gravitational stratification of density, Busse (1986) introduced a topological β-effect by allowing the upper and lower caps of the cylindrical domain to be conical with a tilt designed to mimic the effects of a spherical surface. To illustrate the similarities, we explore Busse's model in the dissipationless limit. Busse (1986) derived dispersion relations for linear waves for systems with and without diffusion. For an inviscid fluid without thermal conduction, the two wave solutions found by Busse (1986) possess frequencies that can be rewritten in the following form where k x is the azimuthal wavenumber, k z is the radial wavenumber, τ is a free-fall time, and K is the azimuthal wavenumber that corresponds to the margin of stability, The free-fall time τ and the threshold wavenumber K, depend on the thickness of the cylindrical shell D, the coefficient of thermal expansion α T , the temperature difference between the inner and outer cylindrical surfaces ∆T , and the angle χ between the upper and lower conical surfaces and the equatorial plane. Note, we have written the centrifugal buoyancy term that appears in Busse's expressions in the form of a more traditional gravitational buoyancy.
In Figure 6 we plot the real and imaginary parts of the eigenfrequencies for these two modes for the parameter value KD/π = 0.3. All of these solutions have a purely real frequency for k x ≤ K and a complex frequency for k x > K. The vertical gray line indicates the marginal wavenumber k x = K that separates the stable and unstable solutions. The color of each curve indicates the radial order n of the mode with k z = nπ/D. Note, the similarity between Figures 4 and 6. Each model has two gravito-inertial waves, and both modes are stable for low azimuthal wavenumber and become unstable once the azimuthal wavenumber crosses a threshold value. The similarities arise because the waves in both models are controlled by a β effect, the effect being topological in Busse's and compressional in ours.
The solid curves in Figure 6 indicate the wave with the faster phase speed, i.e., the plus sign in Equation (39). When the mode is stable, Busse (1986) refers to this mode as the "hydrodynamic mode" and notes that it behaves like a Rossby wave. In the limit of rapid rotation or weak instability, Ωτ → ∞, this fast wave satisfies the dispersion relation where the topological β effect is given by β = 4 sin(χ)Ω/D. In our model, the same limit of rapid rotation or weak instability is achieved by considering N 2 /Ω 2 → 0. In Hindman & Jain (2022) we demonstrate that the local dispersion relation for the fast mode reduces to  -Dalsgaard et al. 1996), a standard model of the Sun's internal structure. The ordinate axis is scaled logarithmically for both positive and negative values that have magnitudes greater than 10 −13 s −2 . In the region between the two black dotted lines, N 2 < 10 −13 s −2 , the scaling is linear. The two red horizontal lines indicate plus and minus the square of the Sun's Carrington rotation rate, ±Ω 2 , where Ω = 2.87 × 10 −6 s −1 . In the lower half of the convection zone, the rotation is large enough to stabilize longwavelength convective modes (Ω > |N |); thus, convectively unstable modes are converted into overstable gravito-inertial waves.  The four modes with the lowest radial order n are shown with black, green, blue, and red curves, while all higher orders are indicated with a common light blue color. The two panels present, as functions of the azimuthal order, (a) the real part of the complex frequency and (b) the growth rate or the absolute value of the imaginary part of the complex frequency. At each value of the azimuthal order, there are two gravito-inertial wave solutions. For low azimuthal orders, both of the gravito-inertial waves are stable and prograde propagating with purely real frequency. The two modes have different frequencies and hence different longitudinal phase speeds. The faster of these is shown with solid curves, while the slow waves are drawn with dashed curves. For sufficiently large azimuthal order, the two modes become unstable and correspond to prograde-propagating convective modes.  Figure 3. The dotted curves indicate the same eigenfrequencies for the semi-infinite domain (see Figure 3). Only four of the lowest radial order modes can be stabilized by rotation. All higher-order modes are unstable for all azimuthal orders. The eigenfrequencies for the two types of domains converge for sufficiently large wavenumber.
Fig. 5.-The real part of the complex eigenfrequencies for the two-layer model, consisting of a weakly unstable polytropic convection zone that overlies a stably stratified isothermal radiative interior. (a) Two families of solutions are possible (shown with blue curves, the gravito-inertial waves that reside in the convection zone and the g modes that are confined to the radiative zone. The black dotted lines indicate four of the gravito-inertial waves that we obtained previously when a perfectly reflecting boundary was placed at the bottom of the convection zone. These are the same eigenfrequencies that were illustrated in Figure 4. In the two-layer model, the interface between the convection zone and the radiative interior is partially transmitting and the two wave cavities can communicate. Because of this, when a gravito-inertial wave has the same frequency as a g mode there is an avoided crossing. (b) A zoom-in view of the avoided crossing marked by the small black box in panel a.  Busse (1986). The solid and dashed curves indicate the fast and slow gravito-inertial waves, respectively. Busse (1986) referred to the fast waves as the hydrodynamic modes and the slow waves as the thermal modes. The unstable modes (and marginally stable modes) were collectively called thermal Rossby waves. The colors correspond to radial order n of the mode as indicated in the legend. Only the four lowest radial orders are plotted. These modes are completely analogous to the gravito-inertial waves that we have explored here (see Figure 4). The β effect is topological in Busse's model and compressional in ours. Both lead to prograde propagating waves.