Free Access
Issue
A&A
Volume 544, August 2012
Article Number A22
Number of page(s) 13
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201218865
Published online 19 July 2012

© ESO, 2012

1. Introduction

The interstellar medium (ISM) is a complex system: its structure and dynamics are governed by the interaction of many processes that cover a wide range of scales, involving both micro- and macrophysics. To understand how the ISM works is also essential for understanding star and planetary system formation. This field has seen much progress in recent years, both on the observational side, with the results from the Herschel Space Observatory (Pilbratt et al. 2010; de Graauw et al. 2010), and on the theoretical side, with ever-improving numerical simulations (see e.g. Banerjee et al. 2009; Vázquez-Semadeni et al. 2010), which now consistently treat self-gravity, thermodynamics and magnetohydrodynamics (MHD).

The new challenge for numerical simulations of the ISM is now to incorporate some treatment of chemistry and grain physics, to be able to compare the results with observations of atomic and molecular lines. One possibility, followed for instance by Glover et al. (2010), is to perform multi-fluid simulations, each fluid representing a chemical species and being coupled to the others via a network of reactions. This approach has the advantage that it naturally solves for the time-dependent distribution and dynamics of the various species in three-dimensional space, but its computational cost requires making use of simplifying assumptions, most notably a small number of species and reactions. For instance, Glover et al. (2010) treated a network of 218 reactions for 32 species, 13 of which were assumed to be in instantaneous chemical equilibrium, so that only 19 were fully treated as time-dependent quantities.

In this paper, we follow a different approach: we post-processed a single-fluid MHD simulation with the Meudon PDR code. This has the advantage of providing a full chemical network (99 species, 1362 reactions), but on the other hand, it implies a one-dimensional, steady-state treatment. With this approach, our main goal is to discuss the effects of realistic density fluctuations on photon-dominated region (PDR) chemistry, since previous studies have focused on uniform density models, e.g. Le Petit et al. (2006), or on simplified clumpy models, such as those by Wolfire et al. (2010). A second goal of this study is to estimate the amount of “dark molecular gas” that can be expected from observations of the diffuse ISM, i.e. gas where hydrogen is mostly in the form of H2, but where CO is too scarce to be seen in the usual tracer that is the (J = 1 → 0) line, given current sensitivities. We wish to compare this estimate with available observations (Grenier et al. 2005; Leroy et al. 2007; Abdo et al. 2010; Velusamy et al. 2010) and models (Wolfire et al. 2010).

The paper is organised as follows: Sect. 2 describes the MHD simulation used and presents the PDR code in its current state, while Sect. 3 gives an overview of the post-processing method used to combine the two. Section 4 discusses the lines of sight selected for running our analysis. Results and comparison to observations are presented in Sect. 5. Section 6 offers a more detailed discussion and a summary. Details on using the Meudon PDR code with density profiles and our strategy for post-processing results can be found in the appendices at the end of the paper.

2. Tools

2.1. The MHD simulation

As a model of the turbulent diffuse ISM, we used data from a MHD simulation performed by Hennebelle et al. (2008) using the RAMSES code (Teyssier 2002; Fromang et al. 2006). A notable advantage of this code is its adaptive mesh refinement (AMR) capabilities, which enable to locally reach extremely high spatial resolutions. The initial setup of the simulation is a cube of homogeneous warm neutral gas (WNM) of atomic gas (a 10:1 mixture of H and He) with density nH = 1 cm-3 and temperature T = 8000   K. The cube is L = 50   pc on each side, and two converging flows of WNM are injected from opposing faces along the X axis with relative velocity ΔVX ~ 40   km   s-1, which means that the Mach number of each flow with respect to the ambient WNM is ℳ ~ 2. Transverse and longitudinal velocity modulations are imposed, with amplitudes roughly equal to the mean flow velocity (20   km   s-1). Periodic boundary conditions are applied on the remaining four faces. The simulation starts on a 2563 grid, with two extra levels of refinement based on density thresholds, so that the effective resolution of the simulation is  ~0.05 pc. For our purposes, we regridded the data regularly on a 10243 cube, therefore using the maximum resolution over the whole domain. A magnetic field is present in the simulation, and is initially parallel to the X axis with an intensity of about 5 μG, consistent with observational values at these densities (Crutcher et al. 2010). The converging flows collide near the midplane X = 0 about 1   Myr into the simulation, and as the total mass grows from  ~3000   M initially to about 10 times this value at the end of the simulation (t = 13.11   Myr), gravity eventually takes over and, combined with the effects of thermal instability (Field 1965; Hennebelle & Pérault 2000; Koyama & Inutsuka 2002; Heitsch et al. 2005; Hennebelle & Audit 2007; Banerjee et al. 2009), leads to the formation of cold dense clumps (nH > 100 cm-3, while T ~ 10−50   K) within a much more diffuse and warm interclump medium (nH ~ 1−10 cm-3 and T ~ 103−104   K). This occurs at t ~ 12   Myr. See Hennebelle et al. (2008) for a more detailed description. Figure 1 shows the column density of the gas viewed along the X axis, which is that of the incoming flows.

2.2. The Meudon PDR code

The Meudon PDR code (http://pdr.obspm.fr/) is a publicly available set of routines (Le Bourlot et al. 1993; Le Petit et al. 2006; Gonzalez Garcia et al. 2008) whose purpose is to describe the UV-driven chemistry of interstellar clouds, and in particular of PDR. It is a steady-state one-dimensional code in which a plane-parallel slab of gas and dust is illuminated on either or both sides by the light from a star or by the standard interstellar radiation field (ISRF), which is defined using expressions from Mathis et al. (1983) and Black (1994). Usually run with homogeneous gas densities, the code can also accept input of density profiles via a text file. At each point along the line of sight, radiative transfer in the UV is treated to solve for the H/H2 transition, using either the approximations by Federman et al. (1979), or an exact method based on a spherical harmonics expansion of the specific intensity (Goicoechea & Le Bourlot 2007). A number of heating (photoelectric effect on grains, cosmic rays) and cooling (infrared and millimeter emission lines) processes contribute to the computation of thermal balance. Outputs of the code include gas properties such as temperature and ionisation fraction, radiation energy density, chemical abundances and column densities, level populations and line intensities. The code is iterative and therefore requires the user to check whether convergence has been achieved.

thumbnail Fig. 1

Total gas column density along the X axis of the MHD simulation snapshot used here. The box in the upper left side marks the position of the clump selected for running our analysis.

Open with DEXTER

3. Method overview

Naively, one would like to run the PDR code on all lines of sight through the simulated cube to derive a three-dimensional chemical structure, as well as line-of-sight integrated observables (an emission map for the CII [158 μm] line, for instance). However, this is neither feasible nor desirable.

Two computational reasons preclude this brute force approach. Firstly, the cost is prohibitive: just one global iteration of the PDR code for a single run typically completes in a couple hours on a GNU/Linux machine equipped with four dual-core 64-bit x86 processors and 64GB of memory. Because a run usually converges in some ten iterations, it takes about a day to process a single line of sight. The PDR code, however, has been ported on the EGEE grid1, which allows us to perform  ~100 runs simultaneously in the same timeframe. Still, it was not reasonable to consider treating more than  ~103 lines of sight in this work.

The second computational problem to consider is that the PDR code has trouble converging in low-density regions, of which there are many in the MHD simulation. These are the vicinities of X =  ±25 pc, where WNM gas (nH = 1 cm-3) is entering the box, but low-density regions are found everywhere throughout the cube (the volume filling factor of regions where is fv = 0.96). This convergence problem might be related to the fact that Ly-α emission, which is a main cooling process for nH ≲ 5 cm-3, is not included in the code. This renders the computation of thermal balance by the PDR code in these regions unreliable.

Besides these computational hurdles, it should be noted that the code is one-dimensional, and treats a density profile as that of a plane-parallel slab of gas. Imagine then that a line of sight intercepts a dense clump of matter: the gas lying directly behind that clump will be shielded from incoming radiation, since the most energetic UV photons will have been absorbed by the clump. This leads to shadowing artefacts in regions that in reality may well be illuminated from other directions, because the ISM has a complex, fractal-like, and evolving structure.

In this paper, we deal with these artefacts in the following way: the physical conditions and chemical composition at every grid point (X,Y,Z) considered in the analysis may be obtained by running the PDR code in p directions going through that grid point, for instance along the main coordinate axes. Thus, each quantity F(X,Y,Z) output by the code has p possible values f1,f2,...,fp. This is the case, in particular, of the radiation energy density E, which has possible values e1,e2,...,ep. To combine the p runs at this grid point, we selected a direction p0 for which the radiation energy density there is maximum, . This choice is discussed in Sect. 6. Of course p0 is a function of (X,Y,Z). We then chose F(X,Y,Z) = fp0 for every quantity F output by the code. This procedure takes better account of the porosity of the simulated structures to the ISRF, while ensuring element conservation at each grid point. Ideally, the more directions p, the better, but this obviously comes with an increased computational cost. In this paper, we chose p = 2 as a compromise, running the PDR code in two orthogonal directions.

The analysis in this paper focuses on a small subset of a single simulation snapshot, which we discuss in the next section. After applying a density threshold n0 = 20 cm-3, which is the lowest value considered in the grid of models run by Le Petit et al. (2006) and therefore deemed sufficient to ensure convergence, we extracted a number of one-dimensional density profiles from this thresholded subset. We ran the PDR code on these profiles and combined results to derive the chemical structure of the clump.

4. Lines of sight selected for PDR computations

thumbnail Fig. 2

Maximum gas density nH along the lines of sight parallel to the X axis within the selected clump. Contours show the total gas column density NH from 3 × 1021 cm-2 to 1.1 × 1022 cm-2 in steps of 1021 cm-2. The “clump” displayed here has a size roughly 1.5   pc    ×    2.5   pc. The two-dimensional slab of gas under study is seen projected as the single-pixel-wide AB strip. The white triangle marks the position of the maximum value of max(nH) in this region, and the white circle that of the maximum of NH. They are separated by 0.65 pc.

Open with DEXTER

4.1. Simulation snapshot

The snapshot chosen to run our analysis on was timed at 7.35 Myr, when the densest parts of the cloud reach nmax ~ 9 × 103 cm-3. Some of the structures present in the simulation at this time are self-gravitating, but we are confident that they are still diffuse enough that the simulation snapshot is representative of a non-starforming region of the ISM. We therefore ran our analysis without any illuminating star, with the ISRF being the only source of primary UV photons.

4.2. Selected clump

To select a representative subset, we note that observational PDRs such as the Horsehead Nebula (see e.g. Pety et al. 2007) are found at the edge of dense and cold clouds of gas and dust, illuminated by ambient FUV light and possibly nearby young stars. It thus makes sense to focus on a “clump”, defined observationally as a connected structure with a significantly higher column density than its surroundings. We identified clumps via a friend-of-friend algorithm on the column density map along the X axis (Fig. 1), using a threshold N0 = 3 × 1021 cm-2, which corresponds to a mean density ⟨nH⟩ = n0 = 20 cm-3 over 50 pc.

The selected clump, which lies at the top left corner of the simulation’s field, harbours an interesting feature: Fig. 2 shows in colour scale the map of the maximum gas density max(nH) encountered along the X axis for every line of sight within the clump. It so happens that the peak of NH does not match that of max(nH), or even a local maximum of the latter. This is important to note for species that may be sensitive to the local gas density rather than to the total column density. Properties of that selected clump are listed in Table 1.

Table 1

Properties of the selected observational clump.

4.3. Selected lines of sight

The observational clump shown in Fig. 2 has a size ΔY × ΔZ ≃ 1.5   pc    ×    2.5   pc. In the X direction, most of its gas is located in the central ΔX  ≃  15 pc around X = 0. Given the pixel size δ ≃ 0.05   pc, applying the PDR code on a structure of that size along the three coordinate axes requires some 25 000 runs, which is beyond the scope of this work. Consequently, we restricted our study to a two-dimensional slab of gas across the observational clump. It is projected in Fig. 2 as the one-pixel-wide strip AB, which is therefore  ~3.7   pc long,  ~0.05   pc wide and contains 76 lines of sight along the X axis. These sample a wide range of column densities, from NH = 6.11 × 1020 cm-2 (near the A end) to NH = 1.11 × 1022 cm-2, corresponding to visual extinctions AV = 0.33 to AV = 5.93, using the conversion from total hydrogen column density NH = N(H) + 2N(H2) (1)with RV = 3.1 and CD = 5.8 × 1021 cm-2   mag-1 (see Table A.1 and Le Petit et al. 2006). The structure of the gas on these lines of sight (Fig. 3) is complex, with many small, dense and cold regions (nH ≃ 103 cm-3, T ≃ 20 K) interconnected via a filamentary structure and embedded within a much more diffuse and warm medium (nH ≃ 1 cm-3, T ≃ 104 K). The overdense regions () are located near the midplane of the simulation (−9   pc ≲ X ≲ 2   pc), where the flows collide and the gas condenses into cold structures.

It is this subset of the simulated cube, shown in Fig. 3, which is the focus of our study, and from which we extracted one-dimensional density profiles. As Fig. 4 shows, this subset indeed contains most of the gas on the lines of sight within the AB strip: Except in the outermost regions in which , column densities along X over that region represent more than half the total column densities over the full 50 pc lines of sight. Figure 4 also shows that whithin the AB strip the peaks of NH and max(nH) are still separated by about 0.55 pc.

thumbnail Fig. 3

Structure of the gas in the two-dimensional slab under study. This is a close-up view on the region X ≃ 0 in which the incoming flows of WNM collide to form cold structures, and only regions in which are shown. The colour image shows the gas density nH in logarithmic scale, while contours show the gas temperature TMHD at 20 K (higher densities), 30 K and 50 K (lower densities). The two-dimensional projection of the velocity field (VX,VY) is also shown as yellow arrows, the length of which indicate the velocity modulus at that location (at the centre of each arrow). The A and B extremities of the observational one-dimensional strip are indicated for reference, and the dashed line marks an example location for the extracted density profiles on which the PDR code is run (Fig. 5). The grey areas are outside of the domain used for PDR computations.

Open with DEXTER

thumbnail Fig. 4

Total gas column densities along the X axis within the AB strip under study (top plot) and maximum gas density max(nH) on the same lines of sight (bottom plot). Shown are the column densities over the full 50 pc lines of sight along the X axis (dashed line) and the column densities for the overdense regions shown in Fig. 3 (solid line). Grey areas mark lines of sight for which less than 50% of the mass is in the overdense region. The dash-dotted lines mark the positions of the maxima of NH and max(nH), which are separated by 0.55 pc.

Open with DEXTER

thumbnail Fig. 5

Example of a density profile used in the PDR code. This is the profile extracted at the location of the dashed line in Fig. 3.

Open with DEXTER

As explained in Sect. 3, we considered two (p = 2) possible directions for the one-dimensional density profiles extracted from the two-dimensional subset, namely those parallel to the X or Y axis. The dashed line in Fig. 3 marks the location of such an extracted profile, which is shown in Fig. 5. Note that we considered profiles to be connectedly overdense, which means that along each line of sight, we may extract several profiles separated by underdense regions nH < n0. For instance, this is the case of the line of sight parallel to the X axis located at Y = −23   pc.

We thus extracted 447 density profiles, of which 156 were parallel to the X axis and 291 were parallel to the Y axis. Their statistical properties are summarised in Tables 2 and 3, emphasising the wide dynamic range they sample in column densities (~104) and mean densities (~30). The high values in density-weighted temperatures correspond to density profiles that never much deviate from n0 = 20 cm-3. For each of these 447 profiles, we ran the PDR code assuming identical illumination on both sides. Since one of the objectives of this paper is to assess the effect of realistic density distributions along the line of sight on the chemical composition of interstellar clouds, we also applied the PDR code to a homogeneous reference model for each extracted profile. We specify this model, called uniform, as having the same mean density ⟨nH⟩ and total visual extinction AV as the inhomogeneous model, which we dub los from now on. Unless otherwise specified, results presented in the next section refer to these los models. The setup for all runs is detailed in Appendix A and their post-processing is described in Appendix B.

5. Results

5.1. Temperature comparison

The PDR code and the MHD simulation both treat thermal balance, so that we have two estimates of the gas temperature, which we can compare: Fig. 6 shows the ratio r = TPDR/TMHD of the gas temperature TPDR output by the PDR code to the gas temperature TMHD computed in the MHD simulation at every point in the subset under study versus the total gas density nH at that point. Average ratios ⟨r⟩ in selected density bins are also shown. It appears quite clearly that ⟨r⟩ is close to 1, with 0.3 ≲ ⟨r⟩ ≲ 2.0 over the whole range of densities.

The fact that TPDR ~ TMHD comes as a pleasant surprise, considering the differences between the PDR and MHD computations: while the former is one-dimensional, steady-state, and includes cooling via the infrared and submillimeter lines from atomic and molecular species, especially H2 transitions (Le Petit et al. 2006), the latter is three-dimensional, dynamical, and only includes cooling via the fine structure lines of CII [158 μm] and OI [63 μm], as well as the recombination of electrons with ionised PAHs. This suggests that unless a very precise knowledge of the temperature is needed, it is probably not necessary, at least as a first approximation and in the range of densities and temperatures probed here, to refine the details of cooling processes in our MHD simulations, because the simple cooling function currently used already yields gas temperatures close to those found using the more detailed processes of the PDR code. Similar conclusions were reached by Glover et al. (2010).

Table 2

Properties of the 156 one-dimensional profiles extracted parallel to the X axis.

Table 3

Same as Table 2 but for the 291 one-dimensional profiles extracted parallel to the Y axis.

thumbnail Fig. 6

Ratio of the gas temperatures TPDR/TMHD versus total gas density nH. Grey crosses show all points, while black circles are average values over density bins in logarithmic scale, with error bars standing for  ±1σ. The dashed line corresponds to TPDR/TMHD = 1.

Open with DEXTER

5.2. Chemical structure and comparison to observations

The spatial distributions of H, H2, C+, C, CO, CS, CH, and CN in the simulation subset are shown in Figs. 7 and 8. In these figures, points were clipped where the density n(α) of species α was below n0 = Xαn0, with Xα a typical threshold abundance for the detection of that species (see Table 4). Although shadowing effects remain (for instance in the atomic hydrogen map), these figures show how some species (e.g. CO, CN) trace denser gas than others (C+, CH). This appears more clearly when plotting these abundances, averaged over density bins, versus total gas density nH (Fig. 9): C+ traces gas uniformly up to nH ≳ 103 cm-3, while CO starts rising up at nH ≳ 250 cm-3, right about where CH flattens out. This break in the slope for CO occurs after the molecular transition ⟨fH2⟩ = 2⟨n(H2)/nH⟩ = 1/2, which is at nH ≃ 100 cm-3. Considering the abundances of C and C+, this means that a significant fraction of the molecular gas (i.e. where hydrogen is mostly in the form of H2) is better traced by C and C+ than by CO. This “dark molecular gas” fraction is the subject of Sect. 5.4. CH, on the other hand, nicely follows H2 (Sheffer et al. 2008). To complete the picture, C and CN have a similar slope throughout the density range, although CN seems to break away slightly at nH ≳ 103 cm-3, to follow CO.

To be more precise on these apparent correlations (CH vs. H2, CS vs. C and CN vs. CO), we show in Fig. 10 abundance ratios for these similarly distributed species in the region where they are all significantly present, that is the cloudlet at (X ≃ −4.7   pc,Y ≃ −24   pc) (see Figs. 7 and 8). We can see that n(CH)/n(H2) and n(CN)/n(CO) have a similar “ringlike” behaviour, rising to maximum values n(CH)/n(H2) ~ 5.5 × 10-8 and n(CN)/n(CO) ~ 10-2 at total gas densities nH ~ 400−500 cm-3, then falling at higher densities, to n(CH)/n(H2) ~ 4 × 10-8 and n(CN)/n(CO) ~ 2 × 10-3, respectively. For n(CS)/n(C), there is a slight loss of azimuthal symmetry, but the overall trend is the same, although maximum values of  ~10-3 are reached at higher densities nH ~ 1000 cm-3 before falling to  ~3 × 10-4 at the peak.

thumbnail Fig. 7

H (top left), H2 (top right), C+ (bottom left) and C (bottom right) abundances for the los models. Contours mark total gas densities 20, 100, 500, 1000, and 2000 cm-3. C+ and C abundance maps are clipped below 2 × 10-4 cm-3 and 10-4 cm-3, respectively. Lines of sight 0, 1, 2, and 3 in the n(H2) map refer to the discussion in the text.

Open with DEXTER

thumbnail Fig. 8

Same as Fig. 7 but for CO (top left), CS (top right), CH (bottom left), and CN (bottom right). Abundance maps are clipped below 2 × 10-5 cm-3 for CO, 2 × 10-10 cm-3 for CS, 2 × 10-8 cm-3 for CH and 2 × 10-9 cm-3 for CN.

Open with DEXTER

The comparison to observational data requires us to work not with densities but with column densities, which are what observers have access to. To his end, we computed column densities for CO, CH and CN on each line of sight parallel to the X or Y axis and plot them versus those of H2 in Fig. 11 to match the observational data plots in Sheffer et al. (2008). We plot data points that correspond to the long lines of sight parallel to X (squares) and to the shorter lines of sight parallel to Y (circles) using a colour scheme to specify the mean gas density ⟨nH⟩ on the line of sight. Fits to observational data derived by Sheffer et al. (2008) are shown as dashed lines, and in the CO panel we also plot actual data points from that same paper.

Our data shows a significant deficit around N(H2) ~ 1020 cm-2 for CO, which may be a general problem with PDR chemistry computations (Sonnentrucker et al. 2007). However, the agreement with the observational fit improves, both in values and in slope, for N(H2) ≳ 2 × 1020 cm-2. The slope found is d [log N(CO)] /d [log N(H2)]  ≃ 5.2, while the observational fit yields 3.07 ± 0.73. On the longer lines of sight, we can see a “loop” structure that needs to be understood. It is easily identified with lines of sight between Y = −24.5   pc and Y = −23   pc (positions 0 to 3 on the H2 and CO maps from Figs. 7 and 8). To be more accurate, following the loop clockwise corresponds to scanning lines of sight parallel to X from 0 to 3, with turnovers at positions 1 and 2. Considering the H2 and CO maps alongside Fig. 11 helps to understand this “loop”. Firstly, lines of sight from 0 to 1 basically intercept just one dense structure with both H2 and CO, so that we have a very similar behaviour to that of the short lines of sight parallel to Y. Secondly, lines of sight from 1 to 2 pass through many less dense CO structures, with a lot of molecular hydrogen in between, leading to the sharp drop in the N(CO) vs. N(H2) relation. At a given N(CO), the difference ΔN(H2) between H2 column densities in the branches 01 and 12 thus represents the “dark gas” (see Sect. 5.4). Finally, lines of sight from 2 to 3 barely have any CO and contain less and less H2 as we approach position 3, where we return to a situation similar to that at position 0. The split between the two branches is definitely related to the mean gas density on the line of sight, the upper branch having ⟨nH⟩ ~ 250 cm-3, the lower one having ⟨nH⟩ ~ 150 cm-3. The overall deficit in CO suggests that we may not shield it sufficiently from the ambient UV field, a hypothesis that we discuss in Sect. 6.

The behaviour of N(CH) with respect to N(H2) (Fig. 11 – middle panel) agrees remarkably well with the observational data fitted by Sheffer et al. (2008). There is a discrepancy for N(H2) ≲ 1020 cm-2, but it is irrelevant, bacause there are no detections in that range, only upper limits.

We typically have a factor 10 deficit for CN (Fig. 11 – bottom panel) compared to observations. It appears that the reaction rate coefficient for the CN + N → C + N2 reaction in the KIDA database2 may be too high (E. Roueff, priv. comm.), which could partly explain that deficit.

Table 4

Typical threshold abundances Xα used in Figs. 7 and 8.

thumbnail Fig. 9

Abundances of H2, C+, C, CO, CH, CS and CN versus total gas density nH. Data points are averaged in the same nH bins as in Fig. 6. The vertical line marks the position of the average molecular transition where  ⟨ fH2 ⟩  = 2⟨n(H2)/nH⟩ = 1/2.

Open with DEXTER

thumbnail Fig. 10

Abundance ratios n(CH)/n(H2) (top), n(CN)/n(CO) (middle), and n(CS)/n(C) (bottom) in the vicinity of the total gas density peak, located at (X ≃ −4.7   pc,Y ≃ −24   pc). Contours mark total gas densities nH of 100, 200, 300, 400, 500, 1000, and 2000 cm-3.

Open with DEXTER

thumbnail Fig. 11

Column densities of CO (top), CH (middle), and CN (bottom) versus column densities of H2, in the los models. Circles correspond to lines of sight parallel to Y and squares to lines of sight parallel to X. Their colours reflect the mean gas density ⟨nH⟩ on these lines of sight. Plus signs on the top (CO) plot stand for observational data points (Sheffer et al. 2008; Crenny & Federman 2004; Pan et al. 2005; Lacour et al. 2005; Rachford et al. 2002, 2009; Snow et al. 2008). The dashed lines are power-law fits from Sheffer et al. (2008). The lines of sight parallel to X marked 0, 1, 2, and 3 on the top panel are the same as in Figs. 7 and 8.

Open with DEXTER

5.3. Comparison of the los and uniform models

To assess the effects of taking into account density fluctuations, as opposed to the assumption of a homogeneous medium, which is usually made when modelling PDRs, we computed the column densities of CO and H2 derived from the uniform models, and plot them in Fig. 12. The behaviour at low N(H2) is very similar to that in the los models (Fig. 11), and the data suffer from the same deficit in CO around N(H2) ~ 1020 cm-2. However, there is a definite difference at higher H2 column densities: the slope of the relation between both column densities is much steeper, d [log N(CO)] /d [log N(H2)]  ≃ 14, than in the los models and in observational data fits. This break occurs later, at N(H2) ~ 5 × 1020 cm-2, and applies to the short lines of sight (parallel to Y and parallel to X between positions 0 and 1) that show essentially one structure in both H2 and CO. However, the maximum column densities reached are slightly lower than in the los models by a factor  ~3, for both CO and H2. This shows the importance of taking into account density fluctuations along the line of sight when modelling PDRs.

For CH, the behaviour is very similar in the uniform and los models, although here also the maximum column densities reached are slightly lower in the uniform models, by a factor  ~2. The observational fit is recovered at somewhat higher H2 column densities (2 × 1020 cm-2 instead of 1020 cm-2), and the scatter of data points is slightly larger.

thumbnail Fig. 12

Same as Fig. 11 but for the uniform models.

Open with DEXTER

Finally, regarding CN, if we ignore the deficit already seen in the los models, we reach the same conclusions: In the uniform models, a break in the slope occurs at higher H2 column densities N(H2) ~ 6 × 1020 cm-2 (instead of N(H2) ~ 4 × 1020 cm-2 in the los models), the slope is definitely steeper in that high-column-density regime, and the maximum column densities reached are lower, here by a factor  ~4.

5.4. Dark molecular gas fraction

To estimate the amount of molecular gas not seen in CO in our simulation – the so-called “dark gas” (Grenier et al. 2005; Planck Collaboration 2011) – we need to compute the CO line emission and compare its spatial distribution with that of H2. To that effect, we focused on the cloudlet at (X ≃ −4.7   pc,Y ≃ −24   pc), where the gas density peak is found, and assumed a very crude cylindrical cloud model by replicating the density maps along the Z axis over a line-of-sight L = 1   pc, which roughly corresponds to the cloudlet’s extent in the (X,Y) plane. This effectively yields column density maps that we can use with the RADEX radiative transfer code (Van der Tak et al. 2007) to obtain emission maps in the CO (J = 1 → 0) rotational transition line at 115.271 GHz and in the [CII] fine structure transition line at 158 μm. To be precise, at each position (X,Y), we treated radiative transfer along Z in a plane-parallel slab geometry. RADEX works with the escape probability formalism (Sobolev 1960), which requires specification of the line width ΔV. We estimate it to be , where σ3D = 1.8   km   s-1 is the total gas velocity dispersion listed in Table 1. Indeed, for any species α with molecular weight μα (μCO = 28 and μC+ = 12), the ratio of thermal to one-dimensional turbulent velocity dispersions reads In the cloudlet under study, T ≲ 100   K, so the above ratio is typically  ≲0.03 for CO and  ≲0.06 for C+. It is thus reasonable to take for all RADEX runs. The code also requires specification of the gas kinetic temperature and the densities of collisional partners (H2 for CO; H2, H and electrons for [CII]), which we obtained from the PDR code outputs. RADEX is thus run on every line of sight parallel to the Z axis, and results are combined into a CO (J = 1 → 0) emission map and a [CII] emission map, both shown in Fig. 13.

thumbnail Fig. 13

Synthetic emission maps in CO(J = 1 → 0) (left) and [CII] at 158 μm (right). The solid contour marks the assumed 0.4 K km s-1 detection threshold for CO, the dashed contour marks the position of the line centre optical depth τCO = 1, and the dotted contour marks the position of the molecular transition fH2 = 1/2. Grey areas are outside of the computational domain.

Open with DEXTER

Between the solid and dotted contours is the “dark molecular gas” region where hydrogen is predominantly in its molecular form but CO emission fails to detect it. We assumed a detection threshold WCO = 0.4   K   km   s-1 consistent with the noise level in e.g. the CO survey of Taurus by Goldsmith et al. (2008). On the other hand, that same gas can definitely be traced in the [CII] line, because it has a typical integrated emission of  ~0.4−0.8   K   km   s-1, while the sensitivity quoted by Velusamy et al. (2010) for the GOTC+ key program is  ~0.1−0.2   K   km   s-1.

thumbnail Fig. 14

Fraction of “dark gas” (solid line) along one-dimensional cuts parallel to the Y axis going through the cloudlet at (X ≃ −4.7   pc,Y ≃ −24   pc). Also shown are the fraction of “dark gas” computed using the definition by Wolfire et al. (2010) (dashed line), and a normalised profile of the total gas column densities NH along the same cuts (dash-dotted line).

Open with DEXTER

The “dark molecular gas” fraction associated with this cloudlet can be estimated by taking one-dimensional cuts parallel to the Y axis going through the CO emission region. Along such a cut, which is parameterised by X, we defined Y0(X) and Y1(X) as the boundaries of the computational domain3 (sharp transition from grey to white in the panels of Fig. 13), and we denote with YCO(X) the region where the integrated emission of CO (J = 1 → 0) exceeds the detection threshold WCO (region enclosed by the solid contour in Fig. 13). We then define the “dark gas” fraction as Figure 14 shows this fraction as a function of the position X of the one-dimensional cut. Obviously, fDG = 1 when the cut does not pass through the CO emission region, and fDG < 1 when some of the H2 is traced by CO. We found that in this cloudlet, at least 20% of H2 is not traced by CO, even at the peak of the gas density. To obtain a mean fraction of dark gas in this cloudlet, we averaged fDG over the range of X coordinates at which CO is seen (i.e. fDG(X) < 1), weighted by the total gas column density NH. This yields , which is somewhat higher than the findings of Velusamy et al. (2010), who identified 53 “transition clouds” with H i and 12CO emission but no 13CO, and found that  ~25% of H2 in these clouds belong to an H2/C+ layer not seen in CO. However, the scatter in observational values is large (Grenier et al. 2005; Abdo et al. 2010), so the small discrepancy is no cause for alarm. Another possible comparison is with Wolfire et al. (2010), who constructed spherical models of molecular clouds to study the dark gas fraction, which they defined in a similar way, except that the boundary of their CO region is specified by the condition of unit optical thickness at the line centre, τCO = 1. In our clump, that isocontour is very close to our own condition ICO = WCO, as can be seen in Fig. 13. Computing the average dark gas fraction with Wolfire et al. (2010)’s condition yields , which is quite close to their results fDG ~ 0.25−0.33. There is a notable difference between their models and ours, however, since our cloudlet has a mass  ~9.5   M inside the CO region, while Wolfire et al. (2010) studied GMCs with masses in the range 105   M to 3 × 106   M. Their impinging UV field was also notably higher (χ = 3–30).

6. Discussion and summary

6.1. Illumination effects

We bypassed the one-dimensionality of the PDR code by combining runs in two orthogonal directions, taking at each grid point the chemical composition corresponding to maximum radiation energy density E. However, this choice is questionable: if a position is shielded from radiation in almost every direction but for one small hole, the illumination at this location resulting from our procedure is too high. As this means forming less molecules, we estimated if it might account for some of the CO deficit seen in Fig. 11. To do so in a simple way, we computed the chemical composition under the opposite assumption, i.e. based on the criterion of minimum local illumination. The result is plotted in Fig. 15, and shows how this indeed helps recovering observed CO column densities for N(H2) ≳ 1020 cm-2, with a consistent scatter. Below N(H2) ≃ 1020 cm-2, a significant CO deficit remains, however.

Obviously, this choice of minimum local illumination is also unphysical, and the reality must lie somewhere in between. A physically better, but more computationally intensive method is being pursued and will be presented in a future paper.

6.2. FGK approximation

Our study used the Federman et al. (1979, FGK) approximation to compute self-shielding. This may underestimate the shielding of CO by molecular hydrogen lines, therefore we performed a few runs of the PDR code using exact radiative transfer (Goicoechea & Le Bourlot 2007). We did this on some lines of sight for which N(H2) ~ 1020 cm-2, to see whether this helps fill the CO deficit in that region. It turns out that the CO column densities so obtained are indeed higher than those found in the FGK approximation, but by a factor  ≲2, which is not sufficient to explain our CO deficit. Because this increased computational time by a factor  ~5–6, we feel that this approach is not to be recommended.

thumbnail Fig. 15

CO column densities versus H2 column densities in the los models when combining data based on a criterion of minimum local illumination. Symbols are the same as in Fig. 11, top panel.

Open with DEXTER

6.3. Steady-state assumption

We took the simulation cube to be a static background, under the assumption that timescales for chemistry and photoprocesses are much shorter than those of the MHD simulation.

From the analysis performed by Le Petit et al. (2006) on uniform density PDRs, it appears that timescales for H2 photodissociation at the edges of a cloud are  ~1000/χ yr, where χ is the FUV radiation strength in units of the Draine (1978) field. In our analysis, we chose χ = 1 so that the corresponding timescale is about 1000 yr.

Estimating timescales for the MHD simulation is more difficult, because what we are interested in is the time over which structures remain coherent, and we do not have access to this information due to the Eulerian nature of the simulation, which makes it impossible to confidently identify structures. For a rough estimate, we may consider the overall crossing time τcross = L/VX ≃ 2.4   Myr, but this does not correspond to the time over which gas is mixed by turbulence at a given scale. For this, we may use the velocity dispersion σ3D listed in Table 1 within the observational clump, whose size is about 2 pc. This yields a dynamical timescale τdyn ≃ 1.1   Myr, which is very similar to the values quoted by Wolfire et al. (2010) to validate the steady-state assumption in their models.

The chemical timescale is that of the formation of molecular hydrogen. As shown by Glover & Mac Low (2007) using numerical simulations of decaying ISM turbulence that include a simplified chemical network, the formation timescale for H2 in turbulent magnetised molecular clouds is τchem ~ 1−2   Myr. It is therefore of the order of the estimated dynamical timescales in our simulation, so that our steady-state assumption seems only marginally valid. However, H2 is formed in dense regions and transported throughout the entire volume via turbulent motions, so we may safely assume steady state, provided we consider a sufficiently late snapshot. Indeed, if H2 starts forming when the converging WNM flows collide near the midplane (τcoll ≃ 1   Myr), and if it is fully formed and transported throughout the entire volume after τchem + τdyn + τcross, this requires taking a snapshot timed at no earlier than 5.4 Myr, which is the case here (t = 7.35   Myr). We conclude that our steady-state assumption is a legitimate one.

6.4. Warm chemistry

We note that chemistry is here driven by UV radiation only, but that there is an important pathway for the formation of many molecular species, which is warm chemistry in turbulence dissipation regions (TDR), studied by Joulain et al. (1998) and Godard et al. (2009). In particular, the CO abundances found in the models by Godard et al. (2009) are higher than in corresponding PDR models, sometimes by almost an order of magnitude. More generally, Godard et al. (2009) argued that observed chemical abundances are on the whole well reproduced if dissipation is caused by ion-neutral friction in sheared structures  ~100   AU thick. A TDR post-processing of our MHD simulation is in preparation, to compare both types of chemistry.

6.5. Summary

We have presented a first analysis of UV-driven chemistry in a simulation of the diffuse ISM, by post-processing it with the Meudon PDR code. Our results show that assuming a uniform density medium when modelling PDRs leads to significant errors: for instance, the maximum CO column densities found with this simplistic assumption are a factor  ~3 lower than those found using actual density fluctuations. The slope of the H2-CO correlation at N(H2) ≳ 5 × 1020 cm-2 is also a factor  ~3 higher than in the more realistic case, and therefore agrees less well with observations. A second result of our study is that in the densest parts of the simulation (nH ≳ 103 cm-3), some 35% to 40% of the molecular gas is “dark”, in the sense that it it not traced by the CO(J = 1 → 0) line, given current sensitivities. It is detectable via the [CII] fine structure transition line at 158 μm, however. As a side result, we found that the simplified cooling used in the MHD simulation by Hennebelle et al. (2008) yields gas temperatures that agree reasonably well with those found using the more detailed processes included in the PDR code.


2

http://kida.obs.u-bordeaux1.fr/

3

The dependence of fDG on the boundaries of the computational domain is necessarily small, because there is little mass at low densities.

4

We used version 1.4.1 of the PDR code, with a fixed H2 formation rate .

5

I.e. density-weighted averages for temperature and pressure.

6

Namely, for the los models, 16 out of 291 along Y and 3 out of 156 along X; for the uniform models, 8 out of 291 along Y and 6 out of 156 along X.

Acknowledgments

The authors acknowledge support from ANR through the SCHISM project (ANR-09-BLAN-231), and computing resources and services from France Grilles and the EGI e-infrastructure. Some kinetic data have been downloaded from the online KIDA (KInetic Database for Astrochemistry, http://kida.obs.u-bordeaux1.fr) database. Colour figures in this paper use the cubehelix colour map by Green (2011).

References

Appendix A: Using the Meudon PDR code with density profiles

Table A.1

Parameters used in the PDR code as input in the .in files.

thumbnail Fig. A.1

Illumination mask computed by comparing at each point (X,Y) the local radiative energy densities EX and EY output by the PDR code along the X and Y directions, respectively. Regions where are marked in black and regions where EX > EY are marked in white. Contour lines of equal total gas density are overlaid at 20, 100, 500, 1000, and 2000 cm-3. Grey areas are outside of the computational domain.

Open with DEXTER

This Appendix is meant as an introduction to using the Meudon PDR code4 with fluctuating density profiles. For more detailed presentations of the code, the reader is referred to Le Bourlot et al. (1999), Le Petit et al. (2006), and Gonzalez Garcia et al. (2008).

The code requires two input files: a .pfl file listing visual extinction AV, temperature T (in K) and total gas density nH (in cm-3) along the line of sight, and a .in file supplying the parameters of the run to perform. These are listed in Table A.1, and some of them require a short comment:

  • modele is the basename chosen for output files. Forthe los models, we used a name of the generic formlos_zZ_yY_xXm-Xp or los_zZ_xX_yYm-Yp, reflecting the positionof the extracted profile in the cube. Note that Z is a constant throughout this paper, as is evident from Fig. 2. For the uniform models, we used names of the generic form uniform_zZ_yY_xXm-Xp or uniform_zZ_xX_yYm-Yp.

  • ifafm is the number of global iterations to use. As Le Petit et al. (2006) pointed out, for diffuse clouds (AV < 0.5) proper convergence may require up to 20 iterations, therefore we selected ifafm=20 for all our models. For information, 415 of the 447 profiles have total AV < 0.5.

  • Avmax is that same total visual extinction through the cloud, which is simply the last AV value in the .pfl file.

  • We set the density densh, temperature tgaz, and pressure presse=densh × tgaz parameters to the average values5 for each profile. They are not used by the code with a density-temperature profile, but they are used, in the reference uniform models, as initial guesses for thermal balance computation.

  • radm and radp specify the strengths χm and χp of the incident radiation field in units of the ISRF on the left and right sides of the profile, respectively. For the runs described in this paper, we used χm = χp = 1.

  • fprofil specifies the .pfl density-temperature profile file. For consistency, we used the same naming scheme as for modele.

  • vturb is the “turbulent velocity dispersion”. It does not include thermal dispersion, therefore we took it to be the standard deviation of the line-of-sight velocity within each structure, noted σX and σY in Tables 2 and 3, respectively.

  • ifisob is a flag specifying whether to use a density profile. For the los models, we therefore set ifisob=1, to enforce the use of a density-temperature .pfl file. However, since thermal balance is solved (ieqth=1), the temperature values in the file are only used as initial guesses. For the uniform models, we set ifisob=0 to use a constant density (specified by densh).

Appendix B: Post-processing of raw outputs

B.1. Resampling

Outputs of the PDR code are FITS data files and XML description files, and we used dedicated scripts to extract specific quantities into plain text files for the subsequent analysis. Among the quantities retrieved are the distance d from the surface of the structure, visual extinction AV, proton column density NH, temperature TPDR, proton density nH, pressure p, ionization fraction xe, and abundances n(α) of 99 chemical species.

Because the PDR code produces its own mesh refinement to better solve for the H/H2 transition, these quantities are sampled irregularly. Consequently, we resampled outputs on the same regular grid as the MHD simulation, using a simple linear interpolation method. This allowed us to build raw maps for all output quantities from PDR code runs along the X and Y directions.

B.2. Missing data

For reasons that are unclear, a few6 of the 894 PDR runs did not complete successfully on the EGEE grid. To supplement the missing data, we interpolated along the perpendicular direction. Consider the ensemble of runs along the X direction: completed runs yield quantities FX(X,Y), and if a run is missing at coordinate Y = Y0, we supplied FX(X,Y0) by linearly interpolating G(Y) = FX(X0,Y) at constant X = X0. This yields a satisfactory completion of the raw data.

B.3. Combination of X and Y runs

We then proceed to the combination of data from runs along the X and Y directions, as described in 3. Figure A.1 shows the “illumination mask” computed by comparing the local radiation energy densities EX and EY output by the PDR code in los models parallel to the X and Y directions, respectively. This mask is then used to build a single data array for each quantity F output by the PDR code at each grid point (X,Y), according to the rule: This helps to reduce the shadowing artefacts due to the one-dimensionality of the PDR code, while ensuring element conservation in each grid cell, and yields the final maps that are analysed and discussed in the main body of the paper. In the discussion (Sect. 6), we also used the inverse choice:

All Tables

Table 1

Properties of the selected observational clump.

Table 2

Properties of the 156 one-dimensional profiles extracted parallel to the X axis.

Table 3

Same as Table 2 but for the 291 one-dimensional profiles extracted parallel to the Y axis.

Table 4

Typical threshold abundances Xα used in Figs. 7 and 8.

Table A.1

Parameters used in the PDR code as input in the .in files.

All Figures

thumbnail Fig. 1

Total gas column density along the X axis of the MHD simulation snapshot used here. The box in the upper left side marks the position of the clump selected for running our analysis.

Open with DEXTER
In the text
thumbnail Fig. 2

Maximum gas density nH along the lines of sight parallel to the X axis within the selected clump. Contours show the total gas column density NH from 3 × 1021 cm-2 to 1.1 × 1022 cm-2 in steps of 1021 cm-2. The “clump” displayed here has a size roughly 1.5   pc    ×    2.5   pc. The two-dimensional slab of gas under study is seen projected as the single-pixel-wide AB strip. The white triangle marks the position of the maximum value of max(nH) in this region, and the white circle that of the maximum of NH. They are separated by 0.65 pc.

Open with DEXTER
In the text
thumbnail Fig. 3

Structure of the gas in the two-dimensional slab under study. This is a close-up view on the region X ≃ 0 in which the incoming flows of WNM collide to form cold structures, and only regions in which are shown. The colour image shows the gas density nH in logarithmic scale, while contours show the gas temperature TMHD at 20 K (higher densities), 30 K and 50 K (lower densities). The two-dimensional projection of the velocity field (VX,VY) is also shown as yellow arrows, the length of which indicate the velocity modulus at that location (at the centre of each arrow). The A and B extremities of the observational one-dimensional strip are indicated for reference, and the dashed line marks an example location for the extracted density profiles on which the PDR code is run (Fig. 5). The grey areas are outside of the domain used for PDR computations.

Open with DEXTER
In the text
thumbnail Fig. 4

Total gas column densities along the X axis within the AB strip under study (top plot) and maximum gas density max(nH) on the same lines of sight (bottom plot). Shown are the column densities over the full 50 pc lines of sight along the X axis (dashed line) and the column densities for the overdense regions shown in Fig. 3 (solid line). Grey areas mark lines of sight for which less than 50% of the mass is in the overdense region. The dash-dotted lines mark the positions of the maxima of NH and max(nH), which are separated by 0.55 pc.

Open with DEXTER
In the text
thumbnail Fig. 5

Example of a density profile used in the PDR code. This is the profile extracted at the location of the dashed line in Fig. 3.

Open with DEXTER
In the text
thumbnail Fig. 6

Ratio of the gas temperatures TPDR/TMHD versus total gas density nH. Grey crosses show all points, while black circles are average values over density bins in logarithmic scale, with error bars standing for  ±1σ. The dashed line corresponds to TPDR/TMHD = 1.

Open with DEXTER
In the text
thumbnail Fig. 7

H (top left), H2 (top right), C+ (bottom left) and C (bottom right) abundances for the los models. Contours mark total gas densities 20, 100, 500, 1000, and 2000 cm-3. C+ and C abundance maps are clipped below 2 × 10-4 cm-3 and 10-4 cm-3, respectively. Lines of sight 0, 1, 2, and 3 in the n(H2) map refer to the discussion in the text.

Open with DEXTER
In the text
thumbnail Fig. 8

Same as Fig. 7 but for CO (top left), CS (top right), CH (bottom left), and CN (bottom right). Abundance maps are clipped below 2 × 10-5 cm-3 for CO, 2 × 10-10 cm-3 for CS, 2 × 10-8 cm-3 for CH and 2 × 10-9 cm-3 for CN.

Open with DEXTER
In the text
thumbnail Fig. 9

Abundances of H2, C+, C, CO, CH, CS and CN versus total gas density nH. Data points are averaged in the same nH bins as in Fig. 6. The vertical line marks the position of the average molecular transition where  ⟨ fH2 ⟩  = 2⟨n(H2)/nH⟩ = 1/2.

Open with DEXTER
In the text
thumbnail Fig. 10

Abundance ratios n(CH)/n(H2) (top), n(CN)/n(CO) (middle), and n(CS)/n(C) (bottom) in the vicinity of the total gas density peak, located at (X ≃ −4.7   pc,Y ≃ −24   pc). Contours mark total gas densities nH of 100, 200, 300, 400, 500, 1000, and 2000 cm-3.

Open with DEXTER
In the text
thumbnail Fig. 11

Column densities of CO (top), CH (middle), and CN (bottom) versus column densities of H2, in the los models. Circles correspond to lines of sight parallel to Y and squares to lines of sight parallel to X. Their colours reflect the mean gas density ⟨nH⟩ on these lines of sight. Plus signs on the top (CO) plot stand for observational data points (Sheffer et al. 2008; Crenny & Federman 2004; Pan et al. 2005; Lacour et al. 2005; Rachford et al. 2002, 2009; Snow et al. 2008). The dashed lines are power-law fits from Sheffer et al. (2008). The lines of sight parallel to X marked 0, 1, 2, and 3 on the top panel are the same as in Figs. 7 and 8.

Open with DEXTER
In the text
thumbnail Fig. 12

Same as Fig. 11 but for the uniform models.

Open with DEXTER
In the text
thumbnail Fig. 13

Synthetic emission maps in CO(J = 1 → 0) (left) and [CII] at 158 μm (right). The solid contour marks the assumed 0.4 K km s-1 detection threshold for CO, the dashed contour marks the position of the line centre optical depth τCO = 1, and the dotted contour marks the position of the molecular transition fH2 = 1/2. Grey areas are outside of the computational domain.

Open with DEXTER
In the text
thumbnail Fig. 14

Fraction of “dark gas” (solid line) along one-dimensional cuts parallel to the Y axis going through the cloudlet at (X ≃ −4.7   pc,Y ≃ −24   pc). Also shown are the fraction of “dark gas” computed using the definition by Wolfire et al. (2010) (dashed line), and a normalised profile of the total gas column densities NH along the same cuts (dash-dotted line).

Open with DEXTER
In the text
thumbnail Fig. 15

CO column densities versus H2 column densities in the los models when combining data based on a criterion of minimum local illumination. Symbols are the same as in Fig. 11, top panel.

Open with DEXTER
In the text
thumbnail Fig. A.1

Illumination mask computed by comparing at each point (X,Y) the local radiative energy densities EX and EY output by the PDR code along the X and Y directions, respectively. Regions where are marked in black and regions where EX > EY are marked in white. Contour lines of equal total gas density are overlaid at 20, 100, 500, 1000, and 2000 cm-3. Grey areas are outside of the computational domain.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.