2.1 Equilibrium state

In figure 1 we sketch how the equilibrium density profile of a planar, pressure-confined, static, self-gravitating, polytropic structure typically appears. We use the $x$ -direction as the direction of stratification, while it is homogeneous with infinite extent in the $y$ and $z$ directions. Therefore, all equilibrium quantities, marked with a subscript $0$ , depend only on position $x$ . The equilibrium state is assumed to be static, because it is methodologically convenient to start with this case before generalising to more realistic situations in future work, and because it is physically relevant for systems with slow accretion, typically slower than the growth of perturbations. The hydrostatic equilibrium reads

(2.1) $$\begin{eqnarray}-\unicode[STIX]{x1D735}p_{0}+\unicode[STIX]{x1D70C}_{0}\boldsymbol{g}_{0}=\mathbf{0},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{0}$ , $p_{0}$ and $\boldsymbol{g}_{0}$ are respectively the equilibrium density, pressure and gravitational acceleration. The system is self-gravitating, meaning that its potential well is shaped by its own density only, which translates into the Poisson equation

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{g}_{0}=-4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{0},\end{eqnarray}$$

where $G$ is Newton’s constant. For our choice of closure relation, we proceed as follows. In the most general case, pressure and density are related by an equation of state of the form $p=p(\unicode[STIX]{x1D70C},s)$ or $p=p(\unicode[STIX]{x1D70C},T)$ where $s$ and $T$ are respectively the specific entropy and the temperature. However, for a non-magnetised fluid to be at rest in a gravitational field, it must necessarily be barotropic, i.e. $p=p(\unicode[STIX]{x1D70C})$ . Indeed, taking the curl of the hydrostatic equilibrium (2.1), since the curl of a gradient vanishes and $\boldsymbol{g}_{0}$ is the gradient of the gravitational potential, we end up with $\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}_{0}\times \unicode[STIX]{x1D735}p_{0}=\mathbf{0}$ . This means that the density and pressure gradients are aligned everywhere, and surfaces of constant density coincide with surfaces of constant pressure. Here we consider the case of a polytropic equation of state,

(2.3) $$\begin{eqnarray}p_{0}=\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}_{0}^{\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the polytropic exponent and $\unicode[STIX]{x1D705}$ is a constant. In this work we derive our results for an arbitrary polytropic index, which makes them very general.

For a polytrope (2.3), the adiabatic speed of sound $c_{a}^{2}\equiv \unicode[STIX]{x2202}p_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}$ is given by

(2.4) $$\begin{eqnarray}c_{a}^{2}=\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70C}_{0}^{\unicode[STIX]{x1D6FE}-1}.\end{eqnarray}$$

Also, the hydrostatic equilibrium (2.1) may be rewritten as

(2.5) $$\begin{eqnarray}g_{0}=c_{a}^{2}\frac{\unicode[STIX]{x1D70C}_{0}^{\prime }}{\unicode[STIX]{x1D70C}_{0}}.\end{eqnarray}$$

Plugging this into (2.2) gives the equation satisfied by the equilibrium density $\unicode[STIX]{x1D70C}_{0}$ , namely

(2.6) $$\begin{eqnarray}(\unicode[STIX]{x1D70C}_{0}^{\unicode[STIX]{x1D6FE}-1})^{\prime \prime }+\frac{4\unicode[STIX]{x03C0}G}{\unicode[STIX]{x1D705}}\frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D70C}_{0}=0,\end{eqnarray}$$

which is known as the Lane–Emden equation. Strictly speaking, we consider here $\unicode[STIX]{x1D6FE}\neq 1$ , but in fact all the results presented in this paper are valid for the isothermal case as well, by taking the limit $\unicode[STIX]{x1D6FE}\rightarrow 1$ . Doing the analysis taking $\unicode[STIX]{x1D6FE}=1$ from the beginning gives interesting analytical results, that we will describe as part of a separate paper. As detailed in Horedt (Reference Horedt2004), the Lane–Emden equation for planar polytropes has simple analytic solutions only for the special cases $\unicode[STIX]{x1D6FE}=2/3$ , 1, 2 and $\infty$ . However, more can be said for instance on the extent of these structures: $\unicode[STIX]{x1D6FE}\leqslant 1$ polytropes extend a priori to infinity, while for $\unicode[STIX]{x1D6FE}>1$ the solution has negative $\unicode[STIX]{x1D70C}$ values, which are unphysical. For the latter, it is therefore common to define the width of the structure by the position $x_{t}$ at which the solution first vanishes. Multiplying (2.6) by $(\unicode[STIX]{x1D70C}_{0}^{\unicode[STIX]{x1D6FE}-1})^{\prime }$ , integrating twice and using $\unicode[STIX]{x1D70C}_{0}(x_{t})=0$ , for centrally flat profiles $\unicode[STIX]{x1D70C}_{0}^{\prime }(0)=0$ , one can show that it is given by

(2.7) $$\begin{eqnarray}x_{t}=\sqrt{\frac{\unicode[STIX]{x1D705}}{8\unicode[STIX]{x03C0}G}}\unicode[STIX]{x1D70C}_{c}^{(\unicode[STIX]{x1D6FE}/2)-1}B\left(\frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FE}},\frac{1}{2}\right),\end{eqnarray}$$

where $B$ is the beta function and $\unicode[STIX]{x1D70C}_{c}\equiv \unicode[STIX]{x1D70C}_{0}(x=0)$ is the mid-plane value of the density. The above length $x_{t}$ must also be taken with an important caveat: by definition at this position the density vanishes, therefore the speed of sound also reaches zero, so that at this boundary any motion is supersonic and we would expect shock waves (Goldreich & Lynden-Bell Reference Goldreich and Lynden-Bell1965). In addition, the infinite extent of $\unicode[STIX]{x1D6FE}\leqslant 1$ polytropes is not realistic. We may go around these two complications by improving the modelling a little, taking into account the existence of an external pressure $p_{ext}$ , possibly small, which is relevant both in the astrophysical and cosmological contexts and accounts for the fact that no astrophysical structure is totally isolated but embedded in its non-empty surroundings. For any $\unicode[STIX]{x1D6FE}$ , we truncate the solution of the Lane–Emden equation (2.6) at the position $x_{b}$ such that the internal pressure equals the external one, i.e.

(2.8) $$\begin{eqnarray}p_{0}(x_{b})=p_{ext}\end{eqnarray}$$

and this defines the boundary of the structure. We illustrate all this in figure 1.

2.2 Equations of motion

Perturbations, marked with a subscript 1, are described in terms of the Lagrangian displacement vector $\unicode[STIX]{x1D743}$ . In the case of a static equilibrium, the Eulerian velocity perturbation $\boldsymbol{v}_{1}$ is related to $\unicode[STIX]{x1D743}$ by

(2.9) $$\begin{eqnarray}\boldsymbol{v}_{1}=\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D743}\end{eqnarray}$$

so that here the linearised momentum conservation reads

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{t}^{2}\unicode[STIX]{x1D743}=\boldsymbol{F}(\unicode[STIX]{x1D743}),\end{eqnarray}$$

with the so-called force operator (Bernstein et al. Reference Bernstein, Frieman, Kruskal and Kulsrud1958; Keppens & Demaerel Reference Keppens and Demaerel2016)

(2.11) $$\begin{eqnarray}\boldsymbol{F}(\unicode[STIX]{x1D743})=-\unicode[STIX]{x1D735}p_{1}+\unicode[STIX]{x1D70C}_{1}\boldsymbol{g}_{0}+\unicode[STIX]{x1D70C}_{0}\boldsymbol{g}_{1},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{1}$ , $p_{1}$ and $\boldsymbol{g}_{1}$ are, respectively, the density, pressure and gravitational acceleration perturbations. In (2.11), the first term gives rise to sound waves due to the compressibility of the fluid. The second corresponds to the fact that the local density perturbation is evolving in an background potential well dictated by the equilibrium potential. Finally, the third term corresponds to the fact that perturbations generate local potential wells which affect the entire density profile $\unicode[STIX]{x1D70C}_{0}$ of the structure. This term is the source of Jeans’ gravitational instability, which is the focus of our study here.

The linearised mass continuity equation expressed in terms of $\unicode[STIX]{x1D743}$ is

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{1}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D743})\end{eqnarray}$$

and the perturbation of the gravitational acceleration $\boldsymbol{g}_{1}$ satisfies the linearised Poisson equation

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{g}_{1}=-4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{1},\end{eqnarray}$$

keeping in mind that it is a curl-free vector field, namely

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D735}\times \boldsymbol{g}_{1}=\mathbf{0}.\end{eqnarray}$$

It is only with that additional constraint that (2.13) is equivalent to the more usual form $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{1}=4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{1}$ , where $\unicode[STIX]{x1D719}_{1}$ is the perturbed gravitational potential, since $\boldsymbol{g}_{1}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{1}$ .

In order to close this set of equations for the perturbations, we need an additional relation. Let us consider the case in which the time scale of the perturbations, i.e. the oscillation period if stable and growth time if unstable, is shorter than the time scale of the heat transfer between neighbouring fluid elements. The evolution of perturbations may then be considered as adiabatic, and it can be shown (e.g. Thompson Reference Thompson2006) that the equation expressing the absence of heat $\unicode[STIX]{x1D6FF}Q=0$ translates into the following relation between the Lagrangian variation of pressure $\unicode[STIX]{x1D6FF}p$ and of density $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}$ :

(2.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}p}{p_{0}}=\unicode[STIX]{x1D6FE}_{ad}\frac{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}},\end{eqnarray}$$

where, in general, $\unicode[STIX]{x1D6FE}_{ad}$ is not equal to the polytropic exponent $\unicode[STIX]{x1D6FE}$ of the polytropic equation of state of the equilibrium (e.g. Cox Reference Cox1980; Toci & Galli Reference Toci and Galli2015). Here however, we will take them equal because the purpose of this work is to focus on the gravitational instability of the Jeans type. This assumption switches off buoyancy, so that we get rid of convective instability, and as far as the stable part of the spectrum is concerned, since the Brunt–Väisälä frequency is now vanishing, the acoustic oscillations we will obtain are equivalent to the p-modes of stellar physics, without complicating the description with g-modes. Finally, rewriting (2.15) using the Eulerian variables $\unicode[STIX]{x1D70C}_{1}$ and $p_{1}$ (related to the Lagrangian description by $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D743}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}_{0}$ and $\unicode[STIX]{x1D6FF}p=p_{1}+\unicode[STIX]{x1D743}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p_{0}$ ) yields our closure relation

(2.16) $$\begin{eqnarray}p_{1}=c_{a}^{2}\unicode[STIX]{x1D70C}_{1}.\end{eqnarray}$$

We look for solutions by separation of variables, i.e. we consider that all quantities $Q$ (namely $\unicode[STIX]{x1D709}_{x},\unicode[STIX]{x1D709}_{y},\unicode[STIX]{x1D709}_{z},g_{1x},g_{1y},g_{1z},\unicode[STIX]{x1D70C}_{1}$ and $p_{1}$ ) are of the form $Q(x,y,z,t)=\hat{Q}(x)Y(y)Z(z)T(t)$ since we are considering a planar geometry. Naturally, the respectively space- and time-translational invariance along the $y$ , $z$ and $t$ dimensions yields exponential dependences for $Y$ , $Z$ and $T$ . Also, given that we are considering a planar geometry with an infinite extent transverse to the stratified $x$ -direction, the $y$ - and $z$ -directions are equivalent, so that we can rotate the axes in the ( $y$ , $z$ ) plane in order to remove any dynamics in the $z$ -direction. Therefore, without loss of generality we take $k_{z}=0$ and we have

(2.17) $$\begin{eqnarray}Q(x,y,z,t)=\hat{Q}(x)\text{e}^{\text{i}(k_{y}y-\unicode[STIX]{x1D714}t)}.\end{eqnarray}$$

2.3 Boundary conditions and symmetry conditions

The main difficulty in studying gravitational instability is that the differential equations satisfied by the perturbations are of order four. This is due to the term $\boldsymbol{g}_{1}$ which, with (2.13), turns the force operator (2.11) into an integro-differential operator. Dealing with fourth-order equations, we need to impose four conditions. In this paper we consider boundary conditions (BC) and symmetry conditions commonly used in the literature. We briefly summarise them here and adapt them to our notations. If need be, one may consult for example Goldreich & Lynden-Bell (Reference Goldreich and Lynden-Bell1965), Elmegreen & Elmegreen (Reference Elmegreen and Elmegreen1978) or Kim et al. (Reference Kim, Kim, Seo and Hong2012) for more details.

Due to the reflection symmetry of the governing equations with respect to the $x=0$ plane, the general solution is a superposition of symmetrical and antisymmetrical modes, which may be considered separately. It has been shown that antisymmetric modes are stable as far as gravitational instability is concerned. Therefore we will focus on symmetric modes, for which

(2.18) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\hat{\unicode[STIX]{x1D709}}_{x}(x=0)=0,\\ {\hat{g}}_{1x}(x=0)=0,\end{array}\right\}\end{eqnarray}$$

where we recall that quantities with a hat are defined through (2.17). We will refer to (2.18) as the ‘symmetry conditions’.

Another condition is obtained as follows. Apply the divergence theorem to the linearised Poisson equation for an infinitesimally thin shell containing the boundary layer. Compute the gravitational acceleration outside the slab by solving Laplace’s equation (since the external fluid is assumed to remain unperturbed), using the fact that it should not diverge at infinity. This results in the constraint

(2.19) $$\begin{eqnarray}{\hat{g}}_{1x}(x_{b})-\text{i}{\hat{g}}_{1y}(x_{b})=4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{b}\hat{\unicode[STIX]{x1D709}}_{x}(x_{b}).\end{eqnarray}$$

Finally, we consider two complementary types of BC. First, the rigid BC corresponds to the case in which the surface at the boundary does not move, i.e. the Lagrangian displacement vanishes:

(2.20) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D709}}_{x}(x_{b})=0\quad (\text{rigid BC}).\end{eqnarray}$$

Second, on the contrary, we treat the boundary as a free surface, but we require continuity of pressure throughout the dynamics, i.e. that the Lagrangian variation of pressure $\unicode[STIX]{x1D6FF}p=p_{1}+\unicode[STIX]{x1D743}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p_{0}$ vanishes at the boundary. With the hydrostatic equilibrium (2.1), the fourth condition in this free BC case is thus

(2.21) $$\begin{eqnarray}\hat{p}_{1}(x_{b})=-\unicode[STIX]{x1D70C}_{0}(x_{B})g_{0}(x_{b})\hat{\unicode[STIX]{x1D709}}_{x}(x_{b})\quad (\text{free BC}).\end{eqnarray}$$

2.4 The eigenvalue problem

Our setting constitutes the closed set of (2.10) and (2.12)–(2.14) and (2.16), together with the conditions (2.18)–(2.21). With the form (2.17) of the perturbations, this set becomes

(2.22) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\hat{\unicode[STIX]{x1D70C}}_{1}=-(\unicode[STIX]{x1D70C}_{0}\hat{\unicode[STIX]{x1D709}}_{x})^{\prime }-\unicode[STIX]{x1D70C}_{0}\text{i}k_{y}\hat{\unicode[STIX]{x1D709}}_{y},\\ -\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D714}^{2}\hat{\unicode[STIX]{x1D709}}_{x}=-c_{a}^{2}\hat{\unicode[STIX]{x1D70C}}_{1}^{\prime }+({\hat{g}}_{0}-(c_{a}^{2})^{\prime })\hat{\unicode[STIX]{x1D70C}}_{1}+\unicode[STIX]{x1D70C}_{0}{\hat{g}}_{1x},\\ -\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D714}^{2}\hat{\unicode[STIX]{x1D709}}_{y}=-\text{i}k_{y}c_{a}^{2}\hat{\unicode[STIX]{x1D70C}}_{1}+\unicode[STIX]{x1D70C}_{0}{\hat{g}}_{1y},\\ {\hat{g}}_{1x}^{\prime }+\text{i}k_{y}{\hat{g}}_{1y}=-4\unicode[STIX]{x03C0}G\hat{\unicode[STIX]{x1D70C}}_{1},\\ {\hat{g}}_{1x}=\left({\displaystyle \frac{{\hat{g}}_{1y}}{\text{i}k_{y}}}\right)^{\prime },\end{array}\right\}\end{eqnarray}$$

where a prime denotes differentiation with respect to $x$ , $^{\prime }=\text{d}/\text{d}x$ . From top to bottom these are respectively: mass continuity, the $x$ and $y$ components of momentum conservation, Poisson equation and the curl-free condition. Together with the boundary and symmetry conditions, this constitutes an eigenvalue problem on $\unicode[STIX]{x1D714}^{2}$ . The ultimate goal of this work is to compute the negative eigenvalues since they give the growth rate of the gravitational instability.

3.1 Reformulated equations of motion

First of all, let us work with dimensionless variables. From the homogeneous case (e.g. Thompson Reference Thompson2006), it is well known that a key length scale in gravitational instability is given by the critical Jeans wavenumber

(3.1) $$\begin{eqnarray}k_{J}(x)\equiv \sqrt{\frac{4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{0}(x)}{c_{a}^{2}(x)}},\end{eqnarray}$$

which marks the balance between pressure and gravity. In this definition $k_{J}$ is position dependent since the equilibrium density profile $\unicode[STIX]{x1D70C}_{0}$ is not assumed homogeneous here. We use the central value (subscripts $c$ ) of this quantity, noted as

(3.2) $$\begin{eqnarray}k_{Jc}\equiv k_{J}(x=0)=\sqrt{\frac{4\unicode[STIX]{x03C0}G}{\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D70C}_{c}^{2-\unicode[STIX]{x1D6FE}}},\end{eqnarray}$$

where we used definition (2.4) of $c_{a}^{2}$ in the last equality, together with the central value of the density $\unicode[STIX]{x1D70C}_{c}\equiv \unicode[STIX]{x1D70C}_{0}(x=0)$ , to define the following dimensionless parameters:

(3.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\bar{x}\equiv k_{Jc}x,\quad \unicode[STIX]{x1D70C}\equiv {\displaystyle \frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}_{c}}},\quad \bar{\unicode[STIX]{x1D714}}^{2}\equiv {\displaystyle \frac{\unicode[STIX]{x1D714}^{2}}{4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{c}}},\quad \bar{k}_{y}\equiv {\displaystyle \frac{k_{y}}{k_{Jc}}},\\ \unicode[STIX]{x1D713}\equiv \unicode[STIX]{x1D70C}k_{Jc}\hat{\unicode[STIX]{x1D709}}_{x},\quad {\mathcal{R}}\equiv \unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-2}{\displaystyle \frac{\hat{\unicode[STIX]{x1D70C}}_{1}}{\unicode[STIX]{x1D70C}_{c}}},\quad {\mathcal{G}}_{x}\equiv {\displaystyle \frac{k_{Jc}}{4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{c}}}{\hat{g}}_{1x},\quad {\mathcal{G}}_{y}\equiv {\displaystyle \frac{k_{Jc}}{4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{c}}}\text{i}\bar{k}_{y}{\hat{g}}_{1y}.\end{array}\right\}\end{eqnarray}$$

All these definitions are straightforward, except that of the variable ${\mathcal{R}}$ (this letter referring to ‘rho’, from $\unicode[STIX]{x1D70C}_{1}$ ). We choose to work with ${\mathcal{R}}\equiv \unicode[STIX]{x1D70C}_{0}^{\unicode[STIX]{x1D6FE}-2}\hat{\unicode[STIX]{x1D70C}}_{1}/\unicode[STIX]{x1D70C}_{c}^{\unicode[STIX]{x1D6FE}-1}$ rather than with the more natural $\hat{\unicode[STIX]{x1D70C}}_{1}/\unicode[STIX]{x1D70C}_{c}$ , because we notice that, using the definition (2.4) of $c_{a}^{2}$ and the hydrostatic equilibrium (2.5), we have the relation

(3.4) $$\begin{eqnarray}c_{a}^{2}\hat{\unicode[STIX]{x1D70C}}_{1}^{\prime }+(-{\hat{g}}_{0}+(c_{a}^{2})^{\prime })\hat{\unicode[STIX]{x1D70C}}_{1}=\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70C}_{0}(\unicode[STIX]{x1D70C}_{0}^{\unicode[STIX]{x1D6FE}-2}\hat{\unicode[STIX]{x1D70C}}_{1})^{\prime },\end{eqnarray}$$

which turns the second equation in (2.22) into a more compact form, such that we may rewrite the whole system of equations into a matrix form that is by blocks. Indeed, with variables (3.3), the system (2.22) reads

(3.5) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\bar{x}}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\\ {\mathcal{R}}\\ {\mathcal{G}}_{y}\end{array}\right)=\left(\begin{array}{@{}cccc@{}}0 & 0 & {\displaystyle \frac{\bar{k}_{y}^{2}}{\bar{\unicode[STIX]{x1D714}}^{2}}}\unicode[STIX]{x1D70C}(\bar{x})-\unicode[STIX]{x1D70C}(\bar{x})^{2-\unicode[STIX]{x1D6FE}} & {\displaystyle \frac{\unicode[STIX]{x1D70C}(\bar{x})}{\bar{\unicode[STIX]{x1D714}}^{2}}}\\ 0 & 0 & -\unicode[STIX]{x1D70C}(\bar{x})^{2-\unicode[STIX]{x1D6FE}} & -1\\ {\displaystyle \frac{\bar{\unicode[STIX]{x1D714}}^{2}}{\unicode[STIX]{x1D70C}(\bar{x})}} & 1 & 0 & 0\\ 0 & -\bar{k}_{y}^{2} & 0 & 0\end{array}\right)\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\\ {\mathcal{R}}\\ {\mathcal{G}}_{y}\end{array}\right),\end{eqnarray}$$

where we plugged the third equation of (2.22) into its first equation. This is already better than the initial formulation, but it is still very complicated, because we do not know what the coefficients in this matrix equation look like explicitly. Indeed, to do so we would need to solve for $\unicode[STIX]{x1D70C}(\bar{x})$ , i.e. integrate the Lane–Emden equation (2.6), but as stated in § 2.1, this has simple analytic solutions only for the special cases $\unicode[STIX]{x1D6FE}=2/3$ , 1, 2 and $\infty$ . And even in these cases, in which the equilibrium state is simple, the equation above, governing the evolution of the perturbations, is not simple. For example with $\unicode[STIX]{x1D6FE}=2$ we have $\unicode[STIX]{x1D70C}(\bar{x})=\cos (\bar{x})$ , and $\unicode[STIX]{x1D6FE}=1$ we have $\unicode[STIX]{x1D70C}(\bar{x})=\cosh (\bar{x}/\sqrt{2})^{-2}$ , so that the differential equation (3.5) is disappointingly complicated, especially for arbitrary  $\unicode[STIX]{x1D6FE}$ .

To improve this, let us do the change of variables (illustrated in figure 1)

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D70C}=\frac{\unicode[STIX]{x1D70C}_{0}(\bar{x})}{\unicode[STIX]{x1D70C}_{c}},\end{eqnarray}$$

which is possible because the equilibrium density profile $\unicode[STIX]{x1D70C}_{0}$ is monotonic so the change $\bar{x}\leftrightarrow \unicode[STIX]{x1D70C}$ is bijective. Interestingly, this change of variables is implicit since we do not need to know the functional form of $\unicode[STIX]{x1D70C}_{0}(\bar{x})$ . Indeed, using the chain rule, $\unicode[STIX]{x1D70C}^{\prime }$ is simple to compute: with the dimensionless variables (3.3), the Lane–Emden relation (2.6) becomes

(3.7) $$\begin{eqnarray}(\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-1})^{\prime \prime }+(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}=0,\end{eqnarray}$$

and multiplying by $(\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-1})^{\prime }$ , it can be rewritten as

(3.8) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\bar{x}}\left[\frac{1}{2}((\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-1})^{\prime })^{2}+\frac{(\unicode[STIX]{x1D6FE}-1)^{2}}{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}\right]=0.\end{eqnarray}$$

We integrate this, determining the constant of integration by considering profiles that are flat at the centre, i.e. $\unicode[STIX]{x1D70C}^{\prime }(\bar{x}=0)=0$ , which is a natural choice. Then we get

(3.9) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D70C}}{\text{d}\bar{x}}=-\sqrt{\frac{2}{\unicode[STIX]{x1D6FE}}}\unicode[STIX]{x1D70C}^{2-\unicode[STIX]{x1D6FE}}\sqrt{1-\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}},\end{eqnarray}$$

where we took the solution with negative sign because the profile is decreasing with $\bar{x}$ . Equation (3.5) can thus be rewritten as

(3.10) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D70C}}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\\ {\mathcal{R}}\\ {\mathcal{G}}_{y}\end{array}\right)=-\sqrt{\frac{\unicode[STIX]{x1D6FE}}{2}}\frac{\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-2}}{\sqrt{1-\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}}}\left(\begin{array}{@{}cccc@{}}0 & 0 & {\displaystyle \frac{\bar{k}_{y}^{2}}{\bar{\unicode[STIX]{x1D714}}^{2}}}\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}^{2-\unicode[STIX]{x1D6FE}} & {\displaystyle \frac{\unicode[STIX]{x1D70C}}{\bar{\unicode[STIX]{x1D714}}^{2}}}\\ 0 & 0 & -\unicode[STIX]{x1D70C}^{2-\unicode[STIX]{x1D6FE}} & -1\\ {\displaystyle \frac{\bar{\unicode[STIX]{x1D714}}^{2}}{\unicode[STIX]{x1D70C}}} & 1 & 0 & 0\\ 0 & -\bar{k}_{y}^{2} & 0 & 0\end{array}\right)\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\\ {\mathcal{R}}\\ {\mathcal{G}}_{y}\end{array}\right).\end{eqnarray}$$

Comparing (3.5) and (3.10) the modification may not look spectacular, but it is in fact a huge progress: equation (3.5) is very complicated because it has complicated coefficients that are not even explicit (in its variable $\bar{x}$ ), while (3.10) contains simple coefficients that are explicit (in its variable $\unicode[STIX]{x1D70C}$ ).

Finally, the occurrence of the square root $\sqrt{1-\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}}}$ in (3.10) calls for another change of variables (illustrated in figure 1), namely

(3.11) $$\begin{eqnarray}\sin \unicode[STIX]{x1D703}=\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}/2},\end{eqnarray}$$

which makes sense since $\unicode[STIX]{x1D70C}\in [0,1]$ by definition, and as we shall see below, this will indeed facilitate greatly our calculations. Note that, in order to lighten the expressions, we will often use the shorthand notations

(3.12) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}s\equiv \sin \unicode[STIX]{x1D703},\\ c\equiv \cos \unicode[STIX]{x1D703}.\end{array}\right\}\end{eqnarray}$$

Therefore, in this paper we will work with the following form of the system (2.22):

(3.13) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D703}}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\\ {\mathcal{R}}\\ {\mathcal{G}}_{y}\end{array}\right)=-\sqrt{\frac{2}{\unicode[STIX]{x1D6FE}}}\left(\begin{array}{@{}cccc@{}}0 & 0 & {\displaystyle \frac{\bar{k}_{y}^{2}}{\bar{\unicode[STIX]{x1D714}}^{2}}}s-s^{(2/\unicode[STIX]{x1D6FE})-1} & {\displaystyle \frac{s}{\bar{\unicode[STIX]{x1D714}}^{2}}}\\ 0 & 0 & -s^{(2/\unicode[STIX]{x1D6FE})-1} & -s^{1-(2/\unicode[STIX]{x1D6FE})}\\ \bar{\unicode[STIX]{x1D714}}^{2}s^{1-(4/\unicode[STIX]{x1D6FE})} & s^{1-(2/\unicode[STIX]{x1D6FE})} & 0 & 0\\ 0 & -\bar{k}_{y}^{2}s^{1-(2/\unicode[STIX]{x1D6FE})} & 0 & 0\end{array}\right)\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\\ {\mathcal{R}}\\ {\mathcal{G}}_{y}\end{array}\right).\end{eqnarray}$$

Let us now define the matrix $\unicode[STIX]{x1D643}$ such that

(3.14) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}{\mathcal{R}}(\unicode[STIX]{x1D703})\\ {\mathcal{G}}_{y}(\unicode[STIX]{x1D703})\end{array}\right)=\unicode[STIX]{x1D643}(\unicode[STIX]{x1D703})\left(\begin{array}{@{}c@{}}{\mathcal{R}}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\\ {\mathcal{G}}_{y}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\end{array}\right),\end{eqnarray}$$

i.e. $\unicode[STIX]{x1D643}$ describes the variation with position $\unicode[STIX]{x1D703}$ of $({\mathcal{R}},{\mathcal{G}}_{y})$ with respect to their values at the centre ( $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ ). We will use subscripts $a=1,2,3$ and 4 to denote its coefficients, namely

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D643}=\left(\begin{array}{@{}cc@{}}h_{1} & h_{2}\\ h_{3} & h_{4}\end{array}\right)\quad \text{where }h_{a}\equiv h_{a}(\unicode[STIX]{x1D703};\bar{\unicode[STIX]{x1D714}}^{2},\bar{k}_{y}^{2}).\end{eqnarray}$$

The matrix $\unicode[STIX]{x1D643}$ characterises completely the eigenfunctions since, by the definition (3.14), $\unicode[STIX]{x1D643}(\unicode[STIX]{x1D703})$ gives the value of ${\mathcal{R}}$ and ${\mathcal{G}}_{y}$ at a given  $\unicode[STIX]{x1D703}$ , and the two other eigenfunctions can be deduced, either in integral form by integrating the first two rows of (3.13), namely

(3.16) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\end{array}\right)=-\sqrt{\frac{2}{\unicode[STIX]{x1D6FE}}}\int _{\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}}\text{d}\unicode[STIX]{x1D703}_{1}\left(\begin{array}{@{}cc@{}}{\displaystyle \frac{\bar{k}_{y}^{2}}{\bar{\unicode[STIX]{x1D714}}^{2}}}s-s^{(2/\unicode[STIX]{x1D6FE})-1} & {\displaystyle \frac{s}{\bar{\unicode[STIX]{x1D714}}^{2}}}\\ -s^{(2/\unicode[STIX]{x1D6FE})-1} & -s^{1-(2/\unicode[STIX]{x1D6FE})}\end{array}\right)\unicode[STIX]{x1D643}(\unicode[STIX]{x1D703}_{1})\left(\begin{array}{@{}c@{}}{\mathcal{R}}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\\ {\mathcal{G}}_{y}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\end{array}\right)\end{eqnarray}$$

(where we have used the symmetry conditions (3.19) below) or expressed with a derivative by inverting the last two rows of (3.13), namely

(3.17) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\\ {\mathcal{G}}_{x}\end{array}\right)=\sqrt{\frac{\unicode[STIX]{x1D6FE}}{2}}\frac{s^{(4/\unicode[STIX]{x1D6FE})-1}}{\bar{\unicode[STIX]{x1D714}}^{2}\bar{k}_{y}^{2}}\left(\begin{array}{@{}cc@{}}-\bar{k}_{y}^{2} & -1\\ 0 & \bar{\unicode[STIX]{x1D714}}^{2}s^{-(2/\unicode[STIX]{x1D6FE})}\end{array}\right)\unicode[STIX]{x1D643}^{\prime }(\unicode[STIX]{x1D703})\left(\begin{array}{@{}c@{}}{\mathcal{R}}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\\ {\mathcal{G}}_{y}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\end{array}\right).\end{eqnarray}$$

3.2 Reformulated boundary conditions and symmetry conditions

Let us start with two preliminary remarks. First, from now on, subscripts $b$ will indicate values taken at the boundary, notably $\bar{x}_{b},\unicode[STIX]{x1D70C}_{b}$ and $\unicode[STIX]{x1D703}_{b}$ . Secondly, note that an interesting feature of the change of variable (3.6) is that, thanks to it, in this work we will derive the eigenvalues without integrating the Lane–Emden equation (which is possible only in a few cases), i.e. without having to obtain the equilibrium density as a function of position $\unicode[STIX]{x1D70C}_{0}(x)$ . What matters is not where the BC is imposed (position $\bar{x}_{b}$ ), but the density $\unicode[STIX]{x1D70C}_{b}$ at the position where it is imposed. This is why, in the end, the explicit formulæ (4.21) are expressed in terms of $\unicode[STIX]{x1D70C}_{b}$ . Hence, given the importance of $\unicode[STIX]{x1D70C}_{b}$ , let us see explicitly how it is related to the external pressure applied on the structure. With the equation of state (2.3) and from the definition (2.8) of $x_{b}$ we have

(3.18) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{b}=\left(\frac{p_{ext}}{p_{c}}\right)^{1/\unicode[STIX]{x1D6FE}},\end{eqnarray}$$

where $p_{c}\equiv p_{0}(x=0)=\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}_{c}^{\unicode[STIX]{x1D6FE}}$ is the central value of the equilibrium pressure.

We now update our formulation of § 2.3 using the new variables (3.3). As far as symmetry conditions are concerned, it is clear that (2.18) becomes

(3.19) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)=0,\\ {\mathcal{G}}_{x}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)=0.\end{array}\right\}\end{eqnarray}$$

Then, the BC (2.19) reads

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D703}_{b})-{\mathcal{G}}_{x}(\unicode[STIX]{x1D703}_{b})+\bar{k}_{y}^{-1}{\mathcal{G}}_{y}(\unicode[STIX]{x1D703}_{b})=0.\end{eqnarray}$$

We also combine (2.20) and (2.21) into a single expression,

(3.21) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D703}_{b})+\unicode[STIX]{x1D6E5}_{BC}{\mathcal{R}}(\unicode[STIX]{x1D703}_{b})=0,\end{eqnarray}$$

where

(3.22) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{BC}\equiv \unicode[STIX]{x1D6FF}_{BC}\frac{4\unicode[STIX]{x03C0}G\unicode[STIX]{x1D70C}_{c}}{k_{Jc}}\frac{\unicode[STIX]{x1D70C}_{b}}{g_{0}(\bar{x}_{b})},\end{eqnarray}$$

and we have defined the parameter

(3.23) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{BC}\equiv \left\{\begin{array}{@{}ll@{}}1\quad & \text{for free BC},\\ 0\quad & \text{for rigid BC},\end{array}\right.\end{eqnarray}$$

which enables us to switch between the two types of BC and spares us the pain of writing separate equations. With the definition of the sound speed (2.4) and the hydrostatic equilibrium (2.5) we may rewrite $\unicode[STIX]{x1D6E5}_{BC}$ using the $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D703}$ variables, as

(3.24) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{BC}=-\unicode[STIX]{x1D6FF}_{BC}\sqrt{\frac{\unicode[STIX]{x1D6FE}}{2}}\frac{\unicode[STIX]{x1D70C}_{b}}{\sqrt{1-\unicode[STIX]{x1D70C}_{b}^{\unicode[STIX]{x1D6FE}}}}=-\unicode[STIX]{x1D6FF}_{BC}\sqrt{\frac{\unicode[STIX]{x1D6FE}}{2}}\frac{(\sin \unicode[STIX]{x1D703}_{b})^{2/\unicode[STIX]{x1D6FE}}}{\cos \unicode[STIX]{x1D703}_{b}}.\end{eqnarray}$$

3.3 Eigenvalue equation

From the BC above, let us derive a single scalar equation which constrains the eigenvalue parameter $\bar{\unicode[STIX]{x1D714}}^{2}$ . We rewrite the boundary conditions (3.20) and (3.21) in matrix form as

(3.25) $$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}1 & -1\\ 1 & 0\end{array}\right)\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D703}_{b})\\ {\mathcal{G}}_{x}(\unicode[STIX]{x1D703}_{b})\end{array}\right)+\left(\begin{array}{@{}cc@{}}0 & 1/\bar{k}_{y}\\ \unicode[STIX]{x1D6E5}_{BC} & 0\end{array}\right)\left(\begin{array}{@{}c@{}}{\mathcal{R}}(\unicode[STIX]{x1D703}_{b})\\ {\mathcal{G}}_{y}(\unicode[STIX]{x1D703}_{b})\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right),\end{eqnarray}$$

and using (3.14) and (3.17), we rewrite this with the coefficients of matrix $\unicode[STIX]{x1D643}$ and their derivatives, namely

(3.26) $$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}q_{1}(\unicode[STIX]{x1D703}_{b}) & q_{2}(\unicode[STIX]{x1D703}_{b})\\ q_{3}(\unicode[STIX]{x1D703}_{b}) & q_{4}(\unicode[STIX]{x1D703}_{b})\end{array}\right)\left(\begin{array}{@{}c@{}}{\mathcal{R}}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\\ {\mathcal{G}}_{y}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\right)\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right),\end{eqnarray}$$

where (all quantities below are evaluated at $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{b}$ )

(3.27) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}q_{1}=\bar{k}_{y}^{2}h_{1}^{\prime }+[1+s^{-(2/\unicode[STIX]{x1D6FE})}\bar{\unicode[STIX]{x1D714}}^{2}]h_{3}^{\prime }-\sqrt{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}}s^{1-(4/\unicode[STIX]{x1D6FE})}\bar{\unicode[STIX]{x1D714}}^{2}\bar{k}_{y}h_{3},\\ q_{2}=\bar{k}_{y}^{2}h_{2}^{\prime }+[1+s^{-(2/\unicode[STIX]{x1D6FE})}\bar{\unicode[STIX]{x1D714}}^{2}]h_{4}^{\prime }-\sqrt{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}}s^{1-(4/\unicode[STIX]{x1D6FE})}\bar{\unicode[STIX]{x1D714}}^{2}\bar{k}_{y}h_{4},\\ q_{3}=\bar{k}_{y}^{2}h_{1}^{\prime }+h_{3}^{\prime }-\sqrt{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}}s^{1-(4/\unicode[STIX]{x1D6FE})}\bar{\unicode[STIX]{x1D714}}^{2}\unicode[STIX]{x1D6E5}_{BC}\bar{k}_{y}^{2}h_{1},\\ q_{4}=\bar{k}_{y}^{2}h_{2}^{\prime }+h_{4}^{\prime }-\sqrt{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}}s^{1-(4/\unicode[STIX]{x1D6FE})}\bar{\unicode[STIX]{x1D714}}^{2}\unicode[STIX]{x1D6E5}_{BC}\bar{k}_{y}^{2}h_{2}.\end{array}\right\}\end{eqnarray}$$

The matrix in (3.26) must be non-invertible for non-trivial solutions to exist: its determinant must vanish,

(3.28) $$\begin{eqnarray}q_{1}q_{4}-q_{3}q_{2}=0.\end{eqnarray}$$

This relation is the starting point of our analysis. It is an equation on $\bar{\unicode[STIX]{x1D714}}^{2}$ , and its solutions are the eigenvalues that we are looking for. For a given $\bar{k}_{y}$ these solutions constitute the discrete spectrum we are looking for, and we will refer to it as the ‘eigenvalue equation’. In the remainder of this paper, we detail how we solve (3.28) and obtain the growth rate of gravitational instability.

4.1 Step I: computing $\sum _{n=1}^{\infty }1/\bar{\unicode[STIX]{x1D714}}_{n}^{2}$

Let us first focus on the series appearing in formula (4.8). An outline is given here, and details are provided in appendix A. Interestingly, at first sight, relation (4.8) does not look like a strategic path to take at all. Indeed, it seems that in order to solve one fourth-order problem (computing $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ ), we will have to solve an infinite number of fourth-order problems (computing the $\bar{\unicode[STIX]{x1D714}}_{n\geqslant 1}^{2}$ values) and in addition we will need to compute a series, which could in principle be extremely difficult too. However, the key point is to notice that, in fact, computing the $\bar{\unicode[STIX]{x1D714}}_{n\geqslant 1}^{2}$ values constitutes a much easier problem than the initial one because, as in stellar oscillation theory, we may compute these high-order modes with extremely high accuracy in the so-called Cowling approximation (Cowling Reference Cowling1941). (We borrow here the terminology from asteroseismology.) It consists in neglecting the perturbation of the gravitational field $\boldsymbol{g}_{1}$ in the linearised momentum conservation (2.10), in which case the problem becomes only second order, instead of fourth. As a reminder of this, we will use the notation

(4.9) $$\begin{eqnarray}\mathscr{C}\equiv \mathop{\sum }_{n=1}^{\infty }\frac{1}{\bar{\unicode[STIX]{x1D714}}_{n}^{2}},\end{eqnarray}$$

where the letter $\mathscr{C}$ stands for ‘Cowling’.

Although by far much simpler than in the full case, the equations in the Cowling approximation remain difficult to solve exactly. The second key point now is that we know that the $\bar{\unicode[STIX]{x1D714}}_{n\geqslant 1}^{2}$ values are large, with even $\bar{\unicode[STIX]{x1D714}}_{n}^{2}\rightarrow \infty$ as $n\rightarrow \infty$ , so that when we want to solve our eigenvalue problem looking for the high-order modes, we know that the parameter $\bar{\unicode[STIX]{x1D714}}$ will be large. Therefore, we may use the Wentzel–Kramers–Brillouin (WKB) method treating $1/\bar{\unicode[STIX]{x1D714}}$ as a small parameter.

With the above two ideas (Cowling & WKB) we derive in appendix A very precise formulæ for the eigenfunctions ${\mathcal{R}}$ and $\unicode[STIX]{x1D713}$ (the only relevant ones in the Cowling approximation) of (3.5), and the eigenvalues $\bar{\unicode[STIX]{x1D714}}_{n\geqslant 1}^{2}$ . In figure 2(a,b), we show an example of a plot of ${\mathcal{R}}$ and $\unicode[STIX]{x1D713}$ using these analytic expressions, compared to the numerical solution. For a very wide range of the parameters $(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D70C}_{b},n,\bar{k}_{y})$ our WKB solution gives excellent results. In figure 2(c), we show the spectrum $\bar{\unicode[STIX]{x1D714}}_{n\geqslant 1}^{2}$ that we obtain with formula (A 26) and compare it to the spectrum solved numerically using a shooting method. As we can see again, both agree extremely well. Let us stress that in the latter plot, the numerical solution corresponds to the full solution of the initial eigenvalue problem (2.22), without the Cowling approximation. Therefore, this panel shows that it is not only the WKB solution itself that works very well (as illustrated in figure 2(a,b) for the eigenfunctions) but also the Cowling approximation, since the eigenvalues (A 26) were obtained by combining both.

Figure 2.

In appendix A we compute the stable part of the spectrum (i.e. harmonics $n\geqslant 1$ ). Here we compare our analytic expressions (red lines) with the numerical solution (black dots) for the eigenfunction ${\mathcal{R}}(\unicode[STIX]{x1D703})$ in (a), $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D703})$ in (b) and for the eigenvalues $\bar{\unicode[STIX]{x1D714}}_{n}^{2}$ given by (A 26) in (c) (shown are $n=1,2,3$ and 4). In (a,b), the $x$ -axis ranges from $\unicode[STIX]{x1D703}_{b}$ (the boundary) to $\unicode[STIX]{x03C0}/2$ (centre of the slab). In this example we chose the rigid BC (hence the fact that the curve for $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D703})$ reaches zero on the left and right of the plot), with values $\unicode[STIX]{x1D6FE}=1.3$ , $\unicode[STIX]{x1D70C}_{b}=0.2$ , $n=10$ and $\bar{k}_{y}=2.5$ . Normalisation has been chosen such that ${\mathcal{R}}(\unicode[STIX]{x1D703}_{b})=1$ . Grey dashed lines indicate the envelopes, such as the one appearing in front of the cosine in expression (A 17). Note that the numerical solution for the eigenfunctions here is done in the Cowling approximation, to show that the WKB solution gives very satisfying results (such an agreement between the analytic and numerical solutions happens for a very wide range of parameters). However, in (c), the analytic result is compared to the full numerical result, i.e. not in the Cowling approximation. Therefore, (c) shows that the approximate spectrum (A 26) fits the full numerical result very well, which validates our procedure of combining the Cowling approximation with the WKB method, and makes expression (4.10) for $\mathscr{C}$ very efficient.

Then finally, in appendix A, we compute the series (4.9) with the pleasant satisfaction of being able to express it in a simple closed form,

(4.10) $$\begin{eqnarray}\mathscr{C}=\frac{1}{2Ab_{BC}^{2}}[\unicode[STIX]{x03C0}b_{BC}\tanh (\unicode[STIX]{x03C0}b_{BC})^{2\unicode[STIX]{x1D6FF}_{BC}-1}+\unicode[STIX]{x1D6FF}_{BC}-1],\end{eqnarray}$$

where

(4.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}b_{BC}=\sqrt{{\displaystyle \frac{B\bar{k}_{y}^{2}+C+D_{BC}}{A}}},\\ A={\displaystyle \frac{\unicode[STIX]{x1D6FE}\unicode[STIX]{x03C0}^{2}}{2I_{-(1/\unicode[STIX]{x1D6FE})}(\unicode[STIX]{x1D703}_{b})}},\\ B=I_{2-(3/\unicode[STIX]{x1D6FE})}(\unicode[STIX]{x1D703}_{b}),\\ C={\displaystyle \frac{(\unicode[STIX]{x1D6FE}-3)}{8\unicode[STIX]{x1D6FE}}}[(5-\unicode[STIX]{x1D6FE})I_{1/\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D703}_{b})+(3\unicode[STIX]{x1D6FE}-5)I_{(1/\unicode[STIX]{x1D6FE})-2}(\unicode[STIX]{x1D703}_{b})],\\ D_{BC}=-2\sin (\unicode[STIX]{x1D703}_{b})^{(1/\unicode[STIX]{x1D6FE})-1}\cos (\unicode[STIX]{x1D703}_{b})\left(1+(\unicode[STIX]{x1D6FF}_{BC}-1){\displaystyle \frac{1+\unicode[STIX]{x1D6FE}}{4}}\right),\end{array}\right\}\end{eqnarray}$$

the subscripts ‘BC’ denoting quantities that depend on the type of BCs ( $\unicode[STIX]{x1D6FF}_{BC}=1$ for free BC, $0$ for rigid BC), and we use the notation

(4.12) $$\begin{eqnarray}I_{p}(\unicode[STIX]{x1D703})\equiv \int _{\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}}(\sin \unicode[STIX]{x1D703})^{p}\,\text{d}\unicode[STIX]{x1D703}=2^{p}\left[B_{(1-\cos \unicode[STIX]{x1D703})/2}\left(\frac{p+1}{2},\frac{p+1}{2}\right)-B_{1/2}\left(\frac{p+1}{2},\frac{p+1}{2}\right)\right],\end{eqnarray}$$

where the second equality comes from the definition of the incomplete beta function.

4.2 Step II: computing $a_{1}/a_{0}$

Let us now focus on the other term, $a_{1}/a_{0}$ , contributing to $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ in formula (4.8), details being provided in appendix B. First of all, it is interesting to note that thanks to (4.8), we need to compute only two of the coefficients $a_{k}$ of the power series (4.3). What is more, these first two are also the simplest. Then, from definition (4.4), it turns out that the only things needed are the functions $h_{a}^{0}$ and $h_{a}^{1}$ of expansion (4.1). That is, out of an a priori infinite number of functions, we need only eight (since $a=1,2,3$ and 4). Now, as detailed in appendix B, in order to obtain $h_{a}^{0}$ and $h_{a}^{1}$ , we first derive the differential equations they satisfy from the equation of motion (3.13). Then, recalling that for gravitational instability only large scale perturbations are relevant (i.e. small $\bar{k}_{y}$ , otherwise pressure balances gravity), we look for the solutions by expanding $h_{a}^{0}$ and $h_{a}^{1}$ in powers of $\bar{k}_{y}$ (cf. definition (4.13) below). The practical advantage of this expansion is that, while $h_{a}^{0}$ and $h_{a}^{1}$ satisfy fourth-order equations, the problem is actually reduced to a ( $n$ infinite) sequence of second-order problems that are fully solvable, namely (B 23) in appendix B. Therefore, we compute $h_{a}^{0}$ and $h_{a}^{1}$ exactly, but order by order in $\bar{k}_{y}$ . For this reason, in the end we obtain $a_{0}$ , $a_{1}$ , and consequently $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ , as power series in $\bar{k}_{y}$ . We show the general procedure in appendix B, but let us give here only the intermediary expressions necessary to derive $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ up to the second order in $\bar{k}_{y}$ . Denoting the expansions in $\bar{k}_{y}$ by

(4.13a-c ) $$\begin{eqnarray}h_{a}^{i}=\mathop{\sum }_{j=0}^{\infty }h_{a,j}^{i}(\bar{k}_{y}^{2})^{j},\quad a_{0}=\mathop{\sum }_{m=0}^{\infty }a_{0,m}(\bar{k}_{y})^{m}\quad \text{and}\quad a_{1}=\mathop{\sum }_{m=0}^{\infty }a_{1,m}(\bar{k}_{y})^{m},\end{eqnarray}$$

we obtain that the useful $h_{a,j}^{i}$ values are given by

(4.14) $$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}h_{1,0}^{0}(\unicode[STIX]{x1D703})\, & =\, & [Q_{1/\unicode[STIX]{x1D6FE}}^{1-(1/\unicode[STIX]{x1D6FE})}(0)]^{-1}y_{Q}(\unicode[STIX]{x1D703}),\\ h_{1,0}^{1}(\unicode[STIX]{x1D703})\, & =\, & {\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}\right)[Q_{1/\unicode[STIX]{x1D6FE}}^{1-(1/\unicode[STIX]{x1D6FE})}(0)]^{-1}\\ \, & \, & \times \,\left[y_{Q}(\unicode[STIX]{x1D703})\displaystyle \int _{\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}}\left({\displaystyle \frac{c}{s}}y_{Q}^{\prime }+y_{Q}\right){\displaystyle \frac{y_{P}}{s}}\,\text{d}\unicode[STIX]{x1D703}_{1}-y_{P}(\unicode[STIX]{x1D703})\displaystyle \int _{\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}}\left({\displaystyle \frac{c}{s}}y_{Q}^{\prime }+y_{Q}\right){\displaystyle \frac{y_{Q}}{s}}\,\text{d}\unicode[STIX]{x1D703}_{1}\right]\!,\end{array}\right\}\qquad\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}$ is the gamma function and we have defined

(4.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}y_{P}(\unicode[STIX]{x1D703})\equiv (\sin \unicode[STIX]{x1D703})^{1-(1/\unicode[STIX]{x1D6FE})}P_{1/\unicode[STIX]{x1D6FE}}^{1-(1/\unicode[STIX]{x1D6FE})}(\cos \unicode[STIX]{x1D703}),\\ y_{Q}(\unicode[STIX]{x1D703})\equiv (\sin \unicode[STIX]{x1D703})^{1-(1/\unicode[STIX]{x1D6FE})}Q_{1/\unicode[STIX]{x1D6FE}}^{1-(1/\unicode[STIX]{x1D6FE})}(\cos \unicode[STIX]{x1D703}),\end{array}\right\}\end{eqnarray}$$

with $P_{\unicode[STIX]{x1D706}}^{\unicode[STIX]{x1D707}}$ and $Q_{\unicode[STIX]{x1D706}}^{\unicode[STIX]{x1D707}}$ the associated Legendre functions of the first and second kind. The useful $a_{n,m}$ values are given by

(4.16) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}a_{0,1}={\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}cs^{2/\unicode[STIX]{x1D6FE}}\left((h_{1,0}^{0})^{\prime }+\unicode[STIX]{x1D6FF}_{BC}{\displaystyle \frac{s}{c}}h_{1,0}^{0}\right),\\ a_{1,0}=\sqrt{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}}s^{2/\unicode[STIX]{x1D6FE}}\left((h_{1,0}^{0})^{\prime }+\unicode[STIX]{x1D6FF}_{BC}{\displaystyle \frac{s}{c}}h_{1,0}^{0}\right),\end{array}\right\}\end{eqnarray}$$

for the linear expansion of $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ (cf. (4.20a ) below), while those for the quadratic expansion (cf. (4.20b ) below) read

(4.17) $$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}a_{0,2}\, & =\, & -\left({\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}\right)^{3/2}c\,s\,h_{1,0}^{0},\\ a_{1,1}\, & =\, & {\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}s^{2/\unicode[STIX]{x1D6FE}}\left[\vphantom{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}}s^{1-(2/\unicode[STIX]{x1D6FE})}\unicode[STIX]{x1D6FF}_{BC}h_{1,0}^{0}+c\left((h_{1,0}^{1})^{\prime }+\unicode[STIX]{x1D6FF}_{BC}{\displaystyle \frac{s}{c}}h_{1,0}^{1}\right)\right.\\ \, & \, & -\left.{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}(I_{1-(2/\unicode[STIX]{x1D6FE})}-I_{-1-(2/\unicode[STIX]{x1D6FE})})\left((h_{1,0}^{0})^{\prime }+\unicode[STIX]{x1D6FF}_{BC}{\displaystyle \frac{s}{c}}h_{1,0}^{0}\right)\right],\end{array}\right\}\end{eqnarray}$$

where all quantities are to be evaluated at $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{b}$ .

4.3 Step III: the growth rate $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$

Finally, we put together the results from steps I and II, and obtain explicit expressions for the growth rate. Now, the resulting expression for $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ from (4.8) comes out as an expansion in $\bar{k}_{y}$ , but this is in fact physically most relevant since gravitational instability precisely occurs for small $\bar{k}_{y}$ . More precisely, with the expansions (4.13) of $a_{0}$ and $a_{1}$ , and expanding $\mathscr{C}$ with a similar notation as

(4.18) $$\begin{eqnarray}\mathscr{C}=\mathop{\sum }_{m=0}^{\infty }\mathscr{C}_{m}(\bar{k}_{y})^{m},\end{eqnarray}$$

relation (4.8) becomes

(4.19) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D714}}_{0}^{2}(\bar{k}_{y})=\bar{\unicode[STIX]{x1D714}}_{0,1}^{2}\;\bar{k}_{y}+\bar{\unicode[STIX]{x1D714}}_{0,2}^{2}\;\bar{k}_{y}^{2}+\mathit{O}(\bar{k}_{y}^{3})\end{eqnarray}$$

(it turns out that $\bar{\unicode[STIX]{x1D714}}_{0,0}^{2}=0$ ), where

(4.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D714}}_{0,1}^{2}=-{\displaystyle \frac{a_{0,1}}{a_{1,0}}}, & \displaystyle\end{eqnarray}$$

(4.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D714}}_{0,2}^{2}=\left(\frac{a_{0,2}}{a_{0,1}}-\frac{a_{1,1}}{a_{1,0}}\right)\unicode[STIX]{x1D714}_{0,1}^{2}+\mathscr{C}_{0}\unicode[STIX]{x1D714}_{0,1}^{4}. & \displaystyle\end{eqnarray}$$

(4.21) $$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}\bar{\unicode[STIX]{x1D714}}_{0,1}^{2}\, & =\, & -\sqrt{{\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}(1-\unicode[STIX]{x1D70C}_{b}^{\unicode[STIX]{x1D6FE}})},\\ \bar{\unicode[STIX]{x1D714}}_{0,2}^{2}\, & =\, & {\displaystyle \frac{4}{\unicode[STIX]{x1D6FE}^{2}}}{\displaystyle \frac{c^{2}}{{\displaystyle \frac{c}{s}}y_{Q}^{\prime }+\unicode[STIX]{x1D6FF}_{BC}y_{Q}}}\left\{{\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}}(1+\unicode[STIX]{x1D6FF}_{BC})s^{-(2/\unicode[STIX]{x1D6FE})}y_{Q}-\left({\displaystyle \frac{c}{s}}y_{P}^{\prime }+\unicode[STIX]{x1D6FF}_{BC}y_{P}\right)\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}\right)\right.\\ \, & \, & \times \,\displaystyle \int _{\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}_{b}}\!\left({\displaystyle \frac{c}{s}}y_{Q}^{\prime }+y_{Q}\right)\!{\displaystyle \frac{y_{Q}}{s}}\,\text{d}\unicode[STIX]{x1D703}_{1}+\left({\displaystyle \frac{c}{s}}y_{Q}^{\prime }+\unicode[STIX]{x1D6FF}_{BC}y_{Q}\right)\!\left[\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{2}{\unicode[STIX]{x1D6FE}}}\right)\!\displaystyle \int _{\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}_{b}}\left({\displaystyle \frac{c}{s}}y_{Q}^{\prime }+y_{Q}\right)\!{\displaystyle \frac{y_{P}}{s}}\,\text{d}\unicode[STIX]{x1D703}_{1}\right.\\ \, & \, & +\,\left.\left.{\displaystyle \frac{1}{c}}(I_{-1-(2/\unicode[STIX]{x1D6FE})}-I_{1-(2/\unicode[STIX]{x1D6FE})})+{\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}}\mathscr{C}_{0}\right]\right\}.\end{array}\right\}\end{eqnarray}$$

We remind the reader that in this expression all quantities outside integrals are evaluated at the boundary $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{b}$ , and we recall the following definitions: $\unicode[STIX]{x1D70C}_{b}=(p_{ext}/p_{c})^{1/\unicode[STIX]{x1D6FE}}$ where $p_{c}$ is the equilibrium pressure at the centre and $p_{ext}\leqslant p_{c}$ is the external pressure; $\unicode[STIX]{x1D703}_{b}=\arcsin \unicode[STIX]{x1D70C}_{b}^{\unicode[STIX]{x1D6FE}/2}$ ; $s=\sin \unicode[STIX]{x1D703}$ ; $c=\cos \unicode[STIX]{x1D703}$ ; $y_{P}$ and $y_{Q}$ are defined in (4.15); $\unicode[STIX]{x1D6E4}$ is the gamma function; $\unicode[STIX]{x1D6FF}_{BC}=0$ for rigid BC and $\unicode[STIX]{x1D6FF}_{BC}=1$ for free BC; $I_{p}$ is defined in (4.12); and $\mathscr{C}_{0}=\mathscr{C}(\bar{k}_{y}=0)$ is given by formula (4.10). (We exclude $p_{ext}>p_{c}$ because in this case the condition (2.8) cannot be fulfilled and no proper equilibrium state exists.)

Let us now give a couple of comments on this formula, beginning with the linear order (first line above). First, this result is exact, in the sense that we did not need the Cowling approximation to obtain it as $\mathscr{C}$ does not show in the linear expansion. Second, it is interesting to note that this result, at the linear order, is independent of the BC. Indeed, while $a_{0,1}$ and $a_{1,0}$ both depend on the BC, that dependence cancels out when we take their ratio. Third, note that we have the following relations:

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D707}\equiv \int _{0}^{\bar{x}}\unicode[STIX]{x1D70C}(s)\,\text{d}s=\frac{1}{1-\unicode[STIX]{x1D6FE}}\int _{0}^{\bar{x}}(\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-1})^{\prime \prime }\,\text{d}s=\frac{1}{1-\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}-1})^{\prime }=\sqrt{\frac{2}{\unicode[STIX]{x1D6FE}}(1-\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}})},\end{eqnarray}$$

where we used the Lane–Emden equation (3.7) in the first equality, in the second equality we used the fact that the profiles we are considering are flat at $\bar{x}=0$ and in the third equality we used expression (3.9) for the derivative of $\unicode[STIX]{x1D70C}$ . From this, we see that $\bar{\unicode[STIX]{x1D714}}_{0,1}^{2}$ in (4.21) actually corresponds to the dimensionless column density (denoted $\unicode[STIX]{x1D707}$ in the rest of the paper) from the centre to the boundary of the slab, i.e. essentially to the mass (per unit area) of the structure. In other words, we recover the well-known asymptotic result that for long wavelengths (small $\bar{k}_{y}$ )

(4.23) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D714}}_{0}^{2}\sim -\unicode[STIX]{x1D707}\bar{k}_{y}\end{eqnarray}$$

i.e. putting back the physical dimensions

(4.24) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}^{2}\sim -2\unicode[STIX]{x03C0}G\unicode[STIX]{x1D6F4}_{0}k_{y},\end{eqnarray}$$

where

(4.25) $$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{0}=\frac{2\unicode[STIX]{x1D70C}_{c}}{k_{Jc}}\unicode[STIX]{x1D707}=\frac{2\unicode[STIX]{x1D70C}_{c}}{k_{Jc}}\int _{0}^{\bar{x}_{b}}\bar{\unicode[STIX]{x1D70C}}(\bar{x})\,\text{d}\bar{x}=\int _{-x_{b}}^{x_{b}}\unicode[STIX]{x1D70C}_{0}(x)\,\text{d}x\end{eqnarray}$$

is the usual, dimensionful column density.

About the quadratic order, we first point out that this result is exact as long as we take $\mathscr{C}_{0}$ from its definition, namely the value of $\mathscr{C}$ defined in (4.9) at $\bar{k}_{y}=0$ , because we have computed exactly the term $a_{1}/a_{0}$ order by order in $\bar{k}_{y}$ . It is only once we use expression (4.10) to express $\mathscr{C}_{0}$ explicitly that the result becomes approximate, since the latter was deduced using the Cowling approximation and the WKB method. However, we find numerically that $\mathscr{C}_{0}$ is small compared to the other terms, particularly when $\unicode[STIX]{x1D70C}_{b}$ approaches unity. Second, we write $\bar{\unicode[STIX]{x1D714}}_{0,2}^{2}$ in (4.21) in this manner in order to highlight the first term in the brackets, which dominates in the high pressure limit $\unicode[STIX]{x1D70C}_{b}\rightarrow 1$ . Third, in figure 3, we show a typical example: the linear expansion is plotted in green, the quadratic expansion in orange and the numerical solution in black dots. As we can see, the quadratic formula is in excellent agreement with the numerical solution for small $\bar{k}_{y}$ , namely for $\bar{k}_{y}$ from zero to a value close to $\bar{k}_{y,max}$ , the wavenumber of fastest growth (i.e. for $0\leqslant \bar{k}_{y}\lesssim 0.5$ in the example shown). Therefore, from this quadratic formula one can already get good approximations of the maximum growth rate and the corresponding wavenumber, and fair approximations in the full unstable range, i.e. up to the critical wavenumber $\bar{k}_{crit}$ at which the fundamental frequency $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ reaches zero again. We will explore this idea further in § 5.

Figure 3.

At the top of this figure we reproduce the formula for the fundamental eigenfrequency $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ as an expansion in $\bar{k}_{y}$ up to the fourth order. The colour coding used in this expression corresponds to the different curves in the plot: green, orange, red and blue curves are respectively the linear, quadratic, cubic and quartic expansions. The numerical solution is shown in black dots. In this plot we naturally see that, as we increase the order, we increase the range of $\bar{k}_{y}$ over which the analytical result fits the numerical solution. In this example we chose free BC, $\unicode[STIX]{x1D6FE}=12$ and $\unicode[STIX]{x1D70C}_{b}=0.6$ , and the full unstable range (the range of $\bar{k}_{y}$ for which $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ is negative) is perfectly matched by the fourth-order expansion. Note that the lower-order expansions are all the better since $\unicode[STIX]{x1D70C}_{b}$ is close to one (i.e. high external pressure), and the second order (4.19) becomes excellent for $\unicode[STIX]{x1D70C}_{b}\gtrsim 0.8$ in general.

Figure 4.

Results fixing $\unicode[STIX]{x1D6FE}=1$ and varying $\unicode[STIX]{x1D70C}_{b}$ . Solid lines are the analytic solutions at fourth order, and dashed lines are the numerical solutions. In (a) we use rigid BC, while in (b) we use free BC. We indicate with an arrow how the slope at the origin varies as we increase the parameter $\unicode[STIX]{x1D70C}_{b}$ . Here we increase the values of $\unicode[STIX]{x1D70C}_{b}$ from 0.4 to 0.9 by steps of 0.1. Formally this evolution can be clearly seen with the expression of $\bar{\unicode[STIX]{x1D714}}_{0,1}^{2}$ in (4.21). The light grey branch in (b) is explained in § 4.4.

4.4 Beyond the quadratic formula

An important aspect of the method presented in this paper is that it is systematic, in the sense that if need be, it is straightforward to increase the accuracy of expression (4.19) by increasing the order of the expansions. (To be more precise, with formula (4.8) we have an ‘intrinsic’ limitation, given by the accuracy of the Cowling approximation. To increase the accuracy even further, one may consider the first three root-coefficient relations, rather than just the first two as we did here. Doing so, we are limited by the accuracy of the Cowling approximation for one harmonic higher, and since that approximation improves for high harmonics, the total accuracy will indeed be improved. However, this improvement is at the cost of more cumbersome expressions, and given how precise the results are already here (cf. figures 4 and 5) it seems to us that this is not necessary in practice here. It could however turn out to be useful in more complicated systems (with magnetic field, non-static equilibrium, non-self-gravitating, etc.).) As a general rule, we find that, in order to get excellent match between numerical and analytic results, we need to increase the order of the expansion in $\bar{k}_{y}$ as we lower the value of $\unicode[STIX]{x1D70C}_{b}$ . Indeed, for large values of $\unicode[STIX]{x1D70C}_{b}$ , typically larger than $0.8$ , the second order can already be excellent, i.e. as good as the fourth order in figure 3. In figures 4 and 5 we show that we obtain excellent agreement for a very wide range of parameters (varying $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D70C}_{b}$ ), but for this to happen on the full unstable range, we need to expand up to the fourth order. Now, the third order is already much more cumbersome than the second order, while it improves the accuracy by only typically several per cent. Therefore quite a lot of effort is required to improve the accuracy of the quadratic formula, and whether or not the gain is worth it depends on the specific problem at hand. As exemplified by figure 3 and further discussed in § 5, the quadratic formula already gives good results on the full unstable range, so in many cases it may constitute the best compromise.

The relevant scales for gravitational instability are $0\leqslant \bar{k}_{y}\leqslant \bar{k}_{crit}$ , by definition of the critical wavenumber $\bar{k}_{crit}$ . As shown in figures 3–5, over that range, our formulæ are very efficient. However, let us note that they are generally not valid much beyond the unstable range. This is no surprise since, in the quartic formula for instance, the function $\bar{\unicode[STIX]{x1D714}}_{0}^{2}$ is a fourth-order polynomial in $\bar{k}_{y}$ , and it has no lower bound for very large $\bar{k}_{y}$ , which is obviously unphysical. This can be seen in figure 4(b): the curve in light grey is part of the $\unicode[STIX]{x1D70C}_{b}=0.4$ solution, but we coloured that branch in light grey because it is outside the relevant domain. In the other curves of figures 3–5 this is not visible simply because it happens beyond the plotted regions.

Figure 5.

Results fixing $\unicode[STIX]{x1D70C}_{b}=0.6$ and varying $\unicode[STIX]{x1D6FE}$ . Solid lines are the analytic solutions at fourth order, and dashed lines are the numerical solutions. In (a) we use rigid BC, while in (b) we use free BC. We indicate with an arrow how the slope at the origin varies as we increase the parameter $\unicode[STIX]{x1D6FE}$ . Here we plot the curves for $\unicode[STIX]{x1D6FE}=0.5$ , 1.5, 3.5 and 5. Formally this evolution can be clearly seen with the expression of $\bar{\unicode[STIX]{x1D714}}_{0,1}^{2}$ in (4.21).

5.1 Isothermal fluid

First, let us consider the isothermal case. For $\unicode[STIX]{x1D6FE}=1$ , $y_{P}=\cos \unicode[STIX]{x1D703}$ and $y_{Q}=-1-\cos \unicode[STIX]{x1D703}\ln \tan \unicode[STIX]{x1D703}/2$ , so that the integrals in (4.21) become straightforward to perform. But in order to get the simplest expressions possible, let us Taylor expand in the high pressure limit, namely $|1-\unicode[STIX]{x1D70C}_{b}|\ll 1$ , i.e. $|\unicode[STIX]{x1D703}_{b}-\unicode[STIX]{x03C0}/2|\ll 1$ . Doing so we obtain the following formulæ for the growth rate as a function of the column density $\unicode[STIX]{x1D707}$ ,

(5.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{rigid BC}\quad \bar{\unicode[STIX]{x1D714}}_{0}^{2}=-\unicode[STIX]{x1D707}\bar{k}_{y}+\left(1+{\textstyle \frac{5}{6}}c_{1}\unicode[STIX]{x1D707}^{2}\right)\bar{k}_{y}^{2}+\mathit{O}(\bar{k}_{y}^{3}), & \displaystyle\end{eqnarray}$$

(5.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{free BC}\quad \bar{\unicode[STIX]{x1D714}}_{0}^{2}=-\unicode[STIX]{x1D707}\bar{k}_{y}+2\unicode[STIX]{x1D707}^{2}(c_{2}-{\textstyle \frac{1}{6}}\unicode[STIX]{x1D707}^{2})\bar{k}_{y}^{2}+\mathit{O}(\bar{k}_{y}^{3}), & \displaystyle\end{eqnarray}$$

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\sqrt{2(1-p_{ext}/p_{c})},\end{eqnarray}$$

with (3.18) and (4.22). We introduced two coefficients, $c_{1}$ and $c_{2}$ which, by the expansion, are equal to unity. We are now going to adjust them in order to improve the efficiency of formulæ (5.1).

For the free BC case (i.e. for the choice of $c_{2}$ ) we proceed as follows. As shown in Goldreich & Lynden-Bell (Reference Goldreich and Lynden-Bell1965), for an incompressible slab the full dispersion relation can be obtained explicitly very simply. Let us denote by $\bar{\unicode[STIX]{x1D714}}_{inc}^{2}$ the maximal growth rate for this case. It is given by

(5.3) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D714}}_{inc}^{2}=-\frac{x}{4}(1-x+\text{e}^{-x})\quad \text{where }x\text{ satisfies }1-2x+(1-x)\text{e}^{-x}=0,\end{eqnarray}$$

that is $x=0.607$ (cf. equation (109) of Goldreich & Lynden-Bell Reference Goldreich and Lynden-Bell1965), and

(5.4) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D714}}_{inc}^{2}=-0.14.\end{eqnarray}$$

This value is key because, as pointed out by previous authors, in the high pressure limit the dynamics depends weakly on $\unicode[STIX]{x1D6FE}$ (this can also be seen in figures 7 and 8 below, where all the curves converge to the same point as $p_{ext}/p_{c}\rightarrow 1$ ). Therefore we should always recover the value (5.4) for the maximum growth rate in the high pressure limit, in particular in the isothermal case considered in this section. We choose the parameter $c_{2}$ in (5.1b ) such that we indeed get $\bar{\unicode[STIX]{x1D714}}_{max}^{2}(p_{ext}/p_{c}\rightarrow 1)=\bar{\unicode[STIX]{x1D714}}_{inc}^{2}$ . This yields

(5.5) $$\begin{eqnarray}c_{2}=-\frac{1}{8\bar{\unicode[STIX]{x1D714}}_{inc}^{2}}.\end{eqnarray}$$

Now for rigid BC (i.e. for the choice of $c_{1}$ ), we note that in the low pressure regime the boundaries go to infinity, therefore in that limit the boundary conditions should not matter any longer. And indeed, in figure 6 the ‘rigid’ and ‘free’ cases agree asymptotically. We can therefore fix the value of $c_{1}$ by equating the maximal wavenumber $\bar{k}_{y,max}$ at $p_{ext}=0$ deduced from (5.1a ) and (5.1b ) respectively. This gives

(5.6) $$\begin{eqnarray}c_{1}=\frac{12c_{2}-7}{5}=-\frac{1}{5}\left(7+\frac{3}{2}\frac{1}{\bar{\unicode[STIX]{x1D714}}_{inc}^{2}}\right).\end{eqnarray}$$

Putting it all together, we finally get that for an isothermal slab:

(5.7a,b ) $$\begin{eqnarray}\displaystyle \text{rigid BC}\quad \bar{k}_{y,max}=\frac{\unicode[STIX]{x1D707}}{2-\left({\displaystyle \frac{7}{3}}+{\displaystyle \frac{1}{2\bar{\unicode[STIX]{x1D714}}_{inc}^{2}}}\right)\unicode[STIX]{x1D707}^{2}}\quad \text{and}\quad \bar{\unicode[STIX]{x1D714}}_{max}^{2}=-\frac{\unicode[STIX]{x1D707}^{2}}{4-\left({\displaystyle \frac{14}{3}}+{\displaystyle \frac{1}{\bar{\unicode[STIX]{x1D714}}_{inc}^{2}}}\right)\unicode[STIX]{x1D707}^{2}}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(5.7c,d ) $$\begin{eqnarray}\displaystyle \text{free BC}\quad \bar{k}_{y,max}=\frac{2}{\unicode[STIX]{x1D707}\left(-{\displaystyle \frac{1}{\bar{\unicode[STIX]{x1D714}}_{inc}^{2}}}-{\displaystyle \frac{4}{3}}\unicode[STIX]{x1D707}^{2}\right)}\quad \text{and}\quad \bar{\unicode[STIX]{x1D714}}_{max}^{2}=\frac{1}{{\displaystyle \frac{1}{\bar{\unicode[STIX]{x1D714}}_{inc}^{2}}}+{\displaystyle \frac{4}{3}}\unicode[STIX]{x1D707}^{2}}. & & \displaystyle\end{eqnarray}$$

These expressions coincide remarkably well with the numerical solution, as shown in figure 6.

Figure 6.

Maximum wavenumber and maximum growth rate as a function of the external pressure $p_{ext}$ , normalised by the central value of the equilibrium pressure $p_{c}=c_{a}^{2}\unicode[STIX]{x1D70C}_{c}$ , for an isothermal slab with free BC in blue and rigid BC in orange. Dots are numerical results, solid lines correspond to (5.7), and the dashed line in (b) corresponds to the incompressible case (5.4).

Figure 7.

Rigid BC: maximum wavenumber and maximum growth rate as a function of the external pressure $p_{ext}$ , normalised by the central value of the equilibrium pressure $p_{c}=c_{a}^{2}\unicode[STIX]{x1D70C}_{c}/\unicode[STIX]{x1D6FE}$ , for various values of $\unicode[STIX]{x1D6FE}$ . Dots are numerical results (corresponding to $\unicode[STIX]{x1D70C}_{b}$ from 0.1 to 0.9 by steps of 0.1) and solid lines correspond to (5.8). With these simple expressions we recover the correct behaviour for a very wide range of $p_{ext}$ and $\unicode[STIX]{x1D6FE}$ .

Figure 8.

Free BC: maximum wavenumber and maximum growth rate as a function of the external pressure $p_{ext}$ , normalised by the central value of the equilibrium pressure $p_{c}=c_{a}^{2}\unicode[STIX]{x1D70C}_{c}/\unicode[STIX]{x1D6FE}$ , for various values of $\unicode[STIX]{x1D6FE}$ . Dots are numerical results (corresponding to $\unicode[STIX]{x1D70C}_{b}$ from 0.1 to 0.9 by steps of 0.1) and solid lines correspond to (5.9). With these simple expressions we recover the correct behaviour for a very wide range of $p_{ext}$ and $\unicode[STIX]{x1D6FE}$ .

5.2 Polytropic fluid

Let us now derive similar expressions, but beyond the isothermal case. To do so, we consider the regime $\unicode[STIX]{x1D6FE}\gg 1$ in formulæ (4.21), in which case the dominant term is the first one. Then, putting $p\equiv p_{ext}/p_{c}$ , we deduce the expressions

(5.8a,b ) $$\begin{eqnarray}\bar{k}_{y,max}=\left(d_{1}+\frac{\ln p}{\unicode[STIX]{x1D6FE}}+\frac{(\ln p)^{2}}{2\unicode[STIX]{x1D6FE}^{2}}\right)\ln \left(\frac{1+\sqrt{1-p}}{\sqrt{p}}\right)\quad \text{and}\quad \bar{\unicode[STIX]{x1D714}}_{max}^{2}=-d_{2}\sqrt{\frac{1-p}{2\unicode[STIX]{x1D6FE}}}~\bar{k}_{y,max}\end{eqnarray}$$

for rigid BC and

(5.9a,b ) $$\begin{eqnarray}\bar{k}_{y,max}=\frac{p^{e_{1}/\unicode[STIX]{x1D6FE}}}{\sqrt{8\unicode[STIX]{x1D6FE}}}\left[\frac{\unicode[STIX]{x1D6FE}}{2}\frac{e_{2}}{\sqrt{1-p}}+e_{3}\ln \left(\frac{1+\sqrt{1-p}}{\sqrt{p}}\right)\right]\quad \text{and}\quad \bar{\unicode[STIX]{x1D714}}_{max}^{2}=-e_{4}\sqrt{\frac{1-p}{2\unicode[STIX]{x1D6FE}}}~\bar{k}_{y,max}\end{eqnarray}$$

for free BC. In the same spirit as above, we introduced coefficients $d_{i}$ and $e_{i}$ in order to extend the domain of applicability of these formulæ. While they are by construction equal to unity, we adjust them to fit the isothermal numerical solution, yielding $(d_{1},d_{2})=(0.956;0.985)$ and $(e_{1},e_{2},e_{3},e_{4})=(0.106;1.174;0.328;0.942)$ . (It may appear a little disappointing to resort, in this very last step, to numerical fits in order to obtain the most efficient formulæ. However, we stress that it is thanks to the foregoing analytic work that we obtained the explicit non-trivial dependences on $\unicode[STIX]{x1D6FE}$ and $p_{ext}$ .) We compare the results to the numerical solution in figures 7 and 8 for rigid and free BC respectively. As we can see, this little tweaking makes the simple expressions (5.8) and (5.9) quantitatively satisfying over a wide range of $\unicode[STIX]{x1D6FE}$ and $p_{ext}$ (or equivalently of the column density since $\unicode[STIX]{x1D707}=\sqrt{(2/\unicode[STIX]{x1D6FE})(1-p_{ext}/p_{c})}$ , even though they were derived in a particular limit. The divergences exhibited by the solid lines in figures 7 and 8 in the $p\rightarrow 0$ limit simply reflect the fact that, in the weak external pressure limit, the expansion in (4.19) would need to be pushed to higher orders in  $\bar{k}_{y}$ .