© S. Sarkar and C. J. Jog 2022

1. Introduction

The vertical structure of the stellar disk in a galaxy has been studied as a self-gravitating, single-component, isothermal disk in the literature, where its self-consistent vertical density distribution is given by a sech² form (Spitzer 1942). However, a real galactic disk is a multi-component system of gravitationally coupled stars and interstellar gas (HI and H₂) embedded in the potential of the dark matter halo. The self-consistent vertical distribution of stars in such a system is determined by the joint gravitational potential of stars, gas, and the dark matter halo (Narayan & Jog 2002; Sarkar & Jog 2018), instead of its self-gravitational force alone. The joint potential of the system is found to constrain the distribution of stars towards the mid-plane, and thus increases the mid-plane density value and decreases the vertical disk thickness (Sarkar & Jog 2018). In the inner Galaxy, gas plays the dominant role in constraining the distribution of stars (Sarkar & Jog 2018). Although gas contains 10–15% of the disk mass (Young & Scoville 1991; Binney & Merrifield 1998), it forms a thin layer about the mid-plane due to its low vertical velocity dispersion and therefore can affect the vertical distribution of stars significantly in the inner Galaxy. Through gravitational coupling, stars also constrain the gas distribution in a similar way. Stars, being a much more massive component, have a stronger effect on gas. Therefore we expect the vertical distribution of stars and gas to be more strongly bound in a coupled system than in single-component cases, and thus it is more robust against perturbations or external tidal interactions.

In this context, it is interesting to study the gravitational potential energy of coupled stellar and gas distribution. The higher the potential energy of the vertical distribution of stars (gas), the more difficult it should be to distort the disk in external tidal encounters. Motivated by the results of the constraining effect of the joint potential, we expect the potential energy of a component to be higher in the gravitationally coupled multi-component system than in its single-component case. With this aim, we study the potential energy per unit area of the vertical distribution of stars, and gas in the single-component self-gravitating cases as well as in the coupled, multi-component system. We consider a gravitationally coupled stars-plus-gas disk, a two-component system, in the inner Galaxy, and explicitly derive the expression for the potential energy per unit area of the disk, following the method proposed by Camm (1967) for the single-component case. We study how the energy corresponding to each component is affected by the gravitational coupling between them. Surprisingly, we find that the energy of each component remains unchanged, and the physical reason for this is explained in the paper. However, despite this, we found that stars and gas are more strongly bound to the mid-plane in the coupled case, and thus are less susceptible to external tidal distortions. For completeness, we have also studied the potential energy per unit area of the disk for a three-component case, that is, for a gravitationally coupled disk of stars and two gas components.

In the outer Galaxy, on the other hand, the dark matter halo plays the dominant role in constraining the vertical distribution of stars and gas significantly (Sarkar & Jog 2018). For simplicity, we did not include the dark matter halo to study the multi-component system here because the dark matter halo is shown to have a less significant effect than gas (stars) on the vertical distribution of stars (gas) in the inner Galaxy (Narayan & Jog 2002; Sarkar & Jog 2018).

We also note that we ignored the bulge in the inner Galaxy. Our theoretical model is applied from a galactocentric radius of 4 kpc onward. The bulge is not a dominant gravitating component in the region studied here (Ghosh et al. 2016; Blum 1995).

We show the formulation in Sect. 2, results in Sect. 3, and give discussion and conclusions in Sects. 4 and 5, respectively. The formulation of the potential energy, given in Sect. 2, is general and applicable to any two (three)-component disk, for instance, for an n-component stellar disk, even though we apply it to the specific and observationally motivated case of a stars-plus-gas disk. In other words, the formulation of the energy does not in any way involve the specific physical nature and properties of gas, such as dissipation or low dispersion.

2. Formulation of the problem

2.1. Gravitational potential energy of a single-component isothermal galactic disk

First we discuss the formulation of the gravitational potential energy of a single component self-gravitating galactic disk that can be taken to be a stars-alone or a gas-alone disk. For a mass distribution contained in a finite volume of space, the gravitational force decreases as ∼1/r² at a large distance and hence the work done to bring a unit mass from infinity to a certain finite distance is obtained to be a finite quantity. The gravitational potential energy of the mass distribution in such a case is defined as the energy released in assembling the finite system from an infinitely dispersed state. For a galactic disk, however, the mass distribution is stratified in plane parallel layers and is infinite on the x − y plane. The vertical gravitational force for such a stratified mass distribution remains constant at a value of −2πGΣ at large z at any galactocentric radius R, Σ being the surface density of the disk at that R, beyond the vertical extent of the mass distribution (and within the disk approximation limit). This results in the energy released in bringing a unit mass from infinity to a certain finite height to be infinite. Therefore, instead, the state of complete collapse of the disk mass on the z = 0 plane is defined here to be the state of the zero potential energy, and the work required to build the disk from that state is considered to be the potential energy stored in the disk. For a detailed discussion of this point and the derivation of Eq. (1) (see below) for a single-component case, see Camm (1967).

We note that the mass contained in a column of unit cross section perpendicular to the mid-plane is finite. Therefore the potential energy of the disk is defined in terms of the energy contained in a column of unit cross section, that is, as the potential energy per unit area of the disk. We use the galactocentric cylindrical coordinates (R, ϕ, z) and consider the disk to be axisymmetric.

The mathematical expression for the gravitational potential energy per unit area of a stars-alone disk has been derived in Camm (1967) as

$$\begin{array}{c}\hfill W=-{\displaystyle {\int }_{-\infty }^{+\infty }}\rho (z)\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z}z\mathrm{d}z,\end{array}$$(1)

where ρ(z) is the vertical mass density distribution of stars, Φ is the gravitational potential of the disk, and dΦ/dz is taken to be the force per unit mass (K_z) due to the self-gravity of the stellar disk, acting along the negative z direction. This represents the self-gravitational energy per unit area of the disk. A similar expression was used by Garrido Pestaña & Eckhardt (2010) and Fridman et al. (1984) to calculate the gravitational potential energy of a single-component disk, but the expression was not derived.

The above expression in Camm (1967) was obtained by using the Poisson equation in the form of d²Φ/dz² = −4πGρ(z), where the vertical force was defined as K_z = dΦ/dz, which is a negative quantity. However, we adopted the standard notation here that is routinely used in the literature, where the Poisson equation is given by d²Φ/dz² = 4πGρ(z) and the vertical force is defined by K_z = −dΦ/dz, where dΦ/dz is positive. Following the treatment in Camm (1967), we therefore derive the expression of the energy as

$$\begin{array}{cc}\hfill W& ={\displaystyle {\int }_{-\infty }^{+\infty }}\rho (z)\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}z}z\mathrm{d}z\hfill \\ \hfill & =-{\displaystyle {\int }_{-\infty }^{+\infty }}\rho (z)K_{z}z\mathrm{d}z.\hfill \end{array}$$(2)

A detailed derivation for a multi-component case is given in Sect. 2.2. The negative sign in front of Eqs. (1) or (2) indicates that the energy is positive.

We note a few important points here. While deriving Eq. (1), Camm (1967) used only the z-term in the Poisson equation. We show that this treatment is justified for a thin galactic disk. We show that the z-term of the Poisson equation is much greater than the R-term ¹. Thus the density distribution, force, and energy become only z-dependent quantities. We followed the same approach while deriving the energy for the multi-component case in Sect. 2.2. We also note that Camm (1967) implicitly assumed the disk to have a constant radial surface density, whereas we considered realistic stellar and gas disks of radially varying surface density. This does not affect the derivation of the energy or the application of the model because the calculation is local. We used the surface density value at any given radius as a local constraint to obtain ρ(z) at a given R, which was then used in the expression for the energy, as discussed below and in Sect. 2.3. Thus the calculation is local, and the energy is independent of the value of surface density at other radii.

We note that although Eq. (2) is derived in Camm (1967), using explicitly only the Poisson equation, for a disk in vertical hydrostatic equilibrium, ρ(z) and K_z are related to each other and have to be obtained by solving the joint hydrostatic balance-Poisson equation. These solutions are required to obtain a numerical value for W. We show this set of equations below for a single-component self-gravitating isothermal disk.

We assume the vertical velocity dispersion (σ_z) of the component to be isothermal along z. The vertical hydrostatic balance equation for a single-component isothermal self-gravitating disk is given by

$$\begin{array}{c}\hfill \frac{{\sigma }_z^2}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}z}=K_z.\end{array}$$(3)

The Poisson equation for a single-component galactic disk is given as

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

We combine these two equations to obtain the joint hydrostatic balance-Poisson equation

$$\begin{array}{c}\hfill {\sigma }_z^2\frac{\mathrm{d}}{\mathrm{d}z}\left[\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}z}\right]=-4\pi G\rho .\end{array}$$(5)

The analytical solution of this equation in form of sech²(z/z₀) was obtained by Spitzer (1942). The analytical expressions of the density distribution and the force are given as

$$\begin{array}{c}\hfill \rho (z)={\rho }_0{\mathrm{sech}}^2(z/z_0);|K_z|=2\frac{{\sigma }_z^2}{z_0}tanh\left(\frac{z}{z_0}\right),\end{array}$$(6)

respectively, where $z_0={({\sigma }_z^2/2\pi G{\rho }_0)}^{1/2}$. Using the analytical form of the solution ρ(z) versus z, we calculate ρ₀ (mid-plane density) and z₀ (scale-height) by integrating ρ(z) vs. z, using the constraint of the observed surface density, defined as $\mathrm{\Sigma }={\int }_{-\infty }^{+\infty }\rho \mathrm{d}z$, at a given galactocentric radius R. The extent of z in the numerical integration is chosen such that the solution obtained is saturated. The obtained value of z₀ in turn is used to calculate |K_z| up to the same limit of z. Thus using ρ(z) and K_z we calculate W using Eq. (2), which gives us the saturated value of the energy.

However, for a gravitationally coupled two-component system, the joint hydrostatic balance-Poisson equation has to be solved numerically to obtain ρ(z) vs. z and |K_z| vs. z. This is discussed in Sect. 2.3.1.

2.2. Gravitational potential energy of a multi-component gravitationally coupled isothermal galactic disk

In this section, we explicitly derive the expression of the gravitational potential energy per unit area for a realistic model of galactic disk that is a multi-component system of gravitationally coupled stars and gas, taken at a given galactocentric radius R. We consider the disks of stars and gas to be coplanar with the same mid-plane at z = 0. We consider the disk to be a thin disk and therefore consider only the z-term in the corresponding Poisson equation for the multi-component system (see Sect. 2.1). We show the detailed derivation for a two-component system, consisting of stars and one gas component, and discuss the same for a three-component system of stars and two gas components in Sect. 3.4. The formulation in each case is done following the same steps as in the treatment for the stars-alone case in Camm (1967). In the following, the subscripts i = s, g in the quantities ρ(z), Σ, Φ denote stars and gas, respectively.

We consider the galactic disk to be a gravitationally coupled stars-plus-gas system where the vertical distribution of each component is determined by the joint gravitational force from the stars and gas. We take the vertical velocity dispersion of each of the two components to be isothermal. We assume that initially, all the mass of the stars-plus-gas system lies on the z = 0 plane, and the potential energy W of this system is assumed to be zero based on the same arguments as discussed in Sect. 2.1. Now the work done per unit area to build a gravitationally coupled stars-plus-gas disk together, at the same time, is stored as the gravitational potential energy per unit area of the stars-plus-gas disk.

Below we present the detailed derivation of the energy per unit area of this coupled stars-plus- gas disk following all the steps used for the stars-alone case in Camm (1967). As stated earlier in Sect. 2.1, we used d²Φ/dz² = 4πGρ(z) and K_z = −dΦ/dz in each step of the derivation, as shown below.

For a gravitationally coupled stars-plus-gas disk, the Poisson equation is given as

$$\begin{array}{cc}\hfill \frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{s}}}{\mathrm{d}{z}^2}+\frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{g}}}{\mathrm{d}{z}^2}& =4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\hfill \\ \hfill \mathrm{Or},\frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}{z}^2}& =4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}}).\hfill \end{array}$$(7)

In the first step, we integrate Eq. (7) and derive an expression for the vertical force of the coupled system to be (after some algebraic manipulations)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}z}& =2\pi G{\displaystyle {\int }_{-\infty }^z({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z-2\pi G{\int }_z^{\infty }({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z}\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill K_{z,\mathrm{coupled}}& =-\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}z}\hfill \\ \hfill & =2\pi G{\displaystyle {\int }_z^{\infty }({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z-2\pi G{\int }_{-\infty }^z({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z.}\hfill \end{array}$$(9)

Now while building up the gravitationally coupled stars-plus-gas disk from z = 0, at any intermediate step, only a fraction of the total mass of stars-plus-gas is distributed along z, denoted by the density distribution of ϵ(ρ_s + ρ_g) in the region z > 0, where ϵ lies between 0 and 1. The rest of the mass still lies on the z = 0 plane with the mass per unit area as $({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})-ϵ{\int }_0^{\infty }({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z$. The value of ϵ at any intermediate step has no specific physical meaning. It was used only to denote a fraction of the sum of the final vertical distribution of stars-plus-gas (i.e., (ρ_s + ρ_g)), at an intermediate step. We also note that ϵ has been taken to be independent of z in Camm (1967).

At any intermediate step, the force against which the work is being done is due to the joint stars-plus-gas system at that step. For any height x(> 0) above the z = 0 plane, the force (K_{fraction, coupled}), for such a system can be written as

$$\begin{array}{cc}\hfill K_{\mathrm{fraction},\mathrm{coupled}}& =2\pi G{\displaystyle {\int }_x^{\infty }ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z-2\pi G{\int }_0^xϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z}\hfill \\ \hfill & -2\pi G(({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})-{\displaystyle {\int }_0^{\infty }ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z).}\hfill \end{array}$$

Here the first two terms account for the force from the stars-plus-gas distribution along z > 0, and the rest of the terms account for the force due to the rest of the mass lying on the z = 0 plane. Now, we can write

$$\begin{array}{cc}\hfill K_{\mathrm{fraction},\mathrm{coupled}}& =2\pi G{\displaystyle {\int }_x^{\infty }ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z-2\pi G{\int }_0^xϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z}\hfill \\ \hfill & -2\pi G(({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})-{\displaystyle {\int }_0^{\infty }ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z)}\hfill \\ \hfill & -2\pi G{\displaystyle {\int }_{-\infty }^0ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z+2\pi G{\int }_{-\infty }^0ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z}\hfill \\ \hfill & =2\pi G{\displaystyle {\int }_x^{\infty }ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z-2\pi G{\int }_{-\infty }^xϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z}\hfill \\ \hfill & -2\pi G({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})+2\pi ϵG({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})\hfill \\ \hfill & =2\pi G{\displaystyle {\int }_x^{\infty }ϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z-2\pi G{\int }_{-\infty }^xϵ({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z}\hfill \\ \hfill & -2\pi G({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})(1-ϵ).\hfill \end{array}$$

The above equation can be rewritten, using Eq. (8), as

$$\begin{array}{c}\hfill K_{\mathrm{fraction},\mathrm{coupled}}=-ϵ\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}x}-2\pi G({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})(1-ϵ).\end{array}$$(10)

Now the work done per unit area of the z-plane to raise the density of stars plus gas between z and z + δz from ϵ(ρ_s + ρ_g) to (ϵ + δϵ)(ρ_s + ρ_g) is the work done in raising the mass (ρ_s + ρ_g)δϵδz from z = 0 to the chosen z plane, and is given by

$$\begin{array}{c}\hfill -({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\delta ϵ\delta z{\displaystyle {\int }_{x=0}^z-\{ϵ\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}x}+2\pi G({\mathrm{\Sigma }}_{\mathrm{s}}+{\mathrm{\Sigma }}_{\mathrm{g}})(1-ϵ)\}\mathrm{d}x.}\end{array}$$(11)

We now calculate the potential energy per unit area of the two-component disk of ρ_s + ρ_g by integrating the above equation over ϵ (from 0 to 1) and z (from −∞ to ∞), following the steps shown in Camm (1967) for the one-component case.

We obtain the expression for the energy as

$$\begin{array}{c}\hfill W_{\mathrm{coupled}}={\displaystyle {\int }_{-\infty }^{\infty }z\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}z}({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\mathrm{d}z.}\end{array}$$(12)

This is the most important result of this paper. Now the above equation can be further expressed as

$$\begin{array}{cc}\hfill W_{\mathrm{coupled}}& ={\displaystyle {\int }_{-\infty }^{\infty }z\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}z}{\rho }_{\mathrm{s}}\mathrm{d}z+{\int }_{-\infty }^{\infty }z\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}z}{\rho }_{\mathrm{g}}\mathrm{d}z}\hfill \\ \hfill & =-{\displaystyle {\int }_{-\infty }^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{s}}\mathrm{d}z-{\int }_{-\infty }^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{g}}\mathrm{d}z.}\hfill \end{array}$$

Using the symmetry of ρ(z) vs. z about z = 0, this can be written as

$$\begin{array}{c}\hfill W_{\mathrm{coupled}}=-2{\displaystyle {\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{s}}\mathrm{d}z-2{\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{g}}\mathrm{d}z,}\end{array}$$(13)

where the first integration can be considered to represent the potential energy per unit area of the stellar disk in the coupled stars-plus-gas system, and the second integration to represent the potential energy per unit area of the gas disk in the coupled stars-plus-gas system. In the limit of ρ_g → 0 or ρ_s → 0, this goes over to the one-component case (see Eqs. (2) and (9)). The separation of energy into the two components, as can be seen from the above equation, may appear somewhat surprising, but we note that all the expressions starting from Eq. (9) can be written in separable form for stars and gas, except that the components are being built against the same coupled force.

The above formulation can be used for any n-component system (e.g., for n > 2). For illustration, we show this for a three-component disk consisting of stars and two gas components in a similar fashion as discussed in Sect. 3.4.

2.3. Theoretical model for the vertical distribution for a multi-component disk and input parameters

2.3.1. Self-consistent vertical distribution in a multi-component disk

For a gravitationally coupled two-component stars-plus-gas disk, the hydrostatic balance of each of the components is determined by the joint gravitational force of stars and gas, and is given by

$$\begin{array}{c}\hfill \frac{{\sigma }_{z,i}^2}{{\rho }_i}\frac{\mathrm{d}{\rho }_i}{\mathrm{d}z}=K_{z,\mathrm{s}}+K_{z,\mathrm{g}},\end{array}$$(14)

where i represents stars (s) or gas (g) and the right-hand side of the equation represents the vertical force from the coupled stars-plus-gas case. For a gravitationally coupled stars-plus-gas galactic disk, the Poisson equation is given as

$$\begin{array}{cc}\hfill \frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{s}}}{\mathrm{d}{z}^2}+\frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{g}}}{\mathrm{d}{z}^2}& =4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})\hfill \\ \hfill \mathrm{Or},\frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}{z}^2}& =4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}}).\hfill \end{array}$$(15)

We combine these two equations to write the joint hydrostatic balance-Poisson equation as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\rho }_i}{\mathrm{d}{z}^2}=\frac{{\rho }_i}{{\sigma }_{z,i}^2}[-4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{g}})]+\frac{1}{{\rho }_i}{\left(\frac{\mathrm{d}{\rho }_i}{\mathrm{d}z}\right)}^2.\end{array}$$(16)

These coupled equations are solved numerically using the fourth-order Runge-Kutta method to obtain ρ_i(z) vs. z, simultaneously for each ith component in an iterative fashion, as discussed in Narayan & Jog (2002) and Sarkar & Jog (2018), until the fifth decimal convergence in the solutions. To solve the equations at a radius, we use the observed surface density of each component as one boundary condition, and dρ_i/dz = 0 at z = 0 as the other boundary condition, where the latter is true for any realistic distribution that is homogeneous very close to the mid-plane. The vertical distributions for a three-component system can be obtained following a similar method, as discussed in Sect. 3.4.

2.3.2. Input parameters

The formulation presented so far is general. Here we apply it for the Milky Way. We considered HI as the gas component here to study the two-component system of the stars-plus-gas disk, and chose the solar radius, taken to be at R = 8.5 kpc, to illustrate the results. The stellar disk was taken to be exponential with the central surface density Σ₀ = 640.9 M_⊙ pc⁻² and radial scale length R_D = 3.2 kpc (Mera et al. 1998). Thus the surface density value of the stellar disk is 45.0 M_⊙ pc⁻² at the solar radius.

The radial velocity dispersion values of stars on the mid-plane were obtained observationally by Lewis & Freeman (1989) up to R = 16 kpc. It falls off exponentially with radius as σ_R,s = 105 exp (−R/8.7 kpc) km s⁻¹. We calculated the corresponding vertical velocity dispersion value (σ_z, s) on the mid-plane by assuming the vertical to radial dispersion ratio to be 0.45 (Dehnen & Binney 1998; Mignard 2000), as observed in the solar neighbourhood. Thus σ_z, s at the solar radius is calculated to be 17.8 km s⁻¹. The dispersion was taken to be isothermal along z.

The surface density value of the HI disk was taken to be 5.5 M_⊙ pc⁻² (Scoville & Sanders 1987). The vertical velocity dispersion of HI at the solar radius was taken to be 8 km s⁻¹ and isothermal, based on the values given by Spitzer (1978) for the Galaxy, and Lewis (1984) for nearly 200 face-on galaxies. The corresponding input parameters to study a three-component system are discussed in Sect. 3.4.

3. Results

We first show the results for the two-component stars-plus-gas (HI) disk in Sects. 3.1, 3.2, and 3.3 in detail. We then show the results for a three-component disk for completeness, taking stars and two gas components (HI and H₂) in Sect. 3.4.

3.1. Calculation of the potential energy of the two-component stars-plus-gas disk

We calculated the vertical distributions of stars and gas (HI) at R = 8.5 kpc first, taking each as a single-component self-gravitating system (as discussed in Sect. 2.1), and then for the coupled two-component system of stars plus gas (as discussed in Sect. 2.3). We compare the vertical density distributions of stars in Fig. 1a and gas in Fig. 1b in these two cases. We note that the vertical distribution of each of the components solved in the coupled system is constrained towards the mid-plane. It has a higher mid-plane density that falls off more sharply along z and therefore has a smaller disk thickness than those in the corresponding one-component case. This is due to an additional gravitational force from a second component in the coupled system. We show the mid-plane density (up to third decimal place) and the half width at half maximum (HWHM) values of the density distribution that define the disk thickness in Table 1. We note that due to the higher mass content of the stellar distribution, stars affect gas more strongly.

Fig. 1. Stellar and gas (HI) vertical density distribution, ρ(z) vs. z at R = 8.5 kpc (panels a and b, respectively). In each plot, the solid curve represents the vertical distribution in the one-component case. The dashed curve represents the distribution obtained in the gravitationally coupled stars-plus-gas system. The coupled gravitational force of the star- plus-gas system constrains each distribution towards the mid-plane, i.e., increases the mid-plane density value of the distribution, which causes the scale height to become smaller and the curve to become steeper.

We also calculated the self-gravitational forces of the stars-alone and gas-alone cases and compare them with the gravitational force of the coupled stars-plus-gas system in Fig. 2. The self-gravitational force for each single-component case was calculated using Eq. (6), as discussed in Sect. 2.1. For the stars-plus-gas system, the coupled gravitational force is given by the right-hand side of Eq. (14). Each component is kept in hydrostatic equilibrium due to this coupled force. We calculated this force numerically from the left-hand side of this equation, substituting ρ_i and dρ_i/dz for any component, both obtained numerically by solving Eq. (16). Figure 2 shows that at each z distance from the mid-plane, this coupled force that keeps stellar and gas distributions in hydrostatic equilibrium in the coupled system is higher than the one-component self-gravitational forces.

Fig. 2. Vertical force per unit mass, i.e., |Kz| vs. z at R = 8.5 kpc. The solid and dashed curves represent the force due to the self-gravity of stars and gas (HI), respectively. The dash-dotted curve represents the force arising from the gravitationally coupled stars-plus-gas system, which is higher than the forces due to the self-gravity of the individual components at each z.

Motivated by these results, we now aim to compare the potential energy values of the stellar and gas distribution in the coupled case to the corresponding single-component cases. We used Eq. (2) and the procedure outlined in Sect. 2.1 to calculate the potential energy per unit area of a stars-alone disk at the solar radius to be 14257.7 M_⊙ pc⁻²km² s⁻² (3 × 10⁷ erg cm⁻²) and of a gas-alone disk to be 351.4 M_⊙ pc⁻² km² s⁻² (7.4 × 10⁵ erg cm⁻²). The sum of the energies is then given by 14609.1 M_⊙ pc⁻² km² s⁻². The z ranges used in the integration of ρ(z) vs. z and in calculation of K_z vs. z were chosen such that the energy values are numerically saturated (see Sect. 2.1).

Now we calculate the potential energy per unit area of the coupled stars-plus-gas disk using Eq. (13) (using K_z, coupled and ρ_i(z) obtained numerically as described in Sect. 2.3.1) and find it to be 14609.02 M_⊙ pc⁻² km² s⁻². This matches the sum of the energies for the stars-alone and the gas-alone cases within the numerical accuracy. Thus very interestingly, despite being in the gravitational force of the coupled system, the work done required to build up a stars-plus-gas disk turns out to be the same as the sum of the energies that would be required to build separate single-component self-gravitating stellar and gas disks. Importantly, the energy values for both stars and gas, in this case, obtained from Eq. (13) is the same as in the corresponding single-component self-gravitating cases within the numerical accuracy, that is, 14257.02 M_⊙ pc⁻² km² s⁻² and 352.0 M_⊙ pc⁻² km² s⁻², respectively. We plot the energy integrand, that is, −zρ(z)K_z vs. z in the single-component and in the coupled case for stars in Fig. 3a and for gas in Fig. 3b. Twice the area under these curves gives the corresponding energy values. We note that the energy integrand is now redistributed along z, conserving the area under the curve.

Fig. 3. Integrand −zρ(z)Kz vs. z (from Eqs. (2), (13)) of W (potential energy per unit area) for stars and gas (HI) at R = 8.5 kpc. Panela: energy integrand for stars in the stars-alone case (solid curve) and in the coupled stars-plus-gas case (dashed curve). Twice the area under the curves gives the energy per unit area of the stellar disk in the two cases, which is obtained to be the same. Panel b: energy integrand for gas (HI) in the gas-alone case (solid curve) and in the coupled stars-plus-gas case (dashed curve). Twice the area under the curves gives the energy per unit area of the gas disk in the two cases, which are obtained to be the same.

To investigate the physical reason that the same potential energy value was obtained, we simplified the expression of the energy further analytically. We discuss this in the following section.

3.2. Analytical simplification of the expression of gravitational potential energy per unit area of a multi-component disk

First we considered a single-component case, for instance the stars-alone case. We substituted the expression of K_z by the left-hand side of the hydrostatic balance equation (Eq. (3)) for a self-consistent distribution for a stars-alone disk into Eq. (2) and calculated it further as shown below,

$$\begin{array}{cc}\hfill W_{\text{stars-alone}}& =-2{\displaystyle {\int }_0^{\infty }z{K}_{z,\mathrm{s}}{\rho }_{\mathrm{s}}\mathrm{d}z}\hfill \\ \hfill & =-2{\displaystyle {\int }_0^{\infty }z\frac{{\sigma }_{z,\mathrm{s}}^2}{{\rho }_{\mathrm{s}}}\frac{\mathrm{d}{\rho }_{\mathrm{s}}}{\mathrm{d}z}{\rho }_{\mathrm{s}}\mathrm{d}z.}\hfill \end{array}$$

For an isothermal dispersion, we obtain this to be

$$\begin{array}{c}\hfill W_{\text{stars-alone}}=-2{\sigma }_{z,\mathrm{s}}^2{\displaystyle {\int }_0^{\infty }z\frac{\mathrm{d}{\rho }_{\mathrm{s}}}{\mathrm{d}z}dz.}\end{array}$$

Applying the method of integration by parts, we obtain

$$\begin{array}{c}\hfill W_{\text{stars-alone}}=-2{\sigma }_{z,\mathrm{s}}^2[(z{\rho }_{\mathrm{s}}{)}_0^{\infty }-{\displaystyle {\int }_0^{\infty }}{\rho }_{\mathrm{s}}\mathrm{d}z].\end{array}$$(17)

Now at a very large z, theoretically, the density value is zero. Here, because ρ_s falls off faster than 1/z (see Fig. 1), the value of zρ_s at the upper limit will tend to zero. Using the numerically obtained solutions also, we can say that at the edge of the distribution, that is, at large z by which the distribution is saturated and the value of the density is negligible, the product zρ_s becomes very small with respect to the other term. Thus the integration is

$$\begin{array}{cc}\hfill W_{\text{stars-alone}}& \approx -2{\sigma }_{z,\mathrm{s}}^2[-{\displaystyle {\int }_0^{\infty }}{\rho }_{\mathrm{s}}\mathrm{d}z]\hfill \\ \hfill & =-2{\sigma }_{z,\mathrm{s}}^2(-\frac{{\mathrm{\Sigma }}_{\mathrm{s}}}{2})\hfill \\ \hfill & ={\sigma }_{z,\mathrm{s}}^2{\mathrm{\Sigma }}_{\mathrm{s}}.\hfill \end{array}$$(18)

Thus the energy per unit area of the stellar disk is dependent only on the intrinsic parameters of the disk, namely its surface density and the vertical velocity dispersion. This expression is also valid for a gas-alone disk. Now we derive the corresponding analytical expression for the coupled stars-plus-gas disk in a similar way (using Eq. (13)), as given below,

$$\begin{array}{c}\hfill W_{\mathrm{coupled}}=-2{\displaystyle {\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{s}}\mathrm{d}z-2{\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{g}}\mathrm{d}z.}\end{array}$$

Substituting K_z, coupled in terms of hydrostatic equilibrium of each component (Eq. (14)), we obtain

$$\begin{array}{c}\hfill W_{\mathrm{coupled}}=-2{\displaystyle {\int }_0^{\infty }\left[z\frac{{\sigma }_{z,\mathrm{s}}^2}{{\rho }_{\mathrm{s}}}\frac{\mathrm{d}{\rho }_{\mathrm{s}}}{\mathrm{d}z}{\rho }_{\mathrm{s}}\mathrm{d}z\right]-2{\int }_0^{\infty }\left[z\frac{{\sigma }_{z,\mathrm{g}}^2}{{\rho }_{\mathrm{g}}}\frac{\mathrm{d}{\rho }_{\mathrm{g}}}{\mathrm{d}z}{\rho }_{\mathrm{g}}\mathrm{d}z\right].}\end{array}$$

Assuming isothermal dispersion, we obtain

$$\begin{array}{cc}\hfill W_{\mathrm{coupled}}& =-2{\sigma }_{z,\mathrm{s}}^2{\int }_0^{\infty }\left[z\frac{\mathrm{d}{\rho }_{\mathrm{s}}}{\mathrm{d}z}\mathrm{d}z\right]-2{\sigma }_{z,\mathrm{g}}^2{\int }_0^{\infty }\left[z\frac{\mathrm{d}{\rho }_{\mathrm{g}}}{\mathrm{d}z}\mathrm{d}z\right].\hfill \end{array}$$(19)

This is similar to what is obtained for the single-component case. Applying integration by parts to each of the integrations and applying the physical argument at the large z limit (as discussed above for the single component case), we obtain W_coupled to be

$$\begin{array}{c}\hfill W_{\mathrm{coupled}}={\sigma }_{z,\mathrm{s}}^2{\mathrm{\Sigma }}_{\mathrm{s}}+{\sigma }_{z,\mathrm{g}}^2{\mathrm{\Sigma }}_{\mathrm{g}}.\end{array}$$(20)

This shows that the potential energy per unit area of each component (stars or gas) only depends on its intrinsic parameters, that is, the surface density and the vertical velocity dispersion, even within the coupled system. Consequently, the energy of any component remains the same in the coupled case as in the single-component case. After obtaining the general expression of energy of the multi-component system rigorously as in Eq. (13), we could simplify it analytically in a straightforward way, and thus could explain the constancy of the energy that was obtained numerically in Sect. 3.1. In a similar way, the energy of any component in the three-component coupled system of stars and two gas components will also remain the same, as we show in Sect. 3.4.

Based on these results, we can argue that physically, due to the joint gravity of the stars-plus-gas disk, the vertical distribution of stars and gas are now constrained toward the mid-plane. Thus due to the higher vertical force in the coupled case, the self-consistent distribution of each component is now effectively extended to a smaller vertical height, so as to conserve the energy per unit area. We also note that the joint gravity here works like an internal force within the system, and therefore it can just redistribute the energy within each of the two components without changing the total value of the energy.

However, instead of comparing the work done needed to build up the complete vertical mass distribution of the disk, we can compare the work done required to take only a unit, test mass from the z = 0 plane to a certain finite height, discussed in the following section. This could verify whether the disk is more likely to resist distortion for the constrained distribution resulting in a coupled case.

3.3. Work done to raise a unit test mass from the mid-plane to a finite height

We derived the expression of the work done or the energy required to raise a unit test mass from the mid-plane to a certain height h in a single-component (stars or gas) disk and in the coupled two-component stars-plus-gas disk. The work has to be done against the self-gravity of stars (gas) and the joint gravity of stars plus gas disk, respectively. We note that this is precisely the measure of the gravitational potential at any height (Bahcall 1984a,b) in these cases.

For a single-component disk (of stars or gas), this work done is given as

$$\begin{array}{c}\hfill E_{z,i}=-{\displaystyle {\int }_0^h{K}_{z,i}\mathrm{d}z,}\end{array}$$(21)

where K_z, i represents the self-gravity of the disk.

In the gravitationally coupled disk of stars plus gas, both stellar and gas distribution are subject to the same coupled force. Therefore, the work done on a unit mass of the stellar or the gas distribution is

$$\begin{array}{c}\hfill E_{z,i,\mathrm{coupled}}=-{\displaystyle {\int }_0^h{K}_{z,\mathrm{coupled}}\mathrm{d}z.}\end{array}$$(22)

We note that in each case the energy is positive. Due to a higher vertical force at each z, as was shown in Fig. 2, the work done to take the unit mass of the stellar (gas) disk to the same height in the coupled case will be higher than in the corresponding single-component case. We show the work done in these cases as a function of vertical height in Fig. 4.

Fig. 4. Work done (Ez) to raise a unit test mass from the mid-plane to any vertical height as a function of z (shown upto z = 1 kpc here) at R = 8.5 kpc (solar neighbourhood) for the stars-alone case against its self-gravity (solid curve), for the gas-alone case against its self-gravity (dashed curve), and for stars or gas in the stars-plus-gas coupled system against their coupled gravity (dash-dotted curve). The work done required to raise a unit mass at any vertical height in the last case is highest. This shows that stars or gas in the coupled system are more strongly bound to the mid-plane than in the corresponding single-component cases.

This shows that to raise a unit mass from the mid-plane to a certain vertical distance, more work is required in the coupled case than that for a single-component case. Thus stars are more strongly bound to the mid-plane of the Galaxy in the coupled case than in the single-component case, and thus the stellar disk will be able to offer more resistance to a given external tidal encounter. In this case, the stellar disk is therefore less likely to be thickened (Walker et al. 1996). Moreover, due to the constrained distribution, the stellar mass distribution is more concentrated towards the mid-plane. This increases the effective gravity of the stellar disk near the mid-plane and helps it to resist external perturbations, which could have led to the generation of warps (Pranav & Jog 2010). Thus the constraining effect of gas on stars makes the stellar disk less likely to be disturbed. A detailed N-body simulation will be able to show this clearly. This is beyond the scope of this paper.

Similarly, a unit mass of the gas disk constrained by the stellar gravity is more strongly bound in the coupled system than in the gas-alone case, and hence significantly more work is required to raise a unit mass of the gas disk to a certain height than in the gas-alone case. Thus gas disk in the coupled system is less likely to be disturbed than the gas-alone case. Furthermore, we note that far more work is required to raise gas to a certain height in the coupled case compared to gas-alone case than the corresponding work required for stars because stars are the more massive component and have a stronger effect on gas.

3.4. Calculation of the potential energy of a three-component disk

For the sake of completeness, we next studied the gravitational potential energy per unit area of a gravitationally coupled three-component galactic disk, consisting of stars and two gas components. The three components were taken to be coplanar with the same mid-plane at z = 0. For illustration, we added H₂ as the second gas component in addition to HI, as seen in the inner Galaxy. We note that this three-component treatment is essential in the inner Galaxy.

For the three-component system of stars, HI, and H₂, the gravitational potential energy per unit area of the disk is the work done to build a column of gravitationally coupled stars, HI, H₂, of unit cross-section, together from z = 0. The Poisson equation for this system is given as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}{z}^2}=4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{HI}}+{\rho }_{{\mathrm{H}}_2}).\end{array}$$(23)

Following the same procedure as in Sect. 2.2, we derived the potential energy per unit area of the three-component disk as

$$\begin{array}{c}\hfill W_{\mathrm{coupled}}={\displaystyle {\int }_{-\infty }^{\infty }z\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{coupled}}}{\mathrm{d}z}({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{HI}}+{\rho }_{{\mathrm{H}}_2})\mathrm{d}z,}\end{array}$$(24)

which can be further expressed as

$$\begin{array}{cc}\hfill W_{\mathrm{coupled}}& =-2{\displaystyle {\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{s}}\mathrm{d}z-2{\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{\mathrm{HI}}\mathrm{d}z}\hfill \\ \hfill & -2{\displaystyle {\int }_0^{\infty }z{K}_{z,\mathrm{coupled}}{\rho }_{{\mathrm{H}}_2}\mathrm{d}z.}\hfill \end{array}$$(25)

These three integrations can be considered to represent the potential energy per unit area of the stellar disk, HI disk, and H₂ disk in the three-component coupled system.

Now for a three-component system, the hydrostatic balance of each component is determined by the joint gravitational force from stars, HI, and H₂, and is given by

$$\begin{array}{c}\hfill \frac{{\sigma }_{z,i}^2}{{\rho }_i}\frac{\mathrm{d}{\rho }_i}{\mathrm{d}z}=K_{z,s}+K_{z,\mathrm{HI}}+K_{z,{\mathrm{H}}_2},\end{array}$$(26)

where i represents stars (s) or HI or H₂, and the right-hand side of the equation represents the vertical force of the coupled system. Combining this equation with Eq. (23), we write the joint hydrostatic balance-Poisson equation for the three-component system as

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\rho }_i}{\mathrm{d}{z}^2}=\frac{{\rho }_i}{{\sigma }_{z,i}^2}[-4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{HI}}+{\rho }_{{\mathrm{H}}_2})]+\frac{1}{{\rho }_i}{\left(\frac{\mathrm{d}{\rho }_i}{\mathrm{d}z}\right)}^2.\end{array}$$(27)

These coupled equations were solved to obtain ρ_i vs. z for stars, HI, and H₂ following the same method as discussed in Sect. 2.3.1. Here, we discuss the results for the stellar distribution alone for simplicity because the results for the gas components follow a similar trend as for stars, as seen for the two-component case in Sects. 3.1 and 3.3. We chose R = 4.5 kpc to illustrate the results here. The stellar surface density at R = 4.5 kpc was calculated to be 157.06 M_⊙ pc⁻² (Mera et al. 1998). The vertical velocity dispersion of stars was obtained to be 28.2 km s⁻¹ (Lewis & Freeman 1989) and assuming the vertical-to-radial dispersion ratio to be the same as in the solar neighbourhood, see Sect. 2.3.2. We chose this radius to prominently show the effect of the third component, namely H₂, on stars. The surface density of H₂ is 19.7 M_⊙ pc⁻² at this radius, which is significantly higher than that of HI, which is 4.6 M_⊙ pc⁻² (Scoville & Sanders 1987). The vertical velocity dispersion of HI is 8 km s⁻¹, as it was in the two-component case (Sect. 2.3.2), and that of H₂ is 5 km s⁻¹ at this radius (Scoville & Sanders 1987).

We show the results for stellar distribution in Fig. 5 in the three-component system, in the two-component system (stars plus HI), and in the stars-alone case. Figure 5a shows that the addition of a second gas component (H₂ here) constrains the stellar distribution towards the mid-plane by raising the mid-plane density value and reducing the disk thickness value compared to the values in the two-component case. We note that due to a higher surface density of H₂, the constraining effect due to this component on stars is more prominent than that due to HI. Figure 5b shows that the coupled gravitational force per unit mass (calculated numerically using Eq. (26)) that keeps the stellar distribution in the hydrostatic equilibrium in the three-component system is higher than that of the two-component system and the stars-alone case at all z.

Fig. 5. Results for stellar distribution in the three-component, gravitationally coupled stars plus HI plus H2 disk, two-component stars plus HI disk, and stars-alone case at R = 4.5 kpc. Panel a: stellar vertical density distribution, ρ(z) vs. z is shown for stars-alone case (solid curve), gravitationally coupled stars plus HI case (dashed curve), and stars plus HI plus H2 case (dash-dotted curve). The gas gravity constrains the distribution of stars towards the mid-plane by raising its mid-plane density and reducing its scale height value. This constraining effect is highest in the three-component case. At the same time, the relative contribution of H2 is higher than that of HI as the surface density of H2 is higher than HI at R = 4.5 kpc. Panel b: vertical force per unit mass, i.e. |Kz| vs. z, acting on stars is shown due to the self-gravity of stars (solid curve), due to the force from the gravitationally coupled stars plus HI disk (dashed curve), and due to the force from the coupled stars plus HI plus H2 disk (dash-dotted curve). The force at any z is highest in the three-component case. Panel c: work done (Ez) to raise a unit test mass from the mid-plane to any vertical height as a function of z (shown up to z = 1 kpc here) for stars against its self-gravity (solid curve), against the gravitational force from the coupled stars plus HI disk (dashed curve), and against the gravitational force from the coupled stars plus HI plus H2 disk (dash-dotted curve). The work done in the three-component case is highest. This shows that stellar distribution in a three-component system is more strongly bound to the mid-plane than in a two-component system.

Despite the higher constraining effect in the three-component system, the potential energy per unit area of the stellar distribution, calculated numerically from Eq. (25), is found to be 124 900.4 M_⊙ pc⁻² km² s⁻² (2.6 × 10⁸ erg cm⁻²), which is the same as that found for stars in the two-component system and in the stars-alone system at R = 4.5 kpc. This result is expected because we find the analytical expression of the total energy per unit area of the three-component coupled disk to be $W_{\mathrm{coupled}}={\sigma }_{z,s}^2{\mathrm{\Sigma }}_{\mathrm{s}}+{\sigma }_{z,\mathrm{HI}}^2{\mathrm{\Sigma }}_{\mathrm{HI}}+{\sigma }_{z,{\mathrm{H}}_2}^2{\mathrm{\Sigma }}_{{\mathrm{H}}_2}$, derived following a similar method as discussed for the two-component system in Sect. 3.2. Thus, we note that the energy of each component remains unchanged.

However, when we calculate the work required to take a unit test mass of stellar distribution from the mid-plane to a certain height h in the above three cases following the method discussed in Sect. 3.3, we find that the work done is highest in the three-component case at any z. We show the work done as a function of z corresponding to the three cases in Fig. 5c. This shows that stars are more strongly bound to the mid-plane of the Galaxy in the three-component case than in the two-component system, and thus the stellar disk is able to offer more resistance to a given external tidal encounter. Thus the constraining effect of the two gas components will make the stellar disk less likely to be disturbed due to external perturbations. Although we do not show results for gas here for conciseness, the stellar component, being more massive, has a higher effect on making the gas disk less likely to be disturbed due to perturbations than the gas-alone case (as was already seen for the two-component case in Sects. 3.1 and 3.3).

4. Discussion

We discuss a few general implications of the model developed in this paper below. First, the main aim of this paper was to determine how the gravitational potential energy per unit area of the disk components changes in view of the constraining effect in the coupled case. The energy of the components turned out to be the same as in their single-component self-gravitating cases, which is in contrast to our initial expectation. Therefore, the question is whether stars and gas are more strongly bound in the coupled case. Our results have shown that it requires a higher amount of energy to raise the unit mass of a component to a certain vertical height in the multi-component case than in its single component case, even though the component contains the same potential energy per unit area in both cases. Thus each component is more strongly bound to the mid-plane in the multi-component system.

Second, we note that the potential energy per unit area of any disk component is dependent on its surface mass density and its vertical velocity dispersion. This implies that the magnitude of the energy depends on the component chosen at a given galactocentric radius and also on the galactocentric radius for a given component, as the above parameters vary along radius. We note that at the solar radius as well as at R = 4.5 kpc, the surface density and the vertical velocity dispersion of stars are so much higher than gas that the energy value per unit area for the stellar disk is much higher than gas. The stellar disk is more extended vertically than the gas disk, and hence we would expect it to be more disturbed by a given external tidal encounter, for instance from a passing satellite galaxy. However, the stellar disk has a higher potential energy per unit area. Although it is more extended, it is therefore less likely to show the effect of a given tidal disturbance. This trend can be confirmed by numerical simulations of an encounter. This is beyond the scope of this paper.

We also note that we have considered only gravitational interaction for both stellar and gas disks. We did not consider any gas dynamical phenomenon.

Third, in the above cases, we have taken two (three) components with different dispersions to identify them as stars and one (two) gas components. It may be an interesting physical question to ask what happens when the stellar disk is divided artificially into n number of components. Following a similar analysis, we find again that the total energy of the components remain unchanged, as is the case for the components with different dispersions. Interestingly, in this case, the net distribution is not vertically more constrained, as expected physically. For example, following the numerical analysis in Sect. 2.3.1, we checked that if we were to divide the disk of Σ into two sub-components of Σ₁ and Σ₂ with the same dispersions, the net distribution in the coupled case would be identical to the one-component case with Σ.

Fourth, the treatment given in the paper is general. Although we have applied it to stars and gas case, it can be applied to n number of stellar sub-components in the Galactic disk as well. Such sub-components have been identified from recent observed data, for example from Gaia (Bovy 2017; Hagen & Helmi 2018; etc.). If the surface density, dispersions are known for these components, then the potential energy per unit area of these components can be determined quantitatively.

5. Conclusions

It has been shown earlier that in a multi-component gravitationally coupled stars-plus-gas disk, the self-consistent vertical distribution of stars is constrained closer to the mid-plane (Sarkar & Jog 2018). In order to understand the implications of this for the energetics of the disk, we obtained the potential energy per unit area at a given galactocentric radius for a multi-component galactic disk. This was obtained as an integration over the vertical density (ρ(z)), the gravitational force, and vertical distance z. To do this, we followed the method developed by Camm (1967) for a single-component self-gravitating disk and explicitly derived the corresponding expression for the multi-component case. For a self-consistent distribution we obtained the density distribution and force as a function of z by numerically solving the joint hydrostatic balance and Poisson equation for the coupled case.

1. We find that the net gravitational potential energy for the stars and gas remain unchanged to that in the single-component cases. This is a surprising result and can be understood by simplifying analytically the general expression for the potential energy that we obtain for the multi-component system. We noted that the potential energy per unit area of each component depends only on its intrinsic parameters, that is, the vertical velocity dispersion and the surface density, in the single-component as well as in the multi-component case. Physically, the energy values remain unchanged because due to the higher joint gravity in the coupled case, the distribution of each component is constrained closer to the mid-plane and thus has a less effective vertical thickness to conserve the energy per unit area. We note that the joint gravity works here like an internal force within the system, and therefore it can just redistribute the energy within each component itself without changing the total value of energy for each component.

2. However, due to the constrained distribution in the coupled cases, the work required to raise a unit test mass to a given height is higher than that in the single-component case. Thus, while constraining in a coupled case does not correspond to any additional gravitational energy in the system, it does indicate that each component in a disk in the coupled case is more strongly bound to the mid-plane. Furthermore, the stellar disk has a higher potential energy per unit area than the gas disk. For a given tidal encounter, the stellar disk is therefore less likely to be disturbed than a gas disk, even though it is more vertically extended.

Acknowledgments

We thank the anonymous referee for very useful comments. We would like to thank Tarun Deep Saini for discussions which helped in clarifying our ideas for developing the problem. S.S. thanks CSIR for a fellowship, and C.J. thanks the DST for support via J.C. Bose fellowship (SB/S2/JCB-31/2014).

References

All Tables

All Figures

Fig. 1. Stellar and gas (HI) vertical density distribution, ρ(z) vs. z at R = 8.5 kpc (panels a and b, respectively). In each plot, the solid curve represents the vertical distribution in the one-component case. The dashed curve represents the distribution obtained in the gravitationally coupled stars-plus-gas system. The coupled gravitational force of the star- plus-gas system constrains each distribution towards the mid-plane, i.e., increases the mid-plane density value of the distribution, which causes the scale height to become smaller and the curve to become steeper.

In the text

	Fig. 2. Vertical force per unit mass, i.e., |Kz| vs. z at R = 8.5 kpc. The solid and dashed curves represent the force due to the self-gravity of stars and gas (HI), respectively. The dash-dotted curve represents the force arising from the gravitationally coupled stars-plus-gas system, which is higher than the forces due to the self-gravity of the individual components at each z.
In the text

Fig. 3. Integrand −zρ(z)Kz vs. z (from Eqs. (2), (13)) of W (potential energy per unit area) for stars and gas (HI) at R = 8.5 kpc. Panela: energy integrand for stars in the stars-alone case (solid curve) and in the coupled stars-plus-gas case (dashed curve). Twice the area under the curves gives the energy per unit area of the stellar disk in the two cases, which is obtained to be the same. Panel b: energy integrand for gas (HI) in the gas-alone case (solid curve) and in the coupled stars-plus-gas case (dashed curve). Twice the area under the curves gives the energy per unit area of the gas disk in the two cases, which are obtained to be the same.

In the text

Fig. 4. Work done (Ez) to raise a unit test mass from the mid-plane to any vertical height as a function of z (shown upto z = 1 kpc here) at R = 8.5 kpc (solar neighbourhood) for the stars-alone case against its self-gravity (solid curve), for the gas-alone case against its self-gravity (dashed curve), and for stars or gas in the stars-plus-gas coupled system against their coupled gravity (dash-dotted curve). The work done required to raise a unit mass at any vertical height in the last case is highest. This shows that stars or gas in the coupled system are more strongly bound to the mid-plane than in the corresponding single-component cases.

In the text

Fig. 5. Results for stellar distribution in the three-component, gravitationally coupled stars plus HI plus H2 disk, two-component stars plus HI disk, and stars-alone case at R = 4.5 kpc. Panel a: stellar vertical density distribution, ρ(z) vs. z is shown for stars-alone case (solid curve), gravitationally coupled stars plus HI case (dashed curve), and stars plus HI plus H2 case (dash-dotted curve). The gas gravity constrains the distribution of stars towards the mid-plane by raising its mid-plane density and reducing its scale height value. This constraining effect is highest in the three-component case. At the same time, the relative contribution of H2 is higher than that of HI as the surface density of H2 is higher than HI at R = 4.5 kpc. Panel b: vertical force per unit mass, i.e. |Kz| vs. z, acting on stars is shown due to the self-gravity of stars (solid curve), due to the force from the gravitationally coupled stars plus HI disk (dashed curve), and due to the force from the coupled stars plus HI plus H2 disk (dash-dotted curve). The force at any z is highest in the three-component case. Panel c: work done (Ez) to raise a unit test mass from the mid-plane to any vertical height as a function of z (shown up to z = 1 kpc here) for stars against its self-gravity (solid curve), against the gravitational force from the coupled stars plus HI disk (dashed curve), and against the gravitational force from the coupled stars plus HI plus H2 disk (dash-dotted curve). The work done in the three-component case is highest. This shows that stellar distribution in a three-component system is more strongly bound to the mid-plane than in a two-component system.

In the text