2.1 The governing equations

In order to determine the true convective stability threshold, we linearise the fully compressible equations about a hydrostatic background in a reference frame rotating with angular velocity $\unicode[STIX]{x1D734}$ . For simplicity we neglect self-gravity and the centrifugal force, so the gravitational acceleration $\boldsymbol{g}$ points directly downward. We adopt Cartesian coordinates in which $\boldsymbol{g}=-g\boldsymbol{e}_{z}$ and $\unicode[STIX]{x1D734}=\unicode[STIX]{x1D6FA}\boldsymbol{e}_{z}$ , and write the linearised equations in the form

where subscripts 0 and 1 refer to unperturbed quantities (which are functions of $z$ only) and their linear perturbations, respectively. The quantities $\boldsymbol{f}_{1}$ and $q_{1}$ are the viscous force and heating rate per unit mass, whose precise form we specify later. In order to close these equations we require relations between the thermodynamic perturbations. For convenience we introduce the specific enthalpy, $H(p,s)$ , which is a function of pressure, $p$ , and specific entropy, $s$ (e.g. Landau & Lifshitz Reference Landau and Lifshitz1980, §14). We can then define the density, $\unicode[STIX]{x1D70C}=1/H_{p}$ , and temperature, $T=H_{s}$ , where the subscripts on $H$ represent partial derivatives. The linear perturbations to $\unicode[STIX]{x1D70C}$ and $T$ are then given by

with similar expressions applying to the spatial variations of the background state:

(Note that $H_{ps}=H_{sp}$ .) The second derivatives in (2.4)–(2.7) can also be expressed in terms of the sound speed, $c$ , and specific heats, $c_{p}$ and $c_{v}$ , as follows

Finally, the buoyancy frequency, $N$ , is defined by the formula $N^{2}=-H_{ps}(\text{d}p_{0}/\text{d}z)(\text{d}s_{0}/\text{d}z)$ . We are concerned here with convective instability, which requires that $N^{2}<0$ . For notational convenience, we therefore define the parameter $\mathscr{N}=-N^{2}$ as a measure of the degree of superadiabaticity. In what follows we will assume that $\mathscr{N}$ is positive throughout the domain, but we will make no assumption about its magnitude. Note that the right-hand side of the momentum equation (2.1) can be written in many different forms, by using the relations (2.4)–(2.7) and the hydrostatic condition $\unicode[STIX]{x1D735}p_{0}=\unicode[STIX]{x1D70C}_{0}\boldsymbol{g}$ ; the form used here has been chosen for later convenience.

The linearised Boussinesq equations, for comparison, are

where $\unicode[STIX]{x1D70C}_{0}$ is approximated as a constant but $s_{0}$ retains its dependence on $z$ . In the Boussinesq approximation, the viscous force and the heating term take the form

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity and $\unicode[STIX]{x1D705}$ is the thermal diffusivity. The derivation of these Boussinesq equations relies upon the following assumptions (e.g. Spiegel & Veronis Reference Spiegel and Veronis1960; Mihaljan Reference Mihaljan1962):

1. (i) the domain height, $L$ say, is much smaller than the scale height of each thermodynamic variable;

2. (ii) perturbations to the thermodynamic variables are even smaller than the variations in their background values;

3. (iii) fluid motions are slow in comparison with the sound speed, $c$ ;

4. (iv) any time scale of the flow is much longer than the acoustic time scale, $L/c$ .

Often these assumptions are not stated explicitly, but are implicit in the scalings assumed for the various physical quantities. Many derivations incorporate additional assumptions in order to simplify the Boussinesq equations still further, for example by neglecting viscous heating (Mihaljan Reference Mihaljan1962; Veronis Reference Veronis1962; Gray & Giorgini Reference Gray and Giorgini1976). In what follows we examine the onset of oscillatory convective instability in the Boussinesq and compressible systems, which renders such additional assumptions unnecessary, because viscous heating vanishes for linear perturbations to a hydrostatic state. Moreover, conditions (ii) and (iii) are automatically satisfied in this case, because near onset the perturbations to the hydrostatic background state are infinitesimally small. Our aim is therefore to determine whether the Boussinesq and fully compressible systems have the same convective stability threshold in the regime described by conditions (i) and (iv).

2.2 The linear stability of the Boussinesq system

We first derive the stability threshold for the Boussinesq equations, summarising the results of Chandrasekhar (Reference Chandrasekhar1953). Suppose, for simplicity, that $\unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D705}$ and $\mathscr{N}$ are positive and constant throughout the domain. We can then seek solutions of (2.9)–(2.13) in which the perturbations are $\propto \text{exp}(\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}-\text{i}\unicode[STIX]{x1D714}t)$ , and thus obtain a dispersion relation

where $k=|\boldsymbol{k}|$ and $k_{h}=\sqrt{k_{x}^{2}+k_{y}^{2}}$ . The onset of instability can be either direct (i.e.  $\unicode[STIX]{x1D714}=0$ ) or oscillatory (i.e.  $\unicode[STIX]{x1D714}^{2}>0$ ) depending on parameter values. In the former case, a perturbation with given $\boldsymbol{k}$ is convectively unstable if

In the latter case, such a perturbation is unstable if

and, at onset, it oscillates with a frequency given by

In the simplest case of Rayleigh–Bénard convection between horizontal plates separated by distance $L$ , with fixed-temperature and stress-free boundary conditions, the most unstable mode can be found by simply setting $k_{z}=\unicode[STIX]{x03C0}/L$ and finding the minimum unstable value of $\mathscr{N}$ over all $k_{h}$ . From an examination of (2.16) and (2.17) it can be seen that oscillatory instability is favoured in a rapidly rotating fluid with low viscosity, and that in such cases the instability broadly resembles a growing inertial wave. In the double asymptotic limit with $\unicode[STIX]{x1D708}\ll \unicode[STIX]{x1D705}\ll \unicode[STIX]{x1D6FA}L^{2}$ the unstable modes are quasi-geostrophic ( $k_{h}^{2}\gg k_{z}^{2}$ ), and the instability boundary is approximately described by

and

Minimising over $k_{h}$ , we obtain the asymptotic scalings for the instability at onset:

It is sometimes helpful to rewrite these results in terms of dimensionless quantities. For Boussinesq convection, the standard dimensionless numbers are the Rayleigh, Taylor and Prandtl numbers, defined respectively as

After putting $k_{z}=\unicode[STIX]{x03C0}/L$ , equation (2.23) can be rewritten as

(Chandrasekhar Reference Chandrasekhar1953) and the double asymptotic limit in which this result applies can be expressed as $1\ll Pr^{-1}\ll Ta^{1/2}$ . Equations (2.21) and (2.22) imply that $\unicode[STIX]{x1D714}\simeq \sqrt{2}\unicode[STIX]{x1D705}k_{h}^{2}$ , so the horizontal length scale of the instability is such that heat can diffuse between the upflows and downflows on the same time scale as their oscillation period. The instability can be understood physically as an inertial wave that is damped by viscosity but amplified by thermal diffusion, which feeds heat into the upflows, making them buoyant, thereby extracting potential energy from the background stratification.

2.3 The linear stability of the fully compressible system

We now consider the stability properties of the fully compressible equations (2.1)–(2.5). We recall that these equations are valid for an arbitrary equation of state and an arbitrary hydrostatic background, so in general their solutions will have complex, non-sinusoidal dependence on $z$ . Nevertheless we can still seek solutions with a given horizontal wavenumber $k_{h}$ and frequency $\unicode[STIX]{x1D714}$ . Our aim is to determine in what asymptotic limit (if any) the fully compressible equations have the same instability threshold as the Boussinesq equations. Assuming that such a limit does exist, we anticipate that the unstable modes will follow the Boussinesq scaling laws (2.21)–(2.23), and in particular that $k_{h}\gg 1/L$ . In that case, we can neglect vertical diffusion, and approximate the diffusive terms in (2.1) and (2.3) as

(Note that we use the correct definition of $\unicode[STIX]{x1D705}$ given by Spiegel & Veronis (Reference Spiegel and Veronis1960), rather than that of Chandrasekhar (Reference Chandrasekhar1961).) The validity of this approximation will be confirmed later, by solving the full set of linear equations numerically. With this simplification, the linear equations (2.1)–(2.5) can be reduced to a pair of coupled ordinary differential equations for the perturbation quantities $p_{1}$ and $u_{z}$ :

In order to make contact with the Boussinesq results, we must now impose the restrictions built into the Boussinesq approximation. First, we take the limit in which the height of the domain, $L$ , is much smaller than any length scale on which the background state varies. This corresponds to condition (i) given earlier. In this limit, and on the assumption that the scaling $\unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D705}k_{h}^{2}$ is still valid, the left-hand sides of (2.27) and (2.28) are dominated by the derivative terms, and the parenthetical terms on the right-hand sides can be approximated as constant. Under condition (i), it is reasonable to expect growing modes near onset to have a sinusoidal dependence on $z$ , which means that we can replace $(\text{d}/\text{d}z)\rightarrow \text{i}k_{z}$ , thus obtaining an approximate dispersion relation

Next, we assume that the oscillation frequency of the instability is much smaller than the resonant acoustic frequency (i.e.  $\unicode[STIX]{x1D714}\ll c/L$ ), in accordance with condition (iv). Provided that the scalings (2.21)–(2.23) are still correct in order of magnitude, the real and imaginary parts of the dispersion relation (2.29) imply that

and

at the onset of instability. We note that the imaginary part of (2.29) is smaller than the real part by a factor of $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}\ll 1$ . Therefore in order to obtain (2.30)–(2.31) it is not sufficient to consider only the leading-order terms in (2.29).

Comparing (2.19) and (2.30), we see that the Boussinesq results are valid only if $(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705})^{1/2}c\gg \unicode[STIX]{x1D6FA}/k_{h}\sim \unicode[STIX]{x1D714}L$ , which is much more restrictive than the condition $c\gg \unicode[STIX]{x1D714}L$ assumed above. Conditions (i)–(iv) are therefore not sufficient to guarantee the validity of the Boussinesq approximation. The additional term in (2.30), compared with (2.19), renders the fully compressible system subject to conditions (i)–(iv) more stable than the equivalent Boussinesq system, because each wavenumber becomes marginally stable at a larger value of $\mathscr{N}$ .

The discrepancy between the true instability threshold and that predicted by the Boussinesq equations suggests that some effect of compressibility remains significant even under conditions (i)–(iv). In the compressible system, equations (2.3) and (2.26) describe the advection of entropy and the diffusion of temperature. If we rewrite (2.3) in terms of the density and pressure perturbations, then we obtain

The term in parentheses is complex, and its imaginary part quantifies the density perturbations arising from pressure perturbations that are out of phase. This is the source of the extra term in (2.30) compared with (2.19). In the Boussinesq equations, the density, temperature and entropy perturbations are all proportional to one another ((2.12) and (2.13)), so this effect is absent. We note that this discrepancy does not arise in direct (i.e. non-oscillatory) convection, for which all perturbations are necessarily in phase with one another, so our findings apply only to oscillatory instabilities.

2.4 A specific example: the case of an ideal gas

To more precisely illustrate the discrepancy between the Boussinesq and compressible results, we now consider the particular case of an ideal gas. With this simplification, the compressible system can be defined in terms of six dimensionless parameters. It is convenient to choose three of these to be $Ra$ , $Ta$ and $Pr$ , as in the Boussinesq system. For the remaining three, we will choose the ratio of specific heats, $c_{p}/c_{v}$ , the ratio of the temperature gradient to the adiabatic temperature gradient, $(\unicode[STIX]{x1D735}/\unicode[STIX]{x1D735}_{ad})=1+(H_{ss}/H_{sp})(\text{d}s_{0}/\text{d}p_{0})$ and the ratio of the domain height to the temperature scale height, $\unicode[STIX]{x1D703}=L/h_{T}$ , where $h_{T}\equiv -(\text{d}\ln T_{0}/\text{d}z)^{-1}$ . Note that, in a domain of finite height, all of these parameters except $c_{p}/c_{v}$ will generally depend on $z$ . However, if the domain is sufficiently shallow, as we will shortly assume, then they may be approximated as constant. For simplicity we will fix $c_{p}/c_{v}=5/3$ and $\unicode[STIX]{x1D735}/\unicode[STIX]{x1D735}_{ad}=2$ , so that $\unicode[STIX]{x1D703}$ is the only additional variable in our compressible system. We anticipate that the system will become Boussinesq in the limit $\unicode[STIX]{x1D703}\rightarrow 0$ , but we need to determine exactly how small $\unicode[STIX]{x1D703}$ must be for the Boussinesq results to hold. Given our choices for the other two parameters, the scale heights of pressure and density are $h_{p}\equiv -(\text{d}\ln p_{0}/\text{d}z)^{-1}=(4/5)L/\unicode[STIX]{x1D703}$ and $h_{\unicode[STIX]{x1D70C}}\equiv -(\text{d}\ln \unicode[STIX]{x1D70C}_{0}/\text{d}z)^{-1}=4L/\unicode[STIX]{x1D703}$ , and we have

Because all of the thermodynamic scale heights are of order $L/\unicode[STIX]{x1D703}$ , condition (i) is simply $\unicode[STIX]{x1D703}\ll 1$ . Condition (iv), which requires that $\unicode[STIX]{x1D714}\ll c/L$ , turns out to be more strict. Combining (2.33) with (2.22) and (2.23), this condition becomes $\unicode[STIX]{x1D703}\ll (8Pr)^{1/2}$ . But even when this condition is satisfied, the Boussinesq results may not be valid. In fact, the marginal stability criterion (2.30) then becomes, in terms of our chosen dimensionless parameters,

Here $k_{h}$ and $k_{z}$ are measured in units of $1/L$ , and we have put $k_{z}=\unicode[STIX]{x03C0}$ , assuming stress-free, fixed-temperature boundary conditions. This expression can be solved for $Ra$ , and the critical Rayleigh number is then found by minimising over all $k_{h}$ . The result is

where

We recover the Boussinesq result, which corresponds to $X=2$ , only if $\unicode[STIX]{x1D703}\ll Pr$ . Under astrophysical conditions, this is a rather extreme restriction. In the solar interior, for example, the Prandtl number is of order $Pr\simeq 10^{-6}$ and the temperature scale height is of order $10^{5}$  km. So our result implies that the Boussinesq approximation is only valid on vertical scales smaller than 1 km! Of course, convection in the solar interior is well developed, and so the linear stability analysis presented here is not directly applicable. However, there is no obvious reason why the Boussinesq equations should be any more valid in the strongly nonlinear regime than they are in the linear one. Moreover, numerical simulations of the Sun (and, indeed, other convective systems) are never performed very far above the convective threshold, owing to computational constraints, and so our results may well be directly applicable to those simulations.

In order to confirm the validity of these results, we have used a linear eigensolver to compute the critical Rayleigh number as a function of Taylor number in an ideal gas with $\unicode[STIX]{x1D703}=0.02$ for three different values of $Pr$ , with other parameters as given above. The solver, which uses the algorithm originally developed by Gough et al. (Reference Gough, Moore, Spiegel and Weiss1976), solves the exact linearised equations for an ideal gas with constant dynamic viscosity and thermal conductivity. We consider the simplest case of stress-free, fixed-temperature boundary conditions, so that the exact stability threshold can be directly compared with the asymptotic scalings obtained from the Boussinesq (2.25) and fully compressible (2.35) equations. This comparison is shown in figure 1. In the rapidly rotating regime, with $Ta\gg 1/Pr^{2}$ , the stability threshold closely matches that predicted by (2.35) in each case. As the Prandtl number is decreased, the Boussinesq prediction becomes less accurate, and always underestimates the true stability threshold.

Figure 1.

The critical Rayleigh number for onset of oscillatory convection as a function of Taylor number in a compressible fluid (solid curves) for three different Prandtl numbers. The dashed and dotted lines show the asymptotic scalings given by (2.25) and (2.35) respectively.