EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 508 / No 2 (December III 2009) A&A, 508 2 (2009) 725-735 Full HTML
Free Access
Issue
A&A
Volume 508, Number 2, December III 2009
Page(s) 725 - 735
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/200912806
Published online 27 October 2009

A&A 508, 725-735 (2009)

Fragmentation of a dynamically condensing radiative layer

K. Iwasaki - T. Tsuribe

Department of Earth and Space Science, Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan

Received 1 July 2009 / Accepted 14 September 2009

Abstract
In this paper, the stability of a dynamically condensing radiative gas layer is investigated by linear analysis. Our own time-dependent, self-similar solutions describing a dynamical condensing radiative gas layer are used as an unperturbed state. We consider perturbations that are both perpendicular and parallel to the direction of condensation. The transverse wave number of the perturbation is defined by k. For k=0, it is found that the condensing gas layer is unstable. However, the growth rate is too low to become nonlinear during dynamical condensation. For $k\ne 0$, in general, perturbation equations for constant wave number cannot be reduced to an eigenvalue problem due to the unsteady unperturbed state. Therefore, direct numerical integration of the perturbation equations is performed. For comparison, an eigenvalue problem neglecting the time evolution of the unperturbed state is also solved and both results agree well. The gas layer is unstable for all wave numbers, and the growth rate depends a little on wave number. The behaviour of the perturbation is specified by $kL_{\rm cool}$ at the centre, where the cooling length, $L_{\rm cool}$, represents the length that a sound wave can travel during the cooling time. For $kL_{\rm cool}\gg1$, the perturbation grows isobarically. For $kL_{\rm cool}\ll1$, the perturbation grows because each part has a different collapse time without interaction. Since the growth rate is sufficiently high, it is not long before the perturbations become nonlinear during the dynamical condensation. Therefore, according to the linear analysis, the cooling layer is expected to split into fragments with various scales.

Key words: hydrodynamics - instabilities - ISM: kinematics and dynamics - ISM: structure - ISM: clouds

1 Introduction

In the interstellar medium (ISM), it is well known that a clumpy low-temperature phase (cold neutral medium, or CNM) and a diffuse high-temperature phase (warm neutral medium, or WNM) can coexist in pressure equilibrium as a result of the balance of radiative cooling and heating due to external radiation fields and cosmic rays (Field et al. 1969; Wolfire et al. 2003,1995). These two phases are thermally stable. On the other hand, gas is thermally unstable in the temperature range between two stable phases, that is, in the range $300~{\rm K} <T< 6000$ K. The unstable gas spontaneously turns into the mixture of stable CNM and WNM by thermal instability (TI). This instability is expected to play an important role in the structure formation and the dynamics of the ISM, and especially, in the molecular cloud formation and origin of turbulence.

The basic properties of TI was investigated by Field (1965), who performed linear analysis of an uniform gas in thermally equilibrium. He derived a criterion for TI. Focusing on one fluid element, Balbus (1986) generalized the Field criterion when the cooling rate is not equal to the heating rate. The effect of magnetic field on TI has been investigated by Field (1965) and Hennebelle & Passot (2006), and other authors.

Recently, many authors have used multi-dimensional numerical simulations to study the turbulent CNM formation driven by TI. Koyama & Inutsuka (2002) suggest that the turbulent CNM formation is induced by TI in a shock-compressed region. Analogous processes have been studied by many authors for a colliding flow of the WNM (Audit & Hennebelle 2005; Vazquez-Semadeni et al. 2007; Heitsch et al. 2006; Hennebelle & Audit 2007), and using two-fluid MHD simulation (Inoue & Inutsuka 2008). The unbalance between cooling and heating rates causes the shock-compressed gas layer to cool and to condense. During the cooling, these numerical simulations shows that the runaway cooling layer quickly fragments into many CNM clumps whose velocity dispersion is equal to a fraction of the sound speed of WNM, where CNM clumps and WNM are tightly interwoven. This complex structure is regarded as produced by TI and possibly by some other hydrodynamical instabilities, such as the nonlinear thin shell instability (Vishniac 1994), the Kelven-Helmholtz instability, and by corrugation instability of the phase transition layers between CNM and WNM (Inoue et al. 2006).

A fluid element that is compressed by a shock wave tends to be a layer rather than a sphere because it is only compressed in the direction perpendicular to the shock front. Once the fluid element enters the thermally unstable regime, the layer cools in a runaway fashion. In this paper, we focus on the fragmentation of the runaway cooling layer. In previous studies, A detailed physical mechanism of the fragmentation of the runaway cooling layer remains poorly understood even in linear regime. The main reason is that it is difficult to select the unperturbed state since the cooling layer evolves temporarily and spatially. Therefore, in previous works, the unperturbed states were limited to spatially uniform gas that cools isochorically (Burkert & Lin 2000; Schwarz et al. 1972) or isobarically (Koyama & Inutsuka 2000).

Iwasaki & Tsuribe (2008) (hereafter IT08) have recently found a family of self-similar (S-S) solutions describing the dynamical condensation of a radiative gas layer where the cooling rate dominates the heating rate. This S-S solution assumed that the net cooling rate is a power-law function and that the heating rate is explicitly neglected. Although it is still too ideal, they are expected to be a good nonlinear one-dimensional model at least in the phase during the transition from WNM to CNM. In this paper, we adopt the S-S solutions as a more realistic unsteady unperturbed state than those in previous works. We perform linear analysis of the S-S solutions against fluctuations perpendicular, as well as parallel to, the direction of the condensation. By performing the linear analysis, we will have some useful insights when and how the cooling layer fragments. Since we focus on the above S-S unperturbed state, the nonlinear thin shell instability, Kelvin-Helmholtz instability, and the corrugation instability are beyond the scope of this paper.

In Sect. 2, we formulate basic equations using a zooming coordinate. Perturbation equations are derived for the linear analysis, with a brief review of the S-S solutions. In Sect. 3, we investigate the stability of the S-S solutions taking only those fluctuations into account that are parallel to the direction of the condensation. In Sect. 4, we consider the fluctuations that are both perpendicular and parallel to the direction of the condensation. In Sect. 5, we discuss the astrophysical implication of the linear analysis and effects of the thermal conduction. Our study is summarized in Sect. 6.

2 Formulation

We consider a dynamically condensing radiative gas layer where the cooling rate dominates the heating rate. The following formula is adopted as the net cooling rate per unit volume and time:

\begin{displaymath}\rho{\cal L}(\rho,P) = \Lambda_0\rho^{2-\alpha}P^{\alpha-1}\propto \rho^2 T^\alpha\;\;{\rm erg\;cm^{-3}\;s^{-1}}. \end{displaymath} (1)

In ISM, in the temperature range of $T\la 10^7$, the main coolant is Bremsstralung for the solar metallicity (Sutherland & Dopita 1993). The cooling rate of Bremsstralung is approximated well by that with $\alpha =0.5$. In the temperature range of $2\times10^5\;{\rm K}\la T \la 10^7\;{\rm K}$, the dominant coolant is metal lines for the solar metallicity (Sutherland & Dopita 1993). The cooling rate of metal lines corresponds to that with $\alpha<0$. However, the cooling rate cannot be expressed by a single $\alpha$. In the temperature range between CNM and WNM, $300\la T\la 6000$, the dominant coolant is CII (Dalgarno & McCray 1972). In this case, the cooling rate is approximated by that with $\alpha\simeq0.6$ as is shown in Sect. 5.1.

Basic equations for a radiative gas are the continuity equation,

\begin{displaymath}\frac{\partial \rho}{\partial t} + \vec{\nabla}\cdot(\rho \vec{v})=0, \end{displaymath} (2)

the equation of motion,

\begin{displaymath}\frac{{\rm D} \vec{v}}{{\rm D} t} +\frac{1}{\rho}\vec{\nabla} P=0, \end{displaymath} (3)

and the entropy equation,

\begin{displaymath}\frac{1}{\gamma-1}\frac{{\rm D} }{{\rm D} t}\left( \ln P \rho... ...ght) = - \gamma^{\alpha}\Lambda_0\rho^{2-\alpha} P^{\alpha-1}, \end{displaymath} (4)

where ${\rm D}/{\rm D} t=\partial/\partial t+\vec{v}\cdot\vec{\nabla}$ indicates the Lagrangian time derivative.

We take the x-axis as the direction of the condensation driven by the cooling and y-axis as the transverse direction. Since the S-S solutions are time-dependent, it is difficult to perform linear analysis in the ordinary Cartesian coordinate, (t,x,y). Bouquet et al. (1985) introduced a zooming coordinate where S-S solutions appear to be stationary (also see Hanawa & Matsumoto 1999). We introduce the similar zooming coordinate since this transformation makes stability analysis easier as follows:

\begin{displaymath}\left( \begin{array}{c} t \\ x \\ y \\ \end{array}\ri... ... y \\ \end{array}\right), \;\; x_0(t)=at_*^{1/(1-\omega)}, \end{displaymath} (5)

where $\omega$ is a free parameter, a corresponds to the cooling length, and $t_* = 1-t/t_{\rm c}$, $t_{\rm c}$ is an epoch when the central density becomes infinity. In the zooming coordinate, density $\Omega$, velocity $\vec{V}$, pressure $\Pi$, and sound speed X are given by
$\displaystyle \Omega(\tau,\xi,y) = \rho(t,x,y)/\rho_0(t),$   $\displaystyle \vec{V}(\tau,\xi,y)=\vec{v}(t,x,y)/v_0(t),$  
$\displaystyle \Pi(\tau,\xi,y) = P(t,x,y)/P_0(t),$   $\displaystyle X(\tau,\xi,y) = c_{\rm s}(t,x,y)/v_0(t),$ (6)

respectively, where

\begin{displaymath}v_0(t)= -\dot{x}_0(t)=\frac{a}{(1-\omega)t_{\rm c}}t_*^{\omega/(1-\omega)}, \end{displaymath} (7)

\begin{displaymath}\rho_0(t) = \left\{ \frac{1}{(1-\omega)t_{\rm c}} \right\}^{-... ...ha+3} \frac{a^{2(1-\alpha)}}{\Lambda_0}t_*^{\beta/(1-\omega)}, \end{displaymath} (8)

and

\begin{displaymath}P_0(t) = \frac{\rho_0v_0^2}{\gamma}\propto t_*^{(2\omega+\beta)/(1-\omega)}, \end{displaymath} (9)

with $\beta=\omega(3-2\alpha)-1$.

In the zooming coordinate, the basic Eqs. (2)-(4) are rewritten as

\begin{displaymath}\frac{{\rm D} \ln \Omega}{{\rm D} \tau} + \vec{\nabla}\cdot\vec{V} = \beta, \end{displaymath} (10)

\begin{displaymath}\frac{{\rm D} \vec{V}}{{\rm D} \tau} +\frac{1}{\Omega}\vec{\nabla} \Pi = \omega \vec{V}, \end{displaymath} (11)

and

\begin{displaymath}\frac{1}{\gamma-1}\frac{{\rm D} }{{\rm D} \tau} \left( \ln \... ...amma-1} -\beta - \gamma^\alpha\Omega^{2-\alpha}\Pi^{\alpha-1}, \end{displaymath} (12)

respectively, where the operators of time and spatial derivative are defined by

\begin{displaymath}\frac{{\rm D}}{{\rm D}\tau} = \frac{\partial}{\partial \tau... ...tial \xi}, x_0(\tau)\frac{\partial}{\partial y} \right),\;\; \end{displaymath} (13)

respectively, where $\vec{e}_\xi$ indicates the unit vector parallel to $\xi$-direction. Supplements for the derivation of Eqs. (10)-(12) is presented in Appendix A.

We apply the zooming transformation only in the x-direction but not in the y-direction. This is because the gas contracts along x-axis but not along y-axis in the unperturbed state. In the ordinary coordinate, the transverse scale of the perturbation is expected to be constant with time. However, if the zooming transformation is also applied in the y-direction, the transverse scale of the perturbation decreases with time in the ordinary coordinate, although the unperturbed gas does not contract along the y-axis. Therefore, we apply the zooming transformation only in the x-direction.

2.1 Review of self-similar solutions

In the zooming coordinate, steady state solutions correspond to S-S solutions that were derived in IT08. In this section, physical properties of the S-S solutions are reviewed briefly.

The S-S solutions are specified by two parameters, $\alpha$ and $\omega$. For convenience, instead of $\omega$, we can use a parameter $\eta $, which is given by

\begin{displaymath}\eta = \frac{(2-\alpha)\{1-(3-2\alpha)\omega\}}{1-\omega}{\cdot} \end{displaymath} (14)

Using these parameters ( $\alpha,\eta$), the time dependences of the central density, $\rho_{00}$, and pressure, P00, are given by

\begin{displaymath}\rho_{00}(t) \propto t_*^{-\alpha_{\rho}(\eta)}\;\;\;{\rm and}\;\;\; P_{00}(t)\propto t_*^{\alpha_P(\eta)}, \end{displaymath} (15)

respectively, where $\alpha_\rho(\eta) = \eta/(2-\alpha)$ and $\alpha_P(\eta) =(1-\eta)/(1-\alpha)$.

The S-S solutions include two asymptotic solutions. For $\eta\sim0$, the time dependences of the central density and pressure are given by

\begin{displaymath}\rho_{00}(t) \sim {\rm const.}\;\;\;{\rm and}\;\;\; P_{00}(t)\propto t_*^{1/(1-\alpha)}, \end{displaymath} (16)

respectively. This time evolution indicates the isochoric mode. For $\eta\sim1$, the time dependences of the central density and pressure are given by

\begin{displaymath}\rho_{00}(t) \propto t_*^{-1/(2-\alpha)}\;\;\;{\rm and}\;\;\; P_{00}(t)\sim{\rm const.}, \end{displaymath} (17)

respectively. This time evolution corresponds to the isobaric mode. Our S-S solutions exist between the two limits, $0<\eta<1$. With the condition that the increasing rate of $\rho_{00}(t)$ is equal to the decreasing rate of P00(t), or $\alpha_\rho(\eta) = \alpha_P(\eta)$, the critical value of $\eta $ is given by

\begin{displaymath}\eta_{\rm eq} = (2-\alpha)/(3-2\alpha). \end{displaymath} (18)

Therefore, the S-S solutions for $0<\eta<\eta_{\rm eq}$ and $\eta_{\rm eq}<\eta<1$ are close to the isochoric and the isobaric modes, respectively. During S-S condensation, the central pressure should not increase with time. From Eq. (15), the range of the value of $\alpha$ is found to be $\alpha<1$.

2.2 Perturbation equations

Perturbation on the S-S solutions is considered. Perturbed variables are defined by

                          $\displaystyle \Omega = \Omega_0(\xi)\{1+\delta \Omega(\tau,\xi,y)\},$  
    $\displaystyle V_x = V_0(\xi) + \delta V_x(\tau,\xi,y),$  
    $\displaystyle V_y = \delta V_y(\tau,\xi,y),$  
    $\displaystyle {\rm and}$  
    $\displaystyle \Pi = \Pi_0(\xi)\{1 + \delta \Pi(\tau,\xi,y)\},$ (19)

where subscript ``0'' indicates the unperturbed state. As the perturbation, we consider the following Fourier mode with respect to y,

\begin{displaymath}\vec{\delta Q}(\tau,\xi,y) = \vec{\delta Q}(\tau,\xi) \exp(ik... ...\vec{\delta Q} = (\delta \Omega, \delta \vec{V},\delta \Pi), \end{displaymath} (20)

and k indicates the wave number of the plane wave that propagates along y-direction. Substituting Eqs. (19) and (20) into Eqs. (10)-(12) and linearizing, we get the following perturbation equations:

\begin{displaymath}\frac{{\rm D} \delta \Omega}{{\rm D} \tau} + \frac{\partial ... ... - (\ln \Omega_0)' \delta V_x - \kappa(\tau){\rm i}\delta V_y, \end{displaymath} (21)

\begin{displaymath}\frac{{\rm D} \delta V_x}{{\rm D} \tau} + \frac{X_0^2}{\gamm... ...rac{X_0^2}{\gamma}(\ln \Pi_0)'( \delta \Pi - \delta \Omega), \end{displaymath} (22)

\begin{displaymath}\frac{{\rm D} {\rm i}\delta V_y}{{\rm D} \tau} = \omega {\rm i}\delta V_y+ \kappa(\tau)\frac{X_0^2}{\gamma}\delta \Pi, \end{displaymath} (23)

and

$\displaystyle \frac{{\rm D} \delta \Pi }{{\rm D} \tau} -\gamma \frac{{\rm D} \d... ...gamma\epsilon_0\delta \Pi -\left(\ln \Pi_0\Omega_0^{-\gamma}\right)'\delta V_x,$     (24)

where ${\rm D}/{\rm D}\tau=\partial/\partial \tau+V_0 \partial /\partial \xi$,

\begin{displaymath}\epsilon_0 = \gamma^{\alpha-1}(\gamma-1)\Omega_0^{2-\alpha}\Pi_0^{\alpha-1}, \;\;{\rm and}\;\;\kappa(\tau) \equiv kx_0(\tau). \end{displaymath} (25)

The time-dependent factors remain in the form of $\kappa(\tau)=kx_0(\tau)$ in (21)-(24), because the transverse scale is not zoomed (see Eq. (5)) as mentioned above. The leads to a problem that the perturbed variables cannot be expanded in the Fourier mode with respect to $\tau$ in general.

\begin{figure} \par\includegraphics[width=12cm,clip]{12806fg1.eps} \end{figure} Figure 1:

Growth rate, $\Sigma $, as a function of $\eta $ for a)  $\alpha =-1.0$ and b) 0.5 for the case with the perturbation parallel to the condensation. For comparison, the increasing rate of the unperturbed central density, $\alpha _\rho (\eta )$, and the decreasing rate of the unperturbed central pressure, $\alpha _P(\eta )$, are shown by the dashed and dotted lines, respectively.

Open with DEXTER

3 Perturbation for k = 0

In this section, we consider the perturbation parallel to the condensation, or for the case with k=0. In this case, since the time-dependent factor, $\kappa(\tau)=kx_0(\tau)$, vanishes, the perturbed variables can be expanded in the Fourier mode with respect to $\tau$ as

\begin{displaymath}\delta Q(\tau,\xi) = \delta Q(\xi) {\rm e}^{\sigma \tau}. \end{displaymath} (26)

By Eq. (26), the time evolution of the perturbations is given by

\begin{displaymath}\delta Q\propto t_*^{-{\Sigma}},\;\;\;{\rm where}\;\; {\Sigma}=\frac{\sigma}{1-\omega}{\cdot} \end{displaymath} (27)

Substituting Eq. (26) into the perturbation Eqs. (21)-(24), one obtains the following ordinary differential equations:

\begin{displaymath}\frac{{\rm d}\delta Q_i}{{\rm d}\xi} = \frac{1}{V_0^2-X_0^2... ...\; \vec{\delta Q} = (\delta \Omega, \delta V_x, \delta \Pi), \end{displaymath} (28)

where the detailed expression of Aij is shown in Appendix B. Equations (28) are solved as a boundary- and eigenvalue problem.

3.1 Boundary conditions

We impose the boundary conditions at $\xi=0$ and at the critical point, $\xi=\xi_{\rm s}$, where V0=X0. The boundary conditions at $\xi=0$ are obtained by the asymptotic limit of the perturbed variables. From the regularity of the perturbed variables at $\xi=0$, we find that perturbed variables should have the following asymptotic forms:

\begin{displaymath}\lim_{\xi\rightarrow0}\delta \Omega(\xi) \simeq \delta \Omega... ...rightarrow0}\delta V_x(\xi) \simeq -\sigma\delta \Omega_0 \xi, \end{displaymath} (29)

and

\begin{displaymath}\lim_{\xi\rightarrow0} \delta \Pi(\xi) \simeq \gamma\frac{\s... ...00}} {\sigma + (\alpha-1)\gamma\epsilon_{00}}\delta \Omega_0, \end{displaymath} (30)

where $\epsilon_{00}=\{2\omega - \beta(\gamma -1)\}/\gamma$.

The boundary conditions at the critical point, $\xi=\xi_{\rm s}$, are obtained from Eqs. (28). At $\xi=\xi_{\rm s}$, the denominator of the righthand side becomes zero. To obtain a regular solution from $\xi=0$ to $\xi=\infty$, the numerator of the righthand side should vanish. Therefore, the boundary conditions are given by the following three equations,

\begin{displaymath}\sum_{j=1}^3 A_{ij}\delta Q_j\Biggr\vert _{\xi=\xi_{\rm s}}=0,\;\; i=1,2,3. \end{displaymath} (31)

The above three equations give only one independent condition.

3.2 Numerical method

Solutions of Eqs. (28) have three integration constants. Therefore, if we set two constants $(\delta \Omega_0,\sigma)$ and impose the boundary condition at $\xi=\xi_{\rm s}$, the solution is completely fixed. In this section, our numerical method for solving Eqs. (28) is described.

We can set $\delta \Omega_0=1$ without loss of generality. For a given $\sigma$, we integrate Eqs. (28) from $\xi=0$ to the critical point, $\xi=\xi_{\rm s}$, using a fourth-order Runge-Kutta method. Eigenvalue, $\sigma$, is modified until the perturbed variables satisfy the boundary condition at $\xi=\xi_{\rm s}$using the Newton-Raphson method. After that, we integrate Eqs. (28) up to $\xi=10^4$.

3.3 Results

Figure 1 shows the dependence of growth rate, $\Sigma $, on $\eta $ for (a)  $\alpha =-1.0$ and (b) 0.5. From Fig. 1, it is seen that $\sigma>0$ for a wide range of $\alpha$ and $\eta $. Therefore, the perturbation is unstable. However, $\sigma$ is smaller than $\alpha_P$ and $\alpha_\rho$. Therefore, the growth rate is too low to grow sufficiently during runaway condensation. This is consistent with the results of the one-dimensional numerical simulation shown in Sect. 5.1.

4 Perturbation with k $\ne $ 0

For $k\ne 0$, the perturbation equations contain the time-dependent factor, $\kappa(t)$. Here we show that this factor, $\kappa=kx_0(t)$, is related to the ratio of the local cooling length at the centre to the wave length of the perturbation. The cooling timescale at t=0 is given by

\begin{displaymath}t_{\rm cool}^{0} = \frac{P_{00}}{(\gamma-1)\gamma^\alpha\La... ...gamma-1)} \Omega_{00}^{\alpha-2}\Pi_{00}^{1-\alpha}t_{\rm c}, \end{displaymath} (32)

where superscript ``0'' indicates the value at t=0. Therefore, using $t_{\rm cool}^{0}$, the collapse time can be expressed as

\begin{displaymath}t_{\rm c} = \frac{\gamma^\alpha(\gamma-1)}{1-\omega}\Omega_{00}^{2-\alpha}\Pi_{00}^{\alpha-1} t_{\rm cool}^{0}. \end{displaymath} (33)

Using Eqs. (6) and (7), at t=0, the sound speed at the centre is given by

\begin{displaymath}c_{00}^0 = X_{00}v_0^0,\;\;{\rm where}\;\;v_0^0t_{\rm c} = a/(1-\omega) = x_0(0)/(1-\omega). \end{displaymath} (34)

Therefore, x0(0) is given by

\begin{displaymath}x_0(0) = \frac{1-\omega}{X_{00}}c_{00}^0 t_{\rm c}. \end{displaymath} (35)

Using Eq. (33), Eq. (35) can be written as

\begin{displaymath}x_0(0) = M_{00} L_{\rm cool}^0=a, \end{displaymath} (36)

where $M_{00}= \gamma^{\alpha-1/2}(\gamma-1) \Omega_{00}^{5/2-\alpha}\Pi_{00}^{\alpha-3/2}$, and $L_{\rm cool}^{0} = c_{00}^{0} t_{\rm cool}^{0}$is the cooling length at the centre at t=0. Since the S-S solutions are invariant under the transformation $\xi\rightarrow m\xi$, $V_0\rightarrow mV_0$, $\Omega_0\rightarrow m^{-2(\alpha-1)} \Omega_0$, and $\Pi_0\rightarrow m^{-2(\alpha-2)}\Pi_0$, with a parameter m, one can take m in such a way that M00 is equal to 1. Hereafter, M00 is set to unity. Therefore, the time-dependent factor, $\kappa(\tau)$, is expressed as

\begin{displaymath}\kappa(\tau) = k x_0(\tau) = k L_{\rm cool}^{0} t_*^{1/(1-\om... ...= kL_{\rm cool}^0{\rm e}^{\rm -\tau} = k L_{\rm cool}(\tau), \end{displaymath} (37)

where $L_{\rm cool}(\tau)=L_{\rm cool}^{0}{\rm e}^{-\tau}$.

4.1 Static approximation

Before we present the fully time-dependent numerical calculation in Sect. 4.2, we at first consider a special case in which the time evolutions of the unperturbed S-S solutions are slower than the growth of perturbation. In this situation, the time evolution of the unperturbed state is negligible during the growth of the perturbations. Therefore, we set x0 to be a constant in Eqs. (21)-(24). This approximation is also valid for $k\ll1/x_0$ where the term, $\kappa $, is negligibly small in Eqs. (21)-(24). We use the Fourier mode as

\begin{displaymath}\vec{\delta Q}(\xi,y,\tau) = \vec{\delta Q}(\xi) {\rm e}^{{\r... ... Q} = (\delta \Omega, \delta V_x,\delta\Pi,{\rm i}\delta V_y). \end{displaymath} (38)

The condition under which the static approximation is valid is given by
          $\displaystyle \left\vert \frac{{\rm d}\ln \delta Q}{{\rm d}\ln t_*} \right\vert = \Sigma \gg \left\vert \frac{{\rm d}\ln \kappa}{{\rm d}\ln t_*} \right\vert$ = $\displaystyle \frac{(2-\alpha)(3-2\alpha) - \eta}{2(2-\alpha)(1-\alpha)},$  
  = $\displaystyle \left\{\begin{array}{cc} (5-2\alpha)/\left\{ 2(2-\alpha) \right\}... ...)/\left\{ 2(1-\alpha) \right\}& \;\;{\rm for}\;\;\eta=0, \\ \end{array}\right.$ (39)

where the definition of $\Sigma $ is the same as that in Eq. (27). Substituting Eq. (38) into the perturbation Eqs. (21)-(24), we get

\begin{displaymath}\frac{{\rm d}\delta Q_i}{{\rm d}\xi} = \frac{1}{V_0^2 - X_0... ... Q} = (\delta \Omega, \delta V_x,\delta\Pi,{\rm i}\delta V_y), \end{displaymath} (40)

where the detailed expression of Aij is shown in Appendix B.

We impose the boundary conditions at $\xi=0$ and $\xi=\xi_{\rm s}$. Since we are interested in the fragmentation of the cooling layer, only the even mode is investigated. For the even mode, the perturbed variable should have the following asymptotic forms in $\xi\ll1$:

       $\displaystyle \lim_{\xi\rightarrow0}\delta \Omega(\xi)$ $\textstyle \simeq$ $\displaystyle \delta \Omega_0,$  
$\displaystyle \lim_{\xi\rightarrow0}\delta V_x(\xi)$ $\textstyle \simeq$ $\displaystyle \delta V_{x0}\xi,$  
$\displaystyle \lim_{\xi\rightarrow0}\delta V_y(\xi)$ $\textstyle \simeq$ $\displaystyle \delta V_{y0},$  
    $\displaystyle \hspace{-2cm}{\rm and}$  
$\displaystyle \lim_{\xi\rightarrow0}\delta \Pi(\xi)$ $\textstyle \simeq$ $\displaystyle \delta \Pi_0 + \delta \Pi_{01}\xi^2.$ (41)

Substituting Eqs. (41) into Eqs. (40), we obtain the following relations:

\begin{displaymath}\sigma\Omega_0 + \delta V_{x0} + \kappa i\delta V_{y0}=0, \end{displaymath} (42)

\begin{displaymath}(\sigma-\omega)i\delta V_{y0} = \frac{\kappa X_{00}^2}{\gamma}\delta \Pi_0, \end{displaymath} (43)

\begin{displaymath}\left\{\sigma + (\alpha-2)\epsilon_{00}\right\}\gamma\delta \... ...\sigma + (\alpha-1)\gamma\epsilon_{00} \right\}\delta \Pi_0=0, \end{displaymath} (44)

and

\begin{displaymath}(2\beta+1 - \omega + \sigma )\delta V_{x0} + 2\frac{X_{00}^2}... ...01} = 2\beta\omega(1-\alpha)(\delta \Pi_0 - \delta \Omega_0).~ \end{displaymath} (45)

The boundary condition at $\xi=\xi_{\rm s}$ is derived in the same way in Sect. 3.1, and we obtain two independent conditions. Numerical method for solving Eqs. (40) is the same as that in Sect. 3.2.

Figure 2 shows an approximate dispersion relation for $\alpha =0.5$ with (a) $\eta=0.95$, (b) 0.75, and (c) 0.10. In Fig. 2, one can see several branches labelled by numbers. The most unstable branch is labelled by (1) for large $\kappa $ limit, and (2) for small $\kappa $limit. Each branch is explained below.

\begin{figure} \par\includegraphics[width=6cm,clip]{12806fg2a.eps} \vspace*{2mm}... ...ps} \vspace*{2mm} \includegraphics[width=6cm,clip]{12806fg2c.eps} \end{figure} Figure 2:

Approximate dispersion relations for $\alpha =0.5$ with a) $\eta =0.1$, b) 0.75, and c) 0.95, where $\kappa $ is the nondimensional wave number, and $\Sigma $ the growth rate of the perturbation. The labels of number represent the branches of modes.

Open with DEXTER

\begin{figure} \par\includegraphics[width=6cm,height=4cm,clip]{12806fg3.eps} \end{figure} Figure 3:

Eigenfunctions for $\eta =0.75$ and $\kappa =8$. The filled circles indicate the values at the critical point. The eigenvalue, $\Sigma $, is equal to 1.196.

Open with DEXTER


(1) The isobaric mode

\begin{figure} \par\includegraphics[width=12cm,clip]{12806fg4.eps} \end{figure} Figure 4:

Growth rate as a function of $\eta $ for a) $\alpha =0.5$ and b) 0.8 for the case with $k\ne 0$. The thick solid lines correspond to the growth rate of the isobaric mode, $\Sigma _{\rm isobaric}$. For comparison, the evolutionary rate of the unperturbed central density, $\alpha _\rho (\eta )$, and pressure, $\alpha _P(\eta )$ are shown by the dashed and dotted lines, respectively. The thin solid lines pointed by arrows correspond to the growth rate in the noninteractive mode, $\Sigma _{\rm non-int}$.

Open with DEXTER

The branch (1) corresponds to the most unstable mode for $\kappa\gg1$, or $1/k\ll L_{\rm cool}$. Since $1/k\ll L_{\rm cool}$, the sound wave can travel the wave length of the perturbation many times during the runaway cooling of the unperturbed state. Therefore, the perturbation is expected to grow in pressure equilibrium with its surroundings, and the mode corresponds to the isobaric mode. An example of eigenfunctions $(\delta \Omega, \delta V_x, \delta V_y, \delta \Pi)$ for the branch (1) is shown in Fig. 3 for $\eta =0.75$ and $\kappa=8.0$. In Fig. 3, it is clearly seen that $\vert\delta \Pi\vert\ll\vert\delta \Omega\vert$. This behaviour also implies the isobaric mode.

The growth rate in the isobaric mode can also be derived analytically by considering the evolution of a fluid element at the centre. The fluid element is assumed to have an isobaric fluctuation, $\rho=\rho_{00}(t) + \delta \rho_{00}$ and P=P00(t). From Eq. (4), the following perturbation equation is obtained:

\begin{displaymath}\frac{\partial}{\partial t}\left( \frac{\delta \rho_{00}}{\rh... ...a}P_{00}^{\alpha-1} \frac{\delta\rho_{00}}{\rho_{00}}{\cdot} \end{displaymath} (46)

From Eqs. (8) and (9), we have

\begin{displaymath}\gamma^{\alpha-1}(\gamma-1)\Lambda_0 \rho_{00}^{2-\alpha}P_{0... ...a-1} = \frac{\epsilon_{00}}{(1-\omega)(t_{\rm c}-t)}{\cdot} \end{displaymath} (47)

Using Eq. (47), Eq. (46) is rewritten as

\begin{displaymath}\frac{\partial}{\partial t}\left( \frac{\delta \rho_{00}}{\rh... ...mega)(t_{\rm c}-t)} \frac{\delta\rho_{00}}{\rho_{00}}{\cdot} \end{displaymath} (48)

Equation (48) can be integrated to give

\begin{displaymath}\frac{\delta \rho_{00}}{\rho_{00}}\propto t_*^ {-(2-\alpha)\epsilon_{00}/(1-\omega)}. \end{displaymath} (49)

Therefore, the growth rate in isobaric mode is given by
                             $\displaystyle \Sigma_{\rm isobaric}$ = $\displaystyle \frac{(2-\alpha)\epsilon_{00}}{1-\omega}$  
  = $\displaystyle \left[ 1-\frac{2-\alpha}{\gamma(1-\alpha)} \right]\eta + \frac{2-... ...for}\;\;\eta=0.75 \\ 1.04 & \;\;{\rm for}\;\;\eta=0.95. \\ \end{array}\right.$ (50)

In Fig. 2b, the eigen-value, $\Sigma $, is found to be 1.196 for the case with $\eta =0.75$ and $\kappa=8.0$. This is very close to $\Sigma_{\rm isobaric}=1.2$ for $\eta =0.75$. From Fig. 2, it is clearly seen that the growth rate approaches the corresponding $\Sigma _{\rm isobaric}$ in the large $\kappa $ limit. The analytic growth rate is derived by the local argument. Therefore, the growth rate is expected to be independent of a global structure of the system. Burkert & Lin (2000) performed a linear analysis on a spatially uniform and isochorically cooling gas. Their growth rate in the isobaric mode and our $\Sigma _{\rm isobaric}$ for $\eta=0$ are identical, although the spatial structure of the unperturbed state is quite different.

For the density perturbation to grow sufficiently during the runaway cooling, it must grow faster than the condensation of the cooling layer. This condition can be expressed by $\Sigma_{\rm isobaric}>\alpha_\rho$. Figure 4 shows $\Sigma _{\rm isobaric}$, $\alpha_\rho$ and $\alpha_P$ as a function of $\eta $ for (a)  $\alpha =0.5$ and (b) 0.8. From Fig. 4, it is found that $\Sigma _{\rm isobaric}$ is greater than $\alpha_\rho$ for all $\eta $. For other $\alpha$, we can investigate analytically as follows: subtracting $\alpha_\rho$ from $\Sigma _{\rm isobaric}$, one obtains

                    $\displaystyle \Sigma_{\rm isobaric} - \alpha_\rho$ = $\displaystyle \left[ 1-\frac{2-\alpha}{\gamma(1-\alpha)} - \frac{1}{2-\alpha}\right]\eta + \frac{2-\alpha}{\gamma(1-\alpha)}$  
  $\textstyle \ge$ $\displaystyle 1-\frac{1}{2-\alpha}>0\;\;{\rm for\;\;\alpha<1}.$ (51)

Therefore, for $\alpha<1$, $\Sigma _{\rm isobaric}$ is more than $\alpha_\rho$, and the isobaric mode can grow in the runaway cooling layer.

\begin{figure} \par\includegraphics[width=7.4cm,clip]{12806fg5.eps} \end{figure} Figure 5:

Schematic picture of the noninteractive mode.

Open with DEXTER


(2) The noninteractive mode

Branch (2) corresponds to the most unstable mode for $\kappa\ll1$. Because the wave length is larger than the cooling length, each part evolves independently according to the S-S solutions. We call this branch the noninteractive mode. Figure 5 shows the schematic picture of the noninteractive mode. In Fig. 5, similarity variables have an initial fluctuation. For example, we consider two different regions, ``A'' and ``B'', where $\rho_{\rm A}=\rho +\Delta \rho$ and $\rho_{\rm B}=\rho$. Due to the difference of $\Delta \rho$, the regions ``A'' and ``B'' have different collapse time, $t_{\rm c}-\Delta t$ and $t_{\rm c}$, respectively. Omitting any terms that do not grow, we find the time evolution of difference, $\Delta \rho$, to be

\begin{displaymath}\frac{\Delta \rho}{\rho_0(t)} = - \Omega(\xi)\Delta t\left[ \... ...\ln\Omega}{{\rm d}\ln\xi} \right] \frac{1}{t_{\rm c}-t}{\cdot} \end{displaymath} (52)

Therefore, in the zooming coordinate, the density perturbation grows as $\delta \Omega(\xi)\propto (t_{\rm c}-t)^{-1}$. Other perturbed variables also grow in the same power law. Therefore, comparing with Eq. (27), the growth rate is given by

\begin{displaymath}\Sigma_{\rm non-int} =1, \end{displaymath} (53)

which is independent of $\alpha$ and $\eta $. The noninteractive mode arises from the fluctuation of the collapse time, $t_{\rm c}$, due to density and pressure perturbations. From the physical mechanism, the perturbation grows in the same way as the unperturbed state. In other words, the perturbation of the isobaric (isochoric) cooling layer grows isobarically (isochorically).

We investigate whether the noninteractive mode grows sufficiently during the runaway cooling or not. First, we consider the case with $\eta> \eta_{\rm eq}$. Since the perturbation grows isobarically, $\Sigma _{\rm non-int}$ is compared to $\alpha_\rho=\eta/(2-\alpha)$. From Fig. 4, it is seen that the growth rate, $\Sigma _{\rm non-int}$, is higher than $\alpha_\rho$ for all $\eta>\eta_{\rm eq}= 0.75$ and 0.86 for $\alpha =0.5$ and 0.8, respectively. Therefore, the noninteractive mode can grow sufficiently. Next, we consider the cooling layer with $\eta<\eta_{\rm eq}$. In this layer, since the perturbation grows isochorically, $\Sigma _{\rm non-int}$ is compared to $\alpha_P=(1-\eta)/(1-\alpha)$. Analytically, it is found that the pressure perturbation can only grow for $\eta>\alpha$.

Koyama & Inutsuka (2000) performed a linear analysis of a spatially uniform and isobarically cooling gas in their appendix. However, they did not find the noninteractive mode. This is because they fixed the collapse time to be spatially constant, and it was assumed not to be influenced by perturbation. Burkert & Lin (2000) performed linear analysis of a spatially uniform and isochorically cooling gas by taking account of the time evolution of the unperturbed state. They showed that a perturbation cannot grow in the condition for long wave length limit. Our result is consistent with theirs.

\begin{figure} \mbox{\includegraphics[width=5.8cm,clip]{12806fg6a.eps}\hspace*{4mm} \includegraphics[width=5.8cm,clip]{12806fg6b.eps} }\end{figure} Figure 6:

Growth rate obtained by results of numerical linear analysis for $\eta =0.1$ a) and 0.75  b). The thick and dashed lines indicate the results of numerical linear analysis and those of approximate linear analysis.

Open with DEXTER


(3) The shear mode

For $\kappa\ll0$, there is a solution where physical variables are very small except for $\delta V_y$ which is spatially constant, and the eigenvalue is $\sigma=\omega$. The similar mode was found by McNamara (1993), who investigated TI of a uniform granular medium. McNamara (1993) called this mode the shear mode. The growth rate is given by

\begin{displaymath}\Sigma_{\rm shear}=\frac{\omega}{1-\omega} = \left\{ \begin... ... 0.37 & \;\;{\rm for}\;\;\eta=0.95. \\ \end{array} \right. \end{displaymath} (54)

In Fig. 2, it is seen that each growth rate in branch (3) has the corresponding value of $\Sigma_{\rm shear}$ for $\kappa\ll1$. The physical meaning of this mode can be understood as follows: for $\kappa\ll1$, since the effect of the pressure gradient with respect to y is very weak, the gas can freely stream with almost constant velocity, vy, in the y direction. On the other hand, the central sound speed, c00(t), decreases as $\propto (t_{\rm c} - t)^{\omega/(1-\omega)}$. Therefore, the ratio of the dynamical velocity to the thermal velocity, vy/c00(t), grows with time as $\propto (t_{\rm c} - t)^{-\omega/(1-\omega)}$, indicating that the growth rate is given by Eq. (54). For the case with larger wave number, the effect of the pressure gradient becomes important. Therefore, the fluid element cannot stream freely, and the growth rate is lower as shown in Fig. 2.


(4) The free-streaming mode

For large $\kappa $, there is another mode in which the velocity perturbation in the x-direction is much greater than that in the y-direction, $\vert\delta V_{x0}\vert\gg \vert\delta V_{y0}\vert$. We call this mode the free-streaming mode. From Eq. (45) with $\delta \Omega_0=\delta \Pi_0=\delta \Pi_{01}=0$, we obtain

\begin{displaymath}\Sigma_{\rm free} = -1 - \frac{2\beta}{1-\omega} = \left\{ ... ... 0.27 & \;\;{\rm for}\;\;\eta=0.95. \\ \end{array} \right. \end{displaymath} (55)

In the free-streaming mode, the growth of the velocity perturbation in the x-direction is hampered by the pressure gradient of the unperturbed state. Therefore, the growth rate is less than the shear mode.


(5) $\vec{k=0}$ mode

The growth rate in this branch for $\kappa\ll1$ coincides with the case with k=0, which is obtained in Sect. 3.

4.2 Linear analysis considering the time evolution of $\kappa $

The static approximation is valid only if $\Sigma $ is much larger than $\vert{\rm d} \ln\kappa/{\rm d}\ln t_*\vert$. However, from Fig. 2, it is found that the growth rate of the most unstable mode for each $\kappa $ is smaller than $\vert{\rm d} \ln\kappa/{\rm d}\ln t_*\vert$ whose values are 1.93, 1.5, and 1.37 for $\eta =0.1$, 0.75, and 0.95 with $\alpha =0.5$, respectively (see Eq. (39)). Therefore, in this section, we perform a linear analysis considering the time evolution of $\kappa(\tau)$ using direct numerical integration. The upwind difference method is used as the numerical method to solve the perturbation equation. The region of calculation is from $\xi=0$ to $\xi=100$ in the zooming coordinate.

We impose the boundary conditions at $\xi=0$ and $\xi=100$. The even mode is set as the boundary condition at $\xi=0$. The free boundary condition is set at $\xi=100$, but it does not influence the inner region since the gas flows out supersonically from the outer boundary of the zooming coordinate. As an initial state, we adopt the eigenfunction of the isobaric mode for $\kappa=30$ which is obtained in Sect. 4.1.

By solving Eqs. (21)-(24), a time evolution of $\vec{\delta Q}$ is obtained. During the calculation, the growth rate at $\tau$ is evaluated by

\begin{displaymath}\Sigma_{\rm num} = \frac{1}{1-\omega}\frac{{\rm d} }{{\rm d}\tau} \left\{\ln \delta \Omega(\xi=0,\tau)\right\} \end{displaymath} (56)

at each instant of time. For comparison with the result of the static approximation, we focus on a relation between $\kappa(t)$ and the growth rate of the density perturbation at the centre, $\Sigma_{\rm num}$. Figure 6 shows the growth rate $\Sigma_{\rm num}$as a function of $\kappa(\tau)$ at each instant of time for $\alpha =0.5$ with (a) $\eta =0.1$ and (b) 0.75. The nondimensional wave number, $\kappa(\tau)$, decreases with time as shown in Eq. (37). Therefore, in Fig. 6, the direction of time is from the right to the left. For comparison, the approximate dispersion relations of branches (1)-(2) in Fig. 2 are superimposed by the dashed lines. In both of Figs. 6a and 6b, the behaviour of the growth rate, $\Sigma_{\rm num}$, moderately agrees with that of the approximated dispersion relations. For $\kappa\gg1$, or initial phase, the growth rate agrees with $\Sigma _{\rm isobaric}$. This is because the growth rate does not depend on $\kappa $ in the short wave length limit. As $\kappa $ decreases and reaches about 1, $\Sigma_{\rm num}$ begins to decrease. For $\kappa\ll1$, $\Sigma_{\rm num}$ approaches asymptotically $\Sigma _{\rm non-int}$, where the effect of $\kappa $ is negligible. The effect of time-depending $\kappa $ is notable only for $0.1<\kappa<10$. Smoother dependence of the growth rate on $\kappa $ is obtained than the approximate dispersion relation.

5 Discussion

5.1 Astrophysical implication

In this section, we estimate the fragmentation time of the cooling layer formed by a collision of WNM. We consider a head-on collision between two thermally equilibrium gases in which density and pressure are $\rho_{\rm WNM}=0.57m_{\rm H}\;{\rm cm}^{-3}$ and $P_{\rm WNM}=3.5\times10^3k_{\rm B}\;{\rm dyne\;cm^{-2}}$. Two clouds collide along the x-axis at t=0 and x=0 with velocity $20\;{\rm km\;s^{-1}}$, i.e., the mach number is 2.17. The cooling function in this temperature region fitted by Koyama & Inutsuka (2002) as follows:

\begin{displaymath}\rho{\cal L} = n (-\Gamma + n\Lambda)\;\;{\rm erg\;cm^{-3}\;s^{-1}}, \end{displaymath} (57)

where

\begin{displaymath}\Gamma = 2\times10^{-26}, \end{displaymath} (58)

\begin{displaymath}\frac{\Lambda}{\Gamma} = 1.0\times10^7 \exp\left( -\frac{11... ... + 1.4\times10^{-2}\sqrt{T}\exp\left( -\frac{92}{T} \right), \end{displaymath} (59)

where n is the number density of the gas. We perform numerical calculation using the one-dimensional Lagrangian Godunov scheme (van Leer 1997) between x=0 and x=Lx = 3.26 pc. At x=0 and x=Lx, we adopt the reflecting and free boundary conditions, respectively. As the initial condition, we add the following density fluctuation:

\begin{displaymath}\rho(t=0,x)=\frac{\rho_{\rm WNM}} {\displaystyle 1+\frac{A}{... ...\rm max}} \sin\left( \frac{2\pi }{L_x} +\theta_i\right) }, \end{displaymath} (60)

where we set $i_{\rm max}=10,\;A=0.5$, and the phase, $\theta_i$, is given by random number between 0 and $2\pi$. The maximum and minimum number density at the initial state are 0.72 cm-3 and 0.44 cm-3, respectively.

\begin{figure} \par\includegraphics[width=6cm,clip]{12806fg7.eps} \end{figure} Figure 7:

Evolutionary track of the gas at the centre after the head-on collision (the thick solid line). The ordinate and abscissa axes indicate the pressure and the number density, respectively. The thin solid line indicates the thermally equilibrium curve. The filled circle represents the gas at the preshock region.

Open with DEXTER

\begin{figure} \par\includegraphics[width=6cm,clip]{12806fg8.eps} \end{figure} Figure 8:

Time evolution of $\alpha _{\rm num}$ (the solid line) and the cooling length (the dotted line), which are evaluated at the centre . The dashed line indicates the time evolution of the critical wave length defined in Sect. 5.3.

Open with DEXTER

\begin{figure} \par\includegraphics[width=12cm,clip]{12806fg9.eps} \end{figure} Figure 9:

Time evolution of a) the number density and b) the temperature, respectively. The rescaled c) number density, $nt_*^{\eta /(2-\alpha )}$ and d) temperature, $T t_* ^{(1-\alpha )(1-\eta )/\{\eta (2-\alpha )\}}$ as a function of the rescaled coordinate, $xt_*^{-1/(1-\omega )}$, where $t_{\rm c}$ and $\eta $ are set to 0.915 Myr and 0.98, respectively. The thick lines in c) and d) indicate the corresponding S-S solution.

Open with DEXTER

Figure 7 shows the evolutionary track of the gas at the centre, x=0. Due to the shock compression, the temperature of gas increases suddenly, and the postshock region enters the thermally unstable phase from the initial stable state. The TI leads to the cooling layer condensing in a runaway fashion. The cooling length in Fig. 8 rapidly decreases with time until it reaches minimum value $1.5\times10^{-3}$ pc at t=0.917 Myr. During runaway cooling, the gas condenses isobarically until it reaches the stable CNM phase, although some pressure oscillation is seen in Fig. 7. We focus on the property of this runaway condensing layer. First we evaluate $\alpha_{\rm num}(t)$, which is defined by

\begin{displaymath}\alpha_{\rm num}=\left( \frac{ \partial (\rho{\cal L})}{\partial T} \right) _{\rm\rho}\Biggr\vert _{x=0} \end{displaymath} (61)

at each instant of time. Figure 8 shows the time evolution of  $\alpha _{\rm num}$. During $0.4<t/{\rm Myr}<0.9$, it is seen that $\alpha_{\rm num}\sim0.61$ is approximately constant. Therefore, in this period, the flow is expected to be approximated well by the S-S solutions.

Figures 9 show the time evolutions of the distributions of (a) the number density and (b) the temperature at t= 0.72, 0.80, 0.88, and 0.90 Myr, respectively. The rescaled number density $nt_*^{\eta /(2-\alpha )}$ and temperature, $T t_*^{\left\{ \eta - (2-\alpha) \right\}/\{(2-\alpha)(1-\alpha)\}}$are shown in Figs. 9c and 9d as a function of rescaled coordinate, $xt_*^{-1/(1-\omega )}$, where $t_{\rm c}$, $\eta $, and $\alpha$are set to 0.915 Myr, 0.98, and 0.61, respectively. From Figs. 9c and 9d, it is clearly seen that the S-S solution $(\alpha=0.61, \eta=0.98)$ describes the results of the numerical calculation well. There are two reasons they are described by the S-S solution well. One is that $\alpha _{\rm num}$ is almost constant, and the cooling function is approximated well by Eq. (1) during the runaway condensation. The other reason relates to the stability of S-S solutions. The one-dimensional calculation solves the evolution only in the direction parallel to the condensation. From the linear analysis in Sect. 3, such a perturbation with k=0 cannot grow enough during the runaway condensation. That is why the S-S solution is realized even though the density is initially fluctuated.

However, the result is expected to be qualitatively different if the perturbation with $k\ne 0$ exists. The perturbation with $k\ne 0$ on the isobarically runaway condensing layer grows as $\propto t_*^{-1}$ (see Sect. 4). The growth rate is independent of wave numbers. From Figs. 8 and 9, the gas evolves obeying the S-S solution between t= 0.4 and 0.9 Myr. If the layer has a perturbation ($k\ne 0$) with amplitude $\delta \rho_0$ at $t=0.4\;{\rm Myr}$, the perturbation grows as

\begin{displaymath}\frac{\delta \rho(t)}{\delta \rho_0}=\frac{1-0.4~ {\rm Myr}/t_{\rm c}} {1 - t/t_{\rm c}}{\cdot} \end{displaymath} (62)

The epochs when the perturbation grows by a factor of 10 and 100 are t= 0.86 and 0.91 Myr, respectively. At 0.86 Myr, the central number density is still as low as 38 cm-3. Therefore, the condensing layer is expected to fragment quickly, and CNM clumps will form.

Inoue & Inutsuka (2008) investigated TI in a shock compressed region formed by WNM-WNM collision using two-dimensional, two-fluid magnetohydrodynamical simulation. The initial condition there is the same as that in our one-dimensional simulation except for the dimension and the way to add density fluctuation in the WNM. Our linear analysis is applicable to the gas evolution before CNM clouds form in their unmagnetized case. In the initial stage of their simulation, in the postshock region, thermally unstable gas initially condenses in a layer structure. Around t=0.5 Myr when the region of the highest density reaches CNM stable state, both small CNM cloudlets and filamentary CNM clouds are generated (Inoue 2009, priv. comm.). According to our linear analysis, an isobarically condensing gas layer is expected to split into fragments that have a variety of lengths in the transverse direction. Therefore, our linear analysis agrees with their simulation qualitatively.

After the first CNM clouds formation, the mass of these CNM clouds continues to increase by the accretion of unstable gas by the shocks. Moreover, some CNM clouds coalesce into larger clouds and others fragment into smaller ones by the turbulent flow (Inoue 2009, priv. comm.). Therefore, the size and mass of CNM clouds varies with time. To understand this phase, which is beyond the scope of this paper, a statistical approach is probably needed. A statistical theory has been proposed by Hennebelle & Audit (2007) using Press-Schechter formalism (Press & Schechter 1974), assuming that the CNM clumps are generated from density fluctuations within WNM whose power spectrum is assumed to be Kolmogorov.

Comparing timescales, Heitsch et al. (2008) discuss the effect of TI, the nonlinear thin-shell instability, Kelvin-Helmholtz instability, and gravitational instability in various densities, temperatures, and scales of fluctuations. In their paper, it is assumed that the gas is unstable only if the scale of the fluctuation is smaller than the sound crossing scale. However, from our linear analysis, perturbations can grow regardless of their scales as long as they are flat rather than spherical.

5.2 The growth rate for ${1<\alpha <2}$

Although the linear analysis on the S-S cooling layer is limited for $\alpha<1$, the thermal stability of the gas for $\alpha>1$ is also roughly understood from Balbus's criterion. For ${1<\alpha <2}$, the gas is isobarically unstable, but it is isochorically stable. For $\alpha>2$, the gas is thermally stable. In this section, we investigate the stability of the gas for ${1<\alpha <2}$ during cooling within the large and small scale limits.

5.2.1 The isobaric mode
For the case with small wave length, perturbation is expected to grow isobarically. By comparison of our results with previous studies in the literature, it is found that the growth rate in the isobaric mode is independent of the global structure of the unperturbed state. Therefore, from local arguments, the growth rate in the isobaric mode of the gas with ${1<\alpha <2}$can also be estimated.

As an unperturbed state, we adopt a cooling gas element whose scale is assumed to be much smaller than the cooling length. In this case, the element cools isobarically. From Eq. (4), the time evolution of the unperturbed gas is given by

    $\displaystyle \rho(t)=\rho_i\left( 1-\frac{t}{t_{\rm cool}'} \right)^{-1/(2-\alpha)},$  
    $\displaystyle \frac{1}{t_{\rm cool}'} = (2-\alpha)\gamma^{\alpha-1}(\gamma-1) P_i^{\alpha-1}\rho_i^{2-\alpha},$ (63)

where $\rho_i$ and Pi represent the initial density and pressure, respectively. In the above unperturbed state, we consider the following isobaric perturbation:

\begin{displaymath}\rho = \rho_0(t) +\delta \rho(t), \end{displaymath} (64)

and

P = P0, (65)

where subscript ``0'' indicates the unperturbed state, and $\delta \rho$ is the density perturbation. Linearizing Eq. (4), one obtains

\begin{displaymath}\frac{{\rm d} }{{\rm d}t}\left( \frac{\delta\rho}{\rho_0} \ri... ...ho_0^{2-\alpha}P_0^{\alpha-1}\frac{\delta \rho}{\rho_0}{\cdot} \end{displaymath} (66)

Using Eq. (63), Eq. (66) is rewritten as

\begin{displaymath}\frac{{\rm d} }{{\rm d}t}\left( \frac{\delta\rho}{\rho_0} \ri... ..._{\rm cool}'} \right)^{-1} \frac{\delta \rho}{\rho_0}{\cdot} \end{displaymath} (67)

Equation (67) is easily integrated to give

\begin{displaymath}\frac{\delta \rho}{\rho_0} \propto \left( 1-\frac{t}{t_{\rm cool}'} \right)^{-1}{\cdot} \end{displaymath} (68)

Comparing Eq. (68) with Eq. (63), one can see that the perturbation grows more slowly than the unperturbed state for ${1<\alpha <2}$. Therefore, the gas is expected to be difficult to fragment during runaway cooling if ${1<\alpha <2}$.

5.2.2 The noninteractive mode
A cooling layer that evolves isobarically is considered. The time evolution of the central density is the same as Eq. (63). When the scale of perturbation perpendicular to the condensation is too large to interact with other regions, each region evolves independently. Here, we focus on the time evolution of density perturbation at the centre (x=0). Initial fluctuation of the central density, $\delta \rho_i$, creates the fluctuation of the cooling time, $\Delta t$. The relative amplitude of density perturbation at x=0 is given by

\begin{displaymath}\frac{\delta \rho}{\rho_0} = \frac{1}{\rho_0}\left( \rho_i + ... ...frac{t}{t_{\rm cool}' - \Delta t} \right)^{-1/(2-\alpha)}-1. \end{displaymath} (69)

Linearizing Eq. (69) with omitting terms that do not grow, we have

\begin{displaymath}\frac{\delta \rho}{\rho_0} = \frac{1}{2-\alpha}\frac{\Delta t... ... cool}'} \left( 1-\frac{t}{t_{\rm cool}'} \right)^{-1}{\cdot} \end{displaymath} (70)

Comparing Eq. (70) with Eq. (63), we can see that the perturbation grows more slowly than the unperturbed state for ${1<\alpha <2}$. Therefore, the gas is expected to be difficult to fragment for ${1<\alpha <2}$ for the large-scale perturbation, as well as the small scale.

5.3 Effects of thermal conduction

In this paper, the thermal conduction is neglected for simplicity. However, for large wave number, the thermal conduction is expected to stabilize TI in the cooling layer (Field 1965). Therefore, there is a critical wave number, $k_{\rm crit}$, such that perturbation with a larger wave number is stabilized by the thermal conduction.

First, we evaluate $k_{\rm crit}$ using an order estimation. Using the characteristic time scale of the thermal conduction, $t_{\rm diff}$, the diffusion equation is given by

\begin{displaymath}\frac{1}{t_{\rm diff}} \left( \frac{P_{00}}{\gamma-1} \right) \sim k^2K(T_{00}) T_{00}, \end{displaymath} (71)

where K is the thermal conduction coefficient. From Eq. (71), the diffusion timescale is given by

\begin{displaymath}t_{\rm diff} \simeq \frac{P_{00}}{(\gamma-1)K(T_{00})T_{00} } k^{-2}. \end{displaymath} (72)

From Eq. (72), one can see that the diffusion timescale is small for large wave number. If $t_{\rm diff}<t_{\rm cool}$, the thermal conduction is expected to stabilize TI. Therefore, $k_{\rm crit}$ can be derived on the condition $t_{\rm diff}\sim t_{\rm cool}$ as

\begin{displaymath}k_{\rm crit} = \sqrt{\frac{k_{\rm K}}{L_{\rm cool}}},\;\;{\rm... ...}\; k_{\rm K} = \frac{P_{00}c_{00}}{(\gamma-1)KT_{00}}{\cdot} \end{displaymath} (73)

Field (1965) derived similar critical wave number to Eq. (73). Since the unperturbed state is time dependent, $k_{\rm crit}$ also evolves with time. Detailed evolution of $k_{\rm crit}$ depends on K. In T<6000, we adopt $K=2.5\times10^3\sqrt{T}$ ergs cm-1 K-1 s-1(Parker 1953). In this case, from Eq. (73), the time evolution of $k_{\rm crit}$ can be derived analytically as
                    $\displaystyle \frac{{\rm d} \ln k_{\rm crit}}{{\rm d}\ln t_*}$ = $\displaystyle \frac{(2\alpha-1)\eta - (2-\alpha)(3-2\alpha)}{4(2-\alpha)(1-\alpha)}$  
  = $\displaystyle \left\{\begin{array}{cc} \displaystyle - \frac{7-2\alpha}{4(2-\al... ...le - \frac{3-2\alpha}{4(1-\alpha)}<0 & {\rm for}\;\eta=0.\\ \end{array}\right.$ (74)

For $\eta=1$ and $\eta=0$, it is found that Eq. (74) is negative for $\alpha<1$. Since Eq. (74) is the linear function for $\eta $, ${\rm d}\ln k_{\rm crit}/{\rm d}\ln t_*$ is negative for all $\eta $. Therefore, $k_{\rm crit}$ increases with time. This means that an initially stable perturbation with wave number ( $k>k_{\rm crit}$) becomes unstable at a certain epoch when $k_{\rm cirt}$ catches up with k.

Figure 8 shows the time evolution of the critical wave length, $\lambda_{\rm crit}=2\pi/k_{\rm crit}$. In Sect. 5.1, the time evolution of the condensing layer can be described by the S-S solution with ( $\alpha_{\rm num}=0.61$, $\eta=0.98$). In Fig. 8, we see that $\lambda_{\rm crit}$ decreases with time during the runaway condensation. Figure 8 also shows that $\lambda_{\rm crit}<L_{\rm cool}$, or $\kappa_{\rm crit}= k_{\rm crit}L_{\rm cool}>1$. This means that the effect of thermal conduction on the dispersion relation (Fig. 6) always appears in the isobaric regime during the runaway condensation.

6 Summary

In this paper, we have investigated the stability of S-S solutions describing the runaway cooling of a radiative gas by linear analysis. The results of our investigation are summarized as follows,

1.
For the case with perturbation only parallel to the flow (k=0), the S-S solutions are unstable. However, the growth rate is too low to become nonlinear during the runaway cooling. Actually, the S-S solutions are realized in one-dimensional hydrodynamical calculations in IT08 and in Sect. 5.1 in this paper.

2.
For the case with transverse perturbation ($k\ne 0$), there are several unstable modes in the S-S solutions. The most unstable modes are the isobaric mode for $k\gg1$ and the noninteractive mode for $k\ll1$. In the isobaric mode, the perturbation grows in pressure equilibrium with its surroundings. On the other hand, the noninteractive mode is originated from each region in the layer condensing independently. Under a static approximation, we derive the approximated dispersion relation. The results of direct numerical integration of the time evolution agree with those using the static approximation.

3.
The S-S solutions for $\eta>\alpha$ are unstable for any wavelength. Especially, if the unperturbed state is isobaric, the growth rate is independent of wave number. Therefore, fluctuations in various scales grow simultaneously, and the gas layer is expected to split into fragments with various scales. The S-S solutions for $\eta < \alpha$ are only unstable in the isobaric mode.
Our linear analysis predicts that the cooling layer splits into fragments quickly even if the size is greater than the local cooling length. Our linear analysis is qualitatively consistent with the results of recent multi-dimensional numerical simulations until CNM clouds form, but the evolution of CNM clouds is beyond the scope of this paper.

Acknowledgements
We would like to thank Tsuyoshi Inoue for valuable discussions. We also would like to thank the anonymous referee for valuable comments and suggestions that improved this paper significantly. K.I. is supported by grants-in-aid for JSPS Fellow (21-1979).

Appendix A: Supplements for derivation of basic equations in the zooming coordinate

Here, we present preparation for derivation of basic Eqs. (10)-(12) in the zooming coordinate given by Eq. (5). In the zooming coordinate, the physical variables, Q(t,x,y), are given by the following unified form:

\begin{displaymath}Q(t,x,y) = Q_0(t)\Theta(\tau,\xi,y),\;\;\; Q_0(t) \propto t_*^{-q/(1-\omega)}, \end{displaymath} (A.1)

where Q(t,x,y) corresponds to $[\rho, v_x, v_y, P]$ (see Eq. (6)), and $\Theta=[\Omega, V_x, V_y, \Pi]$ are the physical variables in the zooming coordinate. From Eqs. (7)-(9), a parameter, q, is given by

\begin{displaymath}q = \left\{ \begin{array}{ll} -\beta & {\rm for}\;\;Q=\rho ... ..._y\\ -2\omega -\beta& {\rm for}\;\;Q=P.% \end{array} \right. \end{displaymath} (A.2)

The temporal and spatial derivatives of Q(t,x,y) in the ordinary coordinate can be expressed in the zooming coordinate as
                            $\displaystyle \left( \frac{\partial Q}{\partial t} \right)_{x,y}$ = $\displaystyle \frac{{\rm d} Q_0(t)}{{\rm d} t} \Theta + Q_0(t) \frac{{\rm d}\ta... ...\frac{\partial \xi}{\partial t}\right)_{x} \frac{\partial \Theta}{\partial \xi}$  
  = $\displaystyle \frac{v_0(t)}{x_0(t)}Q_0(t) \left[ q \Theta + \frac{\partial \Theta}{\partial \tau} + \xi\frac{\partial \Theta}{\partial \xi} \right],$ (A.3)

\begin{displaymath}\left( \frac{\partial Q}{\partial x} \right)_{t,y} = \frac{Q... ...} \right)_{t,x} = Q_0(t) \frac{ \partial \Theta}{\partial y}, \end{displaymath} (A.4)

respectively, where we use

\begin{displaymath}\frac{v_0(t)}{x_0(t)}= -\frac{\dot{x}_0(t)}{x_0(t)}= \frac{1}{(1-\omega)t_{\rm c}t_*} \end{displaymath} (A.5)

from Eq. (7). Using Eqs. (A.3) and (A.4), the Lagrangian time derivative of Q is given by

    $\displaystyle \left(\frac{\partial }{\partial t} + v_0(t) V_x \frac{\partial }{\partial x} + v_0(t) V_y \frac{\partial }{\partial y} \right)Q(t,x,y) =$  
    $\displaystyle ~\frac{v_0(t)}{x_0(t)}Q_0(t) \left[ q + \frac{\partial}{\partial ... ...rtial \xi} + V_y x_0(t)\frac{\partial }{\partial y} \right] \Theta(\tau,\xi,y).$ (A.6)

Using Eqs. (A.3)-(A.6),

\begin{displaymath}\frac{v_0(t)}{x_0(t)} = \gamma^{\alpha-1} \Lambda_0 \rho_0(... ...},\;\;{\rm and}\;\; P_0(t)= \frac{\rho_0(t)v_0(t)^2}{\gamma}, \end{displaymath} (A.7)

basic Eqs. (10)-(12) are derived in the zooming coordinate.

Appendix B: Detailed expression of A $_\textit{\tiny {ik}}$

In this appendix, we provide the detail expression of Aik as

\begin{displaymath}A_{11}= -\left( \frac{V_0^2-X_0^2}{V_0} \right)\sigma - \fr... ...{\gamma}(\ln\Pi_0)' + \frac{X_0^2}{V_0}(\alpha-2)\epsilon_0 \end{displaymath} (B.1)

\begin{displaymath}A_{12}= -V_0(\ln \Omega_0)' + \sigma - \omega + V_0' - \frac{X_0^2}{\gamma V_0}\left(\ln\Pi_0\Omega_0^{-\gamma}\right)', \end{displaymath} (B.2)

\begin{displaymath}A_{13}= \frac{X_0^2}{\gamma} \left[ (\ln\Pi_0)' - \frac{\sigma + (\alpha-1)\gamma\epsilon_0}{V_0} \right], \end{displaymath} (B.3)

A14= - kx0(t)V0, (B.4)

\begin{displaymath}A_{21} = \frac{X_0^2}{\gamma}\left( V_0(\ln\Pi_0)' - (\alpha-2)\gamma\epsilon_0 \right), \end{displaymath} (B.5)

\begin{displaymath}A_{22} = X_0^2(\ln \Omega_0)' - V_0(\sigma - \omega + V_0') + \frac{X_0^2}{\gamma}\left(\ln\Pi_0\Omega_0^{-\gamma}\right)', \end{displaymath} (B.6)

A23= -V0A14, (B.7)

A24= kx0(t)X02, (B.8)

\begin{displaymath}A_{31} = - X_0^2(\ln\Pi_0)' + (\alpha-2)V_0\gamma\epsilon_0, \end{displaymath} (B.9)

\begin{displaymath}A_{32} = \gamma\left\{-V_0(\ln \Omega_0)' + \sigma - \omega +... ...{V_0}{\gamma}\left(\ln \Pi_0\Omega_0^{-\gamma}\right)\right\}, \end{displaymath} (B.10)

\begin{displaymath}A_{33}= X_0^2(\ln\Pi_0)' - V_0\{\sigma + (\alpha-1)\gamma\epsilon_0\}, \end{displaymath} (B.11)

\begin{displaymath}A_{34}= - \gamma kx_0(t)V_0, \end{displaymath} (B.12)

A41= 0, (B.13)

A42= 0, (B.14)

\begin{displaymath}A_{43}= - \frac{V_0^2 - X_0^2}{V_0}(\omega + \sigma), \end{displaymath} (B.15)

and

\begin{displaymath}A_{44}= kx_0(t)(V_0^2-X_0^2)\frac{X_0^2}{\gamma}. \end{displaymath} (B.16)

References

All Figures

  \begin{figure} \par\includegraphics[width=12cm,clip]{12806fg1.eps} \end{figure} Figure 1:

Growth rate, $\Sigma $, as a function of $\eta $ for a)  $\alpha =-1.0$ and b) 0.5 for the case with the perturbation parallel to the condensation. For comparison, the increasing rate of the unperturbed central density, $\alpha _\rho (\eta )$, and the decreasing rate of the unperturbed central pressure, $\alpha _P(\eta )$, are shown by the dashed and dotted lines, respectively.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=6cm,clip]{12806fg2a.eps} \vspace*{2mm}... ...ps} \vspace*{2mm} \includegraphics[width=6cm,clip]{12806fg2c.eps} \end{figure} Figure 2:

Approximate dispersion relations for $\alpha =0.5$ with a) $\eta =0.1$, b) 0.75, and c) 0.95, where $\kappa $ is the nondimensional wave number, and $\Sigma $ the growth rate of the perturbation. The labels of number represent the branches of modes.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=6cm,height=4cm,clip]{12806fg3.eps} \end{figure} Figure 3:

Eigenfunctions for $\eta =0.75$ and $\kappa =8$. The filled circles indicate the values at the critical point. The eigenvalue, $\Sigma $, is equal to 1.196.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=12cm,clip]{12806fg4.eps} \end{figure} Figure 4:

Growth rate as a function of $\eta $ for a) $\alpha =0.5$ and b) 0.8 for the case with $k\ne 0$. The thick solid lines correspond to the growth rate of the isobaric mode, $\Sigma _{\rm isobaric}$. For comparison, the evolutionary rate of the unperturbed central density, $\alpha _\rho (\eta )$, and pressure, $\alpha _P(\eta )$ are shown by the dashed and dotted lines, respectively. The thin solid lines pointed by arrows correspond to the growth rate in the noninteractive mode, $\Sigma _{\rm non-int}$.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=7.4cm,clip]{12806fg5.eps} \end{figure} Figure 5:

Schematic picture of the noninteractive mode.

Open with DEXTER
In the text

  \begin{figure} \mbox{\includegraphics[width=5.8cm,clip]{12806fg6a.eps}\hspace*{4mm} \includegraphics[width=5.8cm,clip]{12806fg6b.eps} }\end{figure} Figure 6:

Growth rate obtained by results of numerical linear analysis for $\eta =0.1$ a) and 0.75  b). The thick and dashed lines indicate the results of numerical linear analysis and those of approximate linear analysis.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=6cm,clip]{12806fg7.eps} \end{figure} Figure 7:

Evolutionary track of the gas at the centre after the head-on collision (the thick solid line). The ordinate and abscissa axes indicate the pressure and the number density, respectively. The thin solid line indicates the thermally equilibrium curve. The filled circle represents the gas at the preshock region.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=6cm,clip]{12806fg8.eps} \end{figure} Figure 8:

Time evolution of $\alpha _{\rm num}$ (the solid line) and the cooling length (the dotted line), which are evaluated at the centre . The dashed line indicates the time evolution of the critical wave length defined in Sect. 5.3.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=12cm,clip]{12806fg9.eps} \end{figure} Figure 9:

Time evolution of a) the number density and b) the temperature, respectively. The rescaled c) number density, $nt_*^{\eta /(2-\alpha )}$ and d) temperature, $T t_* ^{(1-\alpha )(1-\eta )/\{\eta (2-\alpha )\}}$ as a function of the rescaled coordinate, $xt_*^{-1/(1-\omega )}$, where $t_{\rm c}$ and $\eta $ are set to 0.915 Myr and 0.98, respectively. The thick lines in c) and d) indicate the corresponding S-S solution.

Open with DEXTER
In the text

Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.