© ESO, 2014

1. Introduction

Understanding the star formation process remains one of the main areas of research of modern astrophysics. It is directly linked to the way interstellar gas is organized and how dense and cold structures form. Turbulence plays a fundamental role here. It creates local increases in density that are amplified by gravity and that, in some cases, might lead to gravitational contraction. At the same time, the kinetic energy of turbulent motions needs to be dissipated in order for the gravitational collapse to happen. The overall star formation process is therefore directly linked to the way turbulence shapes the density structure of the ISM and how kinetic energy is transferred from large to small scales.

Where star formation occurs, molecular clouds exhibit suprathermal linewidths and strong contrast of density revealed by a lognormal histogram of column density (Vázquez-Semadeni & García 2001; Goodman et al. 2009; Federrath et al. 2010; Brunt 2010) and by power spectra of column density shallower (–2.5, Bensch et al. 2001) than predicted for subsonic or transonic turbulence (–3.66) and than what is observed in the diffuse H i (Miville-Deschênes et al. 2003a). For these reasons they are often modeled as isothermal and supersonic turbulent flows. On the other hand, because turbulence should dissipate in one dynamical time if it is not maintained on a large scale it is thought that star formation might occur within a dynamical (or free fall) time (Elmegreen 2000; Hartmann et al. 2001). In this scenario the gravitational collapse occurs during the cloud formation. Therefore, the initial conditions of the star formation process depend strongly on the structure of the cold clouds in the atomic phase (H i).

Most of the mass of the Galactic interstellar medium (ISM) is H i and there is clear observational evidence that it is thermally bistable (Dickey et al. 2003); it consists of two thermally stable states at significantly different temperatures (the cold neutral medium – CNM – and the warm neutral medium – WNM). This is a natural result of the density and temperature dependence of the cooling processes in play in the diffuse ISM. Like molecular clouds, the structure of the H i is highly complex with strong density contrasts (the CNM fills only about 1% of the volume but has a density 100 times higher than the WNM), but here the combination of turbulence and thermal instability (TI; Field 1965, 1969) is the main driver of the structure. It is important to note that the nature of the density fluctuations in an isothermal and supersonic flow is conceptually different than in a thermally bistable fluid. An isothermal turbulent flow with supersonic velocity fluctuations will naturally develop strong density fluctuations (, Vázquez-Semadeni 2012). In such a flow the density fluctuations are transient, associated with the compression and depression of the turbulent motions. In a two-phase medium, once a high density structure is created by a local compression it remains at high density even when the gas re-expand because it is in pressure equilibrium with the inter-cloud medium. A two-phase medium will therefore produce long-lived high density structures contrary to supersonic flows in which the fluctuations appear and disappear on dynamical timescales. Therefore the thermally bistable H i offers a natural way of producing highly structured initial conditions for the star formation in molecular clouds.

Because of the complex dynamical processes involved, dedicated numerical simulations are useful to develop insights on the exact scenario that leads to the formation of cold structures in the H i. Several studies based on numerical simulations have explored the thermally bistable physics of the H i since the early one-dimensional work of Hennebelle & Pérault (1999) who showed the effect of compression in a WNM converging flow on the triggering of the phase transition and the production of CNM structures. Similarly Koyama & Inutsuka (2000, 2002) studied the impact of the propagation of a shock wave into WNM. They showed how CNM structures form in the post-shock region through the thermal instability which operates faster than the duration of the shock (galactic spiral shock or supernova). In both cases the resulting density contrast between the CNM and WNM is ~100. To obtain a similar contrast in isothermal supersonic turbulence, Mach numbers of the order of 10 are required. These studies have been extended to two and three dimensions simulations of thermally bistable flows with turbulence induced by a converging flow (Audit & Hennebelle 2005; Hennebelle & Audit 2007; Hennebelle et al. 2007; Audit & Hennebelle 2010; Heitsch et al. 2005, 2006), by the magneto-rotational instability (Piontek & Ostriker 2004, 2005, 2007; Piontek et al. 2009) or by a driving in Fourier space (Gazol et al. 2005; Gazol & Kim 2010; Seifried et al. 2011).

It was established that compression of the WNM in transonic converging flows triggers the thermal instability; the increase in density can be such that the gas cools rapidly to temperatures of a few tens of degrees. The cold gas then fragments into clumps, in particular through the Kelvin-Helmholtz instability (Heitsch et al. 2006). The thermal pressure of the cold gas is in approximate equilibrium with the sum of the ram and thermal pressure of the surrounding warm gas. Once pushed at high enough densities, the cold structures can have densities and temperatures close to values typical of molecular clouds (Vázquez-Semadeni et al. 2006; Banerjee et al. 2009).

In general, previous works showed that the efficiency of the formation of the CNM depends directly on the nature and amplitude of the injected energy. One question we wanted to explore in this study is related to the driving of the compression. Most studies have explored the case of converging flows where the increase of the WNM density needed to trigger the thermal instability is obtained by a sustained ram pressure. The origin of such large scale convergent flows is unclear so here we want to see if compression in turbulent motions alone, driven on large scales, could lead to the production of cold gas. More precisely, we wanted to see if turbulent motions, at the level they are observed in the WNM, can induce the phase transition and if not, what are the required physical conditions of the WNM to trigger the formation of cold structures. To do so, like Gazol et al. (2005), Gazol & Kim (2010), Seifried et al. (2011), we used turbulent forcing in Fourier space. In addition, similarly to the study of Federrath et al. (2010) based on isothermal simulations of molecular clouds, we wanted to study the impact of compressive and solenoidal modes in the injected turbulent field.

In this paper we present a parameter study of WNM initial density and turbulent forcing properties on the formation of the CNM in forced turbulence and thermally bistable flows. In a similar study, Seifried et al. (2011) chose initial densities (1 to 3 cm^-3) that are in the thermally unstable regime, to expect the production of a two-phase medium. Our intent is to expand the parameter space to identify the dynamical conditions that lead to the formation of the cold phase. In addition, one fundamental aspect of the current study is that we guide our parameter search based on all available observational constraints of the local ISM. We present 90 simulations performed on a 128³ grid that allowed us to identify more precisely the physical conditions that favor the formation of cold structures while reproducing in detail the observed properties of the H i in the local ISM. For two specific cases that match the observations we performed two 1024³ simulations in order to study in greater details the properties of the H i: temperature, pressure and density histograms, power spectra, and mass-size and velocity dispersion-size relations of the cold structures.

The paper is organized as follow. HERACLES, the numerical code used in this study is briefly presented in Sect. 2. The observational constrains used to lead the parameter study are detailed in Sect. 3. The methodology is presented in Sect. 4. The results of the parametric study and of the high resolution simulations are presented and discussed in Sects. 5 and 6 respectively. The main conclusions of the study are given in Sect. 7.

2. Numerical methods

Despite their great importance in the ISM, we won’t consider magnetic field or gravity in this study to focus on the impact of vortical and compressive turbulence motions on the formation of CNM. All the following simulations are thus hydrodynamical. We performed our simulations with HERACLES (González et al. 2007). The numerical methods are similar to those described in Audit & Hennebelle (2005).

2.1. Euler equations and Godunov scheme

The Euler equations for a radiatively cooling gas are the classical equations of hydrodynamics used for flows at high Reynolds number, which is the case in the ISM: where ρ is the mass density, u the velocity, P the pressure and E the total energy. The atomic gas is considered perfect and ideal with an adiabatic index γ = 5 / 3 and a mean atomic weight μ = 1.4 m_H (m_H is the mass of the proton). The weight μ takes the cosmic abundance of other elements into account: helium and metals, and so the mass density is ρ = μ × n.

The net cooling function is defined as ℒ = n²Λ − nΓ (expressed in erg cm^-3 s^-1), where Λ and Γ are respectively the cooling and heating terms based on the prescription of Wolfire et al. (1995, 2003). They include the most relevant processes in the diffuse medium for our range of temperature [10 K, 10⁴ K] and are plotted and described in Fig. 1 (for more details, see Audit & Hennebelle 2005). The equations solver in HERACLES is a second order Godunov scheme (see Audit & Hennebelle 2010) for the conservative part. It is a MHD solver of type MUSCL-Handcook (Londrillo & Del Zanna 2000; Fromang et al. 2006). The cooling is applied after the solver (see Audit & Hennebelle 2005).

2.2. Turbulence forcing

An external field f, applied in Fourier space, is added to generate large scale turbulent motions. It has the dimension of an acceleration. A study of a driving consisting of a true force compared to a driving consisting of an acceleration has been performed by Wagner et al. (2012). They demonstrated that modeling the force density instead of the acceleration at low wavenumbers directly couples the external energy injection to all scales of the system, as the density field has Fourier modes over the whole range of wavenumbers.

The magnitude of the applied acceleration is where v_S represents the amplitude given to the field¹ and L_S = L_box/k₀ is the characteristic integration scale with k₀ = 2 the characteristic wave number of the stirring. The auto-correlation time of the process is given by τ_S = (L_S/F_S)^1/2 = L_S/v_S. The stirring is applied on large scales, and thus to the small wave numbers: 0 < | k | < 2k₀.

The modes of the turbulent forcing are computed with the stochastic process of Ornstein-Uhlenbeck (Eswaran & Pope 1988; Schmidt et al. 2006). The method is similar to Schmidt et al. (2009) and Federrath et al. (2010). A Helmholtz decomposition is applied to the random field to project it on its compressive and solenoidal components. The projection operator is defined as (4)Choosing a spectral weight ζ = 1 creates a divergence free field, i.e. purely solenoidal. Taking ζ < 1 generates dilatational components which compress or rarify the gas. Finally, ζ = 0 creates a rotation-free field, i.e. purely compressive. In a natural state for which the energy is naturally distributed between the compressive and solenoidal modes (2:1 mixing), ζ takes the value of 0.5. In this case, the Helmholtz operator is proportional to the identity operator and the projection does not alter the energy repartition between the different modes.

3. Observational constraints

The simulations presented in this study are dedicated to the formation of CNM structures out of pure WNM gas. In this section we present global properties of the Galactic H i as deduced from the observations. These properties guided us in the choice of the initial conditions of the simulations and on the qualification of the results obtained.

3.1. Pressure and thermal bistability of the H I

Because of the density and temperature dependence of the heating and cooling processes in the diffuse ISM (see Fig. 1) it has been shown early on by Field (1965), Field et al. (1969) and updated by Wolfire et al. (1995, 2003) that there is a range of pressure where the neutral atomic gas can be in two thermally stable phases in pressure equilibrium. This is often displayed in diagrams like the ones shown in Fig. 2 where the solid curves in the two panels indicate the locus of thermal equilibrium in the P − n and T − n parameter spaces.

The thermal stability is not preserved everywhere along these curve but only where dP/dn > 0. The two dash-dotted vertical lines mark the boundaries in density where the H i is thermally unstable. Any perturbation of the gas sitting on the equilibrium curve in the thermally unstable range would move toward the cold or warm branches. From Fig. 2 one can appreciate the steepness of the equilibrium curve in T − n space, highlighting the great difference in temperature and density of the two stable phases of the H i.

It is important to point out that the bistability of the H i is only present in a given range of pressure, determined by the maximum density of the warm branch (n ~ 0.8 cm^-3) and the minimum density of the cold branch (n ~ 7 cm^-3). For P < 1000 Kcm^-3 there is no CNM and for P> 6000 Kcm^-3 all the H i is cold.

On the top panel of Fig. 2 the dashed curve gives the average pressure observed in the local ISM. According to Jenkins & Tripp (2011) the pressure of the cold phase follows a log-normal distribution with an average of log ₁₀(P/k) equals to 3.58 dex and minimum and maximum values (at the 99% level) equal to 3.04 and 4.11 dex (dotted lines). The minimum CNM pressure reported by Jenkins & Tripp (2011) corresponds exactly to the minimum CNM pressure (3.06 dex) allowed by the thermal instability (see Fig. 2) as predicted by the modeling of the heating and cooling processes of Wolfire et al. (1995, 2003). The maximum pressure reported by Jenkins & Tripp (2011) is significantly higher than the maximum pressure allowed for the thermally stable WNM (3.74 dex). At such high pressure the gas can only be thermally stable in the cold phase.

3.2. WNM density and temperature

Here we want to estimate what is the true density (as opposed to space-averaged) and temperature of the WNM in the local ISM. Typical numbers cited are n ~ 0.2−0.5 cm^-3 and T ~ 6000 − 10⁴ K (Ferrière 2001), but, as emphasized by Kalberla & Dedes (2008), the density of the WNM varies greatly with galacto-centric radius and galactic height z. It is also subject to local variations due to star formation activity (e.g. supernovae, Joung et al. 2009) and between the arm and inter-arm regions.

According to Dickey & Lockman (1990, see also Celnik et al. 1979; Malhotra 1995; Kalberla et al. 2007), the distribution of the H i density with z close to the Sun can be decomposed in a narrow Gaussian of HWHM = 106 pc, a broader Gaussian of HWHM = 265 pc and an exponential with a scale height of 403 pc. Generally the narrow Gaussian is attributed to the CNM and the sum of the two others to the WNM. In this scheme, and still following Dickey & Lockman (1990), the space-average density of the CNM and WNM at z = 0 are 0.395 cm^-3 and 0.172 cm^-3 respectively. These values are weighted by the volume filling factor φ and thus correspond to nφ, where n is the true volume density. The filling factor of the WNM is difficult to estimate; according to Kalberla & Kerp (2009)φ ~ 40% which implies a true density of n_WNM = 0.43 cm^-3.

This value is close to the one that can be derived based on the knowledge of the heating and cooling processes. For the average pressure of P = 3980 Kcm^-3 found by Jenkins & Tripp (2011), using the heating and cooling curves of Wolfire et al. (1995, 2003) and assuming a spatially uniform radiation field with the spectrum and intensity of the Habing field (G₀/ 1.7 where G₀ is the Draine flux – Habing 1968; Draine 1978), the equilibrium density and temperature of the WNM are n_wnm = 0.5 cm^-3 and T_wnm = 7960 K.

Figure 3 shows how n, T and P vary with z assuming a scale height of HWHM = 265 pc for the WNM disk at R = R_⊙ (Dickey & Lockman 1990; Malhotra 1995; Kalberla et al. 2007). Assuming a true volume density n_WNM = 0.5 cm^-3 at z = 0, the corresponding temperature profile T(z) of the WNM can be estimated using the thermal equilibrium curve based on the heating and cooling processes described in Wolfire et al. (1995, 2003). It shows that the variations with z of the WNM properties due to the hydrostatic equilibrium become important only on scales larger than 100 pc.

3.3. WNM turbulent velocity dispersion and Mach number

Apart from its density and temperature, the other relevant parameter for this study is the level of turbulence in the WNM, or equivalently its Mach number. The direct estimate of the Mach number of the WNM is difficult due to its temperature that cannot be easily measured. The traditional method to estimate the H i temperature is based on a comparison of 21 cm absorption and emission measurements. Due to its low opacity, the warm phase produces almost undetectable absorption once observed in front of a strong radio-continuum source. Only a few detections were reported so far (e.g. Carilli et al. 1998; Kanekar et al. 2003; Begum et al. 2010), with deduced temperatures closer to the thermally unstable regime (T ~ 500−6000 K) than to what is expected for the WNM (T ~ 8000 K).

Nevertheless there is indirect evidence that the warm ISM is subsonic. Recent findings on the warm ionized medium (WIM) by Gaensler et al. (2011) confirmed the work of Armstrong et al. (1995) and Redfield & Linsky (2004) that showed that the WIM has the properties of a low Mach number and incompressible turbulent flow. Regarding the WNM, an indication of its low Mach number comes from the smoothness of the 21 cm profiles. Haud & Kalberla (2007) found that the simplest 21 cm spectra of the LAB data (Kalberla et al. 2005) at high Galactic latitudes are well modeled by a single Gaussian of σ_tot = 10.2 ± 0.3 km s^-1. Here we assume this value to be representative of the WNM total velocity dispersion in the solar neighborhood.

The WNM velocity dispersion is the quadratic sum of the thermal (σ_therm) and turbulent (σ_turb) contribution to the gas motions, integrated along the line of sight. Considering the specific line of sight at the Galactic pole, the thermal contribution to the linewidth of the WNM is an integral over Galactic height z: (5)where n(z) is assumed to be a Gaussian with HWHM = 265 pc (Dickey & Lockman 1990). Using the corresponding T(z) shown in Fig. 3, integrating Eq. (5) gives σ_therm = 8.3 km s^-1, a value essentially identical to the thermal width of an isothermal 8000 K gas (σ_therm = 8.1 km s^-1) which illustrates that the WNM temperature does not vary significantly with Galactic height. Assuming σ_tot = 10.2 km s^-1 (Haud & Kalberla 2007), the contribution of turbulent motions to the observed linewidth is σ_turb = 5.9 km s^-1. That σ_turb < σ_therm is another indication that thermal motions dominates over turbulent motions on any scale in the solar neighborhood.

Due to the nature of the turbulent cascade, the turbulent velocity dispersion increases with scale like σ_turb = σ(1)l^q, where q ~ 1/3 for subsonic turbulence and, following the notation of Wolfire et al. (2003), σ(1) is the velocity dispersion (in km s^-1) on a scale of 1 pc and l is the scale (in pc). Assuming that the turbulent contribution to the observed linewidth is representative of motions on the scale of the H i scale height at the Sun location (l = h_z = 265 pc), one obtains σ(1) = 0.92 km s^-1. This calculation is similar to the one of Wolfire et al. (2003) who estimated σ(1) = 1.4 km s^-1 for the WNM. That we obtain a significantly lower value can be attributed mostly to our having taken the thermal contribution to the linewidth into account.

3.4. CNM mass fraction and volume filling factor

One important quantity we want to reproduce in the simulations is the amount of mass in the cold phase. The CNM mass fraction was estimated to 40% by Heiles & Troland (2003) by comparing 21 cm absorption and emission measurements on 79 random lines-of-sights. This is also compatible with the decomposition of the H i density distribution with z of Dickey & Lockman (1990); their narrow component (h = 106 pc) that can be associated to the CNM contributes to 44% of the total H i surface density of the disk. Because the scale height of the CNM is less than for the WNM, the CNM mass fraction varies with z. The CNM contributes to 60% of the H i mass for z < 200 pc.

Like for the WNM one can compute the expected CNM volume density based on the heating and cooling processes (Wolfire et al. 1995, 2003). For the average pressure of Jenkins & Tripp (2011, P = 3980 Kcm^-3) the equilibrium density and temperature are n_cnm = 111 cm^-3, T_cnm = 35.9 K. Following Dickey & Lockman (1990), the space-average density of the CNM is ~0.4 cm^-3 at z = 0, which implies a volume filling factor of the cold phase lower than 1%.

4. Methodology

Based on the observational constraints described in the previous section, we conclude that the WNM in the solar neighborhood has a mean temperature T ~ 8000 K, a mean density of n ~ 0.2−0.5 cm^-3 and subsonic turbulent motions following σ_turb = 0.92 (l/1pc)^1/3 (km s^-1). The first step of this study is to evaluate what are the plausible dynamical conditions that lead to the formation of CNM gas, of the order of 40% in mass, from WNM gas with such physical properties. We explored those conditions with 90 simulations at low resolution (128³ cells). In a second step, two sets of initial conditions that match the H i observations are used to produce high resolution simulations (1024³) and study in more details the properties of the CNM structures.

In this section we first described the characteristic scales involved in the simulations, then we validate the use of 128³ simulations for the exploration of the parameters space by comparing with higher resolution (256³ and 512³) simulations.

4.1. Characteristic scales

4.1.1. Dynamical time and cooling scale of the WNM

If, during the perturbation of an initial WNM cloud, the cooling processes act slower than the compressive ones, the gas will tend to return to its initial state of temperature and density by decompression. The WNM-CNM transition becomes possible when the radiative time is lower than the dynamical time (Hennebelle & Pérault 1999), allowing the condensed gas to keep its higher density even during the decompression: (6)The dynamical time is defined as the time needed by an atom to cross the cloud (of size d) at the sound speed velocity c_S: (7)The radiative time is the characteristic time needed by the medium to cool. It is the ratio of the energy (U) to the energy transfer (ϵ) applied to the medium: (8)From these equations we can calculate the critical size above which the WNM will be unstable and will condense into the unstable regime or CNM structures. We consider typical conditions of the WNM (T = 8000 K and n = 0.5 cm^-3) and a static medium, assuming that the energy is fully thermal. Using Eq. (6) and the heat capacity per mass unity C_V = k_B/m_H(γ − 1): The cooling time thus becomes: (11)Considering the main cooling processes of the diffuse H i described in Sect. 3.1 to evaluate Λ, and using , one obtains the following WNM cooling scale: (12)Therefore, for the thermal instability to develop, the scale of the dynamically perturbed WNM must be greater than λ_cool ~ 22 pc. Similarly, a perturbation of size lower than λ_cool do not efficiently trigger the phase transition.

4.1.2. Sizes of the simulations

To study the WNM-CNM transition it is important to evaluate the scales of cold structures, in order to adjust the numerical resolution. As illustrated by the 2D simulations of Hennebelle et al. (2007) and the 3D simulations of Audit & Hennebelle (2010), the size of the CNM clumps are distributed in the range 10^-3−0.5 pc. Moreover, the Field length λ_F (Field 1965; Begelman & McKee 1990), which gives the typical scale of the thermal front between the WNM and the CNM, as well as the shock front of CNM structures, are both of the order of 10^-3 pc. Therefore, in principle simulations of the thermally bistable H i would require a resolution of the order of 10^-3 pc. Nevertheless, as suggested by Hennebelle & Audit (2007) and as we will show in the next section, the most important scale to resolve in simulations of thermally bistable and turbulent H i is the “sound crossing scale” (i.e., the sound speed times the cooling time). As described in the previous section, this scale is ~22 pc for typical values of the WNM density, providing a minimum size for the simulation box. We chose to assure the phase transition with a box size equal to 40 pc. In the following we use box sizes of 128³ pixels and 1024³ pixels, corresponding respectively to resolution of 0.3 and 0.04 pc.

4.2. Effect of resolution

In Sect. 5 we present a parametric study based on 128³ pixels simulations. Even though we barely resolve the typical scale of CNM structures with those simulations, we will demonstrate using higher resolution simulations that this set up can be used to explore efficiently the parameter space for the formation of the CNM. Figure 4 shows the evolution over time of the CNM mass fraction, the dispersion of the line of sight turbulent velocity σ_turb and the observational Mach number for two representative simulations ([n₀ = 1 cm^-3, ζ = 0.2, v_S = 12.5 km s^-1] and [n₀ = 2 cm^-3, ζ = 0.5, v_S = 7.5 km s^-1]) done on 128³, 256³ and 512³ grids. We calculate the cold mass fraction f_CNM as the sum of the mass in pixels where T < 200 K, over the total mass in the box. The turbulent velocity dispersion is computed the following way (see Miville-Deschênes & Martin 2007, Eq. (8)): (13)where C(x,y) is the velocity centroid at position (x,y): (14)In numerical simulations, the Mach number is often computed as the ratio of the velocity and the sound speed on each cell: (15)where . Because this parameter cannot be directly deduced from observations we consider an expression for the Mach number calculated from integrated quantities along the line of sight, the turbulent velocity dispersion and the thermal velocity dispersion: (16)where the thermal velocity dispersion is defined as (17)Obviously the total² velocity dispersion is simply (18)All quantities shown in Fig. 4 are averaged over the box at different time steps. For all resolution, the mean values of σ_turb and ℳ_obs are in the same range, and their variations are of the same order. These two quantities are thus statistically well represented even in small resolution simulations.

The evolution with time of f_CNM (Fig. 4) as well as the histograms of temperature and density (Fig. 5) indicate that the fraction of the mass in the cold gas depends on the numerical resolution. Higher resolution runs lead to a slightly higher amount of cold gas. The cold peak of the temperature histogram (Fig. 5) is lower in the 128³ and 256³ cases. The time evolution of f_CNM increases more slowly for the smallest resolution run. Nevertheless, the differences in the mass fraction of gas with T < 200 K after 40 Myr remain small for the different resolutions (see Fig. 4). We verified that these results are independant of the initial conditions.

From this analysis, and like Seifried et al. (2011), we conclude that the resolution does not affect the time evolution of σ_turb and ℳ_obs and has a moderate effect on the formation of CNM. Therefore doing a statistical study on small simulations (128³) is legitimate. Higher resolution simulations are needed to study the CNM properties in details, especially its structure.

4.3. Convergence

In order to evaluate the time needed for the 128³ simulations to reach a stationary state, we evaluated σ_turb, ℳ_obs and f_CNM at each Myr. Most simulations converged in 10 Myr but some needed a much longer time to transit to cold phase and never reached a stationary state for the mass fraction of CNM. However we chose to stop all simulations at 40 Myr because of the typical time between two supernovae explosions estimated to be around 20 Myr (see Ferrière 2001, Eq. (17) at z = 0 in a volume of (40 pc)³).

5. Parametric study

All the simulations presented in this study have periodic boundary conditions. The velocity field is modified at every time step by the turbulent stirring in Fourier space (Sect. 2.2). The initial conditions correspond to a uniform and static WNM at a temperature of T = 8000 K. For a box size of 40 pc centered on z = 0, the profiles shown in Fig. 3 indicate that constant n and T are representative of the WNM.

This section presents the parametric study we did on the initial WNM density, n₀, and on two properties of the stirring velocity field: the spectral weight, ζ, and the large scale velocity, v_S. Overall we performed 90 simulations with 128³ cells, varying n₀ from 0.2 to 10 cm^-3, v_S from 5 to 20 km s^-1 and ζ from 0 (rotation free field) to 1 (divergence free). More specifically we explored the parameter space in two ways. First we looked at 48 simulations with variable n₀ and v_S, and a fixed ζ = 0.5 (Sect. 5.1). Then we investigate the effect of ζ and v_S for a constant initial density of n₀ = 1 cm^-3 (Sect. 5.2). A test of the robustness of the results is presented and discussed in Sects. 5.3 and 5.4, where we compare these driven simulations to simulations with decaying turbulence.

5.1. Effect of the initial density

Table 1 and Fig. 6 present the results for 48 simulations for which the turbulence driving has been kept to a natural state (i.e., the ratio of compressible over solenoidal modes of the velocity field is kept in its initial state) using a fixed ζ = 0.5. For those simulations we only vary v_S and n₀. For each simulation we computed the CNM mass fraction f_CNM, the mean Mach number ℳ_obs (Eq. (16)) and the mean velocity dispersion (Eq. (13)) at each Myr. In addtition we computed the average of those quantities over the time interval [20 Myr, 40 Myr] that corresponds to a period where a stationary state has been reached.

We observe that simulations with n₀ between 0.2 and 1.0 cm^-3 do not transit to cold gas. On the other hand, the efficiency of the CNM formation increases rapidly at higher values of n₀. For n₀ = 0.2 and 1.0 cm^-3, the cold mass fraction is very close to zero while a small increase of the initial density (n₀ = 1.5 cm^-3) leads to f_CNM= 0.3−0.4. Initial densities greater than 3.0 cm^-3 lead to more than half of the mass to transit (see Fig. 6 left part). This behavior can be appreciated by looking at integrated density maps (Fig. 7). The map on the left, that corresponds to n₀ = 1 cm^-3, contains only WNM and is smooth. With the increase of n₀, more and more contrasted structures appear, corresponding to the formation of CNM.

One can notice that n₀ = 1.5 cm^-3, the minimum initial density at which significant CNM is formed, corresponds to a pressure (12 000 K cm^-3) much higher than the maximum pressure allowed for the WNM (see Sect. 3.1). At such high pressure and low density, the gas starts in the thermally unstable regime; it is indeed located above the stable branch of the WNM and in a zone of cooling that allows it to evolve rapidly and efficiently toward the cold phase.

We also note that the final turbulent velocity dispersion, σ_turb, is lower for higher values of n₀ and therefore for higher values of the cold mass fraction f_CNM (see Table 1 and Fig. 6). To first order, this can be understood by the fact that once a CNM structure is formed, it concentrates a fraction of the mass in a small volume. The velocity dispersion of this structure is lower than if the same amount of gas was spread out in the WNM on much larger scales (σ_turb ∝ l^q). Because σ_turb is weighted by density (see Eq. (13)), its average value decreases when the amount of mass located in small scale structures increases, or equivalently when f_CNM increases.

To understand this behavior in more details, we studied the distributions of σ_turb, σ_therm and σ_tot in a simulation at high resolution (1024³ cells, n₀ = 2.0 cm^-3, ζ = 0.5 and v_S = 7.5 km s^-1). All distributions are plotted in Fig. 8, each velocity dispersion was computed using all pixels (black line) and also by selecting pixels corresponding to each thermal phase (CNM: T < 200 K, thermally unstable: 200 <T< 5000 K, WNM: T > 5000 K). As expected σ_therm depends on the temperature of the gas. We note that it is also the case for the distributions of σ_turb, the one computed in the CNM being especially different from the other two. The distribution of σ_turb(CNM) peaks at very small values between 0.1 and 0.2 km s^-1 with a tail reaching 5 km s^-1. We recall that these velocity dispersions are computed by integrating along the line of sight (Eqs. (13), (17) and (18)). To better understand the shape of these distributions, we looked at two distinct lines of sight for which σ_turb(CNM) is equal to 0.14 and 4.7 km s^-1 respectively. The density, temperature and velocity (z-component) profiles are plotted in Fig. 9. The line of sight with a low σ_turb(CNM) (Fig. 9 – top) is composed of four CNM structures with temperatures between 100 and 500 K. These clumps have homogeneous internal velocities and small relative velocities (about 3 km s^-1). On the other hand, for the line of sight with a higher σ_turb(CNM) (Fig. 9 – bottom), the internal velocities of the CNM structures are also homogeneous, but their relative velocities are high, about 15 km s^-1. This analysis indicates that individual clumps have an internal velocity dispersion of about 0.2 km s^-1 (the peak at low values in the histogram of σ_turb(CNM)) and that high values of σ_turb(CNM) are the result of the relative velocities between clumps. Moreover, the similarities of the histograms of σ_tot and σ_therm in Fig. 8 indicates that the velocity dispersion is dominated by thermal motions suggesting a subsonic Mach number. Like Heitsch et al. (2005, 2006) we conclude that the cold gas structures have subsonic internal motions and that the supersonic Mach number generally observed in the CNM (Heiles & Troland 2003) can be the result of their relative motions.

A similar explanation is proposed for the increase of the observed Mach number, ℳ_obs, with n₀ and f_CNM (see Table 1). The quantity ℳ_obs is the ratio of σ_turb to σ_therm. With the increase of the amount of cold gas (i.e., f_CNM), the average temperature of the box decreases proportionaly, as well as σ_therm. Previously we have seen that σ_turb also decreases with f_CNM but it does it more slowly than σ_therm due to the contribution of clump relative motions inherited from the WNM.

We also note that, at a given n₀, the efficiency of the CNM formation decreases when the amplitude of the turbulence driving v_S increases (see Fig. 6 – left) and thus when the Mach number increases, an effect also mentioned by Audit & Hennebelle (2005). Sánchez-Salcedo et al. (2002) and Vázquez-Semadeni et al. (2003) argue that supersonic turbulence acts faster than TI and thus dominates the evolution of the cloud, an explanation confirmed by Walch et al. (2011). As stated in Eq. (6) and as emphasized by Hennebelle & Pérault (1999), the formation of CNM occurs when the radiative time is lower than the dynamical time. Equivalently, the cooling time must be lower than the typical time of the stirring: (19)For fixed box size (L_S) and initial temperature, this offers a natural explanation why the transition is inefficient at the lowest densities and the highest amplitudes v_S.

5.2. Effect of the compressive modes

Using isothermal simulations of molecular clouds, Federrath et al. (2008, 2010) showed that the ratio of compressible to solenoidal modes of the forced velocity field has an impact on the density structure of the gas. To understand the effect of the nature of the turbulent forcing on the CNM formation in the case of a thermally bistable gas, we performed 42 simulations (128³ pixels) with different values of the projection weight ζ and a fixed initial density n₀ = 1.0 cm^-3, close to the typical WNM density. We chose an initial density that we showed does not induce CNM formation in the case of the natural state of turbulent forcing (ζ = 0.5, see Sect. 2.2) in order to see if a different combination of compressive versus solenoidal modes can trigger the phase transition at low density. Like previously, the velocity field amplitude v_S was varied from 5 to 20 km s^-1. The velocity field is either purely compressive (ζ = 0), either purely solenoidal (ζ = 1) or a mixing of both components (0.1 ≤ ζ ≤ 0.9). Table 2 summarizes the values of ℳ_obs, ℳ_theo, f_CNM and σ_turb in the stationary regime for the 42 simulations. The behavior of f_CNM and σ_turb with ζ are also presented in Fig. 10.

These results show that a turbulent field with a majority of solenoidal modes or with an energy equilibrium between compressible and vortical modes (ζ ≥ 0.5) cannot trigger the transition to CNM with n₀ = 1 cm^-3 (Fig. 10 – left). This behavior is obvious in the representations of the integrated density in Fig. 11: the more compressive modes the velocity field has, the more cold gas is formed. For a WNM initial density of n₀ = 1 cm^-3, ζ needs to be lower than 0.3 to reach 30% of the mass in the cold phase. Besides, as previously showed and for the same reason, f_CNM increases at lower values of v_S, for a given value of the spectral weight ζ (Fig. 10, – left).

We also note that σ_turb continually increases with the projection weight ζ (Fig. 10, –right), even when no cold gas is created in the simulation. The turbulent motions estimated with σ_turb are thus affected by their repartition in the compressive and solenoidal modes. σ_turb increases with the energy fraction injected in the solenoidal modes for the same total amount of energy injected. The solenoidal modes are efficient at mixing and uniformizing the density more efficiently than the compressive modes, which create higher density contrasts correlated with the velocity field (Federrath et al. 2010). The σ_turb computation gives more weight to the high density regions and thus to the compression regions where the velocity field is close to uniform, explaining why σ_turb decreases when the compressive modes become dominant. This illustrates that a higher fraction of compressible modes is not equivalent to a higher Mach number and emphasizes a clear difference with isothermal simulations (Federrath et al. 2008). The increase of the fraction of compressive modes is applied without increasing the amplitude of the velocity, and thus, without decreasing the crossing time, so that the gas has time to cool.

5.3. Robustness tests

One concern often raised about simulations done using a sustained turbulent driving in Fourier space is how it affects the dynamics of the flow and, in the specific case of thermally bistable H i, how it affects the formation of cold structures. One might be concerned with the fact that compressive modes of the forcing, applied everywhere in space and at every time step, can create unrealistic expanding or compressive flows that could enhance or inhibit the formation of cold clumps. In addition, as emphasized by Avila-Reese & Vázquez-Semadeni (2001), the driving in the ISM is highly intermittent both in time and space, so another concern is that the continuous driving of turbulence might prevent the formation of cold structures altogether.

In order to test the robustness of our results, we produced low resolution simulations (128³ cells) in decaying mode and with a purely solenoidal driving (ζ = 1) to avoid the effect of compressive modes on the formation/destruction of clumps. To test the impact of the continous forcing on the transition, the forcing was turned off at 40 Myear and, after that point, the simulation continued during 80 Myear in decaying mode. We used three initial densities: n₀ = 0.5, 1.0 and 2.0 cm^-3 and the large scale amplitude of the forcing was fixed to v_S = 7.0 km s^-1. At the end (120 My), the Mach number is stabilized below 0.1, which is much lower than the observed value in the H i but allows the study of a complete decay.

Like for the forced turbulence simulations at n₀ = 0.5 cm^-3, the decaying simulation at n₀ = 0.5 cm^-3 does not create any CNM. This result addresses one main concern: it indicates that sustained turbulence forcing in Fourier space is not responsible for the lack of CNM in simulations with initial conditions representative of the WNM in the local ISM.

The high density case (n₀ = 2.0 cm^-3) produces CNM efficiently; it reaches f_CNM = 0.35 while the turbulence is driven but, interestingly, the fraction of CNM increases significantly, up to f_CNM = 0.8, when the forcing is stopped. This is due to the fact that, the turbulence being shut down, the Mach number decreases below 0.1 and the turbulent motions, becoming negligible, do not prevent the unstable gas to cool down anymore. This effect has also been observed by Walch et al. (2011).

One interesting case is the simulation with n₀ = 1.0 cm^-3. For this simulation, the evolution with time of the mass fraction of the different phases, ℳ_obs, and σ_turb are presented in Fig. 12. We previously showed that at this density, the driving must have a majority of compressible modes to trigger the formation of CNM. As expected, no CNM and no unstable gas are created here during the driven part. However, the mass fraction of CNM reaches 60% during the decaying part of the simulation. The decay of turbulence enables the formation of unstable gas at first, which itself transits to CNM after 20 My. This is very interesting in the context of a turbulence driving that would occur by short burst of energy injection. It suggests that such intermittent driving of turbulence could be efficient at forming CNM. The problem with this picture is that the CNM is formed slowly and late after the driving has been shut down. The amount of CNM reaches 30% only when the turbulent velocity dispersion and therefore the Mach number are much lower (σ_turb = 0.4 km s^-1 and ℳ_obs = 0.06) than the measured values (3 km s^-1 and 1.0).

Finally we wanted to push further the comparison between the two types of simulations (driven turbulence and decay) in a case where they are both in the same dynamical state. We looked at the properties of the decaying simulations at a time where the Mach number is 0.4 which corresponds to a time where the driving was stopped and the turbulence is decaying (see Fig. 12). We then selected sustained turbulence driving simulations with the same n₀ and ζ = 1 that have ℳ_obs~0.4 in the stationary regime. For all values of n₀, there are no significant difference in the CNM mass fraction for decaying and sustained turbulence simulations, indicating that the process responsible for the dynamical state of the gas has limited impact on the formation of cold structures.

From the tests we have made by comparing driven and decaying turbulence simulations, we conclude that using sustained turbulence driving in Fourier space is not biasing the results presented here on the formation of the CNM.

5.4. Discussion

• 1.

  The first main result of the parameter studyon 128³ simulations is that the fiducial conditions of the WNM (n₀ = 0.2–0.5 cm^-3 and T = 8000 K) do not lead to the formation of cold gas, whatever the turbulence properties in the range we explored. With a slight increase of the density (n₀ = 1 cm^-3), a majority of compressive modes of the stirring are needed to trigger the transition. A higher increase of the initial WNM density (above n₀ = 1.5 cm^-3) produces CNM very efficiently. A moderate compression effect is therefore needed to trigger the transition with either a compressive velocity turbulent field or with an increase of the WNM density (or pressure). To first order the required increase in pressure (from P₁ to P₂) for the WNM-CNM transition could happen by the confluence of two flows with a relative velocity Δv: (20)Considering n₁ = 0.5 cm^-3 (typical density in the WNM), n₂ = 2.0 cm^-3 (efficient to trigger the transition) and T = 8000 K, one obtains Δv = 10.3 km s^-1. The pressure increase needed to place the WNM into suitable conditions to trigger the transition can happen with converging flows with a relative velocity of 10 km s^-1, which is reasonable in cases of supernovæ explosions or stellar winds, and compatible with the idea that turbulence on large scales is produced by converging flows induced by the star formation activity (Ferrière 2001; Elmegreen & Scalo 2004; McKee & Ostriker 2007; de Avillez & Breitschwerdt 2007; Kim et al. 2010, 2011). That no CNM is formed at low density, whatever the forcing conditions, can imply that typical turbulence motions only triggers the transition of the WNM to CNM if such a transition can spontaneously occur in the medium, i.e., if the WNM is already in a thermally unstable state. This means that turbulence only plays a role in the shaping of the cold medium. This is due to the density of the WNM that is the result of the segregation of the medium into WNM and CNM. Therefore the natural initial condition for the formation of CNM should not be the WNM density but the average density of the ISM. This result would suggest that the formation of CNM should be more efficient in the spiral arms.

• 2.

  The second main result is that a sub or transonic turbulence is required to get CNM gas with properties in agreement with observations (f_CNM ~ 40% and σ_turb ~ 3 km s^-1). In opposition to an isothermal gas, it is rather inefficient to inject supersonic turbulence to produce high density structure in a thermally bistable turbulent flow. This result ties up with the idea of Vázquez-Semadeni et al. (2000), Sánchez-Salcedo et al. (2002) and Audit & Hennebelle (2005) that when t_dyn is lower than t_cool, turbulent compressions and depressions act quicker than the thermal instability, preventing its development. On the contrary, when t_cool<t_dyn, the gas cools down fast enough while turbulence is compressing it that it can reach the stable branch of the CNM before its reexpansion (Hennebelle & Pérault 1999). This argument is in agreement with previous studies. Walch et al. (2011) observed that effect with similations typical of the solar neighborhood physical conditions; they also found that the two phases appear only for low Mach number. Gazol et al. (2005) also mentioned that, as the Mach number increases, the ratio t_cool/t_dyn increases as well and the structures become transient and move away from the cold branch of the thermal equilibrium. Similarly, based on simulations of H i converging flows, Heitsch et al. (2005) observed that the thermal instability is efficient when the density is high and the velocity of the flow low.

• 3.

  The third important result is related to the impact of the distribution of the turbulent energy between the solenoidal and compressive modes. Simulations with moderate initial WNM density (n₀ = 1.0 cm^-3) and purely compressive forcing field (ζ = 0) lead to a fraction of CNM close to what is deduced from observations, but in a very long time (>20 Myr). Besides, the Mach number and the turbulent velocity dispersion also need more than 10 Myr to converge to the observed values, suggesting that the turbulence takes a long time to be fully developed, namely to distribute the energy of the compressive modes toward the solenoidal modes. The convergence and the conversion of the WNM in CNM is much faster (~1 Myr) in the case of turbulence with a natural ratio of solenoidal to compressive modes but in this case a slightly higher WNM density (n₀ > 1.5 cm^-3) is required to trigger the transition. Given these timescales, we conclude that it is more likely that the formation of CNM structures occurs by an increase of the WNM density rather than purely compressive turbulent energy injection at lower density.

A second goal of this work was to study in details the scaling properties of turbulence as well as the properties of the CNM structures formed in a thermally bistable flow. To do so we produced two high resolution simulations (1024³) with initial conditions selected from the results of the low resolution parametric study. The analysis of these two simulations is described in the next section.

6. High resolution simulations

Two high resolution simulations (1024³) were produced within the same framework as the parametric study (box size of 40 pc, turbulence driving in Fourier) with initial conditions identified previously as leading to H i properties close to the observed ones.

The first one (hereafter 1024N01) has an initial density of n₀ = 1 cm^-3, a majority of compressible modes with ζ = 0.2, and a large scale velocity v_S = 12.5 km s^-1. The second one (hereafter 1024N02) has an initial density of n₀ = 2 cm^-3, high enough to allow us to study the natural state of the turbulent velocity field with ζ = 0.5, and a large scale velocity v_S = 7.5 km s^-1. These two initial conditions lead to the formation of about 40% of cold gas, a subsonic Mach number and a turbulent velocity dispersion between 2 and 3 km s^-1.

The initial densities chosen are higher than the observed density of the WNM (see Sect. 3.2) in accordance with our previous results showing that some compression of the WNM is needed in order to trigger the formation of CNM. On the other hand, as about 40% of the mass is converted to CNM with a volume filling factor of about 1%, the density of the WNM in the stationary regime is about 60% the initial value and therefore closer to the observed WNM density.

All the results presented hereafter concern these two 1024³ cells simulations after 16 Myr. The integrated density along the z-axis of both simulations are shown in Fig. 13. These maps reveal structures in clumps and filaments comparable to what is seen in observations of the diffuse ISM. From this first look we already notice that 1024N02 creates more structures than 1024N01, suggesting that it is more efficient at triggering the transition. This will be discussed later.

6.1. Mach number

Table 3 summarizes the values of ℳ_theo, computed on all the pixels and for pixels with T > 200 K, ℳ_obs, and σ_turb. The values of ℳ_theo indicate that 1024N01 is subsonic while 1024N02 is transonic. Besides, removing the cold pixels (T < 200 K) from the computation does not change the results, suggesting that turbulent motions are dominated by the WNM motions, as we will discuss later. On the other hand, the values of ℳ_obs are clearly subsonic for both simulations, in agreement with observations. Lastly, both simulations reproduce very well the observed properties of σ_turb. In fact the computed values of σ_turb are slightly lower than the estimated value in a WNM cloud of 40 pc (~3 km s^-1) but this is due to the cold structures present in the simulations for which the velocity dispersion is much lower.

6.2. Thermal distribution of the gas

As expected with the thermal instability at low Mach number, the distribution of the gas in the P,n,T space presents the evidence of a multiphasic medium. The histograms of T (Fig. 14) have a clear bimodal shape, with a peak in the warm phase and one in the cold phase. Similarly the mass distribution of the gas in the P − n diagram (Fig. 15) shows two preferential zones: the first one at high density and the second one at low density, suggesting again CNM and WNM.

However, the thermal properties of the warm gas are slightly different in both simulations. If the histograms of T of the cold gas peaks at similar temperatures in both cases (40 K for 1024N01 and 50 K for 1024N02), it is not true for the warm gas: while the WNM of 1024N01 has ⟨ T ⟩ ~ 7000 K, close to the prediction of the two-phase model (Field et al. 1969), the warm gas in 1024N02 has ⟨ T ⟩ ~ 3500 K. This effect can also be seen in the pressure-density diagrams (Fig. 15), where the pressure of 1024N02 is lower than the one of 1024N01 at low densities. Therefore, while the gas at low density of 1024N01 follows the stable branch of the WNM, the gas at low density of 1024N02 is shifted below the equilibrium curve where the heating is dominant. This pressure difference is confirmed by the histograms of P presented in Fig. 17, where the mean pressure of 1024N02 (right) is lower than the mean pressure of 1024N01 (left). On the other hand, the histograms of n presented in Fig. 18 peaks around 0.5−0.6 cm^-3 in both cases, similar to the WNM mean density estimated from observations.

Two factors can be considered to explain the presence of warm gas out of equilibrium. First, a higher initial density decreases the WNM cooling time (Eq. (8)): t_cool goes from 1.75 Myr for n₀ = 1.0 cm^-3 to 1.15 Myr for n₀ = 2.0 cm^-3. Therefore, the warm gas located below the equilibrium in 1024N02 is cooled down faster than in 1024N01 and does not have time to reach the equilibrium curve. Secondly, the thermal equilibrium curve considered in the P − n diagram is established for a static gas and does not take the energy contribution of turbulence into account. Depending on the properties of the turbulence, WNM gas can stabilized at a temperature significantly lower than the prediction of the two-phase model.

The definition of the H i thermal phases in the physical parameter space is not always clear in the litterature; it depends somewhat if they are looked at from the simulation or the observation perspective. Here we used three criteria to estimate the mass fraction of the CNM, WNM and thermally unstable gas, based on T and n thresholds. All results are summarized in Table 4.

Whatever criteria we used to estimate the CNM (T or n), the fraction of the total mass in the CNM is 30% and 50% for 1024N01 and 1024N02 respectively. These quantities are in good agreement with the results obtained with the same conditions but on a 128³ grid, providing another validation of the use of low-resolution simulations for parameter space exploration. These values of f_CNM are also in agreement with observational results: Heiles & Troland (2003) estimated that 40% of the mass of the H i lies in the cold phase using 21 cm emission and absorption measurements, and we found f_CNM~ 44% based on large scale properties of the H i emission (see Sect. 3.4). The CNM volume filling factor found in these high resolution simulations (1% for 1024N01 and 4% for 1024N02) are also representative of what we estimated from observations.

The estimate of the mass in the WNM and thermally unstable gas depends significantly on the criteria used. One way to separate the phases is to use the density thresholds that define the range of the thermally unstable gas, as described in Sect. 3.1 (see the dash-dotted lines in Fig. 2). Using that criteria, both simulations hold the same amount of unstable gas (~30%), the difference between the two simulations only lies in the CNM and WNM phases: the amount of gas that does not transfer to CNM in 1024N01 stays in the WNM phase. We note that in the case of 1024N01, the three thermal phases are distributed like 1/3–1/3–1/3, as observed by Heiles & Troland (2003), while in the 1024N02 case, the mass fraction of WNM is lower of a factor 2 and 50% of the mass is in the CNM.

One can also estimate the mass fraction of the WNM and thermally unstable gas using a threshold in T. In both simulations we found that about 20% of its mass is lying between 2000 and 5000 K. In 1024N02, where the WNM-CNM transition is highly efficient, the unstable regime holds the major part of the mass (46% between 200 and 5000 K) and half of it has a temperature lower than 2000 K; in 1024N01, where the WNM-CNM transition is less efficient, only 10% is lying between 200 and 2000 K and a significant amount of the gas lies in the WNM range beyond 5000 K (40%).

Finally it is interesting to look at the spatial distribution of the three phases. In Fig. 16 we present 2D density slices for each simulation with isocontours placed at 200 K (blue) and 2000 K (red). The highest temperature threshold has been consciously chosen lower than the usual limit at 5000 K in order to see the distribution of the thermally unstable gas. The gas lying between 2000 and 5000 K is widely distributed and does not show any preferential location. On the other hand, the thermally unstable gas below 2000 K lies preferentially around the cold structures as sheets around the denser regions, in agreement with Gazol & Kim (2010).

6.3. Pressure range

Jenkins & Tripp (2011) used UV absorption measurements of neutral carbon observed in front of stars to study the distribution of the thermal pressure in the diffuse, CNM. They found that the pressure distribution, weighted by the mass density, is well fitted by a lognormal. It is also the case for our two high-resolution simulations. Figure 17 shows the histogram of P weighted by n (black line) and the average pressure along the line-of-sight weighted by density in each cell (blue line), this quantity being closer to what is measured with real observations for which the 3-dimensional pressure is not available. For comparison purpose, the pressure distribution found by Jenkins & Tripp (2011) is also shown in these plots (red dotted line).

The pressure histograms of Fig. 17 were fitted with a lognormal, following the definition of Jenkins & Tripp (2011): (21)The values found for our lognormal fits are given in the caption of Fig. 17. The simulation 1024N01 peaks at 3800 K cm^-3 (3.58 dex), the same value found by Jenkins & Tripp (2011). The pressure peak of 1024N02 is slightly lower at 2200 K cm^-3 (3.34 dex). These two values lie in the pressure range allowing a two-phase medium. We point out that the width of the histograms are broader for the 3D pressure histograms than for the averaged pressure in both cases (0.29 and 0.20 dex compared to 0.14 and 0.09 dex respectively) due to the average along the line of sight that compresses the dynamical range. In addition, the distributions of 1024N01 are wider than those of 1024N02, which can be explained by the higher amplitude of the turbulent stirring used in 1024N01. In any case, the widths of the pressure histrograms stay close to the result of Jenkins & Tripp (2011) and are thus in good agreement with the observations.

6.4. Density histograms

The density distribution is of great importance in turbulence studies. Many simulations of molecular clouds have shown that the density distribution of an isothermal gas follows a lognormal law with a widening that depends on the Mach number (Vázquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Federrath et al. 2008). The density distribution of a two-phase gas is however more complex and often presents a bimodal shape (Hennebelle et al. 2008; Audit & Hennebelle 2005, 2010). Vázquez-Semadeni et al. (2000), Piontek & Ostriker (2005), Gazol et al. (2005) and Walch et al. (2011) showed that density histograms are broadened and loose the bimodal shape when the Mach number increases, as most of the gas moves to the thermally unstable regime. The density distributions presented in Fig. 18 point to the evidence of a transonic biphasic gas with two distinct peaks, one at low densities, 0.4–0.5 cm^-3, typical of the WNM, and a second one around 10–20 cm^-3, representative of the CNM.

To compare with studies of molecular clouds, we tried to fit lognormals on the cold peaks for each simulation. Regarding the high density part, it is possible to fit the right tail with a lognormal, but for both simulations, the lognormals are overly broad. To study this further, we separated the cells of the box where T < 200 K and plotted the histogram of the CNM density. They did not show any characteristic shape, and it was impossible to fit lognormals on them. We also didn’t find any evidence of a power law tail as Seifried et al. (2011) did. Therefore, even though a lognormal fit is possible on the high density part of the histogram, it is difficult to make any physical interpretation out of it.

6.5. Power spectra

Compared with converging flow simulations, the stirring in Fourier provides a wider range of scales over which the properties of turbulence can be studied. In particular, such simulations offer an opportunity to better study the velocity and density power spectra and therefore characterize some specific properties of a turbulent and thermally bistable medium. In this section we present a study of the statistical properties of the density and velocity fields of 1024N01 and 1024N02. For comparison purposes, we also made the same analysis on an isothermal simulation of 1024³ cells, obtained with the same turbulent stirring than 1024N01 (ζ = 0.2 and v_S = 12 km s^-1), an initial conditions of n₀ = 1 cm^-3 and T = 8000 K. This simulation is subsonic with a Mach number equal to 0.9.

For each 2D plane perpendicular to the z-axis, we calculated the density contrast σ(n)/⟨ n ⟩, the maximum of the density and the slope of the 2D power spectra of the density and the velocity. All quantities were computed for the set of points on each 2D plane perpendicular to the z-axis. The choice of the axis is not important as the simulations are isotropic due to the turbulent stirring in Fourier space. The advantage of computing quantities for each 2D plane is that it offers a way to quantify their variations across the cube. For a scalar field like the density, the power spectrum of 2D slices is equivalent to the power spectrum of the 3D cube (Miville-Deschênes et al. 2003b)³. We chose to compute the statistics of one component of the velocity (v_z) and not its norm to stay closer to spectro-imagery observations that give only access to the velocity component along the line of sight. We also chose not to do the statistics of 21 cm velocity channel maps as their structure is a mixture of velocity and density fluctuations, especially in high density contrast media (Esquivel et al. 2007).

The density contrast, the maximum density, and the slopes of the power spectra as a function of z, are shown for the isothermal and TI simulations in Fig. 19. We also computed the mean power spectra of the density and the velocity, by taking the mean of the 1024 power spectra computed on each plane perpendicular to the z-axis. Figure 20 shows these compensated mean power spectra k^{8 / 3}P(k) for 1024N01 and for the isothermal simulation (the power spectra for 1024N02 being similar to those of 1024N01, they are not displayed here). The errors bars represent the (1σ) variations on the 1024 slices.

The estimation of the slope of the power spectrum in the inertial zone of numerically simulated turbulence can only be done on a small scale range due to the turbulent forcing on large scales and the numerical dissipation at small scales. We evaluated the inertial range [ k_min,k_max ] for each mean power spectrum by minimizing the root mean square of the error (22)where a and β are the two parameters obtained with the power law fit. The size of the simulations (1024³) allows us to expect at least a decade on which a power law fit is suitable. Therefore we imposed the following condition on k_min and k_max: k_max ≥ 10 × k_min. The slopes β as well as the k-ranges used for the different power spectra are summarized in Table 5.

The density power spectrum of the isothermal and subsonic simulation has a slope of –8/3, which corresponds to the Kolmogorov estimation in 2D (see Elmegreen & Scalo 2004). This Kolmogorov law is in agreement with the previous results of Kim & Ryu (2005) and Kritsuk et al. (2007) who emphasized that subsonic and isothermal simulations have a Kolmogorov behavior. Regarding the TI simulations, the density power spectra are far from Kolmogorov with slopes of –1.3 for 1024N01 and –1.5 for 1024N02. This is due to the much higher density contrast. Its values are indeed around 4 (⟨ C_n01 ⟩ = 4.6 ± 1.4 and ⟨ C_n02 ⟩ = 3.8 ± 0.7) while the density contrast of the isothermal simulation is almost constant from slice to slice with a value of about 0.6 (Fig. 19 – left).

As expected, the velocity power spectrum of the isothermal and subsonic simulation also follows the Kolmogorov estimation. However, the velocity power spectra of the two TI simulations are slightly flatter from Kolmogorov with slopes of –2.4 for 1024N01 and –2.3 for 1024N02.

In supersonic turbulent flows it is generally observed that the slope of the density power spectrum flattens with the increase of the Mach number. This effect has been reported for isothermal (Kim & Ryu 2005) and supersonic, thermally bistable (Gazol & Kim 2010) flows. It was also noted by Kritsuk et al. (2007) that the velocity power spectrum goes from Kolmogorov-like to Burgers-like when the Mach number increases. We observe that the density power spectrum of transonic, thermally bistable turbulence has the same behavior than for supersonic turbulence. On the other hand, while the velocity power spectrum of supersonic flows goes from Kolmogorov to Burgers turbulence (Kritsuk et al. 2007), the velocity power spectra in our case stay similar to the one of a subsonic flow with a slope close to Kolmogorov. In fact the velocity power spectra of the TI simulations are even slightly flatter than Kolmogorov, which is the opposite behavior expected for Burgers turbulence. Hennebelle & Audit (2007) emphasized that the density fluctuations of the H i, and thus the density contrast, are due to the thermal instability at low Mach numbers while supersonic shocks dominate at high Mach numbers. The TI simulations are indeed much more structured on the small scales than subsonic isothermal flows. We also want to point out that the power spectrum of the logarithm of the density, which traces more directly the velocity power spectra through the equation of continuity, has a slope very close to the Kolmogorov law for both 1024N01 and 1024N02.

6.6. CNM structures

Two observational facts often attributed to the self-similarity induced by the turbulent cascade in the ISM are the mass-scale (M ∝ L^α – Larson 1981; Elmegreen & Falgarone 1996; Roman-Duval et al. 2010) and velocity dispersion-scale (σ_turb = σ_{1 pc}L^γ – Larson 1981; Heyer & Brunt 2004) relations observed in molecular clouds. If CNM clumps are the precursors of molecular gas, it is interesting to see if they follow similar laws. To identify CNM clumps in the simulations we selected pixels having a density greater than 5 cm^-3. Structures were then identified as spatially connected pixels in 3D. Following the approach of Hennebelle et al. (2007), the size of a structure is defined using the highest eigenvalue λ₁ of its inertial matrix ℐ: The x_i, i ∈ { 1,2,3 }, are the pixels coordinates and dm the infinitesimal mass. One can thus estimate the size of the structure by computing , where M is its mass. We point out that Audit & Hennebelle (2010) investigated other choices than the highest eigenvalue (the geometrical mean of the eigenvalues for instance) and found very similar results.

Figure 21 – left shows the relation between M and L for the cold structures of 1024N01 (this relation for 1024N02 being similar to that of 1024N01, it is not displayed here). There is significant dispersion in this diagram but the points suggest a power law fit. We found slopes of 2.25 ± 0.05 and 2.28 ± 0.03 for 1024N01 and 1024N02 respectively. These results are in agreement with previous works on isothermal (Kritsuk et al. 2007; Federrath et al. 2009) and TI simulations (Audit & Hennebelle 2010). They are also similar to the results obtained on CO observations by Elmegreen & Falgarone (1996) who found slopes between 2.2 and 2.5 on a large sample of independent molecular clouds. Recently, Roman-Duval et al. (2010) measured an index of 2.36 on 580 molecular clouds. The structure observed here on H i simulations is thus similar to molecular clouds⁴.

Similarly we investigated the relation between the velocity dispersion and the size of the CNM clumps. The velocity dispersion is computed as: (26)with v₀ the mean velocity of each structure weighted by density: (27)Figure 21 (right) shows the values of σ( | v |) versus L for the structures of 1024N01 (1024N02 shows similar results). The power law fits give 0.70 × L^0.41 and 0.51 × L^0.42 for 1024N01 and 1024N02 respectively. We note the presence of a scatter of one order of magnitude around the fitted power law, as it as been previously noticed on simulations by Vázquez-Semadeni et al. (1997) and on observations by Falgarone et al. (1992). Keeping in mind that the slopes obtained here (0.41 and 0.42) may be subject to high uncertainty, they are compatible with the results of Audit & Hennebelle (2010) who found 0.33 < γ < 0.53. They are also close to values measured on observations: 0.37 in H i clouds (Larson 1979) and molecular clouds (Larson 1981), and 0.5 obtained on a large sample of molecular clouds by Heyer & Brunt (2004). On the other hand, the values obtained here for σ_{1 pc}, 0.70 and 0.51 km s^-1, are lower than the results presented by Audit & Hennebelle (2010) who found values between 1.2 and 3.3 km s^-1 but with a different forcing scheme. The turbulence in their simulation is indeed produced by WNM converging flows with a Mach number close to 1.5. This type of forcing is very efficient to create cold structures but induces a Mach number in the WNM higher than the one observed in the diffuse ISM.

In order to compare σ_{1 pc} to the work of Wolfire et al. (2003) and to the value 0.92 km s^-1 estimated in Sect. 3, we computed for each CNM structure using the theoretical value γ = 1/3 expected for subsonic turbulence. The histograms of σ(1) presented in Fig. 22 peak around 0.8 and 0.6 km s^-1 for 1024N01 and 1024N02 respectively. These values are in agreement with the velocity dispersion at 1 pc we wanted to reproduce in these simulations and with the work of Wolfire et al. (2003). That σ_{1 pc} is well reproduced suggests that the internal dynamics of the CNM clumps is related to the dynamics in the WNM. To confirm this assertion, we computed the mean velocity dispersion ⟨ σ( | v |) ⟩ (i.e. the average over all cold structures of their individual turbulent velocity dispersion). We obtained 0.3 and 0.4 km s^-1, much lower than the mean thermal velocities of the clumps: 1.7 km s^-1 for 1024N01 and 1.3 km s^-1 for 1024N02. The individual cold structures are therefore subsonic. On the other hand, the relative motions between the clumps are greater. We computed the cloud-cloud velocity dispersion as (28)where N is the number of clumps. We obtained σ_cc = 2.3 km s^-1 for 1024N01 and σ_cc = 1.9 km s^-1. These values are very close to the total velocity dispersions in the boxes (2.8 km s^-1 for 1024N01 and 2.2 km s^-1 for 1024N02), suggesting again that the clump motions are related to the motions in the WNM. This means that the CNM structures keep the memory of the flow in which they were formed and that little dissipation occurs in cloud-cloud collisions.

Part of this relation between CNM and WNM motions could be attributed to the turbulent forcing being applied as an acceleration instead of a force. In the context of a two-phase medium this might introduce biases, in particular an overestimate of the clump-to-clump velocity, because dense clumps are subject to the same accelerations as diffuse gas in the warm medium. Nevertheless, similar velocity dispersion of CNM and WNM structures were also reported by Heitsch et al. (2005, 2006) who simulated the formation of molecular clouds from H i converging flows and therefore did not use turbulent forcing in Fourier. Like them, we thus conclude that a significant fraction of the observed 21 cm line broadening is due to the relative motions of the clumps and is not the result of internal supersonic motions.

6.7. Summary and discussion

Here we summarize the physical properties of the CNM deduced from the analysis of two TI simulations at high resolution (1024³, box size of 40 pc, resolution of 0.04 pc).

• 1.

  First we showed that these simulations reproducewell all the observed physical quantities used as con-straints for the parametric study: a cold mass fraction ofabout 40%, a CNM volume filling fac-tor of 1 to 4%, a transonicMach number and a turbulent velocity dispersion ofabout 3 km s^-1.Furthermore, like in the observations reported by Heiles &Troland (2003), about 1/3of the gas is in the thermally unstable regime. This re-sult also agrees with the previous works on simulationsof Gazol et al. (2005)and de Avillez & Breitschw-erdt (2005). Finally, the pressures are invery good agreement with the observations of Jenkins &Tripp (2011), both for the mean and the disper-sion.

• 2.

  We noted that the thermally unstable gas, defined with temperatures between 200 and 5000 K behave in two different ways, depending on its temperature. The gas between 200 and 2000 K is located at the interface of the cold clumps and the warm gas while the gas beyond 2000 K is dynamically connected to the WNM and is widely distributed in the box. This wide distribution of the unstable gas temperature has also been noticed by Gazol & Kim (2010) who observed that the boundaries between CNM and WNM are sharper when the Mach number is very low (0.2), due to the unperturbed development of the thermal instability by turbulence. This is a natural consequence of the clumps being dynamically formed by locally converging motions. The gas is cooling while approaching the clumps, it is therefore natural that the coldest unstable gas is located around them.

• 3.

  The density distribution is bimodal, bearing the characteristic of a bistable gas, as noticed by several authors (Vázquez-Semadeni et al. 2000; Audit & Hennebelle 2005; Gazol et al. 2005; Piontek & Ostriker 2004; Audit & Hennebelle 2010). This is a result of the gas being transonic. Indeed, a higher turbulence would maintain the gas in the thermally unstable regime (Gazol et al. 2005; Walch et al. 2011) and break the bimodality of the distribution that tends to be close to a single peak lognormal function. As mentioned by Sánchez-Salcedo et al. (2002) and Vázquez-Semadeni et al. (2003), this is due to the fact that for strong turbulence, the gas is not in the thermal equilibrium locus of the pressure-density phase diagram but is rather mostly in the thermally unstable range. This might not be representative of the solar neighborhood where the WNM is thought to be transonic, implying a bimodal histogram of the density. In this context, we tried to describe the high density part of the density histogram with a lognormal or a power law tail but showed that none of them were convincing.

• 4.

  The 2D velocity power spectra of the simulations with thermal instability have slopes equals to –2.4 and –2.3 for 1024N01 and 1024N02 respectively, close to the Kolmogorov law (–8/3). On the other hand, the density power spectra are much flatter (–1.3 and –1.5) and the density contrasts are high (around 4). These signatures are the direct result of the formation of dense structures by the thermal instability in a low Mach number flow (Hennebelle & Audit 2007). We point out that the density power spectra computed on the TI simulations are similar to the one seen in supersonic isothermal flows where the density power spectrum flattens with the increase of the Mach number (Kim & Ryu 2005) while the velocity power spectra are characteristic of subsonic turbulence with a Kolmogorov slope.

• 5.

  Cold clumps with a density higher than 5 cm^-3 follow a mass-scale relation M ∝ L^α, with α = 2.25 and 2.28, in good agreement with observations (Elmegreen & Falgarone 1996; Roman-Duval et al. 2010) as well as with previous numerical studies, even when different conditions and algorithm are used (Kritsuk et al. 2007; Federrath et al. 2009; Audit & Hennebelle 2010). The cold clumps also follow a velocity dispersion-scale relation (σ ∝ L^γ) with γ = 0.41 and 0.42. These values are similar to the ones obtained on numerical simulations (Audit & Hennebelle 2010), to the observed in molecular clouds (Larson 1981; Heyer & Brunt 2004), but also to the value expected for a subsonic turbulence γ = 1/3. These results suggest that the self-similar density structures observed in molecular clouds could be inherited from the H i from which they were formed.

• 6.

  The comparison of the internal velocity dispersion of the cold structures (0.3–0.4 km s^-1) with their thermal velocity dispersion (1.3–1.7 km s^-1) indicates that they are clearly subsonic. The turbulent velocity dispersion of clumps increases with scale with a normalization at 1 pc of σ(1) = 0.6−0.8 km s^-1, very close to the value found in the WNM. Similarly, the cloud-cloud velocity dispersion is 1.9–2.3 km s^-1, close to the total velocity dispersions in the boxes (2.2–2.8 km s^-1). This implies that the turbulent motions inside clumps and the relative velocities between them are related to the motions of the WNM from which they were formed. Individually, the clouds are subsonic but their relative velocities are supersonic relative to the CNM sound speed (in accordance with Koyama & Inutsuka 2002; Heitsch et al. 2005). This result is in complete agreement with the work of Hennebelle & Audit (2007) who noticed that the CNM velocity dispersion estimated from 21 cm spectra is greater than the mean velocity dispersion of the CNM structures identified in their simulation. This means that a non negligible fraction of the measured velocity dispersion is caused by the relative motions of the clumps along the line of sight, suggesting that the observed line broadening is likely to be due to the relative clumps motions rather than supersonic turbulence.

7. Conclusion

We have presented a set of 90 low resolution simulations (128³ pixels) in order to study the physical conditions that lead to the production of CNM structures out of purely WNM gas, given the known observational constraints (velocity dispersion, density, pressure) of the H i in the local ISM and based on the model of the heating and cooling processes of Wolfire et al. (1995, 2003). These simulations show that WNM gas at typical values of density, temperature and turbulent motions do not transit to CNM gas. The typical volume density of the WNM in the disk at the position of the Sun lies between 0.2 and 0.5 cm^-3. We showed that WNM with initial density under 1 cm^-3 never transit into a cold phase, whatever the amplitude or the nature of turbulent motions. Stronger turbulent motions clearly does not help in producing cold gas, in fact we showed that it is the opposite. On the other hand a moderate increase of the WNM density (or equivalently of its pressure) is very efficient; all simulations with an initial density higher than 1.5 cm^-3 do transit and the amount of cold gas in the simulation increases with the density. When the initial density reaches 3 cm^-3, all simulations achieve to create at least 90% of CNM. These results were obtained with a natural balance between compressive and solenoidal modes in the injected velocity field. In this case the formation of CNM gas is almost instantaneous as soon as the turbulence is fully developed. We also explored the case of an intermediate initial density (1 cm^-3) with a variation of the ratio of solenoidal to compressive mode in the turbulent forcing. We showed that, at this density, the transition into CNM only occurs when there is a majority of compressible modes in the turbulent velocity field but this process is much slower (10–20 Myr).

We conclude that the ISM turbulence cascade on its own cannot induce a transition from WNM into CNM at the average density of the WNM. The production of CNM is likely to be due to turbulent motions of moderate amplitude associated to a compressive event of the WNM that could be led by transient phenomena such as outflows, supernova explosion or spiral density waves. This would trigger the phase transition without delay. Once formed, the amount of CNM gas is dynamically stable. Therefore there is no need for a constant energy injection to maintain the high density contrast of the H i.

We also have presented a detailed study of the H i based on two higher resolution simulations (1024³ pixels). To our knowledge these are the largest simulations to date that include the thermal instability and a pseudo-spectral turbulent stirring. The bimodal temperature and density histograms show the evidence of a bistable medium and the mass fractions of each phase, the Mach number of the WNM, and the pressure distributions are all in agreement with observations. With these simulations we confirm a number of results of previous studies based on converging flows (Heitsch et al. 2006; Hennebelle & Audit 2007; Audit & Hennebelle 2010) that showed that the structure of the CNM and WNM are tightly interwoven: the two media share the same velocity fields and the CNM cloud-cloud velocity dispersions is close to the WNM sound speed. While the individual CNM structures are clearly subsonic, their relative motions are supersonic with respect to their own temperature. In addition, thanks to the pseudo-spectral forcing method, the statistical properties of the turbulent bistable flow could be analyzed on a wider range of scales than could be done in previous studies. We found that for such subsonic or transonic bistable flows, the velocity field keeps the properties of subsonic turbulence with a power spectrum close the Kolmogorov while the density is highly contrasted, reminiscent of what is observed for supersonic flows and in molecular clouds. In the case of the H i simulated here, the joined action of the turbulence and the thermal instability allows the formation of long-lasting cold and dense structures with no need for supersonic motions, providing a favorable terrain for the formation of highly contrasted molecular clouds.

---

Acknowledgments

We thank Enrique Vázquez-Semadeni, the referee, for the critical reading of this paper, his input and his help in clarifying many points of this work. Eléonore Saury acknowledges support from the MISTIC ERC grant #267934. Heracles (González et al. 2007) is available on the following link: http://irfu.cea.fr/Projects/Site_heracles/. The simulations at low resolution have been performed on the cluster Sunnyvale at CITA, University of Toronto (Canada), and the 1024³ cells simulations have been performed on the cluster Jade, Cines, Montpellier (France).

References

All Tables

All Figures