Forced turbulence in thermally bistable gas: a parameter study
^{1}
Institut für Theoretische Astrophysik, Universität Heidelberg,
AlbertUeberleStr. 2,
69120
Heidelberg,
Germany
email: dseifried@ita.uniheidelberg.de
^{2}
Institut für Astrophysik Göttingen, Universität Göttingen,
FriedrichHundPlatz
1, 37077
Göttingen,
Germany
^{3}
Lehrstuhl für Astronomie, Institut für theoretische Physik und
Astrophysik, Universität Würzburg, Am Hubland, 97074
Würzburg,
Germany
Received:
8
March
2010
Accepted:
3
September
2010
Context. Thermal instability is one of the dynamical agents for turbulence in the diffuse interstellar medium, where both turbulence and thermal instability interact in a highly nonlinear manner.
Aims. We study basic properties of turbulence in thermally bistable gas for variable simulation parameters. The resulting cold gas fractions can be applied as parameterisation in simulations on galactic scales.
Methods. Turbulent flow is induced on large scales by means of compressive stochastic forcing in a periodic box. The compressible Euler equations with constant UV heating and a parameterised cooling function are solved on uniform grids. We investigate several values of the mean density of the gas and different magnitudes of the forcing. For comparison with other numerical studies, solenoidal forcing is applied as well.
Results. After a transient phase, we observe that a state of statistically stationary turbulence is approached. Compressive forcing generally produces a twophase medium, with a decreasing fraction of cold gas for increasing forcing strength. This behaviour can be explained on the basis of turbulent mixing. We also find powerlaw tails of probability density functions of the gas density in highresolution runs. Solenoidal forcing, on the other hand, appears to prevent the evolution into a twophasemedium for certain parameter regions.
Conclusions. The dynamics of thermally bistable turbulence show a substantial sensitivity to the initial state and the forcing properties.
Key words: hydrodynamics / instabilities / turbulence / methods: numerical / ISM: kinematics and dynamics
© ESO, 2010
1. Introduction
Observations of the interstellar medium (ISM) show large variations in the thermodynamical properties of the gas, which led to the famous threephasemodel of McKee & Ostriker (1977). In this article, we focus on the warm neutral medium (WNM) and the cold neutral medium (CNM), whose occurrence was first explained consistently by Field (1965) with the help of the thermal instability (TI). The classical conception of a more or less static multiphase medium with different phases coexisting in pressure equilibrium (see, for instance, Wolfire et al. 1995) has nowadays been replaced by a picture of a highly turbulent ISM (e.g. BallesterosParedes et al. 2007; McKee & Ostriker 2007). Several observational results point to turbulence in the ISM, for example, the broadening of spectral lines, chemical mixing of the ISM, or the star formation rate (for an overview, see Mac Low & Klessen 2004; Elmegreen & Scalo 2004). As turbulence decays roughly on one dynamical timescale (e.g. Mac Low et al. 1998; Stone et al. 1998), energy has to be injected into the ISM to maintain turbulent motion. Several different sources were proposed to drive turbulence in the ISM, e.g., blast waves from supernova explosions, the differential rotation of galaxies, but also thermal or gravitational instabilities.
While the properties of forced turbulence in the molecular and nearly isothermal phase have been studied extensively (e.g. Kritsuk et al. 2007; Schmidt et al. 2009; Federrath et al. 2010), there is no corresponding analysis for the diffusive ISM so far. Previous work on turbulence in thermally bistable gas encompass decaying turbulence (Kritsuk & Norman 2002), simulations on kpc scales to study the effect of stellar feedback (Joung & Mac Low 2006; Joung et al. 2009), and colliding flows of atomic gas subject to the thermal instability (Heitsch et al. 2006; VázquezSemadeni et al. 2007; Hennebelle & Audit 2007; Banerjee et al. 2009; Audit & Hennebelle 2010; Folini et al. 2010). The last scenario has attracted a lot of attention as a possible model for the formation of molecular clouds from the diffuse ISM. Furthermore, simulations of turbulence produced by purely solenoidal forcing in thermally bistable gas were performed by VázquezSemadeni et al. (2003), Gazol et al. (2005), Kissmann et al. (2008), and Gazol et al. (2009). In these simulations, large amounts of thermally unstable gas (up to 50%) were found, which is also indicated by observational studies (e.g. Dickey et al. 1977; Spitzer & Fitzpatrick 1995; Heiles 2001; Heiles & Troland 2003; Kanekar et al. 2003).
As was demonstrated by Schmidt et al. (2009) for the isothermal case, compressive stochastic forcing produces a flow structure on large scales that resembles colliding flows, because converging shocks are formed between rarefied gas regions. In this work, we study the properties of compressively driven turbulence in thermally bistable gas, depending on the magnitude of the force and the initial state of the gas, which is given by the equilibrium temperature for the mean gas density. In this regard, our approach is similar to the work of Gazol et al. (2005). The mean density is constant because of the periodic boundary conditions. The cooling function follows the definition of Audit & Hennebelle (2005). As Federrath et al. (2008) reported a substantial sensitivity of the density field in isothermal supersonic turbulence on the forcing (see also Federrath et al. 2009, 2010), we applied solenoidal forcing as well. Because of the threedimensional parameter space (mean density, magnitude, and ratio of solenoidal to compressive modes of the forcing), the aim of the present study is to investigate the general behaviour of forced thermally bistable turbulence by means of a series of moderateresolution simulations, using uniform grids with 256^{3} cells. To avoid a further increase in dimensionality of the parameter space, we omit magnetic fields, selfgravity, and thermal conduction. Nevertheless, our results serve as a useful guide estimating the influence of unresolved thermal processes in galactic disk simulations (Tasker & Bryan 2006; Dobbs et al. 2008; Agertz et al. 2009). In such simulations, the cutoff length at the highest level of refinement is typically of the same order of magnitude as the forcing length scale in the periodicbox simulations. In the various prescriptions of the local star formation efficiency, it is neglected that the fraction of cold gas that can turn into stars via gravitational collapse is less than unity. We therefore investigate the influence of the forcing and the mean density on this cold gas fraction in detail. Future work will concentrate on particular parameter sets to obtain the turbulent scaling and fragmentation properties from highresolution simulations and to study the influence of additional physics.
The outline of this article is as follows. Section 2 briefly describes the numerical techniques. The parameters, for which we run simulations, are specified in Sect. 3. The phenomenology of turbulence in thermally bistable gas is described for a representative parameter set in Sect. 4. For this parameter set, the grid resolution was varied from 128^{3} to 512^{3} to check for the resolution dependency of the results. Next, we analyse the variation of statistical properties for different simulation parameters. In Sect. 5, the results are discussed in the context of other numerical studies and observations, followed by a summary of the main results and an outlook in Sect. 6.
2. Numerical techniques
The simulation presented in this article are performed with the open source code Enzo (O’Shea et al. 2005). We solve the compressible Euler equations using the piecewiseparabolic method (PPM) of Colella & Woodward (1984). Besides an external driving force f, which is responsible for generating turbulent motions, we include different heating and cooling terms, combined into a net cooling rate per unit mass ℒ. The resulting equations for the mass density ρ, the velocity u, and the specific total energy e of the fluid can be written as where the total time derivative is defined by (4)The specific total energy e is given by (5)where γ = 5/3 is the adiabatic exponent of monatomic gas, and the gas pressure P is related to the mass density and the temperature via the ideal gas law: (6)The constants k_{B}, μ, and m_{H} denote the Boltzmann constant, the mean molecular weight, and the mass of a hydrogen atom, respectively.
The cooling function ρℒ was defined by Audit & Hennebelle (2005). It includes the finestructure cooling of C ii and O i, as well as the cooling by H (Lyαline) and by electron recombination onto charged grains. The only heating process considered is the photoelectric effect on small grains and polycyclic aromatic hydrocarbons (PAH) caused by farultraviolet galactic background radiation. For detailed information about the different processes see Wolfire et al. (1995, 2003); Spitzer (1978) and Bakes & Tielens (1994). The cooling is implemented in the form of an explicit scheme into Enzo (Niklaus et al. 2009). After each time step, the state variables are updated in several subcycles by adding the resulting total energy increment. As the considered cooling and heating processes are only well defined in the diffuse ISM, the calculation of is constrained to temperatures [10 K, 10 000 K]. For gas with temperatures above 10 000 K, we use ℒ(ρ,10 000 K), which leads to an underestimate of the emitted energy, hence to warmer gas than expected. In strongly driven simulations we apply a temperature cutoff at 30 000 K, assuming that the gas would cool down very fast beyond this limit. Since the speed of sound increases with the temperature, the cutoff allows us to avoid very small time steps arising from unphysical values of the temperature outside the range of the cooling function. The concerned simulation runs are mentioned in Sect. 4.
The specific force f in Eqs. (2) and (3), which generates largescale motions of the gas, is composed in Fourier space. Each mode of the force field is given by a stochastic process that is based on the OrnsteinUhlenbeck process (Eswaran & Pope 1988; Schmidt et al. 2006). For the details of this method, see Schmidt et al. (2009). The forcing varies on the autocorrelation timescale T, which is given by the ratio of the integral length L to the characteristic velocity V of the largescale motions. We set the integral length to half of the box size. Thus, the wavelengths of the forcing modes range from L/2 to 2L. The free parameter V specifies the magnitude of the forcing, i.e., the specific force is ~V^{2}/L. By applying a Helmholtz decomposition in Fourier space, the ratio of solenoidal (divergencefree) to compressive (rotationfree) modes can be arbitrarily adjusted. We use the decomposition parameter ζ such that ζ = 1 corresponds to purely solenoidal forcing, and the force becomes increasingly compressive as ζ decreases from 1 to 0.
3. Simulation parameters
Since we are interested in the local structure of the diffuse ISM, we treat the simulation domain as a part of a larger region with – in the statistical sense – homogeneous properties. Numerically, this is done by using a cubic simulation domain with periodic boundaries and a linear size of X = 40 pc. At the simulation start, the box is filled with atomic hydrogen of uniform density ρ_{0} and temperature . For the mean particle densities n_{0} (which remains constant because of the periodic boundaries), we chose 1.0 cm^{3}, 1.8 cm^{3} and 3.0 cm^{3}, which are comparable to the mean density of the diffuse ISM. In the initial phase, the nearly homogeneous gas will adjust very fast to the equilibrium temperature that is given by the condition ℒ = 0, regardless of the choice of . For this reason, we set equal to the equilibrium temperature corresponding to n_{0}.
Fig. 1 Pressure equilibrium curve for the used cooling and heating processes. Via the isobaric instability criterion (Hunter 1970) the regions of thermal instability (red) and thermal stability (green) are determined. In addition, the isothermal lines at K and K are shown, which define our region of thermally unstable gas. The dots represent the starting points of the simulations, which all lay in the thermally unstable region. 

Open with DEXTER 
The pressure equilibrium curve P_{eq}(n) following from ℒ = 0 and the ideal gas equation (see Sect. 2), is plotted in Fig. 1. Via the isobaric instability criterion of Hunter (1970), one can determine the regions of the thermal instability (TI) in the phase diagram that is shown in Fig. 1. We define three different phases, a cold phase with 200 K, an unstable phase with 200 K 5000 K, and a warm phase with 5000 K. This classification has also been used in other works (Heiles 2001; Heiles & Troland 2003; Audit & Hennebelle 2005) and is warranted by the positions of the isothermal lines at = 200 K, 5000 K, which – of course only very roughly – describe the region of TI in the intermediate and highdensity range. In the following, the terms cold, unstable, and warm phase always refer to the above definitions. The chosen initial conditions of our simulations are located in the thermally unstable region (Fig. 1). In consequence, one can expect that a twophasemedium will evolve. As explained in Sect. 2, we start the simulations with the gas at rest, and turbulence is generated by stochastic forcing. The forcing strength is regulated by the characteristic velocity V. Normalising V by the initial soundspeed , we obtain a characteristic Mach number, Ma = V/c_{0}, which determines the magnitude of the stochastic force. We performed simulations with Ma between 0.2 and 3.0. However, Ma is not an indicator for the average Mach number of the flow, which depends on the phase separation of the gas.
Performed simulations with the mean density n_{0}, the starting temperature , the Mach number Ma, ζ and the simulation time, normalised to the dynamical timescale T.
In one set of simulations, we applied purely compressive forcing (ζ = 0) for different initial densities. In addition, simulations with n_{0} = 1.0 cm^{3} were also carried out for the solenoidal case (ζ = 1). The parameters of all simulations are listed in Table 1. Each simulation ran until an approximate statistically steady state was found for at least one dynamical timescale T. Although it is not obvious that a welldefined statistically steady state can be reached for nonisothermal turbulence, all simulations presented here show a convergence of statistical quantities related to the flow properties and the thermodynamic properties of the gas as time proceeds. Thus, the statistics described in the next sections are averaged over the last dynamical timescale. Since the computational cost of performing simulations of turbulence in nonisothermal gas is significantly higher than for isothermal gas, all runs were performed on a static grid of N = 256^{3} cells. The only exceptions are the runs D1Ma1.0_128 with a resolution of 128^{3} and D1Ma1.0_512 with 512^{3}, to test the resolution dependency of the results. Because of the moderate resolution of the simulations, we only investigated global statistics and probability density functions, but not twopoint statistics or power spectra.
Fig. 2 Time evolution of the temperature field (left panel) and the Mach number (right panel) for the simulation runs D1Ma1.0_128 (upper panel), D1Ma1.0 (mid panel), and D1Ma1.0_512 (lower panel). Plotted are the volumeweighted mean values, the massweighted mean values, and the maximum and minimum values. 

Open with DEXTER 
Fig. 3 Left: temperature PDFs for the runs D1Ma1.0_128, D1Ma1.0 and D1Ma1.0_512, averaged over the last dynamical timescale. The maxima around 55 K and 6500 K indicate the CNM and the WNM. Right: density PDFs for the same runs as in the left panel. For comparison a power law with an exponent of −3.2 is shown. 

Open with DEXTER 
Fig. 4 Temperature (top left), density (top right), vorticity field (bottom left), and Mach number (bottom right) in a slice through the simulation domain of run D1Ma1.0_512 at t = 6.0 T. 

Open with DEXTER 
4. Results
4.1. Basic properties and resolution study
A crucial question for the following analysis is to what extent statistical properties are affected by the numerical resolution. To investigate the dependence on resolution, we consider the simulations D1Ma1.0_128, D1Ma1.0 and D1Ma1.0_512 with compressive forcing, a mean number density n_{0} = 1.0 cm^{3}, and Ma = 1.0. The resolution varies from 128^{3}, over 256^{3} to 512^{3}. Figure 2 shows the time evolution of the temperature field and the Mach number for these three simulations. For them, we see a similar evolution in time both for the temperature field and the velocity field (represented here by the Mach number). The onset of the TI becomes apparent by a rapidly spreading temperature distribution within the first dynamical timescale, where the temperatures finally range from roughly 10 K to several 10 000 K. For the lower resolutions, we did not apply the cutoff at 30 000 K. As one can see, the temperature rises locally to as much as 10^{5} K. For the simulation with 512^{3} grid cells, these extrema entailed problems with the time stepping. Consequently, we introduced the hightemperature cutoff for t > 3 T. After a few dynamical timescales a statistically steady state is approached, in which the mass and volumeweighted mean values are nearly constant. As can be seen in Fig. 2, the massweighted mean of the temperature, , is much lower than the volumeweighted temperature , because the cold gas is associated with highmass densities. The low temperature of the cold, dense gas results in a supersonic massweighted mean Mach number, ℳ_{mw}, while ℳ_{vw} is transonic.
Averaging the different mean values over the last dynamical timescale, ⟨ ⟩ _{t}, we only find small differences between the timeaveraged mean values for the different resolutions. For the massweighted temperature mean values, which correspond to the mean thermal gas energy, we obtain K, 2333 K, and 2326 K for simulation D1Ma1.0_128, D1Ma1.0, and D1Ma1.0_512, respectively. Thus, the relative deviations are less than 5%. For the volumeweighted values , the deviations are even less than 0.6%. For ⟨ ℳ_{vw} ⟩ _{t} and ⟨ ℳ_{vw} ⟩ _{t}, on the other hand, the differences are below 3.7% and 20%, respectively. The larger deviations of the timeaveraged mean Mach numbers point toward a stronger sensitivity of the turbulent velocity field on the numerical resolution, while the thermodynamic properties appear to be quite robust.
Further indicators for the resolution dependency are probability density functions (PDFs). The PDFs of the number density and the temperature, averaged over the last dynamical timescale in the simulations D1Ma1.0_128, D1Ma1.0, and D1Ma1.0_512, are plotted in Fig. 3. The two maxima of the temperature PDFs show that a cold and a warm gas phase are formed, with some gas in the unstable phase in between. The positions of the two maxima around 55 K and 6500 K are almost independent of the numerical resolution. These temperatures can be associated with the CNM and the WNM. Near the extrema of the temperature, there are slight deviations between the PDFs for different resolutions. The small excess of lowtemperature gas in run D1Ma1.0_512 corresponds to a slightly higher amount of very dense gas. The tiny fraction of gas hotter than 10 000 K supports our assumption that the cutoff at 30 000 K in run D1Ma1.0_512 has no significant impact on the dynamics of the whole system. In the temperature interval from roughly 100 K up to 1000 K, we see a decrease in the amount of gas as the resolution increases. Consequently, the phase separation becomes more pronounced at higher resolution. The density PDFs shown in the right panel of Fig. 3 reveal a small excess of highdensity gas for simulation D1Ma1.0_512, which corresponds to the aforementioned increase in cold gas. For the highest resolution, the density PDFs also show a decrease in the amount of gas in the intermediate density range, which corresponds to thermally unstable gas. A notable feature in the highdensity range is the powerlaw shape, which extends over more than one order of magnitude. We estimate the powerlaw exponent of run D1Ma1.0_512 to roughly −3.2, whereas no powerlaw exponents are calculated for the lower resolution runs. Nevertheless, the fit for run D1Ma1.0_512 seems to apply quite well even for run D1Ma1.0_128 and D1Ma1.0, as shown in the right panel of Fig. 3. The power law at these high densities is caused by an effective equation of state (EOS) in this range, which shows a polytropic index γ_{eff} < 1.
Although the statistics we investigated indicate that the simulations D1Ma1.0_128 and D1Ma1.0 are not fully converged with respect to the density and temperature distributions, we conclude that the approximate behaviour of thermally bistable turbulence can be inferred from simulations with 256^{3} grid cells. This is the aim of the parameter study that is presented in the following.
To illustrate the properties of thermally bistable gas, Fig. 4 shows contour plots of the particle density, the temperature, the vorticity norm and the Mach number in 2D slices for simulation D1Ma1.0_512 in the statistically stationary regime. When comparing the slices in the top panel of Fig. 4, one can see that regions of high gas density correspond to low temperatures and vice versa. The most prominent feature that can be seen in the 2D slices is the big clump of cold, dense gas. While the boundaries of this clump are rather sharp, there are also regions in which the cold and warm gas phases entrain each other and form intricate structures. The rich smallscale structure becomes apparent in projections of the mass density and the temperature, which are shown in Fig. 5. Even in these projections, relatively large clumps of cold gas appear prominently. Remarkably, the vorticity is largely reduced in these clumps in comparison to the surrounding medium, which is filled by more or less homogeneous turbulence (see left bottom panel of Fig. 4). Nevertheless, the twodimensional PDF of the vorticity magnitude versus the density, which is shown in Fig. 6, demonstrates that high vorticity is found for the whole range of gas densities. Also the densitydependend mean value shows no clear trend but is more or less constant over the whole range of densities. This seems to contradict the observation made in the vorticity slice in Fig. 4. A possible explanation could be that the vorticity remains low only in larger clumps of overdense gas, while in smaller clumps or in thin filamentary, but nevertheless overdense structures, the vorticity is about as high as in the surrounding lowdensity gas. Presumably, turbulent eddies that are produced by the TI in the unstable phase rapidly decay inside compact cold gas regions, because of the short sound crossing time. The compressive forcing, on the other hand, does not directly generate vorticity and it fails to induce shocks inside the clumps of extremely dense gas. It is possible, however, that these clumps would fragment if the numerical resolution were increased. The Mach number, on the other hand, indicates supersonic motions in the cold, dense gas and subsonic flow in the WNM (see right bottom panel of Fig. 4). This leads to the effect that the mass weighted mean value ℳ_{mw} is higher than the volume weighted mean value ℳ_{vw} as shown in the right panel of Fig. 2.
Slices of the gas density and the vorticity after 1 and 2 integral timescales, which are shown in Fig. 7, demonstrate that turbulence is very inhomogeneous at these stages of the flow evolution. The patches of high vorticity that can be seen are associated with density fluctuations and shocks in the unstable or warm gas phases. Thus, it appears that the initial production of vorticity mainly results from the forcing. As the TI becomes active and temperature fluctuations grow, this production mechanism is accompanied by baroclinic vorticity generation (see, for instance Schmidt 2007). In the late, nearly stationary phase, the high vorticity is found for the whole range of gas densities (see Fig. 6).
4.2. Influence of the forcing strength
So far, we have considered a characteristic Mach number Ma = 1, which means that the velocity of the flow produced by the forcing is close to the speed of sound for the initial gas temperature. Next we analyse the influence of varying the strength of compressive forcing for n_{0} = 1.0 cm^{3} and N = 256^{3}. For the runs D1Ma1.5 and D1Ma2.0, the temperature cutoff at 30 000 K has been applied to avoid artificially high temperatures. In Fig. 8, the temperature and density PDFs are shown for the characteristic Mach numbers in the range from 0.2 to 2.0. In each case, the PDFs are averaged over the last dynamical timescale of the simulation. Except for the temperature PDF of run D1Ma2.0 all other PDFs in the left panel of Fig. 8 show two maxima located in the temperature range of the CNM and the WNM. The positions of the peaks do not vary significantly with the forcing strength in contrast to the peak width. The clearly visible bimodal temperature distribution for the weakly driven simulations D1Ma0.2 and D1Ma0.5 gradually vanishes with increasing forcing strength, while the amount of gas in the unstable range increases. In run D1Ma2.0, no bimodal temperature distribution can be discerned. The density PDFs of the weakly driven simulations with Ma < 1 in the right panel of Fig. 8 also show a bimodal shape, where the strong peak in the lowdensity range (~0.2 cm^{3}) and the plateau around 10 cm^{3} correspond to high and lowtemperature peaks, respectively. With increasing forcing strength, the maximum of the density PDF broadens substantially. Thus, for high characteristic Mach numbers, the density distribution appears to be increasingly dominated by the nearly adiabatic expansion and compression waves produced by the forcing. Moreover, the highdensity tails become flatter; i.e., stronger forcing tends to induce higher densities. However, even for D1Ma2.0, the peak densities up to several 1000 cm^{3} are much higher than what is observed in isothermal simulations with comparable characteristic Mach number (Schmidt et al. 2009). Consequently, the highdensity tails mainly results from the TI (see also Kritsuk & Norman 2002; Li et al. 2003). Accordingly, all density PDFs exhibit a powerlaw like shape in the highdensity range, although the numerical resolution is too low to draw definite conclusions.
Fig. 5 Projected temperature (left) and density (right) fields for run D1Ma1.0_512 at t = 6.0 T. 

Open with DEXTER 
Fig. 6 Twodimensional PDF of the vorticity vs. the density for run D1Ma1.0_512 at t = 6.0 T with logarithmic scaled contour lines. The thick black line shows the densitydependent mean vorticity. 

Open with DEXTER 
The pressuredensity phase diagrams for D1Ma0.2 and D1Ma1.0, which are shown in Fig. 9, demonstrate how the thermodynamic state of the gas is affected by turbulence. In the case of turbulence with low intensity (left panel), the distribution shows typical properties of a static twophasemedium as described by Field et al. (1969). Here the gas exists in a rough pressure equilibrium with deviations from the mean pressure less than half an order of magnitude. Only the highdensity gas is out of pressure equilibrium, because of the very small cooling times corresponding to these densities. For stronger turbulent motions, as in run D1Ma1.0 (right panel of Fig. 9), the gas deviates largely from pressure equilibrium. Thus, for gas in a strongly turbulent state, the classical idea of two distinct phases, coexisting in a rough pressure equilibrium with only small amounts of gas in the unstable range ceases to apply.
To quantify the increasing amount of thermally unstable gas that is observed for strong forcing amplitudes, we determined the mass fractions of the warm, unstable and cold phases as defined in Sect. 3. The results are shown in Fig. 10, where the different mass fractions are plotted against the total specific kinetic energy for each simulation. In accordance with the PDFs, the mass fractions, as well as the specific kinetic energies, are averaged over the last dynamical timescale. Clearly, the turbulence energy resulting from the forcing has a strong influence on the distribution of the gas into the different phases. The mass fraction of the unstable gas increases at the cost of the mass fraction of the cold gas, whereas the mass fraction of the warm phase is only weakly affected. That the redistribution of the gas with increasing forcing strength mainly occurs between the cold and the unstable phase and, to a lesser extent, between the warm and the unstable phase can be understood as the consequence of the low compressibility of the hot gas, in which only sub or transonic motions are induced, regardless of the characteristic Mach number. The dynamics in the unstable phase and the cold phase, on the other hand, are highly sensitive to the forcing strength. Although compressive forcing produces the seeds for the TI, we note that high turbulence intensity inhibits the gas from settling into the cold phase. However, it appears that the TI overcomes this effect once a critical amount of cold gas is accumulated. This is suggested by the occurrence of compact regions of cold gas with weakly turbulent interiors (see Fig. 4).
Fig. 7 Slices of the density (left) and the vorticity field (right) at time t = 1.0 T (top) and 2.0 T (bottom) for run D1Ma1.0_512. 

Open with DEXTER 
Fig. 8 Left: temperature PDFs of simulations with n_{0} = 1.0 cm^{3} and N = 256^{3}. Except for run D1Ma2.0, the PDFs show two maxima in the warm and cold phase. Right: density PDFs for the same runs as in the left panel. The PDFs show a peak around n = 0.2 cm^{3} and a plateau around 10 cm^{3}, which correspond to the WNM and the CNM, respectively. 

Open with DEXTER 
Fig. 9 Pressuredensity phase diagrams of run D1Ma0.2 (left) and D1Ma1.0 (right) with a logarithmic scaling of the contour lines. The equilibrium between heating and cooling is indicated by the solid black line. 

Open with DEXTER 
4.3. Tracer particle statistics
To examine the cause of the higher mass fraction of gas in the thermally unstable range for stronger forcing, we included 150 000 Lagrangian tracer particles in the simulations D1Ma0.2 and D1Ma1.0. This approach is motivated by Audit & Hennebelle (2005), who also found an increase in the thermally unstable gas with increasing turbulence intensity at the cost of the cold gas mass. The authors explain this effect on the basis of a simple semianalytical model. This model implies that the time interval that is spent by a fluid element in the unstable phase before returning to the stable phase becomes longer as the turbulence intensity increases. On the other hand, the mass of gas in the unstable phase is also controlled by the frequency of strong perturbations that bring the gas out of the equilibrium state. Thus, the larger amount of gas in the unstable phase can also be interpreted as a consequence of enhanced turbulent mixing. With the help of the tracer particles, we are able to test the predictions of this model by explicitly following the time evolution of a fluid element moving with the flow. In particular, we measure the timeinterval t_{UP,i}, defined as the time, a tracer particle is located in the unstable range before returning to the stable state of the CNM or the WNM. We find that the mean duration of the stay of a tracer particle in the unstable phase is reduced by a factor of 2.3 if the characteristic Mach number of the forcing increases from 0.2 to 1.0. However, after analysing how often the particle is perturbed out of the equilibrium state into the unstable phase, we find that the frequency of the perturbations, , rises by a factor of ~4.5 from run D1Ma0.2 to D1Ma1.0. This implies that the thermally unstable gas fraction roughly increases from run D1Ma0.2 to D1Ma1.0 by a factor of .
Fig. 10 Timeaveraged mass fractions of the warm (red points), unstable (green points), and cold phase (blue points) for simulations with n_{0} = 1.0 cm^{3} and N = 256^{3} plotted against the specific kinetic energy. An increasing kinetic energy corresponds to an increase in the forcing Mach number (compare Table 1). 

Open with DEXTER 
4.4. Variation in the mean density
It is also of interest to examine how the gas distribution changes if the mean density of the gas is altered. Interestingly, it turns out that even a small increase of the mean density has a noticeable effect on the distribution of the gas among the different phases. In particular, we performed several simulations with compressive forcing for mean densities of 1.8 cm^{3} and 3.0 cm^{3} (see Table 1). Figure 1 shows that the initial state of the gas in these simulations is located closer to the local minimum of the equilibrium curve. For n_{0} = 3.0 cm^{3}, the initial state is located near the lowtemperature boundary of the central unstable region. The dependence of the mass fractions of the different gas phases on the forcing strength is shown in Fig. 11. While the qualitative behaviour is similar to what is observed for n_{0} = 1.0 cm^{3} (see Fig. 10), the fraction of cold gas is clearly greater for a higher mean density. In the case n_{0} = 3.0 cm^{3}, about 90% of the mass settles into the cold phase, more or less independent of the turbulence energy. This result is reasonable because, as the initial state is chosen closer to the cold phase, an increasing amount of the gas will cool and settle into the cold phase due to the external compressions. At densities n_{0} ≳ 10 cm^{3}, most of the gas resides in the nearly isothermal right branch of the equilibrium curve, and only strong rarefactions can push the gas into the warm phase. In this regime, the dynamical properties are comparable to those of isothermal turbulence.
Fig. 11 Left: time averaged mass fraction of the warm (red points), unstable (green points) and cold phase (blue points) for the simulations with n_{0} = 1.8 cm^{3} plotted against the specific kinetic energy. Right: same as in the left panel but for simulations with n_{0} = 3.0 cm^{3}. 

Open with DEXTER 
Fig. 12 Left: density PDFs of the solenoidally driven runs with mean densities of n_{0} = 1.0 cm^{3}. Right: temperature PDFs for the same runs as in the left panel. 

Open with DEXTER 
4.5. Solenoidal forcing
Schmidt et al. (2009) and Federrath et al. (2010) show in their analysis of supersonic isothermal turbulence that there are significant differences between compressive and solenoidal forcing. In particular, the density PDFs are much broader for compressively driven turbulence (Federrath et al. 2008, 2010). For nonisothermal gas, the differences should be even more pronounced because the TI is affected by the different density distributions, which are generated by distinct forms of the forcing. The PDFs of the number density and the temperature resulting from three simulations with solenoidal forcing and a mean density of 1.0 cm^{3} (listed in Table 1) are plotted in Fig. 12, where the temperature cutoff at 30 000 K was applied for the runs D1Ma1.0s and D1Ma2.0s. As one can see, there is no phase separation at all. Both the density and temperature PDFs show a single peak that roughly corresponds to the initial state. With increasing forcing strength, the distributions become broader, and the gas tends to become hotter. The higher temperatures that are found for Ma = 2.0 in comparison to the other runs are caused by the stronger dissipative heating that balances the higher rate of energy injection by the forcing. The amount of cold gas with temperatures below 200 K is zero. Thus, it appears that the TI is largely suppressed by solenoidal forcing, although a significant fraction of the gas is in the unstable region (see Fig. 1). The density PDFs are similar to those reported by Passot & VázquezSemadeni (1998) for simulations of turbulence in a polytropic gas with γ > 1, particularly, with regard to the lowdensity tails. At least for the runs with Ma ≤ 1.0, there appears to be a powerlaw range for densities lower than 1.0 cm^{3}.
5. Discussion
Most simulations considered in this article were run with a resolution of N = 256^{3} grid cells. For the chosen box size of 40 pc, the cutoff length in these simulations is Δ = 0.16 pc. Consequently, cold gas structures with a typical size in the range of a few parsecs are only marginally resolved and substantially smeared out by numerical diffusion. For this reason, it cannot be expected that the properties of these structures can be described in a detailed quantitative manner. However, the resolution study we performed for a particular parameter set suggests that the turbulent phase separation of the gas can be approximated reasonably well. In particular, we showed that trends can be seen if the forcing is varied, which allows for a qualitative discussion. This will help us selecting parameters for more elaborate numerical studies at higher resolution, using adaptive mesh refinement. In particular, an enhanced numerical resolution will be necessary to confirm the indications for powerlaw tails of the mass density PDFs toward high densities and to determine the slope with higher precision. These power laws can be found in nearly all simulations of compressively driven turbulence in the present study. This clearly supports the results of Scalo et al. (1998), Passot & VázquezSemadeni (1998), and Li et al. (2003), who applied a simple polytropic EOS with γ < 1 as a simple model for cooling. These authors obtained skewed lognormal density PDFs with powerlaw tails at high densities. In our simulations, the thermodynamic state of the gas at high densities roughly follows the equilibrium curve, which corresponds to an effective polytropic index γ_{eff} < 1. Powerlaw tails for nonisothermal turbulence followed also from the simulations of decaying turbulence by Kritsuk & Norman (2002) and from a recent collidingflow simulation of Audit & Hennebelle (2010), in which a polytrope with γ = 0.7 was applied.
The main result of our study is the increase in the mass of the thermally unstable gas for higher turbulence energy (Fig. 10), which is roughly balanced by a decrease in the mass of cold gas. Consequently, stronger compressive forcing does not result in more cold, dense gas. The reason is found in the broadening of the pressure distribution with increasing forcing strength and, thus, kinetic energy of the turbulent flow. While weak forcing produces the initial seeds for the TI by perturbing the density of the gas, the forcing drives the gas away from a state of approximate pressure equilibrium if the turbulent velocity becomes comparable to the speed of sound. This effect is shown clearly for runs D1Ma0.2 and D1Ma1.0 in Fig. 9. Only at high densities is the gas mainly concentrated around the equilibrium curve, independent of forcing strength. This is a consequence of the very short cooling time in comparison to the dynamical timescale and the sound crossing time in this range. At low densities, on the other hand, the gas in run D1Ma1.0 is nearly adiabatic because of the long cooling time in the WNM. The effective polytropic index in this regime is found to be γ_{eff} = 1.52. As one can see in Fig. 8 (right panel), this entails a pronounced broadening of the density PDFs in the underdense range if Ma ≳ 1, which is again similar to the results of Passot & VázquezSemadeni (1998) for their simulations with γ > 1.
In Sect. 4.3, we show with the help of tracer particles that the higher frequency of the perturbations of the particles out of the stable regime, , overcompensates the decrease in the mean duration, , a tracer particle stays in the unstable phase if the forcing strength increases. This results in an increase in the unstable mass fraction roughly by a factor of 2.0 for simulation D1Ma1.0 compared to D1Ma0.2, which approximates the ratio of 2.7 that follows from the mass fractions of unstable gas in these simulations reasonably well. Our observation that decreases with increasing forcing strength is contrary to the assumption of Audit & Hennebelle (2005) that increases with the intensity of turbulence. The discrepancy could appear because their model is 2dimensional and applies to the unperturbed transition of a fluid element from the WNM to the CNM, which might be insufficient to describe 3dimensional turbulence with external forcing. In order to rule out that the discrepancy with the model of Audit & Hennebelle (2005) arises because it only describes transitions from the WNM to the CNM, we also examined phase transitions between the WNM and the CNM with the tracer particles. We found that the mean transition time decreases from simulations D1Ma0.2 to D1Ma1.0 by a factor of ~2.5, which is similar to the behaviour of . The averaged time intervals and , a fluid particle stays in the CNM or in the WNM, become significantly shorter as well. Summarising, the tracer particle statistics imply that phase transitions between the different phases occur more frequently with increasing forcing strength, i.e., turbulent mixing is more effective, which results in a higher mass fraction of the unstable phase.
An important question is for which simulation parameters the mass fractions of the different phases that are plotted in Figs. 10 and 11 are comparable to observational results for the atomic ISM (see, for instance, Heiles 2001; Heiles & Troland 2003). One difficulty in such a comparison is that observational values usually do not refer to a particular mean density of the gas, for which measurements are taken. Moreover, the observations presented by these authors do not belong to certain diffuse H iclouds but are absorption measurements on cold gas along the lines of sight against different background radio sources. Besides the mean Mach number in the cold phase, which the authors determined to be around 3, particular properties of these objects are not known. Nevertheless, we note that Heiles & Troland (2003) find that the average mass fraction of atomic gas in the ISM with temperatures in the range from 500 K to 5000 K is at least 29%. Together with the warm gas fraction, their measurements imply that the mass fraction of cold gas should be less than 40%, which is clearly at the lower end of the range of values following from our numerical study. Thus, strong compressive forcing of turbulence in gas at a mean density of 1.0 cm^{3} appears to be typical of the conditions in the diffusive ISM. However, the sources studied by Heiles & Troland (2003) show a large scatter around the average values, and there exist regions with mass fractions of cold gas up to nearly 90%, which would correspond to weak forcing and a higher mean density in our study.
A surprising result is that we did not find any cold phase at all in the simulations with purely solenoidal forcing, regardless of the magnitude of the force. Although the turbulent velocity fluctuations in compressible gas entail density fluctuations and the initial state of the gas is in the unstable regime, it appears that the phase separation is completely suppressed. One possible explanation is that the TI would set in on length scales that are smaller than the numerically resolved scales in our simulations. The minimum length scale of the TI is the Field length (Field 1965), which is defined by (7)where K is the thermal conductivity. According to (Cox 2000), (8)With this, we get a Field length in the CNM that is at least 2 or 3 orders of magnitude smaller than the spatial grid resolution. Therefore, a cold phase might develop at higher resolution. Hennebelle & Audit (2007) show in a 2dimensional resolution study that a spatial resolution below 2 × 10^{3} pc is necessary to reach numerical convergence for a colliding flow of thermally bistable gas, which is unfeasible for 3dimensional simulations with the current computing power. However, the observed behaviour could also be physical and result from the apparent differences between the turbulent velocity field produced by solenoidal forcing compared to compressively driven turbulence (Schmidt et al. 2008b). There are indications of an extended range of length scales, for which rotationally driven turbulence has not only lower density fluctuations but also exhibits a more spacefilling vorticity field (Federrath et al. 2008, 2010). As a consequence, the turbulent pressure could stabilise the gas effectively against the TI. An analogous effect has been observed for forced selfgravitating turbulence in an isothermal gas (e.g. Klessen et al. 2000; Schmidt et al. 2008a). In addition, it is likely that the shearing effect of the velocity fluctuations produced by solenoidal forcing disrupts incipient small condensations of dense gas. On the other hand, Gazol et al. (2005, 2009) and Kissmann et al. (2008) found a cold phase for turbulence in nonisothermal gas that was produced by solenoidal forcing. These findings do not necessarily contradict our results, because Kissmann et al. (2008) used a mean density of 2.0 cm^{3}, and Gazol et al. (2005, 2009) applied a slightly different cooling function. Thus, the initial state was located deeper inside the unstable regime in both cases. For this reason, the TI can overcome the inhibiting effect of turbulence more easily than in our simulations with n_{0} = 1.0 cm^{3}. Indeed, we also found a cold phase developing in simulations with solenoidal forcing if the mean density was set to 3.0 cm^{3}. Taking these results together, it appears that the evolution of thermally unstable gas is extremely sensitive to the initial conditions and the details of the cooling function if the forcing is solenoidal.
The strong sensitivity on the mean density and on the forcing shown here could in principle be used to draw some conclusions about the forcing modes in observed diffuse H iclouds. In practice, however, it is hard to state something about the underlying forcing because of the difficulty determining the mean gas density observationally, and hence deciding whether the gas is thermally unstable or not. Apart from this difficulty, this work suggests that, for a marginally thermally unstable cloud (n_{0} ≃ 1 cm^{3}) showing a bimodal temperature distribution, the forcing should mainly be compressive. On the other hand, for clouds in the highly unstable regime (n_{0} ≳ 3 cm^{3}), it is hardly possible to distinguish between solenoidal or compressive forcing as both forcing modes would produce a twophasemedium with gas distributions that are mostly determined by the thermal instability.
6. Conclusion
In this paper, we have performed a parameter study of forced turbulence in thermally bistable gas, with a cooling function that was defined by Audit & Hennebelle (2005). Considering both compressive and solenoidal forcing, we find substantial variations in the global statistics and probability density functions of the gas. Because of the limited resolution of the simulations, we did not compute twopoint statistics or power spectra, nor did we analyse fragmentation properties. In the following, we summarise our results.

1.
Although we clearly see a phase separation for compressively driven turbulence already at moderate numerical resolution, the PDFs of the mass density and the temperature suggest that a resolution higher than 512^{3} is required to accurately determine the tails of the distributions and the mass fractions of the different phases. For a forcing length scale around 20 pc, a spatial resolution significantly lower than 0.1 pc is implied.

2.
A similarity of the density PDFs of all compressively driven simulations is the powerlaw like shape in the highdensity range, which is caused by an effective polytropic index γ_{eff} < 1. Exploring the highdensity tails in simulations with enhanced numerical resolution will bear consequences on theoretical predictions for the core mass function (Hennebelle & Chabrier 2009).

3.
For weak forcing, the warm and the unstable phase are nearly isobaric. With increasing forcing strength, we find larger deviations from the pressure equilibrium and a greater amount of gas with temperatures in the intermediate, unstable range, which is roughly balanced by less gas in the cold phase. The Lagrangian statistics inferred from tracer particle suggest that this behaviour is caused by stronger turbulent mixing.

4.
The distribution of the gas among the different phases also depends on the mean density. For mean number densities higher than 1.0 cm^{3}, most of the gas tends to settle into the cold phase.

5.
If the forcing is purely solenoidal, we do not find bimodal density and temperature distributions for a mean density of 1.0 cm^{3}. In this case, all the gas is in the thermally unstable and warm phase.
Qualitatively, our results agree with other numerical studies of turbulence in nonisothermal gas (e.g. Kritsuk & Norman 2002; Gazol et al. 2005; Audit & Hennebelle 2005). The largescale flow structure produced by compressive stochastic forcing resembles the colliding flow scenario, since transient frontal collisions of shocks occur randomly (Schmidt et al. 2009). Given the sensitivity of the phase separation on the mean density, a crucial difference arises from the accumulation of mass in the colliding flow scenario, because of the inflow boundary conditions. In a periodic box, on the other hand, the mass is conserved. Consequently, in colliding flow simulations, there is an implicit drift through the parameter space in the course of time. In this regard, we note that Audit & Hennebelle (2010) do not find clear indications of a powerlaw tail of the density PDF in their simulations with cooling. If selfgravity is included, a caveat of periodicbox simulations is that the amount of mass that is confined in gravitationally collapsing objects cannot indefinitely increase at a constant rate. A net inflow of gas from the boundaries, on the other hand, constantly rejuvenates the supply of gas that can be bound by gravitational collapse (e.g. VázquezSemadeni et al. 2007; Hennebelle et al. 2008; Banerjee et al. 2009), although a statistically stationary state cannot be reached either because of the reasons mentioned above. Whether the physical conditions in the interstellar medium on a length scale of a few 10 pc are approximated more closely in simulations of colliding flows without external forcing compared to forced turbulence simulations at a fixed mean density remains a matter of debate. Regardless of this, the results presented here show a substantial sensitivity on the initial state of the gas and the forcing properties. This illustrates that a restriction to highly specialised conditions possibly leads to misleading conclusions, as even minor changes in the studied conditions can cause significantly different results.
In largescale simulations of the ISM (e.g. Slyz et al. 2005; Joung & Mac Low 2006; Joung et al. 2009), the spatial resolutions are usually above 1 pc, and in simulations of galactic disks, the resolution is even coarser (e.g. Agertz et al. 2009; Tasker & Tan 2009). Thus, as shown here, it is quite likely that the mass fractions of the gas in the different phases would change significantly if smaller length scales could be resolved, particularly in the cold phase and in the lowtemperature segment of the unstable phase. This has a potential impact on calculations of the star formation rate, because only cold gas can turn into stars. A subgridmodel for star formation that takes the multiphase structure of the ISM into account has already been developed and applied by Springel & Hernquist (2003), although these authors omitted the influence of turbulent motions that we focussed on. However, our results indicate that the errors due to unresolved thermodynamic processes will be small if the gas density in the grid cells where sink particles representing starforming clouds are added is much higher than 1.0 cm^{3}. Only in regions that just begin to contract, a correction of the star formation efficiency by the mass fraction of cold gas might slow down the growth of the stellar mass in diskgalaxy simulations. In the dense clumps that are formed subsequently, however, most of the gas would cool down to molecularcloud temperatures if a sufficiently low temperature floor was applied, and it can be assumed that the local star formation efficiency is regulated mainly by gravity and smallscale turbulence (Krumholz & McKee 2005; Padoan & Nordlund 2009; Klessen & Hennebelle 2010).
For further studies of turbulence in the multiphase ISM, adaptive mesh refinement (AMR) is one method to increase the numerical resolution in such simulations. Although Niklaus et al. (2009) show that there are serious difficulties in reproducing the properties of 2dimensional colliding flows by means of AMR even if sophisticated refinement criteria are applied, we expect that the situation will be more favourable to 3dimensional simulations of forced turbulence. If the low enstrophy (i.e., the volumeintegrated square of the vorticity) of compact cold gas regions could be confirmed in AMR simulations, it would point to an additional driving mechanism to sustain supersonic turbulence in molecular cloud interiors. Examples include external driving by gravitational accretion (e.g. Klessen & Hennebelle 2010) or internal driving by massive stars (e.g. Fall et al. 2010). Moreover, there is also room for methodological improvements in the cooling. Firstly, the range of the cooling function should be extended to temperatures higher than 10 000 K. On the other hand, the cooling of the gas at low temperatures will be altered if chemical processes, in particular, the generation and destruction of molecular hydrogen, are explicitly included (Glover et al. 2010). Secondly, it is not clear whether an explicit scheme with subcycling will be sufficient to treat the source terms in the energy equation if the effective numerical resolution becomes very high. For subsequent simulations it would also be interesting to study the effect of magnetic fields.
Acknowledgments
We thank Edouard Audit and Patrick Hennebelle for providing their cooling function. We also thank Christoph Federrath for providing his tracer particle algorithm. The numerical simulations presented in this article were performed on HLRB2 at the Leibniz Supercomputing Centre in Garching, as part of project h0972 using the Enzo code developed by the Laboratory for Computational Astrophysics at the University of California in San Diego (http://lca.ucsd.edu). D. Seifried is currently founded by the DFG via the EmmyNoether grant BA3706.
References
 Agertz, O., Lake, G., Teyssier, R., et al. 2009, MNRAS, 392, 294 [NASA ADS] [CrossRef] [Google Scholar]
 Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Audit, E., & Hennebelle, P. 2010, A&A, 511, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef] [Google Scholar]
 BallesterosParedes, J., Klessen, R. S., Mac Low, M.M., & VazquezSemadeni, E. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 63 [Google Scholar]
 Banerjee, R., VázquezSemadeni, E., Hennebelle, P., & Klessen, R. S. 2009, MNRAS, 398, 1082 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cox, A. N. 2000, Allen’s Astrophysical Quantities, 4th edn. (SpringerVerlag) [Google Scholar]
 Dickey, J. M., Salpeter, E. E., & Terzian, Y. 1977, ApJ, 211, L77 [NASA ADS] [CrossRef] [Google Scholar]
 Dobbs, C. L., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2008, MNRAS, 389, 1097 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Eswaran, V., & Pope, S. B. 1988, Comput. Fluids, 16, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Fall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, ApJ, 710, L142 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Klessen, R. S., & Schmidt, W. 2009, ApJ, 692, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Duval, J., Klessen, R., Schmidt, W., & Mac Low, M. M. 2010, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Field, G. B. 1965, ApJ, 142, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJ, 155, L149 [NASA ADS] [CrossRef] [Google Scholar]
 Folini, D., Walder, R., & Favre, J. M. 2010, in Numerical Modeling of Space Plasma Flows, ASP Conf. Series, 9 [Google Scholar]
 Gazol, A., VázquezSemadeni, E., & Kim, J. 2005, ApJ, 630, 911 [NASA ADS] [CrossRef] [Google Scholar]
 Gazol, A., Luis, L., & Kim, J. 2009, ApJ, 693, 656 [NASA ADS] [CrossRef] [Google Scholar]
 Glover, S. C. O., Federrath, C., Mac Low, M., & Klessen, R. S. 2010, MNRAS, 404, 2 [NASA ADS] [Google Scholar]
 Heiles, C. 2001, ApJ, 551, L105 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C., & Troland, T. H. 2003, ApJ, 586, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Audit, E. 2007, A&A, 465, 431 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Chabrier, G. 2009, ApJ, 702, 1428 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., Banerjee, R., VázquezSemadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hunter, Jr., J. H., 1970, ApJ, 161, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Joung, M. K. R., & Mac Low, M.M. 2006, ApJ, 653, 1266 [NASA ADS] [CrossRef] [Google Scholar]
 Joung, M. R., Mac Low, M.M., & Bryan, G. L. 2009, ApJ, 704, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Kanekar, N., Subrahmanyan, R., Chengalur, J. N., & Safouris, V. 2003, MNRAS, 346, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Kissmann, R., Kleimann, J., Fichtner, H., & Grauer, R. 2008, MNRAS, 391, 1577 [NASA ADS] [CrossRef] [Google Scholar]
 Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klessen, R. S., Heitsch, F., & Mac Low, M. 2000, ApJ, 535, 887 [NASA ADS] [CrossRef] [Google Scholar]
 Kritsuk, A. G., & Norman, M. L. 2002, ApJ, 569, L127 [NASA ADS] [CrossRef] [Google Scholar]
 Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Y., Klessen, R. S., & Mac Low, M.M. 2003, ApJ, 592, 975 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Niklaus, M., Schmidt, W., & Niemeyer, J. C. 2009, A&A, 506, 1065 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 O’Shea, B. W., Bryan, G., Bordner, J., et al. 2005, in Adaptive Mesh Refinement – Theory and Applications, ed. T. Plewa, T. Linde, & V. G. Weirs, Springer Lecture Notes in Computational Science and Engineering, 41 [Google Scholar]
 Padoan, P., & Nordlund, A. 2009, ApJ, accepted [arXiv:0907.0248] [Google Scholar]
 Passot, T., & VázquezSemadeni, E. 1998, Phys. Rev. E, 58, 4501 [NASA ADS] [CrossRef] [Google Scholar]
 Scalo, J., VazquezSemadeni, E., Chappell, D., & Passot, T. 1998, ApJ, 504, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, W. 2007, in Structure Formation in the Universe, ed. G. Chabrier, Cambridge Contemporary Astrophysics (Cambridge University Press), 20 [Google Scholar]
 Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006, Comput. Fluids, 35, 353 [CrossRef] [Google Scholar]
 Schmidt, W., Federrath, C., Hupp, M., Maier, A., & Niemeyer, J. C. 2008a, in High Performance Computing in Science and Engineering, Garching/Munich 2007, ed. S. Wagner, M. Steinmetz, A. Bode, & M. Brehm (Springer) [Google Scholar]
 Schmidt, W., Federrath, C., & Klessen, R. 2008b, Phys. Rev. Lett., 101, 194505 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Slyz, A. D., Devriendt, J. E. G., Bryan, G., & Silk, J. 2005, MNRAS, 356, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. 1978, Physical Processes in the Interstellar Medium, 1st edn (WileyInterscience) [Google Scholar]
 Spitzer, Jr., L., & Fitzpatrick, E. L. 1995, ApJ, 445, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, ApJ, 508, L99 [NASA ADS] [CrossRef] [Google Scholar]
 Tasker, E. J., & Bryan, G. L. 2006, ApJ, 641, 878 [NASA ADS] [CrossRef] [Google Scholar]
 Tasker, E. J., & Tan, J. C. 2009, ApJ, 700, 358 [NASA ADS] [CrossRef] [Google Scholar]
 VázquezSemadeni, E., Gazol, A., Passot, T., et al. 2003, in Turbulence and Magnetic Fields in Astrophysics, ed. E. Falgarone, & T. Passot, Lecture Notes in Physics (Berlin: Springer Verlag), 614, 213 [Google Scholar]
 VázquezSemadeni, E., Gómez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Performed simulations with the mean density n_{0}, the starting temperature , the Mach number Ma, ζ and the simulation time, normalised to the dynamical timescale T.
All Figures
Fig. 1 Pressure equilibrium curve for the used cooling and heating processes. Via the isobaric instability criterion (Hunter 1970) the regions of thermal instability (red) and thermal stability (green) are determined. In addition, the isothermal lines at K and K are shown, which define our region of thermally unstable gas. The dots represent the starting points of the simulations, which all lay in the thermally unstable region. 

Open with DEXTER  
In the text 
Fig. 2 Time evolution of the temperature field (left panel) and the Mach number (right panel) for the simulation runs D1Ma1.0_128 (upper panel), D1Ma1.0 (mid panel), and D1Ma1.0_512 (lower panel). Plotted are the volumeweighted mean values, the massweighted mean values, and the maximum and minimum values. 

Open with DEXTER  
In the text 
Fig. 3 Left: temperature PDFs for the runs D1Ma1.0_128, D1Ma1.0 and D1Ma1.0_512, averaged over the last dynamical timescale. The maxima around 55 K and 6500 K indicate the CNM and the WNM. Right: density PDFs for the same runs as in the left panel. For comparison a power law with an exponent of −3.2 is shown. 

Open with DEXTER  
In the text 
Fig. 4 Temperature (top left), density (top right), vorticity field (bottom left), and Mach number (bottom right) in a slice through the simulation domain of run D1Ma1.0_512 at t = 6.0 T. 

Open with DEXTER  
In the text 
Fig. 5 Projected temperature (left) and density (right) fields for run D1Ma1.0_512 at t = 6.0 T. 

Open with DEXTER  
In the text 
Fig. 6 Twodimensional PDF of the vorticity vs. the density for run D1Ma1.0_512 at t = 6.0 T with logarithmic scaled contour lines. The thick black line shows the densitydependent mean vorticity. 

Open with DEXTER  
In the text 
Fig. 7 Slices of the density (left) and the vorticity field (right) at time t = 1.0 T (top) and 2.0 T (bottom) for run D1Ma1.0_512. 

Open with DEXTER  
In the text 
Fig. 8 Left: temperature PDFs of simulations with n_{0} = 1.0 cm^{3} and N = 256^{3}. Except for run D1Ma2.0, the PDFs show two maxima in the warm and cold phase. Right: density PDFs for the same runs as in the left panel. The PDFs show a peak around n = 0.2 cm^{3} and a plateau around 10 cm^{3}, which correspond to the WNM and the CNM, respectively. 

Open with DEXTER  
In the text 
Fig. 9 Pressuredensity phase diagrams of run D1Ma0.2 (left) and D1Ma1.0 (right) with a logarithmic scaling of the contour lines. The equilibrium between heating and cooling is indicated by the solid black line. 

Open with DEXTER  
In the text 
Fig. 10 Timeaveraged mass fractions of the warm (red points), unstable (green points), and cold phase (blue points) for simulations with n_{0} = 1.0 cm^{3} and N = 256^{3} plotted against the specific kinetic energy. An increasing kinetic energy corresponds to an increase in the forcing Mach number (compare Table 1). 

Open with DEXTER  
In the text 
Fig. 11 Left: time averaged mass fraction of the warm (red points), unstable (green points) and cold phase (blue points) for the simulations with n_{0} = 1.8 cm^{3} plotted against the specific kinetic energy. Right: same as in the left panel but for simulations with n_{0} = 3.0 cm^{3}. 

Open with DEXTER  
In the text 
Fig. 12 Left: density PDFs of the solenoidally driven runs with mean densities of n_{0} = 1.0 cm^{3}. Right: temperature PDFs for the same runs as in the left panel. 

Open with DEXTER  
In the text 