text:
# The nonlinear development of the thermal instability in the atomic ISM and its interaction with random fluctuations

## 1 Introduction

The fact that the neutral atomic interstellar medium (ISM) is most likely thermally bistable (Zeldovich & Pikel’ner 1969; Field, Goldsmith & Habing 1969; Wolfire et al. 1995) has had a great impact on our picture of the formation of interstellar structure. In particular, the unstable cooling process has been thought to be responsible for the collection of mass into HI clouds of sizes $`1`$ pc, densities $`10`$–100 cm-3, and temperatures $`100`$ K (e.g., Goldsmith 1970; Schwarz, McCray & Stein 1972; Parravano 1987; Lioure & Chièze 1990; Koyama & Inutsuka 2002). In the classical model of McKee & Ostriker (1977), the cold and warm neutral media are two distinct phases, with clouds forming in dense shells created by shock compressions, aided by the thermal instability (TI hereafter), and later confined by the warm medium’s thermal pressure. The triggering of TI in the warm medium by strong compressions and the subsequent segregation in two separate phases has been recently studied by Hennebelle & Pérault (1999, 2000) and by Koyama & Inutsuka (2000).

However, over the years, a more complex picture of the ISM has emerged, driven in part by the advent of realistic numerical simulations of the ISM that have highlighted the role of inertia and advection (i.e., transport) by the velocity field. In this context, the role of the thermal instability in the overall dynamics of the ISM is still uncertain. Simulations of the turbulent ISM driven by supernovae or stellar winds suggest that the complex structure of the ISM as well as HI clouds can be the result of random gas motions, at least in regions with densities below $`10`$$`100`$ cm-3, with no need for thermal instabilities (e.g. Bania & Lyon 1980; Rosen & Bregman 1995; Passot, Vázquez-Semadeni & Pouquet 1995; Korpi et al.111The simulations of Korpi et al. actually use a cooling function that imply a thermally unstable range between the cold and warm gas if no heating is present (see, e.g., Murray & Lin 1989; Burkert & Lin 2000). Although no global heating was used in those simulations, the presence of local stellar heating causes the gas to be thermally stable while subject to heating, and unstable while cooling. 1999; see also Vázquez-Semadeni & Passot 1999 and Vázquez-Semadeni 2002 for reviews; but see Wada, Spaans & Kim 2000; Koyama & Inutsuka 2002; Kritsuk & Norman 2002 for the alternative view that TI powers ISM turbulence).

Within this scenario, Vázquez-Semadeni, Gazol & Scalo (2000, hereafter Paper I) studied the probability density distribution (PDF) of mass density in low resolution ($`128^2`$) 2D simulations of the thermally bistable ISM at large scales (1-kpc box), finding that the PDF may not reflect the thermally unstable nature of the gas in the presence of magnetic fields, Galactic shear, the Coriolis force, and turbulence driven by stellar-like sources located at the density maxima (clouds). Moreover, simulations at higher resolution ($`512^2`$) with both the same box size (1 kpc; Gazol et al. 2001, hereafter Paper II), and smaller ones (150 and 10 pc; Vázquez-Semadeni et al. 2002, hereafter Paper III) have temperature histograms that indicate the existence of significant amounts ($`50\%`$) of gas with temperatures in the unstable range. The same result could be seen in the simulations by Korpi et al. (1999), although they did not discuss it, and has later also been found in simulations of the late nonlinear stages of the pure instability at smaller scales (6-pc box, at $`256^3`$ resolution) by Kritsuk & Norman (2002). This is in agreement with observations by Dickey, Salpeter & Terzian (1978), Kalberla, Schwartz, & Goss (1985), Spitzer & Fitzpatrick (1995), Fitzpatrick & Spitzer (1997), and Heiles (2001), as well as previous theoretical suggestions (e.g., Lioure & Chièze 1990; Norman & Ferrara 1996). Thus, the study of TI in the context of fluctuating media becomes important, in particular to understand the origin of the large fractions of gas mass in the unstable regime.

At the large scales studied by Paper I, the condensation222For consistency with the nomenclature in the literature, we maintain the name “condensation” even for highly dynamic cases in which “accretion” or “collapse” might be more appropriate. process leading to the formation of clouds and intercloud “voids” becomes highly dynamic (i.e., develops large velocities in the nonlinear stage). The parameter controlling the degree to which the condensation process is dynamic is the ratio ($`\eta `$) of the cooling to the representative (sound or bulk) dynamical time. Although it is customary in the literature to discuss TI in terms of the perturbation wave lengths (or wavenumbers; e.g. Field 1965), in this paper we prefer a description in terms of $`\eta `$ because it allows the consideration of perturbations of a given physical size (say, in parsecs), but that may correspond to different values of $`\eta `$ because of, say, they occur at different mean densities (i.e., equilibrium cooling rates).

It is well known (e.g., Field 1965; Murray & Lin 1989; Burkert & Lin 2000) that perturbations with $`\eta >1`$ (corresponding to small-scale perturbations at a fixed cooling rate or low cooling rate at a fixed length scale) evolve almost isobarically, and thus quasi-statically. Instead, perturbations with $`\eta <1`$ develop large pressure gradients, and thus evolve highly dynamically even under the isobaric mode of the instability, although they take longer times to grow (e.g., Field 1965; Balbus 1995; Meerson 1996). Thus, Paper I suggested that large clouds formed by TI should be bounded by accretion shocks (rather than by contact discontinuities), in a similar manner to that discussed by David, Bregman & Seab (1988) for TI in $`10^6`$$`10^8`$-K gas, or by Rees (1977) in the case of gravitational contraction. It is then important to determine whether large, static, thermal pressure-confined clouds ever develop and, if yes, how long does it take for the dynamic situation to fade away. If this time is long compared to estimated cloud lifetimes \[see Blitz & Shu (1980) and Ballesteros-Paredes, Hartmann & Vázquez-Semadeni (1999) for two different possibilities in the case of molecular clouds\], or to the time scales for major disturbances to pass through the cloud (e.g., Kornreich & Scalo 2000), then the static configuration should not be expected to occur in general in the ISM. Unfortunately, Paper I did not proceed to give a definitive solution to this question because of limitations of the numerical scheme used there. Most previous discussions of the time scales to reach a developed cloud-intercloud medium have been based on the assumption of isobaric evolution (e.g., Parravano 1987), which, however, is strictly valid only in the limit of $`\eta 1`$. Goldsmith (1970) did consider the limit $`\eta 1`$, and mentioned that very long times are required for the quasi-static regime to be reached, such that his simulations (at very low resolutions) never reached it (see also Schwarz et al. 1972).

Small-scale perturbations, on the other hand, are good candidates for rapid, yet not highly dynamic condensation, and constitute the paradigm of cloud formation via the isobaric mode of TI. This brings up the questions of a) whether, in the presence of a spectrum of density fluctuations, the small-scale perturbations can use up most of the available mass before the large-scale ones grow, and thus impede their development, and b) whether the paradigm is maintained when the perturbations are produced by turbulent velocity fluctuations. In the regime of Paper I, in which the cooling times are short, turbulence forced at large scales cannot prevent the development of TI, and in fact actually enhances it, although stellar-like heating at small spatial scales was shown to erase the signature of TI in the density histogram. The latter result does not strictly constitute an inhibition of TI, because the stars form after the condensations, which in turn formed by a combination of the turbulent compressions and TI. Such condensations, in the absence of the stellar heating, would pile up in the thermally-stable density regime. However, true avoidance of the condensation may occur for random velocity fluctuations at the small scales because in this case the compressions tend to heat the gas; i.e., the perturbations are quasi-adiabatic, and can re-expand before they have time to cool.

The discussion above is related to the question of whether the classical picture of cold, dense clouds in pressure equilibrium with their warm, dilute surroundings is a realistic model of the ISM, even if TI were the sole cloud-forming agent in the ISM (which it is not). Indeed, the equilibrium configuration is a possible solution of the hydrodynamic equations in the presence of TI-inducing cooling functions; it is very often used as an initial condition of both analytical and numerical studies, and even as a textbook case. However, the question is whether it is also a self-consistent final state of previous evolution starting from small-amplitude perturbations of thermally unstable gas. Of course, the pressure equilibrium configuration is valid as an initial condition for numerical simulations or as a case study for theoretical analyses, just as Gaussian random fluctuations are in studies of turbulence or Cosmology, but they are not to be interpreted as realistic representations of the ISM. For example, shock-cloud collision simulations (e.g., Klein, McKee & Woods 1995) often start with a cloud in pressure equilibrium with its environment, but end up with fully developed turbulence, suggesting that the latter is a more likely regime for the ISM than their idealized initial conditions.

In this paper we explore the nonlinear development of TI in the isobaric mode, in the context of the atomic ISM, focusing on the issues mentioned above, using high resolution numerical simulations in one dimension (1D) and a numerical scheme better suited for this type of problems than that used in Paper I. We should emphasize that in this paper we choose to isolate the effects of the instability, and concentrate on a somewhat idealized case (one-dimensional geometry, simple parameterized cooling functions neglecting non-equilibrium chemical evolution, neglect of other physical agents such as the magnetic field, differential disk rotation, self-gravity and realistic, stellar energy injection), not necessarily representative of the full situation in the ISM. We make no attempt to discuss the fate of the clouds formed by TI in the real ISM, as this must be accomplished by means of simulations of the fully turbulent ISM in the presence of the agents mentioned above (see Papers I, II and III), and is likely to involve dynamical instabilities whose study requires two- or three-dimensional simulations (see, e.g., Murray et al. 1993; Vietri, Ferrara & Miniati 1997). Instead, our main goals here are simply to illustrate, even in this simplified case, the following features of the condensation process: a) the dynamical nature of the condensation process, providing quantitative estimates of the perturbation sizes in the atomic ISM required for such dynamical evolution to occur, and of the times required for it to subside; b) the existence of accreting gas, which remains at densities in the unstable range, even after the condensation process has reached a slowly evolving stage; c) the competition between growing modes of different sizes; d) the evolution under the isobaric mode of TI of perturbations caused by random velocity fluctuations, and e) their implications on standard scenarios of a multiphase ISM in which most of the gas must reside in one of the stable phases.

The plan of the paper is as follows. We first present the numerical method and discuss the resolution requirements (§2). Next, we present a qualitative discussion of the timescales involved and the regimes that can appear depending on the ratio $`\eta `$3), to then consider the evolution of Gaussian density perturbations, the size scales at which shocks form at advanced stages of the “condensation” process, and the times required to reach quasi-static states (§4.1). We then discuss the competition between the growth of small- and large-scale perturbations by means of simulations with an initial spectrum of density fluctuations (§4.2), and the evolution of turbulent fluctuations in the presence of TI by means of simulations including random forcing (§4.3). In §5 we compare our results with previous work, and discuss the limitations of the one-dimensional simulations we have used as well as some implications of our results on the classical picture of thermal-pressure confinement of clouds, the existence of gas in the unstable regime, and the compressibility of the warm gas in the ISM. Finally, we summarize our results in §6.

## 2 The model

### 2.1 Numerical method and relevant parameters

We solve the hydrodynamic equations for the logarithmic density $`\mathrm{ln}\rho `$, the velocity $`𝒗`$, and the entropy per unit mass $`s`$. In terms of these variables, the mass, momentum and energy conservation equations used to calculate the complete evolution of a thermally unstable medium are

$$\frac{D\mathrm{ln}\rho }{Dt}=\mathbf{}𝒗,$$
(1)
$$\frac{D𝒗}{Dt}=\frac{\mathbf{}P}{\rho }+𝒇(𝒓,t)+\text{viscous force},$$
(2)
$$T\frac{Ds}{Dt}=\mathrm{\Gamma }\rho \mathrm{\Lambda }+\text{viscous heating},$$
(3)
where $`D/Dt=/t+𝒗\mathbf{}`$ is the advective derivative, $`𝒇`$ is a forcing function to be specified that will be used only in the case of forced simulations (§4.3), and $`\mathrm{\Gamma }`$ and $`\mathrm{\Lambda }`$ are the heating and cooling functions, respectively, depending on $`\rho `$ and $`T`$. The equation of state used is that of ideal gases $`P=(\gamma 1)c_\mathrm{v}\rho T`$, where $`\gamma =5/3`$ is the heat capacity ratio $`c_\mathrm{p}/c_\mathrm{v}`$ for monoatomic gases. The specific entropy is related to the density and pressure by $`s=s_{00}+c_\mathrm{v}\left[\mathrm{ln}(P/P_{00})\gamma \mathrm{ln}(\rho /\rho _{00})\right]`$, where $`\rho _{00}`$, $`P_{00}`$ and $`s_{00}`$ are reference values for the density, pressure and entropy.

We use simplified forms for the cooling and heating functions, consisting of a piecewise power-law fit to the “standard” equilibrium $`P`$$`\rho `$ curve of Wolfire et al. (1995, hereafter WHMTB) and a constant background heating rate per unit mass. This choice is at an intermediate level in a hierarchy of increasing complexity in the treatment of the thermal processes in a gas dynamical problem, that goes from assuming isothermal (e.g., studies of molecular clouds by Gammie & Ostriker 1996; Mac Low et al. 1998; Padoan & Nordlund 1999; Ostriker, Stone & Gammie 2001) or polytropic flows (e.g., Vázquez-Semadeni, Passot & Pouquet 1996), to solving the heat equation using simple fits to the cooling and heating functions (e.g., Chiang & Bregman 1988; Rosen & Bregman 1995; Vázquez-Semadeni, Passot & Pouquet 1995, 1996; Paper I; Burkert & Lin 2000; Koyama & Inutsuka 2002; this paper), to following the detailed non-equilibrium chemistry and the cooling rates (e.g., David et al. 1988; Murray & Lin 1989, 1990; Kang et al. 1990; Koyama & Inutsuka 2000). For the purpose of studying the dynamical development of TI from the equilibrium state in the neutral ISM, it is not clear to what extent it is imperative to follow the detailed non-equilibrium chemistry and the resulting rates (see §5.1 for further discussion). Thus, the present work should be taken mainly as an illustration of the highly dynamical and complex regimes that can already appear within this approximation, and the caveat of the simplified cooling functions we use should be kept in mind.

Our fit to WHMTB’s standard cooling curve has the form

$$\mathrm{\Lambda }=C_{i,i+1}T^{\beta _{i,i+1}}\mathrm{for}T_iT<T_{i+1}.$$
(4)
As mentioned above, we parameterize the heating function by $`\mathrm{\Gamma }=\mathrm{\Gamma }_0(\rho /\rho _0)^\alpha `$, with $`\alpha =0`$, since the dominant heating process (photoelectric heating from small grains and PAHs) in the density range of interest depends only weakly on density, roughly as $`\rho ^{0.2}`$ (WHMTB).

We find the following set of $`\beta `$-values and transition temperatures: $`\beta _{12}=2.12`$, $`\beta _{23}=1.0`$, $`\beta _{34}=0.56`$, $`\beta _{45}=3.67`$, and $`T_1=10`$ K, $`T_2=141`$ K, $`T_3=313`$ K, $`T_4=6102`$ K, $`T_5=10^5`$ K. Note that for $`T_3<T<T_4`$ the gas is thermally unstable under the isobaric criterion, and marginally stable for $`T_2<T<T_3`$. For a constant heating of $`\mathrm{\Gamma }_0=0.015`$ erg s-1g-1, the values of the coefficients of the cooling functions in units of erg s-1g-2cm3K$`^{\beta _{i,i+1}}`$ are $`C_{1,2}=3.42\times 10^{16}`$, $`C_{2,3}=9.10\times 10^{18}`$, $`C_{3,4}=1.11\times 10^{20}`$ and $`C_{4,5}=2.00\times 10^8`$. The fitted curve is displayed in Figure 1 (c.f. fig. 3 of WHMTB).

We will denote by $`P_{\mathrm{eq}}(\rho )`$ the pressure corresponding to thermal equilibrium at density $`\rho `$, by $`P_{\mathrm{max}}`$ the equilibrium pressure corresponding to the high-temperature boundary of the thermally-unstable range, and by $`P_{\mathrm{min}}`$ the low-pressure plateau in the thermal equilibrium curve (see Fig. 1). As “standard” unperturbed initial conditions we take $`\rho _0=1`$ cm-3, representative of the mean density in the ISM, and $`T=2400`$ K, the thermal equilibrium temperature corresponding to $`\rho _0`$. This state falls in the unstable range. Other workers have chosen to apply nonlinear perturbations (strong compressions or shocks; e.g., Hennebelle & Pérault 1999; Koyama & Inutsuka 2000, 2002) to linearly stable gas in order to determine whether TI can be triggered in such regimes by strong perturbations. However, since in this paper we are interested in the conditions necessary for the development of TI to be highly dynamic and in the competition between modes, we prefer to take the initial conditions described above, which, we remark, are close to the mean density of the ISM (see, e.g., Ferrière 2001).

We use a sixth-order finite difference scheme for the spatial derivatives on a one-dimensional, non-uniform Cartesian grid, and a third-order time stepping scheme, as described in Sánchez-Salcedo & Brandenburg (2001). The code includes a shock-capturing viscosity proportional to $`\mathbf{}𝒗`$ wherever $`\mathbf{}𝒗<0`$, and zero otherwise. Shocks are resolved over $`3`$$`4`$ zones. To avoid numerical instability due to amplitude errors at high-wavenumber modes of density and entropy, diffusion terms in the continuity and entropy equations have been introduced. Convergence tests (§4.1.1) have been performed to ensure that mass and thermal diffusion do not affect the growth rate of forming density peaks. In particular, the effective mass diffusion scale implies that below scales $`10`$ pixels the growth rate starts to decrease monotonically (see also Paper III).

Most of the simulations presented here assume a plane-parallel geometry. This is motivated by the need to achieve very high resolutions, but is justified by evidence that suggests the presence of sheet-like and filamentary morphologies in the cold neutral medium (e.g. Kim et al. 1998). Although the origin of this morphology may not be caused by TI, we consider this geometry to be more appropriate than, for example, a spherical one for studying the development of TI in the neutral medium of the Galaxy. We discuss this issue further in §5.2.

The non-uniform grid allows for enough resolution as to resolve the structure of the cloud and the accretion fronts. It also allows us to move the boundaries far away, so no perturbation has time to reach them. For simulations of forced flows, we use a uniform mesh with periodic boundary conditions.

### 2.2 Resolution requirements

One condition on the the minimum resolution necessary in the simulations is obtained by requiring that the final condensations be well resolved. Consider a domain of size $`l`$ per dimension covered with a uniform mesh containing $`n_x`$ grid points, in the $`x`$-direction, separated by a distance $`\delta x`$, and an initial density perturbation with a characteristic size $`l_0=\psi l`$. In 1D, the final size of the condensation, $`l_1`$, will be $`l_1l_0\rho _0/\rho _1=\psi l\rho _0/\rho _1`$, where $`\rho _0`$ and $`\rho _1`$ are respectively the mean and final densities. If we require a minimum number of zones, $`n`$, inside $`l_1`$, the total number of grid points must be larger than $`n\rho _1/(\psi \rho _0)`$. For $`\rho _1/\rho _0=50`$, $`n=10`$, and $`\psi =0.5`$, we need $`n_x>1000`$. If this condition is not met, our tests indicate that a stationary regime is reached in which mass diffusion balances accretion, leading to a very unrealistic situation.

If one wants to resolve the out-of-equilibrium thermal structure within accretion fronts, a second condition must be satisfied, namely that $`\delta x`$ must be smaller than the cooling length. The cooling length is given by $`l_\mathrm{c}=\tau _{\mathrm{cool}}v_{\mathrm{fr}}`$, where $`\tau _{\mathrm{cool}}`$ is the cooling time and $`v_{\mathrm{fr}}`$ the front velocity. The condition reads $`n_x>l/(\tau _{\mathrm{cool}}v_{\mathrm{fr}})`$. This is a highly variable condition, as both $`\tau _{\mathrm{cool}}`$ and $`v_{\mathrm{fr}}`$ decrease with increasing density, so the condition is most stringent at high densities. For $`\tau _{\mathrm{cool}}0.2`$ Myr and $`v_{\mathrm{fr}}=0.1`$ km s-1 we need $`n_x>10^4`$ for $`l=100`$ pc and $`n_x>10^5`$ for $`l=1`$ kpc. The need for huge resolution can be alleviated by using a non-uniform grid, with a maximum density of points where the condensation is expected to form.

An alternative approach is to reduce the heating and cooling coefficients by some factor. This modification of the cooling rates is justified only if $`\eta 1`$ at all scales of interest, and at all times, since, in effect, this is equivalent to reducing the simulation box size. However, in regions where the gas is evacuated, the density decreases and the cooling times increases. Hence, it is not true that the condition $`\eta 1`$ holds at all times in the warm diffuse gas, and we cannot reduce the heating and cooling coefficients by the same factor everywhere, but instead we would have to use a density-dependent function which reduces the cooling rates only at very high densities. For these reasons, we prefer to stick to non-uniform grids in simulations of single condensations, and to very-high-resolution uniform grids in multi-fluctuation simulations. Note, however, that these high resolution requirements arise only in order to resolve such small-scale structures as the cloud formed or the thermal structure within radiative shocks. Nonetheless, much lower resolutions still suffice to describe the global dynamics.

## 3 Time scales and general features

In this Section, we describe some basic concepts useful for understanding the remainder of the paper, including the characteristic time scales involved, as well as some features of the growth of fluctuations. We also define the relevant time scales explicitly in the cases of density and velocity fluctuations.

### 3.1 Modes of the instability and the dynamical evolution of fluctuations

As mentioned in the Introduction, in this paper we are concerned with cooling functions which imply instability only for the isobaric mode. In the presence of constant heating ($`\alpha =0`$), this is satisfied for $`\beta <1`$, in which $`\beta `$ is the exponent of the temperature in the cooling function (see §2). For our fit to the standard cooling function of WHMTB, this occurs in the temperature range $`300`$ K $`T6100`$ K, where $`\beta =0.52`$, i.e., between the so called “cold” and “warm” phases of the ISM. In the presence of constant heating, the isochoric mode of TI is unstable for $`\beta <0`$, and thus for temperatures between $`10^5`$ and $`10^6`$ K for classical cooling functions (e.g., Raymond, Cox & Smith 1976). Such temperatures are not reached by our simulations, since we do not include supernova heating. Finally, it can be easily shown that in our case the adiabatic criterion requires $`\beta <1/(1\gamma )=3/2`$ for instability, and thus our simulations are stable under the latter two criteria in the vicinity of thermal equilibrium. Note that, however, far from thermal equilibrium, the system may cool freely, and then the isochoric mode may also become unstable, as in the cooling case (i.e., in the absence of background heating), in which the instability criteria for the isobaric and isochoric modes are, respectively, $`\beta <2`$ and $`\beta <1`$.

### 3.2 Evolution of isobaric density perturbations

Here we consider the growth of initial states with density perturbations of a given physical size, no fluid motions and constant pressure. Isobaric density perturbations constitute a simple case of entropy perturbations, as they imply a perturbation of the ratio $`P/\rho ^\gamma `$, and can be produced in practice by either locally varying the cooling or heating rates, or by velocity fluctuations well into the subsonic regime (see §§3.3 and 4.3).

Let us consider, for simplicity, the condensation process of a perturbation $`\rho _1=\rho \rho _0>0`$. For a linear perturbation of wavelength $`\lambda `$ there are two relevant timescales at $`t=0`$: the sound crossing time $`\tau _\mathrm{s}(0)\lambda /c_\mathrm{s}(0)`$, where $`c_\mathrm{s}(0)`$ is the adiabatic sound speed, and the cooling time $`\tau _{\mathrm{cool}}(0)c_\mathrm{v}/\left|\left(\rho \mathrm{\Lambda }/T\right)_P\right|`$, which is of the order of $`e_0/\rho _0\mathrm{\Lambda }`$, where $`e_0`$ is the initial internal energy per unit mass. The initial time scale ratio is then $`\eta _0\tau _{\mathrm{cool}}/\tau _\mathrm{s}`$.

It is well known (Field 1965; see also Shu 1992; Balbus 1995; Meerson 1996; Burkert & Lin 2000) that entropy perturbations with $`\eta >1`$ evolve almost quasi-statically, maintaining quasi-isobaricity throughout their evolution. On the other hand, entropy perturbations with $`\eta <1`$ cool before they can equate pressures and therefore develop large pressure gradients, which then cause them to evolve highly dynamically even if only the isobaric instability criterion is satisfied. It is also well known that perturbations with $`\eta >1`$ grow on time scales $`\tau _{\mathrm{cool}}`$, because in this case the perturbation can only evolve as fast as it cools, and thus the growth rate saturates at $`1/\tau _{\mathrm{cool}}`$ in the limit of small perturbation scales. Instead, when only the isobaric criterion is satisfied, perturbations with $`\eta <1`$ grow on a time $`\tau _s`$, because they can only evolve as fast as the pressure gradient pushes the gas surrounding the condensation to equate pressures (see Field 1965 and Meerson 1996 for more details). This implies that the growth rate approaches zero in the limit of long wavelengths. In order for finite growth rates to exist in this limit, it is necessary to also satisfy the isochoric criterion.

If the evolution were strictly isobaric, a condensation (cloud) would end up acquiring the density denoted $`\rho _{\mathrm{isob}}`$ in Figure 1. Adopting the cooling function described in §2.1 and a mean interstellar density of $`1`$ cm-3 ($`T_{\mathrm{eq}}2400`$ K), $`\rho _{\mathrm{isob}}=40`$ cm-3. Furthermore, the condition $`\tau _\mathrm{s}\tau _{\mathrm{cool}}`$ occurs at a size scale $`\lambda 10`$ pc (hereafter denoted by $`l_{\mathrm{eq}}`$). Therefore, isobaric evolution cannot occur at scales larger than a few parsecs. Let us thus consider a perturbation (of large enough size scale) such that $`\eta 1`$. Since $`\tau _{\mathrm{cool}}\rho ^{1/\beta }`$ and $`\tau _\mathrm{s}`$ cannot decrease faster than $`\rho ^{1/D}`$ during the condensation process, where $`D`$ is the dimension of the contraction, we conclude that $`\eta 1`$ in the core of the cloud (i.e. at the vicinity of $`x0`$) at all subsequent times, even as it shrinks. This is in sharp contrast with the fact that smaller scales have larger values of $`\eta `$, under the mean conditions of the unperturbed medium.

Relatively high local Mach numbers can be attained in the condensation process under conditions with $`\eta _01`$. An estimate of the maximum velocities that can be reached due to pure TI can be made from Bernoulli-like arguments as $`v_{\mathrm{inst}}^22(P_{\mathrm{max}}P_{\mathrm{min}})/\rho _0`$, leading to fluid velocities $`v_{\mathrm{inst}}4`$ km s-1, for $`\rho _0=1`$ cm-3, $`P_{\mathrm{max}}/k=3000`$ K cm-3 and $`P_{\mathrm{min}}/k=1000`$ K cm-3. This velocity is clearly supersonic in a medium at the temperature corresponding to $`P_{\mathrm{min}}`$, at which $`c_\mathrm{s}2`$ km s-1. With that typical velocity, a region of size $`L`$ ($`L_0`$ at time $`t=0`$) condenses on a characteristic timescale

$$\tau _{\mathrm{inst}}=\frac{L_0}{v_{\mathrm{inst}}}=1.2\left(\frac{L_0}{10\mathrm{pc}}\right)10^6\mathrm{yr}.$$
(5)
The characteristic time given by Eq. (5) can be compared to the free-fall time of a uniform cloud due to self-gravity,

$$\tau _{\mathrm{ff}}=\left(\frac{3\pi }{32G\rho }\right)^{1/2}=4.4\left(\frac{\rho }{1\mathrm{cm}^3}\right)^{1/2}10^7\mathrm{yr}.$$
(6)
Therefore, self-gravity must be included either at scales larger than $`0.4`$ kpc or at late stages due to the density enhancement. In a plane-parallel condensation $`\tau _{\mathrm{ff}}(t)L^{1/2}`$, so that the ratio $`\tau _{\mathrm{ff}}/\tau _{\mathrm{inst}}`$ increases during the condensation process. Thus, for the fiducial values we have used, the flow may reach the thermally stable regime and re-approach pressure equilibrium, leaving only a slow-accretion regime (see §4.1), without self-gravity ever becoming dynamically important. Consequently, overshooting of density will be unable to prompt gravitational instabilities in 1D experiments, but it may be important in more than 1D (e.g., Schwarz et al. 1972; Kang et al. 2000).

### 3.3 Fluctuations in the velocity field

If, in contrast with the cases described above, an initial state of constant density is perturbed with velocity fluctuations (of amplitude $`v_0`$), the evolution may be quite different. The generalization of the instability analysis to a non-static initial state was performed by Hunter (1970, 1971). Velocity perturbations as triggers of TI in initially stable media have been investigated by a number of authors (e.g., Murray & Lin 1989; Kang et al. 1990; Hennebelle & Pérault 1999, 2000; Koyama & Inutsuka 2000), but here we are interested in the opposite effect, i.e., whether perturbations induced by velocity fluctuations in initially unstable media can actually be stable, as first suggested, but not shown, in Paper I. A thorough discussion has been presented in Paper III, and here we just repeat the main points.

In the presence of velocity fluctuations, a new characteristic time scale appears in the system, namely the bulk-velocity crossing time, $`\tau _\mathrm{u}=\lambda /u`$, where $`u`$ is the characteristic velocity difference across scale $`\lambda `$. Thus, the dynamical time $`\tau _{\mathrm{dyn}}`$ to use in $`\eta `$ should be chosen as the minimum of the sound and the bulk velocity crossing times.

Compressive motions in the absence of radiative cooling cause adiabatic perturbations in which the density and pressure fluctuations have the same sign, because of the $`PdV`$ work done on the affected fluid parcel. In the presence of cooling, the actual behavior depends on (the redefined) $`\eta `$. If the cooling time is much longer than the characteristic time of the compression (i.e., if $`\eta 1`$), the cooling is negligible, and the compressed parcel behaves as a quasi-adiabatic sound wave. This case is known to be unstable only if the adiabatic criterion is satisfied, which for our piecewise power-law cooling function requires $`\beta <1/(1\gamma )=3/2`$. In the opposite case ($`\eta 1`$), the $`PdV`$ heating is negligible in comparison with the cooling, and the pressure is determined by the thermal equilibrium condition. In this limit, the velocity perturbations are unstable when the isobaric condition is satisfied, because they behave as a condensation mode (see, e.g., Shu 1992; Paper III). Thus, velocity fluctuations strong enough that they dominate $`\eta `$ are linearly stable at $`\eta >1`$ (small scales) but unstable at $`\eta <1`$ (large scales) in the atomic ISM. Moreover, when $`\tau _\mathrm{u}<\tau _\mathrm{s}`$, the instability of large-scale perturbations manifests itself in a null resistance of the flow to turbulent compressions rather than in spontaneous growth of the perturbations, because the bulk velocity is larger than the sound speed, and so the driver of the density growth is the turbulent velocity compression rather than the thermal pressure gradient.

It is worth remarking that velocity fluctuations are a natural way of producing the density fluctuations that constitute the initial conditions for the subsequent development of TI, because the gas obeys to the continuity equation (mass conservation), and thus the production of a density enhancement requires a converging velocity field.

The case of strongly nonlinear compressions probably depends on the closeness of $`\eta `$ to unity. If $`\eta 1`$, then even very strong compressions behave adiabatically, and are not able to induce very large density enhancements. If $`\eta `$ is still larger than unity at the higher density, triggering the isobaric mode will be extremely difficult. On the other hand, if $`\eta `$ is not much larger than unity initially, then the density increase induced by the compression may raise the cooling rate to large enough values that $`\eta `$ becomes smaller than unity, and triggering of the isobaric instability ensues.

## 4 Numerical results

### 4.1 Evolution of a single density perturbation and final stages of TI

In this section we attempt to quantify the dynamical nature of the condensation process in the context of the formation of interstellar clouds by TI. In particular, we seek to determine the physical perturbation sizes at which the condensation develops supersonic velocities, the times needed for these velocities to subside, and the state of the leftover diffuse gas, which, as we shall see, contains significant fractions of mass still in the unstable density range.

Let us consider first how the nonlinear evolution of a Gaussian density perturbation depends on $`\eta _0`$. To this end, we numerically follow the evolution of perturbed density fields of the form

$$\rho _0(x)=1+A_0\mathrm{exp}\left(\frac{x^2}{2\sigma ^2}\right),$$
(7)
where the density is in units of cm-3. We take $`A_0=0.2`$ and consider three values of $`\sigma `$ such that the FWHM of the Gaussian profile is $`3`$, $`15`$ and $`75`$ pc. We refer to the corresponding simulations as DEN3, DEN15 and DEN75. According to our discussion in §3.2, it is expected that for model DEN3, the evolution will be closer to isobaric and quasi-static ($`\eta _03.7`$, assuming $`\tau _\mathrm{s}=\mathrm{FWHM}/2c_\mathrm{s}`$), whereas dynamical compression regimes are expected for DEN15 ($`\eta _00.75`$) and DEN75 ($`\eta _00.15`$). All these simulations start with zero velocity.

Model DEN3 starts with constant pressure. The initial pressure is $`P/k=2500`$ K cm-3, which equals the equilibrium pressure $`P_{\mathrm{eq}}`$ at $`1`$ cm-3. Figure 2 shows the density, velocity and pressure profiles at four different temporal snapshots. The maximum of the velocity, $`v_{\mathrm{max}}(t)`$, increases continuously until $`t3.8`$ Myr, at which time $`v_{\mathrm{max}}0.86`$ km/s (Mach number $`0.18`$). As expected, although the pressure is not perfectly constant (see Fig. 2b), the evolution can be considered as quasi-isobaric, as can be seen in Figure 4. Indeed, the maximum pressure variation is $`\left|\delta P\right|/P_00.2`$. As a consequence, there is no appreciable density overshoot beyond $`\rho _{\mathrm{isob}}`$3.2), and the final pressure inside the cloud recovers a value very close to the ambient pressure.

By time $`t4`$ Myr, a nearly stationary “slow-growth” stage is reached, in which most of the mass is in the central condensation at nearly zero velocity, with low density material falling in from the surrounding medium, but decelerating smoothly as it enters the cloud, without the formation of shock. The time for reaching this stage increases significantly when the initial perturbation amplitude is decreased ($`5.5`$ Myr for $`A_0=0.1`$). In order to quantify the amount of mass within the condensation, we integrate the density over the length of the region where $`\rho 10\rho _0`$, leading to a column density333Note that the units of column density are M pc-2 even though the simulations are one-dimensional, because the physical units of the density are g cm-3 anyway. $`N0.045`$ M pc-2 at $`t=4.2`$ Myr. This value is 67% of the column density initially within $`|x|<\sigma `$.

It is important to remark that in the “final” slow-growth state mentioned above, the leftover diffuse gas is far from static, with velocities of up to 0.7 km s-1. Although the accretion rate onto the condensation is low because of the low density of the accreting gas, the amount of mass contained in this regime is still large, $`8`$ times larger than the amount of mass in the condensation. This result is in part an artifact of the fact that we considered a single perturbation, which by itself does not consume a large amount of mass by the time it reaches $`\rho _{\mathrm{isob}}`$. Smaller final fractions of mass still in the unstable range are expected in the case of multiple perturbations, but even in this case, much longer times are needed to evacuate the unstable range (§4.2). Smaller unstable fractions are also expected if the background medium is already in the stable diffuse phase and the condensation is triggered by a strong compression (we thank P. Hennebelle for pointing this out).

Another point worth remarking is that the accreting gas evolves in a nearly isobaric regime in response to the instability, but such that the small existing pressure variations have the same sign as the density variations (Fig. 2), and thus this gas has no further tendency to fragment, even though its density lies in the “unstable” range. Moreover, this implies that this density range does not precisely correspond to the unstable temperature range, as shown in Fig. 3.

Let us now consider the case of models DEN15 and DEN75. In these, we set up the initial pressure to the corresponding $`P_{\mathrm{eq}}`$ values, $`P_0=P_{\mathrm{eq}}\left(\rho _0(x)\right)`$, which is reasonable for models with low $`\eta `$’s, i.e. short cooling times. Figure 5 displays four snapshots of DEN15. We can distinguish three evolutionary stages: (1) An initial compression stage in which all the parcels in the flow are in the thermally unstable range, the fluid is continuously accelerated, and the pressure gradient becomes increasingly steeper, eventually generating transonic infall motions. (2) A repressurization, or “crushing”, stage, starting when the core of the cloud enters the thermally stable temperature range and the pressure starts to increase in the center of the cloud, while the outside material is still unstable and continues to accelerate inwards, generating a shock at the edge of the cloud at the time of maximum compression. Note that this is a radiative shock occurring over distances of the order of the cooling length, and so it is well resolved by the numerical grid. (3) A slow, subsonic accretion stage, which slowly fades away. During the initial compression stage the system is closer to thermal equilibrium than in DEN3 (see Fig. 6). The maximum Mach number achieved in this run was $`0.6`$, and occurred at $`6.1`$ Myr. It is necessary to start with a perturbation of FWHM of $`40`$ pc for the infall velocity to be transonic.

The shock wave that emanates from the cloud boundary is worth discussing. This shock is rather weak, as can be seen from the pressure jump in Figure 5b, but has important effects in the warm medium. As shown in the $`P`$$`\rho `$ diagram (Fig. 6), as the shock propagates outwards, it heats the inflowing gas, bringing it closer to an isobaric regime, similar to the one of run DEN3, while bringing it out of thermal equilibrium. In fact, the relative pressure variation $`\delta P/P_0`$ from $`x=0`$ to the position just behind the shock is only $`0.2`$ in the slow accretion stage, whereas the temperature varies from $`50`$ to $`3000`$ K over this same interval. The column density of the cloud is $`0.26`$ M pc-2, about $`76\%`$ of the initial column density within $`|x|<\sigma `$ ($`0.335`$ M pc-2). Since the ram pressure in the boundary of the cloud is larger than in model DEN3, the cloud density is slightly larger too.

We now consider run DEN75, for which $`\eta _00.15`$. Figure 7 shows the evolution of density, pressure and velocity. The three evolutionary stages discussed in model DEN15 are also present in DEN75. In this run, the gas behaves precisely along the thermal equilibrium $`P`$-$`\rho `$ curve in the initial collapsing stage ($`t21`$ Myr), as we can see in Figure 8. In this situation, the gas acquires the maximum possible values of the infall velocity imposed by the development of pure TI. In fact, at the time when the cloud begins to reach the thermally stable phase ($`T140`$ K), the maximum velocity is $`3.0`$ km s-1, and the maximum local Mach number is $`1.2`$. The maximum of the velocity continues increasing up to $`3.6`$ km s-1, agreeing with our estimate for $`v_{\mathrm{inst}}`$ in §3.2 to within 10%. Consequently, the supersonic crushing process is more violent than in DEN15, causing a large density overshoot, of roughly $`30\times \rho _{\mathrm{isob}}`$. Note that the actual value of the density peak in these simulations is limited by shock viscosity and mass diffusion, although this value is of little interest in 1D contractions because, as mentioned in §3.2, it is not viable for triggering gravitational instabilities, except, perhaps, if the densities are high enough to trigger a non-equilibrium instability in the cold gas due to molecular cooling (Koyama & Inutsuka 2000), which we do not include here. The density of the cloud will subsequently decrease slowly with time, but it will remain significantly larger than $`\rho _{\mathrm{isob}}`$ (see below).

During the compression phase, the small-scale Fourier components of the Gaussian perturbation blow up due to their larger growth rates, forming two (symmetric) prominent density peaks, which ultimately merge with each other (Fig. 7; note that only one peak is seen, since only one side of the simulation is shown). Their formation depends on details of the initial conditions, being nonexistent in the presence of a small initial convergent velocity field. The overall evolution of a simulation (DEN75b hereafter) with FWHM=75 pc, $`A_0=0.7`$ and initial velocity field $`v_0P_0/\rho _0`$, with initial maximum Mach number of $`0.4`$, does not present remarkable differences except for the presence of stronger accretion.

As with run DEN15, at later times, the shock wave propagates outwards from the cloud boundary, causing a temperature jump of $`\mathrm{\Delta }T300`$ K across it, nearly restoring pressure equilibrium, albeit at the expense of throwing the outside gas out of thermal equilibrium. Figure 9a shows the density, the bulk and sound speeds, and the cooling times out to $`75`$ pc from the center at $`27.0`$ Myr. Within the region $`1<|x|<75`$ pc from the center, the density lies in the range $`0.5<\rho <3`$ cm-3, with cooling times between $`0.8`$ and $`2.7`$ Myr. Within $`20<|x|<75`$ pc, the cooling time is roughly constant ($`\tau _{\mathrm{cool}}2.0`$$`2.7`$ Myr), and the typical value of $`\eta `$ with respect to the size of this region (55 pc) is $`0.8`$$`0.9`$. This value of $`\eta `$ explains the approximately isobaric evolution of the flow across this region, which is characterized by relatively low densities ($`\rho 0.7`$ cm-3). At intermediate densities, i.e. within $`3<|x|<20`$ pc, $`\eta `$ also has values $`1`$ (again with respect of the region size, $`17`$ pc), because both $`\tau _{\mathrm{cool}}`$ and $`\tau _\mathrm{s}`$ decrease towards the cloud. Thus, this simulation has evolved towards values of $`\eta `$ closer to unity over most of the medium external to the cloud, except at the cloud core ($`|x|<0.25`$ pc), where the condition $`\eta 1`$ is preserved (§3.2). Note also that for perturbations significantly smaller than the sizes of the regions mentioned above, $`\eta `$ is actually $`1`$ in the warm medium.

The stage of slow accretion at the cloud boundary is established after $`25`$ Myr. In Figure 9b the central part of the simulation, including the accretion front, is shown at $`t=27.0`$ Myr. It is worthwhile to note that there is still an appreciable excess of density with respect to $`\rho _{\mathrm{isob}}`$ at that time ($`\rho >2\rho _{\mathrm{isob}}`$). The column density of the cloud is $`1.65`$ M pc-2 at $`27.0`$ Myr, very similar to the initial column density within $`|x|<\sigma `$ ($`1.68`$ M pc-2), whereas the mass column density remaining in unstable regions with $`\rho <3`$ cm-3 and $`1400<T<6000`$ K (not including the region unaffected by the perturbation) is $`3.2`$ M pc-2. In fact, only a small fraction of the low-density gas has reached densities corresponding to the warm stable phase (Fig. 7) (no gas at all reached that phase in runs DEN3 and DEN15). As mentioned with respect to run DEN3, this is in part a consequence of having considered a single density perturbation in an unstable, homogeneous medium.

The final slow-growth state of this large-scale perturbation is actually very similar to that of DEN3, after the shock has established a quasi-isobaric regime. The gas still accreting has subsonic velocities, is out of thermal equilibrium, and is traversing the “unstable” density range, although with density and pressure fluctuations of the same sign, so that it has no further tendency to fragment. Due to the low accretion rate, long times will be required to evacuate the unstable range, as already pointed out by Goldsmith (1970).

#### 4.1.1 Convergence study

In order to assess the effect of numerical resolution, we have conducted three realizations of each model DEN3, DEN15 and DEN75b, in which the resolutions differed by factors of two444We present here the case of DEN75b because it is the most computationally demanding situation.. For DEN3, the number of grid points, $`n_x`$, were $`800`$, $`1600`$ and $`3200`$, which correspond to a resolution around $`x=0`$ of $`3\times 10^2`$, $`1.5\times 10^2`$ and $`7.6\times 10^3`$ pc, respectively. Almost identical results were found (differences of less than $`1\%`$ in all variables). For DEN15 we used $`1800`$, $`3600`$ and $`7200`$ grid points, with associated resolutions $`7.6\times 10^2`$, $`3.8\times 10^2`$ and $`1.9\times 10^2`$ pc, respectively. Again, the variations were insignificant (less than $`1\%`$). For DEN75b, simulations were conducted with $`2300`$, $`4600`$ and $`9200`$ grid points. In Fig. 10 we compare two snapshots for $`n_x=2300`$ and $`n_x=9200`$ at an advanced time when steep gradients are already present. The density and velocity profiles are similar in both cases, whereas the pressure is slightly less smooth in the low resolution simulation. In the high-resolution simulation the pressure is better resolved but the velocity field is unaffected.

It is possible that the growth of small-size structures consisting of a few grid points are limited by resolution. To be sure that the results are not sensitive to resolution, experiments with increasing resolution have been carried out for all of the simulations presented throughout this paper. High enough resolution has been used to ensure that it does not affect the results. See also Paper III for further discussion.

### 4.2 Multi-scale density perturbations

As mentioned above, interesting phenomena may be expected if the initial density profile contains a spectrum of wavelengths (see also Burkert & Lin 2000). Small scale perturbations ($`\eta >1`$) have larger growth rates than larger-scale ones ($`\eta <1`$) and thus will quickly outgrow them. Suppose now that we have two separate scales, one with $`\eta 1`$ and the other one $`\eta 1`$. If the clouds formed from small scales do not consume an important fraction of the mass, their bulk motions will be dominated by the pressure structure of the larger scales. As a consequence, these clouds will behave as if passively advected by the larger ones and will be accelerated towards pressure minima. In this case the evolution may be envisaged as a “linear” superposition of the evolution of small and large scales. In a more general situation, however, the evolution could become extremely complex because of the competition between the growth of perturbations of different sizes. To quantify the relevance of this competition on the dynamics and mass spectrum of the forming clouds, we consider a superposition of fluctuations with wavenumbers in the range $`k_{\mathrm{min}}<k<k_{\mathrm{max}}`$, where $`k_{\mathrm{min}}=2\pi /l`$, with $`l`$ the box size, and $`k_{\mathrm{max}}`$ is varied to consider various extents of the perturbed range of scales. The amplitudes of the fluctuations follow a power-law of the form $`\rho _k=C_0(k/k_{\mathrm{min}})^\delta `$. We consider two values of the ratio $`k_{\mathrm{max}}/k_{\mathrm{min}}`$ and two values of the spectral slope $`\delta `$. Moreover, we consider another set of simulations in a box of size 10 times smaller. Note that, if the amplitude of the initial perturbations is an increasing function of the wavenumber (i.e. $`\delta <0`$), it is clear that the small scales will be able to consume an important fraction of the gas with no chance for large scale perturbations to survive. Thus, we restrict our discussion to $`\delta >0`$, in particular to the values 3/2 and 1/2. For these values of $`\delta `$ we now discuss the conditions under which the small-scale perturbations do not appreciably affect the overall evolution.

First, we compare the evolution of two runs with $`\delta =3/2`$, labeled SP20a and SP4a (for “superposition”; see Table 1). These runs have ratios $`k_{\mathrm{max}}/k_{\mathrm{min}}=20`$ and 4, respectively. We consider a box size of $`250`$ pc, which implies that, for these runs, $`\eta _0(k_{\mathrm{min}})0.04`$. Four snapshots of their evolution are displayed in Figure 11. The number of condensations at intermediate times is $`16`$ and $`5`$ for SP20a and SP4a, respectively. This number depends almost linearly on $`k_{\mathrm{max}}/k_{\mathrm{min}}`$ for these simulations.

The problem of competition for the available mass can be summarized in the question of whether clouds formed from small scales in run SP20a group together to form a density pattern that at large scales is similar to that of run SP4a. If that happens, we can say that perturbations at small scales are unable to suppress the development of TI at larger scales. Note that we speak in terms of Fourier wavelengths only for the initial spectrum. In contrast, once condensations are growing, we focus on their clustering and distribution in physical space.

From panels (c) and (e) of Figure 11, we see that some of the incipient condensations in run SP20a are disrupted when they are embedded in a diverging flow (for instance the clumps between $`x=25`$ and $`x=20`$). Conversely, structures placed in a compression zone grow very fast. Note that some of the clouds in the group near $`x=60`$ are moving relatively fast ($`1.5`$ km/s) at $`25`$ Myr. At later times, the clouds undergo a series of mergers that ultimately lead to four main condensations at positions $`x=110`$, $`x=95`$, $`40`$ and $`60`$ at $`t38`$ Myr.

For this run, the column density ratios of gas mass between the warm, the thermally unstable and the cold phases are $`1:0.5:1`$ at $`15`$ Myr and quickly evolve to $`0.69:0.025:1`$ at $`20`$ Myr, whereas at $`25`$ Myr they are $`0.46:0.01:1`$.

As expected from their larger sizes, the growth rate of the density structures is smaller in run SP4a, but each condensation accretes more matter. In fact, most of the gas is still in the thermally unstable phase at $`15`$ Myr, with column density ratios $`0.9:9.7:1`$. At times $`25`$ Myr, these ratios are $`0.66:0.034:1`$. Still, the fraction of gas in the thermally unstable phase is higher than in SP20 at $`20`$ Myr. At those times ($`25`$ Myr), four main condensations have already formed at $`x=100`$, $`40`$, $`35`$ and $`70`$. Comparing with SP4a, the evolution of run SP20a is somewhat more dynamical, in the sense that clouds acquire larger velocities and are prone to suffer mergers more easily than in run SP4a, in which clouds are more massive and static. The locations of the condensations are the same except for those at $`x=35`$ and $`x=70`$ which have not had time to merge yet; they will do so in an exceedingly long time, at $`t90`$ Myr. We conclude that the final states are not equivalent in this case.

Next, we consider two more sets of simulations using the same density spectrum, but varying the simulation box size. When the box is five times larger, the differences are even more significant. However, if the box size is reduced a factor ten ($`25`$ pc, runs SP20b and SP4b, respectively), so that for $`\eta (k_{\mathrm{min}})=0.4`$, the final states of the two initial spectra are much more similar, as seen especially from their velocity fields (see Fig. 12). Hence, the evolution is not only sensitive to the initial spectrum of the perturbations, but also to $`\eta (k_{\mathrm{min}})`$.

Let us now consider a set of similar experiments but with $`\delta =1/2`$. First, we take $`\eta _0(k_{\mathrm{min}})=0.04`$ (runs SP20c and SP4c). In Figure 13 we see that run SP20c contains several condensations that cannot be associated to any main cloud in SP4c. The most noticeable difference between SP20c and SP4c is that SP20c is able to group condensations into two clusters whereas the four condensations in SP4c would need an exceedingly long time to coalesce, as can be seen by comparing the velocity gradients between positions $`x=110`$ pc and $`x=50`$ pc for the two runs (Figure 13). Hence, the final states of runs SP20c and SP4c are qualitatively very different.

Next, we consider a set of runs with $`\delta =1/2`$ and $`\eta _0(k_{\mathrm{min}})=0.4`$ (runs SP20d and SP4d) and another set with $`\eta _0(k_{\mathrm{min}})=4.0`$ (runs as SP20e and SP4e). Run SP20d still exhibits noticeable differences with respect to SP4d. In Figure 13 we see that there are three condensations between $`x=1`$ and $`x=7`$ pc, generated from small scale perturbations in SP20d, but they are absent in SP4d. These condensations have positive velocities, whilst the gas between $`x=1`$ and $`x=5`$ has the opposite velocity in SP4d. This implies that in this case the emergence of small-scale condensations produces important modifications on the evolution of the gas at large scales because the late-time mass distribution is significantly altered if they are present. However, the opposite conclusion is reached for SP20e and SP4e, which are seen to have more similar states, again judging in particular from their velocity fields (see Fig. 13). The density peaks between $`x=0`$ and $`x=0.75`$ pc in run SP20e do not grow at all between $`6`$ and $`8`$ Myr.

We conclude that both the spectrum exponent and $`\eta (k_{\mathrm{min}})`$ determine the outcome of a superposition of modes of various scales. It appears that, as the initial density power spectrum is shallower, the perturbations need to extend to smaller scales (larger $`\eta `$’s) in order for the small-scale fluctuations not to strongly disturb the growth of the larger-scale ones. In particular, we have found that, for $`\delta =3/2`$, a value $`\eta (k_{\mathrm{min}})=0.4`$ is required, while for $`\delta =1/2`$, a more stringent condition, $`\eta (k_{\mathrm{min}})=4.0`$ is necessary. These critical values could be different for 2D or 3D systems.

It is worthwhile to point out that even though the development of TI may be suppressed at large scales due to mass consumption by small scales, massive clouds may form by coalescence of small clouds (see e.g., Elmegreen 1990, Murray & Lin 1996). As mentioned before, clouds formed at small scales may present appreciable bulk velocities that facilitate the merger of clouds. This effect has recently been seen also in 2- and 3D simulations by Koyama & Inutsuka (2002) and Kritsuk & Norman (2002). The coalescence is exaggerated in one dimensional systems, because the merger rate does not depend on the filling factor of the clouds. However, in more than 1D, the merger rate may still be large if the small-scale clouds are advected towards the condensation center of the large-scale ones (see §5.2). In general, our 1D simulations show that small scale perturbations do not lead to a “forest” of small clouds because they coalesce, building more massive ones.

### 4.3 Simulations of forced flows

We now consider the case of flows forced continuously in time, representative of a turbulent ISM, to determine whether turbulence may oppose the condensation process (cf., §3.3) under some circumstances, and, if so, at which scales and values of the rms Mach number. In a domain of size $`l`$, we force the flow at scales large compared to the final sizes of the condensations. This implies that, once formed, the condensations cannot be disrupted by the random forcing. Note, however, that in higher dimensionalities dynamical instabilities could disrupt the clouds while they form, or after they have formed; see §5.2 for further discussion.

We adopt a random force $`𝒇`$ of the form

$$𝒇(x,t)=\widehat{𝒙}\mathrm{}\left\{N\mathrm{exp}\left[ik\left(t\right)x+i\varphi \left(t\right)\right]\right\},$$
(8)
where $`\mathrm{}`$ denotes the real part, $`k(t)`$ is a time dependent wavenumber, $`\varphi (t)`$ is the phase and $`x`$ the position. Following Brandenburg (2001), we take $`N=f_0c_\mathrm{s}\left[k(t)c_\mathrm{s}/\delta t\right]^{1/2}`$, where $`f_0`$ is a constant factor and $`\delta t`$ is the length of the timestep. The values of $`k/(2\pi /l)`$, where $`l`$ is the box length, and of $`\varphi `$ are selected at each timestep randomly in the ranges $`[3,10]`$ and $`[0,2\pi ]`$, respectively. This implies that the scale of maximum forcing is $`\lambda _\mathrm{f}=l/10`$. The positive exponent (1/2) in $`N`$ implies that the maximum strength of the forcing occurs at the highest wavenumbers $`k`$, i.e., roughly at 1/10 of the box size. A summary of the forced simulations is given in Table 2. There, $`\tau _{\mathrm{cond}}`$ is defined as the time for a condensation to reach a density $`\rho =10\rho _0`$ for the first time, $`e`$ is the ratio of the time-averaged internal energy (with the average performed over $`\tau _{\mathrm{cond}}`$) to its initial value. $`_{\mathrm{rms}}`$ is the rms Mach number, also averaged over $`\tau _{\mathrm{cond}}`$ to increase statistical confidence. Note that $`_{\mathrm{rms}}`$ fluctuates strongly in time, in contrast to $`e`$, which behaves much more smoothly. The different values of $`\tau _{\mathrm{cond}}`$ in one row in Table 2 correspond to runs with different seeds for the random numbers.

Let us first consider experiments with $`l=3`$ pc, in which $`\lambda _\mathrm{f}=l_{\mathrm{eq}}/30`$, since $`l_{\mathrm{eq}}10`$ pc (cf. §3.2). For low Mach numbers (say $`0.1`$), it then holds that the sound speed is larger than the typical values of the velocity fluctuations, with $`\eta 1`$. Consequently, the condensations behave essentially as entropy perturbations, and evolve approximately isobarically (cf. §3). Interestingly, however, the various realizations of run 1, which have $`l=3`$ pc and $`=0.03`$, present a mean $`\tau _{\mathrm{cond}}12`$ Myr, about three times the condensation time for the case with an initial density perturbation of amplitude $`20`$ % and the same size. A large dispersion of $`\tau _{\mathrm{cond}}`$ is apparent for flows with $`_{\mathrm{rms}}0.15`$$`0.3`$, suggesting that the formation of condensations is sensitive to the past history of the flow, becoming somewhat fortuitous for higher Mach numbers. For even larger values of $`_{\mathrm{rms}}`$ ($`0.3`$), eventually the condensation is suppressed altogether (runs 5–9). In these cases, the internal energy initially increases with time until it starts to fluctuate ($`4\%`$) around a value more than twice its initial value. This occurs when cooling is able to compensate the thermal energy injection by the random forces, since the net cooling rate increases also with time as both the strength of the perturbations in density and temperature increase. Specifically, energy injection and net cooling are in balance with a mean value of $`0.038`$ erg s-1 g-1 in run 7 after a time $`0.8`$ Myr, whereas this occurs after a time $`0.35`$ Myr and for a value of $`0.175`$ erg s-1 g-1 for run 9. Our results thus suggest empirically that, for the cooling functions we use, the transition to stabilization occurs roughly at $`e2`$ ($`_{\mathrm{rms}}0.3`$).

On the other hand, for the case $`l=100`$ pc (run 17), it is clear that the forcing is applied at scales where $`\eta 1`$, and indeed the density fluctuations generated by the random compressions rapidly cool and acquire a temperature and a pressure close the thermal equilibrium values. Since the average density is in the unstable range, TI proceeds unimpeded, manifesting itself not so much as a tendency for the fluctuations to grow, but in the null resistance of the flow to external compression (§3.3). Indeed, the random compressions promote the instability, as indicated by the decrease in $`\tau _{\mathrm{cond}}`$ with increasing $`_{\mathrm{rms}}`$. This is the regime of the randomly forced simulations of Paper I, and corresponds to the case of an effective imaginary sound speed of the wave mode, which consequently does not propagate (e.g., Field 1965; Shu 1992).

The case with $`l=10`$ pc (runs 10–16) has $`\lambda _\mathrm{f}l_{\mathrm{eq}}/10`$. Nevertheless, Table 2 shows that the condensation of the density fluctuations cannot be inhibited anyway, even for $`_{\mathrm{rms}}0.6`$. We see that $`e`$ can hardly be increased by factors of $`1.9`$ in this case. We understand this as a consequence of the fact that, for such large Mach numbers, the reduction of the cooling time by the compression-driven density enhancement can reduce $`\eta `$, allowing the triggering of the instability under the isobaric mode by the nonlinear velocity fluctuations. It appears that this does not happen in the case $`l=3`$ at large $`_{\mathrm{rms}}`$ because $`\eta `$ is large enough in that case that the shocks produced are too close to adiabatic and cannot raise the density enough to destroy the adiabaticity. We should remark that for simulations with $`l=10`$ pc and strong forcing ($`_{\mathrm{rms}}>0.5`$, e.g. run 16) compressions can produce transient density fluctuations $`10\rho _0`$, which are dispersed by their own pressure.

We conclude that, as conjectured in §3.3, the growth of perturbations can be suppressed provided that the forcing is applied at scales small enough that $`\eta 1`$ and produces a large enough rms Mach number so that $`\eta `$ is dominated by the turbulent velocity fluctuations rather than by the sound speed. This is because, under these conditions, the perturbations are nearly adiabatic, and are (at least) linearly stable. The fact that even large rms Mach numbers in runs with $`l=3`$ pc do not induce condensations indicates that the stability of these perturbations goes well beyond the first order.

## 5 Discussion

### 5.1 Comparison with previous work

Some of the issues considered in this paper have been discussed by previous authors in the context of the fragmentation of proto-galactic and proto-globular cluster clouds (e.g., Murray & Lin 1989, 1990), cooling flows in galaxy clusters (e.g., David et al. 1988; Brinkmann et al. 1990; Malagoli et al. 1990; Kang et al. 1990) and the ISM (e.g., Hennebelle & Pérault 1999, 2000; Burkert & Lin 2000; Koyama & Inutsuka 2000). In this section we put our results in perspective with respect to those works. In Paper I we have presented a comparison between our ISM-oriented scenario and previous work considering quasi-static development of TI.

It should be first pointed out that some of those studies have solved the chemical rate equations, allowing for out-of-chemical-equilibrium abundances of the relevant coolants (e.g., Murray & Lin 1989, 1990; David et al. 1988; Kang et al. 1990; Koyama & Inutsuka 2000). This is crucial in situations in which considering out-of-chemical-equilibrium cooling rates can make the difference between stability and instability. However, it is not so crucial in relation to the present-day ISM, which is already generally agreed to be thermally unstable in the range 300 K $`T`$ 6000 K under chemical equilibrium conditions. Once the medium is known to be unstable, it seems reasonable to think that variations in the cooling rates due to deviations from chemical equilibrium should only cause modifications in the details of the condensation process, but that the main qualitative behavior should not change significantly. For example, the development of supersonic motions by the development of TI similar to our results in §4.1 in the context of cooling flows and of proto-globular cluster clouds has been observed by David et al. (1988) and Kang et al. (2000), respectively. Moreover, Murray & Lin (1989, 1990) discussed the effects of various types of perturbations, cases with and without background heating, and the need for a minimum cloud collision velocity (equivalent to a compression) in order to trigger TI in proto-globular cluster clouds subject to external heating. The latter result is similar to that of the more recent study by Hennebelle & Pérault (1999) in the context of the ISM, even though these authors did not explicitly solve for the chemical rates. Thus, similar phenomenology is found with or without non-equilibrium chemistry, once the system is unstable. In any case, the caveat of our neglect of non-equilibrium chemistry and rates should be kept in mind.

Now let us specifically consider recent studies of TI in the context of the ISM. First, as mentioned above, the studies of Hennebelle & Pérault (1999) and Koyama & Inutsuka (2000) were again concerned with the triggering of TI in the stable phases of the ISM by means of strong shocks or compressions. These studies are related to the issues we discuss in §§3.3 and 4.3, although in our case the medium is already isobarically unstable from the start. Interestingly, however, the case of velocity perturbations with $`\eta 1`$ is somewhat resemblant of the case of an initially stable medium, in the sense that a strong enough compression is needed in order to raise the cooling rate to high enough values that condensation can occur. Also, it is possible that the induction of TI in the cold neutral medium due to molecular cooling found by Koyama & Inutsuka could occur during the transient overshoot we found in the condensation of large-scale perturbations if we included such cooling, but in any case this phenomenon would only support our suggestion that the static pressure-equilibrium configuration is hard to achieve, since the diffuse clouds could be subject to a further instability that would transform them all the way to small molecular cloudlets.

On the other hand, Burkert & Lin (2000) did aim their study at the formation of clumpy clouds in an unstable ISM, but they considered a different physical situation than the one we have considered here, a fact which caused them to reach significantly different conclusions from us. Specifically, they considered a case with no background heating, a choice that changes the instability criterion and, more importantly, implies that the whole flow is cooling, so that the development of TI consists of a perturbation cooling faster than the background medium. Thus, perturbations that cool isochorically and grow slowly (i.e., large-scale ones) cannot get “too far ahead” of the background gas before a minimum temperature is reached at which the gas either ceases to cool or exits the unstable range. As a consequence, Burkert & Lin concluded, similarly to Murray & Lin (1990) in the context of protogalactic clouds, that such perturbations cannot reach very large density enhancements over the background, and that only isobarically cooling ones do. This is in sharp contrast to our conclusion that it is precisely the large-scale fluctuations that reach the highest amplitudes (albeit only transiently). The reason is that in our case the background is not cooling, and thus there is no competition between it and the fluctuations to reach the minimum temperature. Which case applies at a given location in the real ISM should depend on the presence or absence of a population of heating sources that maintain a relatively constant background heating rate. It is in general agreed that the ISM is permeated by a roughly uniform UV field and cosmic rays that provide a background heating, although non-negligible variations are certainly expected to occur (Hollenbach, Parravano & McKee 2001).

### 5.2 Adequacy of the one-dimensional approach and two-dimensional simulations

The one-dimensional numerical approach that we have used here clearly has some limitations. It was chosen in order to fully resolve the condensations even at the time of maximum compression, but it has the disadvantage of being unable to capture effects that involve solenoidal (vortical) motions in the flow. In particular, it has been suggested that dynamical instabilities such as Rayleigh-Taylor (RT), Ritchmyer-Meshkov (RM), or Kelvin-Helmholtz (KH) may destroy the clouds during or after their formation in stratified media (e.g., Malagoli et al. 1990; Hattori & Habe 1990; Reale et al. 1991; Murray et al. 1993; but see Vietri et al. 1997 for an opposite view). However, within the scope of the present paper, in which we have chosen to study the development of TI in isolation, we are not concerned about the future fate of the clouds formed, but only about issues of the cloud formation process within the idealized scenario that TI is responsible for this process. This amounts to studying the limiting, most favorable case in which the classic scenario of an equilibrium ISM, with thermal-pressure bounded clouds and sharp phase transition could arise, and our main conclusion is that even within this scenario such a picture of the ISM appears unlikely.

Perhaps the closest point of contact of the issues discussed in this paper with the possibility of cloud destruction by dynamical instabilities occurs in the context of our discussion of the evolution of velocity perturbations. As we have seen, these may avoid condensation for arbitrarily long times if they have moderately large amplitudes and occur at small enough scales. Thus, in a sense, the velocity fluctuations even in our simple 1D case have the effect of counteracting the development of the isobaric mode of TI analogous to that of the dynamical instabilities discussed by other authors.

It is thus of interest to compare the effects of those dynamical instabilities on the condensed clouds with the effects of the random forcing we have considered here. Concerning the KH instability, it is expected to occur as the condensations begin to sink in a stratified medium as they develop inverse buoyancy (Balbus & Soker 1989). Thus, it could only start to act after the condensation has reached a sufficiently advanced evolutionary stage. That is, its effect is one of destroying a cloud after it has formed, rather than obstructing its growth, as we have found here.

On the other hand, the RT instability might indeed occur while the condensation proceeds, as the flow is accelerated inwards of the condensation, with the density gradient pointing also inwards. To test this possibility, we have performed two two-dimensional simulations analogous to run DEN75, but including weak random initial velocity perturbations with rms Mach numbers 0.1 and 0.01 with a half-wavelength of 1/8 of the box. We find no traces of development of RT. Instead, the velocity fluctuations generate density fluctuations that evolve in a very similar way to the simulations of §4.2, with the small-scale fluctuations condensing first, and then merging towards the center of the large-scale perturbation. The certainty that RT does not arise comes from the fact that the velocity field converges onto the fragments rather than shear between them, as would be expected for RT. Thus, we conjecture that TI development inhibits that of RT.

Finally, the RM instability probably does not apply to single clouds formed by TI because it occurs upon the interaction of a shock with a discontinuity. Instead, what we have found in this paper is that the evolution of a single perturbation generates a shock wave that propagates outwards from the condensation. It may be relevant, however, in the case of multiple fluctuations, as in §4.2, in which the shock propagating out of one condensation may interact with another. However, this again refers to the future evolution of the condensations formed by TI, which is out of the scope of this paper.

A final note is that, in any case, the effects of all these instabilities go in the same direction as the issues discussed in this paper: as they tend to destroy the condensations, they tend to restore a more continuum-like ISM, rather than one separated in discrete phases.

### 5.3 Implications

The results of this paper suggest a number of implications. All of them require further investigation, but the possibility of their realization in the ISM is of great interest:

1. The outcome of the development of TI, even if it were the sole cloud-forming agent in the atomic ISM, is not necessarily the production of quiescent thermal pressure-bounded clouds in contact with the warm stable phase. Large-scale density perturbations require long times to condense (over 7 Myr for perturbations of size $`>10`$ pc), during which they are in a highly dynamic state, far from pressure equilibrium. At later stages, pressure equilibrium is approximately restored, but at the expense of disrupting thermal equilibrium of the surrounding medium through the shock propagating away from the condensation. This surrounding gas continues to accrete onto the condensation for much longer times, providing a ram presure that keeps the cloud density above the pressure-equilibrium one by factors of up to 2. Moreover, the accreting gas is not in thermal equilibrium, but in a nearly isobaric regime, in which the (small) pressure fluctuations have the same sign as those of the density, and therefore the medium has no strong tendency to fragment any further. In particular, this non-equilibrium situation implies that the unstable density range as shown in Fig. 1 does not exactly coincide with the unstable temperature range.

Smaller-scale (large-$`\eta `$) fluctuations ($`3`$ pc at $`\rho 1`$ cm-3, or even larger at proportionally lower densities) can condense and reach near thermal-pressure balance in shorter times ($`4`$ Myr) if they are nearly isobaric in nature, but they are stable if they are nearly adiabatic, a condition that can be accomplished if they originate from velocity fluctuations with sufficiently large $`_{\mathrm{rms}}`$ at sufficiently small scales (§4.3). Moreover, small-scale fluctuations, if part of a full spectrum of fluctuations, may continue to form part of the condensation of large-scale ones, if the relative amplitude of the former is small enough and they occur at sufficiently large $`\eta `$4.2). In general, it appears that the pure nonlinear development of TI, even in the absence of other agents like magnetic fields, rotation, stellar energy injection, etc., does not necessarily lead to rapid formation of clouds with sharp phase transition at their boundaries.

2. The existence of accreting gas in the “unstable” range at late stages of the evolution, and the stabilization produced by velocity fluctuations in a certain range of scales suggests that in the turbulent ISM the “unstable” phase may actually be significantly populated, since the density fluctuations are in general expected to have a dynamical origin. This may explain observational (e.g., Dickey, Salpeter & Terzian 1977; Kalberla et al. 1985, Spitzer & Fitzpatrick 1995, Fitzpatrick & Spitzer 1997; Heiles 2001) and recent numerical results by Gazol et al. (2001) and Kritsuk & Norman (2002), suggesting that relatively large amounts of gas are present in the ISM in the unstable temperature range.

3. In summary, both the gas accreting onto the condensations at late times and the “unstable” gas in the presence of stable adiabatic fluctuations, are out of thermal equilibrium. In general, it can be said that thermal and pressure equilibria are mutually exclusive for gas in the unstable range, because each implies a different value of the pressure at densities intermediate between $`\rho _0`$ and $`\rho _{\mathrm{isob}}`$ (fig. 3). This incompatibility between the two equilibria can only be avoided by means of sharp phase segregation, i.e., when no gas in the unstable exists. However, if such gas exists, it should in general be expected to be out of thermal equilibrium. This suggests a possible observational test for determining whether the gas apparently seen at unstable temperatures corresponds to this kind of regime, or at least whether it is in thermal equilibrum or not. This could be done by either a) simultaneously determing two of its thermodynamic variables, or b) comparing directly observed cooling rates (e.g., fine structure lines) with theoretical estimates of the heating rate (e.g., photoelectric heating) in specific regions (C. Heiles, J. Scalo, private communications). If this were confirmed, it would provide strong evidence that turbulent motions populate all regions of the thermodynamic variable space, preventing a sharp segregation of the atomic ISM into the stable phases of TI.

4. The nearly adiabatic response of the gas to velocity perturbations in cases with relatively large values of $`\eta `$, suggests that the warm diffuse ISM, in which this condition is satisfied with respect to perturbations of up to several pc (§4.1), may exhibit a very weakly compressible behavior, and thus be close to a Kolmogorov regime, which applies to incompressible turbulent flows. This may explain why observations of the diffuse gas tend to find density fluctuation spectra with the signature of a Kolmogorov regime (e.g., Minter & Spangler 1996; Stanimirovic & Lazarian 2001; Dickey et al. 2001). The numerical results by Kritsuk & Norman (2002) also suggest a nearly adiabatic regime for the diffuse gas.

## 6 Summary

In this paper we have performed a detailed numerical investigation of the nonlinear development of the isobaric mode of TI in the context of the atomic ISM, and of its interplay with a spectrum of fluctuations in density or velocity, aimed at establishing the feasibility of the classical scenario of a thermal- and pressure-equilibrium ISM, in which the clouds are directly confined by the more dilute, stable warm gas, as suggested in the classical two- and three-phase models of the ISM (Field et al. 1969; McKee & Ostriker 1977), even under the idealization that TI and a spectrum of density or velocity fluctuations are the only processes at play. Most of the discussion is based on one parameter, $`\eta `$, the ratio of the cooling to the dynamical times.

We have found that the evolution for initial Gaussian perturbations of FWHM smaller than $`15`$ pc is approximately isobaric. In that case, a slow-growth stage characterized by quasi-stationary accretion of low-density gas is reached in $`4`$$`6`$ Myr for initial amplitudes of $`10`$$`20`$ $`\%`$ over the mean density, although infall velocities $`0.7`$ km s-1 still exist at those times. The infall velocities remain subsonic throughout the evolution for small perturbation sizes, but infall velocities with Mach number $`>0.5`$ are acquired if initial perturbations are large enough ($`15`$ pc). For the latter type of perturbations, a shock wave, formed during the “crushing” and overshooting stage, propagates outwards from the cloud and throws the surrounding medium out of thermal equilibrium, establishing near pressure balance. For the case with a Gaussian density perturbation with FWHM$`=75`$ pc, the final mass density of the cloud is more than twice $`\rho _{\mathrm{isob}}`$ due to the ram pressure of the infalling gas. The time to reach the slow-accretion stage ranges from 4 to 30 Myr for perturbations with initial FWHM of 3 to 75 pc. These times are long enough that a major disturbance is likely to pass through the condensation.

We have quantified the competition between scales by starting with a power-law spectrum of initial density perturbations, to assess the possibility that the smaller scales, which have the fastest growth rates, could consume most of the mass available before large-scale perturbations have time to grow, thus establishing a “forest” of thermal-pressure bounded small clouds in relatively short times. We have found that small-scale perturbations may or may not use up most of the mass available depending on both the spectrum slope $`\delta `$ and the initial values of $`\eta `$ for the smallest-scale perturbations. We have found that the inclusion of large wavenumbers has a progressively smaller effect on the final state as $`\delta `$ increases and as the perturbations are shifted to smaller scales, increasing $`\eta `$. For example, for $`\delta =3/2`$, a value of $`\eta (k_{\mathrm{min}})>0.4`$ is needed in order for the final states to be roughly equivalent, while for $`\delta =1/2`$ this requires $`\eta (k_{\mathrm{min}})>4.0`$. We looked at this problem to determine whether the dynamical character of the condensation of large-scale perturbations might be preserved even in the presence of small-scale ones. Although our simulations indicate that the evolution is also very dynamical even if small-scale perturbations are added, more work is needed to be conclusive concerning this result.

However, we have found that a much more interesting effect appears in the case of a medium continuously stirred by an external force. We have argued that velocity fluctuations are necessary to produce realistic density fluctuations in a continuum, their evolution thus ranging from a nearly adiabatic to a thermal equilibrium response as $`\eta `$ is varied from $`1`$ to $`1`$, provided that in this case $`\eta `$ is determined by the turbulent, rather than the sound crossing time. When the perturbations evolve nearly adiabatically, they cannot condense, because they re-expand before they can cool – i.e., they are stable. We have confirmed this suggestion numerically, showing that flows with $`\rho =1`$ cm-3 forced at scales such that $`\eta >1`$ ($`\lambda 0.3`$$`1`$ pc) and with $`_{\mathrm{rms}}0.3`$ do not develop any condensations over the entire duration of the simulations (between 12 and 26 Myr). However, flows forced at scales such that $`\eta 1`$ ($`\lambda l_{\mathrm{eq}}10`$ pc) produce condensations regardless of the forcing strength, because in this case either the effective waves have an imaginary sound speed ($`\eta 1`$), and trigger the condensation mode, which is unstable, or the density enhancement they induce causes $`\eta `$ to become smaller than unity by increasing the cooling rate ($`\eta 1`$).

As implications, we have suggested that:

1. The existence of accreting gas in the unstable density range and the stabilization of the flow by small-scale velocity fluctuations may be at the origin of the relatively large amounts of gas mass in the unstable regime found in both observations (Dickey et al. 1978; Kalberla et al. 1985, Spitzer & Fitzpatrick 1995, Fitzpatrick & Spitzer 1997; Heiles 2001) and simulations (Gazol et al. 2001; Kritsuk & Norman 2002) of the ISM.

2. Since transonic velocity fluctuations of scales up to several parsecs in the warm diffuse ISM have $`\eta 1`$, this medium should respond nearly adiabatically to them, and the flow regime should then be only weakly compressible. This may explain recent observational results suggesting that its kinetic energy spectrum is very close to the Kolmogorov one (Minter & Spangler 1996; Dickey et al. 2001; Stanimirovic & Lazarian 2001), since the latter is expected to apply to incompressible or weakly compressible cases only.

3. Simultaneous observational determinations of two thermodynamic variables should allow to distinguish whether the gas apparently seen in the unstable temperature range is in thermal equilibrium or not, thus helping to decide whether equilibrium models of the ISM are applicable or not.

We are happy to acknowledge stimulating exchanges with A. Brandenburg, C. Heiles, P. Hennebelle, A. Raga, J. Scalo and E. Zweibel. A deep and thorough referee report prompted much improvement of the paper. This work has received financial support from CONACYT grants 27752-E to E. V.-S., I32888-E to A. G., 32139-E to Luc Binette, and a CONACYT postdoctoral grant to F.J.S.S.