© The Authors 2024

1. Introduction

The local gravitational instability of rotating fluids is an important phenomenon in several astrophysical systems, ranging from protoplanetary discs to galactic discs. In the case of razor-thin rotating discs, Toomre (1964) derived an analytic instability criterion against small axisymmetric perturbations. However, observations of astrophysical systems have shown that the assumption of infinitesimally thin discs is often unrealistic. For instance, both gaseous and stellar discs in nearby star-forming galaxies have non-negligible thicknesses of ∼0.1 − 1 kpc (e.g. O’Brien et al. 2010; van der Kruit & Freeman 2011; Yim et al. 2014; Mackereth et al. 2017; Marchuk & Sotnikova 2017). Different modifications to Toomre’s 2D analysis that approximately account for the finite thickness of the disc have been proposed in several works (e.g. Toomre 1964; Vandervoort 1970; Romeo 1992; Bertin & Amorisco 2010; Wang et al. 2010; Elmegreen 2011; Griv & Gedalin 2012; Romeo & Falstad 2013; Behrendt et al. 2015).

In many cases, astrophysical discs are not purely self-gravitating, but belong to multi-component systems. A typical case is that of late-type galaxies, in which gaseous discs coexist with stellar discs and with pressure-supported components like bulges and dark-matter halos. When two or more discs coexist, the local gravitational stability of each disc can be influenced by the interplay with the other disc, because any gravitational perturbation affects both discs at the same time (Toomre 1964). The stability of multi-component discs has been widely studied in the literature in the case of infinitesimally thin discs and also in 2D analyses that account approximately for the finite disc thickness (Lin & Shu 1966; Kato 1972; Jog & Solomon 1984a,b; Romeo 1994; Wang & Silk 1994; Elmegreen 1995; Rafikov 2001; Shen & Lou 2003; Romeo & Wiegert 2011).

Nipoti (2023, hereafter N23) presented a three-dimensional (3D) stability analysis of rotating stratified fluids, which generalises previous 3D studies that obtained instability criteria based on more restrictive assumptions (Safronov 1960; Chandrasekhar 1961; Goldreich & Lynden-Bell 1965a,b; Genkin & Safronov 1975; Bertin & Casertano 1982; Mamatsashvili & Rice 2010; Meidt 2022). Here, we extend the work of N23, by studying in 3D the local gravitational instability of two-component discs. We first consider the influence on the 3D instability parameter of a second unresponsive component, and then the case in which both components are responsive. We present the application to two-component discs with isothermal vertical distributions, which can represent, for instance, galactic discs with both stellar and gaseous components. Finally, we relax the assumption of vertical isothermal distribution, by studying self-gravitating discs with vertical polytropic distribution, and thus in general with non-zero vertical temperature or velocity-dispersion gradients. Though the computation of two-component discs with generic polytropic vertical distributions is not more technically difficult than that of of two-component isothermal discs, in this work we present results on the polytropic case only for one-component discs, which is a clean case study that allows us to quantify the effect of non-isothermality, without exploring the larger parameter space of two-component systems.

The paper is organised as follows. Section 2 recalls some relevant equations of previous disc stability analyses. In section 3, we describe the properties of two-component disc models with vertical isothermal distributions. The question of the linear stability in 3D of two-component discs with both components responsive is addressed in Section 4. Discs with vertical polytropic distributions are discussed in Section 5 and Section 6 concludes.

2. Preliminaries

In this work we study axisymmetric rotating fluid disc models, adopting, as it is natural to do, a cylindrical system of coordinates (R, z, ϕ). As we limit our stability analysis to axisymmetric perturbations, none of the quantities considered in this paper depend on ϕ. We considered the properties of a disc at a given radius, R. The disc surface density is

$$\begin{array}{c}\hfill \mathrm{\Sigma }={\int }_{-\infty }^{\infty }\rho (z)\mathrm{d}z,\end{array}$$(1)

where ρ(z) is vertical volume density profile. We have assumed that ρ(z) is such that Σ is finite, and we introduced as a reference vertical scale the half-mass half-height z_half, defined to be such that

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Sigma }}{\int }_{-z_{\mathrm{half}}}^{z_{\mathrm{half}}}\rho (z)\mathrm{d}z=0.5.\end{array}$$(2)

The classical Toomre (1964) 2D criterion for gravitational instability of razor-thin rotating discs is Q < 1, where

$$\begin{array}{c}\hfill Q\equiv \frac{\kappa \sigma }{\pi G\mathrm{\Sigma }}\end{array}$$(3)

is Toomre’s instability parameter at a given radius. Here, κ is the epicycle frequency, defined by

$$\begin{array}{c}\hfill {\kappa }^2\equiv 4{\mathrm{\Omega }}^2+\frac{\mathrm{d}{\mathrm{\Omega }}^2}{\mathrm{d}lnR},\end{array}$$(4)

while Ω is the angular frequency, assumed throughout this paper to depend only on R, and σ is the gas velocity dispersion defined by σ² ≡ P/ρ. The quantity σ can also be interpreted as the ‘isothermal’ sound speed of the fluid, but we prefer to refer to σ generically as the velocity dispersion, because this allows us to have more flexibility in the physical interpretation of the fluid, which can represent not only a standard smooth gas supported by thermal pressure, but also a turbulent gas, a gas composed of discrete gas clouds, or even a stellar distribution (see Section 3.2).

2.1. The three-dimensional gravitational instability criterion

We used as a starting point the results of N23, who has studied in 3D the local gravitational instability of rotating stratified fluids. Among the configurations explored by N23, we considered the simplest: vertically stratified discs (section 3.3 of N23), in which it is assumed that Ω = Ω(R) and that the radial gradients of pressure and density are locally negligible. N23 has shown that a sufficient condition for local gravitational instability at a given (R, z) in such a vertically stratified disc is

$$\begin{array}{c}\hfill Q_{3\mathrm{D}}\equiv \frac{\sqrt{{\kappa }^2+{\nu }^2}+c_{\mathrm{s}}{h}_z^{-1}}{\sqrt{4\pi G\rho }}<1(\mathrm{sufficient}\mathrm{for}\mathrm{instability}),\end{array}$$(5)

where ν, defined by

$$\begin{array}{c}\hfill {\nu }^2\equiv \frac{{\rho }_z^{\prime }{p}_z^{\prime }}{{\rho }^2}=\frac{c_{\mathrm{s}}^2}{\gamma }\frac{{\rho }_z^{\prime }}{\rho }\frac{p_z^{\prime }}{p},\end{array}$$(6)

is a frequency related to vertical pressure and density gradients (${\rho }_z^{\prime }\equiv \partial \rho /\partial z$, $p_z^{\prime }\equiv \partial p/\partial z$), h_z is a measure of the disc thickness, and $c_{\mathrm{s}}=\sqrt{\gamma }\sigma$ is the sound speed¹, with γ ≡ (ρ/p)dp/dρ, in which dp is the pressure variation corresponding to a density variation, dρ, during a transformation. Following N23, throughout the paper we assume that h_z = h_70%, where h_70%, defined by

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Sigma }}{\int }_{-h_{70\%}/2}^{h_{70\%}/2}\rho (z)\mathrm{d}z=0.7,\end{array}$$(7)

is the height of an infinitesimal-width strip centred on the midplane containing 70% of the mass per unit surface. The scale height, h_z, appearing in the definition of Q_3D derives from the approximation used by N23 for the linearised Poisson equation for radial perturbations. In Appendix A, we assess quantitatively the validity of this approximation when h_z = h_70%, by considering one-component discs. In the following, we also adopt h_z = h_70% for two-component discs, assuming that it remains a sufficiently good approximation in the presence of a second component.

2.2. Discs with a self-gravitating isothermal vertical distribution

Before discussing more generic models, it is useful to recall the properties of a disc in which the vertical density distribution is that of the self-gravitating isothermal (SGI) slab (Spitzer 1942),

$$\begin{array}{c}\hfill \rho (z)={\rho }_0{\mathrm{sech}}^2\stackrel{\sim }{z},\end{array}$$(8)

where ρ₀ = ρ(0) is the density in the midplane and $\stackrel{\sim }{z}\equiv z/b$, with $b\equiv \sigma /\sqrt{2\pi G{\rho }_0}$ a scale height, where σ is the z-independent velocity dispersion. The gravitational potential is

$$\begin{array}{c}\hfill \mathrm{\Phi }(z)=-{\sigma }^{2}ln\left[\frac{\rho (z)}{{\rho }_0}\right]\end{array}$$(9)

and the surface density is Σ = 2bρ₀, so ρ₀ = πGΣ²/(2σ²). Using p(z) = σ²ρ(z), Eq. (6), and Eq. (8), we obtain for the SGI disc model

$$\begin{array}{c}\hfill {\nu }^2=8\pi G{\rho }_0{tanh}^2\stackrel{\sim }{z}.\end{array}$$(10)

So, for a disc with an angular velocity Ω = Ω(R) and SGI vertical distribution, at any given R the sufficient condition for instability (5) becomes (N23)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Q_{3\mathrm{D}}=\sqrt{\frac{Q^2}{2}{cosh}^2\stackrel{\sim }{z}+2{sinh}^2\stackrel{\sim }{z}}+\sqrt{\frac{\gamma }{2}}\frac{b}{h_z}cosh\stackrel{\sim }{z}<1\\ \hfill (\mathrm{sufficient}\mathrm{for}\mathrm{instability}).\end{array}\end{array}$$(11)

3. Two-component discs with isothermal vertical distributions

We consider here two-component discs, in which each component has an isothermal vertical density distribution.

3.1. Equations

To describe a two-component disc, we use the index i (for i = 1, 2) to label the generic ith component. At any given R the vertical equilibrium is described by a two-component isothermal slab (see Bertin & Pegoraro 2022) composed of two components with density distributions ρ₁(z) and ρ₂(z) in equilibrium in the total gravitational potential,

$$\begin{array}{c}\hfill \mathrm{\Phi }(z)={\mathrm{\Phi }}_1(z)+{\mathrm{\Phi }}_2(z),\end{array}$$(12)

where Φ_i, satisfying

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\mathrm{\Phi }}_i}{\mathrm{d}{z}^2}=4\pi G{\rho }_i,\end{array}$$(13)

is the gravitational potential generated by ρ_i. Each component is in vertical hydrostatic equilibrium in the total potential Φ so

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}_i}{\mathrm{d}z}=-{\rho }_i\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z},\end{array}$$(14)

where $p_i(z)={\sigma }_i{}^2{\rho }_i(z)$, σ_i, and ρ_i(z) are, respectively, the pressure, velocity dispersion, and density of the ith component. Given that σ_i is independent of z, Eq. (14) can be written as

$$\begin{array}{c}\hfill {{\sigma }_i}^2\frac{\mathrm{d}ln{\rho }_i}{\mathrm{d}z}=-\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z},\end{array}$$(15)

which has the solution

$$\begin{array}{c}\hfill {\rho }_i(z)={\rho }_{0,i}{e}^{-\mathrm{\Phi }/{{\sigma }_i}^2},\end{array}$$(16)

where ρ_0, i is the midplane density of the ith component.

Combining Eqs. (12), (13), and (16), we get the ordinary differential equation

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathrm{\Phi }}{\mathrm{d}{z}^2}=4\pi G({\rho }_{0,1}{e}^{-\mathrm{\Phi }/{{\sigma }_1}^2}+{\rho }_{0,2}{e}^{-\mathrm{\Phi }/{{\sigma }_2}^2}),\end{array}$$(17)

which, for given ρ_0, 1, ρ_0, 2, σ₁, and σ₂, can be solved to obtain Φ(z), and then ρ₁(z) and ρ₂(z) from Eq. (16). Eq. (17) can be written in dimensionless form as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{{\mathrm{d}}^2\stackrel{\sim }{\mathrm{\Phi }}}{\mathrm{d}{\stackrel{\sim }{z}}^2}=2(e^{-\stackrel{\sim }{\mathrm{\Phi }}}+\xi e^{-\stackrel{\sim }{\mathrm{\Phi }}/\mu }),\end{array}\end{array}$$(18)

where $\stackrel{\sim }{\mathrm{\Phi }}\equiv \mathrm{\Phi }/{{\sigma }_1}^2$, $\mu \equiv {\sigma }_2{}^2/{\sigma }_1{}^2$, ξ ≡ ρ_0, 2/ρ_0, 1, and $\stackrel{\sim }{z}\equiv z/b_1$, with $b_1\equiv {\sigma }_1/\sqrt{2\pi G{\rho }_{0,1}}$ a scale height. The two-component isothermal slab models defined above have only two free parameters: ξ (the central density ratio of the two components) and μ (their velocity dispersion ratio squared). Assuming that the two components share the same κ, we define the Q parameter of the ith component as

$$\begin{array}{c}\hfill Q_i=\frac{\kappa {\sigma }_i}{\pi G{\mathrm{\Sigma }}_i},\end{array}$$(19)

where Σ_i is the surface density of the ith component.

3.2. Vertical density profiles

We present here an example of an equilibrium two-component isothermal model obtained by solving numerically Eq. (18) for given values of the parameters μ and ξ, with boundary conditions $\stackrel{\sim }{\mathrm{\Phi }}=0$ and $\mathrm{d}\stackrel{\sim }{\mathrm{\Phi }}/\partial \stackrel{\sim }{z}=0$ at $\stackrel{\sim }{z}=0$. In particular, the numerical solutions were computed using the Mathematica² routine NDsolve.

For the sake of definiteness, we now refer to the specific case of a galaxy with a gaseous disc and a stellar disc, and thus we use the equations of Section 3.1, replacing the subscript ‘1’ with ‘gas’ and the subscript ‘2’ with ‘⋆’. In this context, it is natural to assume that the effective gas pressure is $p_{\text{gas}}={\rho }_{\text{gas}}{\sigma }_{\text{gas}}{}^2$, with ${{\sigma }_{\mathrm{gas}}}^2=k_{\mathrm{B}}T/\overline{m}+{{\sigma }_{\mathrm{turb}}}^2$, where T is the gas temperature, k_B is the Boltzmann constant, $\overline{m}$ is the mean molecular or atomic mass, and σ_turb is the turbulent velocity dispersion. This is motivated by the fact that there is observational evidence that galactic gaseous discs (either atomic or molecular) are turbulent, and thus σ_turb contributes substantially to σ_gas (e.g. section 4.2.2 of Cimatti et al. 2019; Bacchini et al. 2020). Alternatively, the gas disc can be modelled as a discrete distribution of gas clouds (e.g. Jeffreson et al. 2022), in which case σ_gas is just the cloud-cloud velocity dispersion. The vertical structure of a stellar disc is described by the same equations used for a gas, but with σ_⋆ given by the vertical stellar velocity dispersion. Given that in a galaxy the stellar disc is in general dynamically hotter than the gaseous disc (at least for sufficiently old stellar populations; van der Kruit & Freeman 2011), without loss of generality, we limit ourselves to exploring cases with $\mu ={\sigma }_{\star }{}^2/{\sigma }_{\text{gas}}{}^2\ge 1$.

An example of the numerically computed vertical density distributions of a two-component isothermal model with ξ = 0.3 and μ = 9 is given in the upper panel of Fig. 1: for this model, Σ_⋆/Σ_gas ≃ 1.27. For comparison, in the same diagram we also plot the vertical gas density distributions of a gaseous SGI slab with the same σ_gas and Σ_gas as the gaseous component of the the two-component system and of a stellar SGI slab with the same σ_⋆ and Σ_⋆ as the stellar component of the two-component system. The stellar component, having a higher velocity dispersion, is more vertically extended than the gaseous component. Comparing the solid and dashed curves in the upper panel of Fig. 1, it is apparent how, for given surface density and velocity dispersion, the volume density distribution depends on the presence or absence of a second component. Due to the deeper gravitational potential, in the presence of the stellar disc the gaseous disc is more concentrated (i.e. has a higher central volume density) than in the self-gravitating case; vice versa, in the presence of the gas disc, the stellar disc is more concentrated than in the self-gravitating case.

Fig. 1. Upper panel: Vertical density profiles (solid curves) of gas and stars of a two-component isothermal disc with μ = 9 and ξ = 0.3. The dashed curves are the profiles of the corresponding SGI models with the same velocity dispersion and surface density. Lower panel: Vertical Q3D profiles of the same models as in the upper panel, assuming Qgas = 0.6 (and thus Q⋆ = 1.4). The horizontal dashed line indicates the instability threshold, Q3D = 1. Here, ρ0, gas and bgas are, respectively, the midplane density and the scale height of the gaseous component of the two-component disc.

For a given vertical gravitational field, the scale height increases for increasing σ (see also Bacchini et al. 2019a). Computing numerically the ratios, z_half, ⋆/z_half, gas (see Eq. 2) and h_z, ⋆/h_z, gas (see Eq. 7), we find that they are reasonably well described by the function μ^5/8. This is illustrated in Fig. 2, where the ratio h_z, ⋆/h_z, gas is compared with μ^5/8 for ξ = 0.5, ξ = 1, and ξ = 2. A very similar behaviour is found for the scale-height ratio, if one defines the scale height of the ith component as Σ_i/(2ρ_0, i), as was done in Romeo (1992), who finds a slope 3/5 and Bertin & Pegoraro (2022; see their figure 2). Fig. 3 shows how the ratios, h_z, ⋆/z_half, ⋆ and h_z, gas/z_half, gas, depend only weakly on μ and ξ in the explored range of values of these parameters. For μ close to unity, h_z, ⋆/z_half, ⋆ ≈ h_z, gas/z_half, gas ≈ 3.16 independently of ξ. For large values of μ, h_z, ⋆/z_half, ⋆ ≈ 3.22 and h_z, gas/z_half, gas ≈ 3.12, with a weak dependence on ξ.

Fig. 2. Upper panel: Stars-to-gas thickness ratio, hz, ⋆/hz, gas, as a function of the squared velocity-dispersion ratio, μ, for two-component isothermal slabs with different values of the midplane density ratio ξ (hz = h70%, defined by Eq. 7). Overplotted in grey is the power law, μ5/8. Lower panel: Percent residuals between the points and the power law of the upper panel.

Fig. 3. Ratios hz, ⋆/zhalf, ⋆ (filled symbols) and hz, gas/zhalf, gas (empty symbols) as functions of the squared velocity-dispersion ratio, μ, for two-component isothermal slabs with different values of the midplane density ratio, ξ. hz = h70% and zhalf are defined by Eq. (7) and Eq. (2), respectively.

3.3. Stability of each component, assuming that the other component is unresponsive

We study here the local gravitational instability of the two-component discs introduced in Section 3.2, using the 3D instability criterion described in Section 2.1. This is possible for each of the two discs, under the simplifying assumption that the other disc is unresponsive; that is, it acts just as a fixed external gravitational potential and does not react dynamically to the perturbations (we discuss the effects of a responsive disc in Section 4). Though simplified, this approach is useful because, given that the back-reaction of the other component favours the instability (Toomre 1964), Q_3D < 1 is also a sufficient condition for instability in the presence of a responsive second component.

As in Section 3.2, we interpret one of the two components as a stellar disc. Stellar discs are collisionless systems so, strictly speaking, their gravitational stability or instability cannot be assessed with the same methods used for fluid systems. However, it has been shown that treating the stellar disc as a fluid is a good approximation, at least in 2D stability analyses (Bertin & Romeo 1988; Rafikov 2001; Romeo & Falstad 2013). To apply the criterion (5) to the stellar component, we thus assume that the stellar component has an isotropic velocity dispersion tensor and we model it as a fluid (see Section 4.3 for a discussion). Here, we assume that the gas and the stars undergo ‘isothermal’ transformations, that is such that $\text{d}{p}_{\text{gas}}={\sigma }_{\text{gas}}{}^2\text{d}{\rho }_{\text{gas}}$ and $\text{d}{p}_{\star }={\sigma }_{\star }{}^2\text{d}{\rho }_{\star }$, which is in practice equivalent to using γ = 1 in the equations introduced in Section 2. In principle, we could consider full axisymmetric disc models in which Σ_gas, σ_gas, κ, and thus Q_gas = κσ_gas/(πGΣ_gas), as well as the corresponding stellar quantities, vary as functions of R (see Bacchini et al. 2024). Instead, as is done in section 5 of N23, we exploited the mapping between Q_gas and R to study the stability at a given R, as a function of z, simply by fixing a value of Q_gas. Once Q_gas is fixed, Q_⋆ is determined by the condition $Q_{\star }/Q_{\mathrm{gas}}=\sqrt{\mu }{\mathrm{\Sigma }}_{\mathrm{gas}}/{\mathrm{\Sigma }}_{\star }$.

For the gaseous disc in the presence of stars, Q_3D can be computed numerically once ρ_gas(z) and its derivative are computed, given that $p_{\text{gas}}={\sigma }_{\text{gas}}{}^2{\rho }_{\text{gas}}$ (with σ_gas independent of z) and that Q_3D (Eq. 5) can be rewritten for the gas component as

$$\begin{array}{c}\hfill Q_{3\mathrm{D},\mathrm{gas}}(z)=\frac{1}{\sqrt{2{\stackrel{\sim }{\rho }}_{\mathrm{gas}}(z)}}[\sqrt{Q_{\mathrm{gas}}^2{{\stackrel{\sim }{\mathrm{\Sigma }}}_{\mathrm{gas}}}^2+\frac{1}{{{\stackrel{\sim }{\rho }}_{\mathrm{gas}}}^2(z)}{\left(\frac{\mathrm{d}{\stackrel{\sim }{\rho }}_{\mathrm{gas}}}{\mathrm{d}\stackrel{\sim }{z}}\right)}^2}+\frac{1}{{\stackrel{\sim }{h}}_{z,\mathrm{gas}}}].\end{array}$$(20)

Similarly, for the stellar component

$$\begin{array}{c}\hfill Q_{3\mathrm{D},\star }(z)=\frac{1}{\sqrt{2{\stackrel{\sim }{\rho }}_{\star }(z)}}[\sqrt{\frac{Q_{\star }^2{{\stackrel{\sim }{\mathrm{\Sigma }}}_{\star }}^2}{\mu }+\frac{\mu }{{{\stackrel{\sim }{\rho }}_{\star }}^2(z)}{\left(\frac{\mathrm{d}{\stackrel{\sim }{\rho }}_{\star }}{\mathrm{d}\stackrel{\sim }{z}}\right)}^2}+\frac{\sqrt{\mu }}{{\stackrel{\sim }{h}}_{z,\star }}],\end{array}$$(21)

which can be obtained if ρ_⋆(z) and its derivative are known numerically. In the above equations, where we have used γ = 1, ${\stackrel{\sim }{\rho }}_{\mathrm{gas}}={\rho }_{\mathrm{gas}}/{\rho }_{0,\mathrm{gas}}$, ${\stackrel{\sim }{\rho }}_{\star }={\rho }_{\star }/{\rho }_{0,\mathrm{gas}}$, ${\stackrel{\sim }{h}}_{z,\mathrm{gas}}=h_{z,\mathrm{gas}}/b_{\mathrm{gas}}$, ${\stackrel{\sim }{h}}_{z,\star }=h_{z,\star }/b_{\mathrm{gas}}$, ${\stackrel{\sim }{\mathrm{\Sigma }}}_{\mathrm{gas}}\equiv {\mathrm{\Sigma }}_{\mathrm{gas}}/(2b_{\mathrm{gas}}{\rho }_{0,\mathrm{gas}})$, and ${\stackrel{\sim }{\mathrm{\Sigma }}}_{\star }\equiv {\mathrm{\Sigma }}_{\star }/(2b_{\mathrm{gas}}{\rho }_{0,\mathrm{gas}})$, where $b_{\mathrm{gas}}\equiv {\sigma }_{\mathrm{gas}}/\sqrt{2\pi G{\rho }_{0,\mathrm{gas}}}$.

Profiles of Q_3D for the gaseous and stellar components are shown as solid curves in the lower panel of Fig. 1, assuming that Q_gas = 0.6, for the same two-component disc models as in the upper panel of Fig. 1. For comparison, in the lower panel of Fig. 1 we plot as dashed curves Q_3D (Eq. 11) for the corresponding SGI gas and stellar discs with the same surface density, velocity dispersion, and Q. This example illustrates that the presence of a second component can affect Q_3D in different ways: in this case, close to the midplane the gas disc has higher Q_3D, gas in the presence than in the absence of the stellar disc (at fixed Q_gas, Σ_gas, and σ_gas), while the stellar disc has lower Q_3D, ⋆ in the presence than in the absence of the stellar disc (at fixed Q_⋆, Σ_⋆, and σ_⋆). This occurs because Q_3D depends not only on the volume density (which, close to the midplane, is invariably higher in the presence of a second component, for given Σ and σ), but also on the shape of the vertical density profile, through (dρ_i/dz)/ρ_i and h_z, i (see Eqs. 20 and 21).

A systematic study of the family of two-component isothermal disc models with Q_gas = 0.6 for a range of values of μ and ξ is displayed as an illustrative example in Fig. 4, which shows that for the considered values of Q_gas there is always instability in the midplane, in the sense that either Q_3D, gas(0) < 1 or Q_3D, ⋆(0) < 1, and that the behaviour of the system depends weakly on μ. In particular, from the upper panel of Fig. 4, plotting Q_3D, gas(0) and Q_3D, ⋆(0) as functions of ξ, we see that when the gas disc is dominant (low ξ) the instability is driven by the gas disc (Q_3D, gas(0) < 1), while when the stellar disc is dominant (high ξ) the instability is driven by the stars (Q_3D, ⋆(0) < 1). The lower panel of Fig. 4 shows, for each component in units of its z_half, the (above or below midplane) extent z_inst of the unstable region, defined by the condition Q_3D(z) < 1 for |z|< z_inst and Q_3D(z) > 1 for |z|> z_inst. When one of the components is dominant (ξ ≪ 1 or ξ ≫ 1), we have 0.5 ≲ z_inst/z_half ≲ 0.7 for Q_gas = 0.6. When the two components are comparable (ξ ≈ 1), z_inst/z_half ≈ 0.3 for Q_gas = 0.6. Fig. 4 suggests that the Q_3D-based stability properties of two-component discs are essentially independent of μ. At fixed Q_gas, the midplane density ratio, ξ, determines which component drives the instability, but only weakly influences the overall instability properties of the system: the minimum between Q_3D, gas(0) and Q_3D, ⋆(0) varies only by about a factor of two when ξ spans two orders of magnitude (see upper panel of Fig. 4). We recall, however, that the analysis based on Q_3D does not capture the full complexity of the problem of the instability of two-component discs, in which the coupling between stellar and gaseous perturbations plays an important role (e.g. Bertin & Romeo 1988; Romeo & Falstad 2013). The effect of this coupling in 3D two-component discs is discussed in the next section.

Fig. 4. Upper panel: Midplane values of Q3D, ⋆ (filled symbols) and Q3D, gas (empty symbols) as functions of the midplane star-to-gas density ratio, ξ, for different values of the star-to-gas squared velocity dispersion ratio, μ. The horizontal dashed line indicates the instability threshold, Q3D = 1. Lower panel: Extent of the unstable strip above the midplane for the stellar (zinst, ⋆, normalised to zhalf, ⋆; filled symbols) and gaseous (zinst, gas, normalised to zhalf, gas; empty symbols) components as functions of ξ for different values of μ. In both panels, we assume that Qgas = 0.6.

4. Stability of two-component discs when both components are responsive

In Section 3.3, we studied the stability of a 3D gaseous disc in the presence of a stellar disc, but neglecting the back-reaction of the stellar disc, and vice versa. Here, we attempt to address the more realistic but more complicated question of the instability of a 3D two-component disc when both components are responsive. As was done in Section 3.3, we assume for simplicity that both components can be treated as fluids (even in the case in which one of the components is a stellar disc). The limitations of this assumption are discussed in Section 4.3.

4.1. Linear perturbation analysis

The governing equations are, for each component, the same as in section 2.1 of N23, but with the gravitational potential, Φ, given by the sum of the gravitational potentials of the two components (Eq. 12):

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial {\rho }_i}{\partial t}+\frac{1}{R}\frac{\partial (R{\rho }_i{u}_{R,i})}{\partial R}+\frac{\partial ({\rho }_i{u}_{z,i})}{\partial z}=0,\hfill \\ \hfill & \frac{\partial u_{R,i}}{\partial t}+u_{R,i}\frac{\partial u_{R,i}}{\partial R}+u_{z,i}\frac{\partial u_{R,i}}{\partial z}-\frac{u_{\varphi ,i}^2}{R}=-\frac{1}{{\rho }_i}\frac{\partial p_i}{\partial R}-\frac{\partial \mathrm{\Phi }}{\partial R}-\frac{\partial {\mathrm{\Phi }}_{\mathrm{ext}}}{\partial R},\hfill \\ \hfill & \frac{\partial u_{\varphi ,i}}{\partial t}+u_{R,i}\frac{\partial u_{\varphi ,i}}{\partial R}+u_{z,i}\frac{\partial u_{\varphi ,i}}{\partial z}+\frac{u_{R,i}{u}_{\varphi ,i}}{R}=0,\hfill \\ \hfill & \frac{\partial u_{z,i}}{\partial t}+u_{R,i}\frac{\partial u_{z,i}}{\partial R}+u_{z,i}\frac{\partial u_{z,i}}{\partial z}=-\frac{1}{{\rho }_i}\frac{\partial p_i}{\partial z}-\frac{\partial \mathrm{\Phi }}{\partial z}-\frac{\partial {\mathrm{\Phi }}_{\mathrm{ext}}}{\partial z},\hfill \\ \hfill & \frac{p_i}{\gamma -1}(\frac{\partial }{\partial t}+{\mathit{u}}_i·\mathrm{\nabla })ln(p_i{{\rho }_i}^{-\gamma })=0,\hfill \\ \hfill & \frac{1}{R}\frac{\partial }{\partial R}\left(R\frac{\partial {\mathrm{\Phi }}_i}{\partial R}\right)+\frac{{\partial }^2{\mathrm{\Phi }}_i}{\partial z^2}=4\pi G{\rho }_i,\hfill \end{array}\end{array}$$(22)

with i = 1, 2, where u_i = u_R,i,u_ϕ,i,u_z,i) is the velocity of the ith component, Φ = Φ₁ + Φ₂, and Φ_ext is any additional gravitational potential (for instance that of the dark-matter halo) assumed to be fixed³.

Here, for each component, we make the same assumptions described in Section 2.1. The unperturbed system is a stationary rotating (u_ϕ, i ≠ 0) solution of Eqs. (22) with no meridional motions (u_R, i = u_z, i = 0), locally negligible radial pressure and density gradients, and u_ϕ, 1 = u_ϕ, 2 = ΩR, with Ω = Ω(R). We note that the fact that the radial pressure gradient is negligible in the unperturbed disc (i.e. |dp_i/dR|/ρ_i ≪ Ω²R) implies that the epicyclic approximation condition, σ_i ≪ κR (e.g. Bertin 2014), is satisfied. This can be seen by considering the following ordering:

$$\begin{array}{c}\hfill \frac{{{\sigma }_i}^2}{R}\sim {{\sigma }_i}^2\left|\frac{\mathrm{d}ln{p}_i}{\mathrm{d}R}\right|\sim \frac{1}{{\rho }_i}\left|\frac{\mathrm{d}{p}_i}{\mathrm{d}R}\right|\ll {\mathrm{\Omega }}^{2}R\sim {\kappa }^{2}R,\end{array}$$(23)

where we have used |dlnp_i/dR| ∼ 1/R and κ ∼ Ω.

As in N23, we perturbed the system (22), writing a generic quantity q = q(R, z, t) (such as ρ_i, p_i, Φ_i, or any component of u_i) as q = q_unp + δq, where the (time-independent) quantity q_unp describes the stationary unperturbed fluid and the (time-dependent) quantity δq describes the Eulerian perturbation. From now on, without risk of ambiguity, we indicate any unperturbed quantity q_unp simply as q. Limiting ourselves to a linear stability analysis, we consider small (|δq/q|≪1) perturbations. Given that the radial perturbations tend to be more unstable than the vertical ones (Goldreich & Lynden-Bell 1965a, N23), we considered purely radial disturbances with spatial and temporal dependence δq ∝ exp[i(k_RR − ωt)], where ω is the frequency and k_R is the radial component of the wave vector, which we assume to be such that |k_R|R ≫ 1, as it is standard in the short-wavelength approximation. Moreover, given that in the one-component case Q_3D is lowest in the midplane, for simplicity we focused only on the midplane, where we can adopt⁴ ${\rho }_{z,i}^{\prime }=0$ and $p_{z,i}^{\prime }=0$. Under these assumptions, the linearised perturbed equations are

$$\begin{array}{c}\hfill \begin{array}{cc}& -\mathrm{i}\omega \delta {\rho }_i+\mathrm{i}{k}_R{\rho }_i\delta u_{R,i}=0,\hfill \\ \hfill & -\mathrm{i}\omega \delta u_{R,i}-2\mathrm{\Omega }\delta u_{\varphi ,i}=-\mathrm{i}\frac{k_R}{{\rho }_i}\delta p_i-\mathrm{i}{k}_R\delta {\mathrm{\Phi }}_1-\mathrm{i}{k}_R\delta {\mathrm{\Phi }}_2,\hfill \\ \hfill & -\mathrm{i}\omega \delta u_{\varphi ,i}+\frac{\mathrm{d}(\mathrm{\Omega }R)}{\mathrm{d}R}\delta u_{R,i}+\mathrm{\Omega }\delta u_{R,i}=0,\hfill \\ \hfill & -\mathrm{i}\omega \delta u_{z,i}=0,\hfill \\ \hfill & -\mathrm{i}\omega \frac{\delta p_i}{p_i}+\mathrm{i}\gamma \omega \frac{\delta {\rho }_i}{{\rho }_i}=0,\hfill \\ \hfill & -({k_R}^2+{h_{z,i}}^{-2})\delta {\mathrm{\Phi }}_i=4\pi G\delta {\rho }_i,\hfill \end{array}\end{array}$$(24)

with i = 1, 2, where the last equation is the perturbed Poisson equation in the form of equation 14 of N23. As was done in the previous sections, here we assume that h_z = h_70% (Eq. 7; see Appendix A and Section 2.1 for a discussion).

The system (24) leads to the biquadratic dispersion relation

$$\begin{array}{c}\hfill {\omega }^4-({\alpha }_1+{\alpha }_2){\omega }^2+{\alpha }_1{\alpha }_2-{\beta }_1{\beta }_2=0,\end{array}$$(25)

where we have defined

$$\begin{array}{c}\hfill {\alpha }_1\equiv B+s-A\frac{s}{s+E_1},\end{array}$$(26)

$$\begin{array}{c}\hfill {\alpha }_2\equiv B+\mu s-\xi A\frac{\mu s}{\mu s+E_2},\end{array}$$(27)

$$\begin{array}{c}\hfill {\beta }_1=A\frac{s}{s+E_1},\end{array}$$(28)

and

$$\begin{array}{c}\hfill {\beta }_2=\xi A\frac{\mu s}{\mu s+E_2},\end{array}$$(29)

with B ≡ κ², $s\equiv \gamma {\sigma }_1{}^2{k}_R{}^2$, A ≡ 4πGρ₁, $E_1\equiv \gamma {{\sigma }_1}^2{h_{z,1}}^{-2}$, $E_2\equiv \gamma {{\sigma }_2}^2{h_{z,2}}^{-2}$, $\mu \equiv {\sigma }_2{}^2/{\sigma }_1{}^2$, and ξ = ρ_0, 2/ρ_0, 1.

The dispersion relation (25) is formally identical to that found by Jog & Solomon (1984b) in their 2D analysis, but with different definitions and physical meaning of the coefficients, so our stability analysis essentially follows that of section II of Jog & Solomon (1984b, see also Hoffmann & Romeo 2012 for further analysis of the same dispersion relation). Eq. (25) has the roots

$$\begin{array}{c}\hfill {\omega }_{-,+}^2=\frac{1}{2}[({\alpha }_1+{\alpha }_2)\pm \sqrt{{({\alpha }_1+{\alpha }_2)}^2-4({\alpha }_1{\alpha }_2-{\beta }_1{\beta }_2)}].\end{array}$$(30)

It is straightforward to show that the argument of the square root is always positive, so the solutions are always real. Given that ${\omega }_{-}^2\le {\omega }_{+}^2$, we can focus on the root, ${\omega }_{-}^2$, to study the stability. When either α₁ < 0 or α₂ < 0, ${\omega }_{-}^2$ < 0, and thus we have instability. We note that the condition α₁ < 0 is equivalent to

$$\begin{array}{c}\hfill F_1(k_R)\equiv \frac{\mathit{As}}{(s+B)(s+E_1)}>1,\end{array}$$(31)

which, following the analysis reported in appendix B3 of N23 leads to the sufficient condition for instability

$$\begin{array}{c}\hfill \frac{\sqrt{B}+\sqrt{E_1}}{\sqrt{A}}<1,\end{array}$$(32)

which is just Q_3D, 1 < 1 in the special case, ${\rho }_{z,1}^{\prime }=0$, considered in this section. It is straightforward to show that α₂ < 0, that is

$$\begin{array}{c}\hfill F_2(k_R)=\frac{\xi A\mu s}{(s+B)(\mu s+E_2)}>1,\end{array}$$(33)

leads to

$$\begin{array}{c}\hfill \frac{\sqrt{B}+\sqrt{E_2}}{\sqrt{A}}<1,\end{array}$$(34)

that is Q_3D, 2 < 1 when ${\rho }_{z,2}^{\prime }=0$. Thus, as was expected, when one of the components has Q_3D < 1 in the midplane (i.e. is unstable, neglecting the back-reaction of the other component), the two-component system turns out also to be unstable when the back-reaction is accounted for.

When both α₁ > 0 (i.e. F₁ < 1) and α₂ > 0 (i.e. F₂ < 1), the condition to have ${\omega }_{-}^2$ < 0 (and thus instability) is

$$\begin{array}{c}\hfill {\alpha }_1{\alpha }_2-{\beta }_1{\beta }_2<0,\end{array}$$(35)

that is

$$\begin{array}{c}\hfill F_{1+2}(k_R)=F_1(k_R)+F_2(k_R)>1,\end{array}$$(36)

where the destabilising effect of the second responsive component is apparent. This is a manifestation of the fact that, because of the mutual back-reactions of the two components, a two-component system can be gravitationally unstable even when each component would be stable if taken individually.

When Q_3D, 1 > 1 and Q_3D, 2 > 1, and thus F₁ < 1 and F₂ < 1 for all k_R, F_1 + 2 > 1 is a sufficient condition for instability.

4.2. Application to two-component discs with stellar and gaseous components

In order to apply quantitatively the results of Section 4.1, we now specialize to the case of a two-component disc with gas and stellar components. Relabelling the component indices 1 and 2 as ‘gas’ and ‘⋆’, respectively, as was done in Sections 3.2 and 3.3, it follows from the calculations of Section 4.1 that when either Q_3D, gas < 1 or Q_3D, ⋆ < 1 in the midplane, we have instability. We thus focus on the case in which in the midplane Q_3D, gas > 1 and Q_3D, ⋆ > 1, so that for all k_R we have α_gas > 0 and α_⋆ > 0, which, using the definitions (31) and (33) can be written, respectively, as

$$\begin{array}{c}\hfill F_{\mathrm{gas}}(k_R)=\frac{2{{\stackrel{\sim }{k}}_R}^2}{(\gamma {{\stackrel{\sim }{k}}_R}^2+Q_{\mathrm{gas}}^2{{\stackrel{\sim }{\mathrm{\Sigma }}}_{\mathrm{gas}}}^2)({{\stackrel{\sim }{k}}_R}^2+{{\stackrel{\sim }{h}}_{z,\mathrm{gas}}}^{-2})}<1\end{array}$$(37)

and

$$\begin{array}{c}\hfill F_{\star }(k_R)=\frac{2\xi {{\stackrel{\sim }{k}}_R}^2}{(\mu \gamma {{\stackrel{\sim }{k}}_R}^2+Q_{\star }^2{{\stackrel{\sim }{\mathrm{\Sigma }}}_{\star }}^2/\mu )({{\stackrel{\sim }{k}}_R}^2+{{\stackrel{\sim }{h}}_{z,\star }}^{-2})}<1,\end{array}$$(38)

where ${\stackrel{\sim }{k}}_R\equiv k_R{b}_{\mathrm{gas}}$. The condition for a disturbance with a radial wave number k_R to be unstable is Eq. (36), which in this case reads

$$\begin{array}{c}\hfill F_{\mathrm{gas}+\star }(k_R)\equiv F_{\mathrm{gas}}(k_R)+F_{\star }(k_R)>1.\end{array}$$(39)

We then consider two-component isothermal discs in which the vertical distributions are computed numerically, as it is described in Section 3.2. In the example shown in Figs. 5 and 6, where we assume that γ = 1 as in Section 3.3, the values of the parameters (ξ = 0.75, μ = 8, Q_gas = 1.1 and Q_⋆ = 1.2) are such that the two-component isothermal disc would not be considered unstable according to the analysis of Section 3.3, but it turns out to be unstable when the back-reaction of one component onto the other is taken into account. Similar to the model shown in Fig. 1, the model shown in Fig. 5, for which Σ_⋆/Σ_gas ≃ 2.26, has higher gas than stellar density close to the midplane (see upper panel of Fig. 5), but κ is assumed to be such that Q_3D, gas > 1 and Q_3D, ⋆ > 1 (see lower panel of Fig. 5). Nevertheless, the system is gravitationally unstable, because there is a range of values of the wave number such that F_{gas + ⋆} > 1 (Fig. 6).

Fig. 5. Same as Fig. 1, but for a two-component isothermal disc model with ξ = 0.75, μ = 6, Qgas = 1.1, and Q⋆ = 1.2.

Fig. 6. Functions Fgas(kR), F⋆(kR), and Fgas + ⋆(kR) (see Eqs. 37, 38 and 39) for the same model as in Fig. 5. The unstable perturbations are those with Fgas + ⋆ > 1.

We stress that the criterion (39) can be applied to observed discs, because all the coefficients of the functions F_gas and F_⋆ can be estimated from observable quantities. In particular, to fully determine the function F_{gas + ⋆}(k_R) at given R one needs estimates of Σ_gas, Σ_⋆, ρ_0, gas, ρ_0, ⋆, h_z, gas, h_z, ⋆, σ_gas, σ_⋆, and κ, which can all be inferred from observational data (see Bacchini et al. 2024). Then, to draw conclusions on the local gravitational instability of the system at given R in the disc midplane, it is sufficient to evaluate F_{gas + ⋆} over a range of values of |k_R|, with a lower limit between 1/R and 1/h_z (see Section 4.1) and an upper limit ≫1/h_z (because the unstable disturbances have |k_R|h_z of the order of unity; Goldreich & Lynden-Bell 1965a; N23).

4.3. Limitations of the application to stellar discs

In Sections 3.3, 4.1, and 4.2, we assumed that the component interpreted as a stellar disc has an isotropic velocity dispersion tensor and we treated it as a fluid in the perturbation analysis. This approach has some limitations, which we discuss here. In the case of 2D disc models, comparisons between fluid and kinetic stability analyses (e.g. Rafikov 2001) have shown that modelling a stellar disc as a fluid turns out to be a sufficiently good approximation, provided the quantity entering the disc gravitational instability diagnostics is correctly interpreted in terms of the collisionless quantities. For instance, it is the radial stellar velocity dispersion, σ_⋆, R, that enters the definition of Q_⋆. Extending the same line of reasoning to our 3D analysis, when assuming that Q_3D is an approximate instability indicator for collisionless discs, we have to take care of relating the fluid quantities appearing in Q_3D to collisionless quantities. As it is clear from the analysis of N23 (see also Eqs. 22 and 24), the sound speed c_s = γσ, appearing in Eq. (5), derives from a radial derivative of the pressure, so in this case σ must be identified with σ_⋆, R. The quantity $c_{\text{s}}^2/\gamma ={\sigma }^2$, appearing in ν², (see Eq. 6), derives from a vertical derivative of the pressure, so in this case σ must be identified with the vertical velocity dispersion, σ_⋆, z. Finally, the scale height h_z, for a given midplane density, is expected to increase for increasing σ_⋆, z, while being independent of σ_⋆, R.

Focusing for simplicity on the midplane (where we expect ν ≈ 0) and assuming that γ = 1, we can thus rewrite Q_3D (Eq. 5) for a stellar disc as

$$\begin{array}{c}\hfill Q_{3\mathrm{D},\star }\approx \sqrt{\frac{{\kappa }^2}{4\pi G{\rho }_{0,\star }}}+\frac{{\sigma }_{\star ,R}}{{\sigma }_{\star ,z}}{\left(\frac{h_z}{{\sigma }_{\star ,z}/\sqrt{4\pi G{\rho }_{0,\star }}}\right)}^{-1}.\end{array}$$(40)

Given that neither the first term on the right-hand side nor the quantity in parentheses depend on the velocity ellipsoid, this equation shows that, for instance, a stellar disc with σ_⋆, R > σ_⋆, z has in fact higher Q_3D, ⋆ (and thus is more stable) than estimated under the assumption of isotropic (σ_⋆, R/σ_⋆, z = 1) velocity distribution. Similar considerations apply to the instability indicators considered in Sections 4.1 and 4.2, given the relationship between Q_3D, ⋆ and F_⋆.

There is observational evidence (Gerssen & Shapiro Griffin 2012; Marchuk & Sotnikova 2017; Pinna et al. 2018; Mackereth et al. 2019; Mogotsi & Romeo 2019) that in present-day stellar discs the velocity dispersion tensor is in general anisotropic, with typically σ_⋆, R > σ_⋆, z, a finding that is also supported by theoretical models (Rodionov & Sotnikova 2013; Walo-Martín et al. 2021). The velocity ellipsoid of stars in a higher-redshift disc is hard to measure observationally, but it is reasonable to expect that higher-redshift discs are dynamically younger, and thus more isotropic, having inherited the velocity distribution from the gas from which their stars have formed. This picture is supported by the finding that in the Milky Way the vertical scale height of the youngest stellar population is similar to that of the cold gas disc (e.g. Bacchini et al. 2019b)

5. One-component self-gravitating discs with polytropic vertical distributions

The vertical structure of discs is often modelled as isothermal, but in fact the velocity dispersion can in general have a vertical gradient. Here, we explore how the stability properties of a stratified disc depend on such a gradient, by exploring one-component self-gravitating discs with polytropic vertical distributions.

5.1. Equations

For a self-gravitating disc with a polytropic vertical distribution, at a given R

$$\begin{array}{c}\hfill p(z)=p_0{\left[\frac{\rho (z)}{{\rho }_0}\right]}^{{\gamma }^{\prime }},\end{array}$$(41)

where γ′ ≡ 1 + 1/n is the polytropic exponent, n is the polytropic index, ρ₀ = ρ(0), and p₀ = p(0). We assume that the gas is in vertical hydrostatic equilibrium in its own gravitational potential, so the vertical stratification is the same as that of a polytropic slab (chapter 2 of Horedt 2004; see also Ibanez & Sigalotti 1984). Eq. (41), combined with the hydrostatic equilibrium equation,

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}z}=-\rho \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z},\end{array}$$(42)

and with the Poisson equation,

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathrm{\Phi }}{\mathrm{d}{z}^2}=4\pi G\rho ,\end{array}$$(43)

gives for finite n the Lane-Emden equation,

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\theta }{\mathrm{d}{\zeta }^2}=\pm {\theta }^n,\end{array}$$(44)

where θ ≡ (ρ/ρ₀)^1/n and

$$\begin{array}{c}\hfill \zeta \equiv z/{[(n+1)p_0/(4\pi G{{\rho }_0}^2)]}^{1/2}.\end{array}$$(45)

On the right-hand side of Eq. (44), the sign is negative when −1 < n < ∞ (i.e. γ′ < 0  ∪  γ′ > 1), while it is positive when −∞< n < −1 (i.e. 0 < γ′ < 1). When n = ±∞ (γ′ = 1), combining Eqs. (41–43) we get

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\theta }{\mathrm{d}{\zeta }^2}=e^{-\theta },\end{array}$$(46)

whose analytic solution is the SGI slab described in Section 2.2. We assume that the fluid undergoes barotropic transformations of the form p ∝ ρ^γ, so dp = γ(p/ρ)dρ and $c_{\text{s}}^2=\gamma {\sigma }^2$. We limit ourselves to convectively stable polytropic distributions, that is those that satisfy the Schwarzschild stability criterion, γ′ ≤ γ. In particular, in this section we consider two representative specific cases.

1. A monoatomic ideal gas undergoing adiabatic transformations with γ = 5/3, for which $c_{\mathrm{s}}=\sqrt{5/3}\sigma$ and convectively stable distributions are those with γ′ ≤ 5/3 (i.e. n < 0  ∪  n ≥ 3/2).

2. A fluid undergoing isothermal transformations with γ = 1, for which c_s = σ and convectively stable distributions are those with γ′ ≤ 1 (n ≤ 0).

We note however that for γ′ < 0 (−1 < n < 0) the density increases with increasing distance from the midplane (Viala & Horedt 1974), so in all cases we exclude from our analysis γ′ < 0 (−1 < n < 0).

5.2. Vertical density and velocity-dispersion profiles

Given that we focus on the cases γ = 5/3 and γ = 1, and that we require convective stability (see Section 5.1), we present here polytropic distributions with a polytropic exponent in the range 0 < γ′ ≤ 5/3 (i.e. n ≥ 3/2  ∪  n < −1).

We numerically solved Eq. (44) with boundary conditions θ = 0 and dθ/∂ζ = 0 at ζ = 0 using the Mathematica routine NDsolve. We find it convenient to normalize the vertical coordinate, z, to the half-mass half-height, z_half. The vertical density and velocity-dispersion profiles are shown in Fig. 7. The velocity dispersion, σ, is defined by σ² = p/ρ, so the profiles in the lower panel of Fig. 7 can also be interpreted as profiles of $\sqrt{T}$ for a disc vertically supported by thermal pressure. We verified that the profiles shown in Fig. 7 are consistent with those tabulated by Ibanez & Sigalotti (1984) and Horedt (2004) for the values of n that we have in common.

Fig. 7. Vertical density (upper panel) and velocity-dispersion (lower panel) profiles (solid curves) of one-component self-gravitating slabs with polytropic vertical distributions for different values of the polytropic index, n. The case of n = ∞ (dotted curves) is the analytic SGI vertical distribution.

5.3. Stability

We discuss here the stability of the self-gravitating discs based on the 3D instability criterion (Eq. 5). When γ = 5/3, we explore polytropic distributions with a polytropic exponent in the range 0 < γ′ ≤ 5/3 (i.e. n ≥ 3/2  ∪  n < −1). When γ = 1, we explore polytropic distributions with a polytropic exponent in the range 0 < γ′ ≤ 1 (i.e. n < −1). With these choices, the systems are guaranteed to be convectively stable. Fig. 8 shows vertical Q_3D profiles for models with Q = 0.4 for a selection of values of n for γ = 5/3 (upper panel) and γ = 1 (lower panel). The Q_3D profiles are qualitatively similar for all values of γ and n, with instability close to the midplane and Q_3D increasing with z. The instability strip, |z|< z_inst (see Section 3.3), tends to be slightly wider for higher n, though z_inst/z_half spans a small range: 0.35 ≲ z_inst/z_half ≲ 0.54 for Q = 0.4. The value Q = 0.4 adopted in Fig. 8 is such to have an unstable region around the midplane for all the explored values of n. The effect of increasing Q is to ‘shift upwards’ the curves in Fig. 8, and thus decrease the extent of the unstable regions, which gradually shrink to zero, starting from polytropic models with the lowest positive values of n.

Fig. 8. Vertical profiles of the local gravitational instability parameter, Q3D, for discs with polytropic vertical density distributions assuming Q = 0.4, when γ = 5/3 (upper panel) or γ = 1 (lower panel). In each panel, the selected values of n are such that the distribution is convectively stable for the assumed γ. The case n = ∞ is the SGI, for which Q3D(z) is analytic. The horizontal dashed line indicates the instability threshold, Q3D = 1.

6. Conclusions

Building on the 3D instability analysis of N23, we have studied the local gravitational instability in two-component 3D axisymmetric discs. We have focused on two-component isothermal discs (without vertical velocity dispersion gradients), but we have complemented our analysis with one-component self-gravitating polytropic discs (with vertical velocity dispersion gradients). The main results of this paper are the following.

• Given that the effect of a second responsive disc is to favour the instability, for each component of a two-component disc Q_3D < 1 is a sufficient condition for instability. Under the assumption that the vertical distributions are isothermal, Q_3D can be computed for each component at a given radius, R, as a function of distance from the midplane, z, using observational estimates of the surface densities, velocity dispersions, and epicycle frequency, κ.

• At a given R and z, and at fixed surface density, velocity dispersion, and κ (thus at a fixed 2D Q parameter), the 3D Q_3D parameter of a disc can be higher, but it can also be lower in the presence of a second component than in the absence of it (see lower panel of Fig. 1). When present, the instability occurs in a strip enclosing the midplane with half-height z_inst that can be computed numerically from observationally inferred quantities.

• We derived a sufficient condition for local gravitational instability of two-component vertically stratified discs that takes into account the mutual back-reactions of the two discs. For instance, a disc consisting of a gaseous disc (index ‘gas’) and a stellar disc (index ‘⋆’), with Q_3D, gas > 1 and Q_3D, ⋆ > 1 at a given R, is unstable in the midplane against radial perturbations with wave number k_R, if F_{gas + ⋆}(k_R) > 1, where F_{gas + ⋆} is a function that can be expressed in terms of observationally inferred quantities (see Section 4.2).

• We have computed Q_3D(z), at a given R, for one-component self-gravitating discs with vertical polytropic distributions, and thus vertical temperature or velocity-dispersion gradients. For a range of values of the polytropic index, n, corresponding to convectively stable configurations, we have found a behaviour qualitatively similar to discs with vertical isothermal distributions (n = ∞). For unstable models, the height of the instability region increases only slightly with n, at fixed Q.

In conclusion, the explorations in this paper provide support to the proposal that Q_3D < 1 is a robust sufficient condition for local gravitational instability of discs, which depends only weakly on the presence of a second component and on the vertical velocity-dispersion gradient. When Q_3D > 1, the local gravitational instability of a two-component disc can be tested by applying the wave-number-dependent criterion (39). In conclusion, whenever the observational data allow one to reconstruct the 3D properties of discs (see Bacchini et al. 2024), 3D local gravitational instability criteria such as those analyzed in this work can be successfully employed.

Acknowledgments

We thank the referee Alessandro Romeo for useful comments that helped improve the paper. The research activities described in this paper have been co-funded by the European Union – NextGenerationEU within PRIN 2022 project n.20229YBSAN – Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch. CB acknowledges support from the Carlsberg Foundation Fellowship Programme by Carlsbergfondet.

References

Appendix A: Gravitational potential of radial perturbations in thick discs

When studying radial axisymmetric linear perturbations in vertically stratified discs, N23, following Goldreich & Lynden-Bell (1965a), assumed a perturbed Poisson equation in the form

$$\begin{array}{c}\hfill -({k_R}^2+h_z^{-2})\delta \mathrm{\Phi }=4\pi G\delta \rho ,\end{array}$$(A.1)

which approximately accounts for the finite vertical extent of the disc. In Eq. (A.1), δΦ is the perturbed potential, δρ is the perturbed density, k_R is the perturbation wave number, and h_z is the disc thickness. Following N23, we assume that h_z = h_70% (Eq. 7), so the above equation becomes

$$\begin{array}{c}\hfill -({k_R}^2+h_{70\%}^{-2})\delta \mathrm{\Phi }=4\pi G\delta \rho .\end{array}$$(A.2)

In this appendix, we assess quantitatively the validity of this approximation of the perturbed Poisson equation. We consider an axisymmetric gravitational potential perturbation, which in cylindrical coordinates has the form

$$\begin{array}{c}\hfill \delta \widehat{\mathrm{\Phi }}(R,z,t)=f(z)exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\mathrm{\Phi }}}_0,\end{array}$$(A.3)

with $\delta {\widehat{\mathrm{\Phi }}}_0=\delta \widehat{\mathrm{\Phi }}(0,0,0)$, and f(z) such that f(0) = 1 and lim_|z|→∞f(z) = 0, and a density perturbation in the form

$$\begin{array}{c}\hfill \delta \widehat{\rho }(R,z,t)=g(z)exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\rho }}_0,\end{array}$$(A.4)

with $\delta {\widehat{\rho }}_0=\delta \rho (0,0,0)$, and g(z) such that g(0) = 1 and lim_|z|→∞g(z) = 0. From the perturbed Poisson equation,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\delta \widehat{\mathrm{\Phi }}=4\pi G\delta \widehat{\rho },\end{array}$$(A.5)

assuming (as in the analysis of N23) that 1/R is negligible compared to |k_R|, we get

$$\begin{array}{c}\hfill \frac{1}{A_0}(-{k_R}^{2}f+\frac{{\mathrm{d}}^{2}f}{\mathrm{d}{z}^2})exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\mathrm{\Phi }}}_0=\frac{4\pi G}{A_0}g(z)exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\rho }}_0,\end{array}$$(A.6)

which is satisfied for

$$\begin{array}{c}\hfill g(z)=\frac{1}{A_0}(-{k_R}^{2}f+\frac{{\mathrm{d}}^{2}f}{\mathrm{d}{z}^2})\end{array}$$(A.7)

and

$$\begin{array}{c}\hfill \delta {\widehat{\mathrm{\Phi }}}_0=\frac{4\pi G}{A_0}\delta {\widehat{\rho }}_0,\end{array}$$(A.8)

where

$$\begin{array}{c}\hfill A_0\equiv -{k_R}^2+{\left(\frac{{\mathrm{d}}^{2}f}{\mathrm{d}{z}^2}\right)}_{z=0}\end{array}$$(A.9)

is a normalization factor that ensures that g(0) = 1. Indicating with ⟨⋯⟩ a weighted average in z, in our perturbed hydrodynamics equations we can approximate $\delta \widehat{\mathrm{\Phi }}$ with

$$\begin{array}{c}\hfill \delta \mathrm{\Phi }(R,t)\equiv ⟨\delta \widehat{\mathrm{\Phi }}⟩=⟨f⟩exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\mathrm{\Phi }}}_0.\end{array}$$(A.10)

Similarly, we can approximate $\delta \widehat{\rho }(R,z,t)$ with

$$\begin{array}{c}\hfill \delta \rho (R,t)\equiv ⟨\delta \widehat{\rho }⟩=⟨g⟩exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\rho }}_0,\end{array}$$(A.11)

where (see Eq. A.7)

$$\begin{array}{c}\hfill ⟨g⟩=\frac{1}{A_0}(-{k_R}^2⟨f⟩+\langle \frac{{\mathrm{d}}^{2}f}{\mathrm{d}{z}^2}\rangle ).\end{array}$$(A.12)

We specialize to the case in which

$$\begin{array}{c}\hfill f(z)=e^{-\frac{z^2}{2H^2}}\end{array}$$(A.13)

and adopt as the weighted average of a generic function F(z)

$$\begin{array}{c}\hfill ⟨F(z)⟩=\frac{{\int }_{-\infty }^{\infty }F(z)f(z)\mathrm{d}z}{{\int }_{-\infty }^{\infty }f(z)\mathrm{d}z}.\end{array}$$(A.14)

For f given by Eq. (A.13) $A_0=-(k_R{}^2+H^{-2})$ (Eq. A.9), so, combining Eq. (A.7) and Eq. (A.13), we get

$$\begin{array}{c}\hfill g(z)=\frac{{k_R}^2+H^{-2}-z^2{H}^{-4}}{{k_R}^2+H^{-2}}{e}^{-\frac{z^2}{2H^2}}.\end{array}$$(A.15)

Calculating explicitly the weighted averages over z, we get

$$\begin{array}{c}\hfill ⟨f⟩=\int f^2(z)\mathrm{d}z/\int f(z)\mathrm{d}z=\frac{1}{\sqrt{2}}\end{array}$$(A.16)

and

$$\begin{array}{c}\hfill ⟨g⟩=\int g(z)f(z)\mathrm{d}z/\int f(z)\mathrm{d}z=\frac{1}{\sqrt{2}({k_R}^2+H^{-2})}({k_R}^2+\frac{1}{2H^2}).\end{array}$$(A.17)

It follows that

$$\begin{array}{c}\hfill \delta \mathrm{\Phi }=\frac{1}{\sqrt{2}}exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\mathrm{\Phi }}}_0\end{array}$$(A.18)

and

$$\begin{array}{c}\hfill \delta \rho =\frac{{k_R}^2+\frac{1}{2}{H}^{-2}}{\sqrt{2}({k_R}^2+H^{-2})}exp[\mathrm{i}(k_{R}R-\omega t)]\delta {\widehat{\rho }}_0.\end{array}$$(A.19)

Combining the two above equations we get

$$\begin{array}{c}\hfill \delta \rho =-\frac{{k_R}^2+\frac{1}{2}{H}^{-2}}{{k_R}^2+H^{-2}}\frac{\delta {\widehat{\rho }}_0}{\delta {\widehat{\mathrm{\Phi }}}_0}\delta \mathrm{\Phi },\end{array}$$(A.20)

which, using Eq. (A.8), gives

$$\begin{array}{c}\hfill -({k_R}^2+\frac{1}{2}{H}^{-2})\delta \mathrm{\Phi }=4\pi G\delta \rho .\end{array}$$(A.21)

The above equation coincides with Eq. (A.2) if we assume that

$$\begin{array}{c}\hfill 2H^2=h_{70\%}^2.\end{array}$$(A.22)

We conclude that, approximating the gravitational potential of the perturbation with its vertical weighted average, the perturbed Poisson equation in the form (A.2) is valid for a density perturbation with a vertical profile

$$\begin{array}{c}\hfill \frac{\delta \widehat{\rho }(R,z,t)}{\delta \widehat{\rho }(R,0,t)}=g(z)=\frac{{k_R}^2+H^{-2}-z^2{H}^{-4}}{{k_R}^2+H^{-2}}{e}^{-\frac{z^2}{2H^2}},\end{array}$$(A.23)

with $H^2=h_{70\%}^2/2$.

We now want to verify whether assuming a perturbation with this vertical structure is consistent with the vertical structure of the unperturbed system. We consider, for instance, an unperturbed SGI vertical density profile (8): in this case h_70% ≃ 1.7346b, so the assumption (A.22) becomes

$$\begin{array}{c}\hfill H=\frac{1.7346}{\sqrt{2}}b=1.2265b.\end{array}$$(A.24)

Assuming that $H^2=h_{70\%}^2/2$, in the left panel of Fig. A.1 we compare, for different values of k_RH, the normalised vertical density profile of the perturbation, $\delta \widehat{\rho }(R,z,t)/\delta \widehat{\rho }(R,0,t)=g(z)$ (with g given by Eq. A.15), with the normalised unperturbed vertical density profile, ρ(z)/ρ(0). The fact that the shape of the density profile of the perturbation is similar to that of the unperturbed density distribution suggests that the considered disturbance can be consistent with the hypothesis of small perturbations |δρ|/ρ ≪ 1. This is illustrated in the right panel of Fig. A.1, showing the ratio $|\delta \widehat{\rho }|/\rho$ as a function of z, assuming, for instance, that $\delta \widehat{\rho }/\rho =0.01$ in the midplane. We note that in Fig. A.1 the plotted range in z is such that it contains 90% of the mass per unit surface of the unperturbed distribution.

Fig. A.1. Left panel: Normalised vertical profiles of the density perturbation $\delta \widehat{\rho }$ (Eqs. A.4 and A.15) with $H=h_{70\%}/\sqrt{2}$ for different values of the radial wave number kR (solid curves), compared with the normalised unperturbed vertical density profile (Eq. 8; dashed curve). Right panel: Vertical profiles of the ratio between the density perturbation and the unperturbed density profile for the same values of kR as in the left panel, assuming that $\delta \widehat{\rho }/\rho =0.01$ at z = 0. The plotted range in z encloses 90% of the mass per unit surface of the unperturbed distribution.

All Figures

Fig. 1. Upper panel: Vertical density profiles (solid curves) of gas and stars of a two-component isothermal disc with μ = 9 and ξ = 0.3. The dashed curves are the profiles of the corresponding SGI models with the same velocity dispersion and surface density. Lower panel: Vertical Q3D profiles of the same models as in the upper panel, assuming Qgas = 0.6 (and thus Q⋆ = 1.4). The horizontal dashed line indicates the instability threshold, Q3D = 1. Here, ρ0, gas and bgas are, respectively, the midplane density and the scale height of the gaseous component of the two-component disc.

In the text

	Fig. 2. Upper panel: Stars-to-gas thickness ratio, hz, ⋆/hz, gas, as a function of the squared velocity-dispersion ratio, μ, for two-component isothermal slabs with different values of the midplane density ratio ξ (hz = h70%, defined by Eq. 7). Overplotted in grey is the power law, μ5/8. Lower panel: Percent residuals between the points and the power law of the upper panel.
In the text

	Fig. 3. Ratios hz, ⋆/zhalf, ⋆ (filled symbols) and hz, gas/zhalf, gas (empty symbols) as functions of the squared velocity-dispersion ratio, μ, for two-component isothermal slabs with different values of the midplane density ratio, ξ. hz = h70% and zhalf are defined by Eq. (7) and Eq. (2), respectively.
In the text

Fig. 4. Upper panel: Midplane values of Q3D, ⋆ (filled symbols) and Q3D, gas (empty symbols) as functions of the midplane star-to-gas density ratio, ξ, for different values of the star-to-gas squared velocity dispersion ratio, μ. The horizontal dashed line indicates the instability threshold, Q3D = 1. Lower panel: Extent of the unstable strip above the midplane for the stellar (zinst, ⋆, normalised to zhalf, ⋆; filled symbols) and gaseous (zinst, gas, normalised to zhalf, gas; empty symbols) components as functions of ξ for different values of μ. In both panels, we assume that Qgas = 0.6.

In the text

	Fig. 5. Same as Fig. 1, but for a two-component isothermal disc model with ξ = 0.75, μ = 6, Qgas = 1.1, and Q⋆ = 1.2.
In the text

	Fig. 6. Functions Fgas(kR), F⋆(kR), and Fgas + ⋆(kR) (see Eqs. 37, 38 and 39) for the same model as in Fig. 5. The unstable perturbations are those with Fgas + ⋆ > 1.
In the text

	Fig. 7. Vertical density (upper panel) and velocity-dispersion (lower panel) profiles (solid curves) of one-component self-gravitating slabs with polytropic vertical distributions for different values of the polytropic index, n. The case of n = ∞ (dotted curves) is the analytic SGI vertical distribution.
In the text

Fig. 8. Vertical profiles of the local gravitational instability parameter, Q3D, for discs with polytropic vertical density distributions assuming Q = 0.4, when γ = 5/3 (upper panel) or γ = 1 (lower panel). In each panel, the selected values of n are such that the distribution is convectively stable for the assumed γ. The case n = ∞ is the SGI, for which Q3D(z) is analytic. The horizontal dashed line indicates the instability threshold, Q3D = 1.

In the text

Fig. A.1. Left panel: Normalised vertical profiles of the density perturbation $\delta \widehat{\rho }$ (Eqs. A.4 and A.15) with $H=h_{70\%}/\sqrt{2}$ for different values of the radial wave number kR (solid curves), compared with the normalised unperturbed vertical density profile (Eq. 8; dashed curve). Right panel: Vertical profiles of the ratio between the density perturbation and the unperturbed density profile for the same values of kR as in the left panel, assuming that $\delta \widehat{\rho }/\rho =0.01$ at z = 0. The plotted range in z encloses 90% of the mass per unit surface of the unperturbed distribution.

In the text