COLD AND WARM ATOMIC GAS AROUND THE PERSEUS MOLECULAR CLOUD. II. THE IMPACT OF HIGH OPTICAL DEPTH ON THE HI COLUMN DENSITY DISTRIBUTION AND ITS IMPLICATION FOR THE HI-TO-H₂ TRANSITION

We use the HI emission and absorption observations from Paper I. The observations were performed with the Arecibo telescope⁷ using the L-band wide receiver and were made toward 27 radio continuum sources located behind Perseus. The target sources were selected from the NRAO VLA Sky Survey (Condon et al. 1998) based on flux densities at 1.4 GHz greater than ∼0.8 Jy, and they are distributed over a large area of ∼500 deg² centered on the cloud (Figure 1).⁸ The angular resolution of the Arecibo telescope at 1.4 GHz is 35. For the observations, a special procedure was adopted to make a "17-point pattern," which includes 1 on-source measurement and 16 off-source measurements (HT03a; Stanimirović & Heiles 2005). This procedure was designed to consider HI intensity variations across the sky and instrumental effects involving telescope gains. The data were processed using the reduction software developed by HT03a, and the final products for each source include an HI absorption spectrum (${e}^{-\tau (v)}$), an "expected" HI emission spectrum (T_exp(v); HI profile that we would observe at the source position if the continuum source were not present), and their uncertainty profiles. Among the 27 sources, 4C +32.14 was excluded from further analyses because of its saturated absorption spectrum. With an average integration time of 1 hr, the rms noise level in the optical depth profiles was ∼1 × 10⁻³ per 1 km s⁻¹ velocity channel. Finally, the derived optical depth and "expected" emission spectra were decomposed into separate CNM and WNM components using the technique of HT03a, and physical properties (optical depth, spin temperature, column density, etc.) were computed for the individual components. We refer to Paper I for details on the observations, data reduction, line fitting,⁹ and CNM/WNM properties.

Zoom In Zoom Out Reset image size

In order to evaluate different methods for deriving the correction for high optical depth, we also use the HI emission data from the GALFA-HI survey (Stanimirović et al. 2006; Peek et al. 2011). GALFA-HI uses ALFA, a seven-beam array of receivers at the focal plane of the Arecibo telescope, to map the HI emission in the Galaxy. Each of the seven dual-polarization beams has an effective beam size of 39 × 41.

For Perseus, Lee et al. (2012) produced an HI cube centered at (R.A., decl.) = (03^h29^m52^s, +30°34'1'') in J2000¹⁰ with a size of ∼15° × 9° by combining a number of individual GALFA-HI projects. We use the same data here, but extend the HI cube up to ∼60° × 18° to include all radio continuum sources in Paper I. The HI column density image derived from the extended HI cube is shown in Figure 1 along with our continuum sources (4C +32.14 excluded).

We use the N(HI) and N(H₂) images from Lee et al. (2012). To derive the N(HI) image, we integrated the HI emission from ${v}_{\mathrm{LSR}}$ = −5 to +15 km s⁻¹ in the optically thin approximation. This velocity range was determined based on the maximum correlation between the N(HI) image and the Two Micron All Sky Survey (2MASS) based A_V data from the COMPLETE Survey (Ridge et al. 2006). In order to construct the N(H₂) image of Perseus with a large sky coverage, we used the 60 and 100 μm data from the IRIS survey (Miville-Deschênes & Lagache 2005) and derived the dust optical depth at 100 μm (${\tau }_{100}$) by assuming that dust grains are in thermal equilibrium. For this purpose, the emissivity spectral index of $\beta =2$ was adopted, and the contribution from very small grains (VSGs) to the intensity at 60 μm (${I}_{60}$) was removed by calibrating the derived ${T}_{\mathrm{dust}}$ image with DIRBE-based ${T}_{\mathrm{dust}}$ data from Schlegel et al. (1998). We then converted the ${\tau }_{100}$ image into the A_V image by finding the conversion factor X for ${A}_{V}=X{\tau }_{100}$ that results in the minimum difference between the derived A_V and the COMPLETE A_V. This calibration of ${\tau }_{100}$ to the COMPLETE A_V was motivated by Goodman et al. (2009), who showed that dust extinction at near-infrared wavelengths is the best probe of total gas column density. Finally, we measured a local dust-to-gas ratio (D/G) by examining the A_V–N(HI) relation for diffuse regions and derived the N(H₂) image by

$$\begin{eqnarray}&&N({{\rm{H}}}_{2})=\displaystyle \frac{1}{2}\left[\displaystyle \frac{{A}_{V}}{{\rm{D}}/{\rm{G}}}-N(\mathrm{HI})\right].\end{eqnarray} \tag{ 2 }$$

The derived N(HI) and N(H₂) images are at 43 resolution (corresponding to ∼0.4 pc at the distance of 300 pc), and their median 1σ uncertainties are ∼5.6 × 10¹⁹ cm⁻² and ∼3.6 × 10¹⁹ cm⁻². See Sections 3 and 4 of Lee et al. (2012) for details on the derivation of the N(HI) and N(H₂) images and their uncertainties.

We use CO integrated intensity (${I}_{\mathrm{CO}}$) data from Dame et al. (2001). Dame et al. (2001) produced a composite CO survey of the Galaxy at 84 resolution by combining individual observations of the Galactic plane and local molecular clouds. The observations were conducted with the Harvard-Smithsonian Center for Astrophysics (CfA) telescope. The final cube had a uniform rms noise of 0.25 K per 0.65 km s⁻¹ velocity channel. To estimate ${I}_{\mathrm{CO}}$ for Perseus, Dame et al. (2001) integrated the CO emission from ${v}_{\mathrm{LSR}}$ = −15 to +15 km s⁻¹. See Section 2 of Dame et al. (2001) for details on the observations, data reduction, and analyses.

As the HI emission spectrum in the direction of a radio continuum source is affected by absorption, it is not possible to obtain an emission profile along exactly the same line of sight as probed by the HI absorption observation. For this reason, we instead estimated the "expected" HI emission spectrum ${T}_{\mathrm{exp}}(v)$, which is the profile we would observe if the continuum source turned off, by modeling the "17-point pattern" measurements. In this modeling, spatial derivatives of the HI emission (up to the second order) were carefully taken into account (see Section 2.1 of Paper I for details). Then the HI column density in the optically thin approximation, ${N}_{\mathrm{exp}}$, can be calculated by

$$\begin{eqnarray}&&{N}_{\mathrm{exp}}\;({\mathrm{cm}}^{-2})=1.823\times {10}^{18}\int {T}_{\mathrm{exp}}(v){dv}\;({\rm{K}}\;\mathrm{km}\;{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 3 }$$

We compute ${N}_{\mathrm{exp}}$ over the velocity range where ${T}_{\mathrm{exp}}(v)$ is above its 3σ noise level. This velocity range covers all CNM and WNM components for each source, and the median velocity range for all 26 sources is ${v}_{\mathrm{LSR}}$ = −39 to +20 km s⁻¹. We then estimate the uncertainty in ${N}_{\mathrm{exp}}$ by propagating the ${T}_{\mathrm{exp}}(v)$ error spectrum through Equation (3), finding a median of ∼3.0 × 10¹⁸ cm⁻².

Additionally, we use spectra from the GALFA-HI cube (pixel size = 1') in the direction of our radio continuum sources to derive the "expected" emission profiles. This simpler approach has been employed when HI absorption spectra were obtained without any special strategy such as the "17-point pattern" (e.g., Dickey et al. 2000; McClure-Griffiths et al. 2001; Dickey et al. 2003). We extract HI emission spectra from a 9 × 9 pixel region (roughly the "17-point pattern" grid size; Figure 1 of HT03a) centered on each continuum source. As the HI spectra right around the continuum source are likely affected by absorption of the background emission, we exclude the HI spectra from the central 3 × 3 pixel region (roughly one Arecibo beamwidth across). By averaging the remaining 72 spectra, we then compute an average emission spectrum ${T}_{\mathrm{avg}}(v)$ and estimate the corresponding HI column density ${N}_{\mathrm{cube}}$ by

$$\begin{eqnarray}&&{N}_{\mathrm{cube}}\;({\mathrm{cm}}^{-2})=1.823\times {10}^{18}\int {T}_{\mathrm{avg}}(v){dv}\;({\rm{K}}\;\mathrm{km}\;{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 4 }$$

The uncertainty in ${N}_{\mathrm{cube}}$ is estimated by calculating the standard deviation of the extracted 72 spectra, ${\sigma }_{T\mathrm{avg}}(v)$, and propagating it through Equation (4). The median value for all 26 sources is ∼1.4 × 10¹⁹ cm⁻².

In Figure 2, we compare the HI column densities calculated using the two methods. Most sources probe the HI column density of ∼(5–16) × 10²⁰ cm⁻², and the last five sources (3C 132, 3C 131, 3C 133, 4C +27.14, and 4C +33.10) extend this range by a factor of ∼2. We find an excellent agreement between the two methods up to ∼3 × 10²¹ cm⁻². The ratio of ${N}_{\mathrm{exp}}$ to ${N}_{\mathrm{cube}}$ ranges from ∼0.9 to ∼1.1 with a median of ∼1.0, suggesting that the two HI column density estimates are consistent within 10%. Considering the more careful examination of spatial variations in the HI emission, we continue by using ${N}_{\mathrm{exp}}$ as ${N}_{\mathrm{low}-\tau }$ in the following sections.

Zoom In Zoom Out Reset image size

In Paper I, we performed Gaussian decomposition of the optical depth and "expected" emission spectra and calculated the properties of individual CNM and WNM components while considering self-absorption of both the CNM and the WNM by the CNM. All Gaussian decomposition results (peak brightness temperature, peak optical depth, spin temperature, etc.) for each component are presented in Table 2 of Paper I. These results enable us to derive the true total HI column density along a line of sight by

$$\begin{eqnarray}{N}_{\mathrm{tot}}\;({\mathrm{cm}}^{-2}) & = & {N}_{\mathrm{CNM}}+{N}_{\mathrm{WNM}}\\ & = & 1.823\times {10}^{18}\displaystyle \int \left(\displaystyle \sum _{0}^{N-1}{T}_{s,n}{\tau }_{0,n}{e}^{-{[(v-{v}_{0,n})/\delta {v}_{n}]}^{2}}\right.\\ & & +\left.\displaystyle \sum _{0}^{K-1}{T}_{0,k}{e}^{-{[(v-{v}_{0,k})/\delta {v}_{k}]}^{2}}\right){dv}\;({\rm{K}}\;\mathrm{km}\;{{\rm{s}}}^{-1}),\end{eqnarray} \tag{ 5 }$$

where the components with subscript n refer to the CNM, the components with subscript k refer to the WNM, ${\tau }_{0}$ is the peak optical depth, v₀ is the central velocity, T₀ is the peak brightness temperature, and $\delta v$ is the $1/e$ width of the component. Here ${N}_{\mathrm{tot}}$ is calculated over the velocity range determined as having ${T}_{\mathrm{exp}}(v)$ higher than its 3σ noise. For the uncertainty in ${N}_{\mathrm{tot}}$, we use errors of the fitted parameters provided by Gaussian decomposition to perform a Monte Carlo simulation where 1000 ${N}_{\mathrm{CNM}}$ and ${N}_{\mathrm{WNM}}$ values are computed from normally distributed parameters. The standard deviations of the ${N}_{\mathrm{CNM}}$ and ${N}_{\mathrm{WNM}}$ distributions are then added in quadrature to estimate the uncertainty in ${N}_{\mathrm{tot}}$. This method of deriving ${N}_{\mathrm{tot}}$ was used by HT03b for their HI absorption measurements toward background sources randomly located over the whole Arecibo sky. In this study, we focus on a localized group of background sources in the direction of Perseus.

The (integrated) correction factor, $f={N}_{\mathrm{tot}}{/N}_{\mathrm{low}-\tau }={N}_{\mathrm{tot}}/{N}_{\mathrm{exp}}$, is shown in Figure 3 (left) as a function of ${N}_{\mathrm{exp}}$. Clearly, the correction factor increases with ${N}_{\mathrm{exp}}$. We then present the (1/${\sigma }^{2}$)-weighted mean values as the blue squares and the linear fit to all 26 data points as the green solid line¹¹ :

$$\begin{eqnarray}f & = & {\mathrm{log}}_{10}({N}_{\mathrm{exp}}/{10}^{20})\times a+b\\ & = & {\mathrm{log}}_{10}({N}_{\mathrm{exp}}/{10}^{20})(0.32\pm 0.06)+(0.81\pm 0.05).\end{eqnarray} \tag{ 6 }$$

In general, the correction factor ranges from ∼1.0 at ∼3.9 × 10²⁰ cm⁻² to ∼1.2 at ∼1.3 × 10²¹ cm⁻² (maximum uncorrected HI column density in Perseus). While f and ${N}_{\mathrm{exp}}$ show a good correlation (Spearman's rank correlation coefficient of 0.80), there are two sources with relatively high correction factors at ∼10²¹ cm⁻², 3C 092 and 3C 093.1. Interestingly, the two are located behind the main body of Perseus (Figure 1). Their high f values of ∼1.5–1.6 could result from an increased amount of the cold HI in the molecular cloud relative to the surrounding diffuse ISM. The CNM fraction is indeed ∼0.4 for both sources, which is higher than the median value of ∼0.3 for all 26 sources (Section 4.3 of Paper I). However, this is not the maximum CNM fraction in our measurements (∼0.6). Observing a denser grid of radio continuum sources behind Perseus and repeating the calculations would be an interesting way to test the cold HI hypothesis. Finally, the correction factor is also presented as a function of the integrated optical depth in Figure 3 (right). As expected, there is a clear correlation (Spearman's rank correlation coefficient of 0.94).

Zoom In Zoom Out Reset image size

We note that our results are not sensitive to HI components at ${v}_{\mathrm{LSR}}$ < −20 km s⁻¹, which are likely unassociated with Perseus (Section 4.2 of Paper I): limiting the calculation of ${N}_{\mathrm{exp}}$ and ${N}_{\mathrm{tot}}$ to ${v}_{\mathrm{LSR}}$ > −20 km s⁻¹ or excluding the five sources showing such HI components at large negative velocities (corresponding to the sources with ${\mathrm{log}}_{10}({N}_{\mathrm{exp}}/{10}^{20})$ > 1.2 in Figure 3; 3C 132, 3C 131, 3C 133, 4C +27.14, and 4C +33.10) results in linear fit coefficients that are consistent with what we present here within uncertainties.¹²

By assuming that each velocity channel represents gas at a single temperature, Dickey et al. (2000) showed that the correction factor per velocity channel can be written as

$$\begin{eqnarray}&&{f}_{\mathrm{chan}}(v)=\displaystyle \frac{{C}_{0}{T}_{{\rm{s}}}(v)\tau (v)}{{C}_{0}{T}_{\mathrm{exp}}(v)},\end{eqnarray} \tag{ 7 }$$

where ${C}_{0}=1.823\times {10}^{18}$ cm⁻²/(K km s⁻¹). In addition, ${T}_{\mathrm{exp}}(v)$ was expressed as

$$\begin{eqnarray}&&{T}_{\mathrm{exp}}(v)={T}_{{\rm{s}}}(v)(1-{e}^{-\tau (v)}).\end{eqnarray} \tag{ 8 }$$

This equation assumes the absence of any radio continuum source behind the absorbing HI cloud. As a result, ${f}_{\mathrm{chan}}(v)$ simply becomes

$$\begin{eqnarray}&&{f}_{\mathrm{chan}}(v)=\displaystyle \frac{\tau (v)}{1-{e}^{-\tau (v)}}.\end{eqnarray} \tag{ 9 }$$

However, while the radio continuum source is absent in Equation (8), some diffuse radio continuum emission is always present and should not be ignored. This emission includes the CMB and the Galactic synchrotron emission that varies across the sky and becomes strong toward the Galactic plane. We refer to a combination of these contributions as ${T}_{\mathrm{sky}}$, and Equation (8) then has to be rewritten as

$$\begin{eqnarray}&&{T}_{\mathrm{exp}}^{*}(v)={T}_{{\rm{s}}}(v)(1-{e}^{-\tau (v)})+{T}_{\mathrm{sky}}{e}^{-\tau (v)}.\end{eqnarray} \tag{ 10 }$$

Considering that HI emission spectra are generally baseline subtracted during the reduction process, ${T}_{\mathrm{sky}}$ can be removed from both sides of Equation (10). Then ${T}_{\mathrm{exp}}(v)={T}_{\mathrm{exp}}^{*}(v)-{T}_{\mathrm{sky}}$, the quantity we have been working with so far, can be expressed as

$$\begin{eqnarray}{T}_{\mathrm{exp}}(v) & = & {T}_{\mathrm{exp}}^{*}(v)-{T}_{\mathrm{sky}}\\ & = & ({T}_{{\rm{s}}}(v)-{T}_{\mathrm{sky}})(1-{e}^{-\tau (v)}).\end{eqnarray} \tag{ 11 }$$

As a consequence, the correction factor becomes

$$\begin{eqnarray}&&{f}_{\mathrm{chan}}(v)=\displaystyle \frac{{T}_{{\rm{s}}}(v)}{{T}_{{\rm{s}}}(v)-{T}_{\mathrm{sky}}}\displaystyle \frac{\tau (v)}{1-{e}^{-\tau (v)}},\end{eqnarray} \tag{ 12 }$$

or with direct observables,

$$\begin{eqnarray}&&{f}_{\mathrm{chan}}(v)={T}_{\mathrm{sky}}\displaystyle \frac{\tau (v)}{{T}_{\mathrm{exp}}(v)}+\displaystyle \frac{\tau (v)}{1-{e}^{-\tau (v)}}.\end{eqnarray} \tag{ 13 }$$

In order to estimate the contribution from the Galactic synchrotron emission, we use the Haslam et al. (1982) 408 MHz survey of the Galaxy. The brightness temperature at 408 MHz is converted to 1.4 GHz using the spectral index of −2.7. As the absolute Galactic latitude of our continuum sources is generally higher than 10°, the synchrotron contribution is small, with ${T}_{\mathrm{sky}}$ ranging from 2.78 to 2.80 K (Table 1 of Paper I). Based on the histogram of ${T}_{{\rm{s}}}$ for the individual CNM components (Figure 5(b) of Paper I), we can then provide a rough estimate of ${T}_{{\rm{s}}}(v)/[{T}_{{\rm{s}}}(v)-{T}_{\mathrm{sky}}]$: the expected range is narrow, from ∼1.0 to ∼1.2. Clearly, for molecular clouds located closer to the Galactic plane the contribution from the diffuse radio continuum emission will be more significant.

In Equation (7), ${f}_{\mathrm{chan}}(v)$ essentially represents the correction that needs to be applied to ${T}_{\mathrm{exp}}(v)$ to calculate the true brightness temperature profile. As a result, the true total HI column density can be obtained by

$$\begin{eqnarray}{N}_{\mathrm{tot}}\;({\mathrm{cm}}^{-2}) & = & 1.823\times {10}^{18}\\ & & \times \displaystyle \int {f}_{\mathrm{chan}}(v){T}_{\mathrm{exp}}(v){dv}\;({\rm{K}}\;\mathrm{km}\;{{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 14 }$$

We derive ${f}_{\mathrm{chan}}(v)$ for all 26 sources (Equation (13)) and present the results in Figure 4. About 91% of the ${f}_{\mathrm{chan}}(v)$ values are between 1 and 2, and the fraction of velocity channels with ${f}_{\mathrm{chan}}(v)\gt 2$ is very small. We then calculate ${N}_{\mathrm{tot}}$ using Equation (14) and show the (integrated) correction factor, $f={N}_{\mathrm{tot}}/{N}_{\mathrm{exp}}$, as a function of ${N}_{\mathrm{exp}}$ in Figure 5 (left). Here the integration is done over the velocity range where ${T}_{\mathrm{exp}}(v)$ is higher than its 3σ noise. The uncertainty in f is derived by running a Monte Carlo simulation where the optical depth and "expected" emission error spectra are propagated through Equations (13) and (14) to compute 1000 f values. The standard deviation of the f distribution is used as the final uncertainty in f.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Similar to the Gaussian decomposition method, we find a good correlation between f and ${N}_{\mathrm{exp}}$ (Spearman's rank correlation coefficient of 0.84). The linear fit determined using all 26 data points is

$$\begin{eqnarray}&&f={\mathrm{log}}_{10}({N}_{\mathrm{exp}}/{10}^{20})(0.25\pm 0.02)+(0.87\pm 0.02).\end{eqnarray} \tag{ 15 }$$

Additionally, f is plotted as a function of the integrated optical depth in Figure 5 (right), again showing a clear correlation (Spearman's rank correlation coefficient of 0.97).

Both graphs in Figure 5 are very similar to those in Figure 3 for the Gaussian decomposition method. Specifically, the linear fit coefficients are consistent within uncertainties. This is surprising considering that the two methods are very different. In particular, the isothermal method assigns a single spin temperature to each velocity channel, while the Gaussian decomposition method allows a single velocity channel to have contributions from several HI components with different spin temperatures.

In Dickey et al. (2000), the authors updated the isothermal method by incorporating the two-phase approximation. As input parameters, this method then required the spin temperature of the cold HI and the fraction of the warm HI that is in front of the cold HI, the quantity they referred to as q. Dickey et al. (2000) showed that for their SMC data there is no difference between the one- and two-phase approximations regarding the correction factor if $q\gtrsim 0.5$, while the difference becomes more pronounced when $q\lesssim 0.25$. In the Gaussian decomposition method, this fraction (F_k in Paper I) is important as well, but it is difficult to constrain. Thus, the fitting process was repeated for F_k = 0, 0.5, and 1, and these results were used to estimate the final uncertainty in ${T}_{{\rm{s}}}$ (see Section 3 of Paper I and HT03a for details).

To understand why the two different methods result in similar correction factors, we compare Equations (5) and (14) (Appendix A) and find that they become comparable regardless of F_k when $\tau \ll 1$ and ${T}_{\mathrm{sky}}\ll {T}_{{\rm{s}}}$. In our observations of Perseus, the median peak optical depth for all CNM components is ∼0.2, and only a small number of the components has the peak optical depth higher than 1 (10 out of 107; Section 4.1 of Paper I), satisfying the condition. In addition, we already showed that ${T}_{\mathrm{sky}}$ is small, with ∼2.8 K for Perseus. The difference between the Gaussian decomposition and isothermal methods, however, will be more significant for molecular clouds that have a large amount of the cold, optically thick HI and/or a substantial contribution from the diffuse radio continuum emission. Due to the more self-consistent way to derive ${T}_{\mathrm{exp}}(v)$, we continue by using the Gaussian decomposition results for further analyses.

Finally, we note that Equations (6) and (15) could be biased against very high optical depths as they result in saturated absorption spectra, e.g., 4C +32.14, the source we had to exclude from our analyses due to its highly saturated absorption profile. In addition, the equations are based on our explicit assumption of the linear relation between f and log₁₀(${N}_{\mathrm{exp}}/{10}^{20}$). The HI emission and absorption measurements obtained by HT03a and HT03b along many random lines of sight through the Galaxy show that the linear relation indeed describes the observations well up to log₁₀(${N}_{\mathrm{exp}}/{10}^{20}$) ∼ 1.5, and the fitted coefficients are consistent with Equations (6) and (15) within uncertainties (Appendix B). There is some interesting deviation from the linear relation, however, particularly for a few sources with log₁₀(${N}_{\mathrm{exp}}$/10²⁰) ≳ 1. This deviation could suggest a non-linear relation at high column densities, and its significance needs to be further examined with more HI absorption measurements. In the future, it will also be important to study the dense CNM using alternative tracers such as CI and CII fine-structure lines (e.g., Pineda et al. 2013).

Dickey et al. (2000). The correction factor calculated by Dickey et al. (2000) for the SMC using full line-of-sight information can be rewritten as $f={\mathrm{log}}_{10}({N}_{\mathrm{exp}}/{10}^{20})0.667+0.066$ for ${N}_{\mathrm{exp}}\gt 2.5\times {10}^{21}$ cm⁻². The range of the HI column density we probe for Perseus barely overlaps with that in Dickey et al. (2000), as the HI column density in the low-metallicity SMC is significantly higher compared to the Galaxy. While the difference with Dickey et al. (2000) in f values depends on ${N}_{\mathrm{exp}}$, our correction factor is only ∼4% higher than what Dickey et al. (2000) suggest when extrapolated to their maximum HI column density of 10²² cm⁻². Similarly, Dickey et al. (2003) used HI emission and absorption data from the Southern Galactic Plane Survey (McClure-Griffiths et al. 2005) in combination with the isothermal method and estimated f = ∼1.4–1.6 for sources located at 326° $\lt \;l\;\lt$ 333° and $| b|$ ≲ 1°. This correction factor is comparable with what we find for Perseus.

Heiles & Troland (2003a, 2003b). HT03a and HT03bperformed Gaussian decomposition of 79 HI emission/absorption spectral line pairs and derived the correction factor, which they called ${R}_{\mathrm{raw}}=1/f$. They found that ${R}_{\mathrm{raw}}$ ranges from ∼0.3 to ∼1.0, corresponding to f = ∼1.0–3.0 (Appendix B). In particular, they estimated $f\sim 1.3$ for the Taurus/Perseus region, which is similar to what we find for Perseus.

Braun et al. (2009). Our correction factor is smaller than what Braun et al. (2009) and Braun (2012) claimed for M31, M33, and the LMC based on the modeling of HI emission spectra. They found that the correction exceeds an order of magnitude in many cases and increases the global HI mass by ∼30%. Even when considering the correction factor per velocity channel, we find the maximum ${f}_{\mathrm{chan}}(v)$ ∼ 4 for only one source and ${f}_{\mathrm{chan}}(v)\lesssim 3$ for the rest of our sources.

Chengalur et al. (2013). Our correction factor can also be compared with predictions from Chengalur et al. (2013). This study performed Monte Carlo simulations where observationally motivated input parameters such as the column density, the spin temperatures of the CNM and the WNM, and the fraction of gas in each of the different phases were provided for ISM models. Our correction factor versus integrated optical depth plots, Figures 3 (right) and 5 (right), can be directly compared with Figure 1(A) in Chengalur et al. (2013). We find that the correction factor by Chengalur et al. (2013) is significantly higher than our estimate, although the general trend of increasing correction with the integrated optical depth is similar. For example, Chengalur et al. (2013) expect $f\sim$20 when $\int \tau {dv}$ ∼ 10 km s⁻¹, while we find only $f\lt 1.5$. Similarly, our Figures 3 (left) and 5 (left) can be compared with Figure 2(A) in Chengalur et al. (2013). We find that our estimate is consistent with the correction factor by Chengalur et al. (2013) for the column density less than 10²¹ cm⁻², while Chengalur et al. (2013) overestimate at the high end of our column density range. If we extrapolate our relation to 10²² cm⁻², we expect a ∼10 times lower correction factor than what Chengalur et al. (2013) suggest. The reason for their very high correction factor could be the inclusion of extremely high column densities (10²³–10²⁴ cm⁻²) in their ISM models, although it is not clear why it less affects the isothermal estimate of the HI column density. Another possible cause for the high correction factor is the neglect of galactic dynamics (e.g., differential rotation) in their Monte Carlo simulations. As recently shown by Kim et al. (2014), consideration of large-scale gas motions is critical to model an overlap between CNM clouds in velocity space. Without proper modeling of dynamics, therefore, CNM clouds could be overcrowded along the velocity axis, resulting in high correction factors.

Liszt (2014). Using the HI absorption data compiled by Liszt et al. (2010), Liszt (2014) estimated the correction factor for radio continuum sources located at high Galactic latitudes. While they did not provide detailed information about how exactly the derivation was done, their correction factor was small with f less than 1.2 for $E(B-V)$ ≲ 0.5 mag. This is comparable to our finding.

Fukui et al. (2014, 2015). Finally, Fukui et al. (2015) estimated the correction factor for the Galaxy at $| b|$ > 15° by exploring the relation between ${\tau }_{353}$ and N(HI) at 33' resolution. Their Figure 13 shows that the correction factor distribution ranges from ∼1.0 to ∼3.0 and peaks at ∼1.5. Using the same methodology, Fukui et al. (2014) found that the correction increases the total HI mass of the high-latitude molecular clouds MBM 53, MBM 54, MBM 55, and HLCG 92-35 by a factor of ∼2. In general, the correction factor by Fukui et al. (2014, 2015) appears higher than our estimate for Perseus. While we do not perform a detailed comparison with Fukui et al. (2014, 2015), we test their claim of the optically thick HI as an alternative of the "CO-dark" gas in Section 7.

To investigate the impact of high optical depth on the observed HI saturation in Perseus, we first rederive the N(H₂) image using the corrected N(HI). In essence, we use the same methodology as Lee et al. (2012): the A_V image is derived using the IRIS 60/100 μm and 2MASS A_V data, and a local D/G is adopted to estimate N(H₂). We refer to Section 4 of Lee et al. (2012) for details on the method for deriving N(H₂) and its limitation, and we summarize here main results.

• 1.  
  Contamination from VSGs: to estimate ${T}_{\mathrm{dust}}$ from the ratio of I₆₀ to I₁₀₀, the contribution from stochastically heated VSGs to I₆₀ must be removed. For this purpose, we compare our IRIS-based ${T}_{\mathrm{dust}}$ with the DIRBE-based ${T}_{\mathrm{dust}}$ from Schlegel et al. (1998) and find that the contribution from VSGs to I₆₀ is 78%. This is the same as what Lee et al. (2012) found.
• 2.  
  Zero-point calibration for ${\tau }_{100}$: we refine the zero point of the ${\tau }_{100}$ image by assuming that the dust column density traced by ${\tau }_{100}$ is proportional to N(HI) for atomic-dominated regions. Based on the zero point of the ${\tau }_{100}$–N(HI) relation, we add 1.1 × 10⁻⁴ to the ${\tau }_{100}$ image. This is slightly smaller than what Lee et al. (2012) added, i.e., 1.8 × 10⁻⁴.
• 3.  
  Conversion from ${\tau }_{100}$ to A_V: we convert ${\tau }_{100}$ into A_V by adopting X = 740 for A_V = $X{\tau }_{100}$ that results in the best agreement between our IRIS-based A_V and the 2MASS-based A_V from Ridge et al. (2006). This is slightly higher than X = 720 used by Lee et al. (2012). We compare our rederived A_V with the A_V image from Lee et al. (2012) and find that the ratio of the new A_V to the old A_V ranges from ∼0.93 to ∼1.02.
• 4.  
  Deriving a local D/G and N(H₂): we examine the A_V–N(HI) relation and find that the slope of A_V/N(HI) = 1 × 10⁻²¹ mag cm² is a good measure of D/G for Perseus. This is slightly lower than what Lee et al. (2012) estimated, i.e., 1.1 × 10⁻²¹ mag cm², and makes sense considering that our rederived A_V is essentially the same with the A_V image from Lee et al. (2012), while the corrected N(HI) is slightly higher than the uncorrected N(HI). Finally, we derive N(H₂) using Equation (2) and mask pixels with possible contaminations (point sources, the Taurus molecular cloud, and the "warm dust ring"), following what Lee et al. (2012) did. We show the final N(H₂) image in Figure 9, as well as the normalized histograms of the new (this study) and old (Lee et al. 2012) N(H₂) images in Figure 8 (right).

Zoom In Zoom Out Reset image size

In conclusion, in our rederivation of N(H₂), all parameters are identical to or only slightly different from what Lee et al. (2012) used. As a result, the new N(H₂) is comparable to the old N(H₂), shown as a good agreement between the two histograms in Figure 8 (right). The rederived N(H₂) ranges from −1.5 × 10²⁰ cm⁻² to 5.1 × 10²¹ cm⁻² with a median of 4.6 × 10¹⁹ cm⁻², and ∼83% of the pixels whose signal-to-noise ratio is higher than 1 have the new N(H₂) differing from the old N(H₂) by only 10%. Finally, we note that ∼30% of all finite pixels have negative N(H₂) values of mostly around −(1–5) × 10¹⁹ cm⁻², which are very close to zero considering the median uncertainty in N(H₂) of ∼5 × 10¹⁹ cm⁻² (Section 6.2).

As Lee et al. (2012) did, we perform a series of Monte Carlo simulations to estimate the uncertainty in N(H₂). In these simulations, we assess the errors in N(HI) and A_V and propagate them together.

For the uncertainty in N(HI), we combine two terms in quadrature: the error from using a fixed velocity width (${\sigma }_{\mathrm{HI},1}$), and the error from the correction for high optical depth (${\sigma }_{\mathrm{HI},2}$). To estimate ${\sigma }_{\mathrm{HI},1}$, we produce 1000 N(HI) images using 1000 velocity widths randomly drawn from a Gaussian distribution that peaks at 20 km s⁻¹ with 1σ of 4 km s⁻¹. The standard deviation of the simulated N(HI) is then computed for ${\sigma }_{\mathrm{HI},1}$. This ${\sigma }_{\mathrm{HI},1}$ is what Lee et al. (2012) used as their uncertainty in N(HI). Similarly, to derive ${\sigma }_{\mathrm{HI},2}$, we generate 1000 N(HI) images by applying the correction to the N(HI) image from Lee et al. (2012) using 1000 combinations of a and b in Equation (6). These a and b values are again drawn from Gaussian distributions whose peaks and widths correspond to the fitted a and b values and their uncertainties. We find that the median 1σ of N(HI) in our study is ∼8.1 × 10¹⁹ cm⁻².

For the uncertainty in A_V, we repeat the exercise done by Lee et al. (2012): deriving a number of A_V images by changing input conditions (1σ noises of the IRIS 60/100 μm images, β, zero-point calibration for ${\tau }_{100}$), and estimating the minimum and maximum A_V values for each pixel. In this exercise, we find that the contribution from VSGs to I₆₀ varies from 76% to 88%, while X varies from 655 to 855.

Finally, we propagate the uncertainty in N(HI) and the minimum/maximum A_V values through a Monte Carlo simulation in order to produce 1000 N(H₂) images. The distribution of the simulated N(H₂) is then used to estimate the uncertainty in N(H₂) on a pixel-by-pixel basis. We find that the median 1σ of N(H₂) in our study is ∼4.8 × 10¹⁹ cm⁻².

From the rederived N(HI) and N(H₂) images, we estimate ${{\rm{\Sigma }}}_{\mathrm{HI}}$ and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ by

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\mathrm{HI}}\;({M}_{\odot }\;{\mathrm{pc}}^{-2}) & = & \displaystyle \frac{N(\mathrm{HI})\;({\mathrm{cm}}^{-2})}{1.25\times {10}^{20}}\\ {{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}\;({M}_{\odot }\;{\mathrm{pc}}^{-2}) & = & \displaystyle \frac{N({{\rm{H}}}_{2})\;({\mathrm{cm}}^{-2})}{6.25\times {10}^{19}}.\end{eqnarray} \tag{ 16 }$$

We find that ${{\rm{\Sigma }}}_{\mathrm{HI}}$ varies by only a factor of ∼2.6 from ∼4.8 to ∼12.7 ${M}_{\odot }$ pc⁻². On the other hand, ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ ranges from −2.4 to 81.8 ${M}_{\odot }$ pc⁻², although ∼98% of the pixels have ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ < 15 ${M}_{\odot }$ pc⁻². As a result, ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ has a small dynamic range of ∼8–30 ${M}_{\odot }$ pc⁻² across most of the cloud.

To compare with the KMT09 predictions aiming at revisiting the HI saturation in Perseus, we focus on the individual dark (B5, B1E, and B1) and star-forming (IC 348 and NGC 1333) regions. The boundaries of each region were determined based on the ¹³CO emission (see Section 5 of Lee et al. 2012 for details) and are shown as the red rectangular boxes in Figure 9. Using all finite pixels in the rectangular boxes, we plot ${{\rm{\Sigma }}}_{\mathrm{HI}}$ and ${R}_{{{\rm{H}}}_{2}}$ as a function of ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ for each region and present the results for IC 348 and B1E in Figures 10 and 11. Additionally, we show the same plots from Lee et al. (2012) for comparison. Note that both this study and Lee et al. (2012) include negative ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ values by using all finite pixels in the rectangular boxes. Almost all (∼90%) of these negative ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ values fluctuate around zero within uncertainties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As in Lee et al. (2012), the following KMT09 predictions are used to fit the observed ${R}_{{{\rm{H}}}_{2}}$ versus ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ profiles:

$$\begin{eqnarray}&&{R}_{{{\rm{H}}}_{2}}=\displaystyle \frac{4{\tau }_{{\rm{c}}}}{3\psi }\left[1+\displaystyle \frac{0.8\psi {\phi }_{\mathrm{mol}}}{4{\tau }_{{\rm{c}}}+3({\phi }_{\mathrm{mol}}-1)\psi }\right]-1\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{3}{4}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{comp}}{\sigma }_{{\rm{d}}}}{{\mu }_{{\rm{H}}}}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\psi =\chi \displaystyle \frac{2.5+\chi }{2.5+\chi e},\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&\chi =2.3\;\displaystyle \frac{1+3.1{Z}^{\prime 0.365}}{{\phi }_{\mathrm{CNM}}}.\end{eqnarray} \tag{ 20 }$$

Here ${\tau }_{{\rm{c}}}$ is the dust optical depth a spherical cloud would have if its HI and H₂ are uniformly mixed, and χ is the ratio of the rate at which LW photons are absorbed by dust grains to the rate at which they are absorbed by H₂. In addition, ${{\rm{\Sigma }}}_{\mathrm{comp}}$ is the total gas column density, ${\sigma }_{{\rm{d}}}$ is the dust absorption cross section per hydrogen nucleus in the LW band, ${\mu }_{{\rm{H}}}$ is the mean mass per hydrogen nucleus, Z' is the metallicity normalized to the value in the solar neighborhood, ${\phi }_{\mathrm{CNM}}$ is the ratio of the CNM density to the minimum CNM density at which the CNM can be in pressure balance with the WNM, and finally ${\phi }_{\mathrm{mol}}$ is the ratio of the H₂ density to the CNM density. We refer to Section 6 of Lee et al. (2012) for a detailed summary of the KMT09 model.

As Equations (17)–(20) suggest, ${R}_{{{\rm{H}}}_{2}}$ in the KMT09 model is simply a function of total gas column density, metallicity, ${\phi }_{\mathrm{CNM}}$, and ${\phi }_{\mathrm{mol}}$ and is independent of the strength of the radiation field in which the cloud is embedded. Following Lee et al. (2012), we adopt Z' = 1 and ${\phi }_{\mathrm{mol}}=10$ (fiducial value by KMT09)¹³  and constrain ${\phi }_{\mathrm{CNM}}$ by finding the best-fit model for the observed ${R}_{{{\rm{H}}}_{2}}$ versus ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$. For this purpose, we perform Monte Carlo simulations where the uncertainties in ${R}_{{{\rm{H}}}_{2}}$ and ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ are taken into account for model fitting. In these simulations, we add random offsets to ${R}_{{{\rm{H}}}_{2}}$ and ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ based on their uncertainties and determine the best-fit curve by setting Z' = 1 and ${\phi }_{\mathrm{mol}}=10$ and finding ${\phi }_{\mathrm{CNM}}$ that results in the minimum sum of squared residuals. We repeat this process 1000 times and estimate the best-fit ${\phi }_{\mathrm{CNM}}$ by calculating the median ${\phi }_{\mathrm{CNM}}$ among the simulated ${\phi }_{\mathrm{CNM}}$. The derived ${\phi }_{\mathrm{CNM}}$ for each region is summarized in Table 1, and the best-fit curves are shown in red in Figures 10 and 11.

Table 1.  Fitting Results for ${R}_{{{\rm{H}}}_{2}}$ vs. ${{\rm{\Sigma }}}_{{\rm{H}}\;{\rm{I}}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$

Region	Best-fit ${\phi }_{\mathrm{CNM}}$
B5	8.75 ± 1.35
IC 348	7.40 ± 0.94
B1E	7.28 ± 1.13
B1	6.93 ± 1.00
NGC 1333	5.28 ± 0.85

Note. The uncertainty in ${\phi }_{\mathrm{CNM}}$ is estimated from the distribution of the simulated 1000 ${\phi }_{\mathrm{CNM}}$ values.

Download table as:  ASCIITypeset image

In the KMT09 model, χ measures the relative importance of dust shielding and H₂ self-shielding and is predicted to be ∼1 for a wide range of galactic environments. In this case, a certain amount of ${{\rm{\Sigma }}}_{\mathrm{HI}}$ is required to shield H₂ againt photodissociation, and H₂ forms out of HI once this minimum shielding column density is obtained. The KMT09 model predicts the minimum shielding column density of ∼10 ${M}_{\odot }$ pc⁻² for solar metallicity, and this is indeed consistent with what we observe in Perseus: ${{\rm{\Sigma }}}_{\mathrm{HI}}$ saturates at ∼7–9 ${M}_{\odot }$ pc⁻² for all five regions. The level of the HI saturation changes between the regions, though, from ∼7 ${M}_{\odot }$ pc⁻² for B5 to ∼9 ${M}_{\odot }$ pc⁻² for NGC 1333. In comparison with Lee et al. (2012), our ${{\rm{\Sigma }}}_{\mathrm{HI}}$ values are slightly higher due to the correction for high optical depth. This correction brings a closer agreement with the KMT09 model.

The excellent agreement with the KMT09 model is also evident from ${R}_{{{\rm{H}}}_{2}}$ versus ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ in Figure 11. For all five regions, we find that ${R}_{{{\rm{H}}}_{2}}$ steeply rises at small ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$, turns over at ${R}_{{{\rm{H}}}_{2}}$ ∼ 1, and then slowly increases at large ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$. In fact, this common trend on a log-linear scale is entirely driven by the almost constant ${{\rm{\Sigma }}}_{\mathrm{HI}}$, and therefore it is not surprising to find such a good agreement with the KMT09 model where the ${{\rm{\Sigma }}}_{\mathrm{HI}}$ saturation is predicted. By fitting the KMT09 predictions to the observed ${R}_{{{\rm{H}}}_{2}}$ versus ${{\rm{\Sigma }}}_{\mathrm{HI}}+{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ profiles, we constrain ${\phi }_{\mathrm{CNM}}$ ∼ 5–9, which is consistent with what Lee et al. (2012) estimated within uncertainties. In the KMT09 model, ${\phi }_{\mathrm{CNM}}$ determines the CNM density by ${n}_{\mathrm{CNM}}={\phi }_{\mathrm{CNM}}{n}_{\mathrm{min}}$, where ${n}_{\mathrm{min}}$ is the minimum CNM density at which the CNM can be in pressure balance with the WNM, and ${\phi }_{\mathrm{CNM}}=5-9$ translates into ${T}_{\mathrm{CNM}}=60-80$ K (Equation (19) of KMT09). This ${T}_{\mathrm{CNM}}$ range is consistent with what we observationally constrained for the HI environment of Perseus via the HI absorption measurements (Paper I): ${T}_{\mathrm{CNM}}$ mostly ranges from ∼10 to 200 K, and its distribution peaks at ∼50 K. Finally, we note that ${\phi }_{\mathrm{CNM}}$ systematically decreases toward the southwest of Perseus, reflecting the observed region-to-region variations in ${{\rm{\Sigma }}}_{\mathrm{HI}}$. The difference in ${\phi }_{\mathrm{CNM}}$, however, is not significant, and this suggests similar χ values for all dark and star-forming regions (Equation (20)). We indeed find χ = ∼1.1–1.8, in agreement with the KMT09 prediction of comparable dust shielding and H₂ self-shielding for H₂ formation.