Structure of dynamical condensation fronts in the interstellar medium

Abstract

1 INTRODUCTION

It is well known that in the interstellar medium (ISM), a clumpy low-temperature phase [(cold neutral medium (CNM)] and a diffuse high-temperature phase [warm neutral medium (WNM)] can coexist in pressure equilibrium as a result of the balance of radiative cooling and heating owing to external radiation fields and cosmic rays (Field, Goldsmith & Habing 1969; Wolfire et al. 1995, 2003). The CNM of atomic gas is observed as H i cloud (n∼ 10–100 cm⁻³, T∼ 10² K). These two phases are thermally stable. On the other hand, thermal instability (TI) arises in the temperature range between these two phases. Thus, the ISM in atomic phase can be interpreted as a bistable fluid.

Pioneering work on the TI has been done by Field (1965) who performed a linear analysis of the thermal equilibrium gas and derived a simple criterion for the TI. Focusing on a fluid element, Balbus (1986) derived a criterion for the TI of the thermal non-equilibrium gas. In the non-linear evolution, Iwasaki & Tsuribe (2008) found a family of self-similar solutions describing the condensation of radiative gas layer assuming a simple power-law cooling rate. Iwasaki & Tsuribe (2009) investigated the linear stability of the self-similar solutions, and they suggested that the condensing layer will fragment in various scales as long as the transverse scale is larger than the Field length. The physics of the bistable fluid has been investigated by many authors. Zel’dovich & Pikel’ner (1969, hereafter ZP69) and Penston & Brown (1970) investigated the structure of the transition front connecting the CNM and WNM under the plane-parallel geometry. The thickness of the transition front is characterized by the Field length below which the TI is stabilized by heat conduction (Field 1965). They found that a static solution is obtained at a so-called saturation pressure. The transition front becomes a condensation (evaporation) front if the surrounding pressure is larger (less) than the saturation pressure. These properties of the transition front depend on its geometry. In the spherical symmetric geometry, Graham & Langer (1973) investigated isobaric flows. They found that spherical clouds have a minimum size below which they inevitably evaporate (a more general description is found in Nagashima, Koyama & Inutsuka 2005). Linear stability of the transition layers has been investigated by Inoue, Inutsuka & Koyama (2006) in the plane-parallel geometry, and they found that the evaporation front is unstable against corrugation-type fluctuation, while the condensation front is stable. Stone & Zweibel (2009) have considered the effect of magnetic field on the instability. The physical mechanism of the corrugation instability is analogous to the Darriues–Landau instability (Landau & Lifshitz 1987). Taking into account the magnetic field perpendicular to the normal of the front, Stone & Zweibel (2010) investigated linear stability of transition layers.

In the above-mentioned theoretical works with respect to the bistable fluid, the phase transition proceeds in a quasi-static manner through heat conduction. However, recent multidimensional numerical simulations show a more dynamical turbulent structure. Kritsuk & Norman (2002) have investigated the non-linear development of the TI in three-dimensional hydrodynamical simulations with periodic boundary condition without any external forcing. As the initial condition, they set the hot ionized medium (T= 2 × 10⁶ K) where cooling dominates heating. The gas is quickly decomposed into two stable phases (CNM and WNM) and the intermediate unstable phase. In the multiphase medium, the supersonic turbulence is developed by conversion from the thermal energy to the kinetic energy. They found that the turbulence decays on a dynamical time-scale. This decaying timescale is longer than that in supersonic isothermal turbulence of one-phase medium where it is smaller or comparable to the flow crossing time (Stone, Ostriker & Gammie 1998). Koyama & Inutsuka (2006) have done similar calculations on a much longer time-scale. As the initial condition, they set an unstable gas in thermal equilibrium state with density fluctuations. They found a self-sustained turbulence with velocity dispersion of ∼0.2–0.4 km s⁻¹ for a period of at least their simulation time (several 10 Myr) in the bistable fluid after the temporal decaying found in Kritsuk & Norman (2002). External forcing such as shock compression can drive stronger long-lived turbulence. Koyama & Inutsuka (2002) investigated development of the TI in a shock-compressed region by using a two-dimensional hydrodynamical simulation. They found that the velocity dispersion is as large as several km s⁻¹, which is larger than the sound speed of the CNM. The supersonic turbulent motion is maintained as long as the shock wave continues to propagate and provide post-shock gas. They suggested that this supersonic translational motion of the CNM is observed as the supersonic turbulence in the ISM. The cloudlets are precursor of molecular clouds because they are as dense as molecular clouds. After their work, many authors have investigated the formation of the CNM in shocked gases in two-dimensional calculations (Audit & Hennebelle 2005; Heitsch et al. 2005, 2006; Hennebelle & Audit 2007) and in three-dimensional calculations (Gazol, Vázquez-Semadeni & Kim 2005; Vázquez-Semadeni et al. 2006; Audit & Hennebelle 2010). The influence of the magnetic field on the development of the TI has been investigated by Inoue & Inutsuka (2008, 2009), Hennebelle et al. (2008) and Heitsch, Stone & Hartmann (2009).

Although the physics of the bistable fluid is developing, the physical mechanism of driving turbulence is not fully understood. As mentioned above, recent numerical simulations have found that the turbulence involving TI is dynamical compared with the ZP69 picture that predicts a very slow transition between CNM and WNM through heat conduction. Moreover, in the shock-compressed region, because of its high pressure, the bistable fluid cannot exist. The ZP69 picture cannot be applied directly to the TI in the shock-compressed region.

In this paper, as a first step to understand the dynamical turbulent structure, we develop the works of ZP69 into a more dynamical transition front. In the steady solutions, there are two parameters: the pressure of the CNM and the mass flux across the front. ZP69 considered solutions only on a line in the parameter space. Elphick, Regev & Shaviv (1992) found steady solutions with various mass fluxes. However, they used a simple cubic function of the cooling rate that enables them to investigate analytically, and they assumed spatially constant pressure. We systematically search steady solutions in larger parameter space even in the high-pressure range where the bistable fluid does not exist by using a realistic cooling rate.

This paper is organized as follows. Basic equations and numerical methods are described in Section 2. In Section 3, we briefly review the solutions connecting the CNM and WNM derived by ZP69. In Section 4, we present solutions describing the dynamical transition layer that has larger mass flux than ZP69. The results are discussed in Section 5 and are summarized in Section 6.

2 BASIC EQUATIONS AND NUMERICAL METHODS

Basic equations for radiative gas under the plane-parallel geometry are the continuity equation,

1

the momentum conservation,

2

the energy equation,

3

and the equation of state,

4

where E=ρv²/2 +P/(γ− 1) is the total energy, γ= 5/3 is the ratio of specific heats, κ is the coefficient of heat conductivity, is the net cooling rate per unit mass, is the gas constant and m_H is the hydrogen mass. Throughout this paper we assume that a gas consists of atomic hydrogen. For the range of temperatures considered, since the gas is almost neutral, we adopt cm⁻¹ K⁻¹ s⁻¹ (Parker 1953). In this paper, we adopt the following fitting formulae of the net cooling rate (Koyama & Inutsuka 2002):

5

6

7

Fig. 1 shows the thermal equilibrium state of this net cooling rate in the (n, P) plane. The gas is subject to the cooling (heating) above (below) the thermal equilibrium curve. The bistable fluid consisting of the WNM and CNM can exist in the pressure range P_min < P < P_max, where P_min/k_B= 1596 K cm⁻³ and P_max/k_B= 5012 K cm⁻³. The saturation pressure where the static front is realized is as large as P_sat/k_B= 2823 K cm⁻³ in the plane-parallel geometry.

Figure 1

Thermal equilibrium state of the net cooling rate in (P, n) plane.

Open in new tabDownload slide

In steady state ⁠, equations (1) and (2) can be integrated with respect to x to give

8

and

9

respectively. The energy equation (3) becomes

10

where j and M denote the mass and momentum fluxes, respectively, that are spatially constant and is the specific heat at constant pressure. From equations (4), (8) and (9), density, pressure and velocity can be expressed by using temperature, mass flux and momentum flux as follows:

11

12

13

For convenience, instead of T, we introduce a new variable θ defined by ⁠. Using θ and equations (11)–(13), equation (10) can be rewritten as

14

Shchekinov & Ibáñez (2001) and Inoue et al. (2006) described a detailed numerical method to solve equation (14). Since equation (14) is the second-order differential equation, we need two boundary conditions. At x=−∞, we consider a uniform CNM with thermal equilibrium, ⁠, where physical variables in the CNM are denoted by the subscript ‘c’. The boundary conditions are given by

15

We have two free parameters: P_c and j. Given P_c, the density and temperature can be obtained from equation (4) and the equilibrium condition ⁠. The momentum flux is determined by P_c and j from equations (8) and (9).

The properties of differential equation (14) can be well understood as trajectories in the phase diagram (θ, dθ/dx). From equation (14), one can see that the spatially uniform state with thermal equilibrium (⁠⁠) corresponds to a stationary point because dθ/dx= d²θ/dx²= 0. The topological property of the CNM at x=−∞ corresponds to saddle point (Elphick et al. 1992; Ferrara & Shchekinov 1993). Given P_c and j, we numerically integrate equation (14) from x=−∞ along one of the eigenvectors obtained from the Jacobi matrix of equation (14) around the stationary point of the CNM.

3 ZEL’DOVICH & PIKEL’NER SOLUTIONS

First, we review the solutions connecting two thermal equilibrium states (CNM and WNM). We call the solutions derived by ZP69 the ZP solutions. The ZP solutions exist only in the pressure range P_min < P_c < P_max (see Fig. 1). Given P_c, Shchekinov & Ibáñez (2001) determine j=j_ZP(P_c) as the eigenvalue problem so that the boundary condition dθ/dx= 0 is satisfied in the WNM at x=∞.

For the case of a static solution, i.e. j=v= 0, equation (14) can be rewritten as

16

(ZP69), where the subscript ‘w’ denotes the physical quantities in the WNM. Equation (16) shows the balance between cooling and heating inside the transition layer. The pressure satisfying equation (16) is called saturation pressure.

In steady solutions with non-zero mass flux j≠ 0, integrating equation (10) from x=−∞ to +∞, one obtains

17

and we neglect the second term of equation (10) since the flow is subsonic, and hence almost isobaric. In the saturation pressure, q= 0 since j= 0. For P > P_sat, q is positive because cooling dominates heating inside the front. Thus, j < 0 from equation (17), i.e. the front corresponds to condensation. In contrast, for P < P_sat, the front describes the evaporation. In the ZP solutions, the mass flux j is very small. Thus, the energy is mainly transferred by heat conduction. From equation (10), one can see that the thickness of the front is characterized by the Field length (Belgelman & McKee 1990) defined by

18

The typical values of l_F is as large as 10⁻³ pc in the CNM, and 10⁻¹ pc in the WNM. From equation (17), the flow speed across the transition layer can be evaluated by |v| ≃l_F/t_c, where t_c is the cooling time-scale (Ferrara & Shchekinov 1993). Since this value is very small in the actual ISM, ZP69 suggested that the transition front is almost static.

Fig. 2(a) represents properties of the ZP solutions for P_c/k_B= 3000 K cm⁻³. The corresponding mass flux is as large as j_ZP=−3.2 × 10⁻³m_H km s⁻¹ (v_w= 4.7 × 10⁻³ km s⁻¹). The first row of Fig. 2 shows trajectories of solutions with various boundary conditions in the phase diagram (θ, dθ/dx). There are three stationary points, which are indicated by filled circles. The CNM and WNM are denoted by ‘C’ and ‘W’, respectively. As mentioned in Section 2, the CNM and WNM correspond to saddle points. On the other hand, the stationary point at the intermediate θ corresponds to the unstable phase denoted by ‘U’ in Fig. 2(a), and it is the spiral point (Elphick et al. 1992; Ferrara & Shchekinov 1993). The solutions starting from the CNM are denoted by the thick solid line. The ZP solutions connect the two saddle points ‘C’ and ‘W’. The second row of Fig. 2(a) represents the temperature profile. The trajectory of the ZP solutions in the (n, P) plane is shown in the third row of Fig. 2(a).

Figure 2

Properties of solutions for (a) j/j_ZP= 1, (b) j/j_ZP= 3 and (c) j/j_ZP= 25. The first row indicates the trajectories of solutions with various boundary conditions. The obtained solutions satisfying equation (15) are shown by thick lines. The labels ‘C’, ‘U’ and ‘W’ denote the CNM, unstable phase and WNM, respectively. The second row shows the temperature profile in each mass flux. The third row represents trajectory of the obtained solution in the (n, P) plane.

Open in new tabDownload slide

4 DYNAMICAL TRANSITION LAYERS

4.1 Transition layers for P_min < P_c < P_max

In the ZP solutions, the mass flux is determined so that the CNM is connected with the WNM. Here, we consider the condensation front (j < 0). What happens if the mass flux is larger than |j_ZP|?

4.1.1 Classification of solutions

It is found that solutions for larger mass flux can be divided into four types of solutions: C-U-o, C-U, C-c and C-C. We describe the property of each solution below.

C-U-o solutions. For the case with j= 3j_ZP, the solution does not reach the saddle point of the WNM but approaches the spiral point of the unstable phase (see the first row of Fig. 2b). From the second row of Fig. 2(b), one can see the oscillation of the temperature profile that corresponds to the spiral motion in the phase diagram. The third row shows the trajectory of the solution in the (n, P) plane. The solution truncates at the heating-dominated region before arriving at the WNM. The truncation point corresponds to the peak of the temperature in the second row of Fig. 2(b). Since the mass flux is still low, the pressure is almost constant. This type of solutions has already been found in Elphick et al. (1992) under the isobaric approximation and a more simplified cooling rate. Inoue et al. (2006) also found them as the finite extent solutions, and they truncated the solution at the first peak of the temperature profile near the transition front. We call this type of solutions the ‘C-U-o’ solutions, where ‘C-U’ means solutions connecting CNM and unstable phase, and ‘o’ means the oscillation.

C-U solutions. For larger mass flux j= 25j_ZP, the topological properties of the stationary point of the unstable phase changes from spiral to node (see the first row of Fig. 2c).1 Thus, from the second row of Fig. 2(c), there is no oscillation in the temperature profile, and the CNM is monotonically connected with the unstable phase. Fig. 3 shows solutions for much larger mass flux, j/j_ZP= 25, 1000, 2000 and 3000. One can see that as the mass flux increases, the density (pressure) of the unstable phase at x=∞ increases (decreases) along the equilibrium curve. Because of the large mass flux, the pressure increases significantly towards the CNM. This type of solutions with lower mass flux has already been found in Elphick et al. (1992) under isobaric approximation. We call this type of solutions the ‘C-U’ solutions.

Figure 3

Trajectories of solutions for j/j_ZP= 1 (solid line), j/j_ZP= 1000 (dashed line), j/j_ZP= 2000 (dotted line) and j/j_ZP= 3000 (dot–dashed line). The CNM pressure is set to P_c= 3000k_B K cm⁻³. The filled circle shows the state at x=∞.

Open in new tabDownload slide

C-c solutions. Fig. 4 shows the solutions for j/j_ZP= 300, 500, 1000 and 1500. Here, we set P_c/k_B= 4000 K cm⁻³. The corresponding mass flux is j_ZP=−2.15 × 10⁻²m_H km s⁻¹. The solution for j= 300j_ZP belongs to the C-U solution. For larger mass flux j= 500j_ZP, the solution truncates at a certain point before arriving at the equilibrium state. This critical point corresponds to above which solutions do not exist (see equations 11–13). The physical variables at the critical point is denoted by the subscript ‘crit’. From equations (8) and (9), the condition T=T_crit is equivalent to

19

Note that equation (19) is slightly smaller than the adiabatic sound speed ⁠. These critical velocity and temperature are also seen in the steady-state structure of the shock front with heat conduction without physical viscosity2 (Zel’dovich & Raizer 1967). At the critical point, one can see that equation (14) becomes infinity. Thus, in order to obtain a finite solution, dθ/dx must vanish at the critical point. From equations (13) and (14), the pressure gradient becomes finite as follows:

20

The pressure at the critical point can be derived from the momentum flux conservation. At the CNM (x=−∞), since the velocity is still subsonic because of its high density, the momentum flux becomes

21

On the other hand, from equation (13), the pressure at the critical point is given by P_crit=M/2. Therefore, the pressure at the critical point can be expressed by P_c as follows:

22

The density at the critical point is given by

23

One can see that P_crit for j= 500j_ZP and 1000j_ZP is roughly equal to P_c/2 = 2000k_B K cm⁻³, and its density increases for larger mass flux as seen in equation (23). Around the critical point, since dT/dx= dθ/dx/κ= 0, the trajectories in the (n, P) plane approach constant temperature lines of T=T_crit in Fig. 4.

Figure 4

Trajectories of solutions for j/j_ZP= 300 (solid line), j/j_ZP= 500 (dashed line), j/j_ZP= 1000 (dotted line) and j/j_ZP= 1500 (dot–dashed line). The CNM pressure is set to P_c= 4000k_B K cm⁻³. The filled circle shows the state at x=∞.

Open in new tabDownload slide

We investigate the topological properties at the critical point in the phase diagram. The Jacobi matrix of equation (14) around the critical point becomes infinite. Thus, there is no solution passing the critical point smoothly. In order to obtain infinite extent solutions, shock front is required. Solutions with shock front are described in Appendix A. We call this type of solutions the ‘C-c’ solutions, where ‘C-c’ means solutions connecting the CNM and the critical point.

C-C solutions. For much larger mass flux j= 1500j_ZP, since the solution passes the equilibrium curve before arriving at the critical point, the solution becomes the infinite extent solution connecting the CNM and CNM (see Fig. 4). We call this type of solutions the C-C solutions, where ‘C-C’ means solutions connecting the CNM and CNM.

Fig. 5 shows classification of solutions in the parameter space (P_c, j) for P_min < P_c < P_max. The white dashed line corresponds to j=j_ZP(P_c) on which the ZP solutions exist. The black area where |j| < |j_ZP| is not focused in this paper. From Fig. 5, one can see that the C-c solutions exist only for P_c > 2P_min. This is because if P_c < 2P_min, the solutions meet the equilibrium state before arriving at the critical point (see equation 22).

Figure 5

Classification diagram of obtained solutions in parameter space (P_c, j). The white dashed line corresponds to j_ZP. The solid line indicates the critical mass flux j_crit(P_c). The black area where |j| < |j_ZP| is not focused in this paper.

Open in new tabDownload slide

4.1.2 Structure of transition layers

We investigate how the structure of transition layers depends on the mass flux. As mentioned in Section 3, in the ZP solutions, the thickness of the front is characterized by the Field length because the mass flux is very small. As |j| increases, the effect of advection (the first term of equation 10) becomes important compared with the heat conduction (the third term of equation 10) in the energy transfer. We can define a critical mass flux j_crit where the effect of advection becomes comparable to the heat conduction as follows:

24

where the subscript ‘m’ denotes a reference state in the solutions, l_{c, m} is the cooling length defined by

25

and is the sound speed. As the reference state, we consider the state where has the maximum value for T > T_c under the constant pressure P=P_c. This state roughly gives minimum Field length and cooling length in the solution. Thus, the critical mass flux is a function of P_c. The actual value of j_crit is shown by the solid line in Fig. 7. In most ranges of P_c, j_crit lies in the region occupied by the C-U solutions.

Figure 7

The same as Fig. 5 but including the pressure range P_c > P_max. The dashed line corresponds to T_crit= 10⁴ K.

Open in new tabDownload slide

For |j| > |j_crit|, advection dominates heat conduction during energy transfer. In this case, the thickness of the front is characterized by the cooling length instead of the Field length. In typical ISM, l_F/l_c is roughly as small as 10⁻². Therefore, as the mass flux increases, the thickness of the front increases. Fig. 6 shows the temperature profiles for j=j_ZP, j_crit, 5j_crit and 10j_crit. The CNM pressure is set to P_c/k_B= 3000 K cm⁻³. One can see that the thickness of the layer increases with the mass flux.

Figure 6

Temperature profiles for j=j_ZP (solid line), j=j_crit (dashed line), j= 5j_crit (dotted line) and j= 10j_crit (dot–dashed line). The CNM pressure is set to P_c/k_B= 300 K cm⁻³.

Open in new tabDownload slide

4.2 Transition layers for P_c > P_max

In the previous section, we found solutions only for P_min < P_c < P_max. However, in actual ISM, the CNM pressure can be greater than P_max, for example, in shock-compressed regions. Our solutions can be directly extended for P_c > P_max.

Fig. 7 is the same as Fig. 5, but it includes the pressure range P_c > P_max. Fig. 7 shows that most of the regions are covered by the C-c solutions for P_c > P_max. The area of the C-U solutions spreads to the pressure range P_max < P_c < 2P_max. This is because solutions in this pressure range can connect with thermal equilibrium states before reaching at the critical point at a certain range of j. Fig. 8 shows trajectories of the solutions for j/j_crit= 1, 2, 5 and 20. The CNM pressure is set to 2.65 × 10⁴ K cm⁻³. One can see that the solution terminates at the critical point where pressure is equal to P_c/2 and where density increases with the mass flux.

Figure 8

Trajectories of solutions for j/j_crit= 1, 2, 5 and 20. The dashed line corresponds to P=P_crit/2.

Open in new tabDownload slide

5 DISCUSSION

5.1 Implications

Koyama & Inutsuka (2006) have investigated multidimensional non-linear development of the TI without external forcing. They found that the turbulence is self-sustained by the TI. The velocity dispersion for sufficiently large simulation box is as large as ∼0.2–0.4 km s⁻¹. They suggested that this velocity is much larger than the prediction in ZP69. Our solutions exist in such a large mass flux. From their simulation, the mass flux becomes |j| ∼ 0.2 m_H km ⁻¹(n/0.5 cm⁻³)(v/0.4 km s⁻¹), where we assume that the accretion velocity to the CNM is comparable to this velocity dispersion. From Fig. 5, the corresponding solutions lie in the C-U or C-U-o solutions.

Koyama & Inutsuka (2002) have investigated the TI in the shock-compressed region. They found that the velocity dispersion is as large as several km s⁻¹. The corresponding mass flux is |j| ∼ 20 m_H km s⁻¹(n/10 cm⁻³)(v/2 km s⁻¹), which is larger than the results in Koyama & Inutsuka (2006). The pressure of the shock-compressed region is larger than P_max, indicting that the two stable phases cannot coexist. In shock-compressed regions, the CNM clouds are embedded not by the stable WNM but by the unstable gas that is supplied continuously by shock compression. The C-c solutions in this pressure range (P > P_max) may describe dynamical condensation from the unstable gas to the CNM in shock compression regions.

5.2 Curvature effect

Our solutions are derived by assuming the plane-parallel geometry. Nagashima et al. (2005) investigated the curvature effect of the transition front by assuming the quasi-steady approximation. They consider the cylindrical and spherical CNM cloud at the centre. In contrast to the plane-parallel solutions, the solutions have three parameters, the pressure of the CNM, the velocity of the front with respect to the centre and the cloud radius. The second parameter corresponds to the mass flux in our solutions. They solved the equations as the eigenvalue and boundary value problems so that the solutions connect the CNM and WNM. Thus, the velocity of the front, v_f(P_c, R_c), is a function of P_c and cloud radius R_c. The curvature effect enhances the heat flux by thermal conduction from the CNM to WNM. Thus, smaller cloud that has larger curvature tends to evaporate quickly. On the other hand, cold gas that has a negative curvature tends to gain mass through condensation. Thus, the evolution of CNM cloud strongly depends on its shape.

By using the quasi-steady approximation, we can easily take the curvature effect into account in our solutions. It is expected that solutions with larger evaporation rate connect the unstable phase and the WNM, while solutions with larger condensation rate connect the CNM and the unstable phase.

6 SUMMARY

In this paper, we have investigated steady condensation solutions of phase transition layers in larger parameter space (P_c, j) than that in previous works. We summarize our results as follows.

• In the pressure range where the three thermal equilibrium phases can coexist under a constant pressure (P_min < P_c < P_max), we find solutions that connect the CNM and the unstable phase. The solutions can be classified into four types of solutions, C-U-o, C-U, C-c and C-C, in order of the mass flux. The C-U-o and C-U solutions connect the CNM and the unstable phase with and without the oscillation of physical quantities in the low-density side. The C-c solution is truncated at the maximum temperature above which there are no steady solutions. Combination of a C-c solution and a shock front provides an infinite extent solution. The C-C solutions connect the CNM and CNM.

• We also find steady solutions in the pressure range (P_c > P_max) where only CNM can be in equilibrium. This pressure range is realized in the shock-compressed region. In this parameter space (P_c, j), most of the solutions are of C-c type.

• We derive a critical mass flux j_crit above which the effect of advection becomes important compared with the heat conduction during energy transfer. For large mass flux, |j| > |j_crit|, the thickness of the front is characterized by the cooling length instead of the Field length, so that the thickness becomes larger.

We thank the referee for valuable comments. We thank Drs Tsuyoshi Inoue and Jennifer M. Stone for valuable discussions. This work was supported by Grants-in-Aid for Scientific Research from the MEXT of Japan (KI: 22864006; SI: 18540238 and 16077202).

Footnotes

REFERENCES

Appendix

APPENDIX A: SOLUTIONS WITH SHOCK FRONT

In this appendix, we present steady solutions with shock front. Here, we assume that pre-shock gas is in thermal equilibrium state. The physical variables in the pre-shock gas are denoted by using the subscript ‘E’. The position of the shock front x=x_sh is determined so that the Rankine–Hugoniot relation is satisfied between the solutions and the thermal equilibrium state by the following method. Given x=x_sh, the specific energy flux is given by

(A1)

where the subscript ‘sh’ denotes physical variables at x=x_sh and is the effective sound speed. The simple relation among v_sh, v_E and c_* is given by

(A2)

From equation (A2), one can get v_E, and the Mach number is obtained. From the Rankine–Hugoniot relations, the physical variables in the pre-shock gas (ρ_E and T_E) are obtained. In general, the pre-shock gas is not in thermal equilibrium state. The position of the shock front x_sh is determined iteratively until is satisfied by using the bisection method.

As examples, we consider two cases of j= 2 and 24j_crit for P_c/k_B= 3.114 × 10⁴K cm⁻³. For the j= 2j_c case, the corresponding physical quantities of the pre-shock gas are n_E= 0.57 cm⁻³ and ⁠. For j= 24j_c, they are n_E= 75 cm⁻³ and ⁠. Thus, the pre-shock gases belong to WNM for j= 2j_c and to CNM for j= 24j_crit. The upper panels of Figs A1(a) and (b) show the density, temperature and pressure distributions for j= 2 and 24j_crit, respectively. The lower panels of Fig. A1 show the trajectories of the obtained solutions in the (n, P) plane.

Figure A1

Upper panels: density (solid line), temperature (dashed line) and pressure (dotted line) profiles of steady-state solutions at which P_c= 3.113 × 10⁴ K cm⁻³ for (a) j= 2j_crit and (b) j= 24j_crit. Lower panels: trajectories of steady-state solutions in the (n, P) plane. The thick grey line in each panel shows the results of one-dimensional hydrodynamical simulation.

Open in new tabDownload slide

In addition, we perform one-dimensional hydrodynamical simulation in order to confirm whether the steady-state solutions with shock front are realized. We consider the situation where the gas flows collide at x= 0 with |v_E−v_c| from x=±∞. Two shock waves propagate from x= 0 outwards. Shock-heated gas quickly cools until the gas reaches the thermal equilibrium state. Then the cold layer forms and its thickness grows by gas accretion. The thick grey line in each panel of Fig. A1 indicates the results of one-dimensional simulations in the growing phase of the cold layer. One can see that the results of one-dimensional simulations are well described by the steady-state solutions.