THE CO-TO-H₂ CONVERSION FACTOR ACROSS THE PERSEUS MOLECULAR CLOUD

To understand how X_CO varies with A_V, we begin by plotting I_CO as a function of A_V in Figure 8 for all five regions defined in Section 4. We use the A_V image at 43 angular resolution derived by Lee et al. (2012). Even though there is a large amount of scatter, several important features are noticeable.

Zoom In Zoom Out Reset image size

First, there appears to be some threshold A_{V, th} ∼ 1 mag below which no CO emission is detected. The sharp increase of I_CO with A_V found from the individual regions (Section 6.2.2) strongly supports the existence of such threshold. This may suggest that CO becomes shielded against photodissociation at A_V ∼ 1 mag in Perseus. Previous observations of molecular clouds have found a similar A_{V, th} ∼ 1 mag (e.g., Lombardi et al. 2006; Pineda et al. 2008; Leroy et al. 2009). Note that a lack of CO detection at A_V ≲ 1 mag is not due to our sensitivity, considering that our mean 1σ uncertainty of A_V is ∼0.2 mag. In addition, the threshold is not the result of the limited spatial coverage of the COMPLETE I_CO image. We made a comparison between our A_V image and the I_CO image with a large spatial area of ∼10° × 7° from the Center for Astrophysics CO Survey (CfA; Dame et al. 2001) at the common angular resolution of 84 and found essentially the same threshold.

Second, I_CO significantly increases from ∼0.1 K km s⁻¹ to ∼70 K km s⁻¹ for a narrow range of A_V ∼ 1–3 mag. This steep increase of I_CO may suggest that the transition from C ii/C i to CO is sharp once shielding becomes sufficiently strong to prevent photodissociation (e.g., Taylor et al. 1993; Bell et al. 2006).

Third, I_CO gradually increases and saturates to ∼50–80 K km s⁻¹ at A_{V, sat} ≳ 3 mag. This is consistent with Pineda et al. (2008), who found A_{V, sat} ∼ 4 mag for Perseus. Similarly, Lombardi et al. (2006) found the saturation of I_CO ∼ 30 K km s⁻¹ for the Pipe nebula at A_{V, sat} ∼ 6 mag (with their adopted relation A_V = A_K/0.112). The saturation of I_CO is expected on the basis of the relation between I_CO and τ, I_CO ∝ 1 − e^−τ, where τ ∝ A_V. Therefore, I_CO does not faithfully trace A_V once it becomes optically thick. The presence of optically thick CO emission in Perseus was hinted by Pineda et al. (2008), who performed the curve of growth analysis for the CO and ¹³CO(J = 1 → 0) observations.

In agreement with Pineda et al. (2008), we find that the relation between I_CO and A_V has significant region-to-region variations across Perseus, contributing to the large scatter in Figure 8. We therefore show I_CO versus A_V for each dark and star-forming region in Figure 9. To emphasize the steep increase and saturation of I_CO, we plot I_CO as a function of A_V on a log–log scale.

Zoom In Zoom Out Reset image size

Among the five regions, B5 and L1448 have the narrowest range of A_V ∼ 1–3 mag, simply reflecting their smaller N(H) range on average. On the other hand, IC 348 has the largest ranges of A_V ∼ 1–11 mag and I_CO ∼ 0.2–50 K km s⁻¹. I_CO steeply increases from ∼0.2 K km s⁻¹ to ∼35 K km s⁻¹ at A_V ∼ 1–3 mag and then saturates to ∼50 K km s⁻¹ at A_V ≳ 3 mag. In the case of B1E/B1, two components are apparent. The first component corresponds to the relatively steep increase of I_CO from ∼1 K km s⁻¹ to ∼20 K km s⁻¹ at A_V ∼ 1.5–3 mag. The second component corresponds to the gradual increase of I_CO from ∼20 K km s⁻¹ to ∼60 K km s⁻¹ at A_V ∼ 1.5–5 mag. Considering the two components together, I_CO saturates to ∼60 K km s⁻¹ at A_V ≳ 3 mag. Lastly, NGC 1333 has the majority of the data points (∼90%) at A_V ≲ 3 mag. I_CO increases from ∼0.5 K km s⁻¹ to ∼70 K km s⁻¹ at A_V ∼ 1–3 mag and then shows a hint of the saturation to ∼80 K km s⁻¹ at A_V ≳ 3 mag. Note that NGC 1333 is the region where I_CO saturates to the largest value in Perseus.

In summary, the most important properties we find from the individual regions are the abrupt increase of I_CO at A_V ≲ 3 mag and the saturation of I_CO at A_V ≳ 3 mag. However, I_CO saturates to different values, from ∼50 K km s⁻¹ for IC 348 to ∼80 K km s⁻¹ for NGC 1333.

We use a modified form of the model in Wolfire et al. (2010) to calculate H₂ and CO abundances and CO line emission. The model in Wolfire et al. (2010) uses a plane-parallel PDR code with one-sided illumination to estimate the distributions of atomic and molecular species as a function of A_V into a cloud. The density distribution is taken to be the median density as expected from turbulence, and the distribution is converted into a spherical geometry. In our modified W10 model, a plane-parallel slab of gas is illuminated by UV photons on two sides and has either a uniform density distribution or a distribution described with a simple step function. The gas temperature and the abundances of atomic and molecular species are calculated as a function of A_V under the assumptions of thermal balance and chemical equilibrium. For details on the chemical and thermal processes, we refer the reader to Tielens & Hollenbach (1985), Kaufman et al. (2006), Wolfire et al. (2010), and Hollenbach et al. (2012).

The input parameters for the modified W10 model are n, G, ζ, v_turb, Z, and DGR. Considering the constraints on the physical parameters of Perseus (Section 2.2), we use a set of the modified W10 models with the following inputs: G = 0.5 $G_{0}^{\prime }$, ζ = 10⁻¹⁶ s⁻¹, v_turb = 4 km s⁻¹, Z = 1 Z_☉, and DGR = 1 × 10⁻²¹ mag cm². For the density distribution, we use both a uniform density distribution with n = 10³, 5 × 10³, and 10⁴ cm⁻³ and a "core–halo" density distribution. The "halo" consists of H i with a fixed density n_halo = 40 cm⁻³, comparable to diffuse cold neutral medium (CNM) clouds (Wolfire et al. 2003), and has N(H i) = 4.5 × 10²⁰ cm⁻² on each side of the slab. In the "core," on the other hand, n abruptly increases to a large density n_core = 10³, 5 × 10³, or 10⁴ cm⁻³. This "core–halo" structure is motivated by observations of molecular clouds that have found H i envelopes with N(H i) ∼ 10²¹ cm⁻² (e.g., Imara & Blitz 2011; Imara et al. 2011; Lee et al. 2012). As the minimum density of the densest regions for both the uniform and "core–halo" density distributions (∼10³ cm⁻³) has already been constrained by previous comparisons between CO observations and LVG and PDR models (Section 2.2), we expect that the modified W10 model with a density much smaller than 10³ cm⁻³ would not reproduce the observed I_CO in Perseus and therefore does not demonstrate the effect of n, n_core < 10³ cm⁻³ in this paper. In addition, we note that the modified W10 model is not sensitive to the exact value of n_halo, as long as this is small enough to contain a small amount of H₂ and CO in the halo (Section 7.1.2 for details).

We run the model for A_V = 0.6, 0.8, 1, 1.5, 2, 2.8, 4.8, 7.2, and 10 mag (uniform density) and A_V = 1.25, 1.3, 1.5, 1.7, 2, 2.8, 4.8, 7.2, and 10 mag ("core–halo"), and the output quantities are N(H i), N(H₂), and I_CO for a given A_V. We summarize the ranges of the output quantities in Tables 2 ("core–halo") and 3 (uniform density). Note that for both the uniform and "core–halo" density distributions, an increase in A_V can be thought of as an increase in size of the dense region. For example, A_V = DGR × N(H) = DGR(n_coreL_core + n_haloL_halo) = 3.1 × 10⁻³ L_coren_core + 0.9 mag, with L_core in units of pc and n_core in units of cm⁻³ for the "core–halo" density distribution. The "core" has a typical size of L_core ≲ 1 pc, while the "halo" is significantly more extended with L_halo ∼ 7 pc. For the uniform density distribution, the size of the slab is generally L_uniform ≲ 1 pc. We note that in the most extreme case the size of the dense region is much smaller than our spatial resolution (L_core ∼ 0.01 pc), implying a considerably small filling factor of the "core" relative to the "halo," but is comparable to the size of small-scale clumps observed in the CO emission (e.g., Heithausen et al. 1998; Kramer et al. 1998).

Table 2. Predictions from the Modified W10 Model with a "Core–Halo" Density Distribution^a

	ncore = 103 cm−3	ncore = 5 × 103 cm−3	ncore = 104 cm−3
N(H i)b (cm−2)	9.00 × 1020–9.60 × 1020	9.00 × 1020–9.20 × 1020	9.00 × 1020–9.00 × 1020
N(H2)b (cm−2)	1.75 × 1020–4.52 × 1021	1.75 × 1020–4.54 × 1021	1.75 × 1020–4.55 × 1021
$R_{\rm H_2}$b	0.39–9.42	0.39–9.87	0.39–10.11
ICOb (K km s−1)	0.075–36.90	0.33–41.40	0.39–33.90
Lhaloc (pc)	7.28	7.28	7.28
Lcoreb, d (pc)	0.11–2.94	0.023–0.59	0.011–0.29
Lhalo-coreb, e (pc)	7.39–10.22	7.30–7.87	7.29–7.57
〈n〉b, f (cm−3)	54.74–316.66	55.42–411.21	55.49–427.51

Notes. ^aFor all three models, n_halo = 40 cm⁻³. ^bThe values are provided for the minimum and maximum column densities (A_V = 1.25 mag and 10 mag). ^cThe size of the diffuse halo L_halo = 9 × 10²⁰ cm⁻²/40 cm⁻³. ^dThe size of the dense core L_core = (N(H) − 9 × 10²⁰)/n_core. ^eThe total size of the slab L_halo-core = L_halo + L_core. ^fThe average density 〈n〉 = N(H)/L_halo-core.

Download table as:  ASCIITypeset image

Table 3. Predictions from the Modified W10 Model with a Uniform Density Distribution

	n = 103 cm−3	n = 5 × 103 cm−3	n = 104 cm−3
N(H i)a (cm−2)	5.71 × 1019–1.30 × 1020	1.15 × 1019–2.66 × 1019	5.59 × 1018–1.33 × 1019
N(H2)a (cm−2)	2.72 × 1020–4.94 × 1021	2.95 × 1020–4.99 × 1021	2.98 × 1020–5.00 × 1021
$R_{\rm H_2}$a	9.53–76.00	51.30–375.19	106.62–751.88
ICOa (K km s−1)	0.84–46.90	9.79–62.40	18.00–60.20
Luniforma, b (pc)	0.19–3.24	0.039–0.65	0.019–0.32

Notes. ^aThe values are provided for the minimum and maximum column densities (A_V = 0.6 mag and 10 mag). ^bThe total size of the slab L_uniform = N(H)/n.

Download table as:  ASCIITypeset image

We compare X_CO versus A_V with predictions from the modified W10 model ("core–halo") in Figure 11(a). While B5 and L1448 probe too narrow ranges of A_V for significant comparisons, the model curves with n_core = 10^{3 − 4} cm⁻³ follow the observed trends for IC 348 and B1E/B1. The situation is more complicated for NGC 1333, where the model matches the observed X_CO only for a partial range of A_V and has difficulties in reproducing the observations at A_V ≲ 3 mag and X_CO ≲ 10¹⁹. In addition, NGC 1333 lacks the decreasing portion of the X_CO versus A_V profile because of the missing data points with small I_CO ≲ 10 K km s⁻¹. Here we provide a description of the detailed comparison between our data of IC 348, B1E/B1, and NGC 1333 and the modified W10 model.

• 1.  
  For IC 348, the model with n_core = 10³ cm⁻³ reproduces well the observed shape of the X_CO versus A_V profile (decreasing X_CO at A_V ≲ 3 mag and increasing X_CO at A_V ≳ 3 mag).
• 2.  
  For B1E/B1, the model with n_core varying from 10³ cm⁻³ to 10⁴ cm⁻³ can reproduce the observed shape of the X_CO and A_V profile.
• 3.  
  For NGC 1333, the observed scatter at small A_V calls for a range of n_core ∼ 10^{3 − 4} cm⁻³. Considering that the models with n_core = 5 × 10³ cm⁻³ and 10⁴ cm⁻³ are essentially identical, however, the data points at A_V ≲ 3 mag with X_CO ≲ 10¹⁹ would not be reproduced by the model with n_core > 10⁴ cm⁻³. In addition, our observational data lack the decreasing portion of the X_CO versus A_V profile. We suspect that this is due to the limited spatial coverage of the COMPLETE I_CO image, which does not adequately sample low column density regions for NGC 1333 (only ∼10% of the data points have I_CO < 10 K km s⁻¹).

Zoom In Zoom Out Reset image size

In Figure 11(b), we compare the observed X_CO versus N(H₂) profiles with the model and find similar results. In summary, the modified W10 model with the "core–halo" structure and the input parameters appropriate for Perseus predicts the ranges of I_CO and N(H₂) in good agreement with our data. IC 348 and B1E/B1 are the best cases where the shape of the X_CO versus A_V profiles and the location of the minimum X_CO are well described by the model. We note that there are some discrepancies at low column densities in NGC 1333, where the data points with X_CO ≲ 10¹⁹ are not reproduced by the model and at the same time the observed data with X_CO ≳ 10²⁰ are missing because of the limited observational coverage.

Next, we plot N(H i) as a function of N(H) in Figure 11(c) and compare the profiles with the modified W10 model. As summarized in Section 2.1, Lee et al. (2012) found a relatively uniform N(H i) distribution across Perseus with ∼(8–10) × 10²⁰ cm⁻². Here we use the same N(H i) image as in Lee et al. (2012) and apply the same boundaries for the five regions as in Section 4. We find that the mean N(H i) varies from ∼7.4 × 10²⁰ cm⁻² (B5) to ∼9.6 × 10²⁰ cm⁻² (NGC 1333 and L1448). The model predicts N(H i) ∼ (9–9.6) × 10²⁰ cm⁻², with essentially no difference between n_core = 10³ cm⁻³ and 10⁴ cm⁻³ models. The predicted N(H i) distribution with ∼9 × 10²⁰ cm⁻² and its uniformity are consistent with what we observe in Perseus. This agreement will persist even if the N(H i) distribution is corrected for high optical depth H i. Our preliminary work on the effect of high optical depth H i that is missing in the H i emission image of Perseus shows that N(H i) increases by a factor of ∼1.5 at most because of the optical depth correction (the corrected N(H i) ∼ (8–18) × 10²⁰ cm⁻²; S. Stanimirović et al., in preparation). The ranges of the predicted N(H i) and N(H₂) distributions are comparable to what we find in Perseus. In Figure 11(d), we plot $R_{\rm H_2}$ against N(H) and indeed find that the model matches well our observations. In particular, the linearly increasing $R_{\rm H_2}$ with N(H) is reproduced well by the model, mainly driven by the uniform N(H i) distribution.

So far we made comparisons between the observations of Perseus and the modified W10 model with the "core–halo" structure. To investigate the role of the diffuse halo in determining H₂ and CO distributions, we show our data for IC 348 and predictions from the modified W10 model both with the "core–halo" structure and the uniform density distribution in Figure 12. The uniform density distribution simply assumes a dense core with n = 10³, 5 × 10³, or 10⁴ cm⁻³. Clearly, the "core–halo" model describes our data better. For example, the uniform density model underestimates the N(H i) distribution compared with the observed one across the cloud. In addition, it predicts the decreasing portion of the X_CO versus A_V profile shallower than our data, while reproducing the observed range of X_CO reasonably well. We compare the "core–halo" model with the uniform density model in detail as follows.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

N(H i) versus N(H): The uniform density model predicts N(H i) significantly smaller than what we measure across Perseus, N(H i) ∼ 9 × 10²⁰ cm⁻². The discrepancy ranges from a factor of ∼10–20 for n = 10³ cm⁻³ to a factor of ∼70–160 for n = 10⁴ cm⁻³. This large discrepancy results from the fact that H₂ self-shielding is so strong that almost all hydrogen is converted into H₂. On the other hand, the density of the halo is small enough that dust shielding is more important than H₂ self-shielding. To provide the sufficient dust shielding for H₂ formation, the entire halo remains atomic with its initial N(H i) ∼ 9 × 10²⁰ cm⁻², resulting in the uniform N(H i) distribution. We expect that if the density of the halo is significantly larger than the current n_halo = 40 cm⁻³, the halo will no longer be purely atomic because of the increased H₂ self-shielding.

N(H₂) versus N(H): All models predict the N(H₂) versus N(H) profile in good agreement with our data, even though the uniform density model slightly overestimates N(H₂) at small N(H). For example, the uniform density model with n = 10⁴ cm⁻³ predicts N(H₂) = 9.96 × 10²⁰ cm⁻² at N(H) = 2 × 10²¹ cm⁻², larger than our data by less than a factor of 2. However, this discrepancy is significant at such small N(H) and results in the small amount of N(H i) ≲ 10¹⁹ cm⁻². In addition, models with different densities predict essentially the same N(H₂) for a given N(H). All these results imply that neither density nor its distribution is critical for the H₂ abundance. Instead, N(H) primarily determines N(H₂).

$R_{\rm H_2}$ versus N(H): While the "core–halo" model reproduces both the range of $R_{\rm H_2}$ and the linear increase of $R_{\rm H_2}$ with N(H), the uniform density model overestimates $R_{\rm H_2}$ for a given N(H) by up to a factor of ∼300. This discrepancy mainly results from the significantly underestimated N(H i) in the uniform density model.

I_CO versus N(H₂): All models reproduce the observed I_CO versus N(H₂) profile reasonably well. In particular, both the "core–halo" and uniform density models with the smallest density show an excellent agreement with our data for IC 348. While the models with n ≳ 5 × 10³ cm⁻³ and n_core ≳ 5 × 10³ cm⁻³ predict larger I_CO at small N(H₂) (up to a factor of ∼10), the difference between the models with different densities becomes negligible at N(H₂) ≳ 1 × 10²¹ cm⁻², where I_CO saturates to ∼45–60 K km s⁻¹ for the uniform density model and ∼30–40 K km s⁻¹ for the "core–halo" model. All these results suggest that I_CO depends on density but only at small N(H₂) and changes in physical parameters other than density (e.g., v_turb) will be required to produce larger I_CO values once I_CO becomes optically thick.

I_CO versus A_V: While the "core–halo" model reproduces the sharp increase of I_CO observed at A_V ≳ 1 mag, the uniform density model predicts the increase of I_CO at A_V ≳ 0.6 mag much more gradually than our data. This difference comes from the fact that the uniform density model has larger density than the "core–halo" model, resulting in the larger I_CO for a given A_V ≲ 3 mag (n ≳ 10³ cm⁻³ for the uniform density model versus 〈n〉 ∼ 55–125 cm⁻³ for the "core–halo" model; Table 2). On the other hand, all models predict the saturation of I_CO to similar values at A_V ≳ 3 mag, suggesting that the larger density in the uniform density model no longer has a significant impact on I_CO because of the large optical depth of I_CO (n ≳ 10³ cm⁻³ for the uniform density model versus 〈n〉 ∼ 180–430 cm⁻³ for the "core–halo" model; Table 2).

X_CO versus A_V: All models reproduce the observed increase of X_CO at A_V ≳ 3 mag because they predict both the range of N(H₂) and the saturation of I_CO comparable to our data. On the other hand, the uniform density model shows the decrease of X_CO at A_V ≲ 3 mag much shallower than our data. This discrepancy mainly results from the less steep increase of I_CO predicted by the model at A_V ≲ 3 mag.

Summary: While we do not perform a full parameter space search, our comparison between the "core–halo" and uniform density models is illustrative and demonstrates that the diffuse halo is essential for reproducing the following observed properties: the uniform N(H i) distribution, the H₂-to-H i ratio for a given N(H), and the sharp increase of I_CO and decrease of X_CO at 1 mag ≲ A_V ≲ 3 mag. Considering that the uniform density model predicts the I_CO distribution extended toward smaller A_V, while producing the N(H₂) distribution in reasonably good agreement with our data (Figures 12(d) and (j)), we expect that the neglect of the diffuse halo will result in the underestimation of the size of "CO-free" H₂ envelope.

The S11 model is essentially composed of two parts. The first part is a modified version of the zeus–mp MHD code (Stone & Norman 1992; Norman 2000). Gas in a periodic box is set to have a uniform density distribution and is driven by a turbulent velocity field with uniform power 1 ⩽ k ⩽ 2, where k is the wavenumber. In addition, the magnetic field has initially orientation parallel to the z-axis, with a strength of 1.95 μG. To model the chemical evolution of the gas, Glover & Mac Low (2007a, 2007b), Glover et al. (2010), and Glover & Clark (2012) updated the zeus-mp MHD code with chemical reactions of several atomic and molecular species. The photodissociation of molecules by a radiation field is treated by the "six-ray approximation" method developed by Glover & Mac Low (2007a). The effect of self-gravity is not included. We refer to Glover & Mac Low (2007a, 2007b), Glover et al. (2010), and Glover & Clark (2012) for details on MHD, thermodynamics, and chemistry included in the S11 model. The second part is a three-dimensional radiative transfer code radmc–3d (C. P. Dullemond et al., in preparation).⁸ Once the simulated molecular cloud reaches a statistically steady state, radmc–3d is executed to model molecular line emission (e.g., CO). To solve the population levels of atomic and molecular species, radmc–3d implements the LVG method (Sobolev 1957), which has been shown to be a good approximation for molecular clouds (e.g., Ossenkopf 1997). We refer to Shetty et al. (2011a) for details on radmc–3d.

The MHD simulation follows the evolution of an initially atomic gas in a (20 pc)³ box with a numerical resolution of 512³. In this paper, we use the S11 model with the following input parameters: initial n = 100 cm⁻³, G = 1 $G_{0}^{\prime }$, ζ = 10⁻¹⁷ s⁻¹, Z = 1 Z_☉, and DGR = 5.3 × 10⁻²² mag cm². This simulation is essentially the same as the "n100 model" in S11 but has a higher numerical resolution and a simpler CO formation model based on Nelson & Langer (1999). We choose this particular simulation because it has a mass of ∼2 × 10⁴ M_☉, consistent with that of Perseus. The input parameters for the S11 model are reasonably close to what we expect for Perseus but not exactly the same as what we used for the modified W10 model. As it has been shown in S11 and Glover & Mac Low (2007b) that the simulated H₂ and CO column densities do not depend on small changes in G and ζ, this simulation would be appropriate for the comparison with our observations (Section 8.4.1 for details).

Compared with the modified W10 model, the final density distribution in the S11 model has a majority of the data points (∼99%) with n < 10³ cm⁻³, resulting in the small median density of ∼30 cm⁻³. Another important difference between the modified W10 model and the S11 model is that H₂ formation in the S11 model does not achieve chemical equilibrium until the end of the simulation. For example, Glover et al. (2010) found from their MHD simulations that the H₂ abundance primarily depends on the time available for H₂ formation and shows no indication of chemical equilibrium up to t ∼ 20 Myr. The gas will eventually become fully H₂ unless the molecular cloud is destroyed by stellar feedback such as photoevaporation by H ii regions and protostellar outflows. On the other hand, the CO abundance is controlled by photodissociation and reaches chemical equilibrium within t ∼ 2 Myr.

The final products of the S11 model include the N(H i), N(H₂), and I_CO images obtained at t ∼ 5.7 Myr. We smooth and regrid the simulated N(H i), N(H₂), and I_CO images so that they have both a spatial resolution of 0.4 pc and a pixel size of 0.4 pc. Recently, Beaumont et al. (2013) compared the COMPLETE data of Perseus with the S11 model and found that the S11 model systematically overestimates N(H₂) (e.g., Figure 5 of Beaumont et al. 2013). One of the possible explanations for this discrepancy is the different size between the simulation box and the individual regions in Perseus. Because the simulation box is larger than the individual regions in Perseus (20 pc versus ∼5–7 pc), the integrated quantities N(H i), N(H₂), and I_CO would need to be scaled. In the case of N(H i) and N(H₂), the scaling is straightforward under the assumption of isotropic density distribution, which is appropriate for the S11 model,⁹ and we simply need to account for the difference between the box and region sizes. However, estimating a proper scaling for I_CO is much more complicated because of the following reasons. First, the I_CO image was produced from the S11 model by integrating the CO brightness temperature, which was estimated by three-dimensional radiative transfer calculations, along a full radial velocity range. Second, the CO emission is optically thick in some parts of the simulation (∼10% of the volume). Rerunning the simulation with a smaller box does not solve the problem as molecular cores and clouds form out of initially larger-scale diffuse ISM. We therefore take an approach of determining the optimal LOS depth that minimizes the difference between our observations and the S11 model by investigating the N(H i) and N(H₂) images simultaneously. For the simulated I_CO image, on the other hand, we do not apply any scaling.

To do this, we estimate the difference between the observed mean and the simulated mean for each of N(H i) and N(H₂) with varying LOS depths. For example, we divide the simulated N(H i) and N(H₂) images by two to calculate the mean N(H i) and N(H₂) for the simulation with the LOS depth of 10 pc. We then normalize the difference by the observed mean of each quantity and calculate the sum of the two normalized differences in quadrature. The results are shown in Figure 13, and we find that the LOS depth that minimizes the difference between our data and the simulation products is 7 pc (Figure 13(c)). While the final quantity in Figure 13(c) has a broad minimum, it is encouraging that the estimated scale length is comparable to both the characteristic size of the five regions in Perseus and the total size of the slab for the "core–halo" model (Tables 1 and 2). As a double check that this scale length is appropriate, we use Larson's law established for turbulent molecular clouds from both observations and MHD simulations: $\sigma _{\rm CO}=(0.96 \pm 0.17) L_{\rm pc}^{0.59 \pm 0.07}$ km s⁻¹ (Heyer & Brunt 2004). For a region size of 20 pc we expect σ_CO ∼ 6 km s⁻¹, while for a region size of 7 pc we expect σ_CO ∼ 3 km s⁻¹. This level of CO velocity dispersion is in agreement with what is shown in Figure 6, confirming that scaling the simulation products to the LOS depth of 7 pc is reasonable. In summary, when we compare our observations with the S11 model, we scale the simulated N(H i), N(H₂), and N(H) images by multiplying them by 7/20 (Figures 14(a)–(c) and Figures 15(a) and (b)). On the other hand, because of the uncertainty in I_CO scaling, we use the original I_CO image produced by the S11 model (Figure 14(d) and Figures 15(c) and (d)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we apply the following thresholds to the simulated data to mimic the sensitivity limits of our observational data: N(H₂) > 3.3 × 10¹⁹ cm⁻² and I_CO > 0.09 K km s⁻¹ (our mean 1σ uncertainties calculated for the data points with N(H₂) > 0 cm⁻² and I_CO > 0 K km s⁻¹). This application of the thresholds to the S11 model is reasonable, considering the minimum N(H₂) ∼ 3.7 × 10¹⁹ cm⁻² and I_CO ∼ 0.2 K km s⁻¹ for the five regions in Perseus.

We first compare our data with the S11 model by constructing normalized histograms of N(H i), N(H₂), N(H), and I_CO in Figure 14. To construct the histograms, we use the data points with N(H₂) > 0 cm⁻² in Figure 2 ("All" histograms in black), as well as those shown in Figure 4 ("Subset" histograms in gray). While the gray histograms are limited to the regions where the CO emission is detected, the black histograms represent the whole Perseus cloud. The simulated data from the S11 model (smoothed, regridded, scaled for 7 pc, and the thresholds applied) are shown as green histograms. Note that the I_CO values from the S11 model are not scaled, and therefore the green I_CO histogram likely represents the upper limit of the actual histogram for sub-regions with a size of ∼7 pc (indicated as an arrow). Because the simulated data (except for I_CO) are scaled to match the properties of the five regions and the thresholds applied to the S11 model are comparable to the minimum N(H₂) and I_CO values of the five regions, the green histograms can be directly compared with the gray histograms. In comparison between our data and the S11 model, we find the following.

First, the black and gray N(H i) histograms are nearly identical. This results from the small variation in N(H i) across the whole Perseus cloud, as discussed in Section 2.1. The green histogram, on the other hand, has a peak at a factor of ∼2 smaller N(H i) and even more importantly a factor of ∼6 broader distribution than the observed data (the black and gray histograms).

Second, the gray and green N(H₂) histograms agree very well: both peak at a similar N(H₂), have a similar width, and show a lognormal-like distribution. The black histogram, on the other hand, is broader and has a tail toward small N(H₂). The difference between the black and gray histograms results from the existence of H₂ beyond the CO spatial coverage (e.g., "CO-dark" H₂ discussed in Lee et al. 2012).

Third, the green N(H) histogram peaks at a similar N(H) compared with the gray histogram, while showing a broader (a factor of ∼2) and lognormal-like distribution. The simulated distribution is broader mainly because the simulated N(H i) has a greater range than what is observed. The black and gray histograms, on the other hand, have a tail toward N(H) ≳ 10^21.4 ∼ 2.5 × 10²¹ cm⁻². This tail is consistent with Kainulainen et al. (2009), who found a deviation from the lognormal distribution at A_V ≳ 3 mag for Perseus (corresponding to N(H) ∼ 2.7 × 10²¹ cm⁻² with DGR = 1.1 × 10⁻²¹ mag cm²) and interpreted it as a result of self-gravity.

Lastly, because the simulated I_CO is not scaled for the LOS depth of 7 pc, we do not compare the exact shapes of the green and gray histograms but emphasize that the simulated I_CO becomes comparable to the observed I_CO only if we use the whole simulation box of 20 pc.

In summary, we find that the scaled S11 model reproduces the observed range of N(H₂) very well. While the predicted N(H i) has a relatively similar mean value compared with the observed N(H i), it has a broader distribution, and this leads to a broader range of N(H) in the simulation. The I_CO values from the S11 model, on the other hand, cannot be properly compared with our observations because of the nontrivial scaling of I_CO with the LOS depth. However, we find that the simulated I_CO is similar to the observed I_CO only when the CO emission is integrated for the full simulation box of 20 pc.

We plot N(H i) against N(H) for each dark and star-forming region and show predictions from the S11 model (smoothed, regridded, scaled for 7 pc, and the thresholds applied) in Figure 15(a). While the observed N(H i) ∼ 9 × 10²⁰ cm⁻² is in the range of the predicted N(H i), the relation between N(H i) and N(H) in the S11 model is different from what we find in Perseus: not only does the simulated N(H i) have a broader distribution, but the S11 model predicts a factor of ∼7 increase of N(H i) for the range of N(H) in Perseus, where we observe less than a factor of 2 variation in N(H i). This suggests that N(H i) linearly correlates with N(H) in the S11 model, and we indeed estimate Pearson's linear correlation coefficient r_p ∼ 0.8.

In addition, the S11 model predicts a factor of ∼2 smaller N(H i) for a given N(H) on average. As a result, $R_{\rm H_2}$ is slightly larger in the S11 model for a given N(H) and increases with N(H) with a slope smaller than what we observe (Figure 15(b)). While our observations show $R_{\rm H_2}$ < 1 for the outskirts of the five regions, the simulation has $R_{\rm H_2}$ > 1 everywhere, even for the regions with small n < 10² cm⁻³.

Next, we plot the observed I_CO as a function of A_V and show the S11 model in Figure 15(c). As discussed in Section 7.2.1, the simulated N(H i) and N(H₂) data can be scaled for the five regions in Perseus, while the simulated I_CO data cannot. To properly examine the relation between I_CO and A_V in the S11 model, therefore, we show the predicted I_CO versus A_V profile without applying the scaling and the thresholds and focus on only the general shape of the profile. We find that the S11 model describes the relation between I_CO and A_V reasonably well: a steep increase of I_CO at small A_V and a hint of the saturation of I_CO at large A_V. Interestingly, the S11 model predicts that I_CO increases with a large scatter at small A_V.

Finally, we show the X_CO versus A_V profile for each dark and star-forming region in Figure 15(d) with the S11 model. As in Figure 15(c), the unscaled N(H i), N(H₂), and I_CO data are used for this comparison. We find that the S11 model predicts a sharp decrease of X_CO at small A_V and a gradual increase of X_CO at large A_V. While a quantitative comparison is not possible without scaling, the simulated data show the characteristic relation between X_CO and A_V in broad agreement with the observational data (particularly for IC 348 and B1E/B1). This is consistent with Shetty et al. (2011a), who performed a number of MHD simulations (n = 100, 300, 1000 cm⁻³ and Z = 0.1, 0.3, 1 Z_☉) and found a steep decrease of X_CO at A_V ≲ 7 mag and a steady increase of X_CO at A_V ≳ 7 mag for all simulations probing a large range of interstellar environments. Relative to the observations, we find that the simulated X_CO shows a significantly larger scatter at small A_V, while the scatter becomes more comparable to what is found in the observations at the high end of the A_V range.

The modified W10 model that is comparable to the observations of Perseus uses the "core–halo" structure motivated by previous studies of molecular clouds (Section 7.1.1). We showed that the model with a uniform density distribution predicts N(H i) much smaller than the uniform N(H i) ∼ 9 × 10²⁰ cm⁻² measured across Perseus. The uniform density model with the largest density n = 10⁴ cm⁻³ predicts the smallest N(H i), up to a factor of ∼160 smaller than what is observed. The main reason is that in the uniform density model, H₂ self-shielding alone counteracts H₂ photodissociation by LW photons. Traditionally, it has been known that G/n determines whether H₂ self-shielding or dust shielding is more important for H₂ formation and controls the location of the transition from H i to H₂ in a PDR (e.g., Hollenbach & Tielens 1997). With G = 0.5 $G_{0}^{\prime }$ and n = 10³ cm⁻³ in the uniform density model, G/n is 5 × 10⁻⁴ cm³, small enough that dust shielding is negligible. In this case, most of the H i is converted into H₂ because of the strong H₂ self-shielding. On the other hand, the "core–halo" model with n_core = 10³ cm⁻³ and n_halo = 40 cm⁻³ has G/n_halo ∼ 0.01 cm³ in the cloud outskirts. This increased G/n makes H₂ self-shielding less important for H₂ formation, and as a result, the gas remains atomic with N(H i) ∼ 9 × 10²⁰ cm⁻².

The fact that the modified W10 model needs a diffuse H i halo to reproduce the observed N(H i) suggests that dust shielding is important for H₂ formation in Perseus. This importance of dust shielding is consistent with what Lee et al. (2012) found from their comparison with the K09 model. The K09 model investigates the structure of a PDR in a spherical cloud on the basis of H₂ formation in chemical equilibrium and predicts the following variable as one of the key parameters that determine the location of the transition from H i to H₂:

$$\begin{equation} \chi = 2.3\frac{(1 + 3.1 Z^{\prime 0.365})}{\phi _{\rm CNM}}, \end{equation} \tag{ 3 }$$

where Z' is the metallicity normalized to the solar neighborhood value and ϕ_CNM is the ratio of the actual CNM density to the minimum CNM density at which the CNM exists in pressure balance with the warm neutral medium (WNM). This χ is the ratio of the rate at which LW photons are absorbed by dust grains (dust shielding) to the rate at which they are absorbed by H₂ (H₂ self-shielding) and is conceptually similar to G/n. K09 predicts χ ∼ 1 in all galaxies where the pressure balance between the CNM and the WNM is valid, suggesting that dust shielding and H₂ self-shielding are equally important for H₂ formation. By fitting the K09 model to the observed $R_{\rm H_2}$ versus $\Sigma _{{\rm H}\,\scriptsize{I}}$+$\Sigma _{\rm H_2}$ profiles, Lee et al. (2012) indeed found χ ∼ 1 for Perseus.

In the modified W10 model, a diffuse H i halo is also required to reproduce the observed steep increase of I_CO at A_V ≳ 1 mag and the sharp decrease of X_CO at A_V ≲ 3 mag (Section 7.1.3). The uniform density model predicts the shallower increase of I_CO at smaller A_V ≳ 0.6 mag, suggesting a less sharp transition from C ii/C i to CO located closer to the surface of the gas slab. The more extended CO distribution eventually results in the reduced "CO-free" H₂ envelope, and therefore the uniform density model with n = 10⁴ cm⁻³ would have the smallest amount of "CO-dark" H₂. The CO distribution deep inside of the gas slab, on the other hand, does not appear to be affected by the presence of the diffuse H i halo because of the saturation of I_CO.

Even though the modified W10 model with the "core–halo" structure reproduces the observed N(H i), N(H₂), and I_CO distributions, the agreement is likely to remain only if the halo density is not significantly larger than 40 cm⁻³. The current density n_halo = 40 cm⁻³ originates from the theoretical (e.g., Wolfire et al. 2003) and observational (e.g., Heiles & Troland 2003) properties of the CNM. While large H i envelopes associated with molecular clouds have been frequently observed (e.g., Knapp 1974; Wannier et al. 1983, 1991; Reach et al. 1994; Rogers et al. 1995; Williams & Maddalena 1996; Imara & Blitz 2011; Lee et al. 2012), a number of fundamental questions still remain to be answered. For example, what are the physical properties of the H i halos, such as density, temperature, and pressure? What is the ratio of the CNM to the WNM in the halos? Is there any correlation between the ratio and the H₂ abundance and star formation? Are the halos expanding or infalling? The high-resolution H i data from the GALFA-H i survey will be valuable for future studies of the extended H i halos around Galactic molecular clouds in a wide range of interstellar environments. Finally, further comparisons between observations and theoretical models will be important to fully constrain the parameter space and density structure of the H i halos.

The timescale of H₂ formation on dust grains, $t_{\rm H_2}$, dominates chemical timescales of PDRs (e.g., Hollenbach & Tielens 1997). For the modified W10 model with the "core–halo" structure, dense regions have n_core ≳ 10³ cm⁻³ where gas is completely molecular ($n_{\rm H_2}$ ∼ 0.5n). In this case, $t_{\rm H_2}$ = 0.5/$\mathcal {R}$n_core ≲ 0.5 Myr, where $\mathcal {R}$ = 3 × 10⁻¹⁷ cm³ s⁻¹ is the rate coefficient for H₂ formation (Wolfire et al. 2008). In diffuse regions with n_halo = 40 cm⁻³, on the other hand, gas is mostly atomic ($n_{\rm H_2}$ ∼ 0.1n) and therefore $t_{\rm H_2}$ = 0.1/$\mathcal {R}$n_halo ∼ 2.6 Myr. Because $t_{\rm H_2}$ of the model is well within the expected age of Perseus, t_age ∼ 10 Myr, the assumption of chemical equilibrium is valid. In other words, Perseus is old enough to reach chemical equilibrium, and therefore it is not surprising that the equilibrium chemistry model (W10) fits our observations very well.

However, for steady state chemistry to be valid, $t_{\rm H_2}$ ≲ t_age is not enough: $t_{\rm H_2}$ should be short compared with the dynamical timescale of a molecular cloud, t_dyn. For Perseus, this requires t_dyn ≳ 3 Myr. As a rough estimate, we calculate a crossing timescale, t_cross = L/σ ∼ 10 pc/1.8 km s⁻¹ ∼ 6 Myr, where we choose L as the characteristic size of the individual regions in Perseus and σ as the mean CO velocity dispersion. This t_cross ∼ 6 Myr satisfies the condition for t_dyn ≳ 3 Myr. However, many dynamical processes are involved with the formation and evolution of molecular clouds, (e.g., cloud–cloud collisions, spiral shocks, stellar feedback; Mac Low & Klessen 2004; McKee & Ostriker 2007) and therefore it is difficult to pin down the exact process that is most relevant for the formation of molecular gas. The good agreement between our data and the modified W10 model with the "core–halo" structure suggests that the characteristic t_dyn for the formation of molecular gas in Perseus should be ≳3 Myr.

In Section 7.2.2, we found that the scaled S11 model predicts N(H₂) comparable to the estimated N(H₂) in Perseus. This excellent agreement will likely hold even if some of the input parameters slightly change. For example, the S11 model was run with G = 1 $G_{0}^{\prime }$, and this is a factor of ∼2 stronger than what we measure across Perseus. Considering that S11 found no considerable difference in N(H₂) for their models with G = 1 $G_{0}^{\prime }$ and 10 $G_{0}^{\prime }$ (Section 3.1 of S11), however, decreasing G from 1 $G_{0}^{\prime }$ to 0.5 $G_{0}^{\prime }$ to match the property of Perseus will not make a significant change in N(H₂). In addition, increasing ζ from 10⁻¹⁷ s⁻¹ to 10⁻¹⁶ s⁻¹ to be consistent with the modified W10 model will not affect N(H₂) very much on the basis of the fact that Glover & Mac Low (2007b) found a negligible change in N(H₂) when ζ increased from 10⁻¹⁷ s⁻¹ to 10⁻¹⁵ s⁻¹ in their MHD simulation with initial n = 100 cm⁻³ (Section 6.3 of Glover & Mac Low 2007b). Increasing DGR from 5.3 × 10⁻²² mag cm² to 1.1 × 10⁻²¹ mag cm² for Perseus will lead to more rapid H₂ formation, but the model with the increased DGR will not be substantially different from the current S11 model since the S11 model becomes H₂-dominated rapidly by t ∼ 3 Myr (Figure 7 of Glover & Mac Low 2011). Finally, the extension of the simulation run up to t ∼ 10 Myr, comparable to the age of Perseus, will not significantly increase N(H₂), considering that Glover & Mac Low (2011) found only a factor of ∼1.3 increase of the mass-weighted mean H₂ abundance from t ∼ 5 Myr to t ∼ 10 Myr for their MHD simulation with initial n = 100 cm⁻³ (Section 3.3 of Glover & Mac Low 2011).

Similarly, small changes in the model parameters will likely make no substantial difference in I_CO. For example, S11 showed that increasing G from 1 $G_{0}^{\prime }$ to 10 $G_{0}^{\prime }$ does not change I_CO for those regions where CO is well shielded against the radiation field (Section 3.1 of S11). Therefore, decreasing G from 1 $G_{0}^{\prime }$ to 0.5 $G_{0}^{\prime }$ will make only a minor change in I_CO at large column densities. Increasing the current DGR of 5.3 × 10⁻²² mag cm² by a factor of ∼2 will cause more rapid CO formation, but I_CO will not be significantly influenced because CO formation in the S11 model reaches chemical equilibrium rapidly by t ∼ 2 Myr. Lastly, we do not expect that running the S11 model up to t ∼ 10 Myr drastically increases I_CO, considering that the MHD simulation with initial n = 100 cm⁻² in Glover & Mac Low (2011) predicts only a factor of ∼2 increase of the mass-weighted mean CO abundance from t ∼ 5 Myr to t ∼ 10 Myr (Section 3.3 of Glover & Mac Low 2011). Note that changes in CO abundance at t > 2 Myr are stochastic fluctuations after chemical equilibrium is achieved.

We therefore conclude that the input parameters used in the S11 model are reasonable for the comparison with the observations of Perseus and small (a factor of few) changes in the input parameters will not result in significant changes in N(H i), N(H₂), and I_CO. Considering that Perseus has most likely reached chemical equilibrium, it provides a suitable testbed for investigating whether results from the time-dependent MHD simulation converge to the time-independent PDR model for molecular clouds that are evolved enough.

As shown in Section 7.2.2, the scaled S11 model produces the N(H₂) distribution in excellent agreement with our observations as well as the modified W10 model. This suggests that the time-dependent H₂ formation model (S11) is consistent with the time-independent H₂ formation model (W10) for a low-mass, old molecular cloud such as Perseus. Our result agrees with Krumholz & Gnedin (2011), who found that time-dependent effects on H₂ formation become important only at extremely low metallicities Z ≲ 10⁻² Z_☉. While the median N(H i) in the S11 model is also in reasonably good agreement with the observations, the simulated N(H i) distribution is a factor of ∼6 broader than the observed one and particularly shows a more extended tail toward small N(H i) ≲ 3 × 10²⁰ cm⁻². This broad N(H i) distribution in the MHD simulation likely results from strong compressions and rarefactions by turbulence, and the predicted N(H i) is on average a factor of ∼2 smaller than the observed N(H i) for a given N(H). The discrepancy becomes significant at small N(H) ∼ 10²¹ cm⁻², where the S11 model underestimates N(H i) by up to a factor of ∼10. Finally, the S11 model predicts that N(H i) increases with N(H), suggesting no minimum N(H i) beyond which the rest of hydrogen is converted into H₂. In the modified W10 model with the "core–halo" structure, on the other hand, the diffuse halo remains atomic with N(H i) ∼ 9 × 10²⁰ cm⁻², and the dense core is fully converted into H₂. Clearly, this discrepancy in N(H i) between the simulation and the observations is significant and interesting. One potential avenue in exploring this in the future is by using a mixture of neutral phases for initial conditions, mimicking in some way the "core–halo" structure in the modified W10 model.

In the case of I_CO, we could not properly compare the S11 model with our observations because of the nontrivial scaling of I_CO for different line of sight depths. Instead, we found that the simulated I_CO becomes comparable to the observed I_CO only if the CO emission is integrated for the full simulation box of 20 pc. This suggests that the S11 model likely underestimates I_CO for the conditions relevant to Perseus. Interestingly, we estimate N(CO) ∼ 1 × 10¹⁷ cm⁻² for B5, IC 348, B1, and NGC 1333 by using the ¹³CO(J = 1 → 0) excitation temperatures, optical depths, and integrated intensities provided by Pineda et al. (2008) and assuming N(CO) = 76N(¹³CO) (Lequeux 2005). This value is in reasonably good agreement with the simulated mean N(CO) ∼ 5 × 10¹⁶ cm⁻² (calculated from the smoothed, regridded, scaled, and thresholds applied S11 model).¹⁰ This comparison suggests that the potential discrepancy in I_CO between the observations and the S11 model would result from the radiative transfer calculations and/or the difference in velocity range. The velocity range of the CO emission, Δv, directly affects I_CO via I_CO = ∫T_Bdv, and a smaller Δv would result in a smaller I_CO for the same T_B.

While we could not compare specific I_CO values predicted by the S11 model with our observations, we found that the S11 model reproduces the observed shape of the I_CO versus A_V profiles reasonably well. This suggests that penetration of UV photons into the ISM, dust shielding, and self-shielding against the ISRF are relatively well captured in the CO formation process by S11. In addition, we noticed that I_CO has a much larger scatter at small A_V. This could result from large density fluctuations in the turbulent medium. The gas in the S11 model would be strongly compressed and rarefied by turbulence, and the gas density at a given A_V can vary over several orders of magnitude (e.g., Figure 14 of Glover et al. 2010). In this case, CO can form in dense clumps even at small A_V, and I_CO therefore shows a large scatter. This scatter is reduced at large A_V where I_CO eventually saturates. Finally, turbulent mixing could spread the CO distribution, contributing to the large scatter of I_CO.

In general, our study shows that the scaled MHD simulation by S11 is successful in reproducing N(H₂) in Perseus, which is a low-mass, old molecular cloud most likely in chemical equilibrium. On the other hand, future model adjustments are required to better match the observed N(H i) and I_CO. We have revealed two important areas of future attention: (1) the role of diffuse halos in the formation of molecular gas and (2) the effect of density fluctuations and turbulent mixing in the spatial distribution of molecular gas. To investigate these two issues, we plan to compare observations of several Galactic molecular clouds with MHD simulations that explore different fractions of neutral phases and a varying degree of turbulence as initial conditions. In particular, our future work will include molecular clouds less evolved and/or forming more massive stars (and therefore more turbulent) than Perseus, where the difference between the MHD and PDR models is likely to be more pronounced.