Interstellar Cloud Formation

# Formation of Interstellar Clouds: Parker Instability with Phase Transitions

Telemachos Ch. Mouschovias, Matthew W. Kunz & Duncan A. Christie
Departments of Physics and Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801
Accepted 2009 January 07. Received 2009 January 06; in original form 2008 October 08
###### Abstract

We follow numerically the nonlinear evolution of the Parker instability in the presence of phase transitions from a warm to a cold Hi interstellar medium in two spatial dimensions. The nonlinear evolution of the system favors modes that allow the magnetic field lines to cross the galactic plane. Cold Hi clouds form with typical masses M, mean densities cm, mean magnetic field strengths G (rms field strengths G), mass-to-flux ratios relative to critical, temperatures K, (two-dimensional) turbulent velocity dispersions km s, and separations pc, in agreement with observations. The maximum density and magnetic field strength are cm and G, respectively. Approximately 60% of all Hi mass is in the warm neutral medium. The cold neutral medium is arranged into sheet-like structures both perpendicular and parallel to the galactic plane, but it is also found almost everywhere in the galactic plane, with the density being highest in valleys of the magnetic field lines. ‘Cloudlets’ also form whose physical properties are in quantitative agreement with those observed for such objects by Heiles (1967). The nonlinear phase of the evolution takes Myr, so that, if the instability is triggered by a nonlinear perturbation such as a spiral density shock wave, interstellar clouds can form within a time suggested by observations.

###### keywords:
instabilities — ISM: clouds — ISM: magnetic fields — ISM: structure — MHD
pagerange: Formation of Interstellar Clouds: Parker Instability with Phase TransitionsReferencespubyear: 2009

## 1 Introduction

Understanding the formation of interstellar clouds and cloud complexes is a problem of longstanding importance, one which has received increased attention in recent years. Since Hi clouds are the progenitors of molecular clouds, the stellar birthplaces, the process by which they form is of fundamental importance to galactic evolution as well as to the early stages of star formation.

Proposed mechanisms for interstellar cloud formation fall into three categories: stochastic coagulation of smaller clouds (Oort, 1954; Field & Saslaw, 1965; Kwan, 1979; Scoville & Hersh, 1979), colliding streams of turbulent flows (Ballesteros-Paredes, Hartmann & Vázquez-Semadeni, 1999; Hartmann, Ballesteros-Paredes & Bergin, 2001; Heitsch et al., 2006), or collective effects involving instabilities (Mouschovias, 1974; Mouschovias, Shu & Woodward, 1974; Blitz & Shu, 1980). Stochastic coagulation is not feasible because collisional agglomeration proceeds too slowly and, also, collisions are more likely to result in disruption rather than merger. Turbulence-driven formation of atomic-hydrogen clouds is problematic because no realistic mechanism has yet been demonstrated to sustain the required large-scale, ordered, focused motions for Myr. In addition, the short lifetimes and highly supercritical mass-to-flux ratios of molecular clouds thought to form within these turbulent ‘atomic precursors’ are in conflict with observations (Mouschovias, Tassis & Kunz, 2006). 111The typical separation between the first appearance of CO and the peak of H in spiral arms corresponds to a timescale yr. This effect cannot possibly be due to dust obscuration, as proposed in García-Burillo, Guélin & Cernicharo (1993), Kim & Ostriker (2006), and Elmegreen (2007). A similar shift occurs between the first appearance of CO and the peak of Pa, which does not suffer from dust obscuration (Scoville et al., 2001). By contrast, large-scale dynamical instabilities can amass material over lengthscales much greater than the sizes of diffuse atomic-hydrogen clouds in a time consistent with observational constraints.

In the Parker (1966) instability, undulations of initially horizontal (i.e., parallel to the galactic plane) magnetic field lines allow matter to slide down along field lines under the action of the vertical component of the galactic gravitational field. The raised portions of the field lines are thereby relieved of the confining weight of the gas and rise even higher, until the tension of the field lines increases sufficiently for a final equilibrium state to be established (Mouschovias, 1974). Mouschovias et al. (1974) have argued that the Parker instability is triggered in spiral arms behind spiral density shock waves, leading to the formation of cloud complexes in the valleys of the curved magnetic field lines.

Several consequences of this mechanism for the formation of interstellar clouds are in excellent, in fact unique, agreement with observations:

1. The implosion of individual clouds by shock waves within these complexes leads to the formation of OB associations and giant Hii regions aligned along spiral arms ‘like beads on a string’ and separated by regular intervals of pc, as observed both in our galaxy and in external galaxies (Westerhout, 1963; Kerr, 1963; Hodge, 1969; Morgan, 1970).

2. The observed moderate enhancement in synchrotron radiation from the interarm to the arm region of spiral galaxies (Mathewson, van der Kruit & Brouw, 1972; van der Kruit, 1973a, b, c) found a natural explanation in the context of the Parker instability (Mouschovias et al. 1974; Mouschovias 1975a). The relief of magnetic stresses by buckling of the field lines in the vertical direction in spiral arms leads to an increase in the magnetic field strength that is much smaller than that expected from one-dimensional compression. Moreover, the tendency for the cosmic-ray pressure to distribute itself uniformly along magnetic flux tubes means that, after the instability develops, most of the relativistic particles are found in the regions where the field is weak, i.e., in the inflated portions of the flux tubes. In other words, the cosmic-ray and magnetic-energy densities are not in equipartition on scales kpc; there is an anticorrelation between them.

3. The anticorrelation between the cosmic-ray density on the one hand and the magnetic-field strength and gas density on the other hand also implies that the intensity of synchrotron radiation along spiral arms should not vary significantly, as observed. This is in sharp contrast with expectations based on correlation between the field strength and the cosmic-ray density (e.g., equipartition between magnetic and cosmic-ray energy densities).

Numerical simulations of the Parker instability with applications to the interstellar medium have been carried out by Basu, Mouschovias & Paleologou (1997; hereafter BMP97), Kim et al. (1998), Kim, Ryu & Jones (2001), and Kim, Ostriker & Stone (2002). Several questions concerning the nonlinear evolution of the instability were addressed. For example, which mode dominates the nonlinear growth of the instability on timescales of interest in the interstellar medium (ISM); how fast the instability grows in the nonlinear regime; and, how great a density enhancement is generated in the galactic plane by downflow along the magnetic arches. A common finding in these papers is that the Parker instability by itself cannot be the agent responsible for interstellar cloud formation; the maximum density enhancement during the nonlinear evolution is only a factor of . Instead, the Parker instability may be thought of as setting the stage upon which smaller-scale processes, such as cloud formation, star formation, and subsequent supernova explosions, act out their individual roles and collectively contribute to the structure and evolution of the ISM. However, Mouschovias (1974) and BMP97 also noted that the Parker instability may be directly responsible for the formation of GMCs, or cloud complexes, if it can generate a sufficient density enhancement near the galactic plane for a critical pressure (Field, Goldsmith & Habing, 1969; Kaplan & Pikel’ner, 1970) to be exceeded, thereby initiating a phase transition to a cooler, denser cloud phase. This point has been re-emphasized by Parker & Jokipii (2000).

We follow numerically the nonlinear evolution of the Parker instability in two spatial dimensions accounting for phase transitions. In Section 2, we present the basic equations and physics of the instability. In Section 3, the dimensionless equations are given and the numerical method of solution is described. The results of the nonlinear development of the instability with phase transitions from a warm neutral medium (WNM) to a cold neutral medium (CNM) are presented in Section 4, followed by a summary and discussion in Section 5. Contact between theory and observations is emphasized.

## 2 Formulation of the Problem

### 2.1 Basic Equations

We consider a conducting fluid of mass density , thermal pressure , and temperature , threaded by a frozen-in magnetic field , in an external gravitational field . The magnetohydrodynamic (MHD) equations describing this system are

 ∂ρ∂t+\boldmath{∇}\boldmath{⋅}(ρ\boldmath{v})=0, (1a)
 ∂(ρ\boldmath{v})∂t+\boldmath{∇}\boldmath{⋅}(ρ\boldmath{v}\boldmath{v}% )=−\boldmath{∇}P+ρ\boldmath{g}+14π(%\boldmath$∇$\boldmath{×}\boldmath{B})\boldmath{×}\boldmath{B}, (1b)
 ∂\boldmath{B}∂t=\boldmath{∇}\boldmath{×}(\boldmath{v}\boldmath{×}\boldmath{B}), (1c)
 ∂u∂t+\boldmath{∇}\boldmath{⋅}(u\boldmath{v})=−P\boldmath{∇}% \boldmath{⋅}\boldmath{v}−ρL+\boldmath{∇}\boldmath{⋅}(κ\boldmath{∇}T), (1d)

where is the velocity and is the internal energy density. The thermal pressure is a function of both density and temperature, , and is determined by a balance of heating and cooling mechanisms in the ISM. It is obtained from the internal energy density of the gas through the ideal gas relation, , where is the ratio of specific heats (). In equation (1d), is the loss function (i.e., net cooling rate per unit mass), where and are the cooling and heating rates, respectively; is the thermal conductivity (see eq. 12 below).

We use an equation of state that allows the existence of two phases, a warm one with K and density cm and a cold one with K and cm. Following Nagashima, Inutsuka & Koyama (2006), we adopt the following heating and cooling rates:222The functional form of these heating and cooling rates has its origin in earlier work by Penston & Brown (1970) and references therein.

 μmHΓ = 2×10−26ergs−1, (2a) μmHΛ(T)Γ = 107exp(−1.184×105T+1000) (2b) + 1.4×10−2√Texp(−92T)cm3,

where is the mean mass per particle in units of the atomic-hydrogen mass, accounting for 10% He abundance by number. The minimum and maximum pressures for the coexistence of the warm/cold phases in a static equilibrium are K cm and K cm. The corresponding transition temperatures K and K define cold (), warm (), and intermediate temperature () phases. Since we do not include the additional phase transition to molecular hydrogen in our loss function, we enforce isothermality for any gas element whose number density exceeds cm. The corresponding temperature of this isothermal branch is K. The thermal equilibrium curve implied by these heating and cooling functions is shown in Figure 1. The pressure (ordinate) and density (abscissa) are normalized to their initial midplane values (see § 3.1 below). The dashed lines denote the -, -, and -isotherms.

Note that an investigation of the Parker instability accounting for phase transitions is not the same as an investigation of the thermal instability occurring alongside the Parker instability (Ames, 1973). The lengthscales involved in the two instabilities are so radically different ( pc for the thermal instability and kpc for the Parker instability) that the two instabilities proceed in parallel without much effect of one on the other. The key element in the possibility of the Parker instability leading directly to the formation of massive Hi clouds and Hi cloud complexes is not so much the detailed nature of the thermal instability but the fact that a transition from the warm to the cold phase can take place by some means (e.g., even if the equation of state on a - diagram had only a horizontal branch joining the two phases, instead of the branch with a negative slope).

## 3 Method of Solution

### 3.1 Initial and Boundary Conditions

We consider a local Cartesian coordinate system corresponding to cylindrical polar coordinates in a galactic disk, with the range of the (radial) coordinate being (the distance from the galactic center) so that curvature terms can be ignored. The magnetic field is initially in the azimuthal direction (), which we take to be along a spiral arm. Because the motions involved in the Parker instability are predominantly azimuthal, rotation can be ignored (see Shu 1974; also Mouschovias et al. 1974).

Observations show that the interstellar magnetic field is predominantly parallel to the galactic plane, although arches of magnetic field lines rising high above the plane are also revealed (Mathewson & Ford, 1970). Parker (1966) assumed a zeroth-order initial (reference) equilibrium state, which, for simplicity, has field lines exactly parallel to the galactic plane (), so that . The vertical galactic gravitational field (due to stars) is taken to be constant and to reverse its direction across the galactic plane; i.e., , where constant . More complicated gravitational fields have been considered analytically by Kellman (1972), Giz & Shu (1993), Kim, Hong & Ryu (1997), Kim & Hong (1998), and Kim et al. (2000) and numerically by Santillán et al. (2000) and Franco et al. (2002). The main effect of a non-constant gravitational acceleration is to increase the growth rate of the instability and shorten the separation of the magnetic valleys.

Although we regard the Parker instability as a driven (or forced) instability in galactic spiral arms and, therefore, the ‘initial’ state is not an equilibrium one, we nevertheless use Parker’s initial equilibrium state as a reference state, from which we obtain the input parameters of the problem. This approach ensures that (i) any ‘final’ outcome of the evolution of the gas-field-gravity system can be reached from Parker’s stratified equilibrium state through continuous deformation of the field lines; and, (ii) the results of the present calculation can be compared and contrasted with those of past calculations on the nonlinear evolution of the Parker instability in the ISM. We thus take the ratio of magnetic and thermal pressures,

 α≡B2ref/8πPref, (3)

to be constant in the reference equilibrium state, in which the density is given by

 ρref(z)=ρref(0)exp(−|z|H), (4)

where

 H≡(1+α)C2ref(0)/g (5)

is the total scale height of the gas in the reference state and is the isothermal sound speed in the reference state at the galactic plane. With K, cm s, and , a value consistent with observations, one finds that pc. The magnetic field strength in this reference state is given by

 Bref(z)=Bref(0)exp(−|z|2H). (6)

The corresponding pressure and temperature distributions are determined by requiring the internal energy to satisfy local thermal equilibrium, based on the heating and cooling functions given in Section 2.1. For simplicity, we ignore the effect of cosmic rays in this study, although their effects on the evolution of the Parker instability are well known by both analytic arguments (Mouschovias, 1975b, 1996; Hanasz & Lesch, 2000) and numerical work (Hanasz & Lesch, 2003; Ryu et al., 2003; Kuwabara, Nakamura & Ko, 2004). Mainly, they tend to speed up the instability.

The solution of the linear instability problem can either have field lines with reflection symmetry about (midplane-symmetric or ‘odd’ mode) or allow field lines to cross (midplane-crossing or ‘even’ mode). In the latter case, the quantities , , and are invariant with respect to a reflection about followed by a translation in the -direction by an amount , where is the horizontal (along the field lines) wavelength of the perturbation. (The terms ‘even’ and ‘odd’ refer to the algebraic form of the velocity perturbation imposed on the reference equilibrium state.) The ‘odd’ mode in an isothermal system does not lead to any appreciable density enhancement, and is therefore less likely than the ‘even’ mode to produce interstellar clouds in a nonisothermal system. Moreover, BMP97 found that the ‘even’ mode is naturally selected by the nonlinear evolution of the system. Hence, we study a model with an initial perturbation that has, and at later times leads to evolution that preserves, the symmetry of the even mode.333The simulations by Kim et al. (1998, 2001) could only study the ‘odd’ mode, due to their assumed symmetry about the galactic plane. BMP97, Santillán et al. (2000), and Kim et al. (2002) considered both the ‘even’ and ‘odd’ modes. This symmetry allows us to study a region of horizontal extent .

The equations are put in dimensionless form by choosing units natural to the problem. The units of velocity , density , and magnetic field strength [] are, respectively, the values of the isothermal sound speed , density , and magnetic field at the galactic plane.444One could choose the unit of to be , but it is more convenient to use ; the two are identical for . The unit of length [] is the thermal-pressure scale height . The implied unit of time is , the sound-crossing time across a thermal-pressure scale height. While the evolution of the Parker instability depends on only one free parameter, , specific values of , , and must be chosen to put the loss function in dimensionless form. We choose g cm (corresponding to a number density of 1 cm), cm s, and K. (This corresponds to a thermally stable state in the warm phase.) Choosing the free parameter , we find the following values for the units: [] = 5.73 km s, [] = 35.45 pc, [] = 6.05 Myr, and [] = 4.17 G.

We choose , so that (the horizontal wavelength corresponding to the maximum growth rate of the linear instability). The density and velocity are assumed to have reflection symmetry about the planes and , where the field lines remain locally horizontal. The upper and lower boundaries, where the field lines remain undeformed, act as ‘lids’ to the system, and are placed at , so that they are far enough from the galactic plane not to affect the evolution of most of the matter between them.555Several numerical studies of the Parker instability miss this crucial point. For example, the simulations by Kim et al. (2002) and Kim & Ostriker (2006) have a total vertical box size of only (with ). Our total vertical box size, by contrast, is kpc. A small vertical size of the computational box suppresses the instability [see eq. (32) of Mouschovias (1974); also, eq. (22) of the review by Mouschovias (1996)]. (In the reference state the density at decreases to of its value at the galactic plane.) The velocity perturbation is taken to have the form

 δvz(y,z)=−εCref(0)cos(πyY)cos(πz2Z), (7)

where , so that the kinetic energy of the perturbation is less than 1% of any other form of energy (gravitational, thermal, or magnetic) in the system. This choice of the amplitude of the perturbation does not, of course, represent the highly nonlinear disturbance introduced in nature by a spiral density shock wave, and it also has the consequence that the instability takes a longer time to enter the nonlinear regime than it does in nature. This choice, nevertheless, allows us to follow the development of the linear phase of the instability as well as avoid the introduction of an arbitrariness in the problem by using such large perturbations that dominate from the outset the system’s evolution. To mimic the natural phenomenon of the instability driven by a spiral density shock wave without having to resort to a three-dimensional calculation, we set the origin of time () at the instant the linear phase of the evolution has led to an increase of the maximum density in the system by a factor of 2. Times earlier than this instant are assigned negative values in what follows.

### 3.2 The Dimensionless Equations

We denote the dimensionless time by and the dimensionless coordinates by , so that

 ∂∂t→gCref(0)∂∂τand\boldmath{∇}→gC2ref(0)\boldmath{∇}′.

All other dimensionless quantities are adorned with a tilde. The dimensionless equations are then

 ∂˜ρ∂τ+\boldmath{∇}′\boldmath{⋅}(˜ρ˜% \boldmath{v})=0, (8a)
 ∂(˜ρ˜\boldmath{v})∂τ+\boldmath{∇}′\boldmath{⋅}(˜ρ˜\boldmath{v}˜\boldmath{v})=−\boldmath{∇}′˜P−˜ρ^\boldmath{e}z−α\boldmath{∇}′˜B2+2α˜\boldmath{B}\boldmath{⋅}\boldmath{∇}′˜\boldmath{B}, (8b)
 ∂˜\boldmath{B}∂τ=\boldmath% {∇}′\boldmath{×}(˜\boldmath{% v}\boldmath{×}˜\boldmath{B}), (8c)
 ∂˜u∂τ+\boldmath{∇}′\boldmath{⋅}(˜u˜\boldmath{v})=−˜P\boldmath{∇}′\boldmath{⋅}˜\boldmath{v}−˜ρ˜L+\boldmath{∇}′\boldmath{⋅}(˜κ\boldmath{∇}′˜T), (8d)

where is the dimensionless pressure, and

 ˜u≡uρref(0)C2ref(0),˜L≡LgCref(0),˜κ≡κgTref(0)ρref(0)C5ref(0), (9)

are the dimensionless internal energy, loss function, and conductivity, respectively. The dimensionless density, magnetic field, and velocity perturbations in the reference state are, respectively,

 ˜ρref(z′)=exp(−|z′|1+α), (10a)
 ˜Bref(z′)=exp(−|z′|2(1+α)), (10b)

and

 ˜vz(y′,z′)=−εcos(πy′Y′)cos(πz′2Z′). (10c)

These equations are subject to the boundary conditions

 ˜vz(y′,±Z′)=0, ˜vy(0,z′)=0, ˜vy(Y′,z′)=0, (11a) ˜Bz(y′,±Z′)=0, ˜Bz(0,z′)=0, ˜Bz(Y′,z′)=0, (11b)
 ∂˜ρ(y′,z′)∂y′∣∣∣y′=0,Y′=0, ∂˜ρ(y′,z′)∂z′∣∣∣z′=±Z′=0, (11c) ∂˜vz(y′,z′)∂y′∣∣∣y′=0,Y′=0, ∂˜vy(y′,z′)∂z′∣∣∣z′=±Z′=0, (11d) ∂˜By(y′,z′)∂y′∣∣ ∣∣y′=0,Y′=0, ∂˜By(y′,z′)∂z′∣∣ ∣∣z′=±Z′=0. (11e)

### 3.3 Numerical Issues

We integrate equations (8a) - (8d) using a version of the publicly-available Zeus-MP code (Hayes et al., 2006) distributed by the Laboratory for Computational Astrophysics of the University of California at San Diego. Zeus-MP uses a time-explicit, operator-split, finite-difference method for solving the MHD equations on a staggered mesh, capturing shocks via quadratic and linear artificial viscosities. Details concerning the algorithms used can be found in Hayes et al. (2006) and references therein.

We have made modifications to this code in order to follow phase transitions. Because cooling times can be very short compared to the timescales relevant to cloud formation, the energy equation update from the net cooling terms is solved implicitly using Newton-Raphson iteration. At the start of each iteration, the timestep is initially computed from the Courant-Friedrichs-Levy condition using the sound speed, Alfvén speed, and conduction parameter (see below). This is followed by a call to the cooling subroutine. The change in temperature within each zone is limited to 10% of its initial value. If this requirement is not met for all cells in the computational domain, the timestep is reduced by a factor of 2, and the implicit energy update is recalculated. As in Piontek & Ostriker (2004), we have found that the timestep is reduced at most once or twice per convergent timestep.

The update from the conduction operator is solved explicitly, using a second-order accurate differencing of . As Koyama & Inutsuka (2004) pointed out, the importance of incorporating conduction in simulations that contain thermally unstable gas has sometimes been overlooked in past work. Without conduction, the growth rates for thermal instability are greatest at the smallest lengthscales, and unresolved growth at the grid scale may occur. The inclusion of conduction has a stabilizing effect on the thermal instability at small scales, and the conduction parameter can be adjusted to allow spatial resolution of the thermal instability on the computational grid. We tune the thermal conduction coefficient so that the number of zones in a Field length satisfies in the thermally unstable regime:666The Field length is when is a function of temperature. It is essentially the critical lengthscale inside which conduction stabilizes the thermal instability.

 κ=ρ2ΛT(2Δyπ)2max[1−(∂lnΛ∂lnT),0]. (12)

As a parcel of gas transits (thermodynamically) from an unstable state to a stable state, the thermal conductivity vanishes in a continuous fashion.

The timestep is necessarily limited not only by the usual sound and Alfvén cell-crossing times, but also by a diffusion timescale that is set by the thermal conduction parameter:

 (Δt)cond=ρ4(γ−1)(Δy)2κ. (13)

The quadratic dependence on the grid spacing is a general feature of any numerical problem that includes diffusion. This dependence is particularly restrictive when high resolution is required to accurately capture important physical processes. In this work, a high resolution is needed so that the thermal conduction parameter is restricted to reasonable values.

The problem is solved on an orthogonal, nonuniform grid. The horizontal () grid is spatially uniform, with 2048 zones and a resolution pc (i.e., ). The vertical () grid contains 4096 zones and is in general nonuniform, with the resolution being greatest near the galactic plane. There are 2048 uniformly-spaced zones in the 15 initial thermal-pressure scale heights straddling the galactic plane, so that pc (i.e., ); the remaining 2048 zones are logarithmically-spaced above and below this uniform-grid region. The transition between the uniform and nonuniform regions is chosen so that the grid spacing varies smoothly.

## 4 Results

In Figure 2, we show the time evolution of the central (dimensionless) (a) density, (b) temperature, and (c) magnetic field strength at , which is located at the bottom of a magnetic valley. As found by previous studies of the isothermal Parker instability, there are three distinct phases of evolution: linear, nonlinear, and relaxation into a ‘final’ equilibrium state.

During the linear phase (), the velocities remain subsonic and the central density, temperature, and magnetic-field strength vary only slightly. Soon after the onset of the nonlinear phase (at ), the central density, temperature and magnetic-field strength change very rapidly. After the density has increased by a factor of , the pressure of many fluid elements has also increased above the critical value for a phase transition to occur to a cold (cloud) phase. This leads to a rapid density enhancement by more than two orders of magnitude in a time (i.e., Myr). During the phase transition, approximately isobaric evolution causes the temperature to decrease by about two orders of magnitude. The magnetic field strength increases only slightly during the early part of the nonlinear phase, since gas motions are primarily along field lines. However, after enough gas has accumulated in the valleys of the field lines, its weight compresses the field lines and the maximum magnetic-field strength increases by a factor of three (from G to G). During the subsequent relaxation phase, the compressed gas rebounds because it had overshot its equilibrium state, and eventually settles in a state that has a central density cm. The magnetic field strength continues to increase as more mass elements enter the cloud phase and the galactic gravitational field (i) compresses the cloud perpendicular to the galactic plane and (ii) flattens it in the galactic plane. The resulting physical conditions are such that, if we had included the gas’ self-gravity in the calculation, individual self-gravitating clouds would separate out. The density oscillations evident after occur about a ‘final’ equilibrium state; they are the result of motions left over by the cloud formation process.

As explained in Section 3.1, the linear phase of evolution is expected to be very short-lived in nature because the instability is externally driven by a spiral density shock wave (see Mouschovias et al. 1974). For this reason, the ‘origin’ of time is chosen to be the instant by which the central density has increased by a factor of 2. It is thus seen that the time required to form cold interstellar cloud complexes from the WNM does not exceed Myr.

Figure 3 exhibits contour plots of the density (color) and magnetic field lines (solid lines) at times , , , 2, 7, and 12. For ease in visualization, we show a region that covers one full horizontal wavelength, , of the instability. Regions of high (low) density are colored red (blue) with a color continuum between them (as shown in the colorbar to the right of the figures). By the end of the linear phase of the evolution (), the magnetic field lines have begun to deform significantly and the gas has begun its rapid descent into the magnetic valleys. The maximum velocity has exceeded the sound speed by a factor of 1.61. The evolution up to this point is nearly indistinguishable from the evolution of the isothermal Parker instability (determined by BMP97), because few gas elements have transitioned into the stable CNM phase by this time (see below). By , the system is well into the nonlinear stage of its evolution. The field lines have become considerably arched, even near . There is a strong flow of gas [maximum velocity ] into the magnetic valleys, leading to the formation of shocks. By , the rising portions of the field lines have almost stopped their expansion and are nearly at their ‘final’ equilibrium positions.

In Figure 4, we display the large-scale velocity (arrows) and density (solid lines) fields at time . There is considerable density structure in the galactic plane. The velocity field shows prominent and distinctive S-shaped shocks extending far above and below the galactic plane. Shocks of this character have also been seen in the BMP97 simulations of the isothermal Parker instability. Mouschovias (1996) noted that, if such shocks are observed, they cannot possibly be mistaken for supernova remnants and can therefore be taken as evidence for the nonlinear development of the Parker instability in the interstellar medium.

Figure 5a shows the innermost , part of the system shown in Figure 3 at time . The CNM has a clear sheet-like morphology extending high above the galactic plane, but it also has a significant extent in the galactic plane, in agreement with observations by Heiles & Jenkins (1976) and Heiles & Troland (2003). The maximum velocity in this region is km s, also in agreement with observations. Gas elements slide down the deformed magnetic field lines and undergo a phase transition, from the WNM to the CNM. A shock front marks the CNM-WNM interface, across which the density changes by a factor in pc. In addition to the large structures discussed above, there are small cold ‘cloudlets’ of size pc strung along field lines near the bottom of the panel. These resemble the cloudlets discovered by Heiles (1967) and further discussed in Heiles & Troland (2003). Figure 5b shows a typical cloudlet; the local magnetic field lines are also shown. The maximum number density of this cloudlet is cm; its major (minor) axis is (1.2) pc, giving an axis ratio of . The median column density is cm. These numbers are in good quantitative agreement with the observations discussed in § 8.3 of Heiles & Troland (2003).

There is also a great deal of CNM away from the magnetic valleys. In Figure 6, we show a cold Hi cloud elongated in the galactic plane at times (a) and (b) . There are notable density variations and ripples along the CNM-WNM interface, where there is significant velocity shear. The major (minor) axis of the cloud at is (1.35) pc, giving an axis ratio of . The magnetic field makes an angle of with respect to the minor axis of the cloud.

In light of the recent Hi surveys by Heiles & Troland (2003), it is of interest to calculate the fractions of interstellar mass in each thermal phase. In Figure 7, we show histograms of the fraction of the total mass in each thermal phase (CNM, unstable, and WNM) at times , , , 2, 7, and 12. Mass is binned by temperature in intervals of width . The vertical dashed lines denote the temperature boundaries separating the different thermal phases (, ). The mass fraction in each phase is shown in each figure as a percentage. The reference state () has all the mass in the WNM. As the Parker instability progresses, phase transitions are triggered, leading to a transfer of mass from the WNM, through the thermally-unstable phase, into the CNM. At the time when cold HI clouds have been formed and have settled into their final quasi-equilibrium states, the percentages of mass in the CNM and WNM are approximately 40% and 60%, respectively. (This number is found by simply averaging the results between and 12.) Enough time has elapsed for there to be an insignificant amount of mass in the thermally-unstable phase.

Taken at face value, this calculation predicts the existence of cold Hi clouds with maximum number density cm, maximum magnetic field strength G, and temperatures K. The mean values of the density and field strength in the cold clouds are found to be cm and G, respectively. The rms value of the magnetic field strength in the CNM is G, while the maximum line-of-sight value of the field occurs for a line-of-sight in the azimuthal direction in the galactic plane and is G. The (two-dimensional) turbulent velocity dispersion inside the CNM clouds is km s; this corresponds to a (thermal) Mach number . The mass of a typical cold cloud is calculated by assuming a thickness in the third (radial) direction corresponding to that of a galactic shock ( pc); we thus find the typical mass to be M. The spatial separation between cloud complexes is pc. An estimate of the maximum mass-to-flux ratio of typical cold Hi clouds yields , where is the critical mass-to-flux ratio calculated by Mouschovias & Spitzer (1976).

## 5 Summary and Discussion

We have performed two-dimensional simulations of the nonlinear evolution of the Parker instability properly accounting for phase transitions between the warm and cold neutral media. The evolution is initiated by small-amplitude velocity perturbations that have, and at later times lead to, evolution that preserves the symmetry of the dominant midplane-crossing (‘even’) mode — the mode that allows field lines to cross the galactic plane. This mode has been shown by BMP97 to be naturally selected when starting from a spectrum of random velocity perturbations. The qualitative features of the nonlinear evolution are similar to those found by calculations that ignored phase transitions: arches of magnetic field lines rise high above the galactic plane, accompanied by accumulations of mass in the resulting magnetic valleys. However, there are important quantitative differences between the two sets of calculations. Phase transitions, a consequence of realistic heating and cooling included in the present calculation, lead to structures (cold clouds) with typical densities much greater and temperatures much smaller than those achieved by isothermal calculations, and which are in excellent quantitative agreement with observations of cold Hi clouds.

The main result of this work is that the Parker instability is directly responsible for the formation of interstellar clouds and cloud complexes, despite recent claims to the contrary by Kim et al. (1998, 2001, 2002) and Santillán et al. (2000). Cold Hi clouds are seen to form whose masses ( M), mean magnetic field strengths (G), rms magnetic field strengths (G), mean densities ( cm), mass-to-flux ratios ( relative to critical), temperatures ( K), turbulent velocity dispersions ( km s), and spatial separations ( pc) are all in agreement with observations. The maximum number density is cm and the maximum magnetic field strength G. These physical conditions are such that, had we included the additional phase transition from atomic to molecular hydrogen, molecular clouds would have formed very rapidly within these giant Hi clouds.

We find that the fraction of ISM mass in the WNM is approximately 60%, in excellent agreement with observations by Heiles & Troland (2003). The bulk of the CNM has a sheet-like morphology. We also note the presence of cold ‘cloudlets’ of size pc strung along field lines. These cloudlets are both qualitatively and quantitatively similar to those first discovered by Heiles (1967). The median column density of a typical cloudlet is cm.

An apparent discrepancy between this work and the observational results of Heiles & Troland (2003) is the amount of WNM in the thermally-unstable phase. In the simulation, the maximum percentage of WNM in the unstable phase is . By contrast, the observed percentage is by Heiles & Troland (2003). The discrepancy is easily understood as a consequence of the fact that the observations see the ISM as it is today, after clouds not only have formed but also some of them have given birth to stars, including OB associations, which alter the physical conditions of their parent clouds as well as the amount of matter in the warm unstable phase. In other words, to account for the observed amount of gas in the thermally unstable phase, a calculation must include feedback from at least some of the energetic events that occur on relatively short timescales after clouds and stars have formed (e.g., stellar winds, UV radiation, or supernovae). Including such effects is beyond the scope of the present study, whose focus is cloud formation, not events that follow star formation.

Our results differ substantially from those of the recent three-dimensional simulations by Kosiński & Hanasz (2007). Important differences also exist in the formulation of the problem. Their neglect of thermal conductivity, together with their very low resolution ( compared to our ), results in thermally-unstable modes showing unresolved growth on their grid size of 10 pc (our grid resolution is pc). In addition, their assumed symmetry about the galactic plane limits the problem to the ‘odd’ mode, which is not favored during the nonlinear evolution over the ‘even’ mode we have followed. Their heating and cooling functions imply a stable CNM phase at an unrealistically-high temperature of K (depending on the model). For the realistic heating and cooling functions we have employed, these temperatures lie in the thermally-unstable regime; we find that the thermally stable CNM phase exists below a temperature of 185 K.

The results of this work lend strong support and quantify the scenario of cloud formation envisioned early on by Mouschovias et al. (1974). Taken in conjunction with the many successes of the theory of magnetic support of self-gravitating clouds and the ambipolar-diffusion–initiated fragmentation and star formation in such clouds (e.g., see § 4.2 of Mouschovias et al. 2006), the present work shows that magnetic fields play a crucial role not only in the formation of stars, but also in the formation and early evolution of their birthplaces, the interstellar clouds.

## Acknowledgments

We thank Shantanu Basu for providing the code developed for BMP97, which was used in an early phase of this work, as well as Carl Heiles and Konstantinos Tassis for useful discussions. TM acknowledges partial support from the National Science Foundation under grant NSF ASF-07-09206. All computations were performed on the Turing cluster (a 1536-processor Apple G5 X-serve cluster devoted to high-performance computing in engineering and science) maintained and operated by the Computational Science and Engineering Program of the University of Illinois at Urbana-Champaign.

## References

• Ames (1973) Ames S., 1973, ApJ, 182, 387
• Ballesteros-Paredes, Hartmann & Vázquez-Semadeni (1999) Ballesteros-Paredes J., Hartmann L., Vázquez-Semadeni E., 1999, ApJ, 527, 285
• Basu, Mouschovias & Paleologou (1997; hereafter BMP97) Basu S., Mouschovias T. Ch., Paleologou E. V., 1997, ApJ, 480, L55 (BMP97)
• Blitz & Shu (1980) Blitz L., Shu F. H., 1980, ApJ, 238, 148
• Elmegreen (2007) Elmegreen B. G., 2007, ApJ, 668, 1064
• Field & Saslaw (1965) Field G. B., Saslaw W. C., 1965, ApJ, 142, 568
• Field, Goldsmith & Habing (1969) Field G. B., Goldsmith D. W., Habing H. J., 1969, ApJ, L149
• Franco et al. (2002) Franco J., Kim J., Alfaro E. J., Hong S. S., 2002, ApJ, 570, 647
• García-Burillo, Guélin & Cernicharo (1993) García-Burillo S., Guélin M., Cernicharo, J., 1993, A&A, 274, 123
• Giz & Shu (1993) Giz A. T., Shu F. H., 1993, ApJ, 404, 185
• Hanasz & Lesch (2000) Hanasz M., Lesch H., 2000, 543, 235
• Hanasz & Lesch (2003) Hanasz M., Lesch H., 2003, A&A, 412, 331
• Hartmann, Ballesteros-Paredes & Bergin (2001) Hartmann L., Ballesteros­Paredes J., Bergin E. A., 2001, ApJ, 562, 852
• Hayes et al. (2006) Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188
• Heiles (1967) Heiles, C. 1967, ApJS, 15, 97
• Heiles & Jenkins (1976) Heiles, C., & Jenkins, E. B. 1976, A&A, 46, 333
• Heiles & Troland (2003) Heiles C., Troland T. H., 2003, ApJ, 586, 1067
• Heitsch et al. (2006) Heitsch F., Slyz A. D., Devriendt J. E. G., Hartmann L. W., Burkert A., 2006, ApJ, 648, 1052
• Hodge (1969) Hodge P. W., 1969, ApJ, 156, 847
• Kaplan & Pikel’ner (1970) Kaplan S. A., Pikel’ner S. B., 1970, The Interstellar Medium. Cambridge, Harvard University Press
• Kellman (1972) Kellman S. A., 1972, ApJ, 175, 363
• Kerr (1963) Kerr F. J., 1963, in Kerr, F. J., ed., The Galaxy and the Magellanic Clouds. Dordrecht, Reidel, 81
• Kim & Hong (1998) Kim J., Hong S. S., 1998, ApJ, 507, 254
• Kim & Ostriker (2006) Kim W.-T., Ostriker E. C., 2006, ApJ, 646, 213
• Kim, Hong & Ryu (1997) Kim J., Hong S. S., Ryu D., 1997, ApJ, 485, 228
• Kim et al. (1998) Kim J., Hong S. S., Ryu D., Jones T. W., 1998, ApJ, 506, L139
• Kim et al. (2000) Kim J., Franco J., Hong S. S., Santillán A., Martos M. A., 2000, ApJ, 531, 873
• Kim, Ryu & Jones (2001) Kim J., Ryu D., Jones T. W., 2001, ApJ, 557, 464
• Kim, Ostriker & Stone (2002) Kim W.-T., Ostriker E. C., Stone J. M., 2002, ApJ, 581, 1080
• Kosiński & Hanasz (2007) Kosiński R., Hanasz M., 2007, MNRAS, 376, 861
• Koyama & Inutsuka (2004) Koyama H., Inutsuka S., 2004, ApJ, 602, L25
• Kuwabara, Nakamura & Ko (2004) Kuwabara T., Nakamura K., Ko C. M., 2004, ApJ, 607, 828
• Kwan (1979) Kwan J., 1979, ApJ, 229, 567
• Mathewson & Ford (1970) Mathewson D. S., Ford V. L., 1970, Mem. RAS, 74, 143
• Mathewson, van der Kruit & Brouw (1972) Mathewson D. S., van der Kruit P. C., Brouw W. N., 1972, A&A, 17, 468
• Morgan (1970) Morgan W. W., 1970, in Becker W., Contropoulos G., eds, The Spiral Structure of Our Galaxy. Dordrecht, Reidel, 9
• Mouschovias (1974) Mouschovias T. Ch., 1974, ApJ, 192, 37
• Mouschovias (1975a) Mouschovias T. Ch., 1975a, PhD Thesis, Univ. Calif., Berkeley
• Mouschovias (1975b) Mouschovias T. Ch., 1975b, A&A, 40, 191
• Mouschovias (1996) Mouschovias T. Ch., 1996, in Tsinganos K. C., ed., NATO ASI Ser. C, Vol. 281, Solar and Astrophysical Magnetohydrodynamic Flows. Dordrecht, Kluwer, 475
• Mouschovias & Spitzer (1976) Mouschovias T. Ch., Spitzer L. Jr. 1976, ApJ, 210, 326
• Mouschovias, Shu & Woodward (1974) Mouschovias T. Ch., Shu F. H., Woodward P. R., 1974, A&A, 33, 73
• Mouschovias, Tassis & Kunz (2006) Mouschovias T. Ch., Tassis K., Kunz M. W., 2006, ApJ, 646, 1043
• Nagashima, Inutsuka & Koyama (2006) Nagashima M., Inutsuka S.-i., Koyama H., 2006, ApJ, 652, L41
• Oort (1954) Oort J., 1954, B.A.N., 12, 177
• Parker (1966) Parker E. N., 1966, ApJ, 145, 811
• Parker & Jokipii (2000) Parker E. N., Jokipii J. R., 2000, ApJ, 536, 331
• Penston & Brown (1970) Penston M. V., Brown F. E., 1970, MNRAS, 150, 373
• Piontek & Ostriker (2004) Piontek R. A., Ostriker E. C., 2004, ApJ, 601, 905
• Ryu et al. (2003) Ryu D., Kim J., Hong S. S., Jones T. W., 2003, ApJ, 589, 338
• Santillán et al. (2000) Santillán A., Kim J., Franco J., Martos M., Hong S. S., Ryu D., 2000, ApJ, 545, 353
• Scoville & Hersh (1979) Scoville N. Z., Hersh K., 1979, ApJ, 229, 578
• Scoville et al. (2001) Scoville N. Z., Polletta M., Ewald S., Stolovy S. R., Thompson R., Rieke M., 2001, ApJ, 122, 3017
• Shu (1974) Shu F. H., 1974, A&A, 33, 55
• van der Kruit (1973a) van der Kruit P. C., 1973a, A&A, 29, 231
• van der Kruit (1973b) van der Kruit P. C., 1973b, A&A, 29, 249
• van der Kruit (1973c) van der Kruit P. C., 1973c, A&A, 29, 263
• Westerhout (1963) Westerhout G., 1963, in F. J. Kerr, ed., The Galaxy and the Magellanic Clouds. Dordrecht, Reidel, 78
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters