A&A 508, 725735 (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 5600043, 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 timedependent, selfsimilar 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 ,
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
at the centre, where
the cooling length,
,
represents the length that a sound wave can travel
during the cooling time. For
,
the perturbation grows isobarically.
For
,
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 lowtemperature phase (cold neutral medium, or CNM) and a diffuse hightemperature 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 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 multidimensional 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 shockcompressed region. Analogous processes have been studied by many authors for a colliding flow of the WNM (Audit & Hennebelle 2005; VazquezSemadeni et al. 2007; Heitsch et al. 2006; Hennebelle & Audit 2007), and using twofluid MHD simulation (Inoue & Inutsuka 2008). The unbalance between cooling and heating rates causes the shockcompressed 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 KelvenHelmholtz 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 selfsimilar (SS) solutions describing the dynamical condensation of a radiative gas layer where the cooling rate dominates the heating rate. This SS solution assumed that the net cooling rate is a powerlaw function and that the heating rate is explicitly neglected. Although it is still too ideal, they are expected to be a good nonlinear onedimensional model at least in the phase during the transition from WNM to CNM. In this paper, we adopt the SS solutions as a more realistic unsteady unperturbed state than those in previous works. We perform linear analysis of the SS 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 SS unperturbed state, the nonlinear thin shell instability, KelvinHelmholtz 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 SS solutions. In Sect. 3, we investigate the stability of the SS 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:In ISM, in the temperature range of , the main coolant is Bremsstralung for the solar metallicity (Sutherland & Dopita 1993). The cooling rate of Bremsstralung is approximated well by that with . In the temperature range of , the dominant coolant is metal lines for the solar metallicity (Sutherland & Dopita 1993). The cooling rate of metal lines corresponds to that with . However, the cooling rate cannot be expressed by a single . In the temperature range between CNM and WNM, , the dominant coolant is CII (Dalgarno & McCray 1972). In this case, the cooling rate is approximated by that with as is shown in Sect. 5.1.
Basic equations for a radiative gas are
the continuity equation,
the equation of motion,
and the entropy equation,
where indicates the Lagrangian time derivative.
We take the xaxis as the direction of
the condensation driven by the cooling and
yaxis as the transverse direction.
Since the SS solutions are timedependent, it
is difficult to perform linear analysis
in the ordinary Cartesian coordinate, (t,x,y).
Bouquet et al. (1985) introduced a zooming coordinate
where SS 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:
where is a free parameter, a corresponds to the cooling length, and , is an epoch when the central density becomes infinity. In the zooming coordinate, density , velocity , pressure , and sound speed X are given by
respectively, where
and
with .
In the zooming coordinate,
the basic Eqs. (2)(4) are rewritten as
and
respectively, where the operators of time and spatial derivative are defined by
respectively, where indicates the unit vector parallel to direction. Supplements for the derivation of Eqs. (10)(12) is presented in Appendix A.
We apply the zooming transformation only in the xdirection but not in the ydirection. This is because the gas contracts along xaxis but not along yaxis 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 ydirection, the transverse scale of the perturbation decreases with time in the ordinary coordinate, although the unperturbed gas does not contract along the yaxis. Therefore, we apply the zooming transformation only in the xdirection.
2.1 Review of selfsimilar solutions
In the zooming coordinate, steady state solutions correspond to SS solutions that were derived in IT08. In this section, physical properties of the SS solutions are reviewed briefly.The SS solutions are specified by two parameters,
and .
For convenience, instead of ,
we can use a parameter ,
which
is given by
(14) 
Using these parameters ( ), the time dependences of the central density, , and pressure, P_{00}, are given by
respectively, where and .
The SS solutions include two asymptotic solutions.
For ,
the time dependences of the central density and pressure are
given by
(16) 
respectively. This time evolution indicates the isochoric mode. For , the time dependences of the central density and pressure are given by
(17) 
respectively. This time evolution corresponds to the isobaric mode. Our SS solutions exist between the two limits, . With the condition that the increasing rate of is equal to the decreasing rate of P_{00}(t), or , the critical value of is given by
(18) 
Therefore, the SS solutions for and are close to the isochoric and the isobaric modes, respectively. During SS condensation, the central pressure should not increase with time. From Eq. (15), the range of the value of is found to be .
2.2 Perturbation equations
Perturbation on the SS solutions is considered. Perturbed variables are defined bywhere subscript ``0'' indicates the unperturbed state.
As the perturbation, we consider the following Fourier mode with respect to y,
and k indicates the wave number of the plane wave that propagates along ydirection. Substituting Eqs. (19) and (20) into Eqs. (10)(12) and linearizing, we get the following perturbation equations:
and
where
,
(25) 
The timedependent factors remain in the form of 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 in general.
Figure 1: Growth rate, , as a function of for a) and b) 0.5 for the case with the perturbation parallel to the condensation. For comparison, the increasing rate of the unperturbed central density, , and the decreasing rate of the unperturbed central pressure, , 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 timedependent factor,
,
vanishes,
the perturbed variables can be expanded in the Fourier mode with respect to as
By Eq. (26), the time evolution of the perturbations is given by
Substituting Eq. (26) into the perturbation Eqs. (21)(24), one obtains the following ordinary differential equations:
where the detailed expression of A_{ij} 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 and at the critical point, , where V_{0}=X_{0}. The boundary conditions at are obtained by the asymptotic limit of the perturbed variables. From the regularity of the perturbed variables at , we find that perturbed variables should have the following asymptotic forms:and
where .
The boundary conditions at the critical point,
,
are obtained
from Eqs. (28).
At
,
the denominator of the righthand side becomes zero.
To obtain a regular solution from
to
,
the numerator of the righthand side should vanish.
Therefore, the boundary conditions are given by the following three equations,
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 and impose the boundary condition at , the solution is completely fixed. In this section, our numerical method for solving Eqs. (28) is described.We can set without loss of generality. For a given , we integrate Eqs. (28) from to the critical point, , using a fourthorder RungeKutta method. Eigenvalue, , is modified until the perturbed variables satisfy the boundary condition at using the NewtonRaphson method. After that, we integrate Eqs. (28) up to .
3.3 Results
Figure 1 shows the dependence of growth rate, , on for (a) and (b) 0.5. From Fig. 1, it is seen that for a wide range of and . Therefore, the perturbation is unstable. However, is smaller than and . Therefore, the growth rate is too low to grow sufficiently during runaway condensation. This is consistent with the results of the onedimensional numerical simulation shown in Sect. 5.1.4 Perturbation with k 0
For , the perturbation equations contain the timedependent factor, . Here we show that this factor, , 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(32) 
where superscript ``0'' indicates the value at t=0. Therefore, using , the collapse time can be expressed as
Using Eqs. (6) and (7), at t=0, the sound speed at the centre is given by
(34) 
Therefore, x_{0}(0) is given by
Using Eq. (33), Eq. (35) can be written as
(36) 
where , and is the cooling length at the centre at t=0. Since the SS solutions are invariant under the transformation , , , and , with a parameter m, one can take m in such a way that M_{00} is equal to 1. Hereafter, M_{00} is set to unity. Therefore, the timedependent factor, , is expressed as
where .
4.1 Static approximation
Before we present the fully timedependent numerical calculation
in Sect. 4.2,
we at first consider a special case in which the time evolutions of the unperturbed
SS 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 x_{0} to be a constant in Eqs. (21)(24).
This approximation is also valid for
where the term, ,
is negligibly small
in Eqs. (21)(24).
We use the Fourier mode as
The condition under which the static approximation is valid is given by
where the definition of is the same as that in Eq. (27). Substituting Eq. (38) into the perturbation Eqs. (21)(24), we get
where the detailed expression of A_{ij} is shown in Appendix B.
We impose the boundary conditions at
and
.
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 :
Substituting Eqs. (41) into Eqs. (40), we obtain the following relations:
and
The boundary condition at 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 with (a) , (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 limit, and (2) for small limit. Each branch is explained below.
Figure 2: Approximate dispersion relations for with a) , b) 0.75, and c) 0.95, where is the nondimensional wave number, and the growth rate of the perturbation. The labels of number represent the branches of modes. 

Open with DEXTER 
Figure 3: Eigenfunctions for and . The filled circles indicate the values at the critical point. The eigenvalue, , is equal to 1.196. 

Open with DEXTER 
(1) The isobaric mode
Figure 4: Growth rate as a function of for a) and b) 0.8 for the case with . The thick solid lines correspond to the growth rate of the isobaric mode, . For comparison, the evolutionary rate of the unperturbed central density, , and pressure, 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, . 

Open with DEXTER 
The branch (1) corresponds to the most unstable mode for , or . Since , 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 for the branch (1) is shown in Fig. 3 for and . In Fig. 3, it is clearly seen that . 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,
and
P=P_{00}(t).
From Eq. (4),
the following perturbation equation is obtained:
From Eqs. (8) and (9), we have
Using Eq. (47), Eq. (46) is rewritten as
Equation (48) can be integrated to give
(49) 
Therefore, the growth rate in isobaric mode is given by
In Fig. 2b, the eigenvalue, , is found to be 1.196 for the case with and . This is very close to for . From Fig. 2, it is clearly seen that the growth rate approaches the corresponding in the large 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 for 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
.
Figure 4 shows
,
and
as a function of
for (a)
and (b) 0.8.
From Fig. 4, it is found that
is greater than
for all .
For other ,
we can investigate analytically as follows:
subtracting
from
,
one obtains
=  
(51) 
Therefore, for , is more than , and the isobaric mode can grow in the runaway cooling layer.
Figure 5: Schematic picture of the noninteractive mode. 

Open with DEXTER 
(2) The noninteractive mode
Branch (2) corresponds to the most unstable mode for
.
Because the wave length is larger than the cooling length, each part evolves independently
according to the SS 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
and
.
Due to the difference of
,
the regions ``A'' and ``B'' have different collapse time,
and ,
respectively.
Omitting any terms that do not grow, we find the time evolution of difference,
,
to be
(52) 
Therefore, in the zooming coordinate, the density perturbation grows as . Other perturbed variables also grow in the same power law. Therefore, comparing with Eq. (27), the growth rate is given by
(53) 
which is independent of and . The noninteractive mode arises from the fluctuation of the collapse time, , 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 . Since the perturbation grows isobarically, is compared to . From Fig. 4, it is seen that the growth rate, , is higher than for all and 0.86 for and 0.8, respectively. Therefore, the noninteractive mode can grow sufficiently. Next, we consider the cooling layer with . In this layer, since the perturbation grows isochorically, is compared to . Analytically, it is found that the pressure perturbation can only grow for .
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.
Figure 6: Growth rate obtained by results of numerical linear analysis for 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
,
there is a solution where physical variables are very small except for
which is spatially constant, and the eigenvalue is
.
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
In Fig. 2, it is seen that each growth rate in branch (3) has the corresponding value of for . The physical meaning of this mode can be understood as follows: for , since the effect of the pressure gradient with respect to y is very weak, the gas can freely stream with almost constant velocity, v_{y}, in the y direction. On the other hand, the central sound speed, c_{00}(t), decreases as . Therefore, the ratio of the dynamical velocity to the thermal velocity, v_{y}/c_{00}(t), grows with time as , 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 freestreaming mode
For large ,
there is another mode in which
the velocity perturbation in the xdirection is
much greater than that in the ydirection,
.
We call this mode the freestreaming mode.
From Eq. (45) with
,
we obtain
(55) 
In the freestreaming mode, the growth of the velocity perturbation in the xdirection is hampered by the pressure gradient of the unperturbed state. Therefore, the growth rate is less than the shear mode.
(5) mode
The growth rate in this branch for coincides with the case with k=0, which is obtained in Sect. 3.
4.2 Linear analysis considering the time evolution of
The static approximation is valid only if is much larger than . However, from Fig. 2, it is found that the growth rate of the most unstable mode for each is smaller than whose values are 1.93, 1.5, and 1.37 for , 0.75, and 0.95 with , respectively (see Eq. (39)). Therefore, in this section, we perform a linear analysis considering the time evolution of 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 to in the zooming coordinate.We impose the boundary conditions at and . The even mode is set as the boundary condition at . The free boundary condition is set at , 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 which is obtained in Sect. 4.1.
By solving Eqs. (21)(24),
a time evolution of
is obtained.
During the calculation, the growth rate at
is evaluated by
(56) 
at each instant of time. For comparison with the result of the static approximation, we focus on a relation between and the growth rate of the density perturbation at the centre, . Figure 6 shows the growth rate as a function of at each instant of time for with (a) and (b) 0.75. The nondimensional wave number, , 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, , moderately agrees with that of the approximated dispersion relations. For , or initial phase, the growth rate agrees with . This is because the growth rate does not depend on in the short wave length limit. As decreases and reaches about 1, begins to decrease. For , approaches asymptotically , where the effect of is negligible. The effect of timedepending is notable only for . Smoother dependence of the growth rate on 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 headon collision between two thermally equilibrium gases in which density and pressure are and . Two clouds collide along the xaxis at t=0 and x=0 with velocity , i.e., the mach number is 2.17. The cooling function in this temperature region fitted by Koyama & Inutsuka (2002) as follows:(57) 
where
(58) 
(59) 
where n is the number density of the gas. We perform numerical calculation using the onedimensional Lagrangian Godunov scheme (van Leer 1997) between x=0 and x=L_{x} = 3.26 pc. At x=0 and x=L_{x}, we adopt the reflecting and free boundary conditions, respectively. As the initial condition, we add the following density fluctuation:
(60) 
where we set , and the phase, , is given by random number between 0 and . The maximum and minimum number density at the initial state are 0.72 cm^{3} and 0.44 cm^{3}, respectively.
Figure 7: Evolutionary track of the gas at the centre after the headon 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 
Figure 8: Time evolution of (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 
Figure 9: Time evolution of a) the number density and b) the temperature, respectively. The rescaled c) number density, and d) temperature, as a function of the rescaled coordinate, , where and are set to 0.915 Myr and 0.98, respectively. The thick lines in c) and d) indicate the corresponding SS 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
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
,
which is defined by
(61) 
at each instant of time. Figure 8 shows the time evolution of . During , it is seen that is approximately constant. Therefore, in this period, the flow is expected to be approximated well by the SS 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 and temperature, are shown in Figs. 9c and 9d as a function of rescaled coordinate, , where , , and are set to 0.915 Myr, 0.98, and 0.61, respectively. From Figs. 9c and 9d, it is clearly seen that the SS solution describes the results of the numerical calculation well. There are two reasons they are described by the SS solution well. One is that 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 SS solutions. The onedimensional 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 SS solution is realized even though the density is initially fluctuated.
However, the result is expected to be qualitatively different if the perturbation with
exists.
The perturbation with
on the isobarically runaway condensing layer
grows as
(see Sect. 4). The growth rate is independent of wave numbers.
From Figs. 8 and 9,
the gas evolves obeying the SS solution between t= 0.4 and 0.9 Myr.
If the layer has a perturbation ()
with
amplitude
at
,
the perturbation grows as
(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 WNMWNM collision using twodimensional, twofluid magnetohydrodynamical simulation. The initial condition there is the same as that in our onedimensional 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 PressSchechter 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 thinshell instability, KelvinHelmholtz 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
Although the linear analysis on the SS cooling layer is limited for , the thermal stability of the gas for is also roughly understood from Balbus's criterion. For , the gas is isobarically unstable, but it is isochorically stable. For , the gas is thermally stable. In this section, we investigate the stability of the gas for 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 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
where
and P_{i} represent the initial density and pressure, respectively.
In the above unperturbed state, we consider the following isobaric perturbation:
(64) 
and
P = P_{0},  (65) 
where subscript ``0'' indicates the unperturbed state, and is the density perturbation. Linearizing Eq. (4), one obtains
Using Eq. (63), Eq. (66) is rewritten as
Equation (67) is easily integrated to give
Comparing Eq. (68) with Eq. (63), one can see that the perturbation grows more slowly than the unperturbed state for . Therefore, the gas is expected to be difficult to fragment during runaway cooling if .
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, , creates the fluctuation of the cooling time, . The relative amplitude of density perturbation at x=0 is given byLinearizing Eq. (69) with omitting terms that do not grow, we have
Comparing Eq. (70) with Eq. (63), we can see that the perturbation grows more slowly than the unperturbed state for . Therefore, the gas is expected to be difficult to fragment for for the largescale 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, , such that perturbation with a larger wave number is stabilized by the thermal conduction.First, we evaluate
using an order estimation.
Using the characteristic time scale of the thermal conduction,
,
the diffusion equation is given by
where K is the thermal conduction coefficient. From Eq. (71), the diffusion timescale is given by
From Eq. (72), one can see that the diffusion timescale is small for large wave number. If , the thermal conduction is expected to stabilize TI. Therefore, can be derived on the condition as
Field (1965) derived similar critical wave number to Eq. (73). Since the unperturbed state is time dependent, also evolves with time. Detailed evolution of depends on K. In T<6000, we adopt ergs cm^{1} K^{1} s^{1}(Parker 1953). In this case, from Eq. (73), the time evolution of can be derived analytically as
For and , it is found that Eq. (74) is negative for . Since Eq. (74) is the linear function for , is negative for all . Therefore, increases with time. This means that an initially stable perturbation with wave number ( ) becomes unstable at a certain epoch when catches up with k.
Figure 8 shows the time evolution of the critical wave length, . In Sect. 5.1, the time evolution of the condensing layer can be described by the SS solution with ( , ). In Fig. 8, we see that decreases with time during the runaway condensation. Figure 8 also shows that , or . 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 SS 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 SS solutions are unstable. However, the growth rate is too low to become nonlinear during the runaway cooling. Actually, the SS solutions are realized in onedimensional hydrodynamical calculations in IT08 and in Sect. 5.1 in this paper.
 2.
 For the case with transverse perturbation (), there are several unstable modes in the SS solutions. The most unstable modes are the isobaric mode for and the noninteractive mode for . 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 SS solutions for 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 SS solutions for are only unstable in the isobaric mode.
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 grantsinaid for JSPS Fellow (211979).
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:(A.1) 
where Q(t,x,y) corresponds to (see Eq. (6)), and are the physical variables in the zooming coordinate. From Eqs. (7)(9), a parameter, q, is given by
(A.2) 
The temporal and spatial derivatives of Q(t,x,y) in the ordinary coordinate can be expressed in the zooming coordinate as
respectively, where we use
from Eq. (7). Using Eqs. (A.3) and (A.4), the Lagrangian time derivative of Q is given by
(A.7) 
basic Eqs. (10)(12) are derived in the zooming coordinate.
Appendix B: Detailed expression of A
In this appendix, we provide the detail expression of A_{ik} as(B.1) 
(B.2) 
(B.3) 
A_{14}=  kx_{0}(t)V_{0},  (B.4) 
(B.5) 
(B.6) 
A_{23}= V_{0}A_{14},  (B.7) 
A_{24}= kx_{0}(t)X_{0}^{2},  (B.8) 
(B.9) 
(B.10) 
(B.11) 
(B.12) 
A_{41}= 0,  (B.13) 
A_{42}= 0,  (B.14) 
(B.15) 
and
(B.16) 
References
 Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [EDP Sciences] [CrossRef]
 Balbus, S. A. 1986, ApJ, 303, L79 [NASA ADS] [CrossRef]
 Bouquet, S., Feix, R., & Fijalkow, E. 1985, ApJ, 293, 494 [NASA ADS] [CrossRef]
 Burkert, A., & Lin, D. N. C. 2000, ApJ, 537, 270 [NASA ADS] [CrossRef]
 Dalgarno, A., & MaCray, R. 1972, ARA&A, 10, 375 [NASA ADS] [CrossRef]
 Field, G. B. 1965, ApJ, 142, 531 [NASA ADS] [CrossRef]
 Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJ, 155, L149 [NASA ADS] [CrossRef]
 Hanawa, T., & Matsumoto, T. 1999, ApJ, 521, 703 [NASA ADS] [CrossRef]
 Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052 [NASA ADS] [CrossRef]
 Heitsch, F., Hartmann, L. W., & Burkert, A. 2008, ApJ, 683, 786 [NASA ADS] [CrossRef]
 Hennebelle, P., & Audit, E. 2007, A&A, 465, 431 [NASA ADS] [EDP Sciences] [CrossRef]
 Hennebelle, P., & Passot, T. 2006, A&A, 448, 1083 [NASA ADS] [EDP Sciences] [CrossRef]
 Inoue, T., & Inutsuka, S. 2008, ApJ, 687, 303 [NASA ADS] [CrossRef]
 Inoue, T., Inutsuka, S., & Koyama, H. 2006, ApJ, 652, 1331 [NASA ADS] [CrossRef]
 Iwasaki, K., & Tsuribe, T. 2008, MNRAS, 384, 1554 [NASA ADS] [CrossRef]
 Koyama, H., & Inutsuka, S. 2000, ApJ, 321, 980 [NASA ADS] [CrossRef]
 Koyama, H., & Inutsuka, S. 2002, ApJ, 564, L97 [NASA ADS] [CrossRef]
 McNamara, S. 1993, Phys. Fluids A, 5, 3056 [NASA ADS] [CrossRef]
 Parker, E. N. 1953, ApJ, 117, 431 [NASA ADS] [CrossRef]
 Press, W., & Schechter, P. 1974, ApJ, 187, 425 [NASA ADS] [CrossRef]
 Schwarz, J., McCray, R., & Stein, R. F. 1972, ApJ, 175, 673 [NASA ADS] [CrossRef]
 Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef]
 van Leer, B. 1997, J. Comput. Phys., 135, 229 [NASA ADS] [CrossRef]
 VazquezSemadeni, E., Gomez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870 [NASA ADS] [CrossRef]
 Vishniac, E. 1994, ApJ, 428, 186 [NASA ADS] [CrossRef]
 Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152 [NASA ADS] [CrossRef]
 Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 [NASA ADS] [CrossRef]
All Figures
Figure 1: Growth rate, , as a function of for a) and b) 0.5 for the case with the perturbation parallel to the condensation. For comparison, the increasing rate of the unperturbed central density, , and the decreasing rate of the unperturbed central pressure, , are shown by the dashed and dotted lines, respectively. 

Open with DEXTER  
In the text 
Figure 2: Approximate dispersion relations for with a) , b) 0.75, and c) 0.95, where is the nondimensional wave number, and the growth rate of the perturbation. The labels of number represent the branches of modes. 

Open with DEXTER  
In the text 
Figure 3: Eigenfunctions for and . The filled circles indicate the values at the critical point. The eigenvalue, , is equal to 1.196. 

Open with DEXTER  
In the text 
Figure 4: Growth rate as a function of for a) and b) 0.8 for the case with . The thick solid lines correspond to the growth rate of the isobaric mode, . For comparison, the evolutionary rate of the unperturbed central density, , and pressure, 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, . 

Open with DEXTER  
In the text 
Figure 5: Schematic picture of the noninteractive mode. 

Open with DEXTER  
In the text 
Figure 6: Growth rate obtained by results of numerical linear analysis for 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 
Figure 7: Evolutionary track of the gas at the centre after the headon 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 
Figure 8: Time evolution of (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 
Figure 9: Time evolution of a) the number density and b) the temperature, respectively. The rescaled c) number density, and d) temperature, as a function of the rescaled coordinate, , where and are set to 0.915 Myr and 0.98, respectively. The thick lines in c) and d) indicate the corresponding SS solution. 

Open with DEXTER  
In the text 
Copyright ESO 2009