MA Emilio Schmidt

The isothermal self-gravitating cylinder

Deriving the ODE

Sought is the stationary density profile \(\rho\) of an infinitely long, isothermal, self-gravitating gas cylinder. To do this, the hydrostatic Euler equation has to be solved, which is given by:

\begin{align} \boldsymbol{\nabla} P = \boldsymbol{f}_\text{ext} \end{align}

Here, \(P\) is the pressure field and \(\boldsymbol{f}_\text{ext}\) is the force density acting on the gas. In the considered case, this is the gas's own gravity. Since the gravitational field is a conservative force field, it can be expressed by a potential \(\phi_\text{ext}\):

\begin{align} \boldsymbol{f}_\text{ext} = - \rho \boldsymbol{\nabla} \phi_\text{ext} \end{align}

Inserting yields:

\begin{align} \boldsymbol{\nabla} P = - \rho \boldsymbol{\nabla} \phi_\text{ext} \end{align}

The pressure field of the gas under consideration can be described by the equation of state of an isothermal gas, which is given by:

\begin{align} P = \frac{R T}{M} \rho \end{align}

Here, \(R\) denotes the universal gas constant, \(T\) the spatially homogeneous temperature distribution, and \(M\) the molar mass of the gas. Inserting into the hydrostatic Euler equation yields:

\begin{align} \phantom{\Rightarrow}\;& \frac{R T}{M} \boldsymbol{\nabla} \rho = - \rho \boldsymbol{\nabla} \phi_\text{ext} \\[6pt] \Leftrightarrow\;& \frac{1}{\rho} \boldsymbol{\nabla} \rho = - \frac{M}{R T} \boldsymbol{\nabla} \phi_\text{ext} \\[6pt] \Leftrightarrow\;& \boldsymbol{\nabla} \ln(\rho) = - \frac{M}{R T} \boldsymbol{\nabla} \phi_\text{ext} \\[6pt] \Rightarrow\;& \Delta \ln(\rho) = - \frac{M}{R T} \Delta \phi_\text{ext} \end{align}

It proves useful to use the Poisson equation of Newtonian gravity, which is given by:

\begin{align} \Delta \phi_\text{ext} = 4 \pi G \rho \end{align}

Inserting yields:

\begin{align} \phantom{\Rightarrow}\;& \Delta \ln(\rho) = - \frac{4 \pi G M}{R T} \rho \, , \quad \text{def.:} \, a := \frac{4 \pi G M}{R T} \\[6pt] \Leftrightarrow\;& \Delta \ln(\rho) = - a \rho \end{align}

For an infinitely long, straight, self-gravitating gas cylinder without an externally imposed azimuthal or axial structure, there is symmetry with respect to rotation about the cylinder axis and translation along the axis. The underlying equations share this symmetry, which is why corresponding derivatives disappear in the considered equation:

\begin{align} \phantom{\Rightarrow}\;& \frac{1}{r} \partial_r \big(r \partial_r \ln(\rho)\big) = - a \rho \\[6pt] \Leftrightarrow\;& \partial_r \bigg(r \frac{\partial_r \rho}{\rho}\bigg) = - a r \rho \\[6pt] \Leftrightarrow\;& \frac{\partial_r \rho}{\rho} - r \bigg(\frac{\partial_r \rho}{\rho}\bigg)^2 + \frac{\partial_r^2 \rho}{\rho} = - a r \rho \end{align}

Solving the ODE

Due to the symmetry, the quantities considered depend only on the radius \(r\). Therefore, in the differential equation considered, the partial derivatives can be replaced by total derivatives:

\begin{align} \frac{\text{d}_r \rho}{\rho} - r \bigg(\frac{\text{d}_r \rho}{\rho}\bigg)^2 + r\frac{\text{d}_r^2 \rho}{\rho} = - a r \rho \end{align}

Since the density is positive by definition, the following substitution can be made without restriction:

\begin{align} u(r) = \ln\left(\frac{\rho(r)}{\rho_0}\right). \end{align}

Here, \(\rho_0 = \text{const.} > 0\) is a constant reference density. Its specific value is irrelevant for further consideration. It only serves to ensure that the argument of the logarithm is dimensionless. The first two derivatives of the new quantity \(u\) are then given by:

\begin{align} \text{d}_r u = \frac{\text{d}_r \rho}{\rho} \quad \text{und} \quad \text{d}_r^2 u = \frac{\text{d}_r^2 \rho}{\rho} - \left(\frac{\text{d}_r \rho}{\rho}\right)^2 \end{align}

Inserting this into the differential equation under consideration yields:

\begin{align} \frac{\text{d}_r u}{r} + \text{d}_r^2 u = - a \rho_0 \exp(u) \end{align}

To simplify this equation further, another substitution is performed:

\begin{align} z = u + \sqrt{2} t \, , \quad \text{with} \, t = \sqrt{2} \ln\left(\sqrt{a \rho_0} r\right) \end{align}

With \(\text{d}_r t = \frac{\sqrt{2}}{r}\), the following applies for the first two derivatives of \(u\) with respect to \(r\):

\begin{align} \text{d}_r u = \frac{\sqrt{2}}{r} \text{d}_t u \quad \text{and} \quad \text{d}_r^2 u = - \frac{\sqrt{2}}{r^2} \text{d}_t u + \frac{2}{r^2} \text{d}_t^2 u \end{align}

Inserting yields:

\begin{align} \phantom{\Rightarrow}\;& 2 \text{d}_t^2 u = - a \rho_0 r^2 \exp(u) \\[6pt] \Leftrightarrow\;& 2 \text{d}_t^2 u = - \exp(z) \end{align}

Furthermore, the following applies:

\begin{align} \phantom{\Rightarrow}\;& \text{d}_t^2 z = \text{d}_t^2 u + \sqrt{2} \text{d}_t^2 t \\[6pt] \Leftrightarrow\;& \text{d}_t^2 z = \text{d}_t^2 u \end{align}

Thus, it follows that:

\begin{align} \phantom{\Rightarrow}\;& 2 \text{d}_t^2 z = - \exp(z) \\[6pt] \Rightarrow\;& 2 \left(\text{d}_t z\right) \left(\text{d}_t^2 z\right) = - \left(\text{d}_t z\right) \exp(z) \\[6pt] \Leftrightarrow\;& \text{d}_t\left(\left(\text{d}_t z\right)^2 + \exp(z)\right) = 0 \\[6pt] \Rightarrow\;& \left(\text{d}_t z\right)^2 + \exp(z) = c_1 \, , \quad c_1 \in \mathbb{R} \\[6pt] \Leftrightarrow\;& \text{d}_t z = \pm \sqrt{c_1 - \exp(z)} \end{align}

At first, it might seem that both signs need to be treated separately. However, it can actually be shown that both equations lead to the same solutions. In the following, therefore, only the equation with a positive sign will be considered. To solve this, the following substitution is used:

\begin{align} z = \ln(y) \end{align}

This results in the following for the differential equation under consideration:

\begin{align} \frac{\text{d}_t y}{y} = \sqrt{c_1 - y} \end{align}

Separation of variables yields:

\begin{align} \phantom{\Rightarrow}\;& \frac{\text{d}y}{y\sqrt{c_1 - y}} = \text{d}t \\[6pt] \Rightarrow\;& \frac{1}{\sqrt{c_1}} \ln\left(\bigg|\frac{\sqrt{c_1 - y} - \sqrt{c_1}}{\sqrt{c_1 - y} + \sqrt{c_1}}\bigg|\right) = t + C_2 \\[6pt] \Leftrightarrow\;& \bigg|\frac{\sqrt{c_1 - y} - \sqrt{c_1}}{\sqrt{c_1 - y} + \sqrt{c_1}}\bigg| = \exp\left(\sqrt{c_1} t + \sqrt{c_1} C_2\right) \, , \quad \text{with} \, c_2 = \exp\left(\sqrt{c_1} C_2\right) \\[6pt] \Leftrightarrow\;& \bigg|\frac{\sqrt{c_1 - y} - \sqrt{c_1}}{\sqrt{c_1 - y} + \sqrt{c_1}}\bigg| = c_2 \exp\left(\sqrt{c_1} t\right) \, , \quad \text{with} \, A = c_2 \exp\left(\sqrt{c_1} t\right) \\[6pt] \Leftrightarrow\;& \bigg|\frac{\sqrt{c_1 - y} - \sqrt{c_1}}{\sqrt{c_1 - y} + \sqrt{c_1}}\bigg| = A \, , \quad \text{with} \, y = \exp(z) \Rightarrow \sqrt{c_1 - y} > \sqrt{c_1} \; \forall y \\[6pt] \Leftrightarrow\;& \frac{\sqrt{c_1} - \sqrt{c_1 - y}}{\sqrt{c_1} + \sqrt{c_1 - y}} = A \\[6pt] \Leftrightarrow\;& y = 4 c_1 \frac{A}{(1 + A)^2} \\[6pt] \Leftrightarrow\;& \exp(u) = \frac{4 c_1}{a \rho_0 r^2} \frac{A}{(1 + A)^2} \\[6pt] \Leftrightarrow\;& \rho = \frac{4 c_1 \rho_0}{a \rho_0 r^2} \frac{A}{(1 + A)^2} \end{align}

The quantity \(A\) can be written as follows:

\begin{align} \phantom{\Rightarrow}\;& A = c_2 \exp\left(\sqrt{c_1} t\right) \\[6pt] \Leftrightarrow\;& A = c_2 \exp\left(\sqrt{2 c_1} \ln(\sqrt{a \rho_0} r)\right) \\[6pt] \Leftrightarrow\;& A = c_2 (\sqrt{a \rho_0} r)^{\sqrt{2 c_1}} \end{align}

Define additionally:

\begin{align} b := \sqrt{2 c_1} \quad \text{und} \quad r_0 := \left(c_2^\frac{1}{b} \sqrt{a \rho_0}\right)^{-1} \end{align}

This yields the following general solution for the differential equation under consideration:

\begin{align} \phantom{\Rightarrow}\;& \rho(r) = \frac{2 b^2}{a r^2} \frac{\left(\frac{r}{r_0}\right)^b} {\left(1 + \left(\frac{r}{r_0}\right)^b\right)^2} \\[6pt] \Leftrightarrow\;& \rho(r) = \frac{2 b^2}{a r_0^2} \frac{\left(\frac{r}{r_0}\right)^{b-2}} {\left(1 + \left(\frac{r}{r_0}\right)^b\right)^2} \end{align}

Boundary conditions

Gas only

The integration constants \(r_0\) and \(b\) must then be determined using appropriate boundary conditions. For the problem considered, it makes sense to specify boundary conditions for the center of the cylinder, i.e., for \(r = 0\). In the considered case, the density at the center of the cylinder must be finite and nonzero, since both a vanishing and a diverging density at the origin would lead to unphysical behavior:

\begin{align} \rho(r=0) = \rho_c \, , \quad \text{with} \quad 0 < \rho_c < \infty \end{align}

Furthermore, rotational symmetry and the requirement for a smooth, regular solution imply that the density profile at the origin should not have a gradient. Consequently, \(\rho(r)\) has an extremum at \(r=0\), and the following applies:

\begin{align} \text{d}_r \rho(r) \big|_{r = 0} = 0 \end{align}

It should be noted that the substitution \(t=\sqrt{2}\ln\big(\sqrt{a\rho_0}r\big)\) is only defined for \(r>0\). The differential equation was thus initially solved in the domain \(r \in \mathbb{R}\). In order to formulate boundary conditions at \(r=0\), the solution must be extended continuously to \(r=0\). Accordingly, the central values are understood as right-sided limits:

\begin{align} \rho(r = 0) = \lim\limits_{r \rightarrow 0} \rho(r) \quad \text{and} \quad \text{d}_r \rho(r) \big|_{r = 0} = \lim\limits_{r \rightarrow 0} \text{d}_r \rho(r) \end{align}

With the first boundary condition, the following then applies:

\begin{align} \phantom{\Rightarrow}\;& \lim\limits_{r \rightarrow 0} \rho(r) = \rho_c \\[6pt] \Leftrightarrow\;& \lim\limits_{r \rightarrow 0} \frac{2 b^2}{a r_0^2} \frac{\left(\frac{r}{r_0}\right)^{b-2}} {\left(1 + \left(\frac{r}{r_0}\right)^b\right)^2} = \rho_c \end{align}

Since the denominator converges to \(1\) in the limit case \(r\rightarrow 0\), the local behavior of the density in the immediate vicinity of the origin is determined exclusively by the numerator. This, together with the fact that the density at the origin must not vanish or diverge, directly implies that the following must hold:

\begin{align} b = 2 \end{align}

Inserting into the solution for the density profile provides:

\begin{align} \rho(r) = \frac{8}{a r_0^2} \frac{1} {\left(1 + \left(\frac{r}{r_0}\right)^2\right)^2} \end{align}

Furthermore, the first boundary condition implies that:

\begin{align} \phantom{\Rightarrow}\;& \rho_c = \frac{8}{a r_0^2} \\[6pt] \Leftrightarrow\;& r_0^2 = \frac{8}{a \rho_c} \end{align}

Inserting this gives the stationary density profile of an infinitely long, isothermal, self-gravitating gas cylinder:

\begin{align} \rho(r) = \frac{\rho_c}{\left(1 + \frac{1}{8} a \rho_c r^2\right)^2} \end{align}

It is noticeable that the second boundary condition is not explicitly included in the determination of the integration constants \(r_0\) and \(b\). The reason for this is that the requirement for a finite central density different from zero already restricts the solution so much that only the case \(b=2\) remains. Thus, the density profile takes the above form and automatically satisfies the second boundary condition. It is therefore not redundant, but is already implicitly contained in the regularity and symmetry requirements.

With the determined density distribution, the gravitational potential induced by the density distribution can now be determined using Poisson's equation. The Poisson equation, using the given symmetry, is given by:

\begin{align} \frac{1}{r} \partial_r \left(r \partial_r \phi_\text{ext}(r)\right) = 4 \pi G \rho(r) \end{align}

Also here, due to the symmetry, the partial derivatives can be replaced by total derivatives:

\begin{align} \phantom{\Rightarrow}\;& \frac{1}{r} \text{d}_r \left(r \text{d}_r \phi_\text{ext}(r)\right) = 4 \pi G \rho(r) \\[6pt] \Leftrightarrow\;& \text{d}_r \left(r \text{d}_r \phi_\text{ext}(r)\right) = 4 \pi G \rho(r) r \end{align}

Integration yields:

\begin{align} \phantom{\Rightarrow}\;& r \text{d}_r \phi_\text{ext}(r) = 4 \pi G \int_{0}^{r} \rho(r^\prime) r^\prime \text{d} r \\[6pt] \Leftrightarrow\;& r \text{d}_r \phi_\text{ext}(r) = \frac{32 \pi G}{a r_0^2} \int_{0}^{r} \frac{r^\prime}{\left(1 + \left(\frac{r^\prime}{r_0}\right)^2\right)^2} \text{d} r \\[6pt] \Leftrightarrow\;& r \text{d}_r \phi_\text{ext}(r) = \frac{16 \pi G}{a} \frac{\left(\frac{r}{r_0}\right)^2} {1 + \left(\frac{r}{r_0}\right)^2} \\[6pt] \Leftrightarrow\;& \text{d}_r \phi_\text{ext}(r) = \frac{16 \pi G}{a r_0} \frac{\frac{r}{r_0}}{1 + \left(\frac{r}{r_0}\right)^2} \end{align}

Another integration yields the gravitational potential of of an infinitely long, isothermal, self-gravitating gas cylinder:

\begin{align} \phantom{\Rightarrow}\;& \phi_\text{ext}(r) = \frac{16 \pi G}{a r_0} \int_{0}^{r} \frac{\frac{r^\prime}{r_0}}{1 + \left(\frac{r^\prime}{r_0}\right)^2} \text{d} r \\[6pt] \Leftrightarrow\;& \phi_\text{ext}(r) = \frac{8 \pi G}{a} \ln\left(1 + \left(\frac{r}{r_0}\right)^2\right) \end{align}

The case with a core

The solution to the ODE reads

\begin{align} \rho &= \frac{\lambda b^2(r/r_0)^{b-2}}{4\pi r_0^2(1+(r/r_0)^b)^2}\\[2ex] \text{with}\quad \lambda &= \frac{2k_\text{B}T}{Gm} = \int_0^\infty 2\pi r\rho(r,r_0,b{=}2)\,\text dr \end{align} being the only possible total mass/length in the case of gas only, where \(b=2\) is necessary for zero force at \(r=0\). Notably, \(\lambda\) does not depend on \(r_0\), i.e. every width of \(\rho\) is equally allowed. The reason is that no length scale can be constructed from the input parameters \(k_\text{B}T\), \(G\) and \(m\).

This changes in the presence of a core with mass/length \(\lambda_\text{c}\) and radius \(R\). Vanishing force at \(r=0\) is replaced by force balance at \(r=R\), which leads to the condition \begin{align} \frac{r_0}{R} &= \left(\frac{b+2-4\lambda_\text{c}/\lambda}{b-2+4\lambda_\text{c}/\lambda}\right)^{1/b} ~, \end{align} which in turn introduces the constraint \(b>\vert 4\lambda_\text{c}/\lambda-2\vert\).

The resulting gas mass is \begin{gather} \int_R^\infty 2\pi r\rho(r,r_0,b)\,\text dr = \lambda\frac{b+2}{4}-\lambda_\text{c}\\ \Leftrightarrow\quad \frac{\lambda_\text{c}+\lambda_\text{gas}}{\lambda}=\frac{b+2}{4} ~. \end{gather} That is, the amount of gas can be tuned via \(b\). Note that for \(\lambda_\text{c}<\lambda/2\), it cannot be chosen arbitrarily small.

If the core is actually a numerical necessity and only meant to emulate the gas contained inside, it's \(b=2\) again, and the most convenient choice is \(\lambda_\text{c}=\lambda_\text{gas}=\lambda/2\) which yields \(r_0=R\).