H2 distribution during formation of multiphase molecular clouds
Valeska Valdivia, Patrick Hennebelle, Maryvonne Gérin, Pierre Lesaffre
AAstronomy & Astrophysics manuscript no. Valdivia_27325 c (cid:13)
ESO 2018November 5, 2018 H distribution during the formation of multiphase molecularclouds Valeska Valdivia , , , Patrick Hennebelle , , Maryvonne Gérin , , and Pierre Lesa ff re , Laboratoire de radioastronomie, LERMA, Observatoire de Paris, École Normale Supérieure, PSL Research University, CNRS,UMR 8112, 75005, Paris, Francee-mail: [email protected] Laboratoire AIM, Paris-Saclay, CEA / IRFU / SAp - CNRS - Université Paris Diderot, 91191 Gif-sur-Yvette Cedex, Francee-mail: [email protected] Sorbonne Universités, UPMC Univ. Paris 06, UMR 8112, LERMA, Paris, France, F-75005Received Month dd, yyyy; accepted Month dd, yyyy
ABSTRACT
Context. H is the simplest and the most abundant molecule in the interstellar medium (ISM), and its formation precedes the formationof other molecules. Aims.
Understanding the dynamical influence of the environment and the interplay between the thermal processes related to theformation and destruction of H and the structure of the cloud is mandatory to understand correctly the observations of H . Methods.
We performed high-resolution magnetohydrodynamical colliding-flow simulations with the adaptive mesh refinement codeRAMSES in which the physics of H has been included. We compared the simulation results with various observations of the H molecule, including the column densities of excited rotational levels. Results.
As a result of a combination of thermal pressure, ram pressure, and gravity, the clouds produced at the converging point ofHI streams are highly inhomogeneous. H molecules quickly form in relatively dense clumps and spread into the di ff use interclumpgas. This in particular leads to the existence of significant abundances of H in the di ff use and warm gas that lies in between clumps.Simulations and observations show similar trends, especially for the HI-to-H transition (H fraction vs total hydrogen column den-sity). Moreover, the abundances of excited rotational levels, calculated at equilibrium in the simulations, turn out to be very similarto the observed abundances inferred from FUSE results. This is a direct consequence of the presence of the H enriched di ff use andwarm gas. Conclusions.
Our simulations, which self-consistently form molecular clouds out of the di ff use atomic gas, show that H rapidlyforms in the dense clumps and, due to the complex structure of molecular clouds, quickly spreads at lower densities. Consequently, asignificant fraction of warm H exists in the low-density gas. This warm H leads to column densities of excited rotational levels closeto the observed ones and probably reveals the complex intermix between the warm and cold gas in molecular clouds. This suggeststhat the two-phase structure of molecular clouds is an essential ingredient for fully understanding molecular hydrogen in these objects. Key words.
H2– molecular clouds – ISM – column density – star formation
1. Introduction
It is well known that stars form in dense and self-gravitatingmolecular clouds and that the star formation rate per unit area( Σ SFR ) is on average relatively well correlated with the H sur-face density ( Σ H ) (Leroy et al. 2013; Lada et al. 2012; Wong &Blitz 2002; Bigiel et al. 2008).Although H is the simplest and most abundant moleculein the interstellar medium (ISM), it is very hard to observe di-rectly. Because of its homonuclear nature, H lacks a permamentdipole and only weak quadrupolar transitions are allowed. More-over, excitation energies are very high and require very hightemperatures or strong ultraviolet (UV) fields to excite its ro-tational levels (Kennicutt & Evans 2012). H can be observed inemission by infrared (IR) rovibrational transitions (Burton et al.1992; Santangelo et al. 2014; Habart et al. 2011), or in absorp-tion at far-ultraviolet (far-UV) wavelengths from the Lyman andWerner bands. These bands were first observed by Spitzer &Jenkins (1975) using the Copernicus satellite. The
Far Ultravi-olet Spectroscopic Explorer (
FUSE ) (Moos et al. 2000) o ff erednew perspectives on the study of H in the ISM through its sen- sitivity, which is 10 times higher than Copernicus in the far-UVpart of the spectrum, providing measurements of H column den-sities (of the total column density of H and for several rotationallevels J ) along translucent lines of sight (Wakker 2006; She ff eret al. 2008).While it seems to be clear that molecular gas is well corre-lated with star formation at di ff erent scales, it is still an openquestion how the gas becomes molecular. Since understand-ing the atomic-to-molecular hydrogen (HI-to-H ) transition isof the highest importance for understanding the star-formingprocess, numerous models have been developed (e.g. Krumholzet al. 2008; Sternberg et al. 2014) and seem able to reproducemany of the observed constraints. While extremely useful, thesemodels leave aside the dynamical aspects of H formation. Inparticular, the question of the relatively long timescale that isneeded to form the H molecules has for many years been dif-ficult to reconcile with short-lived and quickly forming molecu-lar clouds. The typical timescale for H formation is of the or-der of 10 yr / n (Hollenbach et al. 1971), which would lead toages older than 10 Myr for molecular clouds of mean densi- Article number, page 1 of 16 a r X i v : . [ a s t r o - ph . GA ] D ec & A proofs: manuscript no. Valdivia_27325 ties of the order of 10-100 cm − (Blitz & Shu 1980). However,Glover & Mac Low (2007a) and Glover & Mac Low (2007b),who have simulated the formation of molecular hydrogen in su-personic clouds, showed that H can form relatively quickly inthe dense clumps induced by the shocks, leading to a signifi-cantly shorter timescale (typically of about 3 Myr). Recently,Micic et al. (2012) explored the influence of the turbulent forc-ing and showed that it has a significant influence on the timescaleof H formation. This is because compressible forcing leads todenser clumps than solenoidal forcing, for which the motions areless compressible. This clearly shows that dynamics is playingan important role in the process of HI-to-H transition, at least atthe scale of a molecular cloud.In the Milky Way most of the molecular hydrogen is in theform of low-temperature gas, but di ff erent observations havebrought to light the existence of large amounts of fairly warmH in the ISM (Valentijn & van der Werf 1999; Verstraete et al.1999). Gry et al. (2002) and Lacour et al. (2005) have shownthat the H excitation resulting from UV pumping and from H formation cannot account for the observed population of H ex-cited levels ( J > is likely explained by the pres-ence of a warm and turbulent layer associated with the molecu-lar cloud. But in such warm gas, which is characterised by verylow densities ( n H ≈ −
10 cm − ) and very high temperatures( T ≈ − K), H formation on grain surfaces becomesnegligible, which accordingly leads to an apparent contradiction.Recently, Godard et al. (2009) have investigated the possibilitythat warm H could form during intermittent high-energy dissi-pation events such as vortices and shocks. They showed in par-ticular that under plausible assumptions for the distributions ofthese events, the observed abundance of excited H moleculescould be reproduced. Another possibility to explain the abun-dance of these excited H molecules is that they could be formedin dense gas and then be transported in more di ff use and warmermedium. This possibility has been investigated by Lesa ff re et al.(2007) using a 1D model and a prescription to take into ac-count the turbulent di ff usion between the phases. In particular,they found that the abundance of H molecules at low densityand high temperature increases with the turbulent di ff usion e ffi -ciency.While di ff erent in nature, previous works (Glover & MacLow 2007a,b; Lesa ff re et al. 2007; Godard et al. 2009) there-fore consistently found that the formation of H is significantlyinfluenced by the dynamics of the flow in which it forms. Al-though the exact mechanism that leads to the formation ofmolecular clouds is still under investigation, it seems unavoid-able to consider converging streams of di ff use gas (e.g. Hen-nebelle & Falgarone 2012; Dobbs et al. 2014). What exactlytriggers these flows is not fully elucidated yet but is most likelydue a combination of gravity, large-scale turbulence, and large-scale shell expansion. Another important issue is the thermody-namical state of the gas. The di ff use gas out of which molec-ular clouds form is most likely a mixture of phases made ofwarm and cold atomic hydrogen (HI) and even possibly some-what di ff use molecular gas. Such a multiphase medium presentsstrong temperature variations (typically from 8000 K to lowerthan 50 K). This causes the dynamics of the flow to be signif-icantly di ff erent from isothermal or nearly isothermal flow (seeAudit & Hennebelle 2010, for a comparison between barotropicand two-phase flows). Moreover, the two-phase nature of thedi ff use atomic medium has been found to persist in molecularclouds (Hennebelle et al. 2008; Heitsch et al. 2008a; Banerjeeet al. 2009; Vázquez-Semadeni et al. 2010; Inoue & Inutsuka2012). This in particular produces a medium in which super- sonic clumps of cold gas are embedded in a di ff use phase ofsubsonic warm gas. Consequently, the question about the conse-quences of the multiphase nature of a molecular cloud are for theformation of H . Clearly, various processes may occur. First ofall, the dense clumps are continuously forming and accreting outof the di ff use gas because molecular clouds continue to accrete.Second, some of the H formed at high density is likely to bespread in low-density high-temperature gas because of the phaseexchanges.We here present numerical simulations using a simple modelfor H formation within a dynamically evolving turbulent molec-ular cloud that is formed through colliding streams of atomicgas. Such flows, sometimes called colliding flows, have beenwidely used to study the formation of molecular clouds becausethey represent a good compromise between the need to assem-ble the di ff use material from the large scales and the need to ac-curately describe the small scales within the molecular clouds(Hennebelle & Pérault 1999; Ballesteros-Paredes et al. 1999;Heitsch et al. 2006; Vázquez-Semadeni et al. 2006; Clark et al.2012). The content of this paper is as follows: in Sect. 2 wepresent the governing equations and the physical processes, suchas H formation. In Sect. 3 we describe the numerical setup thatwe use to perform the simulations. In Sect. 4 we present the re-sults of our simulations, first describing the general structure ofthe cloud and then focussing on the H distribution. We performa few complementary calculations that aim at better understand-ing the physical mechanisms at play in the simulations. Then inSect. 5 we compare our results with various observations thathave quantified the abundance of H and find very reasonableagreements. Finally, we summarise and discuss the implicationsof our work in Sect. 6.
2. Physical processes
We consider the usual compressive magnetohydrodynamicalequations that govern the behaviour of the gas. These equationswritten in their conservative form are ∂ρ∂ t + ∇ · ( ρ v ) = , (1) ∂ρ v ∂ t + ∇ · ( ρ vv − BB ) + ∇ P = − ρ ∇ φ, (2) ∂ E ∂ t + ∇ · [( E + P ) v − B ( Bv )] = − ρ L , (3) ∂ B ∂ t + ∇ · ( vB − Bv ) = , (4) ∇ φ = π G ρ, (5)where ρ , v , B , P , E , and φ are the mass density, velocity field,magnetic field, total energy, and the gravitational potential of thegas, respectively, and L is the net loss function that describes gascooling and heating. The fluid equations (Eqs. 1-5) are complemented by an equationthat describes the formation of the H molecules, ∂ n H ∂ t + ∇ · ( n H v ) = k form n ( n − n H ) − k ph n H , (6) Article number, page 2 of 16aldivia et al.: H distribution during formation of molecular clouds Fig. 1.
Column density maps. Total hydrogen column density map (left), H column density map (centre), and molecular fraction (right). The gasenters the box from both sides ( x -axis), and the conditions for the other sides are periodic. The gas is compressed along the x -direction and thecloud forms toward the middle of the box. y [ p c ] l o g n [ c m − ]
20 Myr y [ p c ] l o g n H [ c m − ]
20 Myr y [ p c ] l o g χ s h i e l d
20 Myr
Fig. 2.
Slices cut through the mid-plane. The panel on the left shows the total number density of hydrogen nucleons, the panel in the centreshows the number density of molecular hydrogen (on the same scale), and the panel on the right shows the mean shielding factor for the H photodissociation rate coe ffi cient, calculated as (cid:104) exp( − τ d , ) f shield (cid:105) . where n is the total hydrogen density, n H represents the densityof H , and k form and k ph represent the formation and destructionrates of H . formation When two hydrogen atoms encounter each other in the gasphase, they cannot radiate the excess of energy because of thelack of an electric dipole moment, and therefore H forma-tion in the gas phase is negligible. It is widely accepted to-day that H is formed through grain catalysis (Hollenbach &Salpeter 1971; Gould & Salpeter 1963). Hydrogen atoms canadsorb onto the grain surfaces and encounter another hydro-gen atom to form a H molecule through two mechanisms:the Langmuir-Hinshelwood mechanism, where atoms are ph-ysisorbed (e ffi cient in shielded gas), and the Eley-Rideal mech-anism, or chemisorption (e ffi cient on warm grains) (Le Bourlotet al. 2012; Bron et al. 2014). The detailed physical descriptionof these two mechanisms is complex and the numerical treat-ment is computationally expensive. For this reason we adopted asimpler description for H formation on grain surfaces, given bythe following mean formation rate: k form , = × − cm s − . (7)This rate was first derived by Jura (1974), using Copernicus ob-servations, and was confirmed by Gry et al. (2002) using
FUSE observations. The formation rate depends on the local gas tem-perature and on the adsorption properties of the grain, therefore this value is corrected for by the dependence on temperature andby a sticking coe ffi cient: k form = k form , (cid:114) T
100 K × S ( T ) . (8) S ( T ) is the empirical expression for the sticking coe ffi cient ofhydrogen atoms on the grain surface as described in Le Bourlotet al. (2012), S ( T ) = + (cid:16) TT (cid:17) β , (9)where we use the same fitting values as Bron et al. (2014),namely T =
464 K, and β = . destruction The main mechanism that destroys the H molecule is photodis-sociation by absorption of UV photons in the 912 - 1100 Å range(Lyman and Werner transitions). In a UV field of strength G ,H is photodissociated in optically thin gas at a rate (Draine &Bertoldi 1996) of k ph , = . × − G s − , (10)but the H gas can protect itself against photodissociation by twoshielding e ff ects. The first shielding e ff ect is the continuous dustabsorption, while the second e ff ect is the line absorption due Article number, page 3 of 16 & A proofs: manuscript no. Valdivia_27325 to other H molecules, called self-shielding. Draine & Bertoldi(1996) showed that the photodissociation rate can be written as k ph = e − τ d , f shield ( N H ) k ph , . (11)The first term is the e ff ect of the dust. Here τ d , = σ d , N tot is the optical depth along a line of sight ( σ d , = × − cm − is the e ff ective attenuation cross section for dust grains at λ = N tot is the total column density of hydrogen). Thesecond term is the self-shielding factor, and we use the sameapproximation, f shield = . + x / b ) + . + x ) / exp( − . × − (1 + x ) / ) , (12)where x = N H / × cm − , b = b / cm s − , where b is theDoppler-broadening parameter, which is typically 2 km s − , butit can reach values as high as 10 km s − (Shull et al. 2000). Weassume that the turbulent contribution dominates, then we use b = − (see Appendix B for the influence of this choice onthe shielding coe ffi cient). For the heating and cooling of the gas, we used the same treat-ment as we did previously (Valdivia & Hennebelle 2014) (seealso Audit & Hennebelle 2005), which we describe in this sec-tion and call the standard heating and cooling, to which we haveadded the thermal feedback from H , described in Sect. 2.3.1.The dominant heating process in the gas is due to the ultravi-olet flux from the interstellar radiation field (ISRF) through thephotoelectric e ff ect on grains (Bakes & Tielens 1994; Wolfireet al. 1995), where we use the e ff ective UV field strength, cal-culated using our attenuation factor χ for the UV field. We alsoinclude the heating by cosmic rays (Goldsmith 2001), which isimportant in well shielded regions.The primary coolant at low temperature is the P / → P / [CII] fine-structure transition at 158 µ m (Launay & Roue ff µ m and146 µ m line emission (Flower et al. 1986; Tielens 2005; Wolfireet al. 2003). For the Lyman α emission, which becomes dom-inant at high temperatures, we use the classical expression ofSpitzer (1978). We also include the cooling by electron recom-bination onto dust grains using the prescription of Wolfire et al.(1995, 2003) and Bakes & Tielens (1994). We note that strictlyspeaking, this cooling function does not include molecular cool-ing in spite of the high densities reached in the simulation. How-ever, it leads to temperatures that are entirely reasonable even athigh densities and therefore appears to be su ffi ciently accurate,in particular because we are not explicitly solving for the forma-tion of molecules in addition to H (see Levrier et al. 2012, for aquantitative estimate). To the atomic cooling, described before, we added the molecularcooling by H lines and the heating by H formation and destruc-tion.During the H formation process, about 4 . internalenergy (rotational and vibrational excitation), and into the grainheating is not well known, and the derived fraction of this en-ergy that actually heats the gas varies from one author to another [cm − ] l o g V / V t o t Volume-weighted PDF t =5 . t =10 . t =15 . t =20 . Fig. 3.
Density PDF evolution. Colour lines show the density PDF at t = (Le Bourlot et al. 2012; Glover & Mac Low 2007a). We consideran equipartition of the energy, which means that one-third of theenergy released goes into heating the gas, Γ form = . × − k form n H n erg s − cm − . (13)H photodissociation provides an additional heating source. Fol-lowing Black & Dalgarno (1977) and Glover & Mac Low(2007a), we assumed 0 . Γ ph = . × − k ph n H erg s − cm − . (14)H contributes to the cooling of the gas through line emis-sion and can become strong for the di ff use medium under certainconditions (Glover & Clark 2014; Gnedin & Kravtsov 2011). AsH is a homonuclear molecule, only quadrupolar transitions areallowed, and thus, the ortho and para states can be treated inde-pendently. In our simulations we adopted the cooling functionfrom Le Bourlot et al. (1999). This is a function of temperature,total density, n (HI) / n (H ) relative abundance, and of the ortho-to-para-H ratio (OPR). It considers transitions between the first51 rovibrational energy levels. As the cooling function is quiteinsensitive to the OPR, we fixed it to 3 for simplicity.
3. Numerical setup and initial conditions
As already discussed, colliding flows are a practical ansatz togather matter to form a molecular cloud (Inoue & Inutsuka2012; Audit & Hennebelle 2005; Vázquez-Semadeni et al. 2007;Heitsch et al. 2005, 2006, 2008b; Micic et al. 2013) to mimiclarge-scale converging flows, as in Galactic spiral arms or insuper-bubble collisions.
The setup is very similar to the one that Valdivia & Hennebelle(2014) used. The size of the simulation box is L =
50 pc andit is initially uniformly filled with atomic gas, with initial num-ber density n tot = − and initial temperature T = V =
15 km s − , aligned with the x -axis, whichis modulated by a function of amplitude (cid:15) = .
5, producing
Article number, page 4 of 16aldivia et al.: H distribution during formation of molecular clouds − ]76543210 l o g V / V t o t Volume-weighted PDF, t =20 . V tot =1 . , x ∈ [0 , V tot =0 . , x ∈ [17 , V tot =0 . , x ∈ [22 , V tot =0 . , x ∈ [27 , Fig. 4.
Volume filling factor as a function of density in di ff erent regions. V tot is the volume fraction of each region. a slightly turbulent profile, as in Audit & Hennebelle (2005).We use periodic boundary conditions for the remaining sides.The gas is initially uniformly magnetised, with a moderate mag-netic field ( ∼ . µ G) aligned with the inflowing gas. This con-figuration avoids boundary issues and facilitates the buildingof the cloud. Introducing an angle between the magnetic andthe velocity fields would be more realistic, but this very sim-plified framework requires relatively small angles (Hennebelle& Pérault 2000; Körtgen & Banerjee 2015; Inoue & Inutsuka2012).
Molecular clouds are embedded in the interstellar radiation field(ISRF), which is assumed to be isotropic and constant through-out the space (Habing 1968). The whole simulation is embed-ded in an external and isotropic UV field that is assumed tobe monochromatic, of strength G = . χ , which is calculated at eachpoint of the simulation box for each timestep using our tree-based method, which is fully described in Valdivia & Hennebelle(2014). Hence, the strength of the local UV field, after dustshielding, can be written G = χ G , where χ = (cid:104) e − . V (cid:105) = (cid:104) e − . × − N tot (cid:105) , is the mean shielding factor due to dust. Thismean value is calculated using a fixed number of directions. Inour previous work (Valdivia & Hennebelle 2014) we have shownthat the number of directions is not crucial for the dynamical evo-lution of the cloud. For simplicity and numerical e ffi ciency, wetherefore used 12 directions ( M = N = ffi ciency.We adapted the code to also compute the self-shielding(stated by Eq. 12). More precisely, along each direction we com-puted not only the total gas density, but also the H column den-sity from which we computed the shielding coe ffi cient, f shield ,and the extinction due to dust, τ d , . Since the UV flux is as-sumed to be isotropic, the final photodissociation rate is obtained [cm − ] l o g P t h [ g c m − s − ] Thermal pressure vs density, t =20Myr [cm − ] l o g P r a m [ g c m − s − ] Dynamic pressure vs density, t =20Myr [cm − ] l o g P m a g [ g c m − s − ] Magnetic pressure vs density, t =20Myr Fig. 5. Di ff erent pressures in the ISM vs density. The top panel showsthe thermal pressure, P th = nkT , the middle panel the ram pressure, P ram = ρ V , and the bottom panel the magnetic pressure, P ma g = B / π in the simulation box. The red line shows the mean value per densitybin. The colour scale indicates the number of points. by taking the mean value over all directions. We define χ shield as χ shield = k ph k ph , = (cid:104) e − τ d , f shield ( N H ) (cid:105) . (15) To perform our simulations we employed the adaptive mesh re-finement code RAMSES (Teyssier 2002; Fromang et al. 2006).RAMSES solves the MHD equations using the HLLD Riemannsolvers. It preserves the nullity of the divergence of the magneticfield thanks to the use of a staggered mesh. To solve Eq. (6),
Article number, page 5 of 16 & A proofs: manuscript no. Valdivia_27325
Fig. 6.
Statistics of clumps. The top panel shows the mass spectrum, thesecond panel their internal velocity dispersion as a function of size, thethird panel the kinetic β = P ram / P ma g , and the bottom panel displaysthe H fraction as a function of mass. we used operator splitting, solving first for the advection (whichis identical to the conservation equation) and then subcyclingto solve for the right-hand side. In AMR codes the refinementis done on a cell-by-cell basis. In our simulations, the refine-ment criterion is the density. When a cell reaches a given den-sity threshold, it is split into eight smaller cells, each one hav-ing the same mass and volume. The process is repeated recur-sively until the maximum resolution is reached. In our fiducialrun, we allowed two AMR levels (cid:96) min = (cid:96) max =
10, lead-ing to an equivalent numerical resolution of 1024 cells and an [cm − ] l o g χ UV screening vs density, t= Myr log n [cm − ] l o g χ s h i e l d Total shielding for H2 vs density, t= Myr Fig. 7.
UV screening factor vs density (top panel) and total shieldingcoe ffi cient for H vs density (bottom panel) at t =
20 Myr. The red lineshows the mean value per density bin, and the colour scale indicates thenumber of points. e ff ective spatial resolution of about 0 .
05 pc. For the first refine-ment level (from level (cid:96) = (cid:96) =
9) the density thresholdis n thresh =
50 cm − and for the next refinement, the densitythreshold is n thresh =
100 cm − . To investigate numerical conver-gence, a particularly crucial issue when chemistry is considered,we also performed a high-resolution run for which the resolu-tion was doubled (that is to say, using levels from 9 to 11 withthe same refinement criterion).We ran our simulations for about 20 Myr. We note thatthese calculations are somewhat demanding because of the shorttimesteps induced by high temperatures, therefore it was not pos-sible to run the high-resolution simulations up to this point. Thisrun reaches up to 15 Myr. To check further for convergence, wealso performed low-resolution runs (128 and 256 , see Ap-pendix) that we compared with our highest resolution simula-tions. We also performed complementary runs in which we mod-ified the physics of H formation to better understand how H isformed. More precisely, in these runs we suppressed the H for-mation above a certain density threshold,
4. Results
We now present the results of our calculations. We start by dis-cussing the general properties of the clouds because it is essen-
Article number, page 6 of 16aldivia et al.: H distribution during formation of molecular clouds t = t =
10 Myr t =
15 Myr t =
20 Myr . f ( H l o g m i n H f o r m [ M fl ] n ∈ [0 . , .
0] cm − simulationat equilibrium f ( H -4-3-2-1 l o g m i n H f o r m [ M fl ] n ∈ [0 . , .
0] cm − f ( H -5-4-3-2-1 l o g m i n H f o r m [ M fl ] n ∈ [0 . , .
0] cm − f ( H -3-2-1 l o g m i n H f o r m [ M fl ] n ∈ [0 . , .
0] cm − f ( H -3-2-10 l o g m i n H f o r m [ M fl ] n ∈ [1 . , .
0] cm − f ( H -4-3-2-101 l o g m i n H f o r m [ M fl ] n ∈ [1 . , .
0] cm − f ( H l o g m i n H f o r m [ M fl ] n ∈ [1 . , .
0] cm − f ( H l o g m i n H f o r m [ M fl ] n ∈ [1 . , .
0] cm − f ( H -1012 l o g m i n H f o r m [ M fl ] n ∈ [10 . , .
0] cm − f ( H l o g m i n H f o r m [ M fl ] n ∈ [10 . , .
0] cm − f ( H l o g m i n H f o r m [ M fl ] n ∈ [10 . , .
0] cm − f ( H l o g m i n H f o r m [ M fl ] n ∈ [10 . , .
0] cm − f ( H -2-1012 l o g m i n H f o r m [ M fl ] n ∈ [100 . , .
0] cm − f ( H -3-2-10123 l o g m i n H f o r m [ M fl ] n ∈ [100 . , .
0] cm − f ( H -10123 l o g m i n H f o r m [ M fl ] n ∈ [100 . , .
0] cm − f ( H -3-2-101234 l o g m i n H f o r m [ M fl ] n ∈ [100 . , .
0] cm − f ( H -101234 l o g m i n H f o r m [ M fl ] n ∈ [1000 . , .
0] cm − f ( H -101234 l o g m i n H f o r m [ M fl ] n ∈ [1000 . , .
0] cm − f ( H -101234 l o g m i n H f o r m [ M fl ] n ∈ [1000 . , .
0] cm − f ( H -101234 l o g m i n H f o r m [ M fl ] n ∈ [1000 . , .
0] cm − Fig. 10.
Mass in H form per density bin. Comparison between the simulation (blue) and the calculation at equilib-rium (dashed) for di ff erent timesteps. From left to right: t = , , , and 20 Myr. From top to bottom: density bins n ∈ (0 . , , (1 , , (10 , , (100 , , (1000 , − ]. Note the change of vertical scale. tial to understand in which context the hydrogen molecules form.Then we discuss the abundance and distribution of H itself. Fi-nally, we perform various complementary calculations aiming atbetter understanding in which conditions H forms. As many colliding flow calculations have been performed duringthe past decade, we do not attempt here to display all the steps that the flow experiences (see references above for accurate de-scriptions), so we quickly summarise the main steps. First, theWNM flows collide, which triggers the transition for a fractionof the gas into moderately dense gas , that is to say, densitiesof the order of 100-1000 cm − as a consequence of the thermalinstability and ram pressure. Second, when enough gas has been If the incoming flow is fast enough, a collision is not even neces-sary, as it is the case for example in the study presented by Koyama &Inutsuka (2002). Article number, page 7 of 16 & A proofs: manuscript no. Valdivia_27325 [Myr] l o g f ( H ) Mass fraction in H2 form vs t
Fig. 8.
Evolution of the molecular fraction in the simulation box as afunction of time. accumulated, gravity becomes important and triggers infall, firstat small scales and then at larger ones. This leads to a continuousincrease of the dense gas fraction.The left panels of Figs. 1 and 2 show the column densityalong the z -direction and a density cut in the x y plane at time 20Myr, respectively. As we show below, this corresponds to a phasewhere gravity already plays a significant role. The density cutshows that the flow is very fragmented. The dense gas is organ-ised into dense clumps that are embedded in di ff use and warmgas. The column density typically spans 2-3 orders of magnitudefrom a few 10 cm − to at least 10 cm − in some very compactregions, the dense cores in which gravity triggers strong infall.The column density map is also clumpy, but less obviously struc-tured than the density cut. In many respects the density distribution is an essential cloudproperty that reveals the dynamical state of the cloud andstrongly influences the formation of H .Figure 3 shows the density probability distribution function(PDF) at time 5, 10, 15, and 20 Myr. The peak at low density( n (cid:39) − ) is simply a consequence of the initial conditionsand boundary conditions (which inject WNM from the bound-aries), the rest, at higher densities, represents the cold phase,whose formation has been triggered by the thermal instability.Because of the supersonic turbulence that develops in the coldgas, the density PDF is broad and presents a lognormal shape(Kritsuk et al. 2007; Federrath et al. 2008). At later times, morecold gas accumulates and the PDF broadens. The high-densitytail tends to become less steep. At intermediate densities ( ρ ∼ -10 cm − ) the shape is compatible with a power law withan exponent between -1 and -2. This is typical of what has beenfound in super-sonic turbulent simulations that include gravity(Kritsuk et al. 2011). At very high densities ( > − cm − ),the density field flattens. This occurs because we did not use sinkparticles, and the very dense gas therefore piles up and accumu-lates within a few grid cells. This feature is thus a numericalartefact and represents a limit of the simulations.Figure 4 shows the density PDF for di ff erent computationalbox regions at a time of 20 Myr. The four lines show the resultsfor four di ff erent box regions (as shown in the label). The blackline shows the whole box, the red line is limited to the densest [Myr] l o g f ( H ) Molecular fraction vs t n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ] [Myr] l o g f ( H ) Molecular fraction vs t T ∈ ( . , . ) [K]T ∈ ( . , . ) [K]T ∈ ( . , . ) [K]T ∈ ( . , . ) [K]T ∈ ( . , . ) [K]T ∈ ( . , . ) [K] Fig. 9.
Molecular fraction evolution. In the top panel we show the evolu-tion per density bin (the density increases from purple to red), in the bot-tom panel the evolution per temperature bin (the temperature increasesfrom purple to salmon). regions of the computational box. The volume filling factor isclearly dominated by the warm and di ff use gas. Selecting gasof densities n < − , we find that it typically occupies 70%for the region located between x =
27 and x =
37 (85% for thewhole computational box). The dense gas, even in the densest re-gion, occupies only a tiny fraction, for example we find that thegas denser than 100 cm − has a filling factor of about 3% (respec-tively 1.5% for the whole box). The remaining 26% (13%) arefilled with gas of densities between 3 and 100 cm − . We there-fore conclude that the interclump medium that occupies most ofthe volume, is itself made of two components: a warm gas thatis similar to the standard WNM, but can be slightly denser, anda cold gas that is similar to the CNM, but contains, as we showbelow, a significant fraction of H . They typically occupy about70% and 25% percent of the volume, respectively. Another important diagnostic for characterizing the dynamics ofa medium are the di ff erent pressures. Figure 5 shows the ther-mal, P th , dynamical, P ram , and magnetic, P ma g , pressures equalto P th = nkT , P ram = ρ V and P ma g = B / π , respectively.The dynamical pressure clearly dominates the thermal pres-sure by typically one to two orders of magnitude in the dense Article number, page 8 of 16aldivia et al.: H distribution during formation of molecular clouds [Myr] l o g f ( H ) Molecular fraction vs t n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ] [Myr] l o g f ( H ) Molecular fraction vs t n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ] [Myr] l o g f ( H ) Molecular fraction vs t n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ]n ∈ ( . , . ) [cm − ] Fig. 11.
Evolution of the molecular fraction, where the formation of H has been suppressed for gas denser than a fixed threshold. The top panelcorresponds to the standard case, with no suppression of H formation.The middle panel shows the evolution for a density threshold of n = − ], while the panel at the bottom corresponds to a densitythreshold of n =
100 [cm − ] gas ( n >
10 cm − ), while they are more similar in the di ff usegas. The magnetic pressure lies in between the two, with val-ues a few times higher than that of the thermal pressure. P ram and P ma g increase with density, and at densities of the order of n (cid:39) cm − , they are about one order of magnitude higherthan their mean values at n (cid:39)
10 cm − . This means that whilethe low-density gas that fills the volume provides some confin-ing pressure, it has a limited influence on the clumps. − ]333231302928272625 l o g Γ , Λ [ e r g s s − H − ] Heating and cooling vs density Γ std Γ H2 f Γ H2 d Λ std Λ H2 Fig. 12.
Various cooling and heating contributions as a function of den-sity. As can be seen, H does not have a dominant influence except atdensities of about 4 cm − where H cooling (green line) becomes com-parable to the standard ISM cooling (blue line) and at high densitieswhere the heating by H formation (red line) is significant (but withoutmodifying the temperature substantially). Altogether, these results indicate that the molecular cloudproduced in this calculation can be described by density fluc-tuations or clumps, induced by both ram pressure and gravity.The clumps occupy a tiny fraction of the volume, which is filledby a mixture of warm di ff use and dense HI gas (occupying afraction >
70% and >
20% of the volume, respectively). This low-density material feeds the clumps, which grow in mass. This pic-ture agrees well with the measurement performed by Williamset al. (1995), who found that the interclump medium typicallyhas a density of a few particles per cm − and a high velocitydispersion of several km s − . Because the cloud is organised into clumps, we further quanti-fied the simulations by providing some statistics. To identify theclumps, we used a density threshold of 1000 cm − and a friends-of-friends algorithm. This structure is important for the chem-istry evolution because very significant density and temperaturegradients arise at the edge of clumps . Figure 6 shows the massspectrum, the velocity dispersion, and the kinetic β parameter(that is to say, P ram / P ma g ) of the clumps. The mass spectrum formasses above a few 0.1 M (cid:12) presents a power law with an in-dex of about − .
7, which is very similar to what has been shownin related works (e.g. Audit & Hennebelle 2010; Heitsch et al.2008b). The velocity dispersion within the clumps as a functionof their size is about σ (cid:39) − ( L / β parameter is typically ofthe order of, or slightly below 1, showing that the magnetic fieldplays an important role within the clumps. We now discuss the H abundance. The middle and right panelsof Fig. 1 show the H column density and molecular abundance, f (H ) = n (H ) / n . As expected, the high column density re-gions are dominated by the H molecules, and values of f (H ) Article number, page 9 of 16 & A proofs: manuscript no. Valdivia_27325 close to 1 are obtained there. The values of f (H ) are obviouslysignificantly lower in the outer parts of the cloud. The middle and right panels of Fig. 2 show a cut of n (H ) andthe value of the shielding parameter, χ shield = k ph / k ph , . This lat-ter steeply decreases when entering the cloud, where it takes val-ues of the order of 10 − . A comparison between the density cutwith k ph / k ph , reveals that there are regions of di ff use materialin which the shielding parameter is low. This is because theseregions are surrounded by dense gas in which H is abundant,which provides an e ffi cient self-shielding. As a consequence,there are low-density regions that contain a relatively high abun-dance of H (see left and middle panels of Fig. 2).A more quantitative statement can be drawn from Fig. 7,which displays the distribution of the dust shielding, χ (toppanel), and of the total shielding χ shield as a function of den-sity (bottom panel). As expected, most low-density cells are as-sociated with χ values that are close to 1, that is to say, theyare weakly shielded. There is a fraction of them, however, thatpresents values of χ shield as low as 10 − . At higher densities, themean value of χ decreases with a steep drop between n (cid:39) and 10 cm − , which corresponds to the point where the dustsignificantly absorbs the UV external field. For the total shield-ing, a steep drop is observed at about 10 cm − . There is a broadscatter, however, which increases from n = n = cm − for the dust shielding and tends to decrease for the total shield-ing χ shield . This indicates that the mean value of χ shield is notsu ffi cient information to quantify the abundance of H expectedat a given density. This is a clear consequence of the complexcloud structure. Since χ shield plays a key role in the formationof H , this complex distribution certainly makes molecular hy-drogen formation in a multiphase turbulent cloud a complicatedproblem. Figure 8 shows f (H ), the mean molecular fraction within thewhole cloud, as a function of time. During the first 5 Myr, f (H )exponentially increases from nearly 0 to about (cid:39) .
1. After thisphase, f (H ) continues to increase in a more steady almost lin-ear way. Near 15 Myr, about 40% of the gas is in H molecules.Since, as discussed before, the cloud is rather inhomogeneousin density and temperature, we also plot the time evolution perbin of density and temperature. The corresponding curves aredisplayed in Fig. 9. As expected, the total distribution closelyfollows the values of the densest bins, that is to say, it corre-sponds to densities higher than 100-1000 cm − . For these den-sities, the timescale for H formation is thus of a few ( (cid:39) − / n yr. The timescale we observe in the simulation therefore agreeswell with this scaling. Altogether, this agrees well with the re-sults of Glover & Mac Low (2007a,b), who found that H couldbe formed quickly in molecular clouds because it preferentiallyforms into clumps that are significantly denser than the rest ofthe cloud. In the present case, the mean density of our cloud attime 5-10 Myr would be about 10-100 cm − (depending on ex-actly which area is selected), and this would lead to a timescalefor H formation of the order of 10 to 100 Myr.The relatively short amount of time ( (cid:39)
10 Myr) after whichthere is a significant amount of H in low density gas ( n (cid:39) − − ) is somewhat more surprising. The expected timescale for H formation in this density range would be about 10 Myr. Fig-ure 9 shows that the behaviour of f (H ) is similar, although notidentical, to the behaviour at higher densities. There is a fast in-crease followed by a slowly increasing phase. The slope changesat about 4 Myr however, and this is di ff erent from the higher den-sity evolution. The steady evolution starts at about 8 Myr insteadof 5-6 Myr for the denser regions. This suggests that the pres-ence of H at low density is triggered by formation of H in thedense gas. This can occur in various non-exclusive ways. First ofall, some dense gas can expand back and mix with di ff use gas.This may occur either through a sonic expansion, through evapo-ration, or through numerical di ff usion. Second of all, some of thedi ff use gas may be surrounded by dense gas in which f (H ) islarge and therefore has a low χ shield . Below we attempt to analysethese mechanisms in more detail.Finally, we note that in the bottom panel of Fig. 9 thereis a small (1-10%) but nevertheless non-negligible fraction ofH within gas of relatively high temperature, that is, T > is the first molecule involved in producing othermolecules and since some of them, such as CH + , are formedthrough reactions with high activation barriers (of the order of4300 K for CH + , see Agúndez et al. 2010), the presence of H could have consequences for the production of these species, asproposed in Lesa ff re et al. (2007). In the same way, as discussedbelow, this warm H contains molecules in excited high J rota-tional levels, which can therefore contribute to the gas cooling. To further quantify the distribution of H in the cloud, we dis-play histograms of f (H ) per density bin (Fig. 10). We also drawthe distribution of f (H ) obtained at equilibrium in the samepanel. Knowing the density, n , the temperature, T and the shield-ing factor, χ shield , we solve Eq. (6) at equilibrium. We note thatthis distribution is not fully self-consistent since the value of χ shield should in principle be recalculated to be consistent withthis equilibrium distribution (this would imply performing sev-eral iterations). In practice, the goal here is simply to comparewith the time-dependent distribution to gain insight into the H production mechanisms, and it is therefore easier to have exactlythe same formation and destruction rates.In particular, the equilibrium distribution is expected to bedi ff erent from the time-dependent distribution for at least twomain reasons. First, since the H formation timescale is some-what long, if the time-dependent f (H ) lies below the equilib-rium timescale, it will indicate that the H abundance has beenlimited by the time to form the molecule. Second, if the time-dependent f (H ) is above the equilibrium value, it will indicatethat f (H ) has increased because of an enrichment coming fromdenser gas. As discussed above, this could operate through theexpansion or the evaporation of cold clumps, a process that wecall turbulent di ff usion, or through numerical di ff usion. The lat-ter is quantified in the appendix by performing convergence stud-ies that suggest that it remains limited. Another possibility isthat the H excess has been produced during a phase where thefluid parcel was more shielded by the surrounding material andtherefore the UV field was lower. The two e ff ects probably actsimultaneously and are di ffi cult to separate.At a time of 5 Myr (left panels), all density bins show anexcess of the equilibrium distribution with respect to the time-dependent distribution. This is clearly due to the long timescaleneeded to form H . The di ff erences between the equilibrium andtime-dependent distributions become eventually less important.For example, the two distributions are obviously closer at a time Article number, page 10 of 16aldivia et al.: H distribution during formation of molecular clouds of 20 Myr (right panels) than at a time of 5 Myr. However, theyare not identical. Since all distributions at a time of 15 Myr and20 Myr are similar, the persistence of the di ff erences betweenthe equilibrium and time-dependent distributions indicates thatthis is most likely the result of a stationary situation. It mostlikely reflects the accretion process, that is to say, HI gas (possi-bly mixed with a fraction of H ) is continuously being accretedwithin denser clumps and therefore the dense gas does not be-come fully molecular. This interpretation is also consistent withthe fact that the e ff ect is more pronounced for the fourth densitybin (line 4) than for the fifth and highest one (line 5). Consideringnow the lowest density bin (first line), we find the reverse e ff ectat a time of 10 and 15 Myr: the time-dependent distribution dom-inates the equilibrium distribution. This is most likely an e ff ectof the turbulent di ff usion. Because of the low mass contained inthis low-density bin, this remains a limited e ff ect, however.The second density bin (line 2), which corresponds to n be-tween 1 and 10 cm − , is slightly puzzling. No significant di ff er-ences are seen between the two distributions, which is surprising.One possibility is that various e ff ects compensate for each other.That is to say, the time delay to form H (clearly visible for thethree denser density bins) may be compensated for by an en-richment from the denser gas. This might possibly occur for thesecond density bin at time 20 Myr (fourth column, second line)where overall a small excess of time-dependent H is found, butfor f (H ) > . As a complement to the time evolution of f (H ), it is also worth-while to investigate the spatial fluctuations at a given time. Inthis respect, the clumps constitute natural entities for studyingthe spatial variations. The bottom panel of Fig. 6 shows f (H )within the clumps. For the most massive clumps, the value of f (H ) is close to 1 and the dispersion remains weak. In contrast,the dispersion becomes very high for the less massive clumps,and the value of f (H ) can in some circumstances be rather low.This shows that clump histories, that is to say, their ages andthe local UV flux in which they grow, have a major influence on f (H ). formation To better quantify the interdependence of the di ff erent densitybins, we conducted three complementary low-resolution runs forwhich the formation of H was suppressed ( k form was set to zero)above a given density threshold. Figure 11 displays the result. Inthe top panel, no threshold is applied. In the middle and bottompanels a threshold of 1000 cm − and 100 cm − , respectively, isapplied. For the threshold 1000 cm − , the values of f (H ) arenot significantly changed except for the highest density range(for n > − ). On the other hand, with a threshold 100cm − , the values of f (H ) decrease by a factor of about 3 for alldensity bins. This clearly shows that most of the H molecules inthe low-density gas form at a density of a few 100 cm − . Whilethe shielding provided by molecules in gas denser than this valuecould contribute to enhance the more di ff use gas, the filling fac-tor of the dense gas is too small to strongly a ff ect the di ff use gas,which has a much larger filling factor (see Fig. 4). Thereforewe conclude that the H abundance within the low-density gas( <
100 cm − ) is very likely a consequence of turbulent mixingand gas exchange between di ff use and dense gas. ( N HI +2 N H ) [cm − ] l o g f ( H ) Copernicus SurveyFUSE Halo SurveyFUSE Rachford+02FUSE Rachford+09 Av [ mag ] Fig. 13.
Molecular fraction as a function of total hydrogen column den-sity. Comparison with the
Copernicus survey (Savage et al. 1977) and
FUSE (Gillmon et al. 2006; Rachford et al. 2002, 2009).
To quantify the influence of the H molecule on thermodynam-ics, we present the contributions of the various heating and cool-ing terms to the thermal balance as a function of density. Fig-ure 12 shows the standard ISM cooling (blue curve) and thecooling by H (green curve). The latter clearly remains belowthe former throughout, except at a density of about 3-5 cm − ,where they become similar. At densities above 3 × cm − the"standard" heating curve (magenta) is dominated by the heat-ing caused by cosmic rays and reaches the intermediate valuefor shielded gas proposed by Goldsmith (2001). While the heat-ing due to H destruction remains negligible at any density (or-ange curve), the heating due to its formation (red curve) is thedominant heating mechanism at high densities. This does not af-fect the temperature very strongly, however. It remains equal toabout 10 K. Moreover, at higher densities (above ∼ cm − )other processes that are not included in our simulations can takeover, such as collisions with dust grains and the cooling by othermolecular lines, such as CO.Finally, as can be seen for densities between 1 to ∼
50 cm − ,the cooling dominates the heating by a factor of a few. This in-dicates that there is another source of heating equal to the dif-ference, which is due to the mechanical energy dissipation. Thislatter therefore appears to have a contribution similar to the UVheating. This explains why warm gas is actually able to survivewithin molecular clouds. In the same way that density is muchhigher than in the rest of the ISM, kinetic energy is also higherand provides a significant heating (e.g. Hennebelle & Inutsuka2006).
5. Comparison to observations
We now compare our results with various observations that fallinto two categories. The first observations have attempted tomeasure the molecular fraction, that is, the total column den-sity of H with respect to the total column density. The secondobservations have measured the excited levels of H , thereforepresumably tracing high-temperature gas. These two sets of ob-servations are therefore complementary and very informative. Article number, page 11 of 16 & A proofs: manuscript no. Valdivia_27325
The H column density has been estimated along several linesof sight mainly (though not exclusively) with the Copernicus satellite (Savage et al. 1977) and the
FUSE observatory (Gill-mon et al. 2006; Rachford et al. 2002, 2009). These observationsprobe lines of sight spanning a wide range of column densitiesand therefore constitute a good test for our simulations, althougha di ff erence to keep in mind is that the lines of sight extractedfrom our simulations are all taken from a 50 pc region and there-fore are more homogeneous than the observed lines of sight.Figure 13 shows the molecular fraction f (H ) as a functionof total column density for all lines of sight of the simulations(taken along the z -direction) and the available lines of sight re-ported in Gillmon et al. (2006), and Rachford et al. (2002, 2009).The simulation results and the observations agree rather well,many observational data directly fall into the same regions assimulated data. In particular, the two regimes, that is, the ver-tical transition branch at column densities of between 10 and3 × cm − as well as the higher column density region, areclearly seen both in observations and in the simulation.This said, there are also data points that are not reproducedby any of the lines of sight from the simulations. This is particu-larly the case for column densities higher than 3 × cm − and for the Copernicus survey at column densities of around (cid:39) cm − . The most likely explanation is that the UV fluxis di ff erent from the standard value we assumed in our study.We stress in particular that the measurements were made in ab-sorption toward massive stars. Because these are strong emittersof UV radiation, it may not be too surprising that our measure-ments lead to higher values of f (H ). A more quantitative es-timate should entail a detailed modelling of every line of sight,including the UV flux in the regions of interest and specific cloudparameters, such as column densities. This is beyond the scopeof the present paper. We now compare our results with the observations of rotationallevels of the H molecule, which require high temperatures to beexcited.As investigated by previous authors (Lacour et al. 2005; Go-dard et al. 2009), these data cannot be explained by pure UV ex-citation. While Lacour et al. (2005) concluded that to explain theobservations, a warm layer associated with the cold gas shouldbe present, Godard et al. (2009) performed detailed modellingthat entailed shocks or vortices. The general idea is that in suchdissipative small -scale structures the temperature reaches highvalues because of the intense mechanical heating. Interestingly,Godard et al. (2009) obtained a nice agreement between theirmodel and the data. We stress that the small dissipative scalesthat are determinant in these models cannot be described inour simulations, which have a limited resolution of the order of (cid:39) . , it is worthwhile to investigate whetherthe inferred populations of excited H can be reproduced quan-titatively. We calculated the population of the first six rotational levels ( J )of H as in Flower et al. (1986), based on the approach of Elitzur& Watson (1978). This calculation was made in post-processingfor all grid cells. It needs the total hydrogen number density, the E J [K]121416182022 l o g N J / g J [ c m − ] Excited levels of H2, O/P= 0.7 l o g N t o t [ c m − ] E J [K]121416182022 l o g N J / g J [ c m − ] Excited levels of H2, O/P= OPR(T) l o g N t o t [ c m − ] Fig. 14. H population distribution along di ff erent lines of sight at atime t =
20 Myr, using an ortho-to-para ratio of 0 . H number density, the He number density (we assumed 10% ofthe total H), the temperature, and the ortho-to-para ratio (OPR).It assumes thermal equilibrium and includes spontaneous de-excitation and collisional excitation by HI, by He, and by otherH molecules.The OPR is not very well known in the ISM (e.g. Le Bour-lot 2000) and may need a specific time-dependent treatment,which is beyond the scope of this paper. Even though the OPRof newly formed H molecules is thought to take values close to3 (Takahashi 2001), observations di ff er. In high-density regions( n (H ) ∼ cm − ) observations seem to favour values as lowas 0 .
25 (Neufeld et al. 2006; Pagani et al. 2011; Dislaire et al.2012), while in translucent clouds intermediate values close to0 . , weassumed an OPR = .
7, suitable for translucent clouds such asthose produced in our simulations. We also considered for com-parison an OPR given by thermal equilibrium that reaches thestatistical weight of 3 at high temperatures ( T (cid:38)
300 K).
Article number, page 12 of 16aldivia et al.: H distribution during formation of molecular clouds Fig. 15. H level population distribution for individual cells in the sim-ulation calculated at 20 Myr: using an ortho-to-para ratio OPR = . After we obtained the populations of the H rotational levelsand column density, we integrated along several lines of sight inthe z -direction. Figure 14 shows the distribution of the column densities corre-sponding to the various rotational levels in the simulation. Thetop panel shows results for an OPR equal to 0 . , while the bottompanel shows results for the thermal equilibrium assumption. Thecolour coding shows the total column density of the correspond-ing line of sight. The symbols correspond to the observationaldata of Wakker (2006) (cross), Rachford et al. (2002) (circle),Gry et al. (2002) (square), and Lacour et al. (2005) (diamond).The simulation results are divided into two groups of points,very low column densities coming from the WNM that is lo-cated outside the molecular cloud (blue points) and higher col-umn densities (green and yellow points) that come from the gasbelonging to the molecular cloud. For the J = moleculesare in one of these two levels. Their respective abundance sig-nificantly depends on the assumption used for the OPR. Whilethere is about a factor 10 between the first two rotational lev-els for N J /g J when the OPR is equal to 0.7, it is almost a factorthousand when the OPR is assumed to be at thermal equilibrium. E J [K]0246810121416 b = ( b t h + b t u r b ) / [ k m s − ] b parameter, t =20 [Myr] log N tot [cm − ]0 500 1000 1500 2000 2500 E J [K]02468101214 b = ( b t h + b t u r b ) / [ k m s − ] b parameter, t =20 [Myr] log N tot [cm − ]0 500 1000 1500 2000 2500 E J [K]02468101214 b = ( b t h + b t u r b ) / [ k m s − ] b parameter, t =20 [Myr] log N tot [cm − ] Fig. 16. b -parameter along several lines of sight for excited levels of H and comparison with the data from Lacour et al. (2005). The top panelshows lines of sight along the x -axis, the middle panel those along the y -axis, and the bottom panel the lines along the z -axis. The situation is di ff erent for the higher levels. The highestcolumn densities of the excited rotational levels do not corre-spond to the highest total hydrogen column densities, and thereis generally no obvious correlation between the two. This is be-cause the high J levels come from the warm gas (with temper-atures of between a few hundred and a few thousand Kelvin), Article number, page 13 of 16 & A proofs: manuscript no. Valdivia_27325 which itself has a low column density and is largely independentof the cold gas.The comparison with the observational data is very enlight-ening. The agreement with the simulation results is very goodin general. The observational data points have column densitiesthat are very similar to the simulation data. The best agreementis found for OPR = . .
7. This might be in conflict with some of thedata of Gry et al. (2002) and Lacour et al. (2005) for J = .To help understand the origin of the excited rotational lev-els in the simulation, Fig. 15 displays the mean density of theH excited levels for five density bins and the complete distri-bution for the two OPR used in this work. Figure 15 shows thatthe choice of the OPR mainly a ff ects the dense, and thus coldergas, which is the main contributor to the populations of the twofirst levels ( J = J = ff ects the ratio between the populations of these two levels. Atlower densities in warmer gas, the OPR becomes closer to 3,but the populations do not di ff er much from those calculated us-ing OPR = .
7. For the excited level the highest contributions arefound in gas of total density between 1 and 10 cm − . This meansthat very di ff use gas with n < − and moderately dense gas(with n between 10 and 100 cm − ) contribute in roughly equalproportions to the excited levels. We recall that the peak of thecooling contribution of H lies precisely in this density domain(see Fig. 12). The conclusion we can draw so far is that the population of theexcited levels of H ( J ≥
2) is dominated by the warm H thatis interspersed between the cold and dense molecular clumps.This "layer" of warm gas is directly associated with the molec-ular cloud and agrees well with the proposition made by Lacouret al. (2005). These authors made another interesting observa-tion. They found that the width of the lines associated with theexcited levels increases with J . To investigate whether our sim-ulation exhibits the same trends, we calculated the mean ve-locity dispersion along several (125 along each axis) lines ofsight for each excited level. Figure 16 shows the mean Doppler-broadening parameter, b , linked to the velocity dispersion σ by b = / σ , and where σ = (cid:82) ρ (cid:16) ( v i − (cid:104) v i (cid:105) ) + C s (cid:17) dx (cid:82) ρ dx , (16)where C s is the local sound speed, v i is the i component of thelocal velocity, and (cid:104) v i (cid:105) is the mean velocity along the i -direction.Figure 16 displays the results. Since the colliding flow configu-ration is highly non-isotropic, we calculated these velocity dis-persions along the x -axis (top panel), the y -axis (middle panel),and the z -axis (bottom panel). While a more accurate treatmentwould consist of simulating the line profiles and then applyingthe same algorithm as the authors used, this simple estimate isalready illustrative. From σ , we can infer the width of the lineand the b parameter. The trends in the simulation and in the observations are sim-ilar. Higher J levels tends to be associated with larger b . Theagreement is quantiatively satisfying for the y − and z -directions.Similar to the observations, the simulated lines of sight presenta slightly higher velocity dispersion for higher J . The agreementis poor for the lines of sight along the x -direction because theypresent a dispersion that is too high with respect to the observa-tions. This is most likely an artefact of the colliding flow config-uration.The trend of larger b for higher J stems for the fact thathigher levels need warmer gas to be excited. Therefore the fluidelements, which are enriched in high rotational levels, havehigher temperatures and therefore higher velocity dispersions(since both are usually correlated).Altogether, these results suggest that the excited rotationallevel abundances reveal the complex structure of molecularclouds that is two-phase in nature and entails warm gas deeplyinterspersed with the cold gas due to the mixing induced by tur-bulence.
6. Conclusions
We have performed high-resolution magneto-hydrodynamicalsimulations to describe the formation of molecular clouds out ofdi ff use atomic hydrogen streams. We particularly focused on theformation of the H molecule itself using a tree-based approachto evaluate UV shielding (see Valdivia & Hennebelle 2014).In accordance with previous works (Glover & Mac Low2007a), we found that H is able to form much faster than sim-ple estimates, based on cloud mean density, would predict. Thereason for this is that because the clouds are supersonic and havea two-phase structure, H is produced in clumps that are muchdenser than the clouds on average.As a result of a combination of phase exchanges and highUV screening deep inside the multiphase molecular clouds (nu-merical convergence tests suggest that numerical di ff usion is notresponsible of this process), a significant fraction of H devel-ops even in the low-density and warm interclump medium. Thiswarm H contributes to the thermal balance of the gas, and inthe range of density 3 −
10 cm − its cooling rate is similar to thestandard cooling rate of the ISM.Detailed comparisons with Copernicus and
FUSE observa-tions showed a good agreement overall. In particular, the fractionof H varies with the total gas column density in a very similarway, showing a steep increase between column densities 10 and 3 × cm − and a slow increase at higher column den-sities. There is a trend for the high column density regions topresent H fractions significantly below the values inferred fromthe simulations, however. This is a possible consequence of theconstant UV flux assumed in this work. Interestingly, the col-umn densities of the excited rotational levels obtained at thermalequilibrium reproduce the observations fairly well. This is a di-rect consequence of the presence of H molecules in the warminterclump medium and suggests that the H populations in ex-cited levels reveal the two-phase structure of molecular clouds. Acknowledgements.
We thank the anonymous referee for the critical reading andvaluable comments. We thank J. Le Bourlot and B. Godard for the insightful dis-cussions.V.V. acknowledges support from a CNRS-CONICYT scholarship. This researchhas been partially funded by CONICYT and CNRS, according to the December11, 2007 agreement.P.H. acknowledges the financial support of the Agence National pour laRecherche through the COSMIS project. This research has received fundingfrom the European Research Council under the European Community’s SeventhFramework Program (FP7 / Article number, page 14 of 16aldivia et al.: H distribution during formation of molecular clouds M.G. and P.L. thank the French Program "Physique Chimie du Milieu Interstel-laire" (PCMI)This work was granted access to the HPC resources of MesoPSL financed by theRegion Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the AgenceNationale pour la Recherche
References
Agúndez, M., Goicoechea, J. R., Cernicharo, J., Faure, A., & Roue ff , E. 2010,ApJ, 713, 662Audit, E. & Hennebelle, P. 2005, A&A, 433, 1Audit, E. & Hennebelle, P. 2010, A&A, 511, A76Bakes, E. L. O. & Tielens, A. G. G. M. 1994, ApJ, 427, 822Ballesteros-Paredes, J., Vázquez-Semadeni, E., & Scalo, J. 1999, ApJ, 515, 286Banerjee, R., Vázquez-Semadeni, E., Hennebelle, P., & Klessen, R. S. 2009,MNRAS, 398, 1082Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846Black, J. H. & Dalgarno, A. 1977, ApJS, 34, 405Blitz, L. & Shu, F. H. 1980, ApJ, 238, 148Bron, E., Le Bourlot, J., & Le Petit, F. 2014, A&A, 569, A100Burton, M. G., Hollenbach, D. J., & Tielens, A. G. G. 1992, ApJ, 399, 563Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bonnell, I. A. 2012, MNRAS,424, 2599Dislaire, V., Hily-Blant, P., Faure, A., et al. 2012, A&A, 537, A20Dobbs, C. L., Krumholz, M. R., Ballesteros-Paredes, J., et al. 2014, Protostarsand Planets VI, 3Draine, B. T. & Bertoldi, F. 1996, ApJ, 468, 269Elitzur, M. & Watson, W. D. 1978, A&A, 70, 443Federrath, C. 2013, MNRAS, 436, 1245Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79Flower, D. R., Pineau-Des-Forets, G., & Hartquist, T. W. 1986, MNRAS, 218,729Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371Gillmon, K., Shull, J. M., Tumlinson, J., & Danforth, C. 2006, ApJ, 636, 891Glover, S. C. O. & Clark, P. C. 2014, MNRAS, 437, 9Glover, S. C. O. & Mac Low, M.-M. 2007a, ApJS, 169, 239Glover, S. C. O. & Mac Low, M.-M. 2007b, ApJ, 659, 1317Gnedin, N. Y. & Kravtsov, A. V. 2011, ApJ, 728, 88Godard, B., Falgarone, E., & Pineau Des Forêts, G. 2009, A&A, 495, 847Goldsmith, P. F. 2001, ApJ, 557, 736Gould, R. J. & Salpeter, E. E. 1963, ApJ, 138, 393Gry, C., Boulanger, F., Nehmé, C., et al. 2002, A&A, 391, 675Habart, E., Abergel, A., Boulanger, F., et al. 2011, A&A, 527, A122Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421Hayes, M. A. & Nussbaumer, H. 1984, A&A, 134, 193Heitsch, F., Burkert, A., Hartmann, L. W., Slyz, A. D., & Devriendt, J. E. G.2005, ApJ, 633, L113Heitsch, F., Hartmann, L. W., & Burkert, A. 2008a, ApJ, 683, 786Heitsch, F., Hartmann, L. W., Slyz, A. D., Devriendt, J. E. G., & Burkert, A.2008b, ApJ, 674, 316Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A.2006, ApJ, 648, 1052Hennebelle, P., Banerjee, R., Vázquez-Semadeni, E., Klessen, R. S., & Audit, E.2008, A&A, 486, L43Hennebelle, P. & Falgarone, E. 2012, A&A Rev., 20, 55Hennebelle, P. & Inutsuka, S.-i. 2006, ApJ, 647, 404Hennebelle, P. & Pérault, M. 1999, A&A, 351, 309Hennebelle, P. & Pérault, M. 2000, A&A, 359, 1124Hollenbach, D. & Salpeter, E. E. 1971, ApJ, 163, 155Hollenbach, D. J., Werner, M. W., & Salpeter, E. E. 1971, ApJ, 163, 165Ingalls, J. G., Bania, T. M., Boulanger, F., et al. 2011, ApJ, 743, 174Inoue, T. & Inutsuka, S.-i. 2012, ApJ, 759, 35Jura, M. 1974, ApJ, 191, 375Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531Körtgen, B. & Banerjee, R. 2015, ArXiv e-printsKoyama, H. & Inutsuka, S.-i. 2002, ApJ, 564, L97Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2008, ApJ, 689, 865Lacour, S., Ziskin, V., Hébrard, G., et al. 2005, ApJ, 627, 251Lada, C. J., Forbrich, J., Lombardi, M., & Alves, J. F. 2012, ApJ, 745, 190Larson, R. B. 1981, MNRAS, 194, 809Launay, J.-M. & Roue ff , E. 1977, Journal of Physics B Atomic MolecularPhysics, 10, 879Le Bourlot, J. 2000, A&A, 360, 656Le Bourlot, J., Le Petit, F., Pinto, C., Roue ff , E., & Roy, F. 2012, A&A, 541, A76Le Bourlot, J., Pineau des Forêts, G., & Flower, D. R. 1999, MNRAS, 305, 802 Leroy, A. K., Walter, F., Sandstrom, K., et al. 2013, AJ, 146, 19Lesa ff re, P., Gerin, M., & Hennebelle, P. 2007, A&A, 469, 949Levrier, F., Le Petit, F., Hennebelle, P., et al. 2012, A&A, 544, A22Micic, M., Glover, S. C. O., Banerjee, R., & Klessen, R. S. 2013, MNRAS, 432,626Micic, M., Glover, S. C. O., Federrath, C., & Klessen, R. S. 2012, MNRAS, 421,2531Moos, H. W., Cash, W. C., Cowie, L. L., et al. 2000, ApJ, 538, L1Myers, A. T., McKee, C. F., & Li, P. S. 2015, MNRAS, 453, 2747Nehmé, C., Le Bourlot, J., Boulanger, F., Pineau Des Forêts, G., & Gry, C. 2008,A&A, 483, 485Neufeld, D. A., Melnick, G. J., Sonnentrucker, P., et al. 2006, ApJ, 649, 816Pagani, L., Roue ff , E., & Lesa ff re, P. 2011, ApJ, 739, L35Rachford, B. L., Snow, T. P., Destree, J. D., et al. 2009, ApJS, 180, 125Rachford, B. L., Snow, T. P., Tumlinson, J., et al. 2002, ApJ, 577, 221Santangelo, G., Antoniucci, S., Nisini, B., et al. 2014, A&A, 569, L8Savage, B. D., Bohlin, R. C., Drake, J. F., & Budich, W. 1977, ApJ, 216, 291She ff er, Y., Rogers, M., Federman, S. R., et al. 2008, ApJ, 687, 1075Shull, J. M., Tumlinson, J., Jenkins, E. B., et al. 2000, ApJ, 538, L73Spitzer, L. 1978, Physical processes in the interstellar mediumSpitzer, Jr., L. & Jenkins, E. B. 1975, ARA&A, 13, 133Sternberg, A., Le Petit, F., Roue ff , E., & Le Bourlot, J. 2014, ApJ, 790, 10Takahashi, J. 2001, ApJ, 561, 254Teyssier, R. 2002, A&A, 385, 337Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar MediumValdivia, V. & Hennebelle, P. 2014, A&A, 571, A46Valentijn, E. A. & van der Werf, P. P. 1999, ApJ, 522, L29Vázquez-Semadeni, E., Colín, P., Gómez, G. C., Ballesteros-Paredes, J., & Wat-son, A. W. 2010, ApJ, 715, 1302Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870Vázquez-Semadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006,ApJ, 643, 245Verstraete, L., Falgarone, E., Pineau des Forets, G., Flower, D., & Puget, J. L.1999, in ESA Special Publication, Vol. 427, The Universe as Seen by ISO,ed. P. Cox & M. Kessler, 779Wakker, B. P. 2006, ApJS, 163, 282Williams, J. P., Blitz, L., & Stark, A. A. 1995, ApJ, 451, 252Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes,E. L. O. 1995, ApJ, 443, 152Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ,587, 278Wong, T. & Blitz, L. 2002, ApJ, 569, 157 Appendix A: Numerical convergence: resolutionstudy log n [cm − ] l o g M [ M o ] Mass-weighted PDF, t =15 Myr Fig. A.1.
Mass-weighted density PDF for a series of numerical resolu-tions. The discrepancy between low- and high-resolution runs is verysignificant. The di ff erence between the resolution of the standard andhigh resolution runs is low, except for the highest density bins, whichcontain little mass. Numerical resolution is a crucial issue particularly in thecontext of chemistry. Here we compare runs of various resolu-tions.
Article number, page 15 of 16 & A proofs: manuscript no. Valdivia_27325 t =
10 Myr t =
15 Myr l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − l o g m i n H f o r m [ M fl ] n ∈ [ . , . ] cm − Fig. A.2.
Mass in H per density bin. Comparison betweenthe standard resolution simulations (levels 8-10, in blue)and the high-resolution simulations (levels 9-11, dashed) ata time of 10 and 15 Myr. From top to bottom: density bins n ∈ (0 . , , (1 , , (10 , , (100 , , (1000 , − ].The di ff erences between the two remain fairly limited, showing thatnumerical convergence is not a great problem. Figure A.1 shows the mass-weighted density PDF at a timeof 15 Myr for four di ff erent resolutions, 128 , 256 , 1024 , and2048 (for the two last and higher resolutions this represents thehighest e ff ective resolution since AMR was used, as describedabove). The density PDF of low-resolution runs is quite di ff erentfrom the higher ones, which shows the importance of the numer-
14 16 18 20 22 24log N H2 l o g χ s h i e l d χ shield dependence on the b parameter b = 1 km/sb = 2 km/sb = 5 km/sb = 10 km/s Fig. B.1.
Total shielding coe ffi cient for H , χ shield , as a function of totalH column density for several values of the b Doppler parameter. Thedashed lines show the expected value in absence of shielding by dust. ical resolution. The two highest resolution runs present similarbut not identical PDFs. This suggests that for the highest resolu-tion runs, numerical convergence for the density PDF is nearlyreached, although strictly speaking, it would be necessary to per-form runs with an even higher resolution. This conclusion agreeswell with the results of Federrath (2013), who found that the den-sity PDF converges for a resolution inbetween 1024 and 2048 for isothermal self-gravitating simulations.Figure A.2 shows f (H ) distribution in five density bins attimes of 10 and 15 Myr for the standard and high-resolution runs(e ff ective 1024 and 2048 resolution). Overall, the two simula-tions agree well. The di ff erences are similar to the di ff erencesfound in the density PDF. In particular, this suggests that numer-ical di ff usion is not primarily responsible of the fraction of warmH that we observe at low density (as discussed in Sect. 4.2). Appendix B: Influence of the b Doppler parameter
Figure B.1 shows the dependence of the total shielding coe ffi -cient for H (shielding by H and by dust) on the b Dopplerparameter as a function of total H column density. The totalcolumn density was approximated as N tot ≈ × N H , which setsa lower limit to the influence of dust shielding. The χ shield param-eter varies within less than one order of magnitude for di ff erent b parameters ranging from 1 to 10 km s − . In the H column den-sity range of interest ( N H ∼ − cm − ) all values are verysimilar. At higher column densities the shielding is dominated bydust.) all values are verysimilar. At higher column densities the shielding is dominated bydust.