GALACTIC SPIRAL SHOCKS WITH THERMAL INSTABILITY IN VERTICALLY STRATIFIED GALACTIC DISKS

Spiral arms in disk galaxies are regions of ongoing active star formation, sharply outlined by bright young star complexes. They usually span the entire optical disks and sometimes extend even to outer gaseous disks (e.g., Dickey et al. 1990; Boomsma et al. 2008; Bertin & Amorisco 2010, and references therein). Such large-scale spiral patterns may be the manifestation of spiral density waves which propagate with a constant pattern speed through stellar disks (Lin & Shu 1964, 1966), or may be transient features driven, for example, by tidal interactions with companion galaxies (e.g., Toomre & Toomre 1972; Hernquist 1990; Salo & Laurikainen 2000; Oh et al. 2008; Dobbs et al. 2010). Regardless of the origin of spiral features, it is widely accepted that the interstellar medium (ISM) is strongly compressed when it encounters stellar arms, forming narrow dust lanes in optical images (e.g., Elmegreen et al. 2006; Shetty et al. 2007). The densest parts of the shocked layers subsequently undergo gravitational collapse and produce downstream secondary structures, such as OB star complexes and giant H ii regions (e.g., Baade 1963; Elmegreen & Elmegreen 1983), filamentary gaseous spurs (also referred to as "feathers;" e.g., Scoville et al. 2001; Kennicutt 2004; Willner et al. 2004; La Vigne et al. 2006; Gordon 2007; Corder et al. 2008), and giant molecular associations or atomic superclouds (e.g., Elmegreen & Elmegreen 1983; Vogel et al. 1988; Rand & Kulkarni 1990; Koda et al. 2009).

Since shock compression within the arms is the first step toward star formation in grand-design spiral galaxies, understanding structural and dynamical evolution of these gas flows is essential to a host of fundamental problems, such as global star formation rates, the nature of the Hubble sequence, galaxy evolution, etc. Since the pioneering work of Roberts (1969) who obtained one-dimensional stationary shock profiles, there have been numerous studies of the structure of galactic spiral shocks under the simplifying assumption that the gas remains isothermal (e.g., Woodward 1975, 1976; Lubow et al. 1986; Martos & Cox 1998; Kim & Ostriker 2002; Gómez & Cox 2002, 2004; Wada & Koda 2004; Boley & Durisen 2006; Kim & Ostriker 2006; Dobbs & Bonnell 2006). In particular, Woodward (1975) and Kim & Ostriker (2002) showed that the one-dimensional shock profiles found by Roberts (1969) represent stable equilibria when the fluid quantities are allowed to vary only in the direction perpendicular to the arms. The growth of axisymmetric self-gravitating modes is limited by postshock expansion (Lubow et al. 1986).

When the direction parallel to the arm is included in models, on the other hand, isothermal spiral shocks in two dimensions are prone to various kinds of instabilities. Balbus (1988) showed that when self-gravity is included, spiral arms are unstable to transient swing instability. When magnetic fields are included, spiral arms are subject to magneto-Jeans instability, in which embedded parallel magnetic fields that exchange angular momentum limit the stabilizing effect of galaxy rotation, encouraging non-axisymmetric perturbations to grow into giant clouds and other arm substructures (Elmegreen 1994; Kim & Ostriker 2002, 2006; Shetty & Ostriker 2006). Wada & Koda (2004) showed that spiral shocks in two-dimensional thin disk models are unstable to a vorticity-generating wiggle instability and develop arm substructures (see also Dobbs & Bonnell 2006), although these in-plane modes appear to be suppressed by embedded magnetic fields (Shetty & Ostriker 2006; Dobbs & Price 2008) or by vertical shear and mixing when all the three dimensions are included in models (Kim & Ostriker 2006).

While steady in-plane shock solutions are subject to instability, shock models that include the vertical degree of freedom do not even have steady solutions. Instead, the shock front in vertically stratified disks moves back and forth relative to the mean position (Kim & Ostriker 2006). These shock "flapping" motions arise mainly because the vertical oscillation period of the gas is, in general, incommensurate with the arm-to-arm crossing time, so that the gas streamlines are not closed. Kim et al. (2006, hereafter Paper I) showed that the shock flapping is able to feed random gas motions on the scale of disk scale height that persist despite dissipation. The induced gas velocity dispersions reach a sonic level, suggesting that spiral shocks may be a considerable source of the ISM turbulence. Motions driven by shock flapping motions destroy any coherent vortical structures that would otherwise grow near the spiral shocks, suppressing the wiggle instability. Since gravity is a long-range force and insensitive to small-scale density structure, however, magneto-Jeans instabilities still grow within the arms in three-dimensional disk models, in spite of non-steady motions induced by shock flapping (Kim & Ostriker 2006).

Phase transitions caused by thermal instability (TI) create a multi-phase ISM, with important consequences for galactic star formation. In the classical picture of the ISM, TI changes an otherwise uniform ISM into warm rarefied material and cold dense clouds in a rough pressure balance (e.g., Field 1965; Field et al. 1969; Meerson 1996; Heiles 2001; Wolfire et al. 2003), while there also exists significant mass in the thermally unstable temperature range (e.g., Heiles & Troland 2003). Supernova blast waves create an additional, hot component that is organized into bubbles or cavities (Cox & Smith 1974; McKee & Ostriker 1977), although the total mass contained in the hot phase is much smaller than that in the cold and warm phases (e.g., Cox 2005). Cold atomic clouds transform to molecular clouds if their volume density (to produce molecules fast enough) and column density (to self-shield molecules against photodissociation) are sufficiently high (e.g., Elmegreen 1993; Draine & Bertoldi 1996), as in, for instance, shocks in turbulent flows (Glover & Mac Low 2007; Glover et al. 2010). That the star-forming molecular clouds strongly correlate with spiral arms (e.g., Stark 1979; Solomon et al. 1985; Kenney 1997; Zimmer et al. 2004; Shetty et al. 2007) suggests that spiral shocks too should trigger phase transitions from warm and diffuse to cold and dense conditions.

Effects of TI on spiral shocks were first studied by Shu et al. (1973), who calculated one-dimensional shock profiles consisting of two stable phases in equilibrium. Although they allowed for phase transitions, they assumed instantaneous thermal equilibrium, which precluded the existence of transitory thermally unstable gas in their calculations. Using direct time-dependent numerical simulations including ISM heating and cooling, Tubbs (1980) and Marochnik et al. (1983) found that spiral shocks trigger phase transitions if the initial density is large enough. Because of strong numerical diffusion associated with insufficient resolution, however, they were unable to capture TI in the postshock transition zone, which is the thermally unstable regime.

In Kim et al. (2008, hereafter Paper II), we used high-resolution one-dimensional simulations to study dynamical and thermodynamical evolution of gas flows across spiral arms with ISM heating, cooling, and thermal conduction. We found that even with TI, a quasi-steady state develops with the following recurring cycle: both warm and cold phases in the interarm region are shocked and immediately transform to denser cold gas in the arm, which subsequently evolves to be TI-unstable due to postshock expansion in a transition zone, and separates back into warm and cold phases. For a model with the initial number density of 2 cm^-3, the gas stays in the arm, transition, and interarm zones for 14%, 22%, and 64% of the arm-to-arm crossing time, respectively. The gas mass in the TI-unstable temperature range was ∼25%–30% of the total, and mostly located in the transition zone, suggesting that postshock expanding flows are important for producing intermediate-temperature gas. Paper II also found that TI in association with one-dimensional spiral shocks can drive random gas motions at ∼1.5 km s^-1 in the interarm and transition zones; these velocities are ∼5–7 times larger than those from pure TI alone (e.g., Kritsuk & Norman 2002; Piontek & Ostriker 2004).

In this paper, we extend the one-dimensional models of Paper II into two dimensions, in order to study the effect of vertical disk stratification on the dynamics and structure of multi-phase galactic spiral shocks. The current work also extends the vertically stratified isothermal models considered in Paper I to include the effects of the ISM heating and cooling. Our key objective is to find how the flapping motions of spiral shocks inherent in stratified disks interact with multi-phase gas produced by TI, to change the shock structure and drive random gas motions in each phase. We also study the internal properties of clouds that form in self-gravitating models. Although Dobbs & Bonnell (2007, 2008) and Dobbs & Price (2008) ran smoothed particle hydrodynamics (SPH) simulations to study shock structure and cloud formation in three dimensions, they used pre-determined cold and warm particles and did not allow the transitions between them. Dobbs et al. (2008) included ISM heating and cooling in the energy equation and thus handled TI self-consistently, focusing on the formation of molecular clouds in spiral shocks. Using grid-based three-dimensional simulations, Wada (2008) studied dynamics of galactic gas flows under the influence of self-gravity, a spiral-arm potential, radiative cooling, star formation, and energy feedback from supernova explosions. Although these three-dimensional global models are of course more realistic, our local models are useful for studying the detailed dynamics of spiral shocks at high spatial resolution, and allow us to isolate each effect of the physical processes involved.

This paper is organized as follows. In Section 2, we describe the basic equations we solve and specify the model parameters. In Section 3, we present the results of one-dimensional solutions for vertical disk equilibria including heating and cooling, also providing simple analytic expressions for the mass fractions and scale heights. In Section 4, we present the overall evolution, structure, and statistical properties of spiral shock flows with TI in stratified disks, based on the results of two-dimensional simulations. The effect of self-gravity is discussed in Section 5. In Section 6, we summarize our results and discuss their implications.

The local formulation used in the present study is similar to that in Papers I and II. In this section, we explain our numerical method and model parameters.

2.1. Basic Equations

We consider a local region centered on a spiral arm that is assumed to be tightly wound with a pitch angle sin i ≪ 1 and rotating at a constant pattern speed Ω_p. We set up a local Cartesian frame, centered at the position (R, ϕ, z) = (R₀, Ω_pt, 0), that corotates with the spiral arm. The x- and y-axes of the local frame are aligned in the midplane parallel and perpendicular to the local arm segment, while the z-axis points the direction perpendicular to the galactic plane (Roberts 1969; Paper I). Our simulation domain is a two-dimensional rectangular region with size L_x × L_z in the x–z plane (hereafter XZ plane). We assume that all physical variables are independent of y (quasi-axisymmetric), while allowing nonzero velocity in the y-direction in order to handle epicycle motions associated with galactic rotation self-consistently.

In this local frame, the galactic differential rotation is translated into the background velocity as

$$\begin{equation} \mathbf {v}_0=R_0(\Omega _0-\Omega _p)\sin i\hat{\mathbf {x}} + [R_0(\Omega _0-\Omega _p)-q_0\Omega _0 x]\hat{\mathbf {y}}, \end{equation} \tag{ 1 }$$

where Ω₀ = Ω(R₀) and $q_0\equiv -(d\ln \Omega /d\ln R)|_{R_0}$ denotes the local shear rate in the absence of the spiral potential (Kim & Ostriker 2002, 2006). Under the local approximation (i.e., L_x ≪ R₀ and |v| ≪ R₀Ω₀), the equations of hydrodynamics expanded in the local frame are

$$\begin{equation} \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf {v}_T)=0, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial \mathbf {v}_T}{\partial t}+\mathbf {v}_T\cdot \nabla \mathbf {v}_T= -\frac{1}{\rho }\nabla P -q_0\Omega _0 v_{0x}\hat{\mathbf {y}} -2\mathbf {\Omega }_0\times \mathbf {v}-\nabla (\Phi _s +\Phi _{\rm ext}), \end{equation} \tag{ 3 }$$

$$\begin{equation} \frac{\partial e}{\partial t} +\mathbf {v}_T\cdot \nabla e = -\frac{\gamma }{\gamma -1}P\nabla \cdot \mathbf {v}_T -\rho \mathcal {L}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \nabla ^2\Phi _s=4\pi G\rho, \end{equation} \tag{ 5 }$$

(e.g., Roberts 1969; Shu et al. 1973; Kim & Ostriker 2006), where v_T ≡ v₀ + v is the total velocity in the local frame, Φ_s is the self-gravitational potential of the gas, Φ_ext is the external stellar potential, and $\rho \mathcal {L}(\rho, T)$ is the net cooling function. Other symbols have their usual meanings. We adopt an ideal gas-law P = (γ − 1)e with γ = 5/3.

The external stellar potential Φ_ext varies in both the x- and z-directions and is separable into two parts as

$$\begin{equation} \Phi _{\rm ext}=2\pi G \rho _* z^2 +\Phi _{\rm sp}\cos \left(\frac{2\pi x}{L_x}\right), \end{equation} \tag{ 6 }$$

where ρ_* is the unperturbed midplane stellar density, Φ_sp is the amplitude of the spiral potential, and L_x = 2πR₀sin i/m is the arm-to-arm separation for an m-armed spiral. The first term in Equation (6) adopts a constant density ρ_* for the stellar disk vertically; this is a reasonable assumption given that most of the gas is located within one stellar scale height of the midplane. The second term is a local analog of a logarithmic spiral potential considered in Roberts (1969) and Shu et al. (1973). To parameterize the spiral arm strength, we employ the dimensionless parameter

$$\begin{equation} F\equiv \frac{m}{\sin i}\left(\frac{|\Phi _{\rm sp}|}{R_0^2\Omega _0^2}\right), \end{equation} \tag{ 7 }$$

corresponding to the ratio of the maximum force due to the spiral potential to the mean radial gravitational force (Roberts 1969).

The interstellar gas is subject to the net cooling $\rho \mathcal {L} \equiv n^2\Lambda [T]-n\Gamma$, where n = ρ/(μm_H) is the gas number density with μ = 1.27 being the mean molecular weight per particle. For the cooling rate of the diffuse ISM, we take the fitting formula suggested by Koyama & Inutsuka (2002),

$$\begin{eqnarray} \Lambda (T) &=& 2\times 10^{-19}\exp \left(\frac{-1.184\times 10^5}{T+1000}\right) \nonumber\\ && + 2.8\times 10^{-28}\sqrt{T}\exp \left(\frac{-92}{T}\right)\,{\rm erg} \,{\rm cm}^3 \,{\rm s^{-1}}, \end{eqnarray} \tag{ 8 }$$

with temperature T in degrees Kelvin.

For the gas heating function, we consider two different cases: (1) a constant heating rate Γ = Γ₀ = 2.0 × 10⁻²⁶ erg s^-1 and (2) a density-dependent heating rate

$$\begin{equation} \Gamma = \Gamma _0\exp (n/ n_0)^3, \end{equation} \tag{ 9 }$$

with n₀ = 20 cm^-3. The first, uniform heating rate corresponds to the mean heating rate due to the photoelectric effect on small grains and polycyclic aromatic hydrocarbons (PAHs), photodissociations of hydrogen molecules, and ionizations by cosmic rays and X-rays (e.g., Koyama & Inutsuka 2002). The second, density-modified heating rate that increases stiffly with n is to treat the effect of star formation feedback in a very simple way. Without such a feedback, high-density clouds produced inside spiral arms in our self-gravitating models would undergo catastrophic collapse, preventing further computation. In the real ISM, gravitational collapse leads to new-born stars which will in turn disperse the parent clouds by injecting radiative and mechanical energies. Investigating the details of star formation and consequent feedback processes is a very active current research area. Most previous works have adopted simplified feedback prescriptions that depend on gas-consumption rate, star formation efficiency, type of energy injection, etc., with considerable uncertainties in the parameters (e.g., Springel et al. 2005; Joung & Mac Low 2006; Shetty & Ostriker 2008; Koyama & Ostriker 2009a). More realistic feedback prescriptions will be considered in a subsequent paper.

Figure 1 plots the equilibrium cooling curves in the density versus pressure plane. The solid line corresponds to the uniform heating rate, while the dashed curve is for the modified heating rate. The dotted lines indicate isotherms. The modified heating rate changes the equilibrium curve dramatically only for high-density gas, while making a negligible difference for low-density material. The equilibrium pressure has a local maximum P_max/k_B = 5.0 × 10³ K cm^-3 at n_crit,1 = 1.0 cm^-3 for both heating rates, while attaining a local minimum at n_crit,2 = 8.7 cm^-3 with P_min/k_B = 1.6 × 10³ K cm^-3 for the constant heating rate, and at n_crit,2 = 6.9 cm^-3 with P_min/k_B = 1.7 × 10³ K cm^-3 for the modified heating rate. Under the constant heating rate, the gas temperature along the equilibrium curve is a monotonically decreasing function of density, although it is insensitive to n at the low-density end with n < n_crit,1. This is not the case under the density-modified heating rate where the equilibrium temperature increases with density when n>n_crit,3 = (1/3)^1/3n₀ = 13.9 cm^-3 in order to model feedback. We thus classify the gas into two components based on its density rather than temperature: rarefied component if n < n_crit,1 and dense component if n>n_crit,1. Note that thermally unstable gas with n_crit,1 < n < n_crit,2 belongs to the dense phase according to our classification.

Zoom In Zoom Out Reset image size

As Equation (4) indicates, we explicitly ignore the effect of thermal conduction in the present work. Paper I found that large translational motions in a finite difference scheme give rise to numerical diffusion that tends to suppress the growth of TI, similarly to thermal conduction. The amount of numerical conductivity in our models is typically $\mathcal {K}_n = 10^9 \,{\rm erg\,s^{-1}}\,{\rm cm}^{-1}\,{\rm K}^{-1}$ for the background velocity v_0x = 13 km s^-1, grid spacing Δx = 2.5  pc, and the perturbation wavelength λ = 20  pc. Inclusion of physical conductivity larger than $\mathcal {K}_n$ would resolve the wavelengths of the most unstable TI. But, this would in turn reduce the time step greatly, making computation essentially unpractical.³ We note that by neglecting the thermal conduction term in the energy equation, some of our results may depend on the numerical resolution, although the mass fractions appear to be insensitive to the resolution (Paper I).

2.2. Model Parameters and Numerical Methods

We consider a simulation box located near the Solar neighborhood at galactocentric radius of R₀ = 8 kpc. We adopt the galactic rotational velocity of R₀Ω₀ = 208 km s^-1 with a flat rotation curve (q₀ = 1). The corresponding angular velocity is Ω₀ = 26 km s^-1 kpc⁻¹, and orbital period is t_orb ≡ 2π/Ω₀ = 2.4 × 10⁸ yr, which we use as the time unit in our presentation. For spiral arm parameters, we take pattern speed Ω_p = 0.5 Ω₀, pitch angle sin i = 0.1, and azimuthal wavenumber m = 2. The corresponding arm-to-arm separation is L_x = 2πR₀sin i/m = 2.5 kpc, which is set equal to the size of the simulation box along the x–direction. We fix the spiral arm strength to F₀ = 5%.

Our initial gaseous disks, in the absence of the spiral-arm perturbations, are taken to be isothermal and in vertical hydrostatic equilibrium under the linear stellar gravity g_z = −4πGρ_*z (see Equation (6)). The corresponding density distribution is a Gaussian profile

$$\begin{equation} \rho (z)=\rho _0 \exp \left(-\frac{z^2}{2h_g^2}\right), \end{equation} \tag{ 10 }$$

with a scale height

$${\selectfont\fontsize{9}{10}{\begin{equation} h_g=\frac{c_R}{(4\pi G\rho _*)^{1/2}}= 128{\,\rm \,pc}\left(\frac{c_R}{7\,{\rm km\,s^{-1}}}\right)\left(\frac{\rho _*}{0.056\,{M_\odot }{\,\rm \,pc}^{-3}}\right)^{-1/2}, \end{equation}}} \tag{ 11 }$$

where c_R = 7 km s^-1 is the isothermal sound speed of the initial disks and ρ_* = 0.056 M_☉  pc⁻³ is the stellar density near the solar neighborhood (Holmberg & Flynn 2000). We take L_z = 7.5 h_g = 960  pc as the vertical size of the simulation domain (i.e., |z| ⩽ L_z/2).

Tables 1 and 2 summarize the model parameters and some simulation outcomes for models with and without spiral potential perturbations, respectively. Column 1 labels each run. The prefixes NU and NM stand for non-self-gravitating models ("N") with the uniform heating rate ("U") and the modified heating rate ("M"), respectively, while the prefix SM indicates self-gravitating models ("S") with the modified heating rate ("M"). As will be discussed below, Column 2 gives the initial gas surface density Σ₀. Columns 3 and 4 give the mass fractions, f_a ≡ 〈∫ρ_adxdz/∫ρdxdz〉 (with a = D or R), of dense and rarefied components, respectively. Here, the angle brackets 〈〉 denote a time average over t/t_orb = 5–8 for non-self-gravitating models and over t/t_orb = 8–11 for self-gravitating models. Columns 5 and 6 give the scale heights, H_a ≡ 〈∫ρ_az²dxdz/∫ρ_adxdz〉^1/2, of the dense and rarefied components, respectively. Column 7 gives the average scale height of the whole gas H_ave ≡ (f_DH²_D + f_RH²_R)^1/2.

Table 1. Models Without Spiral Arms (F = 0%)

Modela	$\Sigma _0\; ({ M_\odot \,\rm pc^{-2}})$	fD(%)	fR(%)	$H_{ D}\;(\rm pc)$	$H_{ R}\; (\rm pc)$	$H_{\rm ave}\; (\rm pc)$
(1)	(2)	(3)	(4)	(5)	(6)	(7)
NU.S02	2	0	100	0	126	126
NU.S05	5	0	100	0	119	119
NU.S10	10	71	29	2	125	67
NU.S20	20	86	14	4	127	48
NM.S02	2	0	100	0	126	126
NM.S05	5	0	100	0	119	119
NM.S10	10	69	31	4	125	70
NM.S20	20	85	15	9	130	50
SM.S02	2	0	100	0	121	121
SM.S05	5	0	100	0	107	107
SM.S10	10	82	18	4	100	43
SM.S20	20	94	6	7	84	21

Note. ^aThe prefixes NU refers to the non-self-gravitating models with the uniform heating rate, NM for the non-self-gravitating models with the modified heating rate, and SM for the self-gravitating models with the modified heating rate.

Download table as:  ASCIITypeset image

Table 2. Models With Spiral Arms (F = 5%)

Model	$\Sigma _0\; ({ M _\odot \,\rm pc^{-2}})$	fD(%)	fR(%)	$H_{ D}\; (\rm pc)$	$H_{ R}\; (\rm pc)$	$H_{\rm ave}\; (\rm pc)$
(1)	(2)	(3)	(4)	(5)	(6)	(7)
NU.S02	2	11	89	25	129	122
NU.S05	5	62	38	10	130	81
NU.S10	10	81	19	7	130	57
NM.S02	2	12	88	28	130	122
NM.S05	5	60	40	20	132	84
NM.S10	10	81	19	25	135	64
SM.S02	2	26	74	14	124	107
SM.S05	5	91	9	21	121	42
SM.S10	10	95	5	43	123	51

Note. Model name prefixes are as in Table 1.

Download table as:  ASCIITypeset image

We integrate the time-dependent partial differential Equations (2)–(5) using a modified version of the Athena code (Gardiner & Stone 2005, 2008; Stone et al. 2008; Stone & Gardiner 2009). Athena employs a single-step, directionally unsplit Godunov scheme for magnetohydrodynamics in multi-spatial dimensions. Among the various schemes contained in Athena, we take a piecewise linear method for spatial reconstruction, HLLE Riemann solver to compute the fluxes (Harten et al. 1983; Einfeldt et al. 1991), and van Leer algorithm for directionally unsplit integration (Stone & Gardiner 2009). Since our simulations involve strong shocks for the dense medium (with typical Mach numbers ∼7–10), we adopt the first-order flux correction when the net mass flux out of a cell exceeds the initial mass of the cell in order to avoid an occurrence of negative density (see, e.g., Lemaster & Stone 2009). Our models employ a 1024 × 512 zones over the simulation box, corresponding to the resolution of Δx = 2.4  pc and Δz = 1.9  pc.

We adopt the shearing-periodic boundary condition at the x-boundaries (Hawley et al. 1995). In the z-direction, we use the outflow condition for the velocity and the vacuum condition for the gravitational potential (e.g., Koyama & Ostriker 2009a). For the density and pressure at the z-boundaries, we linearly extrapolate the logarithmic density, while keeping temperature fixed, whenever dρ/dz < 0. This produces a balance between the vertical pressure gradient and the gravitational source term at the boundaries, similarly to the "conducting" boundary in Parrish & Stone (2005). When dρ/dz>0, on the other hand, we switch to the continuous boundary condition for the density and pressure to reduce artificial mass inflow due to the extrapolation. Under our boundary conditions, the gas can freely escape from the vertical boundary; we have checked that the total mass is nonetheless conserved within 2% for all models. Because of the very short cooling time, energy updates from the net cooling functions are made implicitly based on Newton–Raphson iterations (Piontek & Ostriker 2004). To solve for the gravitational potential in our simulation domain, we adopt a method introduced by Koyama & Ostriker (2009a) which, by using the fast Fourier transform technique, is much more efficient than a hybrid method involving Green's functions (e.g., Kim et al. 2002).

While our main objective is to investigate the overall dynamics and structure of spiral shocks in vertically stratified disks under the influence of TI, in this section we focus on the quasi-static vertical equilibria with heating and cooling in the absence of the spiral arm potential (i.e., F = 0). This allows us to study the effect of TI on vertical disk structure. We run one-dimensional simulations with physical quantities varying only with z. We consider an initially isothermal disk with Σ₀ = 2, 5, 10, or 20 M_☉ pc⁻², and evolve it subject to TI.

For disks with large surface density (models with Σ₀ ⩾ 10 M_☉ pc⁻²), TI grows rapidly (≪t_orb), transforming the initially constant-temperature gas into thermally bistable phases. The cold, dense gas falls toward the midplane to form a dense slab, while the warm, rarefied gas rises up buoyantly. The infall is supersonic relative to the dense medium. At early time, the dense slab surrounded by the upper rarefied gas undergoes vertical expansions and contractions a few times. As the kinetic energy dissipates through shocks at the interfaces, the whole configuration evolves toward vertical hydrostatic equilibrium typically within ∼0.6t_orb. Figure 2 shows density profiles for S10 models with Σ₀ = 10 M_☉ pc⁻². Solid and dotted lines are for non-self-gravitating NU and NM models, respectively, while the dashed line is for the self-gravitating SM models. The difference between models NU.S10 and NM.S10 is not significant since the maximum midplane density is not much larger than n₀ = 20 cm^-3, below which the heating rate is almost density independent. For model SM.S10, self-gravity compresses the midplane slab further at the expense of the rarefied medium at |z|>H_D. Nevertheless, the phase transition between dense and rarefied components turns out to occur at almost the same density n_trans ∼ 0.25–0.35 cm^-3, corresponding to the transition pressure P_trans/k_B ∼ 2000–2200 cm^-3 K for all models that are unstable to TI. Note that P_trans/k_B is above the minimum pressure for a cold medium with our adopted cooling and heating functions, P_min/k_B = 1600–1700 cm^-3 K.

Zoom In Zoom Out Reset image size

Once vertical hydrostatic equilibrium is attained, we measure the mass fractions f_D and f_R, and the scale heights H_D and H_R of the dense and rarefied components, respectively; these values are listed in Table 1. Figure 3 plots f_R and H_R as functions of the initial surface density Σ₀. The results of NU and NM models are denoted by asterisks and diamonds, respectively, while open squares are for SM models. Models with low surface density (Σ₀ = 2 M_☉ pc⁻²) do not experience TI and thus establish a single-phase equilibrium consisting only of the rarefied medium. Since the warm gas is nearly isothermal at c_R ≈ 7 km s^-1 and self-gravity is weak in these models, the equilibrium density profiles in low Σ models are approximately given by Equation (10), with surface density $\Sigma = \rho _0 h_g \sqrt{2 \pi }=\rho _0 c_R /\sqrt{2 G \rho _*}$. Since the midplane pressure P_R(0) = c²_Rρ₀ of the rarefied gas cannot exceed P_max along the thermal equilibrium curve shown in Figure 1, the surface density for a single-phase equilibrium with only a rarefied component is constrained to be less than Σ_max = P_max/(2Gρ_*c²_R)^1/2 = 6.7 M_☉ pc⁻². Similarly, the condition P_R(0) = P_min yields Σ_min = 2.1 M_☉ pc⁻² as the minimum surface density for a two-phase equilibrium in which dense and rarefied components coexist. The two vertical dot-dashed lines in Figure 3 mark Σ_min and Σ_max. Below Σ_min, only a rarefied phase is possible, whereas above Σ_max, both dense and rarefied phases must be present.

Zoom In Zoom Out Reset image size

For Σ_min < Σ₀ < Σ_max, both single (rarefied) phase and two-phase equilibria can be realized. Which type of equilibrium emerges depends, of course, on the initial disk conditions. In the case of our models with Σ₀ = 5 M_☉ pc⁻², the initial midplane density and pressure are n(0) = 0.5 cm^-3 and P(0)/k_B = 3770 K cm^-3, smaller than than n_crit,1 and P_max/k_B. Since cooling and heating occur almost isobarically, even the densest gas in these models is unable to overcome P_max to turn into the dense component, for this case.

Figure 3 also shows that for the models that reach a two-phase equilibrium, self-gravity reduces the rarefied-gas fraction in mass as well as its scale height compared to those in non-self-gravitating counterparts. Self-gravity also makes the density profile of the rarefied component deviate significantly from a Gaussian profile. A thin midplane dense slab, containing the majority of the gas mass, exerts a uniform gravity on the rarefied gas lying above it, providing an additional confining force. In the Appendix, we describe a simple way to estimate f_R and H_R as functions of the total gas density, assuming that the rarefied component can be characterized by a fixed sound speed c_R and that its self-gravity is negligible. The resulting theoretical predictions, with c_R = 7 km s^-1, for self-gravitating and non-self-gravitating cases are plotted in Figure 3 as dashed and solid curves, respectively. These are overall in good agreement with the numerical results. Small discrepancies between the theoretical and numerical values of H_R for disks with Σ₀ = 5 M_☉ pc⁻² arise from the fact that the rarefied gas in these models has larger midplane pressure than in any other models.⁴ In view of the thermal equilibrium curve shown in Figure 1, this implies that the rarefied medium in S05 models is coldest, corresponding to c_R ≃ 6.3 km s^-1, making the scale height smaller than the theoretical estimate based on c_R = 7 km s^-1.

Now we turn to our main theme, nonlinear gas flows with TI across spiral arms in a stratified disk. In this section, we study overall evolution, structure, and statistical properties such as density and random velocity distributions of spiral shocks for non-self-gravitating models. Effects of self-gravity will be discussed in the next section.

4.1. Overall Evolution

We begin by describing evolution of our fiducial models NU.S10 and NM.S10 with Σ₀ = 10 M_☉  pc⁻² that employ the uniform and density-modified heating rates, respectively. We slowly turn on the spiral potential amplitude such that it attains full strength F₀ = 5% at t/t_orb = 1.5. Snapshots of volume density in logarithmic color scale at early epochs t/t_orb = 1.50, 2.24, and 2.40 are shown in Figures 4 and 5, respectively. Figure 6 plots the gas distribution in the n-P/k_B plane for both models at t/t_orb = 2.40. Initially, the disk is in hydrostatic equilibrium with a constant sound speed of c_R = 7 km s^-1. Since the initial disk is out of thermal equilibrium, it quickly evolves into a two-phase equilibrium configuration, as explained in Section 3. As F increases, both the dense gas near the midplane and the rarefied gas at high |z| respond to the growing spiral potential and form a shock front near the potential minimum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the gas flows at this time are fairly horizontal without much vertical mixing, as evidenced by the instantaneous streamlines shown in Figures 4(a) and 5(a), the shock profile at each height is very similar to those of the one-dimensional cases studied in Paper II. The shock strength and gas phase in the postshock regions depend on the mean density and temperature at that height. Near the midplane at |z| < H_D (=7 and 25  pc for models NU.S10 and NM.S10, respectively), the dense slab is so cold that the shock is very strong with a typical Mach number $\mathcal {M}\sim 10$, resulting in a far denser postshock region. In the high-|z| regions (|z|>130  pc), on the other hand, the gas is warm and has a low mean density (<0.1 cm^-3) enough to remain warm even after the shock compression. It is the mid-altitude rarefied medium (at H_D < |z| < 130  pc) that is able to achieve a postshock pressure larger than P_max and thus undergoes a phase transition to the dense component after experiencing isobaric cooling (Mufson 1974; Inoue & Inutsuka 2008; Paper II). Since the shock is stronger at lower |z| in a stratified disk and since a stronger shock tends to move downstream (e.g., Kim & Ostriker 2002), the shock front when it first develops is naturally curved in the XZ plane. Figures 4(a) and 5(a) show that the shock front is concave when seen from the upstream direction, with mean slopes of |dx_sp(z)/dz| ≈ 0.83 at |z| < 130  pc and 0.13 at |z|>130  pc, where x_sp(z) is the shock location at z.

The dense gas produced at the shock at moderate z begins falling toward the midplane under the influence of the external gravity as it moves downstream. The reduction of the velocity in the direction normal to the concave shock front also helps the downward motion of the gas. On the other hand, the dense gas near the midplane has a large pressure and thus slightly expands vertically after the shock. The vertical expansion is more extreme in NM models than in NU models. The falling gas collides with the expanding gas, reducing the rising motion of the latter. The streamlines shown in Figures 4(a) and 5(a) illustrate these motions at early time.

The rarefied gas which crosses the shock at high |z| also falls toward the midplane as it follows galaxy rotation. This builds up thermal pressure at low |z|, so the flow rebounds to high-altitude regions. Since the period of vertical oscillation, ∼(Gρ_*)^−1/2, is in general not commensurate with the interval between arm crossings, the streamlines of the rarefied gas are not closed. This causes the shock front to sway back and forth around its mean position in the direction perpendicular to the arm (e.g., Kim & Ostriker 2006; Paper I). During the course of the flapping motions, the shock front has a convex shape (seen from upstream) when the postshock regions are maximally compressed (Figure 4(b)), while it becomes concave when the gas in the postshock regions is in full vertical expansion (Figure 4(c)). These flapping motions of the shock front, alternating between convex and concave shapes, occur quasi-periodically with a period of ∼0.5 t_orb and have an amplitude of Δx/H_R ∼ 1 at |z| = H_R (=130  pc in model NU.S10). The shock flapping motions are able to tap some of the kinetic energy in galaxy rotation to supply random kinetic energy for the gas. We will quantify the amplitudes of random gas motions driven by flapping in Section 4.3.

One of the special features of galactic spiral shocks is that gas experiences acceleration after the maximum shock compression, forming a postshock expansion zone (e.g., Balbus 1988; Kim & Ostriker 2006; Papers I & II). Any parcel of gas becomes gradually less dense as it moves downstream in the expansion zone. In model NU.S10, the shock compression and subsequent cooling is so strong that the shocked dense gas in the midplane can reach n>10³ cm^-3 (see also model SC20 in Paper II). With such a large postshock density, this gas can still remain dense, with n>n_crit,2, even after emerging from the expansion zone located at x/L_x ∼ 0–0.3. This TI-stable dense gas travels almost ballistically in the interarm region, reenters the arm, and combines with other dense gas to produce a few condensations. Figures 4(b) at time t/t_orb = 2.24 shows three large condensations located at x/L_x ∼ −0.3, 0, and 0.14, which are stretched horizontally due to the expanding background velocity. The condensations move nearly horizontally to appear at x/L_x = −0.12, 0.14, and 0.35 when t/t_orb = 2.40 (Figure 4(c)). In model NU.S10, the dense condensations stay in the arm (−0.05 < x/L_x < 0.05) for about ∼0.15 t_orb, in the expansion zone (0.05 < x/L_x < 0.25) for about ∼0.30 t_orb, and in the interarm region for the remainder (∼0.55 t_orb) of the cycle. The width and residence time of each zone are insensitive to the model parameters.

With the density-modified heating rate, on the other hand, the postshock dense gas in model NM.S10 has a moderate density (∼30–40 cm^-3), so that the postshock expansion is able to take it to the thermally unstable regime (n_crit,1 < n < n_crit,2). Subsequently, the expanding dense gas suffers from TI and separates back into dense and rarefied gas in the interarm region. The large thermal pressure also prevents the formation of dense condensations in this model.

Figure 7 plots the temporal evolution of the mass fractions of dense phase (solid lines) and rarefied phase (dashed lines), respectively, for models NU.S10 (thick lines) and NM.S10 (thin lines). At early time, f_D increases rapidly as the gas cools and collapses toward the midplane to form a dense layer that bounces appreciably at t/t_orb ∼ 0.1. The mass fractions flatten at t/t_orb ∼ 0.6 when vertical hydrostatic equilibrium is established, well before the effect of the spiral potential becomes substantial. As the spiral potential attains its full strength at t/t_orb = 1.5, f_D increases slightly due to the phase transition of the rarefied to dense phases occurring at the shock fronts. Although the flows are fully nonlinear with strong unsteady motions and phase transitions, there is no noticeable secular variation in the mass fractions, which remains at f_D ∼ 0.8 after t/t_orb = 5; the associated temporal fluctuation amplitudes are about 6%–9% relative to the mean values. We thus conclude that in a statistical sense, the spiral shocks in our models have reached a quasi-steady state at t/t_orb>5. Compared with models without spiral arms discussed in Section 3, the shock compression and associated phase transitions decrease the rarefied gas fraction by 46% for S10 models. In fact, all of the non-self-gravitating models with spiral arms have comparable total surface density of rarefied gas, ∼1.9 M_☉  pc⁻², lower than the value $\Sigma _R= P_{\rm trans}/(c_R \sqrt{2 G \rho _*})\approx 2.8 \,{M_\odot }{\,{\rm \,pc}}^{-2}$ that would be predicted using a uniform surface density. Note that both NU and NM models have almost the same dense and rarefied mass fractions since the modified heating rate does not affect the rarefied medium much.

Zoom In Zoom Out Reset image size

The evolution of S02 and S05 models is qualitatively similar to that of S10 models in that phase transitions occur at the shock and in the postshock expansion zone, although the former with low postshock density do not produce much dense gas even under the uniform heating rate. When the spiral potential is absent, the equilibrium disks of these models consist entirely of the rarefied gas with the midplane pressure P(0)/k_B ∼ 3500–4000 cm^-3 K for S05 models and P(0)/k_B ∼ 1500–2000 cm^-3 K for S02 models. But, the shock compression increases the pressure by about a factor of 3, corresponding to a typical Mach number $\mathcal {M}\approx v_{0x}/c_R\sim 2$ for the rarefied medium,⁵ making the midplane postshock pressure larger than P_max. As a result, the dense medium in S05 models comprises about 60% of the total mass, undergoing TI in the postshock expansion zone. In S02 models, the postshock pressure barely exceeds P_max, so that the shocked dense gas, comprising about 10% of the total, easily disperses to return to the rarefied gas in the interarm region. Flapping motions of spiral shocks are correspondingly weaker in these models, with Δx/H_R ∼ 0.7 and 0.3 for S05 and S02 models, respectively.

4.2. Time-averaged Shock Structure

To visualize spiral shock structure in each model, we construct a time-averaged distribution of number density 〈n〉. Here, the angle brackets 〈〉 denote a time average over t/t_orb = 5–8. Figure 8 displays 〈n〉 for all the non-self-gravitating models in logarithmic color scale. It is apparent that S05 and S10 models possess a thin dense layer everywhere near the midplane, while the dense gas is found only inside the arm regions in S02 models. The shock compression factors in the time-averaged density profiles are ∼7–10, which is larger than the adiabatic shock jump due to enhanced radiative cooling in the shocked gas (cf., Mufson 1974; Inoue & Inutsuka 2008; Paper II). The shock transition layer in S05 and S10 models is relatively broad because of rather strong flapping motions of the shocks, while S02 models exhibit relatively sharp discontinuities. Compared to NU models, arms in NM models have larger pressure and are more expanded vertically, similar to "hydraulic jumps" that occur when the equation of state is stiffer than isothermal (e.g., Martos & Cox 1998). Table 2 lists the time-averaged values of the mass fractions and scale heights of dense and rarefied components, as well as the overall average scale height.

Zoom In Zoom Out Reset image size

Figure 9 plots the mass-weighted probability distribution functions (PDFs), averaged over t/t_orb = 5–8, of gas density and temperature for models NU.S10 and NM.S10. The PDFs are in general bimodal, as is expected for a bistable cooling function. For model NU.S10, the dense and rarefied peaks are centered at (n, T) ∼ (200 cm^-3, 30 K) and ∼(0.2 cm^-3, 7100 K), respectively, mostly distributed near the thermal equilibrium curves. Only a small fraction of the gas is in the thermally unstable regime. The dense portion of the PDF in model NU.S10 is well fitted by a lognormal distribution (thin solid line) with a standard deviation of Δ(ln n) = 1.2, which is one of the characteristics of near-isothermal turbulence (cf., Wada & Norman 2007; Wada 2008). With enhanced heating, on the other hand, model NM.S10 shows a sharp density cutoff in the density PDF at n ∼ 50 cm^-3 and has a dense peak shifted to n ∼ 20 cm^-3. Because of the stiff equation of state, the dense gas in model NM.S10 is not as cold as in model NM.S10. This not only thickens the midplane dense slab, but also sets an upper limit on the gas density, which in turn prevents the formation of dense condensations. In model NM.S10, all the postshock dense gas becomes thermally unstable in the expansion zone and separates into dense and rarefied gas in the interarm region.

Zoom In Zoom Out Reset image size

Paper II showed that for one-dimensional models, the density profiles of multi-phase spiral shocks are more symmetric and have a wider arms than isothermal counterparts. This is because the strength of spiral shocks in the multi-phase models fluctuates depending on whether the incoming gas is warm or cold, resulting in slight oscillations of the shock fronts in the direction perpendicular to the arm. In addition, spiral shocks in the XZ plane undergo flapping motions, which can further widen the arms. To see this, we plot in Figure 10 the time-averaged surface density profiles Σ(x) = ∫^∞_−∞〈ρ〉dz after taking a boxcar average with window of 8  pc. The solid and dotted lines are for NU and NM models, respectively. Shown also as dashed lines are the density profiles n(x)/n₀ from one-dimensional simulations (i.e., without vertical stratification) under the uniform heating rate; the initial number density n₀ of the one-dimensional counterpart was chosen equal to the density-weighted mean density n_ave = Σ₀/(2π^1/2μm_HH_ave), with the average disk thickness H_ave listed in Column 7 of Table 2. For S05 and S10 models for which the shock flapping motions are appreciable, the arms are considerably wider and less centrally peaked than in the one-dimensional models. Due to the flapping motions, dense condensates formed in the NU.S10 model oscillate slightly in the x-direction when they pass through the shock, resulting in broader arms than in the NM.S10 model. For S02 models, the one-dimensional shock consists only of the warm rarefied gas, while the two-dimensional shocks contain a small amount (∼10%) of dense gas due to an additional compression in the z-direction; the difference in profiles is therefore slight.

Zoom In Zoom Out Reset image size

In studies of galactic disk structure, it has been the customary to assume effective hydrostatic equilibrium in the vertical direction. Using numerical simulations without spiral arms, Koyama & Ostriker (2009b) explicitly demonstrated that turbulent, multi-phase disks are in effective hydrostatic equilibrium, provided that the turbulent pressure arising from random gas motions is taken into account. When a spiral potential is present, the gas surface density and velocity dispersions depend upon the distance x from the minimum of the spiral potential. It is interesting to study whether "local" effective hydrostatic equilibrium is still established at each x.

From the time-averaged density distribution, we measure the density-weighted vertical scale height H(x), sound speed c_s(x), and vertical velocity dispersion δu_z(x) via

$$\begin{eqnarray} H^2(x) &=& \frac{\int \langle \rho \rangle z^2dz}{\int \langle \rho \rangle dz}, \quad c_{s}^2(x) =\frac{\int \langle P \rangle dz}{\int \langle \rho \rangle dz}, \quad {\rm and} \nonumber\\ \delta u_{z}^2(x) &=& \frac{\int \langle \rho [v_z-\langle v_z\rangle ]^2\rangle dz}{\int \langle \rho \rangle dz}. \end{eqnarray} \tag{ 12 }$$

In the absence of self-gravity, the "estimated" vertical scale height is given by

$$\begin{equation} H_{\rm est}^2(x)=\frac{c_s^2 + \delta u_{z}^2}{4\pi G\rho _*} \end{equation} \tag{ 13 }$$

for effective hydrostatic equilibrium (Koyama & Ostriker 2009b). Figure 11 plots H(x) (solid lines) and H_est(x) (dotted lines) for NU and NM models as functions of x. The measured vertical scale height is overall in excellent agreement with the estimated value at all horizontal locations. For all models, c_s is about 5–7 times larger than δu_z. This implies that the disks with spiral arms, in time-averaged sense, are effectively in vertical hydrostatic equilibrium, with the main support provided by thermal pressure (in these models without stellar feedback). Unsteady motions associated with shock flapping and movements of dense condensations are mostly horizontal, affecting the vertical force balance relatively little. A small mismatch between H and H_est near x/L_x = −0.1 in S02 models arises from the fact that shocks in these models exhibit weak flapping motions and retain sharp discontinuities in the time-averaged configurations. In this case, the ∂〈ρv_xv_z〉/∂x term in the momentum equation has a non-negligible contribution to the vertical force balance, which was ignored in the derivation of Equation (13).

Zoom In Zoom Out Reset image size

4.3. Velocity Dispersions

Paper I showed that two-dimensional (isothermal) spiral shocks exhibit strong flapping motions in the XZ plane and are able to generate a sonic level of random gas motions in the arm regions. On the other hand, Paper II showed that in one-dimensional spiral shocks with TI, random gas motions amount to only ∼1–2 km s^-1. In this subsection, we quantify the level of random gas motions driven by shock flapping motions and TI in our two-dimensional models.

The velocity field of gas moving across spiral arms is a combination of several different components including streaming motions, oscillations of the shock fronts themselves, and random motions. Since streaming velocities that are ordered and vary with x are much larger than the true random motions of the gas, it is important to subtract the former from the total velocity as cleanly as possible. For this purpose, we first construct time-averaged templates of the velocity field 〈v_i〉 (with i = x, y, or z) for the dense and rarefied components separately. We then calculate the density-weighted velocity dispersions using

$$\begin{equation} \delta v_{i}(t)=\left[\frac{\int \rho [v_{i}-\langle v_{i}\rangle ]^2 dxdz}{\int \rho dx dz}\right]^{1/2}. \end{equation} \tag{ 14 }$$

Figure 12 plots δv_i(t) for the dense and rarefied components in models NU.S10 and NM.S10 as solid and dashed curves, respectively, for a time span of t/t_orb = 5–8. Figure 13 plots the mean values 〈δv_i〉 along with their standard deviations Δδv_i = (〈δv²_i〉 − 〈δv_i〉²)^1/2 for all the non-self-gravitating models; the values of 〈δv_i〉 and Δδv_i are listed in Table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Induced Random Velocity Dispersions

Model	Σ0	Dense Component				Rarefied Component
		〈δvx〉	〈δvy〉	〈δvz〉		〈δvx〉	〈δvy〉	〈δvz〉
	(M☉  pc−2)	(km s−1)	(km s−1)	(km s−1)		(km s−1)	(km s−1)	(km s−1)
NU.S02	2	3.04 ± 1.53	3.28 ± 1.23	1.44 ± 0.34		2.62 ± 0.58	2.25 ± 0.29	1.61 ± 0.17
NU.S05	5	3.12 ± 0.87	3.20 ± 0.92	0.61 ± 0.14		3.04 ± 0.34	2.75 ± 0.33	1.64 ± 0.14
NU.S10	10	4.26 ± 1.75	3.60 ± 1.89	0.41 ± 0.08		3.47 ± 0.81	2.94 ± 0.41	1.72 ± 0.14
NM.S02	2	3.07 ± 1.24	3.24 ± 1.20	1.52 ± 0.43		2.79 ± 0.40	2.45 ± 0.18	1.82 ± 0.11
NM.S05	5	2.99 ± 0.66	3.20 ± 0.81	1.03 ± 0.29		3.17 ± 0.59	2.69 ± 0.46	1.74 ± 0.11
NM.S10	10	3.18 ± 0.62	3.51 ± 1.08	1.11 ± 0.40		3.07 ± 0.42	2.63 ± 0.33	1.69 ± 0.12
SM.S02	2	4.47 ± 1.33	3.97 ± 1.49	0.72 ± 0.18		3.04 ± 0.50	2.49 ± 0.19	1.85 ± 0.14
SM.S05	5	6.36 ± 1.79	6.95 ± 2.11	0.75 ± 0.10		6.07 ± 1.44	5.50 ± 1.45	2.20 ± 0.14
SM.S10	10	10.52 ± 5.55	8.41 ± 3.91	3.60 ± 1.71		10.60 ± 2.94	7.92 ± 1.78	4.96 ± 1.26

Note. Model name prefixes are as in Table 1.

Download table as:  ASCIITypeset image

Figure 12 shows that the density-weighted velocity dispersions for the dense component exhibit large-amplitude fluctuations, with periods roughly of ∼0.5 t_orb. The standard deviations of the fluctuations are Δδv_i ∼ (0.2–0.5)〈δv_i〉 for the dense component; deviations are only Δδv_i ∼ (0.1–0.2)〈δv_i〉 for the rarefied phase. These variations of δv_x and δv_y are caused mostly by oscillations of the shock front relative to the mean position. With large spatial variations of streaming velocities across the arm, the small offset of the shock position as well as the instantaneous locations of dense condensates result in large values of Δδv_i. We thus regard the local minima of δv_i, approximately equal to σ_i ≡ 〈δv_i〉 − Δδv_i, as the upper limits to the level of random gas motions.

Figure 13 shows that for NU models, 〈δv_x〉 and 〈δv_y〉 increase with Σ₀. This is mainly because the shock compression and associated phase transition are stronger with larger Σ₀, leading to stronger flapping motions. Nevertheless, σ_x ∼ σ_y ∼ 2–3 km s^-1 for both dense and rarefied components, insensitive to Σ₀. This indicates that the portion of kinetic energy in the shock flapping motions that goes into random gas motions is quite limited. The remaining portion is simply associated with the horizontal shock oscillations near the midplane. Since the shock flapping motions at low |z| are mostly horizontal, the random vertical motions of the dense gas in NU models are forced predominantly by the impact of rarefied gas arriving from high altitudes. The fact that the rarefied gas has σ_z ∼ 1.7 km s^-1, almost independent of Σ₀, suggests the flapping motions drives more-or-less constant vertical motions at high |z|; this is because the total mass of rarefied gas is almost the same in all models, equivalent to a surface density of 1.9 M_☉  pc⁻². Since the fraction of the rarefied component decreases with increasing Σ₀ (see Table 2), the ratio of vertical kinetic energy in the rarefied gas to the mass of dense gas decreases with Σ₀. This causes σ_z of the dense medium to decrease with increasing Σ₀, roughly as σ_z ∝ Σ^−0.8₀ in NU models. For NM models, the dense gas in the immediate postshock region is overpressured due to the strong heating and thus expands vertically, enhancing σ_z compared to those in NU models.

We now explore the formation of self-gravitating clumps and their internal properties. For NU models, we find that the inclusion of self-gravity always results in catastrophic collapses of self-gravitating clouds that form in the postshock region, preventing us from continuing simulations further. For this reason, we present the results of self-gravitating models only with the density-modified heating rate ("SM" models). Instead of running self-gravitating models from t = 0, we make use of the "saturated-state" data sets of NM models at t/t_orb ∼ 4.8 and restart them by slowly turning on the gaseous self-gravity over a time interval of 1.5 t_orb.

Neglecting the effect of the rarefied medium, the gravitational susceptibility of a midplane dense layer in NM models can be measured by the average Toomre stability parameter

$$\begin{equation} Q_D = \frac{1}{f_D}\frac{\kappa c_D}{\pi G \Sigma _0} = \frac{0.27}{f_D} \left(\frac{\Sigma _0}{10\,{M_\odot \,{\rm pc}^{-2}}}\right)^{-1}, \end{equation} \tag{ 15 }$$

where c_D = 1 km s^-1 is the mean sound speed of the dense gas. With the inclusion of self-gravity, the dense layer in model SM.S10 has Q_D = 0.34 and thus is quite unstable, initiating the collapse of high-density regions. Due to the stiff equation of state, however, the collapsing clouds soon reach an equilibrium state with only a moderate central density 30 cm^-3, which is only 1.5 times larger than the average density of the dense gas in model NM.S10. As these clouds follow galaxy rotation, they merge together in regions of converging streaming velocities, eventually resulting in two big condensations. Figure 14 shows the density snapshots in logarithmic color scale of model SM.S10 at t/t_orb = 7.32, 7.48, and 7.69. Two clouds are widely separated from each other during traversal of the interarm region (Figure 14(a)). After one cloud enters the spiral shock, it loses most of its x-momentum, and the second cloud can then collide with it when it enters the shock at high speed (Figure 14(b)). Since the two clouds are on their own epicyclic orbits before the collision, however, they retain quite different v_y even after the collision, preventing them from merging into a single entity. Due to the Coriolis force, the two clouds subsequently have different v_x, so that the merged entity elongates in the expansion zone (Figure 14(c)), and separates back into two pieces in the interarm region.

Zoom In Zoom Out Reset image size

The dense layer in model SM.S05 has Q_D = 0.9 and is thus marginally unstable. With the aid of shock compression, the dense gas in the postshock region collapses and eventually forms two self-gravitating clumps. In this model, the collision of these clumps at the shock and subsequent break up in the expansion zone is similar to model SM.S10. With Q_D = 12, on the other hand, the dense layer in model SM.S02 is gravitationally stable and does not form dense clumps. Compared with model NM.S02, self-gravity in model SM.S02 increases the fraction of the dense phase by about a factor of 2, which in turn decreases its vertical velocity dispersion by a similar factor.

The presence of self-gravity leads to stronger shock flapping motions than in the NU models, increasing the velocity dispersion. In model SM.S02, the in-plane velocity dispersions of the dense component increases by a factor of ∼1.2–1.5 compared to model NU.S02; for low surface density, the self-gravity is not sufficiently strong to have a major effect. For the SM.S05 and SM.S10 models, however, the stronger self-gravity make a larger difference to shock flapping, in turn driving larger velocity dispersions. Correcting for streaming, the velocity dispersions of both dense and rarefied phases in model SM.S05 reach σ_x ∼ σ_y ∼ 4–5 km s^-1, about twice larger than those in model NM.S05. In model SM.S10, the dense gas velocity dispersions are σ_x ∼ σ_y ∼ 4–5 km s^-1, while the rarefied gas has the velocity dispersions up to ∼7 km s^-1. Note that these in-plane velocity dispersions in multi-phase, self-gravitating models are similar to those in the isothermal self-gravitating models with F = 5% studied in Paper I.

Since the self-gravitating clumps produced in models SM.S05 and SM.S10 move almost ballistically, the position of the largest surface density at a given time does not always correspond to the minimum of the spiral potential. In fact, the clumps are near the potential minimum (|x|/L_x < 0.1) only for Δt/t_orb = 0.35, making the definition of the arm regions rather ambiguous. In addition, due to accretion onto the clumps, the rarefied medium in these models amounts to less than 10% of the total mass, much smaller than the observed mass fraction of the warm gas near the solar neighborhood (Heiles 2001; Heiles & Troland 2003). For these reasons, these clumps are unlike real self-gravitating clouds in spiral galaxies. Nevertheless, we believe these model clouds may provide clues to the internal properties and virial balance of real interstellar clouds, in that they represent a limiting situation in which internal turbulent feedback from star formation is absent.

Keeping in mind the caveats mentioned above, we proceed as follows to calculate the cloud properties. To define the boundary of the dense clumps, we choose the threshold density n_th = 21 cm^-3, corresponding to P_max in the thermal equilibrium curve. Using a CLUMPFIND algorithm (e.g., Williams et al. 1994), we find the interior of each cloud with n>n_th. We then measure the mean density $\rho _{\rm cl}\equiv \overline{\rho }$, the averaged sound speed $c_{\rm cl}\equiv (\overline{P/\rho })^{1/2}$ and the mean one-dimensional internal velocity dispersion $\sigma _{\rm cl}\equiv \sum _i (\overline{v_{i}^2}-\overline{v_{i}}^2)^{1/2} /3^{1/2}$ of each cloud, where the overlines denote the mass-weighted average. We then count the total number of pixels N_cl on the XZ plane occupied by each cloud, and calculate the cloud size R_cl ≡ (N_clΔxΔz/π)^1/2. Assuming a spherical shape, we calculate the total mass M_cl ≡ 4πρ_clR³_cl/3 and the virial ratio of each cloud via

$$\begin{equation} \alpha \equiv \frac{5\big(\sigma _{\rm cl}^2+c_{\rm cl}^2\big)R_{\rm cl}}{GM_{\rm cl}}, \end{equation} \tag{ 16 }$$

which is the ratio of the total kinetic energy to the gravitational potential energy for a uniform spherical cloud (e.g., Bertoldi & McKee 1992; McKee & Ostriker 2007); note that central concentration would decrease α. Figure 15 gives an example of the shapes and internal velocities of two clumps shown in Figure 14(a). Note that the velocity vectors are distributed quite randomly, indicative of (subsonic) turbulent motions within the clouds.

Zoom In Zoom Out Reset image size

Table 4 summarizes the average properties of the clouds that form in models SM.S05 and SM.S10. The typical size and mass of the clouds are found to be R_cl ∼ 45–60  pc and M_cl ∼ (4–9) × 10⁵ M_☉, respectively, with the clouds in model SM.S05 somewhat smaller and less massive than in model SM.S10. Overall, α ∼ 2 for all the clouds, suggesting that they are (marginally) gravitationally bound. The mean sound speed inside the clouds is c_cl = 3.7–5.1 km s^-1, ∼3–4 times larger than the one-dimensional internal velocity dispersions σ_cl = 0.9–1.8 km s^-1. This indicates that the major support against self-gravity comes predominantly from the thermal energy, a consequence of the density-modified heating rate we adopt. The relatively low value of the internal turbulent velocity dispersion suggests that the interaction of a dense, gravitationally-bound cloud with its surroundings can drive only a modest level of internal turbulence.

While stellar spiral arms in disk galaxies provide smoothly varying low-amplitude gravitational potential perturbations, the response of the interstellar gas to them is quite dramatic. Spiral shocks compress the ISM and the high post-shock densities may trigger growth of arm substructures and star formation. In addition, radiative cooling and heating of the gas makes the ISM inherently inhomogeneous, producing two phases that differ in density and temperature by about two orders of magnitude. Moreover, the vertical stellar gravity tends to produce stratification of the cold and warm gas due to differential buoyancy; this stratification can be modified by vertical turbulence, however. Interactions of these processes may significantly affect the gas flows and shock structures, compared with results from our previous work (and that of other groups), which employed an isothermal approximation (Paper I) or neglected the vertical degree of freedom (Paper II). In this paper, we have conducted nonlinear hydrodynamic simulations in a two-dimensional slice perpendicular to a local segment of a spiral arm that is tightly wound (with a pitch angle sin i = 0.1) and rotates rigidly (with a pattern speed Ω_p/Ω₀ = 0.5). To handle the Coriolis force arising from the galaxy rotation self-consistently, we allow for gas motions parallel to the arm (i.e., perpendicular to the domain of the simulation). We consider two different forms of gas heating; the usual constant heating rate (for NU models) and the density-modified heating rate (for NM and SM models), which mimics the effect of star formation feedback in a very simple way, to limit runaway collapse. We start from initially isothermal disks that are in vertical hydrostatic equilibrium but out of thermal equilibrium. We slowly turn on the amplitude of the spiral arm potential such that it attains a full strength at 1.5 orbital times. Magnetic fields are neglected in the present work.

Our main results and their astronomical implications are as follows.

• 1.  
  Two-phase disk equilibria without spiral arms. In the absence of spiral-arm potential perturbations (and other sources of turbulence), the vertical structure of equilibrium disks depends on the disk surface density Σ₀. When Σ₀>Σ_max ≡ P_max/(2Gρ_*c²_R)^1/2 → 6.7 M_☉ pc⁻² (for Solar neighborhood conditions), the disk experiences TI and evolves toward an equilibrium configuration with a thin slab of dense gas (n>1 cm^-3) near the midplane sandwiched between layers of rarefied gas (with n < 1 cm^-3). Here, P_max/k_B = 5000 cm^-3 K is the maximum pressure allowed for the thermally-stable rarefied gas with our adopted heating and cooling functions, c_R ≈ 7 km s^-1 is its density-weighted sound speed, and ρ_* = 0.056 M_☉  pc⁻³ is the stellar density near the Solar neighborhood. In our models, the transition between the dense and rarefied phases occurs approximately at P_trans/k_B ≈ 2100 cm^-3 K, insensitive to Σ₀, as long as the disk is unstable to TI. Without self-gravity, the vertical distribution of the rarefied gas in equilibrium is well fitted by a Gaussian profile whose surface density is fixed to Σ_R = Σ_NG ≡ P_trans/(2Gρ_*c²_R)^1/2 → 2.8 M_☉ pc⁻². When self-gravity is included, the gravity from the midplane dense layer compresses the overlying rarefied component further, forcing overpressured rarefied gas to transform to the dense phase. The resulting surface density of the rarefied gas is reduced to $\Sigma _R ={\mathcal F}\Sigma _{\rm NG}$, with the reduction factor ${\mathcal F}$ defined by Equation (A5). When Σ₀ < Σ_min ≡ P_min/(2Gρ_*c²_R)^1/2 → 2.1 M_☉ pc⁻², the disk has too low a pressure to produce the dense component; equilibrium consists only of the rarefied gas. When Σ_min < Σ₀ < Σ_max, either a single rarefied disk or a two-phase disk is possible, depending on the initial conditions. Our S02 and S05 models with Σ₀ = 2 and 5 M_☉ pc⁻², respectively, that start from a warm isothermal configurations all end up with a single-component rarefied disk, in the absence of spiral perturbations.
• 2.  
  Shock flapping motions in vertically stratified disks. Spiral-arm potential perturbations lead to spiral shocks in the gas, which are vertically curved and non-stationary, showing strong flapping motions perpendicular to the arms. Similarly to the one-dimensional cases studied in Paper II, the shock compression and postshock expansion in two dimensions allow phase transitions, but only if the gas density at the shock and/or the postshock expansion zone reaches the thermally-unstable range (1 cm^-3 ≳ n ≳ 7–9 cm^-3). In model NU.S10 with Σ₀ = 10 M_☉ pc⁻² and the uniform heating rate, the shocked dense gas has large enough density that the postshock expansion is unable to return it to the thermally unstable regime. As a consequence, there is a large amount of interarm dense gas entering the shock in this model, which collides with other dense gas in the arm, producing dense condensations. In other models with lower surface density or the density-modified heating rate, the shocked gas re-expands and becomes thermally unstable, returning to either the dense or the rarefied phase in the interarm region.The shock flapping motions in our models arise due to the fact that the arm crossing time of gas is incommensurate with the vertical oscillation period, so that steady flows are not possible. Seen from the upstream side, the shock is convex when the postshock regions are maximally compressed, and concave when the postshock vertical expansion is strongest. These non-steady motions of shock fronts are commonly seen in numerical simulations with sufficient resolution (e.g., Martos & Cox 1998; Gómez & Cox 2002, 2004; Kim & Ostriker 2006; Paper I; Wada 2008), although simulations with low resolution (Tubbs 1980) or particles (e.g., Dobbs et al. 2008) do not capture the flapping motions clearly.It is interesting to note that radio continuum images of the 5 kpc arm (or the Scutum arm) toward the galactic longitude l = 30°–32° in the Milky Way shows a bow shock feature in the warm ionized medium with temperature ∼10⁴ K (Sofue 1985), similarly to a convex shock front seen in Figure 4(b). The radio emission from the bow shock is presumably thermal radiations from ionized gas, with emission measure of ∼7000  pc cm⁻⁶. The curvature of the observed bow shock, as measured by the longitudinal offset Δl of the shock at latitude b relative to the shock front at the midplane, is Δl/b ∼ 0.5 for b = 05. This value is about a half of the maximum curvature of the shock front |x_sp(H_R) − x_sp(0)|/H_R ∼ 0.85 in our S10 models, where x_sp(z) denotes the shock position at height z. This strongly suggests that the bow shock associated with the 5 kpc arm is most likely a cross section of a galactic spiral shock that is undergoing flapping motions. Velocity information is needed to determine whether the 5 kpc arm regions are currently being compressed or expanding in the course of the flapping motions.In this paper, for consistency with our local approximation we have considered tightly wound arms with a very small pitch angle i ∼ 58, and a pattern speed Ω_p/Ω₀ = 0.5. Observed spiral arms of external galaxies are often more loosely wound with i ∼ 10°–30° (e.g., Seigar et al. 2008) and span a wide range of galactocentric radii with differing Ω_p/Ω₀. For fixed F, a larger arm pitch angle would imply a larger streaming velocity v_0x perpendicular to the arm (see Equation (1)). Spiral shocks would then become correspondingly stronger and shifted farther downstream (e.g., Roberts 1969; Shu et al. 1973; Kim & Ostriker 2002), exhibiting larger amplitude flapping motions (Paper I). On the other hand, |v_0x| is increasingly small as Ω_p approaches Ω₀. Consequently, the spiral shock as well as associated flapping motions would become weaker as one approaches corotation, where the gas would simply concentrate near the potential minimum, without involving a shock.
• 3.  
  Time-averaged shock structure. Within a few orbital times after the development of spiral shocks, gas flows reach a quasi-steady state in the sense that the mass fractions of dense and rarefied gas do not change appreciably with time. For models with Σ₀ = 10 M_☉ pc⁻², the quasi-steady mass fraction of the rarefied gas is f_R ∼ 19%, which can be compared to f_R ∼ 30% when the spiral potential is absent. Despite the shock flapping motions, most of the gas is found close to thermal equilibrium, with a small fraction thermally unstable. The density and temperature PDFs are thus bimodal. For model NU.S10, the dense and rarefied peaks are located at (n, T) ∼ (200 cm^-3, 30 K) and ∼(0.2 cm^-3, 7100 K), respectively, and the dense part of the density PDF is described by a lognormal distribution. The time-averaged structure can be well represented by local vertical hydrostatic equilibrium, supported mainly by the thermal pressure rather than gas random motions. This indicates that the vertical hydrostatic balance is a reasonable approximation even in the presence of spiral shocks. The profiles of surface density perpendicular to the arm are more-or-less symmetric with a shock compression factor of ∼7–10, and have broad arm regions whose width correlates with the strength of the shock flapping motions. The fractional widths of the arm, postshock expansion zone, and interarm region are typically 10%, 20%, and 70% of the arm-to-arm distance, where the gas stays for 15%, 30%, and 55% of the arm-to-arm crossing time, respectively. The shock flapping motions in the XZ plane make the arm wider than in one-dimensional spiral shocks where the arm takes up only 1% of the arm-to-arm distance (Paper II).The dense gas produced from TI and shock compression tends to sink toward the midplane to form a thin slab, while high-altitude regions are dominated by warm rarefied gas. The thickness of the dense slab is H_D ∼ 10–40  pc, depending on the total gas content, heating rate, and presence/absence of self-gravity, while the scale height of the rarefied gas is $H_R\sim 130{\,{\rm \,pc}}\approx c_R/\sqrt{4 \pi G \rho _*}$ insensitive to the parameters. For model NM.S10, the thickness ratio of the dense to rarefied components is about 5, which is not much different from the results of Dobbs et al. (2008) who reported that the warm gas extends vertically up four times more than the cold gas. With high density and pressure, the dense slab would transform to molecular clouds if the appropriate chemical reactions for molecule formation were included (e.g., Dobbs & Bonnell 2007; Dobbs & Price 2008). Thin distributions of the cold dense gas are in fact common in numerical simulations of galactic disks with TI where turbulence is driven by magnetorotational instability (Piontek & Ostriker 2007), stellar feedback via H ii regions (Koyama & Ostriker 2009a, 2009b), or supernovae explosions (Korpi et al. 1999; de Avillez & Berry 2001; Joung & Mac Low 2006; Joung et al. 2009). The observed molecular distribution the Milky Way has a scale height of ∼35  pc within the Solar circle, somewhat reduced for the most massive clouds (e.g., Malhotra 1994; Bronfman et al. 2000; Stark & Lee 2005). Galactic H ii regions are also within about 30  pc of the midplane (Lockman 1977). In the inner Galaxy and near the solar circle, the scale height of the cold H i layer is about 1.5 times smaller than the warm H i gas (e.g., Falgarone & Lequeux 1973; Crovisier 1978; see also Ferriére 2001), although the cold and warm phases appear to have a similar scale height in the outer Galaxy (e.g., Dickey et al. 2009).
• 4.  
  Random gas motions driven by shock flapping motions. The flapping motions of spiral shocks stir the gas and supply random kinetic energy. Allowing for incomplete subtraction of streaming motions in the arm region, the induced density-weighted velocity dispersions are σ_x ∼ σ_y ∼ 2–3 km s^-1 for both dense and rarefied components for the non-self-gravitating models, with larger values corresponding to disks with larger Σ₀. Compared with the results of Paper II where the vertical coordinate was suppressed, these values are similar to those for the dense gas in the arms and larger by a factor of ∼2–3 for the interarm rarefied gas. This implies that it is the rarefied gas that is more efficiently stirred by the shock flapping motions. The self-gravitating models have larger velocity dispersions, in the range σ_x ∼ σ_y ∼ 4–5 km s^-1 for the dense and σ_x ∼ σ_y ∼ 4–7 km s^-1 for the rarefied gas, indicating that self-gravity enhances shock flapping and velocity dispersions, especially for rarefied gas. These in-plane velocity dispersions in the current multi-phase models are similar to those in the isothermal models considered in Paper I.The vertical velocity dispersions of the rarefied gas in NU and NM models are σ_z ∼ 1.7 km s^-1, insensitive to Σ₀. In NU models, the vertical motions of the dense gas are excited preferentially by vertical motions of the rarefied gas. Since the mass fraction of the rarefied gas decreases with Σ₀, the vertical velocity dispersions of the dense gas in NU models is a decreasing function of Σ₀, varying roughly as σ_z ∝ Σ^−0.8₀. In NM models, the postshock gas is overpressured due to enhanced heating and thus expands vertically, increasing σ_z compared to NU models.The level of random gas motions in our models are generally smaller than the observed velocity dispersions ∼7–10 km s^-1 for atomic gas in the solar neighborhood (e.g., Heiles & Troland 2003) and for the larger molecular clouds in the Milky Way (e.g., Stark & Brand 1989; Solomon et al. 1987; Heyer et al. 2009). Thus, we conclude that other sources of the interstellar turbulence (e.g., Mac Low & Klessen 2004; Elmegreen & Scalo 2004) must exceed that provided by spiral shocks. One of the dominant mechanisms is of course supernova explosions (e.g., Korpi et al. 1999; de Avillez & Breitschwerdt 2005; Joung & Mac Low 2006; Joung et al. 2009). In outer regions of galaxies where star formation activity is low, the magnetorotational instability (e.g., Sellwood & Balbus 1999; Piontek & Ostriker 2005, 2007) and/or cosmic infall of gas (e.g., Santillán et al. 2007; Klessen & Hennebelle 2010) may play an important role in driving the ISM turbulence. H ii region expansion and radiation pressure are important in injecting energy into the ISM as GMCs are dispersed (e.g., Matzner 2002; Murray et al. 2010). At large scales, self-gravitating instability with galactic rotation and shear can drive turbulence at near-sonic levels (e.g., Kim & Ostriker 2002; Wada & Norman 2002; Kim & Ostriker 2007; Agertz et al. 2009).
• 5.  
  Effect of self-gravity and properties of self-gravitating clouds. When self-gravity is included, dense gas in SM models with Σ₀ ⩾ 5 M_☉ pc⁻² suffers from gravitational instability, eventually forming two large clouds in each model. These are separate in the interarm regions, temporarily merge in the arm, and then break up into two pieces in the postshock expansion zone. These clouds have a radii ∼45–60  pc and mass ∼(4–9) × 10⁵ M_☉ each, and are gravitationally bound with a virial parameter of α ∼ 2. In our present models, we have not attempted to include realistic star formation feedback, but instead increase the heating rate at high density to prevent collapse. As a consequence, the main support against self-gravity comes from thermal pressure. The mean thermal sound speed and internal velocity dispersion of the clouds are c_cl ∼ 3.7–5.1 km s^-1 and σ_cl ∼ 0.9–1.8 km s^-1, respectively. For models with Σ₀ = 2 M_☉ pc⁻², self-gravity is insufficient to form bound clouds. Nevertheless, self-gravity increases the dense gas fraction by a factor of ∼2 compared to the non-self-gravitating counterpart of this model, which in turn decreases the vertical velocity dispersion of the dense gas by a similar factor.Formation of self-gravitating clouds in our two-dimensional models requires the production of the dense gas due to TI, and then additional shock compression. Although our present models do not capture the cloud destruction process, bound clouds created inside spiral arms may be disrupted before they leave the arms if feedback from star formation is sufficiently strong (Shetty & Ostriker 2008; Wada 2008). Nevertheless, the presence of high-density, self-gravitating clouds in the interarm regions opens an interesting possibility that the spiral shocks—where the diffuse gas is strongly compressed—do not necessarily coincide with the regions of highest gas density (in gravitationally bound clouds). For example, Patrikeev et al. (2006) found strongly polarized nonthermal radio emission that may trace magnetic arms, detected preferentially upstream of the CO arms in the inner disk of the Whirlpool galaxy M51 (see also, e.g., Fletcher et al. 2010). We note, however, that the current unmagnetized models are not yet able to provide clues to the relation between gaseous and magnetic arms. It will be interesting to see how TI, spiral shocks, and realistic star formation feedback conspire with magnetized self-gravitating instabilities to create bound clouds and arm substructures (possibly including separate magnetic arms), and to generate turbulence in the gas.

The authors are grateful to B.-C. Koo and G. Park for information on galactic spiral arms, and also to the referee for an insightful report. This work was supported in part by the Korea Research Foundation Grant funded by the Koran Government (MOEHRD) (KRF–2007–313–C00328) and in part by by KICOS through the grant K20702020016-07E0200-01610 provided by MOST. Simulations were performed by using the supercomputing resource of the Korea Institute of Science and Technology Information through the grant KSC-2009-S02-0008. The work of ECO is supported by the U.S. National Science foundation under grant AST-0908185.

We consider thermally bistable two-phase disks in which a cold, dense layer with surface density Σ_D with thickness H_D is surrounded by a warm, rarefied medium. Since the scale height of the dense layer is very small compared to that of the rarefied gas, we approximate the former as razor-thin (H_D ≈ 0). We further assume that the mass fraction of the rarefied gas is small, so that its self-gravity is unimportant. Let ρ_R(z) denote the density distribution of the rarefied gas. In the presence of the external gravity from a stellar disk of uniform density ρ_*, the condition of vertical hydrostatic equilibrium for the warm gas reads

$$\begin{equation} c_R^2\frac{d \ln \rho _R}{d z}= -4\pi G\rho _*z-2\pi G \Sigma _D {\rm sign}(z), \end{equation} \tag{ A1 }$$

where c_R is the isothermal sound speed of the rarefied gas, assumed to be independent of z. Integrating Equation (A1) over z, one obtains

$$\begin{equation} \rho _R(z)=\rho _R(0)\exp \left[-\frac{1}{2h_g^2}\left(z^2+\frac{\Sigma _D}{\rho _*}|z|\right)\right] \end{equation} \tag{ A2 }$$

where h²_g = c²_R/(4πGρ_*) is the Gaussian scale height which the rarefied gas would have in the absence of self-gravity. The surface density, Σ_R, of the rarefied medium is then given by

$$\begin{equation} \Sigma _R=2\int _{H_D\simeq 0}^\infty \rho _R(z)dz = \Sigma _{\rm NG}\mathcal {F}(s_0), \end{equation} \tag{ A3 }$$

where

$$\begin{equation} \Sigma _{\rm NG}\equiv (2\pi)^{1/2}h_g \rho _R(0) =\frac{P_R(0)}{\sqrt{ 2G\rho _*c_R^2}} \end{equation} \tag{ A4 }$$

is the surface density without gas self-gravity,

$$\begin{equation} \mathcal {F}(s_0)\equiv \exp (s_0){\rm erfc}\big(s_0^{1/2}\big) \end{equation} \tag{ A5 }$$

is the reduction factor, and

$$\begin{equation} s_0\equiv \frac{\pi G\Sigma _D^2}{2c_R^2 \rho _*}, \end{equation} \tag{ A6 }$$

measures the strength of gravity due to the dense gaseous slab relative to the external vertical gravity (see, e.g., Kim et al. 2002). Note that the results of Section 3 suggest that when two-phase equilibria are established, the interface between dense and rarefied media has a constant pressure P_trans, so that we may take ρ_R(0) = P_R(0)/c²_R = P_trans/c²_R since the dense medium has a very small scale height (e.g., Piontek & Ostriker 2007), such that

$$\begin{equation} \Sigma _R=\frac{ P_{\rm trans} \mathcal {F}(s_0)}{ \sqrt{2 G \rho _* c_R^2}}. \end{equation} \tag{ A7 }$$

The scale height of the rarefied medium is given by

$$\begin{equation} H_R^2=\frac{\int _{-\infty }^\infty \rho _R(z)z^2dz}{\int _{-\infty }^\infty \rho _R(z)dz}= h_g^2\left[(1+2s_0)-\sqrt{\frac{s_0}{\pi }}\frac{2}{\mathcal {F}(s_0)}\right], \end{equation} \tag{ A8 }$$

where Equation (A2) is used.

The condition of mass conservation requires

$$\begin{equation} \Sigma _0=\Sigma _D + \Sigma _R \end{equation} \tag{ A9 }$$

so that

$$\begin{equation} s_0=\frac{\pi G (\Sigma _0-\Sigma _R)^2 }{2 c_R^2 \rho _*}. \end{equation} \tag{ A10 }$$

For self-gravitating cases, we fix ρ_* = 0.056 M_☉  pc⁻³ and c_R = 7 km s^-1, and solve Equations (A7) and (A10) iteratively to find Σ_R and Σ_D = Σ₀ − Σ_R as functions of Σ₀. The resulting values of f_R = Σ_R/Σ₀ and H_R are plotted in Figure 3 as dashed lines. For non-self-gravitating cases, Σ_R = Σ_NG = P_R(0)/(2Gρ_*c²_R)^1/2 and H_R = h_g corresponding to s₀ = 0, plotted as solid lines in Figure 3.

Note that since Σ_D < Σ₀, s₀ ≪ 1 if πGΣ²₀/(2c²_Rρ_*) ≪ 1; for our models with ρ_* = 0.056 M_☉  pc⁻³ and Σ₀ < 10 M_☉  pc⁻², πGΣ²₀/(2c²_Rρ_*) < 0.2. When s₀ ≪ 1, ${\cal F}(s_0) \approx 1$, so that H_R ≈ h_g → 128  pc and Σ_R → 2.8 M_☉ pc⁻² for ρ_* = 0.056 M_☉  pc⁻³ and P_trans/k_B = 2100 K cm^-3.