# Influence of Alfvén waves on thermal instability in the interstellar medium

P. Hennebelle, T. Passot

► To cite this version:

P. Hennebelle, T. Passot. Influence of Alfvén waves on thermal instability in the interstellar medium. Astronomy & Astrophysics - A&A, 2006, 448, pp.1083-1093. 10.1051/0004-6361:20053510. hal-00388614

HAL Id: hal-00388614

https://hal.science/hal-00388614v1

Submitted on 17 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Influence of Alfvén waves on thermal instability in the interstellar medium

P. Hennebelle<fnref n="1" /> and T. Passot<fnref n="2" />

Received 25 May 2005 / Accepted 14 November 2005

## ABSTRACT

The effect of Alfvén waves on the thermal instability of the Interstellar Medium (ISM) is investigated both analytically and numerically. A stability analysis of a finite amplitude circularly polarized Alfvén wave propagating parallel to an ambient magnetic field in a thermally unstable gas at thermal equilibrium is performed, leading to a dispersion relation that depends on 3 parameters, namely the square ratio of the sonic and Alfvén velocities (β), the wave amplitude and the ratio between the wave temporal period and the cooling time. Depending on the values of these 3 parameters, the Alfvén waves can stabilize the large-scale perturbations, destabilize those whose wavelength is a few times the Alfvén wavelength λAW, or leave the growth rate of the short scales unchanged. To investigate the non-linear regime, two different numerical experiments are performed in a slab geometry. The first one deals with the development of an initial density perturbation in a thermally unstable gas in the presence of Alfvén waves. The second one addresses the influence of those waves on the thermal transition induced by a converging flow. The numerical results confirm the trends inferred from the analytic calculations, i.e. the waves prevent the instability at scales larger than λAW and trigger the growth of wavelengths close to λAW, therefore producing a very fragmented cold phase. The second numerical experiments shows that i) the magnetic pressure prevents the merging of the CNM fragments therefore maintaining the complex structure of the flow and organizing it into groups of clouds ii) these groups of CNM clouds have an Alfvénic internal velocity dispersion iii) strong density fluctuations (≈10ρ_cnm) triggered by magnetic compression occur. We note that during this event there is no stiff variation of the longitudinal velocity field. This is unlike the hydrodynamical case for which the clouds are uniform and do not contain significant internal motions except after cloud collisions. In this situation a strong density fluctuation occurs, accompanied by a stationary velocity gradient through the cloud.

**Key words.** ISM: instabilities – magnetohydrodynamics – turbulence – ISM: clouds – ISM: magnetic fields

## 1. Introduction

It is now well established that the atomic interstellar medium is a thermally bistable gas which at a pressure of about 4000 K cm⁻³ can be in two different phases, namely the Warm Neutral Medium (WNM) and the Cold Neutral Medium (CNM), roughly in pressure equilibrium (Kulkarni & Heiles 1987; Field et al. 1969; Wolfire et al. 1995, 2003). It is also believed that the atomic gas is strongly magnetized with a magnetic intensity around 5 µG (Troland & Heiles 1986; Heiles 1987; Heiles & Troland 2005). Although the structure of the magnetic field is poorly known, some observational evidences seem to indicate that the fluctuating part is at least comparable to the uniform one suggesting that magnetized waves may be of great importance for the dynamics of the ISM. Another important observational result is the absence of correlation between the magnetic intensity and the density in the interstellar atomic gas. Various works have investigated the effect of MHD waves on the dynamics of a polytropic or nearly polytropic (non

thermally unstable) magnetized gas (Dewar 1970; Goldstein 1978; McKee & Zweibel 1993; Passot et al. 1995; Falle & Hartquist 2002; Passot & Vázquez-Semadeni 2003).

It is thus of great interest to investigate the simultaneous role of magnetic fields and of the bistable nature of the flow on the physics of the ISM. Field (1965) considers the effect of a transverse magnetic field on the thermal instability and generalizes the isobaric criterion. In the context of cooling flows, Loewenstein (1990) studies the thermal instability in the presence of a static field. Hennebelle & Pérault (2000) investigates the role of an initially uniform magnetic field, analytically and numerically, when the thermal condensation is dynamically triggered. They propose a mechanism based on magnetic tension to explain the thermal collapse in a magnetized flow and argue that the magnetic intensity in the WNM and in the CNM should not be very different. Piontek & Ostriker (2004, 2005) study the magneto-rotational instability in a thermally unstable gas.

<footnote n="1">¹ Laboratoire de Radioastronomie Millimétrique, UMR 8112 du CNRS, École Normale Supérieure et Observatoire de Paris, 24 rue Lhomond, 75231 Paris Cedex 05, France
e-mail: patrick.hennebelle@ens.fr</footnote>

<footnote n="2">² CNRS, Observatoire de la Côte d'Azur, BP 4229, 06304 Nice Cedex 4, France
e-mail: passot@obs-nice.fr</footnote>
Here we investigate the effect of a finite amplitude circularly polarized Alfvén wave on the thermal instability, both analytically and numerically. These Alfvén waves, which are exact solutions of the MHD equations, are non dissipative and are thus very likely to be present in a magnetized gas such as the ISM. Similar waves, although at a smaller scale and in a regime affected by dispersive effects, exist in the solar wind upstream of the earth’s bow shock (Spangler et al. 1988). In-situ satellite observations have clearly identified circularly polarized quasi-monochromatic Alfvén wave packets, probably generated by reflected protons.

In Sect. 2 we present a stability analysis of a circularly polarized Alfvén wave propagating in a thermally unstable gas at thermal equilibrium and obtain a dispersion relation that generalizes the relation obtained by Field (1965). We then solve this equation numerically for various regimes and discuss the consequences for the ISM. In Sect. 3, we perform various numerical experiments to confirm the analytic prediction and to investigate the non-linear regime. Section 4 concludes the paper.

## 2. Analysis

The equations governing, in the magnetohydrodynamic (MHD) limit, the one-dimensional motion of a plasma permeated by a uniform magnetic field $B_0$ in a slab geometry read

$$ \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} = 0 \qquad (1) $$

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = - \frac{1}{\rho} \frac{\partial}{\partial x} \left( P + \frac{|b|^2}{8\pi} \right) \qquad (2) $$

$$ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} = \frac{B_x}{4\pi\rho} \frac{\partial b}{\partial x} \qquad (3) $$

$$ \frac{\partial b}{\partial t} + \frac{\partial}{\partial x}(ub) = B_x \frac{\partial v}{\partial x}, \qquad (4) $$

$$ \left( \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} \right) = -(\gamma - 1) T \frac{\partial u}{\partial x} - \frac{1}{C_v} \mathcal{L} + \frac{1}{\rho C_v} \partial_x (\kappa(T) \partial_x T) \qquad (5) $$

where all fields only depend on the coordinate $x$ and time $t$.

The ambient field $B_0$ is assumed to be oriented in the $x$ direction. The velocity field $V$ has a component $u = V_x$ along the $x$ coordinate and two transverse components combined in the complex number $v = V_y + iV_z$. Similarly, we write the magnetic field as $b = B_y + iB_z$, the component $B_x = B_0$ remaining constant. The mass density and thermal pressure are denoted by $\rho$ and $P$ respectively and we assume a perfect gas law $P = \frac{\rho_0 T}{\mu}$ where $R = \frac{k_B}{m_H}$ is the universal gas constant, $k_B$ the Boltzmann constant and $\mu$ the mean molecular weight in units of the hydrogen mass $m_H$. Heating and cooling processes are combined in a single net cooling function $\mathcal{L}$. The parameter $\gamma$ denotes the ratio of the specific heats at constant pressure $C_p$ and at constant volume $C_v$. Note that $C_v = \frac{R}{(\gamma-1)\mu}$. Thermal diffusivity due to neutrals is the dominant one. It is isotropic in spite of the presence of the magnetic field and equal to $\kappa(T) = 5/3C_v\eta(T)$, where $\eta(T) = 5.7 \times 10^{-5}(T/1\,\text{K})^{1/2} \, \text{g cm}^{-1} \, \text{s}^{-1}$ (Lang 1974).

Throughout most of Sect. 2, thermal diffusivity will be neglected for simplicity, leading to the assumption of a vanishingly small Field length. Its effect will briefly be addressed in

Sect. 2.5 and it is explicitly taken into account in the numerical simulations of Sect. 3.

Two main non-dimensional numbers can be defined, the sonic Mach number $M_s = V_0/c_s$ ratio of a typical velocity $V_0$ with the constant sound speed $c_s = \sqrt{\frac{\gamma k_B T}{\mu m_H}}$ and the Alfvénic Mach number $M_a = V_0/c_a$, where $c_a = B_x/(4\pi\rho_0)^{1/2}$ is the Alfvén speed of the unperturbed system. The plasma beta is here defined by $\beta = \frac{M_a^2}{M_s^2}$.

These equations have exact solutions in the form of circularly polarized plane Alfvén waves of arbitrary amplitude (Ferraro 1955). They read $b_0 = B_\perp \exp[-i\sigma(k_0x - \omega_0t)]$ with constant density $\rho_0$ and temperature $T_0$, and zero longitudinal velocity $u_0 = 0$. The transverse wave velocity is related to the magnetic field perturbation by $v_0 = V_\perp \exp[-i\sigma(k_0x - \omega_0t)] = -\frac{b_0}{\sqrt{4\pi\rho_0}}$. The polarization of the wave is determined by the parameter $\sigma$, with $\sigma = +1$ ($\sigma = -1$) for a right-handed (left-handed) wave. In the absence of dispersive effects we can, without restriction, take $\sigma = -1$. Moreover, $B_\perp$ can be taken real. The dispersion relation reads $\omega_0^2 = c_a^2 k_0^2$.

### 2.1. Ambipolar diffusion and wave steepening

Before considering the idealized configuration of single fluid perfectly coupled to the magnetic field with a circularly polarized monochromatic Alfvén wave train of infinite length, it is important to estimate the timescale of two important processes which are not considered in this study. The first one is the ambipolar diffusion that takes place in weakly ionized gas and that induces energy dissipation due to the friction between the neutrals and the ions. The second one is the steepening of the MHD waves that takes place when the wave polarization is not perfectly circular or when the wave packet is modulated.

#### 2.1.1. Timescale of the ambipolar diffusion

The ambipolar diffusion time is given by (Shu 1992):

$$ t_{\text{da}} = \frac{4\pi\gamma_{\text{da}}\rho_{\text{n}}\rho_{\text{i}}L^2}{B^2}, \qquad (6) $$

where $\gamma_{\text{da}}$ is the friction coefficient between ions of density $\rho_i$ and neutrals of density $\rho_n$, $L$ is the typical scale to be considered and $B$ is the magnetic intensity.

In the case of a molecular cloud, it has been estimated to $3 \times 10^{13} \, \text{cm}^3 \, \text{g}^{-1} \, \text{s}^{-1}$ (Shu 1992). In the case of the atomic gas, since the mass of the neutral and the ions are two times smaller, $\gamma_{\text{da}}$ is about three times larger, leading to

$$ t_{\text{da}} \approx 5 \times 10^{15} \left(\frac{5 \mu \text{G}}{B}\right)^2 \frac{\xi_i}{10^{-3}} \left(\frac{n_n}{100 \, \text{cm}^{-3}}\right)^2 \left(\frac{L}{10^{18} \, \text{cm}}\right)^2 \text{s}, \quad (7) $$

where $\xi_i$ is the gas ionization and $n_n$ the particle density of the neutrals.

In the CNM the ionization is about $4 \times 10^{-4}$, therefore at the scale of say 0.1 pc, the ambipolar diffusion time is about $2 \times 10^{14} \, \text{s}$ or about 10 Myr. In the WNM, the density is about $1 \, \text{cm}^{-3}$, the ionization about 0.1, therefore at a scale of say 1 pc,
the ambipolar diffusion time is $5 \times 10^{14}$ s. In both cases, am-bipolar diffusion operates on a rather long characteristic time and can thus safely be neglected.

### 2.1.2. Timescale of the wave steepening

The equation governing the dynamics of parallel propagating Alfvén waves in the long wavelength, small amplitude limit was derived by Cohen & Kuslrud (1974). A similar equation was derived previously by Rogister (1971) from the Vlasov-Maxwell system in the context of a collisionless plasma, taking into account kinetic effects such as dispersion and Landau damping. The main point is that the evolution of the magnetic field components perpendicular to the ambient field involves a nonlinearity of the form $\alpha \partial_x (|b|^2 b)$, where $\alpha$ is a coefficient depending on $\beta$. From this formula, one sees in particular that circularly polarized Alfvén waves, for which $|b|$ is constant, do not steepen. In the case where the wave amplitude is modulated, the characteristic time of steepening is given by (using Eq. (40) in Cohen & Kuslrud 1974)

$$ \tau = \frac{4(1-\beta)B_0^2}{3c_A \partial_x |b|^2} \quad (8) $$

Denoting by $\lambda$ the wavelength, by $l_c$ the typical length of the amplitude modulation and by $\omega_s$ the characteristic growth rate associated with wave steepening, one gets

$$ \frac{\omega_s}{\omega_0} = \frac{3}{4(1-\beta)} \frac{\lambda}{l_c} \frac{b_0^2}{B_0^2} \quad (9) $$

The Cohen-Kuslrud equation is not valid for $\beta$ close to unity, a point where sound waves have the same phase speed as Alfvén waves and have thus to be retained explicitly in the nonlinear dynamics. In general however, we see that the ratio $\omega_s/\omega_0$ depends quadratically on the wave amplitude and is inversely proportional to the normalized scale of amplitude modulation. As an example, for a wave train of amplitude $b_0/B_0 = 0.5$ and coherence length $l_c/\lambda = 5$, with $\beta = 0.5$, one gets $\omega_s/\omega_0 = 0.075$, a value about five times smaller than the typical growth rates shown in Fig. 1. This process of wave steepening thus appears subdominant with respect to the combined thermal and decay instabilities considered in this paper. Moreover, as will be demonstrated by the numerical simulations, wave steepening does indeed occur but it does not affect our conclusions.

Another process is at play in real three-dimensional situations, namely turbulent cascade, which occurs mainly in transverse directions as a result of nonlinear interaction between counter-propagating wave packets. The nonlinear eddy turnover time is expected to be smaller than the other processes mentioned above for high enough amplitude and situations where Alfvén waves propagate in both directions in roughly equal amounts. When waves propagate in a privileged direction, as in the solar wind, the turbulent cascade is much less efficient, leaving enough time for the other instabilities to develop. The interaction between waves and turbulence in compressible MHD is still a very debated topic (see e.g. Cho & Lazarian 2003). Three-dimensional simulations will have to be performed to address the competition between these phenomena.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Growth rate as a function of the wave number. The first panel is for the adiabatic case, second, third and fourth panels display results for various values of $\omega_c/\omega_0$. The solid line corresponds to the case $A=0$, dotted line to $\beta=0.9$, $A=0.5$, short dashed line to $\beta=0.5$, $A=0.5$, dot-dashed line to $\beta=0.5$, $A=1$, double dot-dashed line to $\beta=0.1$, $A=0.5$ and long dashed line to $\beta=0.1$, $A=1$.</figcaption></figure>

### 2.2. Derivation of the dispersion relation

We shall now consider perturbations about the exact solution described above, namely a circularly polarized Alfvén wave propagating in a uniform gas at thermal equilibrium and determine the dispersion relation for normal modes. Such an analysis has been performed by Goldstein (1978) for an isothermal gas and by Lou (1996) for a self-gravitating gas. Let us write the perturbations as

$$ \rho = \rho_0 + \delta\rho, \quad T = T_0 + \delta T, \quad v = v_0 + \delta v, \quad b = b_0 + \delta b, \quad u = \delta u. $$

The linearized equations are, denoting by a star the complex conjugate,

$$ \partial_t \delta\rho + \rho_0 \partial_x \delta u = 0, $$

$$ \rho_0 \partial_t \delta u = - \partial_x \delta P - \frac{1}{8\pi} \partial_x (b_0 \delta b^* + b_0^* \delta b), $$

$$ \rho_0 \partial_t \delta v + \delta \rho \partial_t v_0 + \rho_0 \delta u \partial_x v_0 = \frac{1}{4\pi} B_x \partial_x \delta b, \quad (10) $$

$$ \partial_t \delta b + \delta u \partial_x b_0 + b_0 \partial_x \delta u - B_x \partial_x \delta v = 0, $$

$$ C_v (\partial_t \delta T + (\gamma - 1) T_0 \partial_x \delta u) = - \left( \left( \frac{\partial \mathcal{L}}{\partial \rho} \right)_T \delta \rho + \left( \frac{\partial \mathcal{L}}{\partial T} \right)_\rho \delta T \right). $$
Taking perturbations in the form

$$
\begin{align}
\delta\rho &= \delta\bar{\rho} \exp(-i\omega t + ikx) + \text{c.c.}, \nonumber \\
\delta T &= \delta\bar{T} \exp(-i\omega t + ikx) + \text{c.c.}, \nonumber \\
\delta u &= \delta\bar{u} \exp(-i\omega t + ikx) + \text{c.c.}, \nonumber \\
\delta v &= \delta\bar{v}^+ \exp(-i(\omega + \omega_0)t + i(k + k_0)x) \nonumber \\
&\quad + (\delta\bar{v}^-)^* \exp(i(\omega - \omega_0)t - i(k - k_0)x), \tag{11} \\
\delta b &= \delta\bar{b}^+ \exp(-i(\omega + \omega_0)t + i(k + k_0)x) \nonumber \\
&\quad + (\delta\bar{b}^-)^* \exp(i(\omega - \omega_0)t - i(k - k_0)x), \nonumber
\end{align}
$$

we obtain the following equations

$$
\begin{align}
& -\omega \hat{\delta \rho} + k \rho_0 \delta \hat{u} = 0, \nonumber \\
& -\rho_0(\omega + \omega_0) \hat{\delta v}^+ + (-\omega_0 \hat{\delta \rho} + k_0 \rho_0 \delta \hat{u}) V_\perp = \frac{B_x}{4\pi} (k + k_0) \hat{\delta b}^+, \nonumber \\
& -\rho_0(\omega - \omega_0) \hat{\delta v}^- + (\omega_0 \hat{\delta \rho} - k_0 \rho_0 \delta \hat{u}) V_\perp = \frac{B_x}{4\pi} (k - k_0) \hat{\delta b}^-, \nonumber \\
& -\omega \rho_0 \delta \hat{u} + k \frac{k_B}{\mu m_H} (\rho_0 \hat{\delta T} + T_0 \delta \hat{\rho}) + \frac{B_\perp}{8\pi} k (\hat{\delta b}^+ + \hat{\delta b}^-) = 0, \tag{12} \\
& -(\omega + \omega_0) \hat{\delta b}^+ + B_\perp (k + k_0) \delta \hat{u} - B_x (k + k_0) \delta \hat{v}^+ = 0, \nonumber \\
& -(\omega - \omega_0) \hat{\delta b}^- + B_\perp (k - k_0) \delta \hat{u} - B_x (k - k_0) \delta \hat{v}^- = 0, \nonumber \\
& i C_v (-\omega \hat{\delta T} + (\gamma - 1) T_0 k \delta \hat{u}) = - (\partial_\rho \mathcal{L} \delta \hat{\rho} + \partial_\tau \mathcal{L} \delta \hat{T}). \nonumber
\end{align}
$$

The fluctuations $\hat{\delta T}$, $\hat{\delta \rho}$, and $\hat{\delta b}^{\pm}$ can easily be obtained as a function of $\hat{\delta u}$ in the form

$$
\begin{align}
\omega \delta \rho &= k \rho_0 \delta u, \nonumber \\
\delta T &= \frac{\delta u}{-\partial_T \mathcal{L} + i C_v \omega} \left( i k C_v (\gamma - 1) T_0 + \frac{k}{\omega} \rho_0 \partial_\rho \mathcal{L} \right), \tag{13} \\
\delta b^{\pm} &= \frac{\delta u}{\omega_{\pm}^2 - c_a^2 k_{\pm}^2} B_{\perp} k_{\pm} \frac{\omega_0^2 k}{\omega k_0} \left( \frac{\omega_{\pm}}{\omega_0} \frac{\omega}{\omega_0} \frac{k_0}{k} \pm \frac{\omega}{\omega_0} \frac{k_0}{k} - \pm 1 \right). \nonumber
\end{align}
$$

One then gets the dispersion relation

$$
\left( -\omega^2 + k^2 \frac{k_B}{\mu m_H} T_0 + k^2 \frac{k_B}{\mu m_H} \left( \frac{i C_v (\gamma - 1) T_0 \omega + \rho_0 \partial_\rho \mathcal{L}}{-\partial_T \mathcal{L} + i C_v \omega} \right) \right) \\
\times \left( \frac{\omega}{\omega_0} - \frac{k}{k_0} \right) \left( \left( \frac{\omega}{\omega_0} + \frac{k}{k_0} \right)^2 - 4 \right) = \\
- \omega_0^2 \left( \frac{B_\perp}{B_x} \right)^2 \left( \frac{k}{k_0} \right)^2 \left( \left( \frac{\omega}{\omega_0} \right)^3 + \frac{k}{k_0} \left( \frac{\omega}{\omega_0} \right)^2 - 3 \frac{\omega}{\omega_0} + \frac{k}{k_0} \right), \quad (14)
$$

which can be rewritten

$$
\tilde{\omega}^2 - \beta \tilde{k}^2 \frac{\tilde{\omega} + i \tilde{\omega}_b}{\tilde{\omega} + i \tilde{\omega}_c} = \tilde{k}^2 A^2 \frac{\tilde{\omega}^3 + \tilde{k} \tilde{\omega}^2 - 3 \tilde{\omega} + \tilde{k}}{(\tilde{\omega} - \tilde{k})((\tilde{\omega} + \tilde{k})^2 - 4)}, \quad (15)
$$

where $\tilde{\omega}_c = \frac{\omega_c}{\omega_0} = \frac{\partial_T \mathcal{L}}{C_c \omega_0}$, $\tilde{\omega}_b = \frac{\omega_b}{\omega_0} = \frac{T_0 \partial_T \mathcal{L} - \rho_0 \partial_\rho \mathcal{L}}{\gamma C_c T_0 \omega_0} = \frac{1}{\gamma C_c \omega_0} \left(\frac{\partial \mathcal{L}}{\partial T}\right)_P$, $\tilde{\omega} = \omega/\omega_0$, $\tilde{k} = k/k_0$ and $A = B_\perp/B_x$.

Note that in the limit $A = 0$, Eq. (15) becomes identical to the dispersion relation obtained by Field (1965) whereas if $\mathcal{L} = 0$, it reduces to the dispersion relation obtained by Goldstein (1978).

## 2.3. Asymptotic behaviors

Before numerically solving Eq. (15), we consider various asymptotic limits.

### 2.3.1. Static magnetic field

It is possible to recover the dispersion relation in the absence of waves for a situation where the ambient field is oblique, making an angle $\theta$ with the x-axis. In this case $B_x = B_0 \cos \theta$ and one must take the limit $\omega_0 \to 0$ with $\frac{\omega_0}{k_0} = c_A = \frac{B_0}{\sqrt{4\pi\rho_0}} \cos \theta$ and $A = \tan \theta$. It follows that

$$
\omega^2 - k^2 c_s^2 \frac{\omega + i \omega_b}{\omega + i \omega_c} = \frac{k^2 \omega^2 v_A^2 \sin^2 \theta}{\omega^2 - k^2 v_A^2 \cos^2 \theta}, \quad (16)
$$

where we denote $v_A^2 = \frac{B_0^2}{4\pi\rho_0}$.

In the case $\theta = \pi/2$ and close to the threshold ($\omega \approx 0 \ll kc_s$), we have

$$
\omega \approx -i \frac{v_A^2 \omega_c + c_s^2 \omega_b}{v_A^2 + c_s^2}. \quad (17)
$$

Since the magnetic field is purely transverse, the magnetic tension vanishes and the magnetic pressure adds up to the thermal pressure, making the gas more stable. The criterion for thermal instability is simply $v_A^2 \omega_c + c_s^2 \omega_b \le 0$ which was first obtained by Field (1965). It shows that a purely transverse magnetic field can suppress the thermal instability. This is because, in this geometry, the magnetic pressure is proportional to the square of the density and therefore $\partial_\rho P_{\text{tot}}$ can be positive even if $\partial_\rho P_{\text{therm}}$ is negative.

### 2.3.2. Instability thresholds

This section addresses the neighborhood of the instability, a situation where $\Re(\omega) \to 0$.

We shall first discuss the case $\Re(\omega) = 0$, corresponding to the so-called condensation or entropy mode. In the situation where $|\tilde{\omega}_b|$ and $|\tilde{\omega}_c|$ are smaller than $\tilde{k}$, the characteristic cooling time is longer than the Alfvén crossing time and the dispersion relation reduces to

$$
\tilde{\omega} \approx -i \frac{\beta (4 - \tilde{k}^2) \tilde{\omega}_b + A^2 \tilde{\omega}_c}{\beta (4 - \tilde{k}^2) + A^2}. \quad (18)
$$

When $\beta$ is not too small, sound waves have time to restore pressure equilibrium while the gas cools and one expects that the growth rate will be close to the isobaric one. In the opposite case, where $\beta$ is very small, the growth rate should be close to the isochoric one since Alfvén waves are not accompanied by pressure or density perturbations. These conclusions are easily recovered from Eq. (18) which for very small amplitude ($A^2/\beta(4-\tilde{k}^2) \ll 1$) gives

$$
\tilde{\omega} \approx -i \left( \tilde{\omega}_b + \frac{A^2}{\beta(4-\tilde{k}^2)} (\tilde{\omega}_c - \tilde{\omega}_b) \right), \quad (19)
$$

while, in the case where $\beta(4 - \tilde{k}^2)/A^2 \ll 1$ it rewrites

$$
\tilde{\omega} \approx -i \left( \tilde{\omega}_c + \frac{\beta(4 - \tilde{k}^2)}{A^2} (\tilde{\omega}_b - \tilde{\omega}_c) \right). \quad (20)
$$

In the case where $A^2/\beta \ll 1$, more likely to be met in the ISM, the effect of the Alfvén wave depends both on the sign
of $\tilde{\omega}_c - \tilde{\omega}_b$ and of that of $\tilde{k} - 2$. In a typical region of the ISM with $T \approx 1000$ K, $\tilde{\omega}_c - \tilde{\omega}_b \ge 0$ and therefore the waves stabilize (destabilize) the gas if $\tilde{k} < 2$, (respectively $\tilde{k} \ge 2$).

In the limit $\tilde{k} \gg 1$ and/or for $A = 0$, one recovers the isobaric growth rate, the destabilizing effect of the Alfvén waves becoming asymptotically small as $\tilde{k}$ increases.

If $\tilde{\omega}_b \ll \tilde{k} \ll 1$, we find a growth rate similar to that given by Eq. (17) except for the factor 2 that divides $v_A^2$. This is due to the magnetic tension that tends to unbend the magnetic field lines, making the stabilization of the magnetic pressure less efficient.

Two different limits are obtained when $\tilde{k} \ll \tilde{\omega}_b$. When the isochoric criterion is not verified and $\tilde{\omega}/\tilde{k} \equiv \alpha$ remains finite, a situation where the growth rate vanishes with $\tilde{k}$, one gets

$$ (1 - \alpha) \left( \alpha^2 - \beta \frac{\tilde{\omega}_b}{\tilde{\omega}_c} \right) = (1 - 3\alpha) \frac{A^2}{4}. \quad (21) $$

In the absence of waves, the growth rate asymptotically approaches $\tilde{\omega} = \pm \beta^{1/2} \tilde{k} (\tilde{\omega}_b / \tilde{\omega}_c)^{1/2}$ (Meerson 1996). Equation (21) shows in general, assuming $\alpha$ and $A$ small and thus $\alpha^2 \approx \beta \tilde{\omega}_b / \tilde{\omega}_c + A^2/4$, that the Alfvén waves have a stabilizing effect.

When $\tilde{\omega}_c < 0$, the growth rate for $\tilde{k} \ll |\tilde{\omega}_c|$ remains finite, equal to $\tilde{\omega} = -i\tilde{\omega}_c$, independently of the presence of Alfvén waves.

For completeness, we now consider the case of adiabatic perturbations, corresponding to a wave mode with $\Re(\omega) \ne 0$. For this purpose, we restrict ourselves to the case of low amplitude waves, i.e. $A \ll 1$ and we set $\tilde{\omega} = \tilde{\omega}_r + i\tilde{\omega}_i$. We therefore have $\tilde{\omega}_r \gg \tilde{\omega}_i$. Moreover, it is assumed that $\tilde{\omega}_c \ll 1$ and $\tilde{\omega}_b \ll 1$ so that they can be neglected when multiplied by $\tilde{\omega}_i$. With these assumptions, one finds that to the first order

$$ \tilde{\omega} = \beta^{1/2} \tilde{k} + \frac{i}{2} (\tilde{\omega}_b - \tilde{\omega}_c) + \frac{A^2}{2\beta} \left( \frac{\beta^{1/2} \tilde{k} + i \tilde{\omega}_c}{\beta^{1/2} - 1} \right) \left( \frac{(\beta^{1/2} + \beta)\tilde{k}^2 + 1 - 3\beta^{1/2}}{(1 + \beta^{1/2})^2 \tilde{k}^2 - 4} \right). \quad (22) $$

The real part, $\omega_r \approx C_s k$ simply corresponds to a sonic wave.

When $A = 0$, one finds that the conditions for thermal stability of a sonic wave is $\tilde{w}_c - \tilde{w}_b \ge 0$, a criterion already obtained by Field (1965).

When $A \ne 0$, the stability criterion is modified according to Eq. (22), showing that the effect of the waves depends in a complex way on $\beta$ and $\tilde{k}$.

## 2.4. Growth rate and physical discussion

We now numerically solve Eq. (15) which can be rewritten as a sixth order polynomial containing 4 parameters to be specified, namely $\beta$ and $A$ (characterizing the Alfvén wave) and $\tilde{\omega}_b$ and $\tilde{\omega}_c$ (function of the thermal processes). The value of $\tilde{\omega}_c = \omega_c / \omega_0$ represents the ratio of the temporal period of the Alfvén wave divided by the cooling time and can be arbitrarily chosen. However, once $\tilde{\omega}_c$ is specified, $\tilde{\omega}_b$ depends on the thermal function. In order to estimate this parameter we use the standard cooling function of the neutral atomic ISM (Wolffre et al. 1995, 2003) which is used in Audit & Hennebelle (2005)

in the thermally unstable regime ($n = 3$ cm$^{-3}$ and $T \approx 500$ K). To obtain the dispersion relation we integrate Eq. (14) for $k/k_0$ between 0.01 and 100 using logarithmic spacing. The roots of the polynomials are obtained using the zroots subroutine (Press et al. 1992). Here, we restrict our attention to the unstable branch only. In order to verify our method, we have reproduced the dispersion relation for the decay instability of an isothermal gas presented in Goldstein (1978).

Figure 1 displays the results of the numerical integration of Eq. (14). First panel shows the adiabatic case for $\beta = 0.9$, $A = 0.5$ (dotted line), $\beta = 0.5$, $A = 0.5$ (short dashed line), $\beta = 0.5$, $A = 1$ (dot-dashed line), $\beta = 0.1$, $A = 0.5$ (double dot-dashed line), $\beta = 0.1$, $A = 1$ (long dashed line). As expected the circularly polarized Alfvén wave is unstable (decay instability) in a range of $k$ extending to a few times $k_0$. Both the growth rate and the largest unstable value of $k$ increase with $\beta^{-1}$ and $A$.

The second panel shows results for $\tilde{\omega}_c = \omega_c / \omega_0 = 1$. The full line corresponds to the hydrodynamical case whereas the others correspond to the same values of $\beta$ and $A$ as in the first panel. Various interesting features can be seen. i) When $k \to 0$ the effect of the Alfvén wave is to decrease the growth rate and therefore to stabilize these modes with respect to thermal instability. For $\beta = 0.5$, $A = 1$ and for $\beta = 0.1$, the modes whose wavenumber is smaller than $k \approx 1$ are perfectly stable. ii) The intermediate modes (i.e. $k \approx k_0$) are more unstable when $\beta$ is smaller and $A$ is higher. This is due to the decay instability that the wave undergoes for these values of $k$. iii) When $k \to \infty$, the growth rate is independent of $\beta$ and $A$. These features are in good agreement with the asymptotic limit $\omega \approx 0$ studied in the previous section. Cases with a larger $\beta$, namely 1 and 1.5, have also been explored but are not described here since they are not directly relevant for the regions of the ISM we consider in this paper. Although the decay instability disappears for $\beta \ge 1$ (for small enough amplitude), we find no qualitative difference with the cases $\beta \le 1$. In particular the intermediate wavelengths are still destabilized by the waves. This is in good agreement with the analytical study of the instability threshold (see Eq. (18)).

The third panel shows results for $\tilde{\omega}_c = 0.1$. In this case the cooling time is 10 times larger than the period of the Alfvén waves. The dispersion relation is qualitatively similar to the previous case. Quantitatively, however, the intermediate wavelengths are much more unstable than the small wavelengths. This is due to the fact that in this range of parameters, the cooling time being larger than the dynamical time of the waves, the fastest growing instability is the decay instability. Note that since $\omega_0$ is 10 times larger in this case than in the case displayed in panel 2, the value of $\omega/\omega_0$ in the hydrodynamical case (full line) is 10 times lower.

The fourth panel shows results for $\tilde{\omega}_c = 10$. In this case, the growth rate of the thermal instability is much lower than the growth rate of the decay instability. Therefore, the only effect of the Alfvén wave is to stabilize the gas with respect to thermal instability. This effect increases when $\beta$ decreases and when $A$ increases.

These results suggest that the presence of non-linear circularly polarized Alfvén waves in a thermally unstable medium like the neutral interstellar atomic gas can have two main effects. If the cooling time is short with respect to the temporal
<figure><img src="image_3.png" /><figcaption>Fig. 2. Growth rate as a function of the wave number when the effect of the thermal diffusivity is taken into account. The curve styles are associated with the same values of β and A as in Fig. 1.</figcaption></figure>

period of the waves, then the waves stabilize the gas. Therefore the gas can survive longer in the thermally unstable domain possibly leading in the ISM to a larger fraction of thermally unstable gas. Since the short wavelengths are the most unstable, the trend is that the CNM is fragmented into several clouds. If the cooling time is larger than the wave period, then the decay instability makes the intermediate modes ($k \approx k_0$) more unstable. In that case, the fraction of thermally unstable gas is not necessarily larger (depending on β and A) but the CNM should be fragmented in structures having a size of about $1/(k_0 \times \xi)$, where $\xi$ is the density ratio between the CNM and the WNM.

## 2.5. Effect of thermal diffusivity

Here we briefly consider the effect of thermal diffusivity. When this term is taken into account, the dispersion relation (15) becomes:

$$ \bar{\omega}^2 - \beta \bar{k}^2 \frac{\bar{\omega} + i(\bar{\omega}_b + \bar{\kappa}_b \bar{k}^2)}{\bar{\omega} + i(\bar{\omega}_c + \bar{\kappa}_c \bar{k}^2)} = \bar{k}^2 A^2 \frac{\bar{\omega}^3 + \bar{k}\bar{\omega}^2 - 3\bar{\omega} + \bar{k}}{(\bar{\omega} - \bar{k})((\bar{\omega} + \bar{k})^2 - 4)}, \quad (23) $$

$$ \text{where } \bar{\kappa}_c = \kappa(T_0)\omega_0\beta/(C_v\rho_0C_s^2) \text{ and } \bar{\kappa}_b = \bar{\kappa}_c/\gamma. $$

Figure 2 displays the growth rate obtained from Eq. (23) using the fiducial value of $\kappa(T_0)$ given at the beginning of Sect. 2 and for $\omega_c = 0.3, 1, 3$. The different curves are associated with the same values of β and A as in Fig. 1. For small and intermediate values of $\bar{k}$ the shape and the values of $\bar{w}$ are very similar to the case of vanishing thermal conductivity. As expected however, thermal conduction introduces a cut-off at small

scales (Field 1965). The value of $\bar{k}$ for which $\bar{\omega}$ vanishes depends on β and increases (by a factor of 2 to 3 for the values considered here) when β decreases, confirming the trends inferred previously, i.e. the CNM should be more fragmented in the presence of Alfvén waves leading to smaller CNM structures.

## 3. Numerical study

In order to test the analytic results presented in the previous section and to investigate the non-linear regime, numerical simulations are performed in a slab geometry.

For this purpose we use the 1D adaptive mesh refinement (AMR) code presented in Hennebelle & Pérault (1999, 2000). The AMR technique is very helpful to simultaneously resolve the sharp thermal fronts ($\approx 10^{-3}$ pc) and the larger scale ($\ge 10$ pc) involved in the problem. The code has been extensively tested and the results in the hydrodynamical case have been closely compared with high resolution simulations using a second order Godunov scheme. The growth rates for the parametric instability of an Alfvén wave in an adiabatic gas have been calculated and shown to match, within an accuracy of a few percent, the results of the first panel of Fig. 1.

Two numerical experiments are carried out. First, we consider a situation for which the gas is initially thermally unstable and we study the development of the thermal instability in the presence of Alfvén waves. Although these initial conditions are somewhat artificial, they are simple and close to the assumption of the analytic analysis, making comparison easier. Second, we set up more realistic initial conditions corresponding to a converging flow of thermally stable WNM. In this case the ram pressure of the flow drives the thermal collapse dynamically.

### 3.1. Case of initially thermally unstable gas

In order to study the effect of the circularly polarized Alfvén waves on the development of the thermal instability, we start the simulation with thermally unstable gas ($n = 3 \text{ cm}^{-3}$ and $T \approx 500$ K) with density fluctuations of amplitude 0.5. The cooling function is described in Audit & Hennebelle (2005). The computational domain, which initially contains 5000 pixels, has periodic boundary conditions and a length of $3 \times 10^{19}$ cm, corresponding to 20 wavelengths of the initial Alfvén wave.

Figure 3 shows three snapshots of the density field in a purely hydrodynamic run (no MHD wave is present). The initial perturbation grows as a result of the thermal instability. At time $t = 2.55$ Myr the density is about three times its initial value. At time $t = 3.64$ Myr a cloud of CNM having a size of about 0.3 pc has formed.

Figure 4 displays the density and the y-component of the magnetic field for 3 snapshots showing the development of the thermal instability in the presence of Alfvén waves of amplitude $B_\perp = 1 \mu\text{G}$ and for $B_x = 5 \mu\text{G}$. In this situation the value of $\omega_c/\omega_0$ is about 0.03, $A = 0.2$ and $\beta = 0.17$. The peak density at time $t = 1.82$ Myr is about 4.9 whereas the initial peak density is 4.5. In the hydrodynamic case the peak density at
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Density field for 3 snapshots illustrating the development of the thermal instability in the hydrodynamical case. The density of the perturbation increases until the gas reaches thermal equilibrium.</figcaption></figure>

the same time is about 6.2 showing that the waves have signifi-
cantly slowed down the growth of the perturbation, by a factor
of about (6.2–4.5)/(4.9–4.5) ≈ 4. This factor is significantly
larger than what is predicted by the linear theory which pre-
dicts a difference of about 10% (in the case of a perturbation
having an initial amplitude of 10⁻² we verified that the growth
agrees with the linear theory). After the central density has in-
creased by a factor of about 2 (panel 2), the waves drastically
change the structure of the gas and create significant density
contrasts (≈10%). The resulting cloud (panel 3) contains den-
sity fluctuations of about ≈50% of its maximum value, and
is therefore very different from the uniform cloud formed in
the hydrodynamical case. These large fluctuations are due to
magnetic pressure variations. Due to the contraction, the waves
inside the growing perturbation (panel 2 and 3) have a larger
amplitude and a shorter wavelength than the waves in the sur-
rounding medium. According to the analysis of the preceeding
sections, this effect tends to increase the influence of the waves.

Figure 5 shows results for $B_{\perp} = 2.5$ µG and for $B_x = 5$ µG. The waves strongly influence the gas evolution. The initial perturbation is stabilized by the waves and does not develop (second panel). On the contrary, the waves trigger the formation of structures having a wave number $k \approx k_0$ (first and second panels). These structures keep condensing (panel 3) and finally about 12 CNM clouds form. Their size is about

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Density and y-component of the magnetic field for 3 snapshots illustrating the development of the thermal instability in the presence of Alfvén waves of amplitude $B_{\perp} = 1$ µG and a longitudinal magnetic field $B_x = 5$ µG.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Same as Fig. 4 for $B_{\perp} = 2.5 \mu\text{G}$ and $B_x = 5 \mu\text{G}$.</figcaption></figure>

0.03 pc and therefore about 10 times smaller than the size of the cloud that forms in the hydrodynamical case. The formation of these clouds is significantly (50%–100%) faster than in

the two previous cases. Therefore, in this range of parameters, the waves destabilize the gas and accelerate the formation of CNM structures.

Figure 6 shows results for $B_{\perp} = 2.5 \mu\text{G}$ and for $B_x = 2.5 \mu\text{G}$. In this situation the value of $\omega_c/\omega_0$ is about 0.06, $A = 1$ and $\beta = 0.34$. As can be seen in the first panel, the initial perturbation does not grow (the peak density is smaller at time $t = 1.82$ Myr than the initial value). Unlike the previous case, the Alfvén waves do not efficiently trigger the formation of structures at $k \approx k_0$ (see first panel). However since the waves are unable to stabilize small wavelengths, the initial perturbation breaks down in several structures that finally develop (panel 2). This nevertheless occurs after a significant delay showing that in this situation, the gas spent a longer amount of time in the thermally unstable state. At time $t = 5.41$ Myr about 6 small structures of size $\approx 0.03$ pc and one larger structure of size $\approx 0.3$ pc have formed.

## 3.2. Case of dynamically induced thermal condensation

We now consider a converging flow of WNM (Hennebelle & Pérault 1999) in a simulation box of length 150 pc. The peak velocities are $2.2 C_{\text{unm}}$ and $-2 C_{\text{unm}}$ and the peak to peak distance is about 60 pc. Two simulations are performed. The first one is purely hydrodynamical, while the second one starts with an Alfvén wave of amplitude $2.5 \mu\text{G}$ with 100 spatial periods in the integration domain. The total magnetic intensity is $5 \mu\text{G}$ and the longitudinal one $B_x \approx 4 \mu\text{G}$. An initial resolution of 25 000 pixels is used to ensure an accurate description of the Alfvén wave (it corresponds to 250 pixels per period).

Figure 7 shows the hydrodynamical result. At time $t = 6.09$ Myr, two clouds with a size of about 0.1 pc have formed. They present weak density gradients and have no significant internal velocity. The two clouds have a relative velocity of about $0.15 \text{ km s}^{-1}$ and undergo a collision at time $t = 6.3$ Myr. A shock-compressed layer forms with a density of about $\approx 700 \text{ cm}^{-3}$ and a length $\approx 0.02$ pc. During the time of the collision, the structure presents a stiff velocity gradient.

Figure 8 displays four snapshots of the longitudinal velocity field, the density and the y-component of the magnetic field in the MHD case. The results displayed in the first panel confirm the trends observed in the numerical experiments of Sect. 3.1. The gas fragments into few small CNM structures having a physical length as small as a few 0.01 pc (note that due to the AMR scheme, these structures are well described). This situation is very different from the hydrodynamical case for which the structures are much larger and uniform. A comparison between the two times displayed in first panel of Fig. 8 reveals that all the structures do not form at the same time. The structure of the magnetic field varies very rapidly in the new condensations ($x = 242$ cm and $x = 243$ cm), therefore compressing them. It varies less steeply near older structures ($x \approx 244$ cm) since the field lines have time to unbend. Between the clouds, the magnetic field is much more uniform. This is a consequence of the fact that the Alfvén speed is about 10 times larger in the WNM than in the CNM. It is also clear that the
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same as Fig. 4 for $B_{\perp} = 2.5 \mu\text{G}$ and $B_{\chi} = 2.5 \mu\text{G}$.</figcaption></figure>

intercloud magnetic pressure plays an important role in preventing the merging of the CNM structures therefore maintaining the complexity of the flow and organizing it into groups

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Spatial zoom showing the thermal condensations induced by the large scale converging flow in the hydrodynamical case. Two snapshots are displayed.</figcaption></figure>

of structures rather than into a single cloud. In these groups of structures the longitudinal velocity dispersion is not negligible, unlike in the hydrodynamical case. It is about 0.1–0.2 km s⁻¹, i.e. a few times the sound speed or the Alfvén speed of the CNM. This velocity dispersion is due to the transfer of magnetic energy into longitudinal motions because of the magnetic pressure fluctuations. Observing such group of structures with a low spatial resolution (≈0.2 pc) may lead to a rather different picture, namely a broad, uniform and turbulent (having a Mach number $M \approx 1–2$) CNM structure.

The second panel of Fig. 8 shows two later snapshots of the same numerical experiment. The group of structures seen in the first panel is now located at $x \approx 251.5$ pc. Large fluctuations (≈100%) of density and magnetic field are still present as well as a longitudinal velocity dispersion of about 0.1–0.2 km s⁻¹. Another smaller group of structures ($x \approx 253.5$ pc) has formed. At time $t = 9.11$ Myr (dotted line), the first cloud undergoes a large density fluctuation ($n_{\text{max}} \approx 600$ cm⁻³) on a scale of about 0.01 pc. At the same time the cloud is compressed, because of the (magnetic) interaction with the other cloud, so that its length is divided by a factor $\approx 2$. An interpretation based on the previous analytic results is that since the amplitude and wavenumber of the Alfvén wave in the cloud increases as a result of the compression, the decay instability is triggered leading to a larger density fluctuations. Interestingly, in contrast to the large fluctuations undergone by $B_y$, the longitudinal velocity field remains relatively smooth.

## 4. Discussion and conclusion

The analysis presented in Sect. 2 as well as the numerical experiments discussed in the previous section show that even modest amplitude Alfvén waves may have a strong impact on the structure of the multiphase ISM. This study focused on circularly polarized parallel propagating Alfvén waves, mainly to allow analytic calculations. These waves, which are exact
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as Fig. 7 in the presence of circularly polarized Alfvén waves (see text for detail). Four snapshots are displayed.</figcaption></figure>

solutions of the MHD equations, are very weakly dissipative and are therefore very likely to be present in the ISM.

Their effects depend on $\beta$, $\tilde{\omega}_c$ and $A$, respectively the square ratio of the sound to the Alfvén speeds, the ratio between the wave temporal period and the cooling time and the wave amplitude. Depending on the values of these parameters, these waves may: i) stabilize the wavelengths larger than that of

the Alfvén wave, $\lambda_{AW}$, that would otherwise be thermally unstable, therefore enhancing the fraction of thermally unstable gas in the ISM; ii) destabilize the wavelengths comparable to $\lambda_{AW}$ and thus fragment the CNM into several spatially correlated small clouds; iii) induce strong density fluctuations within preexisting CNM structures (up to 10 times the mean density); iv) maintain an Alfvénic velocity dispersion within the CNM structures by pumping their energy into longitudinal motions. Finally magnetic pressure tends to prevent the merging between CNM clouds.

These effects are not observed in one-dimensional purely hydrodynamical simulations leading to larger and almost uniform structures with weak internal motions. In two dimensions (Koyama & Inutsuka 2002; Audit & Hennebelle 2005) the situation is more complex due to the role of turbulence, with the coexistence of thermally unstable gas and small-scale structures, but the latter are still relatively uniform and do not contain a significant velocity dispersion. According to the present analysis, MHD waves may enhance the thermal fragmentation found in these simulations, maintain an Alfvénic velocity dispersion and generate large density fluctuations within the CNM structures.

Recent observational progress has revealed the presence of interesting features, and of a large quantity of thermally unstable gas (e.g. Heiles 2001; Heiles & Troland 2003; Miville-Deschênes et al. 2003). Heiles (1997) summarizes and discusses observations of tiny small-scale structures and more recently Braun & Kanekar (2004, 2005) and Staminorovic & Heiles (2005) report the detection of very low column density CNM clouds.

Qualitatively, the present study as well as the 2D simulations of Koyama & Inutsuka (2002) and Audit & Hennebelle (2005) show similar features including thermally unstable gas (also observed by Gazol et al. 2001 in simulations at larger scale), low column density structures (down to and even smaller than $10^{18}$ cm$^{-2}$) and large density fluctuations at small scale (up to $10^3$ cm$^{-3}$). Small-scale structures and large density fluctuations appear to be a natural outcome of the two-phase nature of the flow. In the absence of turbulence or a magnetic field, CNM structures have a typical length of about 0.1 pc (or a column density of $\approx 10^{19}$ cm$^{-2}$) (Audit & Hennebelle 2005). In the presence of waves or turbulence, a growing structure can however fragment into smaller pieces. Moreover, since both the sound and the Alfvén speed change by a factor $\approx 10$ between the WNM and the CNM, with a transition occurring over a short distance of the order of a Field length, supersonic or super-Alfvénic CNM motions are expected to be more frequent.

No systematic and quantitative study attempting to closely compare simulations and observations has been carried out yet. This is clearly a major challenge for the future.

***Acknowledgements.*** This work has received partial financial support from the French national program PCMI.

## References

Audit, E., & Hennebelle, P. 2005, A&A, 433, Braun, R., & Kanekar, N. 2004, in the IMF at 50, ed. E. Corbelli, F. Palla, & H. Zinnecker (Kluwer) [arXiv:astro-ph/0409427]
Braun, R., & Kanekar, N. 2005, A&A, 436, L53

Cho, J., & Lazarian, A. 2003, Mon. Not. R. Astron. Soc., 345, 325

Cohen, R. H., & Kulsrud, R. M. 1974, Phys. Fluid, 17, 2215

Dewar, R. L. 1970, Phys. Fluid, 13, 2710

Falle, S. A., & Hartquist, T. W. 2002, MNRAS, 329, 195

Ferraro, V. C. A. 1955, Proc. R. Soc. London, Ser. A, 223, 310

Field, G. 1965, ApJ, 142, 531

Field, G., Goldsmith, D., & Habing, H. 1969, ApJ Lett., 155, 149

Gazol, A., Vázquez-Semadeni, E., Sánchez-Salcedo, F., & Scalo, J. 2001, ApJ, 557, L124

Goldstein, M. 1978, ApJ, 219, 700

Heiles, C. 1987, Interstellar processes, ed. D. Hollenbach, & H. Thronson (Reidel)

Heiles, C. 1997, ApJ, 481, 193

Heiles, C. 2001, ApJ, 551, 105

Heiles, C., & Troland, T. 2003, ApJ, 586, 1067

Heiles, C., & Troland, T. 2005, ApJ, in press
[arXiv:astro-ph/0501482]

Hennebelle, P., & Pérault, M. 1999, A&A, 351, 309

Hennebelle, P., & Pérault, M. 2000, A&A, 359, 1124

Koyoma, H., & Inutsuka, S. 2002, ApJ, 564, L97

Kulkarni, S. R., & Heiles, C. 1987, Interstellar processes, ed. D. Hollenbach, & H. Thronson (Reidel)

Loewenstein, M. 1990, ApJ, 349, 471

Lou, Y. Q. 1996, MNRAS, 279, L67

McKee, C. F., & Zweibel, E. G. 1995, ApJ, 440, 686

Meerson, B. 1996, Rev. Mod. Phys., 68, 215

Miville-Deschênes, M.-A., Joncas, G., Falgarone, E., & Boulanger, F. 2003, A&A, 411, 109

Passot, T., Vázquez-Semadeni, E., & Pouquet, A. 1995, ApJ, 441, 702

Passot, T., & Vázquez-Semadeni, E. 2003, A&A, 398, 845

Piontek, R. A., & Ostriker, E. C. 2004, ApJ, 601, 905

Piontek, R. A., & Ostriker, E. C. 2005, ApJ, in press
[arXiv:astro-ph/0504669]

Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 1994, Numerical Recipies (Cambridge University Press)

Rogister, A. 1971, Phys. Fluids, 12, 2733

Shu, F. 1992, Gas Dynamics (Mill Valley CA: University Science Books)

Spangler, S., Fuselier, S., Fey, A., & Anderson, G. 1988, J. Geophys. Res., 93(A2), 845-857, 10.1029/88JA01034

Stanimirović, S., & Heiles, C. 2005, ApJ, in press

Troland, T., & Heiles, C. 1986, ApJ, 301, 339

Wolfire, M. G., Hollenbach, D., & McKee, C. F. 1995, ApJ, 443, 152

Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2003, ApJ, 587, 278