Low-frequency Internal Gravity Waves Are Pseudo-incompressible

Starting from the fully compressible fluid equations in a plane-parallel atmosphere, we demonstrate that linear internal gravity waves are naturally pseudo-incompressible in the limit that the wave frequency ω is much less than that of surface gravity waves, i.e., ω≪gkh , where g is the gravitational acceleration and k h is the horizontal wavenumber. We accomplish this by performing a formal expansion of the wave functions and the local dispersion relation in terms of a dimensionless frequency ε=ω/gkh . Further, we show that, in this same low-frequency limit, several forms of the anelastic approximation, including the Lantz–Braginsky–Roberts formulation, poorly reproduce the correct behavior of internal gravity waves. The pseudo-incompressible approximation is achieved by assuming that Eulerian fluctuations of the pressure are small in the continuity equation—whereas, in the anelastic approximation, Eulerian density fluctuations are ignored. In an adiabatic stratification, such as occurs in a convection zone, the two approximations become identical. However, in a stable stratification, the differences between the two approximations are stark and only the pseudo-incompressible approximation remains valid.


INTRODUCTION
Numerical simulations of convection in low-mass stars, the Earth's atmosphere, giant planets, and many other astrophysical objects all must face the tyranny of sound.Generally, sound waves propagate quickly and have high frequencies; thus, the typical timescale associated with acoustics is far shorter than those arising from convection and large-scale circulations.In a numerical simulation, this short timescale ensures through the CFL condition that sound waves control the size of the timestep that can be taken while still maintaining numerical stability.The difference can be dramatic.For example, at the base of the Sun's convection zone, the speed of sound is roughly 200 km s −1 while the convective flow speed is on the order of 20 m s −1 (e.g.Miesch et al. 2012).A numerical simulation that is forced to track sound waves for stability will need to take 10 4 times as many time steps to evolve the solution for the same duration as a simulation that could ignore the acoustic wave field.This inflation of the necessary computational work is particular onerous since the immense timescale difference between the deep convection and the sound waves indicates that the two phenomena are essentially decoupled.
A variety of methods have been proposed to mitigate this dilemma; almost all involve modifications to the fluid equations to either temper the impact of sound waves or to remove sound altogether.One way to reduce the influence of sound on the time step is to artificially lower the speed at which sound waves propagate (e.g., Rempel 2005Rempel , 2006;;Hotta et al. 2012;Käpylä et al. 2016;Iijima et al. 2019).Successful application of such Reduced Speed of Sound Techniques (RSST) requires that the sound speed be reduced sufficiently to make sound waves tractable, but to maintain enough celerity in the sound waves such that they do not interact strongly with the convective motions.
A more common solution is to surgically remove terms from the continuity equation such that sound waves are no longer a permissible solution to the fluid equations.These "sound-proofed" equation sets typically apply to low-Mach number motions with small thermodynamic fluctuations about a hydrostatic background atmosphere.The most venerable of these techniques is the Boussinesq approximation, whereby the fluid is assumed to be incompressible with constant density.In the highly stratified atmospheres of stars and giant planets where the mass density can vary by orders of magnitude, treatments that can account for the stratification are necessary.In these stratified systems, the fundamental presumption is that for sedate motions a displaced parcel of fluid quickly equilibrates thermodynamically with its new surroundings.In astrophysics the most common of these extensions to the Boussinesq framework is the anelastic approximation (e.g., Batchelor 1953;Ogura & Phillips 1962;Gough 1969;Gilman & Glatzmaier 1981;Bannon 1996), which removes all density fluctuations that appear in the continuity equation.A similar technique called the pseudo-incompressible approximation is a bit subtler, removing only the influence of Eulerian pressure fluctuations from the continuity equation (e.g., Durran 1989;Klein 2009;Vasil et al. 2013).
Such sound-proofing techniques have been used extensively in stellar and planetary convection simulations where the convecting layer spans many density scale heights.In regions of efficient convection, where the redistribution of heat and mass by the convective motions efficiently drives the atmosphere towards an adiabatic stratification, the most common forms of the anelastic and pseudo-incompressible equations are identical and either approximation works well.However, in a stably stratified fluid, the two approximations differ to the extent that they may violate their underlying assumptions, leading to different dynamics.Specifically, Klein et al. (2010), Brown et al. (2012) and Vasil et al. (2013) have demonstrated that anelastic formulations do a disservice to internal gravity waves leading to a loss of energy conservation and to large errors in the wave frequencies.Further, Klein et al. (2010) and Vasil et al. (2013) have demonstrated that although the pseudo-incompressible approximation does far better in preserving the properties of internal gravity waves, it too evinces discrepancies from the fully compressible wave forms.
Here, we demonstrate that internal gravity waves naturally approach the pseudo-incompressible limit as their frequency becomes very low.The discrepancies noted by Klein et al. (2010) and Vasil et al. (2013) arise only when the wave frequencies become large and the assumption of sedate motions in a state of pressure-equilibrium is lost.We accomplish this by deriving internal gravity waves in a plane-parallel atmosphere with a general stratification and subsequently performing a low-frequency expansion of the local dispersion relation and of the wave functions.We find that, to lowest-order in the frequency, internal gravity waves are incompressive.To the next order in the frequency, they become pseudo-incompressible.All forms of the anelastic approximation fail to produce the correct behavior for both the dispersion relation and the wave functions.
In the next section we formulate the anelastic and pseudo-incompressible approximations.Section 3 derives the governing equation for internal gravity waves in a general stratification for a fully compressible fluid.We explore the low-frequency limit of these waves in Section 4, deriving the magnitude and ordering of terms in the continuity and momentum equations.In Section 5 we rederive internal gravity waves using three different sound-proofed equation sets and discuss the integrity of each approximation.Finally, in Section 6 we summarize and discuss the implications of our results.

The Anelastic Approximation
The anelastic condition is a relatively simple replacement for the continuity equation that captures significant density variation in the mean properties of the fluid.For instance, in a gravitationally stratified fluid with velocity, u, and time-averaged density that varies with height, ρ 0 (z), the continuity equation is replaced with (1) This expression can be derived from the full continuity equation, by making two assumptions that are often appropriate for flows of low Mach number: 1) the time derivative of the mass density ρ is inconsequential and 2) the fractional fluctuations of the density around the background density are small, i.e., |ρ 1 /ρ 0 | ≪ 1 where ρ = ρ 0 + ρ 1 .
The popularity of the anelastic approximation arises from two important properties.When the continuity equation is replaced by the anelastic condition, Equation (1), sound waves are removed as a permissible solution to the fluid equations and the mass flux ρ 0 u can be written using stream functions.Brown et al. (2012) and Vasil et al. (2013) both remarked that when the anelastic form of the continuity equation is employed, the fluid equations are no longer energy conserving without modifications to the momentum equation.To enforce conservation of energy, an otherwise unmotivated change to the buoyancy force is required.For an inviscid fluid, the vertical momentum equation can be written in the following form, with the pressure P and specific entropy density s decomposed into a steady hydrostatic background and a fluctuation about that background, P = P 0 + P 1 and s = s 0 + s 1 .The vertical velocity is w, c p is the specific heat capacity at constant pressure, and z is the height within the atmosphere with concomitant unit vector ẑ anti-aligned with gravity, g = −g ẑ.Further, the quantity N 2 = gc −1 p ds 0 /dz is the square of the atmosphere's buoyancy or Brunt-Väisälä frequency.In Equation (3), we have ignored the density fluctuation in the inertial term on the left-hand side, subtracted the steady hydrostatic component from the force balance, and used the ideal gas law to rewrite the density fluctuation in terms of the pressure and entropy fluctuations.To ensure energy conservation, the term involving the buoyancy frequency must be discarded or be physically subdominant.In a convection zone, where efficient heat transport drives the atmosphere towards an adiabatic gradient with N 2 ≈ 0, this approximation is completely justified and has been coined the Lantz-Braginsky-Roberts (LBR) formulation of the anelastic approximation (Lantz 1992;Braginsky & Roberts 1995).Conversely, in a stably stratified region, the term is not small and cannot generally be self-consistently ignored.
We will examine two distinct formulations of the anelastic approximation.Both replace the continuity equation with the anelastic condition (1).One of these approximations-which we will dub the "fiducial" anelastic approximation-will make no further assumptions, leaving the momentum equation unmodified.The other formulation will be the LBR anelastic approximation as discussed above, which ensures energy conservation by excising a specific term from the momentum equation.

The Pseudo-Incompressible Approximation
The pseudo-incompressible approximation as proposed by Durran (1989) modifies the continuity equation under the assumption that Eulerian fluctuations of the gas pressure can be ignored.Following Durran (2008), we start by defining the potential density ρ * for an ideal gas, ρ * ≡ ρ e s/cp . (4) If we take the convective derivative of the potential density and utilize the continuity equation ( 2) and the thermal energy equation, we obtain a prognostic equation for the potential density where T is the temperature and Q represents all irreversible thermodynamic processes, such as thermal diffusion, viscous heating, radiative transfer, etc.Finally, by invoking Equation ( 4) and the equation of state for an ideal gas, we replace the time derivative of the potential density with the time derivative of the gas pressure, In the preceding equations, γ is the gas's adiabatic exponent and c is the sound speed given by c 2 = γP/ρ.Equation ( 8) is an exact form of the continuity equation for which no approximation has been made other than the gas being ideal.The pseudo-incompressible approximation is achieved by assuming that the term involving the time derivative of the gas pressure is negligible, Such an approximation is valid in the limit of infinite sound speed and is consistent with slow motions of low Mach number for which a displaced parcel of fluid rapidly reaches pressure equilibration with its new surroundings.Most importantly, making this approximation removes sound waves from the fluid equations in the same way that anelasticity does.Durran's form of the pseudo-incompressible approximation (Durran 2008) involves replacing the continuity equation by the preceding equation, but otherwise leaving the other fluid equations unmodified-specifically, the momentum equation remains the same.For isentropic motion, the pseudo-incompressible condition reduces to a form that is reminiscent of the anelastic relation with the mass density replaced by the potential density.However, for flows with low Mach number, thermodynamic fluctuations are small and we can safely linearize Equation ( 10), replacing the potential density by the potential density of the hydrostatic background atmosphere (denoted by '0' subscripts), The last equivalency in Equation (11) arises by noting that the potential density is the density that a fluid parcel would possess if displaced adiabatically to a fiducial height in the atmosphere where P 0 = P , ρ 0 = ρ, and s 0 = 0. Like the anelastic approximation, the flow field can be expressed using streamfunctions when Equations (10) and ( 11) are valid, We remind the reader that these two equations were derived using two assumptions: 1) the advective time scales are fast compared to diffusion times-i.e., isentropic motion, and 2) thermodynamic fluctuations are small compared to the background atmosphere.

INTERNAL GRAVITY WAVES IN A GENERAL STRATIFICATION
Consider a plane-parallel atmosphere with a gas pressure P 0 and mass density ρ 0 related through hydrostatic balance, dP 0 /dz = −gρ 0 .Further, let the thermal structure of the atmosphere be general and specified by the vertical variation of the specific entropy density, s 0 .We start with the linearized fluid equations for a fullycompressible ideal gas, We have ignored rotation, magnetism, and all dissipative mechanisms, including viscosity, thermal conduction, and radiative transfer.The thermodynamic variables s 1 , ρ 1 , and P 1 are the Eulerian fluctuations of the specific entropy density, the mass density, and the gas pressure respectively.
Since gravity provides the only preferred direction, internal gravity waves can be treated as a 2D phenomenon that propagates vertically and in a single horizontal direction.Let ẑ be the unit vector that is antiparallel to the constant gravitational acceleration, g = −g ẑ.Further, let x be the horizontal unit vector that is aligned with the wave's horizontal direction of propagation.Finally, seek plane-wave solutions with the form where k h is the horizontal wavenumber, ω is the temporal frequency, and f (z) is a vertical wave function.
The transformed set of equations can be manipulated to express the velocity and its divergence solely in terms of the Lagrangian pressure fluctuation, δP .The resulting equations are a coupled system of ODEs, with the vertical coordinate z as the independent variable and u and w being the horizontal and vertical velocity components, u = ux + wẑ.The Lagrangian pressure fluctuation is related to the Eulerian pressure fluctuation and the vertical velocity, The denominator of Equations ( 18) and ( 19) is spatially constant and will appear later.Therefore for convenience we make the definition, Equations ( 18)-( 20) can be combined to produce a single stand-alone ODE with δP as the dependent variable, where N is the buoyancy frequency and H is the density scale height, In Equation ( 23), the term that involves the sound speed is responsible for the propagation of high-frequency acoustic waves and the term with the buoyancy frequency leads to internal gravity waves.As we will see in the following subsection, the first-derivative term ensures energy conservation for both varieties of wave.
Once one has solved for the Lagrangian pressure fluctuation by applying boundary conditions to Equation (23), the velocity components, u and w, can be found directly through the use of Equations ( 18) and ( 19).Subsequently, all of the thermodynamic fluctuations can then be derived through Equations ( 14), ( 16), and (21), All of the thermodynamic fluctuations appear as linear combinations of the two velocity components.

Energy Conservation and the First Derivative
Here we demonstrate that any viable sound-proofing technique must produce an appropriate coefficient for the first-derivative term that appears in Equation ( 23).This term is crucial for energy conservation.To see this, consider the vertical energy flux for an acoustic-gravity wave, F (z) = ⟨w P 1 ⟩, where angular brackets <> indicate a temporal average over a wave period.Since, the second term on the right-hand side of Equation ( 21) is 90 degrees out of phase with the vertical velocity, in a time average the second term's contribution vanishes and the energy flux can be written just in terms of the Lagrangian pressure fluctuation, where the superscript asterisks denote complex conjugation.By employing Equation ( 19), one can demonstrate that this flux is inversely proportional to the mass density and proportional to the Wronskian of the Lagrangian pressure fluctuation and its complex conjugate, Abel's Identity tells us that to within an unknown multiplicative constant, C, the Wronskian depends only on the coefficient of the first derivative term in the ODE.For the ODE here, the necessary integration is trivial to perform, Hence, the energy flux is constant with height even though the coefficients of the ODE are vertically variable, The constancy of the energy flux with height in the atmosphere is one way to characterize the conservation of energy by acoustic-gravity waves.
From this analysis, we can deduce that any approximation that incorrectly reproduces the first derivative term, may produce wave solutions with energy fluxes that vary with height.Consequently, such approximations will fail to conserve energy.For example, if the first derivative term is artificially set to zero, the flux will be inversely proportional to the mass density and F (z) will spuriously increase with height.This is the fundamental reason why Brown et al. (2012) and Vasil et al. (2013) found a lack of energy conservation when applying a variety of anelastic approximations to an isothermal atmosphere.Those approximations failed to correctly reproduce the first-derivative term of the ODE.Here we show that it is a general property for any stratification, not just an isothermal one.

Local Dispersion Relation
For a general stratification, the coefficients of the ODE ( 23) are functions of height and the solutions will not be sinusoidal.However, by making a change of variable that converts the ODE into standard form (i.e., a Helmholtz equation that lacks a first-derivative term), a local dispersion relation can be generated which is appropriate in a WKB framework (e.g., Bender & Orszag 1999).The required change of variable involves the square root of the mass density, δP = (αρ 0 ) 1/2 ψ.We include the constant α inside the square root purely for the sake of symmetry in later sections when we explore various sound-proofing techniques.Here, its inclusion is unnecessary and only introduces a multiplicative constant which factors out of the resulting ODE, In the preceding equations, k z (z) is a local vertical wavenumber and ω c (z) is the acoustic-cutoff frequency which depends on the stratification through the density scale height H, We denote vertical derivatives of atmospheric quantities using a superscript prime, i.e., the vertical derivative of the density scale height is given by H ′ ≡ dH/dz.
From the preceding analysis, we see that acousticgravity waves vary over two relevant vertical spatial scales: a local vertical wavelength and an envelope scale.The wavelength is given by the local dispersion relation (32) and hence depends on the wave frequency as well as the characteristic frequencies of the atmosphere-i.e., the buoyancy frequency N , the acoustic cut-off frequency ω c , and the Lamb frequency k h c.The envelope scale is associated with vertical variation of the envelope function (αρ 0 ) 1/2 that appears in the change of variable above.This function provides a local amplitude of the wave function (in a WKB sense).
Since the envelope function only depends on the mass density, the envelope scale is solely determined by the atmospheric stratification through the density scale height H.For later convenience, we choose to define the envelope scale Λ as twice the scale length associated with the envelope function such that Λ = H,

INTERNAL GRAVITY WAVES IN THE LOW-FREQUENCY LIMIT
Our primary goal is to see how each wave variable scales with frequency and to therefore determine which terms are important in the fluid equations in the lowfrequency limit.We start by non-dimensionalizing, using the reciprocal of the horizontal wavenumber k −1 h and the frequency of surface gravity waves √ gk h for the characteristic length and frequency.We choose c p and ρ to be typical values of the entropy and mass density, respectively.Of particular importance is the non-dimensional wave frequency, which will serve as a small parameter in our lowfrequency expansions.Thus, when we speak of low frequencies we are considering frequencies that are small compared to those of surface gravity waves, ω 2 ≪ gk h or equivalently ε ≪ 1.This assumption will assure that the acoustic waves and the internal gravity waves decouple cleanly.
In combination, Equations ( 31) and ( 32) indicate that the vertical wavelength of an internal gravity wave becomes very short as the frequency vanishes.To leading order in the frequency, the vertical wavenumber is determined by the ratio of the buoyancy frequency to the wave frequency, Hence, in the low-frequency limit the vertical wavelength becomes a short spatial scale, whereas the envelope or atmospheric scale remains long.This scale separation dictates that we must define a non-dimensional height ζ that appropriately rescales the vertical derivatives in the fluid equations to respect the short scale, If we denote the non-dimensional forms of the wave variables and atmospheric profiles using a tilde, the wave equation ( 23) becomes, where the non-dimensional atmospheric profiles are given by and the non-dimensional Lagrangian pressure fluctuation is defined as follows: Similarly, the non-dimensional form for the local dispersion relation is given by where kc is a nondimensional wavenumber that is the ratio of the acoustic cutoff frequency to the Lamb frequency, .
As expected, the leading order behavior of the local vertical wavenumber in Equation ( 41) demonstrates that the vertical wavelength becomes very short in the lowfrequency limit, Modifications to the vertical wavenumber arising from a finite frequency first appear at order unity, O(1), whereas the term in the dispersion relation responsible for the propagation of highfrequency acoustic waves appears at O(ε 2 ).

Frequency Dependence of the Other Wave Variables
The non-dimensional forms of the other fluid variables can be generated through Equations ( 18), ( 19), and ( 26) and are related to the Lagrangian pressure fluctuation through differential operators, where ρ0 = ρ 0 /ρ is the non-dimensional atmospheric density.
We can immediately see that internal gravity waves possess motions that are nearly horizontal for low frequencies.The vertical velocity component w is small by a factor of ε.Furthermore, while the Lagrangian pressure fluctuation remains order unity in size, δP ∼ O(1), the Eulerian pressure fluctuation becomes small, P 1 ∼ O(ε).Both the entropy and density fluctuations remain order unity.The fact that the Eulerian pressure fluctuation vanishes in the limit of low frequency is consistent with the pseudo-incompressible approximation and ensures that the internal gravity waves and acoustic waves decouple in that limit.However, since the mass density fluctuation does not vanish, these limits further suggest that this decoupling is not accomplished through the anelastic limit.We explore this result fully in the next subsection.
In order to make obvious the relative magnitude of terms in subsequent equations, we define alternate dimensionless variables for the vertical velocity and Eulerian pressure fluctuation, Both W and Θ are order unity because the prefactors in their definitions absorb the leading-order behavior as the frequency becomes small.

Low-Frequency Limit of the Continuity Equation
Consider the dimensional form of the continuity equation (15), where the equation of state ( 16) is used to replace the density fluctuation In order to sound proof the equation set, we need to eliminate the term involving the Eulerian pressure fluc-tuation.This term is responsible for producing the pressure fluctuations that generate the restoring force for acoustic oscillations.
The anelastic approximation does indeed eliminate this pressure term, but it is overkill and removes the entire left-hand side of the continuity equation above.
In particular, the term involving the entropy fluctuation is also thrown away.For low-frequency internal gravity waves, this is inconsistent.If the continuity equation is non-dimensionalized it becomes obvious that the entropy term is the same order as other terms that are retained by the anelastic approximation, The leading-order behavior consists of the two orderunity terms that appear in square brackets on the righthand side of Equation (50).The first correction for nonzero frequency is comprised of the two first-order terms, O(ε); one of these is the foresaid entropy term.
The term involving the Eulerian pressure fluctuation is second order, O(ε 2 ).
The lowest-order self-consistent approximation that one could make would be to keep just the leading-order terms, resulting in an assumption of incompressibility, ∇•u ≈ 0. The next self-consistent approximation would be the retention of all zero-order and first-order terms.As we will show next, this approximation is equivalent to the pseudo-incompressible condition.
We demonstrate pseudo-incompressibility by using the energy equation ( 14) to replace the entropy fluctuation in Equation ( 49) with the vertical velocity and then combining the two first-order terms using the definition of the buoyancy frequency, The last term on the right-hand side is equal to the vertical velocity divided by the scale height for the potential density, i.e., the density scale height for an adiabatic density stratification, Hence, the terms on the right-hand side of Equation ( 51) can be cleanly combined, A self-consistent low-frequency approximation is to discard all second-order terms, leading to the pseudoincompressible approximation, ∇ • (ρ * 0 u) = 0.

Low-Frequency Limit of the Momentum Equation
When transformed into the spectral representation, the vertical component of the momentum equation ( 13) is given by and non-dimensionalization of this equation yields, It is now obvious from the preceding equation that the inertial term on the left-hand side becomes the smallest term in the low-frequency limit; it is a second-order correction.The right-hand side consists solely of terms that are zero order in the dimensionless frequency.Hence, to first order, the balance is simply the hydrostatic relation between the perturbed pressure and the perturbed density, The pseudo-incompressible and fiducial anelastic approximations both leave the vertical momentum equation unmodified.But, the LBR formulation of the anelastic approximation drops a term whose removal is formally valid only in an adiabatic (or near-adiabatic) stratification.The vertical momentum equation ( 54) can be rewritten in the following manner either by linearizing Equation (3) or by dividing the vertical momentum equation ( 54) by the mass density and pulling the density into the gradient operator that appears in the pressure force by use of the chain rule.
The LBR formulation of the anelastic approximation removes the term involving the buoyancy frequency, even in stable stratifications where the buoyancy frequency is not small.We demonstrate that this approximation is inconsistent with low-frequency gravity waves by considering the nondimensional form of Equation (57), The LBR approximation inconsistently ignores the firstorder term while retaining the inertial term (which is second-order).

THE INTEGRITY OF THREE SOUND-PROOFING TREATMENTS
In this section we examine the success or failure of a variety of sound-proofing methods in reproducing the appropriate behavior of low-frequency internal gravity waves.We have already discussed how all anelastic formulations inconsistently reject terms in the continuity equation and how the LBR anelastic formulation is further inconsistent with its treatment of the vertical momentum equation.Here we will examine how these inconsistencies propagate and produce errors in the dispersion relation and wave functions.To ease comparison, here we provide the local dispersion relation for a fully compressible fluid in both its dimensional and nondimensional forms-i.e., Equations ( 32) and (41), Further, in Table 1, we summarize the function α(z), the local wavenumber k z , and the envelope scale Λ in the low-frequency limit for a fully compressible fluid and for all three sound-proofing treatments.We retain terms only up to first-order in the dimensionless frequency ε.

Pseudo-incompressible Approximation
Since the pseudo-incompressible approximation is selfconsistent in its treatment of the continuity equation and correct to first order in the frequency, we expect that this approximation should produce low-frequency internal gravity waves that are correct to first order in the dispersion relation and in the wave functions.To demonstrate that this expectation is true we rederive the wave equation for internal gravity waves but with the continuity equation ( 15) replaced by the pseudoincompressible condition, ∇ • (ρ * 0 u) = 0. We simply present the result, Table 1.Wave properties achieved under various sound-proofing approximations as indicated in the first column.The second column indicates the function α(z).The third and fourth columns provide the square of the local vertical wavenumber k 2 z , and the scale length Λ of the amplitude envelope for internal gravity waves in the low-frequency limit.The wave frequency and horizontal wavenumber are indicated by ω and k h , respectively.The atmosphere is characterized by the vertical profiles of the sound speed c, the density scale height H, the scale height for an adiabatic stratification (i.e., the scale height for the potential density) H * = c 2 /g, the buoyancy frequency N , and the acoustic cutoff frequency ωc.For the vertical wavenumber and envelope scale, all terms with a magnitude O(ε 2 ) or smaller have been neglected.Since, the leading-order terms in the vertical wavenumber are O(ε −2 ), the neglected terms are small by a factor of ε 4 , i.e., they are fourth order.
In this expression, θ PI (z) is a dimensionless function that depends on the temporal frequency ω, the horizontal wavenumber k h , and the potential density ρ * 0 through the following definitions Compared to the fully-compressible equations, the quantity α PI has been augmented by ω 2 g/H * , and is therefore no longer a constant function of height.
A direct comparison of Equation ( 61) with the wave equation for a fully compressible fluid (23) reveals that there are three spurious terms: both of the terms involving θ PI , as well as the term (N 2 /c 2 )δP .To demonstrate that all of these spurious terms are small in magnitude and can be safely ignored in the low-frequency limit, we nondimensionalize Equation (61), and we recognize that the function θ PI is an order-unity quantity for low frequencies, Thus, all of the spurious terms are second-order or higher in the dimensionless frequency ε and the Lagrangian pressure fluctuation that is generated by Equation ( 61) is correct to first-order.
Based on this result we should expect the local dispersion relation to also be correct to first order and this is indeed the case.The transformation that converts the ODE into a Helmholtz equation has the same form as we found for the fully-compressible equations, but now the function α = α PI (z) varies with height.This change of variable leads to the following local dispersion relation, with a nondimensional form given by All of the terms contained by the error term, E PI , in the preceding equations are spurious and do not appear in the local dispersion relation for a fully compressible fluid.However, all spurious terms appear as a correction that is smaller than the leading order behavior by a factor of ε 2 or smaller.Hence, the pseudo-incompressible approximation leads to a local dispersion relation that is correct to first order.Finally, the envelope scale can be read directly from the coefficient of the first-derivative term in the ODE, Λ −1 = H −1 + ω 2 θ PI /g.To first order in the frequency, the envelope scale is simply the density scale height.

Fiducial Anelastic
For the fiducial anelastic approximation, where the only modification to the fully-compressible fluid equations is made to the continuity equation, the resulting ODE for the Lagrangian pressure fluctuation is as follows, where the α and θ functions take on subtly but crucially different forms, Here, α FA and θ FA differ from the pseudo-incompressible case, Equations ( 62) and ( 63), by the appearance of H instead of H * .
A direct comparison of Equation ( 69) with the ODE (23) appropriate for a fully compressible fluid illustrates that fiducial anelastic generates a variety of spurious and incorrect terms.Specifically, the terms in the square brackets are spurious and the entire coefficient of the first-derivative term is incorrect.To ascertain the magnitude of these mistakes, we nondimensionalize, Fiducial anelastic performs rather poorly in reproducing the behavior of low-frequency internal gravity waves.
The ODE is correct only to leading order in ε with inconsistencies appearing at first-order in the coefficient of the first derivative.The first term in this coefficient contains the reciprocal of the scale height of the potential density, where it should instead possess the reciprocal of the density scale height-see Equation ( 38).Interestingly, conversion of the ODE to standard form-via the change of variable δP = (α FA ρ * ) 1/2 ψresults in a local dispersion relation that is correct to first order, or (74) In addition to all of the spurious terms that appear in the square brackets, the acoustic cut-off frequency is incorrect, For ease of comparison, in Table 1 we have reworked the right-hand side of Equation ( 73) to extract the correct form of the acoustic cutoff frequency.Despite these issues, the errors all appear at second order or higher in the dimensionless frequency ε, meaning that the erroneous terms divided by the leading order behavior are small by a factor of ε 2 .The fact that the ODE itself is incorrect at first order manifests in the envelope function, (α FA ρ * ) 1/2 , which is wrong at all orders.As we will see in a subsequent section this results in first-order errors in the wave functions even though the dispersion relation is correct to first order.

LBR Anelastic
In the framework of the LBR anelastic approximation, in addition to the anelastic treatment of the continuity equation, i.e., ∇ • (ρ 0 u) ≈ 0, a term in the vertical momentum equation is removed.When these two modification to the fluid equations are adopted, the resulting ODE that describes internal gravity waves becomes, where α and θ are now The non-dimensional form of the ODE becomes Despite the inconsistent treatment of the vertical momentum equation, the LBR form of the anelastic approximation generates an ODE that is correct to first order in ε.The spurious terms that appear in the square brackets are second order or higher and the coefficient of the first derivative is correct to first order.As expected the local dispersion relation-once again achieved by the change of variable δP = (α LBR ρ 0 ) 1/2 ψ, is correct to first order, and (81)

Comparison of the Vertical Wavelengths
In the previous subsections we demonstrated that the three approximations generate errors to the vertical wavelength of internal gravity waves that are second order in the dimensionless frequency ε.Hence, if the only test of fidelity was to reproduce the local dispersion relation, all of the sound-proofing treatments would fair equally well.This is born out by a comparison of the vertical wavenumber that is achieved in an isothermal atmosphere by each treatment.This type of atmosphere is one of the most lenient of all potential atmospheres because all of the characteristic frequencies, i.e., N , ω c , and k h c, become constant functions of height, as do the scale heights H and H * .As a consequence, the vertical wavenumber k z becomes a constant and the quantity θ vanishes identically for all approximations.When θ is zero, many of the spurious terms disappear from the local dispersion relations.
Figure 1 shows the performance of the three approximations in an isothermal atmosphere.The leftmost panel illustrates the isocontours of the vertical wavenumber achieved in a fully-compressible fluid as a function of horizontal wavenumber k h and temporal frequency ω.The remaining three panels provide the same isocontours for the sound-proofing treatment indicated at the top of the panel.The solid black contours in each panel are for the fully-compressible fluid, while the dashed red curves show the same contours under the relevant approximation.The value of each contour is marked in the left-most panel.In each panel, four isocontours of the nondimensional frequency ε = ω/ √ gk h are overlayed for reference and appear as dotted orange curves.To see how well an approximation reproduces the correct behavior, one should compare the red and black curves within a panel.We would expect that the differences should be small for low values of the nondimensional frequency, i.e., in the lower-right portion of the diagram, and large for high values of ε (upper left).From the four panels, it is clear that all three approximations reproduce the vertical wavenumber well as long as the dimensionless frequency is small, i.e., ε ≲ 0.3.

Comparison of Wave Cavity Boundaries
Since an isothermal atmosphere is so special (because many of the spurious terms in the dispersion relation vanish), it is wise to examine the behavior of the local dispersion relation in a more complicated atmosphere.We have chosen to examine the vertical wavenumber in a polytropically stratified atmosphere.Such atmospheres have thermodynamic profiles that are power laws in the depth, In the expressions above, A is an arbitrary constant and m is the polytropic index.Polytropes can be stably or unstably stratified depending on the values of the adiabatic index γ and the polytropic index m; if m > (γ − 1) −1 , the atmosphere is stable to convective overturning.
A convenient feature of a polytropic atmosphere is that it is self-similar, lacking an intrinsic spatial scale (see Hindman & Jain 2022).Therefore, the local dispersion relation becomes independent of the horizontal wavenumber if we express the frequency in terms of our nondimensional frequency ε = ω/ √ gk h and we write all of the atmospheric profiles using a nondimensional depth −k h z.Because of this property, we can generate a single dispersion diagram that illustrates the vertical wavenumber as a function of dimensionless depth and frequency that is valid for all horizontal wavenumbers.Figure 2 presents the resulting dispersion diagram for each treatment of the fluid equations for a polytropic atmosphere with a polytropic index of m = 3 (which is stably stratified for an adiabatic index of γ = 5/3).The left-most panel is for a fully-compressible fluid and the right three panels are for the three sound-proofing formalisms.The blue region in each diagram corresponds to those depths in the atmosphere where a wave of the given frequency is vertically evanescent and the black and red contours have the same meaning as in Figure 1.The upper panels show a range of dimensionless frequency that is wide enough to contain both the branch of low-frequency internal gravity waves and the branch of high-frequency acoustic waves (if present).The lower panels show a zoom-in view at low-frequencies that focuses on the gravity waves.Note, at a given frequency, there are two turning points where the local vertical wavenumber vanishes.Hence, the internal gravity waves are vertically trapped in a wave cavity for g modes.The turning points are indicated by the thick curves.Similarly, in the fully-compressible fluid, the acoustic waves are trapped in a p mode cavity.
The pseudo-incompressible condition does minimal damage to the g-mode cavity, see Figure 2b.The boundaries move only slightly even for the highest-frequency waves.Further, the vertical wavenumbers within the cavity are weakly affected even for frequencies that  1 for a summary).In each panel, the solid black curves correspond to the isocontours of the square of the dimensionless vertical wavenumber (kzH) 2 for a fully compressible atmosphere (where the density scale height H is a constant function of height for an isothermal atmosphere).The value of each contour is indicated by a black label in panel a.Further, the thick black contour corresponds to the zero contour that separates domains of vertical wave propagation (k 2 z > 0) and evanescence (k 2 z < 0).In panels b-d, the dashed red curves indicate the same contours but for the approximation indicated at the top of the panel.In each panel, the domain of evanescent waves is indicated by the blue shading, while the region of vertical propagation is unshaded.The dotted curves in each panel are isocontours of the dimensionless frequency.Since the dimensionless frequency is a function of wavenumber, ε = ω/ √ gk h , isocontours are curved lines with low values in the lower-right portion of the diagram and high values in the upper left.All approximations reproduce the correct vertical wavenumber when the dimensionless frequency ε is small.Differences between the approximations begin to appear for moderate to large values of the dimensionless frequency ε > 0.3.
are large enough that we might suspect that the lowfrequency limit is invalid.The anelastic models fare poorly, however.The fiducial anelastic approximation does a horrendous job of reproducing the wave cavity boundaries.In fact, there appears to be a residual of the acoustic cavity that is highly distorted and appears at frequencies halfway between the acoustic and gravity wave branches.While, the LBR form of the anelastic approximation does not have spurious wave cavities at high frequency, it fails to reproduce the boundaries of the g-mode cavity with fidelity.The highest-frequency gravity waves that are vertically propagating have frequencies that are too high by a factor of about one-third.Further, errors in the vertical wavenumber become noticably large for relatively low values of the dimensionless frequency, ε > 0.1.

Errors in the Wave Functions
In sections 5.1-5.3,we found that the pseudoincompressible approximation and the LBR formulation of the anelastic approximation both introduced errors in the Lagrangian pressure fluctuation that appeared at second order.The fiducial anelastic approximation produced errors at first-order.So at first glance the LBR approximation seems to fare well.However, as we shall soon see, when we consider other wave vari-ables, such as the fluid velocity components, the pseudoincompressible approximation becomes the clear winner.
In the same manner that one derives equations ( 18) and ( 19) for a fully compressible fluid, similar equations can be derived for each of the approximations.When pseudo-incompressibility is adopted, one obtains the following: To see the magnitude of the spurious terms, we nondimensionalize, If one compares these expressions with Equations ( 43) and ( 44), it is clear that all spurious terms appear at second order in the dimensionless frequency.Since all of the other fluid variables (i.e., ρ 1 , P 1 , and s 1 ) are linear combinations of the two velocity components-see The pseudo-incompressible and LBR anelastic approximations eliminate all such acoustic waves.The fiducial anelastic approximation leaves a highly distorted residual domain of propagating acoustic waves.In general, all three approximations do well in reproducing the correct vertical wavenumber when the dimensionless frequency is small ε ≲ 0.1.However, the pseudo-anelastic approximation has the least distortion to the spatial extent of the wave cavity even for frequencies as large as ε ≈ 0.3.26), the wave functions for all of the fluid variables are correct to first order when the pseudoincompressible approximation is utilized.

Equation (
Both of the anelastic approximations falter.For the fiducial anelastic approximation the nondimensional forms for the two velocity components are and for the LBR anelastic formulation we obtain with αLBR ≡ 1 − ε 4 + ε 2 Ñ 2 + ε 2 H−1 .Both have errors in the horizontal velocity that appear at first order (i.e., the term involving ε Ñ 2 ).The fiducial anelastic approximation has the added shame that the Lagrangian pressure fluctuation itself is only correct to zero order and hence all fluid variables suffer from the same deficiency.For the LBR approximation, the first order error in the horizontal velocity u propagates to errors of similar size in the fluctuations of the Eulerian pressure P 1 and density ρ 1 .

DISCUSSION
We have demonstrated that internal gravity waves within a fully-compressible fluid become pseudoincompressible in the low-frequency limit.Discrepancies from the solutions for a fully compressible fluid appear at second order in the non-dimensional frequency, i.e., the relative errors are O(ω 2 /gk h ).Conversely, the two anelastic approximations that we consider are inconsistent in the terms they neglect or retain in the continuity equation and vertical momentum equation.This inconsistency leads to errors in the wave functions that appear at first order, O(ω/ √ gk h ).A summary of the fractional errors in the vertical wavenumber, envelope scale length, and in the eigenfunctions appears in Table 2.
These errors in the eigenfunctions arise from errors in either the local vertical wavenumber (the short spatial scale) or in the amplitude envelope of the oscillations (the long spatial scale)-see Tables 1 and 2. Many of the errors in the local dispersion relation explicitly require vertical variation in the atmospheric profiles of the density scale height and buoyancy frequency.Both Brown et al. (2012) and Vasil et al. (2013) explicitly considered isothermal atmospheres for which the scale heights and the characteristic frequencies are constants.So many of the errors identified here failed to materialize in those previous studies.Brown et al. (2012) examined the behavior of internal gravity waves under the influence of three distinct anelastic treatments (including the LBR and fiducial anelastic formulations), and found that the LBR formulation suffered from the least deviation from the fully compressible result.Here we have demonstrated that the apparent success of the LBR approximation is only in reproducing the local dispersion relation.If one considers the wave functions directly, the LBR anelastic approximation fails at first order, just like fiducial anelastic.

Conservation of Energy
We can explore conservation of energy under each approximation by computing the vertical energy flux F (z).Using Abel's Identity, as we did for a fully-compressible fluid in section 3.1, we find a general expression for the energy flux that is valid for all three sound-proofing treatments, Each approximation generates a distinct form for α and has a different Wronskian because the coefficients of the first-derivative term in the respective ODEs differ.
For the pseudo-incompressible equations, using Equations ( 62) and ( 63), we find that the vertical energy flux is a constant function of height, Hence, energy is conserved.It is interesting to note that we have not utilized the small parameter in this derivation of the energy flux.So, energy is conserved even when the low-frequency expansions have questionable validity because the dimensionless frequency is not small.Performing the same calculations for the two anelastic treatments reveals that the LBR formulation conserves energy (for the same reasons that the pseudoincompressible equations do) and the fiducial anelastic equations lack energy conservation, The vertical energy flux F FA derived from the fiducial anelastic equations depends on the atmosphere's specific entropy density and, thus, in an atmosphere without adiabatic stratification the wave will deposit or extract energy as it travels.

Applicability in Numerical Simulations
In numerical simulations, it is hard to overstate the utility in converting the continuity equation from a parabolic prognostic equation to an elliptic PDE constraint, as is accomplished by both the anelastic and pseudo-incompressible approximations, In addition to removing sound waves and hence unthrottling the simulation's timestep, the imposition of constraints with this form allow the fluid velocity to be expressed using stream functions.Of course, this reduces the number of variables that must be evolved from one time step to the next.However, this is done at the expense of increasing the spatial order of the now reformulated momentum equations in stream function form that is now devoid of any elliptic constraints.This may demand auxiliary boundary conditions on the streamfunctions that are not readily available.Moreover, if linear coupling in the system is treated as explicit in numerical time-stepping algorithms, it is known, specifically for spectral schemes, that the numerical accuracy of the scheme can be degraded at high resolutions.Fortunately, recent advances have shown that this degradation is avoided if linear couplings remain implicit at Table 2. Magnitude of the fractional errors that are introduced in internal gravity waves by three different sound-proofing techniques.Each column lists the size of the error divided by the leading order behavior for the wave property indicated at the top of the column.The size of each error is presented in terms of the dimensionless frequency ε = ω/ √ gk h .The pseudoincompressible approximation evinces the smallest errors, all appearing at second order.Both of the anelastic approximations have errors that appear at first order or larger.
In the derivation of the pseudo-incompressible condition above, two related assumptions are made.First, the Mach number, Ma, of the flows is small such that the advection timescale is much longer than a soundcrossing time for a typical flow structure.This ensures that fluid motions are in a constant state of pressure equilibration-i.e., the Eulerian pressure fluctuation is small.Second, we have assumed that fluctuations in the potential density are small compared to that of the background state.This later assumption is self-consistent with low-Mach number flows.Notably, unlike the anelastic approximation discussed below, it does not restrict density fluctuations to be small compared to that of the background state.Finally, since we have ignored diffusive effects in the derivation of the pseudo-incompressible constraint, i.e., we have ignored Q in Equation ( 8), we have made the further assumption that the Péclet number is large, Pe ≫ 1, such that the thermal diffusion timescale is long compared to the advective time scale.To summarize, for the pseudoincompressible constraint to be valid, we must have the following ordering of timescales, or equivalently in terms of nondimensional numbers where U is a typical flow speed, L is a typical length scale, and κ is the thermal diffusivity.The validity of the anelastic constraint requires the same assumption of low Mach number, Ma ≪ 1, but makes a different stricture on the effectiveness of thermal diffusion.Since, we must ignore Eulerian fluctuations of the mass density in the continuity equation, the equation of state dictates that, in addition to small pressure fluctuations, we must have small entropy or temperature fluctuations.In the convection zone of a star or planet, where the stratification is essentially adiabatic, entropy fluctuations are naturally small; anelasticity holds; and the anelastic and pseudo-incompressible conditions are equivalent.However, in a region of stable stratification, the only way that the entropy or temperature fluctuations can remain small is if temperature homogeneity is diffusively maintained across flow structures (see Bannon 1996).This requires that the thermal diffusion time is short compared to the advective time scale.Summarizing, anelasticity requires or equivalently The limitation of low Mach number is easily met in many astrophysical and geophysical applications.Convection is sedate in the Jovian planets, in the Earth's interior, and in the deep layers of low-mass stars.Wave motions and circulations in the stably stratified regions of stars and planets are similarly often low Mach number.The requirements on the Péclet number are usually the more restrictive of the two assumptions.For example, the thermal diffusion time in the Sun is typically millions of years; using the solar radius as the length scale, L = R ⊙ ≈ 700 Mm, and a thermal diffusivity appropriate for photon diffusion, κ ∼ 10 7 cm 2 s −1 , we obtain τ diff ∼ 16 Myr.If we consider the meridional circulation at the base of the Sun's convection zone and adopt a typical flow speed of 1 m s −1 , we obtain an advective timescale of 20 years, leading to a Péclet number of Pe ∼ 10 6 .Clearly, these motions are not anelastic; thermal diffusion cannot act rapidly enough to eliminate the temperature fluctions generated by advection.However, the motions do satisfy both of the requirements for pseudo-incompressibility, Ma ≪ 1 and Pe ≫ 1.Although large Péclet numbers are often true from an astrophysical perspective, numerical simulations are often performed in regimes where Pe ∼ O(1).The anelastic approximation offers no resolution to this problem, but the pseudo-incompressible equations do.The restriction on the Péclet number Pe can be relaxed if the irreversible thermodynamic terms are retained, where ρ * 0 is the potential density of the background state.Of course, the retention of Q will usually render a stream function formalism without the requirement of an elliptic constraint impossible.Finally, we wish to note a final advantage of the pseudo-incompressible approximation over anelasticity.While both sound-proofing schemes are well justified in a convection zone where the stratification is nearly adiabatic, if one wishes to simulate both stable and unstable regions in the same computational domain, the pseudo-incompressible approximation allows one to do so smoothly with a uniform treatment.The anelastic approximation will result in flows that violate the underlying assumptions of the approximation.

Figure 1 .
Figure 1.Propagation diagrams for an isothermal atmosphere for four treatments of the fluid equations: (a) a fully compressible fluidi.e., no approximation, (b) the pseudo-incompressible condition, (c) the fiducial anelastic approximation, and (d) the LBR formulation of the anelastic approximation (see Table1for a summary).In each panel, the solid black curves correspond to the isocontours of the square of the dimensionless vertical wavenumber (kzH) 2 for a fully compressible atmosphere (where the density scale height H is a constant function of height for an isothermal atmosphere).The value of each contour is indicated by a black label in panel a.Further, the thick black contour corresponds to the zero contour that separates domains of vertical wave propagation (k 2 z > 0) and evanescence (k 2 z < 0).In panels b-d, the dashed red curves indicate the same contours but for the approximation indicated at the top of the panel.In each panel, the domain of evanescent waves is indicated by the blue shading, while the region of vertical propagation is unshaded.The dotted curves in each panel are isocontours of the dimensionless frequency.Since the dimensionless frequency is a function of wavenumber, ε = ω/ √ gk h , isocontours are curved lines with low values in the lower-right portion of the diagram and high values in the upper left.All approximations reproduce the correct vertical wavenumber when the dimensionless frequency ε is small.Differences between the approximations begin to appear for moderate to large values of the dimensionless frequency ε > 0.3.

Figure 2 .
Figure2.Propagation diagrams for a polytropic atmosphere under different approximations to the fluid equations.In each panel, the solid black curves correspond to the isocontours of the square of the dimensionless vertical wavenumber (kz/k h ) 2 for a fully compressible atmosphere.These contours are plotted versus a non-dimensional depth, −k h z, and the dimensionless frequency, ε = ω/ √ gk h .The thick black contour corresponds to the zero contour that separates domains of vertical propagation (k 2 z > 0) and evanescence (k 2 z < 0).The dashed red curves indicate the same contours but for the approximation indicated at the top of the column.The background colors have the same meaning as in Figure1.The upper panels illustrate a larger range of frequency and capture the high-frequency acoustic branch.The pseudo-incompressible and LBR anelastic approximations eliminate all such acoustic waves.The fiducial anelastic approximation leaves a highly distorted residual domain of propagating acoustic waves.In general, all three approximations do well in reproducing the correct vertical wavenumber when the dimensionless frequency is small ε ≲ 0.1.However, the pseudo-anelastic approximation has the least distortion to the spatial extent of the wave cavity even for frequencies as large as ε ≈ 0.3.

Table 1
Comparison of Sound-Proofing Techniques

Table 2
Fractional Errors Introduced by Sound-Proofing Techniques Errors in the Errors in the Errors in the Equation Set Vertical Wavenumber kz Envelope Scale Λ Wave Functions δP , u