THREE-DIMENSIONAL HYDRODYNAMIC SIMULATIONS OF MULTIPHASE GALACTIC DISKS WITH STAR FORMATION FEEDBACK. II. SYNTHETIC H i 21 cm LINE OBSERVATIONS

In this paper, we utilize recent three-dimensional hydrodynamic simulations of multiphase, turbulent galactic disks presented in Paper I to investigate synthetic 21 cm lines. Our numerical models represent a local patch of a galactic disk including galactic differential rotation, external gravity from stars and dark matter, gaseous self-gravity, interstellar cooling and heating, and star formation feedback. We treat star formation feedback to the neutral ISM in two ways. To represent the radiative stage of expanding supernova (SN) remnants, we inject radial momentum, which drives turbulence at a wide range of scale. We also apply (spatially uniform) heating at a rate proportional to the recent star formation rate (SFR), to model far-ultraviolet (FUV) emission. While simplified, this treatment incorporates both thermal and turbulent energy feedback, both of which are required to regulate the ISM properties and SFR properly. The reader is referred to Paper I for complete descriptions of model setup, numerical methods, and detailed evolution. Here, we briefly summarize the dynamical behavior and physical properties of the disk when it has reached a quasi-equilibrium state.

In equilibrium, the disk models contain highly turbulent gas at a wide range of temperatures. CNM clouds tend to settle to the midplane, where they collect to form more massive gravitationally bound clouds (GBCs). Star formation in these GBCs produces feedback that heats the disk and drives expanding SN shells. Shell expansion and heating puff up the disk vertically, reducing the mean density and rate of GBC formation. Fewer GBCs leads to lower SFR and less feedback, until gas settles to the midplane and begins a new cycle of star formation. After an initial transient, the models reach a state characterized by quasi-periodic fluctuations with period³ of $t_{\rm osc}\sim 0.5\sqrt{\pi /G\rho _{\rm sd}}$, where ρ_sd is the midplane volume density of stars and dark matter.

After a saturated state is reached, the thermal and dynamical equilibrium model proposed by Ostriker et al. (2010) and Ostriker & Shetty (2011) describes average disk properties very well (see also Kim et al. 2011; Shetty & Ostriker 2012; Paper I). The balance between cooling and heating sets the mean midplane thermal pressure P_th to a level between the minimum pressure for the CNM, P_min, and the maximum pressure for the WNM, P_max (see Wolfire et al. 2003). Ostriker et al. (2010) assumed that P_th is approximately equal to the geometric mean value $P_{\rm two}\equiv \sqrt{P_{\rm min}P_{\rm max}}$, which is proportional to the heating rate. Photoelectric heating by FUV radiation, believed to be the dominant heating process for the diffuse ISM (Wolfire et al. 1995, 2003), is proportional to the mean FUV radiation field and hence to the SFR surface density, Σ_SFR. We find the approximation P_th ∼ P_two holds within 60% in our numerical simulations.

Turbulent driving and dissipation are also expected to be balanced, resulting in the midplane turbulent pressure P_turb ∼ P_driv ≡ 0.25(p_*/m_*)Σ_SFR, where p_* and m_* are, respectively, the total radial momentum injected to the ISM by a single SN, and total mass in stars per SN (Ostriker & Shetty 2011). Our numerical simulations confirm that P_turb ∼ P_driv. The total (turbulent plus thermal) midplane pressure is then approximately linearly proportional to Σ_SFR. In addition, the total midplane pressure must match the vertical gravitational weight in order to satisfy vertical dynamical equilibrium (e.g., Piontek & Ostriker 2007; Koyama & Ostriker 2009; Kim et al. 2010; Hill et al. 2012). The level of the SFR in equilibrium is self-regulated such that all these constraints are simultaneously satisfied. When the vertical gravitational field is dominated by that of stars (plus dark matter), the equilibrium model predicts $\Sigma _{\rm SFR}\propto P_{\rm th}\propto P_{\rm tot}\propto \Sigma \sqrt{\rho _{\rm sd}}$, where Σ is the gaseous surface density; our simulations show this is satisfied (Ostriker et al. 2010; Kim et al. 2011; Paper I). This situation is expected to hold for atomic-dominated regions in most galactic disks.

Paper I presents results from three-dimensional simulations at varying Σ. Other model parameters are set as ρ_sd∝Σ² and Ω∝Σ. The midplane volume density of stars and dark matter is ρ_sd = 0.05 M_☉ pc⁻³(Σ/10 M_☉ pc⁻²)², and the angular velocity at the center of the local galactic patch is Ω = 28 km s⁻¹ kpc⁻¹(Σ/10 M_☉ pc⁻²). Table 1 summarizes model input parameters and selected, time-averaged physical quantities from the three-dimensional models of Paper I. Column 1 gives the name of each model from Paper I. Column 2 gives the gas surface density Σ. Columns 3–5 give selected properties measured in the simulations averaged over one orbital time after saturation: the SFR surface density Σ_SFR, the midplane thermal pressure, and the scale height of the WNM, respectively. We list the expected column density of the WNM-only LOS, N_WNM, in Column 6 using Equation (21).

In order to construct synthetic 21 cm lines analogous to observations of the Milky Way Galaxy, we assume an observer sitting at (x, y, z) = (0, 0, 0), the center of the simulation domain. The Galactic center is located toward the −x direction and the z = 0 plane represents the midplane of the Galactic disk. We use the QA10 model, a solar neighborhood analogue. This virtual observer conducts mock observations 10⁴ times toward random LOSs of (l, b) distributed uniformly in dl and d(sin b) to ensure uniform sampling per unit solid angle. Here, the galactic longitude l is measured counterclockwise from the −x axis in the z = 0 plane, and the galactic latitude b is measured from the midplane at z = 0. As we briefly described in Section 2.1 (see also Paper I), the model disk suffers quasi-periodic evolutionary cycles with period t_osc ∼ 0.27t_orb. To represent the range of disk conditions, we sample fairly in time with Δt = 0.1t_osc for a period 2t_osc starting from 1t_orb.

For each LOS, H i 21 cm emission/absorption lines are constructed as follows. We first extract the LOS number density n(s), temperature T_k(s), and velocity v(s) as a function of the path length s along a pencil beam at (l, b). We hereafter omit the functional dependence on s for convenience. The local Cartesian coordinates (x, y, z) can be converted to the Galactic coordinates (s, l, b) using

$$\begin{eqnarray*} x&=&{-}s\cos b \cos l \nonumber \\ y&=&{-}s\cos b \sin l\nonumber \\ z&=&s\sin b. \nonumber \end{eqnarray*}$$

Along each ray, we sample at an interval of Δs = 2 pc, calculating the physical quantities via a trilinear interpolation with the eight nearest grid zones. The background rotational velocity relative to the observer v₀ = (0, −Ωx, 0) is also added to the LOS velocity. The path length is increased until the vertical coordinate z reaches the boundary of the simulation, z = ±L_z/2, by making use of the periodic boundary conditions in the horizontal direction. Since our numerical model is local, it is unphysical to perform mock observations toward very low b, for which the path length would extend beyond the disk scale length R_d ∼ 3.75 kpc (Kalberla & Kerp 2009). For this reason, we limit the path length to be smaller than s_max = 3 kpc, resulting in a minimum latitude for LOSs of b_min = 4°.9.

In order to illustrate our virtual H i sky, we calculate the column density and the harmonic mean temperature as N_H ≡ ∫nds and T_{k, avg} ≡ N_H/∫(n/T_k)ds, respectively. Figure 1 displays aitoff projection images of (a) N_H and (b) T_{k, avg} at t/t_orb = 1. Here, for visualization purposes only, we evenly divide the Galactic longitude and latitude with 1° resolution down to |b| = 0, while mock observations used in this paper are randomly and fairly sampled within an unit solid angle as described above. Note again that the path length is artificially limited to s_max = 3 kpc for |b| < 5°. The overall morphology of the column density map resembles the radio survey maps quite well (e.g., Kalberla et al. 2005).

Zoom In Zoom Out Reset image size

The emissivity and absorption coefficient⁴ of the 21 cm line are respectively given by (e.g., Spitzer 1978; Draine 2011)

$$\begin{equation} j(v_{\rm ch})=\frac{3}{16\pi }h\nu _0A_{10}n\phi (v_{\rm ch}) \end{equation} \tag{ 1 }$$

and

$$\begin{equation} \kappa (v_{\rm ch})=\frac{3}{32\pi }\frac{hc^2}{k\nu _0}A_{10}\frac{n}{T_s}\phi (v_{\rm ch}), \end{equation} \tag{ 2 }$$

where ν₀ = 1.4204 GHz is the frequency at the line center, A₁₀ = 2.8843 × 10⁻¹⁵ s⁻¹ is the Einstein A coefficient of the line (Gould 1994), and ϕ(v_ch) is the normalized line profile as a function of velocity that satisfies ∫ϕ(v_ch)(ν₀/c)dv_ch = 1. We assume a Gaussian velocity distribution

$$\begin{equation} \phi (v_{\rm ch})=\frac{1}{\sqrt{\pi }\Delta v}\frac{c}{\nu _0}e^{-[(v_{\rm ch}-v)/\Delta v]^2}, \end{equation} \tag{ 3 }$$

where Δv ≡ (2kT_k/m_H)^1/2 is the Doppler width of the line. In order to calculate the absorption coefficient, we need the spin temperature, T_s, at each point along the LOS.

The spin temperature is defined such that the Boltzmann distribution at T_s gives the actual level populations. There are three principal mechanisms that determine the level populations of the hyperfine state: collisional transitions, direct radiative transitions by 21 cm photons, and indirect radiative transitions involving intermediate levels mainly by scattered Lyα photons (the Wouthuysen-Field (WF) effect; Wouthuysen 1952; Field 1958). Including all three processes and assuming a balance between excitation and de-excitation, we can calculate the actual level population. The spin temperature is then given by a weighted mean of gas kinetic temperature T_k, the brightness temperature of the background 21 cm radiation field T_R = 3.77 K,⁵ and the effective temperature of the Lyα field T_α (Field 1958; see also Liszt 2001):

$$\begin{equation} T_s = \frac{T_R+y_cT_k+y_\alpha T_\alpha }{1+y_c+y_\alpha }, \end{equation} \tag{ 4 }$$

where the normalized transition probabilities via collisions and Lyα photons are respectively given by

$$\begin{equation} y_c \equiv \frac{T_{0}}{T_k}\frac{R_{10}^c}{A_{10}}\quad \hbox{and}\; y_\alpha \equiv \frac{T_{0}}{T_\alpha }\frac{R_{10}^\alpha }{A_{10}}. \end{equation} \tag{ 5 }$$

Here, $R_{10}^c$ and $R_{10}^\alpha$ are the rates of net downward transition by collision and the WF effect at given temperatures T_k and T_α, respectively, and T₀ ≡ hν₀/k = 0.0681 K. If collisions dominate the transition between levels, T_s = T_k. However, for typical conditions in the WNM, T_s departs from T_k since the rate of collisional transitions is insufficient for "thermalization" (Liszt 2001).

We assume that neutral hydrogen is the only collisional partner and adopt an approximate expression for the collisional de-excitation rate coefficient (Draine 2011, chapter 17; see also Allison & Dalgarno 1969; Zygelman 2005; Furlanetto et al. 2006)

$$\begin{equation} k_{10}\approx \left\lbrace \begin{array}{@{} ll}1.19 \times 10^{-10}T_2^{0.74-0.20\ln T_2} \;{\rm cm}^{-3}{\;\rm s}^{-1} \\ \qquad 20\;{\rm K}<T_k<300\;{\rm K} \\ 2.24\times 10^{-10}T_2^{0.207}e^{-0.876/T_2} \;{\rm cm}^{-3}{\;\rm s}^{-1} \\ \qquad 300\;{\rm K}<T_k<10^3\;{\rm K}\end{array} \right., \end{equation} \tag{ 6 }$$

where T₂ = T_k/100 K. For T_k > 10³ K, we extrapolate the second expression for k₁₀. Collisional de-excitation due to electron impact is less than that from neutral hydrogens provided the ionization fraction remains less than 3%, which is generally the case for the CNM and WNM in the solar neighborhood (Wolfire et al. 1995). The net downward transition rate is $R_{10}^c=nk_{10}$.

We consider the WF effect using the simple parameterized formula for y_α derived by Field (1958):

$$\begin{equation} y_\alpha =5.9\times 10^{11}\frac{n_{\alpha }}{T_\alpha T_k^{1/2}}, \end{equation} \tag{ 7 }$$

where n_α is the Lyα photon number density near the Lyα line center. T_α is determined by the spectrum near Lyα, which involves a difficult scattering problem. To a first approximation, Wouthuysen (1952) and Field (1958) argued that the spectrum near Lyα is given by the Planck curve corresponding to the atomic kinetic temperature so that T_α = T_k (see also Field 1959).⁶ We adopt a fixed value of n_α = 10⁻⁶ cm⁻³ for galactic Lyα radiation (see Liszt 2001).

Since we simply adopt a single, somewhat large value for n_α, the current treatment gives an upper limit of the WF effect (although the WF effect is still insufficient to fully thermalize the WNM—see Liszt 2001). For a lower limit on the WNM spin temperature, we use the spin temperature neglecting the WF effect. Note that the WF effect has little impact for the overall results in this paper in spite of about an order of magnitude difference in the spin temperature of the WNM. This is because we consider channel- or LOS-integrated properties, which are dominated not by the low-optical depth WNM but by the high- or intermediate-optical depth gas. In the remainder of this paper, we mainly present results based calculations of T_s without the WF effect (where we have found it is unimportant), while we compare results with and without the WF effect when it produces any non-negligible differences.

Figure 2 displays the distribution of the spin temperature T_s without the WF effect (top row; (a) and (b)) and with the WF effect (bottom row; (c) and (d)) as a function of the kinetic temperature T_k (left column; (a) and (c)) and the thermal pressure P/k = 1.1nT_k (right column; (b) and (d)) for the snapshot at t/t_orb = 1 of the QA10 model. The gray scale represents the fraction of the total mass in each bin. The red dotted line denotes the value of the spin temperature calculated along the thermal equilibrium density curve, n_eq ≡ Γ/Λ(T_k), for the average heating rate of the QA10 model Γ = 1.6 × 10⁻²⁶ erg s⁻¹ and using Λ given by Equation (6) of Paper I.

Zoom In Zoom Out Reset image size

Figure 2 shows that T_s ∼ T_k for the gas at T_k < 500 K, whereas a range of T_s is possible for the warm and unstable gas with T_k > 1000 K. Note, however, that low values of T_s in warm gas (without the WF effect) occur only when the pressure is quite low, far from the midplane. The WF effect brings T_s closer to T_k at low density and pressure (panels (c) and (d)), but the transition is still not thermalized. Since most of the gas is close to thermal equilibrium (because the cooling time is shorter than the dynamical time), the most populated part of the distribution follows the equilibrium curve. However, a non-negligible fraction of gas is out-of-equilibrium with higher/lower density than n_eq due to dynamical contraction/expansion, resulting in higher/lower spin temperature than the red dotted line (see also Figures 8 and 10 of Paper I). In practice, without the WF effect, the density range for out-of-equilibrium gas is quite wide (more than an order of magnitude) so that T_s spans from a few hundreds of Kelvin to T_k when T_k > 1000 K. With the WF effect, T_s of the WNM is generally larger than 1000 K but still significantly smaller than T_k.

Given T_s, we now can calculate the absorption coefficient from Equation (2), and the contribution to the local optical depth from the ith element along an LOS is then τⁱ(v_ch) = κ(v_ch; s_i)Δs, where s_i = iΔs. The total optical depth along an LOS for a given velocity channel v_ch is then:

$$\begin{equation} \tau (v_{\rm ch})=\sum _{i=1}^N \tau ^i(v_{\rm ch}), \end{equation} \tag{ 8 }$$

where N is the total number of LOS elements in a given (l, b) direction. The synthetic emission line is constructed by summing up the brightness temperature of each element with foreground attenuation:

$$\begin{equation} T_B(v_{\rm ch})=\sum _{i=1}^N\left[{ T_s^i (1-e^{-\tau ^i(v_{\rm ch})})\exp \left({-\sum _{j=1}^{i-1}\tau ^j(v_{\rm ch})}\right)}\right]. \end{equation} \tag{ 9 }$$

We use velocity channels with resolution Δv_ch = 1 km s⁻¹ from −50 km s⁻¹ to 50 km s⁻¹.

Figures 3 and 4 plot examples of high and low optical depth LOSs toward (l, b) = (121°, 12°) and (l, b) = (66°, 33°) at t/t_orb = 1, respectively. In the left column, we plot (top to bottom) the number density, temperature, and velocity versus path length along the LOS. In the right column, we plot the synthetic emission (T_B(v_ch)) and absorption ($1-e^{-\tau (v_{\rm ch})}$) lines as well as the "observed" spin temperature T_{s, obs}(v_ch) (see Equation (12)) as a function of velocity channel. The spin temperature T_s without the WF effect at each point along the LOS, and the harmonic mean spin temperature in each velocity channel, T_{s, avg}(v_ch) (see Equation (11)), are presented as red dotted lines in middle left and bottom right plots, respectively, while the blue dashed lines denote profiles with the WF effect. Note that the optical depth is slightly decreased due to the WF effect (dashed line in middle right plot), while the synthetic emission remains unchanged. We only plot T_{s, avg}(v_ch) for τ(v_ch) > 10⁻³ to show the range of synthetic lines we use for further analysis. From Figure 3, it is evident that density structure is quite complicated along a given LOS, even when T_B(v_ch) and τ(v_ch) are relatively simple.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The column density and harmonic mean spin temperature are the main physical properties that can be inferred directly from the H i 21 cm emission and absorption lines observations (e.g., Heiles & Troland 2003b; Kalberla et al. 2005; Dickey et al. 2009). From Equations (2) and (8), integrating the density along the LOS gives

$$\begin{eqnarray} N_{\rm H} &=& \int n ds = 1.813\times 10^{18}\;{\rm cm}^{-2} \nonumber\\ &&\times \int T_{s, {\rm avg}}(v_{\rm ch})\tau (v_{\rm ch}) dv_{\rm ch} (dv_{\rm ch}\hbox{ in} \; {\rm km}\;{\rm s}^{-1}), \end{eqnarray} \tag{ 10 }$$

where the harmonic mean spin temperature for a given velocity channel is defined by

$$\begin{eqnarray} T_{s, {\rm avg}}(v_{\rm ch}) &\equiv& \frac{\int \kappa (v_{\rm ch})T_s ds}{\tau (v_{\rm ch})}= \frac{\int \kappa (v_{\rm ch}) T_s ds}{\int \kappa (v_{\rm ch}) ds} \nonumber\\ &=& \frac{\int n\phi (v_{\rm ch}) ds}{\int (n/T_s)\phi (v_{\rm ch}) ds}. \end{eqnarray} \tag{ 11 }$$

Note, however, that T_{s, avg}(v_ch) is not a direct observable since the definition in Equation (11) requires full LOS information. For observations, an "observed" spin temperature in a given velocity channel is defined by

$$\begin{equation} T_{s, {\rm obs}}(v_{\rm ch})\equiv \frac{T_B(v_{\rm ch})}{1-e^{-\tau (v_{\rm ch})}}, \end{equation} \tag{ 12 }$$

where T_B(v_ch) is the observed brightness temperature, and $1-e^{-\tau (v_{\rm ch})}$ is obtained by differencing off and on a background source. By inspecting Equation (9), the approximation T_{s, obs}(v_ch) ≈ T_{s, avg}(v_ch) holds for channels with low optical depth and absorption dominated by a single layer (Spitzer 1978). By comparing T_{s, avg}(v_ch) (red and blue) and T_{s, obs}(v_ch) (black solid and dashed) in bottom-right panels of Figures 3 and 4, we can see that the "observed" spin temperatures inferred from the lines in fact agree quite well with the corresponding harmonic mean spin temperatures even near the line center with large optical depth. This is true whether or not WF effect contributes to setting the spin temperature.

In order to quantify the agreement between the "observed" and "true" harmonic-mean spin temperatures, we calculate T_{s, obs}(v_ch)/T_{s, avg}(v_ch) as a function of the channel optical depth τ(v_ch) for all LOSs and channels of a snapshot at t/t_orb = 1. Figure 5 displays the distribution for the case without WF in gray scale. As expected, if the channel is optically thin τ(v_ch) ≲ 0.1, the two spin temperatures are in very good agreement. The difference emerges at τ(v_ch) ≳ 0.1 and gets larger as the channel optical depth increases. T_{s, obs}(v_ch) can either under- or overestimate the "true" harmonic mean spin temperature, but remains within a factor of 1.2 and 1.5 of T_{s, avg}(v_ch) for the channel optical depth of τ(v_ch) ∼ 1 and 10, respectively. Note that T_s is larger when the WF effect is included, which reduces both the contribution to the harmonic-mean average spin temperature and the contribution to the optical depth, from warm regions in each velocity channel. The agreement between T_{s, obs}(v_ch) and T_{s, avg}(v_ch) is best at low optical depth, so the use of lower limit for the WNM spin temperature (omitting WF) provides a conservative conclusion for the agreement between "observed" and "true" values. We have also directly checked that when the WF effect is included, Figure 5 is essentially unchanged.

Zoom In Zoom Out Reset image size

We now explore using T_{s, obs}(v_ch) as a proxy for T_{s, avg}(v_ch) to obtain the "observed" column density from the synthetic lines. Substituting Equation (12) for T_{s, avg}(v_ch) in Equation (10),

$$\begin{eqnarray} N_{\rm H,obs} & \equiv & 1.813\times 10^{18}\;{\rm cm}^{-2} \nonumber\\ &&\times \int \frac{\tau (v_{\rm ch})T_B(v_{\rm ch})}{1-e^{-\tau (v_{\rm ch})}}d v_{\rm ch} (dv_{\rm ch}\hbox{ in} \; {\rm km}\;{\rm s}^{-1}). \end{eqnarray} \tag{ 13 }$$

The quantities on the right-hand side are all observables; Equation (13) is also known as the "isothermal" estimator of the H i column density (Spitzer 1978; Dickey & Benson 1982; Chengalur et al. 2013). In the optically thin limit, we have the "thin" column density

$$\begin{eqnarray} &&N_{\rm H,thin}\equiv 1.813\times 10^{18}\;{\rm cm}^{-2}\int T_B(v_{\rm ch}) d v_{\rm ch} (dv_{\rm ch}\hbox{ in} \; {\rm km}\;{\rm s}^{-1}), \nonumber\\ \end{eqnarray} \tag{ 14 }$$

which is what observers obtain when there is no absorption line information. The "observed" spin temperature obtained from an optical-depth weighted average in a given LOS is:

$$\begin{equation} T_{s, {\rm obs}}\equiv \frac{\int \tau (v_{\rm ch})T_{s, {\rm obs}}(v_{\rm ch})dv_{\rm ch}}{\int \tau (v_{\rm ch})d v_{\rm ch}}. \end{equation} \tag{ 15 }$$

Figure 6 displays distributions of the mock observational data for the ratios of "observed" to "true" (a) column density N_{H, obs}/N_H and (b) harmonic mean spin temperature T_{s, obs}/T_{s, avg} as a function of the integrated optical depth τ_int ≡ ∫τ(v_ch)dv_ch. Figure 6 includes all LOSs and all temporal snapshots over an interval 2t_osc (rather than a single data dump) from the simulation. Here, the "true" harmonic mean spin temperature is defined by T_{s, avg} ≡ ∫nds/∫(n/T_s)ds, which can also be obtained by taking the optical depth weighted average of T_{s, avg}(v_ch) analogous to Equation (15). Since T_{s, obs}(v_ch) reproduces T_{s, avg}(v_ch) very well, N_{H, obs} and T_{s, obs} also give excellent estimates of the true column density and harmonic mean spin temperature, respectively. In particular, 90% and 99% of all LOSs respectively have N_{H, obs} and T_{s, obs} within 5% and 12% of the true values. As in Figure 6, the estimators are best at low τ_int. The bulk of underestimated data points shown near τ_int ∼ 0.5 are due to a specific snapshot at t/t_orb = 1.3, when a cold cloud happens to surround the virtual observer. The cold cloud surrounding the observer provides foreground absorption at moderate opacity for all LOSs, violating the assumption of a single cold layer and resulting in overall underestimation of T_{s, obs}(v_ch) and hence N_{H, obs} and T_{s, obs}.

Zoom In Zoom Out Reset image size

The distribution of the ratio of the "thin" (Equation (14)) to "true" column density is shown as thin contours in Figure 6(a). Evidently, N_{H, thin} significantly underestimates the "true" column density when τ_int > 1. However, the majority of the N_{H, thin}/N_H distribution remains >0.7 up to τ_int < 10, while the distribution extends to values as small as ∼0.2 at high τ_int. In observations, N_{H, thin}/N_{H, obs} ∼ 0.6–0.8 for N_{H, obs} ∼ 4–20 × 10²¹ cm⁻² (Dickey et al. 2003; see also Dickey & Benson 1982).

Another important physical property that can be deduced from emission/absorption line observations is the CNM mass fraction (e.g., Dickey et al. 2009). In the classical two phase picture of the ISM (e.g., Field et al. 1969; Wolfire et al. 1995), the harmonic mean kinetic temperature T_{k, avg} along the LOS is defined by

$$\begin{equation} \frac{N_{\rm H}}{T_{k, {\rm avg}}}\equiv \frac{N_c}{T_c}+\frac{N_w}{T_w}, \end{equation} \tag{ 16 }$$

where N_c and N_w are the column density of the CNM and WNM, respectively, and T_c and T_w are the temperature of the CNM and WNM, respectively. Since T_c ≪ T_w, the CNM mass fraction f_c ≡ N_c/N_H ≈ T_c/T_{k, avg} unless f_c is extremely small.

Similarly, for a two phase ISM, Equation (15) can be simplified (using Equation (10) and T_{s, obs}(v_ch) ≈ T_{s, avg}(v_ch)) as

$$\begin{equation} \frac{N_{\rm H,obs}}{T_{s, {\rm obs}}}\approx \frac{N_c}{T_{s,c}}+\frac{N_w}{T_{s,w}} \approx \frac{N_c}{T_c}+\frac{N_w}{T_{s,w}}, \end{equation} \tag{ 17 }$$

where T_{s, c} and T_{s, w} are the spin temperatures of the CNM and WNM, respectively. The second approximation in Equation (17) utilizes T_{s, c} ≈ T_c from Figure 2. This gives

$$\begin{equation} f_{c, {\rm obs}}\equiv \frac{N_c}{N_{\rm H,obs}}=\frac{T_c}{T_{s, {\rm obs}}}\left({\frac{T_{s,w}-T_{s, {\rm obs}}}{T_{s,w}-T_c}}\right). \end{equation} \tag{ 18 }$$

Here, we have used N_{H, obs} ≈ N_H and we keep the term in parentheses since typical spin temperatures of the WNM are not as high as the kinetic temperature of the WNM (see Figure 2). Thus, for a given T_{s, obs}, one can estimate the CNM mass fraction by assuming a value of T_c and T_{s, w}. If T_{s, w} ≫ T_{s, obs}, then f_{c, obs} ∼ T_c/T_{s, obs}.

In order to test the feasibility of this method, we first calculate the "true" mass fractions of the CNM (f_c; T_k < 184 K), unstable neutral medium (UNM; f_u; 184 K < T_k < 5050 K), and WNM (f_w; T_k > 5050 K) for a given LOS by integrating number density n of each component directly. Figure 7 plots the distribution of the "true" mass fractions of (a) CNM (b) UNM and (c) WNM as a function of the "observed" spin temperature T_{s, obs} without the WF effect for all LOSs and snapshots. For T_{s, obs} < 400 K, the main distribution of f_c follows roughly f_{c, obs} ≈ 80 K/T_{s, obs}, the "observed" CNM mass fraction. Here, we adopt T_c = 80 K based on the mean temperature of the CNM in the simulation. The scatter of the distribution indicates that the temperature of the CNM is not a constant but spans a range of 50 K < T_c < 100 K (see the dotted lines in Figure 7(a)). For f_c < 0.2 with T_{s, obs} > 400 K, however, this simple approximation for f_c employing T_c alone is no longer valid. Instead, by adopting T_{s, w} = 1500 K (the typical spin temperature of the WNM without the WF effect—see Figures 2(a) and (b)), we find that f_{c, obs} from Equation (18) follows the overall trend of the f_c distribution quite well (see the dashed line in Figure 7(a)). In any case, the intrinsic scatter of this mass fraction estimator is as high as a factor of two.

Zoom In Zoom Out Reset image size

Figures 7(b) and (c) show that there is no regime in which either the UNM or the WNM is individually predominant in any range of the spin temperature. For the range of 200 K < T_{s, obs} < 1000 K, the UNM is in the majority, but the WNM mass fraction is nearly comparable. For the range of T_{s, obs} > 1000 K, the WNM becomes the majority, but the UNM fraction is still not negligible. The most probable UNM mass fraction is 30–50% for T_{s, obs} > 1000 K. This implies that any attempt to use the spin temperature to estimate relative UNM and WNM abundances would not be reliable. Observation of spin temperature in the regime T_{s, obs} ≳ 1000 K only implies that there is negligible CNM.

Figure 8 shows the same distributions of phases as a function of the "observed" spin temperature T_{s, obs}, but now with the WF effect included. The distributions of all gas components are unchanged at T_{s, obs} < 1000 K, while f_u and f_w distributions are stretched toward higher T_{s, obs} at T_{s, obs} > 1000 K. When the WF effect is included, T_s spans a narrower range at T_{s, obs} > 1000 K (see Figure 2), and as a result f_w → 1 sharply at around T_{s, obs} ∼ 4000 K. Note however that the value of T_s above 1000 K is quite uncertain, as it depends on poorly constrained parameters controlling the WF effect; the specific values of f_w shown here in that range should therefore be considered cautiously. Taken together, Figures 7 (T_s without WF) and 8 (T_s with WF) imply that the typical "observed" spin temperature of the WNM (including substantial amount of the UNM) would be in a range between 1000 and 5000 K, consistent with previous calculation by Liszt (2001).

Zoom In Zoom Out Reset image size

Classically, H i emission/absorption line observations have reported a negative correlation between the peak optical depth τ_peak and the spin temperature at the peak T_{s, peak} since Lazareff (1975) first noted this correlation. The typical slope between log T_{s, peak} and log τ_peak is found to be −0.35 (e.g., Payne et al. 1983; Kulkarni & Heiles 1987). As noted by Heiles & Troland (2003b; see also Braun & Walterbos 1992), however, not only T_{s, peak} and τ_peak but also N_H and Δv_FWHM are mutually related as

$$\begin{equation} N_{\rm H}=1.93\times 10^{20}\;{\rm cm}^{-2}\;\tau _{\rm peak}\left({\frac{T_{s, {\rm peak}}}{100\;{\rm K}}}\right) \left({\frac{\Delta v_{\rm FWHM}}{\; {\rm km}\;{\rm s}^{-1}}}\right), \end{equation} \tag{ 19 }$$

which can be obtained by integrating Equation (2) over the LOS. Here $\Delta v_{\rm FWHM}=2\sqrt{\ln (2)}\Delta v$ is the FWHM of the line profile. The multivariate analysis performed by Heiles & Troland (2003b) concludes that neither the historical τ_peak–T_{s, peak} relationship nor the τ_peak–N_H–T_{s, peak} relationship has physical significance beyond the relation in Equation (19).

It is more informative to consider the distribution in three-dimensional space of T_B(v_ch)–τ(v_ch)–T_{s, obs}(v_ch) as in Figure 4 of Roy et al. (2013a). Figures 9 and 10 display distributions of the synthetic line data without and with the WF effect, respectively, in the plane of T_B(v_ch) and τ(v_ch) for all LOSs, channels, and simulation snapshots of model QA10 with dv_ch = 1 km s⁻¹. Different color contours in Figures 9(a) and 10(a) represent different spin temperature ranges: T_{s, obs}(v_ch) < 200 K in red, 200 K < T_{s, obs}(v_ch) < 1000 K in green, and 1000 K < T_{s, obs}(v_ch) in blue. Overlaid are the observational data points extracted from Figure 4 of Roy et al. (2013a). The distribution from our simulation follows the observed data points very well. Roy et al. (2013a) pointed out the lack of data points for low or intermediate optical depth measurements (0.01 < τ < 0.1) in high (T_B > 50 K; regime A) and intermediate (1 K < T_B < 10 K; regime B) brightness temperatures; these regions are also weakly populated by our synthetic line data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To understand the paucity of data in the A regime, let us consider an LOS that consists only of static WNM. From Equation (19), we obtain T_{s, peak}τ_peak = 52 K(N_H/10²⁰ cm⁻²)(Δv_FWHM/km s⁻¹)⁻¹, and N_H > 10²¹ cm⁻² is required for Δv_FWHM = 10 km s⁻¹ to have T_B ∼ T_{s, peak}τ_peak > 50 K. Since the WNM number density is about 0.3 cm⁻³ in the solar neighborhood (Wolfire et al. 1995), the path length would need to exceed 1 kpc to have N_H > 10²¹ cm⁻². Note that this is lower limit since we assume a "static" medium. To have such a large column density for the WNM within a narrow velocity range, the path length would need to be much longer. It is thus highly unlikely to have no CNM along such a long path, which would be traversing a low-|b| part of the ISM. Thus, for solar neighborhood conditions, the brightness temperature in a given channel would be as high as in regime A only if there are CNM clouds along the LOS. With a contribution from the CNM, the optical depth of the channel would be higher than that of the WNM-only LOS, and the spin temperature would be lower, moving data points upward from regime A in Figure 9. Since T_B ∼ T_s for an optically thick channel, which is only possible when the CNM dominates, the maximum brightness temperature would be limited to the maximum temperature of the CNM, T_B ≲ 200 K. The majority of both real and mock data points for high optical depth have brightness temperature close to the typical spin temperature of the CNM, T_B ∼ T_c ∼ 80 K.

The reason for the lack of the observational points in regime B is related to the minimum optical depth of the CNM. The column density of a single CNM cloud is ∼n_cl_c, where n_c and l_c are the number density and size of the CNM cloud. The peak optical depth of the cloud is then

$$\begin{equation} \tau _c=0.17\left({\frac{n_c}{10\;{\rm cm}^{-3}}}\right)\left({\frac{l_c}{\;{\rm pc}}}\right) \left({\frac{100\;{\rm K}}{T_s}}\right)\left({\frac{\; {\rm km}\;{\rm s}^{-1}}{\Delta v_{\rm FWHM}}}\right). \end{equation} \tag{ 20 }$$

The brightness temperature of the CNM is then T_B ∼ τ_cT_s ∼ 17 K for a typical small cloud. The high peak optical depth and brightness temperature of a single CNM cloud implies that observational points in regime B with T_B < 10 K and T_s < 200 K would either require a single CNM layer thinner than 1 pc (e.g., tiny scale atomic structure which is not resolved in our simulations; see Heiles 2007 for a review) or far outer wings from multiple CNM clouds. A few data points are observed in this regime in both real and mock observations, in contrast to the relatively sharp limit excluding regime A.

In Figures 9(b) and 10(b), we display the distributions of the synthetic line data without and with the WF effect, respectively, for different ranges of the harmonic mean "kinetic" temperature. The red, green, and blue contours denote CNM (T_{k, avg}(v_ch) < 184 K), UNM (184 K < T_{k, avg}(v_ch) < 5050 K), and WNM (5050 K < T_{k, avg}(v_ch)), respectively. T_{k, avg}(v_ch) is defined analogously to Equation (11) for T_k. As we have seen in Figures 7 and 8, the majority of the UNM and WNM distributions (the innermost green and blue contours) are also completely mixed in T_B(v_ch)–τ(v_ch) plane. This is in contrast to the sharply separated spin temperature distributions in Figures 9(a) and 10(a). The difference between contours of T_{s, obs}(v_ch) and T_{k, avg}(v_ch) again warns against naive use of the spin temperature as a proxy for the gas phase since there is not a one-to-one correspondence between T_{s, obs}(v_ch) and T_{k, avg}(v_ch) at all temperatures.

Using T_{s, obs}(v_ch), however, it is possible to separate the T_B(v_ch)–τ(v_ch) plane into CNM-dominated, UNM-dominated, and UNM–WNM-mixed regimes. The CNM-dominated regime is obviously defined by T_{s, obs}(v_ch) ≲ 200 K. As uncertainty in the WF parameters affects the spin temperature at low pressure (see Figure 2), the upper limit of the UNM-dominated regime is not well defined. From Figure 9(b), with no WF effect, a conservative definition of the UNM-dominated regime is 200 K ≲ T_{s, obs}(v_ch) ≲ 1000 K with T_B(v_ch) ≳ 10 K. From Figure 10(b), with WF effect at a high level included, the UNM-dominated regime is extended to 200 K ≲ T_{s, obs}(v_ch) ≲ 2000 K. In other words, for τ(v_ch) > 10⁻³ the one-to-one correspondence between T_{s, obs}(v_ch) and T_{k, avg}(v_ch) persists up to T_{s, obs}(v_ch) ≲ 1000 K (without the WF effect) and T_{s, obs}(v_ch) ≲ 2000 K (with the WF effect). Note that green points in the observational data are mostly in the UNM-dominated regime regardless of the WF effect. Therefore, the presence of the UNM is evident in the observed data, although the relative proportions of UNM and WNM in observed data is uncertain since the majority of the UNM is expected to be buried in the UNM–WNM–mixed regime.

Recently, Kanekar et al. (2011) have proposed the existence of an H i column density threshold, N_{H, lim} = 2 × 10²⁰ cm⁻² for CNM to form, based on 21 cm emission and absorption line observations detailed in Roy et al. (2013a). They found that observed LOSs have a median spin temperature of T_{s, avg} ∼ 240 K for N_H > N_{H, lim} and T_{s, avg} > 1000 K for N_H ≲ N_{H, lim}. This implies that the LOSs with N_H > N_{H, lim} consist of both CNM and WNM, while the WNM is predominant for the LOSs with N_H ≲ N_{H, lim}.

Here, we explain the observed CNM threshold behavior in the context of thermal and dynamical equilibrium in the ISM. As summarized in Section 2.1, our model disks are in thermal and dynamical equilibrium in an average sense (Kim et al. 2011, Paper I). Vertical dynamical equilibrium demands a balance between vertical pressure support and the weight of the gas. Also, thermal equilibrium for a two-phase neutral medium implies the midplane thermal pressure P_th lies between P_min and P_max ($P_{\rm th}\sim \sqrt{P_{\rm min}P_{\rm max}}$) is a good first order approximation as in Wolfire et al. (2003) and Ostriker et al. (2010). The equilibrium profile for a WNM-only vertical LOS can be approximated by a Gaussian profile with the midplane density of n_w(0) = P_th/(1.1kT_w) and scale height of H_w = σ_{z, w}/(4πGρ_sd)^1/2, where T_w is the thermal equilibrium temperature at P_th for the WNM, $\sigma _{z,w}=(v_{\rm turb}^2+c_w^2)^{1/2}$ is the total (thermal and turbulent) vertical velocity dispersion of the WNM, and ρ_sd is the midplane density of stars and dark matter (which dominate the potential in the solar neighborhood). The one-sided WNM column density along a vertical LOS is

$$\begin{eqnarray} N_{\rm WNM}&=&\int _0^{\infty } n_w(z)dz=\sqrt{\frac{\pi }{2}}n_w(0)H_w =1.14\frac{P_{\rm th}H_w}{kT_w}\nonumber \\ &=&2.0\times 10^{20}\;{\rm cm}^{-2}\left({\frac{P_{\rm th}/k}{3000\;{\rm cm}^{-3}\,{\rm K}}}\right) \left({\frac{T_w}{7000\;{\rm K}}}\right)^{-1}\nonumber\\ &&\times \left({\frac{\sigma _{z,w}}{10\; {\rm km}\;{\rm s}^{-1}}}\right) \left({\frac{\rho _{\rm sd}}{0.1\; {M}_{\odot }\;{\rm pc^{-3}}}}\right)^{-1/2}. \end{eqnarray} \tag{ 21 }$$

The reference values in the last expression are typical of the solar neighborhood; this gives N_WNM comparable to the N_{H, lim} value reported in Kanekar et al. (2011). In Column 7 of Table 1, we list corresponding N_WNM calculated from the numerical outcomes of P_th and H_w in several numerical models representing a range of galactic conditions. For arbitrary galactic latitude b, the column density of a WNM-only LOS would be N_WNM/sin |b|.

Figure 11 displays the distribution of N_Hsin |b| in gray scale as a function of sin |b|, drawn from our model QA10. The observed data points, plotted as filled circles with color denoting the spin temperature (taken from Roy et al. 2013a) follow the same general distribution as the mock observation data (gray scale). In order to quantify where the WNM dominates, we calculate the fraction of no-CNM LOS (f_c = 0) in each bin. From outside to inside, the contours demark no-CNM LOS fraction of >50% (cyan), >70% (green), and >90% (magenta). The innermost magenta contour (more than 90% of LOSs lack CNM) envelopes the observed data points with high spin temperatures (>10³ K). The horizontal dashed line marks N_WNM from Equation (21). Note that the spin temperatures of the observed data change systematically from small (CNM dominated) to large (UNM+WNM dominated) as N_Hsin |b| becomes smaller. This is consistent with Figures 7 and 8, which show log T_{s, obs} > 3 only for LOSs without CNM.

Zoom In Zoom Out Reset image size

Since the disk is highly turbulent and time-variable, the midplane pressure of the WNM can be somewhat larger and smaller than the mean midplane pressure P_th. Thus, it is possible to have no-CNM LOSs with column density larger than N_WNM/sin |b|. However, most LOSs, especially at small |b|, consist of both CNM and WNM, as seen in the region excluded from the contours in Figure 11. Where both WNM and CNM are present, N_Hsin |b| will exceed N_WNM. We note, however, that the lowest values of N_Hsin |b| from Roy et al. (2013a) appear in the regime that we expect will be WNM-dominated. We therefore suggest that the "cold threshold" column density seen by Kanekar et al. (2011) is not the manifestation of a minimum shielding column, but instead represents the vertical H i column with only WNM that is consistent with both thermal and dynamical equilibrium in the local Milky Way disk.

If external galaxies were observed with sufficiently high resolution and sensitivity, the column density of a WNM-only LOS would similarly vary with local disk conditions. In order to address this, we consider what an H i observer would see for different extragalactic ISM conditions as modeled by the set of simulations listed in Table 1. We vertically integrate the gas density through the entire disk to obtain the column density $N_{\rm H}\equiv \int _{-\infty }^{\infty } n dz$. Note that this column density is twice that defined by Equation (10), which assumes an observer at the midplane. We also calculate the column density of the CNM (N_c) and the UNM plus WNM (N_{u + w}) in the same way. Here, we consider the UNM and WNM together since they are not distinguishable in H i line observations (see Figures 7 and 9).

Figure 12 displays the distribution of the UNM+WNM mass fraction N_{u + w}/N_H as a function of the total column density N_H. The color scale shows the distribution of the logarithmic number fraction for each of the simulations listed in Table 1. The symbols and errorbars respectively plot the mean and standard deviation of the distribution for each N_H bin. The vertical dashed line in each panel denotes the column density 2N_WNM for a WNM-only LOS under equilibrium conditions through the full disk thickness. For N_H < 2N_WNM, most of LOSs have N_{u + w}/N_H higher than 90%. Above this column density, LOSs are dominated by the two-phase mixture. Thus, 2N_WNM represents a lower limit on the total column for cold gas to be present.

Zoom In Zoom Out Reset image size

The overall distribution and the threshold column density move toward higher N_H, as the total gas surface density increases from 2.5 M_☉ pc⁻² to 20 M_☉ pc⁻² (from QA02 to QA20 models). The self-regulated equilibrium model for atomic-dominated regions (Ostriker et al. 2010) predicts $P_{\rm th}\propto \Sigma _{\rm SFR}\propto \Sigma \sqrt{\rho _{\rm sd}}$ and $H_w\propto 1/\sqrt{\rho _{\rm sd}}$ if the ratio of thermal to total pressure, and the vertical velocity dispersion are approximately constant, as verified numerically for QA-series (Kim et al. 2011; Paper I). From Equation (21), we thus expect N_WNM∝P_thH_w∝Σ, explaining why N_WNM increases nearly linearly in the QA series (see Table 1).