Radiative cooling of gravitationally collapsing molecular cloud

Colliding flows and the birth of stars

Caption_right

Star-forming region S106 (captured by the Hubble Space Telescope)

Getting from a diffuse gas cloud to a star starts with the gravitational collapse of molecular clouds, the cold, dense pockets of interstellar gas where stars form. The collapse is probably set off by colliding flows, turbulent gas streams that crash into one another and compress the material. Molecular clouds are much denser than their surroundings, so some compression must be involved. The question is which forces drive it, and how they shape the evolution of these structures.

Shock waves

In molecular clouds, shock waves can be set off by colliding gas streams or by the cloud moving through the interstellar medium. A shock compresses the gas it passes through, and the densest regions can then overcome their internal pressure and collapse under their own gravity, eventually forming stars.

Self-similarity

Caption_left

Murakami et al. (2004)

Collapsing clouds look chaotic, but they often follow a self-similar pattern. In physics, self-similarity means that a system looks the same at different scales. For example, you would probably be familiar with the Larson-Penston (LP) solution (Larson 1969; Penston 1969) which uses the continuity and momentum conservation equations to give a self-similar solution for the gravitational collapse of an isothermal gaseous sphere.

The behavior of x (distance or radius) and t (time) in these self-similar solutions could be such that the radius and density are related by a power-law expression that remains constant as the collapse progresses. This means that if you were to plot the density profile of the cloud at different times, the shape of the profile would remain the same, but the scale would change according to a power-law relationship with time.

Planar geometry

Caption_right

Plane-parallel geometry

In this project, we use a plane-parallel geometry, where physical properties vary only along the x-axis and are uniform along y and z. In a self-similar flow with this geometry, the velocity, temperature and density profiles keep the same shape as they evolve, only rescaled.

So we will have:

• $v=v_x,\rho ={\rho }_x$
• $\frac{d}{dy}=\frac{d}{dz}=0$

Radiative cooling

One of our main assumptions is that radiative losses balance out compressional heating, keeping the contracting gas relatively cold. As the gas is compressed, the extra internal energy escapes as radiation, so it heats up less than gas under adiabatic compression would.

Objectives

The objective here is to present a self-similar solution that considers two physical phenomena: self-gravity and radiative cooling.

• We address the complete set of hydrodynamic equations,
• Introduce a cooling term into the energy equation to account for heat loss, and search for self-similar solutions.

The physical system

The compressible Navier-Stokes equations with cooling are formulated as follows:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0$$ $$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot (\rho \mathbf{uu})=-\nabla P-\rho \nabla \varphi$$ $$\Delta \varphi =4\pi G\rho$$ $$\frac{\partial \rho (ϵ+\frac{u^2}{2})}{\partial t}+\nabla \cdot \left[\rho \left(ϵ+\frac{u^2}{2}\right)\mathbf{u}\right]=-\nabla \cdot (P\mathbf{u})-\rho \mathbf{u}\cdot \nabla \varphi -\Lambda$$

where t is time, u is the velocity vector, φ is the gravitational potential, G is the gravitational constant, $ϵ$ is the specific internal energy and Λ is the volumic cooling rate. The net heat exchange includes both heating and cooling, but we only consider cases where cooling dominates, so the net effect is a loss of energy. The gas follows the ideal gas law:

$$\frac{k_B}{\mu m_p}T=\frac{P}{\rho }=(\gamma -1)ϵ$$

where $k_B$ is the Boltzmann constant, µ is the mean particle (molecular) weight, $m_P$ is the proton mass, $\gamma$ is the specific heat ratio, and T is the temperature.

We parameterize the cooling function as:

$$\Lambda ={\Lambda }_0{\rho }^m{T}^n,{\Lambda }_0\in {\mathbb{R}}^{+}$$

For the values of m and n:

• Depending on the cooling mechanism, we can have m = 1 or 2:
  • m = 2 for collisional cooling
  • m = 1 for thermal radiation
• Since it’s a nonlinear steep dependence on temperature, we can approximate with a power law that has n larger than 1.

Next, we define the mass per unit area as:

$$\Sigma ={\int }_0^x\rho dx$$

From Gauss’s law for gravity:

$$\int \text{g}\cdot \text{n}=2g=-4\pi G2\Sigma$$ $$\to g=-2\pi G2\Sigma =-\nabla \varphi$$

Thus, in the x dimension:

$$\frac{\partial \varphi }{\partial x}=2\pi G2\Sigma$$

Now, the one-dimensional system for planar symmetry can be written as follows:

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}=0\text{\hspace{1em}(1)}$$ $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=-\frac{1}{\rho }\frac{\partial P}{\partial x}-4\pi G\Sigma \text{\hspace{1em}(2)}$$ $$\Sigma ={\int }_0^x\rho dx\text{\hspace{1em}(3)}$$ $$\rho \left(\frac{\partial ϵ}{\partial t}+u\frac{\partial ϵ}{\partial x}\right)=-P\frac{\partial u}{\partial x}-{\Lambda }_0{\rho }^m{T}^n\text{\hspace{1em}(4)}$$

Self-similar equations

We will use the above equations and the ideal gas law, along with the following group transformation:

$$x=\lambda \hat{x},t={\lambda }^a\hat{t},u={\lambda }^b\hat{u},T={\lambda }^c\hat{T},\rho ={\lambda }^d\hat{\rho },\varphi ={\lambda }^e\hat{\varphi }$$

Here, the symbols with hats represent the scaled physical quantities in relation to the unscaled parent system by the scaling factor $\lambda$. For the transformed system to have the same form as the original one (so that self-similar solutions can exist), the constants a, b, c, d and e must satisfy:

$$1-a=b=\frac{c}{2}=1+\frac{d}{2}=\frac{e}{2}=\frac{3-2m}{5-2m-2n}\text{\hspace{1em}(5)}$$

The following similarity Ansatz (inspired by Murakami et al. 2004) would allow us to eliminate the temporal dependence from the one-dimensional system equations. Basically, all the scales of the physical quantities are uniquely determined as a function of time only:

$$x=A|t{|}^{1/a}\xi \text{\hspace{1em}(6)}$$ $$u=\frac{A}{a}|t{|}^{b/a}\nu (\xi )\text{\hspace{1em}(7)}$$ $$T=\frac{A^2}{a^2}\frac{\mu m_p}{k_B}|t{|}^{c/a}\tau (\xi )\text{\hspace{1em}(8)}$$ $$\rho =B|t{|}^{-2}g(\xi )\text{\hspace{1em}(9)}$$ $$\frac{\partial \varphi }{\partial x}=2\pi G2\Sigma =4\pi GAB|t{|}^{(e-1)/a}\sigma (\xi )\text{\hspace{1em}(10)}$$ $$\Sigma ={\Sigma }_0+{\int }_0^x\rho dx=AB|t{|}^{1/a-2}\sigma (\xi )\text{\hspace{1em}(11)}$$ $$\sigma (\xi )={\sigma }_0+{\int }_0^{\xi }g(\xi )d\xi \text{\hspace{1em}(12)}$$

Here, A and B are positive constants that establish the scales of the radius and density. Note that d/a = -2 in the density Ansatz (9), whatever the values of m and n.

With this Ansatz, the mass, momentum and energy equations reduce to ordinary differential equations:

$$(\nu \mp \xi )g^{\prime }+({\nu }^{\prime }\pm d)g=0\text{\hspace{1em}(13)}$$ $$\pm b\nu +(\nu \mp \xi ){\nu }^{\prime }+\frac{(g\tau {)}^{\prime }}{g}+K_1\sigma =0\text{\hspace{1em}(14)}$$ $$\frac{\pm c\tau +(\nu \mp \xi ){\tau }^{\prime }}{\gamma -1}+\tau {\nu }^{\prime }+K_2{g}^{m-1}{\tau }^n=0\text{\hspace{1em}(15)}$$

where the prime denotes the derivative with respect to ξ and the plus and minus signs correspond to t>0 and t<0. You may have noticed that we had four equations in the system, where the third equation is Poisson’s equation for gravity. But this equation is automatically satisfied in the above step, so we don’t need to include its reduced form here. These three ODEs are a second-order ODE system for g, τ and v.

Apparently, the system above is characterized by two dimensionless parameters $K_1$ and $K_2$:

$$K_1=4\pi a^{2}GB$$ $$K_2={\Lambda }_{0}a{B}^{m-1}{\left(\frac{A}{a}\right)}^{2(n-1)}{\left(\frac{\mu m_p}{k_B}\right)}^n$$

Additionally, from the mass conservation law:

$$\frac{\partial \Sigma }{\partial t}=-\rho u\text{\hspace{1em}(16)}$$ $$\sigma =\frac{g(\xi \mp \nu )}{1-2a}\text{\hspace{1em}(17)}$$

Asymptotic behaviors

An analytical solution is not available for this nonlinear system, except for a single trivial solution when t < 0. However, we can derive asymptotic solutions in some limits and use numerical integration to get full solutions. We use the terminology introduced by Whitworth & Summers (1985) and identify three regimes: the early interior path (t < 0 and ξ is small), the exterior path (ξ is large), and the late path (t > 0 and ξ is small).

For now, I’ve only managed to investigate the first regime - the early interior path (t<0, ξ <<1).

Early interior path

The time interval t < 0 corresponds to the period preceding the formation of a mass singularity at the core of the collapsing area. During this phase, the central region exhibits a homologous collapse with a flat density and temperature profile. We assume the following polynomial series:

$$g=1+g_1\xi +g_2{\xi }^2+\mathcal{O}({\xi }^3)\text{\hspace{1em}(18)}$$ $$\nu ={\nu }_0+{\nu }_1\xi +{\nu }_2{\xi }^2+{\nu }_3{\xi }^3+\mathcal{O}({\xi }^4)\text{\hspace{1em}(19)}$$ $$\tau =1+{\tau }_1\xi +{\tau }_2{\xi }^2+\mathcal{O}({\xi }^3)\text{\hspace{1em}(20)}$$

As A and B are free parameters, we can choose the scaling so that the leading terms of g and τ are 1. The physical boundary conditions are:

$$g^{\prime }(0)=0,\nu (0)=0,{\tau }^{\prime }(0)=0$$

These mean there is no singularity at the origin, and no density or temperature gradient there. Substituting them into the equations gives $g_1={\nu }_0={\nu }_2={\tau }_1=0$. By symmetry, all odd-order terms of g and τ, and all even-order terms of ν, vanish. This simplifies the polynomial series to:

$$g=1+g_2{\xi }^2,\nu ={\nu }_1\xi +{\nu }_3{\xi }^3,\tau =1+{\tau }_2{\xi }^2\text{\hspace{1em}(21)}$$

Inserting into the ODE system, we can get:

$${\nu }_1=d=\frac{4n-4}{5-2m-2n}\text{\hspace{1em}(22)}$$ $$2g_2({\nu }_1+1)+3{\nu }_3=0\text{\hspace{1em}(23)}$$ $$g_2+{\tau }_2=-\frac{{\nu }_1^2}{4}-\frac{K_1}{2}\text{\hspace{1em}(24)}$$ $${\tau }_2\left[\frac{2({\nu }_1+1)}{\gamma -1}+(n-1)K_2\right]=g_2\left[2({\nu }_1+1)-(m-1)K_2\right]\text{\hspace{1em}(25)}$$ $$g_2=\left(-\frac{{\nu }_1^2}{4}-\frac{K_1}{2}\right)\frac{2\frac{{\nu }_1+1}{\gamma -1}+(n-1)K_2}{2\gamma \frac{{\nu }_1+1}{\gamma -1}+(n-m)K_2}\text{\hspace{1em}(26)}$$ $$K_2=\frac{c}{\gamma -1}-d\text{\hspace{1em}(27)}$$

So $K_1$ sets both the magnitude and the sign of $g_2$. Setting ${\tau }_0=1$ in the series is equivalent to fixing $K_2$: changing ${\tau }_0$ and $K_2$ only rescales the system, and gives exactly the same solution in physical units.

From the constraints above, it is apparent that for $g_2$ and ${\tau }_2$ to be equal to 0, we must have:

$$K_1=K_0=-\frac{{\nu }_1^2}{2}$$

The effective polytropic index can be obtained by substituting the asymptotic values into the equation that represents the effective equation of state (EOS),

$$P\propto {\rho }^{\Gamma }\text{\hspace{1em}(28)}$$

leading to:

$${\tau }_2=g_2({\Gamma }_i-1)\text{\hspace{1em}(29)}$$

Using the constraints above again, we have:

$${\Gamma }_i=\frac{g_2+{\tau }_2}{g_2}=\frac{2\gamma (d+1)+(n-m)(c-(\gamma -1)d)}{2(d+1)+(n-1)(c-(\gamma -1)d)}\text{\hspace{1em}(30)}$$

Apparently, ${\Gamma }_i$ is a weighted mean between $\gamma$ and $(n-m)/(n-1)$.

Sonic point

The sonic point matters because a physical solution has to pass through it smoothly, but the energy equation makes it hard to analyse analytically.

By expressing the set of ODEs as linear combinations of their derivatives, we can derive:

$$g^{\prime }=\frac{{\Delta }_1}{\Delta },{\nu }^{\prime }=\frac{{\Delta }_2}{\Delta },{\tau }^{\prime }=\frac{{\Delta }_3}{\Delta }\text{\hspace{1em}(31)}$$

where:

$$\Delta =(\nu \mp \xi )\left[\tau -\frac{(\nu \mp \xi {)}^2-\tau }{\gamma -1}\right]$$ $${\Delta }_1=-\frac{g}{\gamma -1}\left[\mp c\tau \mp \frac{K_1}{1-2a}g(\nu \mp \xi {)}^2\pm b\nu (\nu \mp \xi )\mp d(\nu \mp \xi {)}^2\right]+K_2{g}^m{\tau }^n\mp dg\tau$$ $${\Delta }_2=\frac{\nu \mp \xi }{\gamma -1}\left[\mp c\tau \mp \frac{K_1}{1-2a}g(\nu \mp \xi {)}^2\pm b\nu (\nu \mp \xi )\mp d\tau \right]-(\nu \mp \xi )K_2{g}^{m-1}{\tau }^n$$ $${\Delta }_3=[(\nu \mp \xi {)}^2-\tau ](K_2{g}^{m-1}{\tau }^n\pm \frac{c\tau }{\gamma -1})\pm \tau \frac{K_1}{1-2a}g(\nu \mp \xi {)}^2\mp b\nu \tau (\nu \mp \xi )\pm d{\tau }^2$$

whew yeah that’s a lot…

So in order for the solutions with t<0 to transition smoothly through the sonic surface, it is necessary for both conditions,

$$\Delta ={\Delta }_2=0\text{\hspace{1em}(32)}$$

to be met simultaneously at the sonic point, which is where the solution intersects the sonic surface.

Two distinct types of solutions can smoothly pass through the sonic point.

• The first type is a trivial case with $K_1=K_0=-\frac{{\nu }_1^2}{2}$, resulting in constant density and temperature. This leads to a linear velocity profile.
• The second valid value of $K_1$ needs to be determined through numerical integration, using the shooting method towards the sonic point. Essentially, we solve the above problem by reducing it to an initial value problem. We find solutions for different initial values of $K_1$ until we find the solution that also satisfies the sonic surface conditions mentioned above.

Results

Trivial solution

When $K_1=K_0=-\frac{{\nu }_1^2}{2}$, a trivial solution is pretty easy to find:

$$g(\xi )=1,\nu (\xi )=d\xi ,\tau (\xi )=1.\text{\hspace{1em}(33)}$$

This solution characterizes the homologous collapse of a flat structure that stretches infinitely and is applicable solely for t < 0.

Physical solution

We solve the system numerically, starting the integration at small ξ and t < 0. With the shooting method, we find the value of K1 that gives a continuous solution, and then extend the t < 0 solution to very large ξ. The figures below show profiles for two pairs of m and n.

Self-similar solution for two cases: (m, n) = (1.5, 3), (2, 5) from top to bottom. The dimensionless density, velocity and temperature are shown in blue, orange and green respectively. The dashed line shows the sonic point.

From the velocity profile (orange line), we can see that the compression is faster than the sonic velocity; thus, we have a supersonic collapse.

Density-temperature diagram for two cases: (m, n) = (1.5, 3), (2, 5) from top to bottom. The dashed line shows the slope from the asymptotic solution.

We recall that ${\tau }_2=g_2({\Gamma }_i-1)$, where ${\Gamma }_i-1$ is shown on the figures above. Obviously, both the numerical and asymptotic solutions give the same slope values.