Valeska Valdivia
email: valeska.valdivia@lra.ens.fr

Patrick Hennebelle
email: patrick.hennebelle@lra.ens.fr

Maryvonne Gérin

Pierre Lesaffre

A fast tree-based method for estimating column densities in Adaptive Mesh Refinement codes Influence of UV radiation field on the structure of molecular clouds

Keywords:

The interstellar medium (ISM) is a highly complex system. It corresponds to an intermediate scale between stars and galaxies. The interstellar gas is present throughout the galaxy, filling the volume between stars. A wide variety of coupled processes, such as gravity, magnetic fields, turbulence and chemistry, participate in its evolution, making the modeling of the ISM a challenging problem. A correct description of the ISM requires a good treatment of the magnetohydrodynamics (MHD) equations, gravity, thermal balance, and chemical evolution within the molecular clouds. This thesis work aims at a better understanding of the formation and evolution of molecular clouds, specially how they become "molecular", paying particular attention to the transition HI-to-H 2 . We have performed ideal MHD simulations of the formation of molecular clouds and the formation of molecular hydrogen under the influence of gravity and turbulence, using accurate estimates for the shielding effects from dust and the self-shielding for H 2 , calculated with a Tree-based method, able to provide fast estimates of column densities.

We find that H 2 is formed faster than predicted by the usual estimates due to local density enhancements created by the gas turbulent motions. Molecular hydrogen, formed at higher densities, could then migrate toward low density warmer regions.

Total H 2 column densities show that the HI-to-H 2 transition occurs at total column densities of a few 10 20 cm -2 . We have calculated the populations of rotational levels of H 2 at thermal equilibrium, and integrated along several lines of sight. These two results reproduce quite well the values observed by Copernicus and FUSE, suggesting that the observed transition and the excited populations could arise as a consequence of the multi-phase structure of molecular clouds. As H 2 formation is prior to further molecule formation, warm H 2 could possibly allow the development of a warm chemistry, and eventually explain some aspects of the molecular richness observed in the ISM.

VIII

Résumé

Le milieu interstellaire (MIS) est un système extrêmement complexe. Il correspond à une échelle intermédiaire entre les étoiles et les galaxies. Le gaz interstellaire est présent dans toute la galaxie, remplissant l'espace entre les étoiles. Une grande diversité de processus couplés, comme la gravité, le champs magnétiques, la turbulence et la chimie, participe à son évolution, faisant de la modélisation du MIS un problème ardu. Une description correcte du MIS nécessite un bon traitement des équations de la magnetohydrodynamique (MHD), de la gravité, du bilan thermique et de l'évolution chimique à l'intérieur du nuage moléculaire.

L'objectif de ce travail de thèse est une meilleure compréhension de la formation et de l'évolution des nuages moléculaires, et plus particulièrement de la transition du gaz atomique en gaz moléculaire. Nous avons réalisé des simulations numériques de la formation des nuages moléculaires et de la formation de l'hydrogène moléculaire sous l'influence de la gravité et de la turbulence MHD, en utilisant des estimations précises de l'écrantage par les poussières et de l'autoécrantage par la molécule H 2 . Ceci a été calculé grâce à une méthode en arbre, à même de fournir une rapide estimation des densités de colonne.

Nous avons trouvé que l'hydrogène moléculaire se forme plus rapidement que prévu par les estimations classiques du fait de l'augmentation de densité locale provoquée par les fluctuations turbulentes du gaz. L'hydrogène moléculaire, formé à des densités plus élevées, peut alors migrer vers les régions plus chaudes et moins denses.

Les densités de colonne totale d'hydrogène moléculaire montrent que la transition HI-H 2 se produit à des densités de colonne de quelques 10 20 cm -2 . Nous avons calculé les populations des niveaux rotationnels de H 2 à l'équilibre thermique et intégré le long de plusieurs lignes de visées. Ces résultats reproduisent bien les valeurs observées par Copernicus et FUSE, suggérant que la transition observée et les populations excitées pourraient être une conséquence de la structure multi-phasique des nuages moléculaires. Comme la formation de H 2 précède la formation des autres molécules, le H 2 chaud pourrait permettre le développement d'espèces endothermiques et éventuellement expliquer certains aspects de la richesse moléculaire observée dans l'ISM. IX X

List of Figures

x = n e /n H , from [START_REF] Dalgarno | Heating and Ionization of HI Regions[END_REF]. The discontinuity toward 10 4 K is due to the cooling by Lyman α emission, while at lower temperatures the cooling is due to CII and OI. . . 1.5 3D internal velocity dispersion of molecular clouds and condensations as a function of its maximum linear dimension L. The fit corresponds to σ v [km s -1 ] = 1.1L [pc] 0.38 . Figure from Larson (1981). Symbols are described in Table 1 in 1.13 Ultraviolet background in the solar neighborhood. Estimates from Habing (1968); Draine (1978) and [START_REF] Mathis | Interstellar radiation field and dust temperatures in the diffuse interstellar matter and in giant molecular clouds[END_REF] and observations from [START_REF] Henry | Far-ultraviolet studies. VII -The spectrum and latitude dependence of the local interstellar radiation field[END_REF] and [START_REF] Gondhalekar | Observations of the interstellar ultraviolet radiation field from the S2/68 sky-survey telescope[END_REF]. Figure from [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF] tional levels v = 0 to 3, and the rotational levels J for ortho and para-H 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Cooling function per hydrogen molecule, for n tot = 10 4 cm -3 and n(H)/n(H 2 ) = 10 -2 . Cooling functions have been calculated for two different OP R = 10 -4 and 3. From Le Bourlot et al. (1999).

3.1 Comparison of Sod test for two AMR codes (Flash and Enzo), and two SPH codes (Gadget-2 and Hydra). Bottom pannels shows two simulations of a galaxy merger performed with two of the most popular simulation techniques. The panel on the left: SPH simulation performed with the code Gadget-2 (image from Gadget-2 website http://www.mpa-garching. mpg.de/gadget). The panel on the right: AMR simulation using the Ramses code [START_REF] Teyssier | The Driving Mechanism of Starbursts in Galaxy Mergers[END_REF]). . . . . . . . . . . . . . . 3.3 Refinement strategy. The figure on the left shows the tree at different levels, while the figure on the right shows the leaf cells. 3.4 Tree-data structure in Ramses. Each oct, at any resolution level points toward its parent cell, to the 6 cells neighbors of its "parent", and it points to its 8 "children" octs. . . . . . . . . . . 3.5 Mesh in space and time showing cells used in the Godunov scheme. The Interstellar Medium

The interstellar medium (ISM) is the dilute material that fills the space between stars. Even if at first sight the space between stars seems to be empty, observations have unveiled a complex and rich structure. During the 17th century natural philosophers Francis Bacon and Robert Boyle discussed the nature of the space between stars, suggesting that interstellar space was not empty. Nebulous objects in the night sky have been observed since antiquity, but without understanding their nature. During the 2nd century A.D., Claudius Ptolemaeus recorded in his Almagest observations of nebulous stars. In the 18th century the french astronomer and comet hunter Charles Messier found fixed, faint and diffuse unknown objects, identified as "nebulae", that could have been mistaken for comets. Today we know that these "nebulae" correspond mainly to galaxies, planetary nebulae and interstellar clouds of gas. During the late 18th and early 19th centuries William and Caroline Herschel worked on a catalogue of nebulae. In On the Construction of the Heavens, [START_REF] Herschel | On the Construction of the Heavens[END_REF] gives the first representation of our galaxy using star counts. They noticed regions where the amount of stars was less, suggesting that the presence of matter in the line of sight could be responsible of the lack of stars. The birth of photography in the 19th century allowed astronomers to observe fainter objects. In 1919 the American astronomer E. E. Barnard found 182 dark nebulae in the sky, interpreted later as dense absorbing regions. The advent of radioastronomy during the 1930s gave a new impulse to the study of the interstellar medium. In particular, the discovery of the 21-cm line of atomic hydrogen during the 1950s [START_REF] Ewen | Observation of a Line in the Galactic Radio Spectrum: Radiation from Galactic Hydrogen at 1,420 Mc./sec[END_REF], and the CO molecular line at millimeter wavelengths, allowed the study of the interstellar gas both atomic and molecular, and the multiple instruments developed in recent years (satellite and ground-based) have opened the gates to wavelengths and spatial resolutions never reached before.

The ISM presents an intermediate scale between large structures and small objects, and the interplay between the different physical processes governs the star formation process. Its study is indeed a key step for understanding the structure formation and the evolution of galaxies.

The following section presents a general description of the interstellar medium, paying particular attention to the different physical processes present in the ISM.

The Cycle of Interstellar Matter

The interstellar medium is a turbulent, multiphase and highly inhomogeneous fluid that is in constant evolution. Stars are made of this insterstellar material and they form in the densest and coolest parts of molecular clouds. During the star formation process, gas undergoes dramatic changes. The total number density varies from a typical value of the order of 0.5 -1 cm -3 (for the diffuse medium) to several 10 6 cm -3 (for dense clumps). Temperature experiences a consequent variation. It drops from temperatures of the order of 10 4 K to some 10 K. When gravity overcomes all kind of supports, such as thermal pressure, the gravitational collapse can lead to the formation of a new star. In our galaxy the mass of gas that is converted into stars is about 1.3 M per year [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF]. During its lifespan, the star obtains its energy through thermonuclear fusion of hydrogen into helium, and modifies its environment through stellar winds and radiation. Once the star has run out its fuel, it dies as a supernova or as a planetary nebula containing a white dwarf in its center, depending on the initial mass of the star. In any case, the stellar evolution influences the ISM, enriching it in heavy elements and injecting kinetic energy, that contributes to support the interstellar turbulence. The stellar feedback in the Milky Way, through stellar and planetary winds and supernovae, send back to the ISM some 0.5 M per year. The intergalactic medium, the highly diluted gas that fills the space between galaxies, contributes to maintain the amount of gas present in the ISM. For the case of our galaxy the infall rate is about 0.5 M .

Composition of the ISM

Nowadays it is well known that the interstellar space is filled with magnetic field, cosmic rays, ions and molecules. The gas in the ISM inherits its chem-ical composition from the primordial nucleosynthesis. In terms of mass, the hydrogen counts for 70 %, the helium 28 %, and the heavy elements (also called metals) count only for 2 %. Apropos of hydrogen, it can be found in three forms. Roughly 60 % of the total mass of hydrogen is neutral atomic (HI), while the ionized hydrogen (HII) represents 23 %. The rest of the hydrogen can be found in dense and cold regions forming molecules (H 2 ).

Another important component of the ISM is the dust. Dust mass is insignificant (only 1 % of total H gas), but it plays a key role in the thermal balance, radiative transfer, and chemistry. It attenuates and scatters radiation (specially in the UV), it contributes to the gas heating through the photoelectric effect, and it acts as a catalyser for the molecule formation (specially H 2 ) on its surface.

Phases in the ISM

Interstellar gas is not homogeneous. It presents a highly complex structure and the different properties of the gas can vary over several orders of magnitude. In a seminal paper [START_REF] Mckee | A theory of the interstellar medium -Three components regulated by supernova explosions in an inhomogeneous substrate[END_REF] proposed a theory of the interstellar medium, where they identify three dominant phases. The phases they proposed have different physical conditions and are regulated by supernova explosions. The coldest and denser phase corresponds to the cold neutral medium (CNM), whith typical densities of 20 -50 [cm -3 ] and a temperature of 50 -100 [K]. The other two phases are warmer (6000 -10000 [K]), and less dense (0.2 -0.5 [cm -3 ]), but essentially different in their chemical compositions. The first of these warm phases is mostly neutral and atomic, known as the warm neutral medium (WNM), while the second, the warm ionized medium (WIM), is mainly ionized. Finally these three phases would be embedded in the hot ionized medium (HIM), or coronal gas, a tenous, fully ionized and extremely hot phase, with densities as low as 0.0065 [cm -3 ] and temperatures of the order of 10 6 [K] [START_REF] Ferrière | The interstellar environment of our galaxy[END_REF]. In this model, the different phases are spatially defined: the CNM occupies the innermost region of the cloud, shielded by layers of WNM and WIM. Nowadays we know that the global structure of the ISM is more complexe than this, but this classical picture is still valid for describing the interstellar gas. This vision is now more accurate. Following [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF], these phases (summarized in Tab. 1.1) can be described as: HIM: Hot Ionized Medium Also known as Coronal Gas (in analogy to the gas in the solar corona). This is the most diffuse phase of the ISM, composed mainly of ionized hydrogen. This is an extremely hot plasma and it has been heated on a relatively short timescale, primarily by supernova blast waves or hot superwinds from massive stars. Even though the heating is not constant, the cooling time in this phase is typically longer than 1 Myr. It fills a large volume of the galaxy, but the filling factor is not well known and it seems that it should be smaller than thought before.

Warm Ionized Medium (WIM) It is a plasma mainly composed of fully ionized hydrogen (HII). It fills about 10 % of the volume, and its extension is more than 1 kpc from the midplane of the galaxy. The gas is continuously ionized by Lyman continuum radiation from hot young stars, so its lifetime is close to that of the ionizing stars (typically ∼ 3 -10 Myr for O stars).

Warm Neutral Medium (WNM) Essentially warm atomic hydrogen, with a very low ionization fraction ( 8%). The gas in this phase presents low density (typically 0.5 cm -3 ) and a typical temperature of about 8000 K.

Recent estimates seems to point out that this phase would fill a volume of at least 30%. This phase contains about the 60% of the neutral atomic hydrogen and provides the bulk of the HI emission.

Cool Neutral Medium (CNM) This phase is composed of compact and cold atomic gas, with a very low volume filling factor (∼ 1%), forming fillamentary and sheet-like structures. This phase is in global pressure equilibrium wih the WNM (Field et al. 1969;Wolfire et al. 1995Wolfire et al. , 2003)). The cooling processes present in the ISM allow the gas to transit from the WNM to the CNM (Audit and Hennebelle 2005). The gas in this phase can be found associated with molecular gas.

Diffuse H 2 This is a molecular phase generally found embedded in the CNM, forming a subphase, where column densities are large enough to shield H 2 molecules against photodissociation. Although the gas is relatively denser than previous phases (∼ 100 cm -3 ), it is not dominated by self-gravity.

Dense H 2 Chemically, this phase is very similar to the previous one, but the main difference is that gravity plays a stronger role in this phase. Densities can reach values as high as 10 6 cm -3 , forming dense clumps and cores.

The main tool to study the molecular gas is CO line observations, but many other molecules are present in these phases. Values of density and temperature for the different components extracted from [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF].

These phases present a transition towards the necessary conditions for star formation. As we are interested on the formation of cold and compact structures, the scope of this work comprehends two fundamental steps in the star forming process. The first one is the transition that the gas experiences from the WNM toward the CNM, while the second one is the transition from the atomic to molecular hydrogen.

Heating and cooling processes in the diffuse ISM

A wide variety of thermal processes are present in the ism. These processes depend on the environmental physical conditions, but in return they determine the physical state of gas. In this section we describe the most important heating and cooling processes.

Heating

Photoelectric effect The photoelectric heating is the most efficient process for heating the diffuse gas. Interstellar gas bathes in a diffuse radiation field (Draine 1978). When the photon energy is higher than 13.6 eV, it will more likely ionize atomic hydrogen. Due to the high abundance of hydrogen, the medium has a large opacity for such photons, and the diffuse radiation field is devoid of photons with energies higher than 13.6 eV, but the remaining UV photons can still interact through photoelectric effect with heavier elements and dust grains that have lower ionization thresholds. Once the electron has been removed, it acquires an important fraction of the energy of the incident photon, and contributes to the heating of the gas through a thermalization process. The heating rate depends on the abundance and on nature of the species. Heavy elements are not abundant enough to significantly contribute to the gas heating through photoelectric effect. On the other hand, ultraviolet photoelectric emission from grains can explain by itself the different temperatures observed in the ISM [START_REF] Watson | Heating of Interstellar H i Clouds by Ultraviolet Photoelectron Emission from Grains[END_REF]. As very small grains (VSG) and polycyclic aromatic hydrocarbons (PAH) dominates in terms of total surface, and electrons can escape easily [START_REF] Puget | A new component of the interstellar matter -Small grains and large aromatic molecules[END_REF], Bakes and Tielens (1994) proposed the following heating rate for VSG : nΓ pe = 10 -24 n G 0 ergs cm -3 s -1 (1.1) where G is the incident radiation field and is the fraction of UV radiation absorbed by grains (with the values of Wolfire et al. (1995)) and it corresponds to: = 4.9 × 10 -2

1.0 + (G 0 T 1/2 /n e )/1925 0.73 + 3.7 × 10 -2 (T /10 4 ) 0.7 1.0 + (G 0 T 1/2 /n e )/5000 (1.2) n e is the electron density. Wolfire et al. (2003) give, in their Appendix C, an analytic expression by simplifying the treatment of PAHs, which reads:

n e = 2.4 × 10 -3 T 100 K 1/4
φ PAH (1.3) for PAHs, using G 0 = 1, a total ionization rate of ζ t = 10 -16 s -1 , and a parameter that scales the PAH collision rate φ PAH , that is equal to 0.5 in their standard model.

Collisional Ionization by cosmic rays Cosmic rays (CR), mostly low energy protons accelerated in supernovas, were the first heating source that has been considered [START_REF] Goldsmith | Thermal Properties of Interstellar Gas Heated by Cosmic Rays[END_REF]). This process is not dominant in the diffuse ISM, but it is the main heating and ionizing process when gas becomes very dense and UV photons cannot penetrate. The heating rate by cosmic rays depends, in principle, on the energy of the ionization rate, but this rate can vary over two orders of magnitude in the galaxy depending on the environment. Goldsmith (2001) has estimated an intermediate value for the cosmic ray heating rate, given by: nΓ cr = 10 -27 n ergs cm -3 s -1 (1.4)

Mechanical heating Macroscopical mouvements in the gas, including turbulence and shocks, can inject energy and contribute to heating through dissipative processes. The mechanical heating is the most difficult to evaluate due to the extremely compressive, non adiabatic and magnetized nature of the ISM, and due to the lack of knowledge about the nature of the dissipative processes themselves. Nevertheless, it has been shown that dissipative processes can contribute to heat the gas. In particular, the works of [START_REF] Scalo | Heating of dense interstellar clouds by magnetic ion slip -A constraint on cloud field strengths[END_REF] (in the dense gas) and [START_REF] Ferriere | Hydromagnetic wave heating of the low-density interstellar medium[END_REF] (in low-density gas) have shown that the dissipation of magnetic waves, through ambipolar diffusion, can contribute substantially to heat the gas.

Similarly [START_REF] Spitzer | Acoustic waves in supernova remnants[END_REF] has shown that the dissipation of acoustic waves produced by shocks in supernova remnants can contribute to heat the diffuse phase. [START_REF] Falgarone | The intermittency of turbulence in interstellar clouds: implications for the gas kinetic temperature and decoupling of heavy particles from the gas motions[END_REF] and [START_REF] Joulain | Non-equilibrium chemistry in the dissipative structures of interstellar turbulence[END_REF] have addressed the question of whether a burst of dissipation of turbulent energy could heat the gas through viscous dissipation and eventually drive a "warm" chemistry. In a more recent work Hennebelle and Inutsuka (2006) have shown that a mechanical heating seems to be crucial to the survival of warm neutral gas inside molecular clouds. They have found that turbulent energy dissipation can allow the existence of warm gas in low-pressure environments, while in high-pressure gas the dissipation of MHD waves can maintain the warm gas.

From a numerical point of view, for large scale simulations hardly ever these dissipative scales are well described, but the injection of mechanical energy entrains an additional heating of gas that can be estimated by knowing the other heating and cooling rates.

Cooling

The cooling in the diffuse medium is mainly due to collisionally excited lines.

The mechanism of the radiative cooling is quite simple: during the collision a part of the kinetic energy of the gas is transferred to the particle (mostly CII and OI in the diffuse gas), and excites its internal energy levels. Radiative decay follows, emitting a photon, that can escape from the gas if it is optically thin. The energy of the photon will be the energy difference between the two levels where the transition occurs. The efficiency of the transition as a coolant depends on the frequency of collisions, on the excitation energy required (comparable or less that thermal kinetic energy), on the probability of excitation during the collision, the emission of a photon before a new collision, and on the ability of the photon to leave the gas before being re-absorbed. For example, if the photon is produced during an allowed transition, it will be more likely absorbed by the gas and will not contribute to the cooling. That is the reason why the main cooling lines are produced by forbidden transitions.

In the low density approximation, where every collisional excitation is followed by the emission of a photon, for a particle is collisionally excited from level l (lower) to an excited level u (upper), the cooling rate can be written

n 2 Λ ul = n X n l C lu hν lu [erg s -1 ] (1.5)
where n X is the density of colliders, n l is the density of particles in the lower level contributing to the cooling, and C lu is the collisional excitation probability. Assuming local thermal equilibrium (LTE), both collisional probabilities of excitation and de-excitation are related by

C lu = C ul g u g l exp - hν kT (1.6)
and writing the Boltzmann-like dependence on temperature for the de-excitation probability, including the collision strength, Ω ul = Ω lu , and assuming that it does not vary with v

C ul = 2π kT 1/2 2 m 3/2 e Ω ul g u (1.7)
then we can derive the following general expression for the rate of energy loss

n Λ ul = n X n l 2π kT 1/2 2 m 3/2 e Ω ul g l exp - hν kT [erg s -1 ] (1.8)
It is important to notice that the collision strength value comes from a complicated and detailed calculation of atomic and molecular physics [START_REF]Oscillator strengths and effective collision strengths for electron excitation of Mg VI[END_REF]. Fine structure lines of CII Figure 1.2 shows, on the left hand side, the internal energy structure of ionized carbon (CII) for the ground state 2 P 1/2 and the first excited state 2 P 3/2 . The energy difference between both levels is ∆E = 1.3 × 10 -14 ergs = 92 K, and the transition from 2 P 3/2 toward 2 P 1/2 of CII at 157.7 µm is the main contribution to the cooling at low temperature. CII is predominantly excited by HI and by electrons. [START_REF] Launay | Fine-structure excitation of ground-state C/+/ ions by hydrogen atoms[END_REF] have calculated the collision rate between CII and HI, while Hayes and Nussbaumer (1984) calculated the collision rate with electrons in the low electron density limit. Then using the expression from Wolfire et al. (1995), the cooling rate by the CII fine-structure transition can be written:

nΛ CII = 2.54 × 10 -14 A C f CII γ HI n HI + γ e n e exp -92 T ergs s -1 (1.9)
where A C is the total abundance of carbon, f CII the fraction of singly ionized carbon. γ HI and γ e are, respectively, the collisional de-excitation rates for collisions with neutral hydrogen, and with electrons, that can be written γ HI = 8 × 10 -10 T 100 0.7

(1.10) γ e = 2.8 × 10 -7 T 100 -0.5 (1.11) Fine structure lines of OI At higher temperatures (T ≥ 100 K) fine levels of other atoms and ions can be collisionally excited. Figure 1.2 displays, on the right hand side, the internal energy structure of neutral oxygen (OI) for the ground state ( 3 P 2 ) and the two first excited levels ( 3 P 1 and 3 P 0 ). For neutral oxygen (OI), two transitions contribute to the cooling. The first one is a suprathermal emission at 145.52 µm ( 3 P 0 → 3 P 1 ). The second one is a subthermal emission at 63.18 µm ( 3 P 1 → 3 P 2 ). The rates of radiative cooling by OI in both transitions, using the rates of collisional excitation for OI from [START_REF] Launay | Fine structure excitation of carbon and oxygen by atomic hydrogen impact[END_REF], as given in Flower et al. (1986) are (1.12)

n 2 Λ OI = 10 -26 T
where n HI and n OI are the number density of neutral hydrogen and neutral oxygen, respectively.

Lyman alpha emission At temperatures of a few thousands degrees neutral hydrogen can be collisionally excited by electrons, and the Lyman α emission line becomes important [START_REF] Spitzer | Physical processes in the interstellar medium[END_REF]. The mechanism that assures the cooling is different from the previous ones. The Lyman α photons are absorbed by dust, that re-emites the absorbed energy in the form of infrared wavelengths photons, that can escape easily from the cloud. The main contribution to the cooling comes from the transition from the first excited level n = 2 toward the ground level, at 121.6 nm, and the cooling rate will be nearly proportional to exp(-∆E12 kT ), where ∆E 12 is the energy difference between the first exited level and the ground level. The cooling rate can be expressed as follows

n 2 Λ Ly = 7.3 × 10 -19 n e n HI exp -118400 T (1.13)
where n e and n HI are the number density of electrons and the number density of neutral atomic hydrogen, respectively.

Electron recombination onto positively charged dust grains Charged particles can stick onto grain surfaces and remove thermal energy from the surrounding gas when the thermal energy of particles that recombine on grain surfaces is higher that the energy of ejected electrons. This electron recombination mechanism can become important at high temperatures, T 10 4 K, due to the increased Coulomb interaction caused by the increase in grain charge (Bakes and Tielens 1994). Wolfire et al. (1995Wolfire et al. ( , 2003)), based on the analytical fit calculated by Bakes and Tielens (1994), give the cooling rate

n 2 Λ rec = 4.65 × 10 -30 T 0.94 G 0 T 1/2 n e φ PAH β n e φ PAH n [erg cm -3 s -1 ] (1.14)
with β = 0.74/T 0.068 , the parameter φ PAH and the electron density, n e , are the same as in Eq. 1.3.

These are the main cooling processes that take place in the atomic gas, but molecules play a key role in the evolution of the ISM. In particular H 2 , the most abundant molecule in the ISM, plays a role in the thermal balance. This will be analized in more detail in Section 2.4. Wolfire et al. (1995). Field et al. (1969) have shown that neutral atomic gas is thermally bistable, and in the solar neighborhood both phases exist at roughly P/k = 10 3 -10 4 cm -3 K [START_REF] Jenkins | Copernicus observations of C I -Pressures and carbon abundances in diffuse interstellar clouds[END_REF]. Wolfire et al. (1995Wolfire et al. ( , 2003) ) have calculated the thermal equilibrium gas temperature for the diffuse neutral medium including several heating and cooling processes. They have found that, under the conditions present in the solar neighborhood, the photoelectric heating produces two stable phases, and that the rapid drop in temperature, due to CII line cooling, leads to a thermal instability between both phases. Figure 1.3 shows the curve of thermal equilibrium for their standard model, where both CNM and WNM coexist.

Thermal bistability in the ISM

The bistable nature of the interstellar gas arises from the shape of the cooling function. Fig. 1.4 (Dalgarno and McCray 1972) shows how the cooling function [START_REF] Dalgarno | Heating and Ionization of HI Regions[END_REF]. The discontinuity toward 10 4 K is due to the cooling by Lyman α emission, while at lower temperatures the cooling is due to CII and OI. is steeper at lower temperatures, and becomes flatter at temperatures between 100 and 8000 K. In the first region, when gas cools and temperature decreases, the efficiency of the cooling decreases strongly and the gas cannot cool more. On the second region, as the cooling is mostly insensitive to temperature, it maintains its efficiency and the energy loss continue. The first case is the stable region, while the second one is unstable.

Turbulence

The first evidences of the turbulence in the ISM were found through the kinematic analysis of molecular clouds. Several works aimed at characterizing the nature of the internal velocity distribution observed in molecular clouds and in HII regions. The velocity fluctuations observed in HII regions were found to be in good agreement with the spectral law of turbulent motions and consistent with the predictions of the theory of Kolmogorov (von Hoerner 1951;[START_REF] Courtes | Astrophysique : vitesses radiales de l'hydrogene interstellaire EN E mission des regions H II lointaines[END_REF][START_REF] Münch | Internal Motions in the Orion Nebula[END_REF], in particular [START_REF] Wilson | Internal Kinematics of the Orion Nebula[END_REF] showed that turbulent gas could be described by the Kolmogorov theory for incompressible turbulence, even though interstellar gas might be highly compressible. Larson (1981) found a powerlaw relation between (Fig. 1.5) the velocity dispersion σ v and the characteristic size L of the structure, σ v ∝ L γ , where γ ≈ 0.38 consistent with the index γ = 1/3 from the Kolmogorov theory. (Solomon 0.38 . Figure from Larson (1981). Symbols are described in Table 1 in the same work, but among the most interesting regions we can recall the Taurus Molecular Cloud Complex (T), the ρ Ophiuchus Complex (ρ), and the Orion Complex (O). et al. 1987) perfomed a study of 273 molecular clouds in the Galactic disk, finding an index of γ ≈ 0.5, steeper that the value found by Larson. More recently [START_REF] Heyer | The Universality of Turbulence in Galactic Molecular Clouds[END_REF] studied 27 GMCs, finding an index of γ = 0.59. The nature of turbulence is essentially a multiscale phenomenon, where energy is transferred from the large scales toward the small dissipative scales. The existence of a scale law relating velocity dispersion and size enforces the idea that the hierarchy of structures present in the ISM has its origin likely in the turbulence.

In recent years the new observational capabilities have permitted us to reach smaller scales, unveiling a highly structured and fragmented medium. Emissionline observations of molecular clouds show clumps and filamentary structures over a large range of scales [START_REF] Falgarone | The small-scale density and velocity structure of quiescent molecular clouds[END_REF][START_REF] Falgarone | Small-Scale Density and Velocity Structure of a Molecular Cloud Edge[END_REF][START_REF] Wiesemeyer | High resolution studies of protostellar condensations in NGC 2024[END_REF]. The mass spectrum of clumps is characterized by a power law, revealing the self-similarity of the distribution of structures over several orders of magnitude from the largest structures, the Giant Molecular Clouds, with scales of few tens to hundred pc, to the smallest structures, the protostellar cores, with scales of the order of 0.01 pc.

But not only kinematical observations highlight the turbulent nature of the ISM. The presence of molecules whose formation is governed by highly endothermic reactions, requires additional sources of energy able to open new formation paths. Turbulent dissipation could supply such energy. This has motivated the work of [START_REF] Godard | Chemical probes of turbulence in the diffuse medium: the TDR model[END_REF], who, by coupling a turbulent dissipation model to a chemical network, have been able to produce the observed abundances of such species. In the same work they have found that CH + and SH + in addition to the excited levels of H 2 are the most sensitive species to turbulent dissipation, and they can be used to estimate the turbulent dissipation rate, as to constrain the ion-neutral velocity drift in such regions.

All these works seem to stress the key role of turbulence in the evolution of molecular clouds and in the regulation of star formation in the galaxy. Recent theoretical and observational works have payed special attention to the role of turbulence in the evolution of molecular clouds and on the star formation process (Mac [START_REF] Low | Control of star formation by supersonic turbulence[END_REF][START_REF] Elmegreen | Interstellar Turbulence I: Observations and Processes[END_REF][START_REF] Scalo | Interstellar Turbulence II: Implications and Effects[END_REF][START_REF] Mckee | Theory of Star Formation[END_REF][START_REF] Kennicutt | Star Formation in the Milky Way and Nearby Galaxies[END_REF]. Large scale turbulent motions may induce density fluctuations that can lead to molecular cloud formation (Ballesteros-Paredes et al. 1999), while turbulence at smaller scales is the main support against selfgravity for molecular clouds, regulating the star formation rate (SFR). Turbulence at intermediate scales alters the amount of gas present in different phases of the ISM, in particular the amount of unstable gas transiting between the WNM and the CNM is correlated to the amount of turbulence (Audit and Hennebelle 2005).

Nowadays we know that the turbulence in the ISM is highly compressible, magnetized, multi-phase and interacts with gravity, but its properties are still controversial.

Incompressible hydrodynamical turbulence

The first attempt to characterize turbulence was done during the late 15th and early 16th centuries by Leonardo da Vinci. He experimented with obstructions in water flows and free water jets falling into pools, being a pioneer in flow visualization (Fig. 1.6). He wrote: "Observe the motion of the surface of the water, which resembles that of hair, which has two motions, of which one is caused by the weight of the hair, the other by the direction of the curls; thus the water has eddying motions, one part of which is due to the principal current, the other to random and reverse motion."

With this simple description da Vinci identified the injection of energy, the presence of vorticity and the turbulent cascade, which tranfers energy from one scale to the smaller ones, by describing the disruption of bigger eddies into smaller ones.

Turbulence is characterized by random motions within the fluid. The first modern definition of turbulence was built experimentally in the late 19th century by Osborne Reynolds, who characterized the transition from laminar to turbulent flow in a pipe flow. This transition corresponds to an instability that develops when the inertial forces, v • ∇ v, overcome viscous forces, ν∇ 2 v, that stabilize the fluid (here ν is the kinematic viscosity of the fluid). A fluid can be characterized by the ratio between these two forces, called the Reynolds number, which is defined as follows

Re = inertial forces viscous forces = ρv 2 L L 2 µv L L = v L L ν (1.15)
the dynamic viscosity µ and the kinematic viscosity ν are related by ν = ν/ρ. v L is the characteristic velocity at scale L. When Re 1 diffusion dominates and the flow is laminar. On the other hand, when Re 1 the flow becomes unstable. Viscosity is no longer able to stabilize the flow, inertial forces dominate, the flow becomes unstable, and turbulence develops. It is important to notice that to the present day there is not an exact treatment of turbulent flows, and that in

Kolmogorov theory

Since the early works of da Vinci and the characterization of a turbulent state by Reynolds, there were no new developments in the theory of turbulence until 1939-1941 with the works of Andrey Kolmogorov and Alexander Obukhov [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF][START_REF] Obukhof | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF]. The theory was completed later by taking into account the objection of Lev Landau (1949) on the universality of the energy transfer [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number[END_REF][START_REF] Obukhov | Some Specific Features of Atmospheric Turbulence[END_REF]. Their works were based on the turbulent cascade mechanism proposed by Lewis Fry Richardson during the 1920s, who described the mechanism in the form of a rhyme: Big whorls have little whorls That feed on their velocity And little whorls have lesser whorls And so on to viscosity

The assumptions of Richarson's idea, that constitute the base of the theory of Kolmogorov, can be summarized as follows:

• There exist an "external scale" L, at which energy is injected • There exist an "internal scale" l D , at which energy is dissipated • There exist in the turbulent flow vortices on all possible scales L < l < l D • There exist a certain uniform mechanism of energy transfer from the coarser-scaled vortices to the finer.

The central idea in Kolmogorov theory consist in the phenomenological description of the turbulence as a hierarchy of vortices set up by the energy transfer rate, which is independent of the scale. Further assumptions are needed: boundary influences (external forcing) are ignored, the fluid is supposed spatially isotropic and homogeneous (statistically identical in any direction).

The turbulence is supposed to be composed of eddies of different sizes, each one characterized by a scale l, a velocity v l , and a turnover timescale τ l = l/v l . Largest eddies have sizes comparable to the flow length scale. At high Re the fluid becomes unstable and eddies decay into smaller ones, transferring a part of the energy to the smaller scales. If at the new scale the fluid is still unstable, eddies will break into smaller ones until viscous forces will be able to stabilize the fluid again. Through this process energy cascades from the largest scale to the smallest scale, where energy can be disipated. At scale l the energy rate transfer can be estimated dimensionally in the following manner

∝ v 2 l τ l = v 3 l l → v l ∝ ( l) 3/2 (1.16)
Due to the isotropy and homogeneity assumptions, we can estimate the kinetic energy at scale l as

E l ∝ 1 2 v • v = 1 2 v 2 x + v 2 y + v 2 z ≈ 1 3 (v l ) 2
(1.17)

The energy spectrum describes the energy distribution among eddies at all scales. Associating a wavenumber k = 2π/l to the scale l, and using dimensional arguments, we can deduce the energy spectrum as

E(k) ∝ E l k ∝ 1 k ( l) 2/3 E(k) = C 2/3 k -5/3 (1.18)
which is the Kolmogorov, or Kolmogorov-Obukhov, spectrum. C is the universal Kolmogorov constant, which value, C = 1.5, has been obtained experimentally. This spectrum is valid for incompressible fluids (purely solenoidal modes) in the inertial range L < l < l D . The power spectrum P( k) is related to the energy spectrum by

E(k) = P( k)δ(| k| -k)d k (1.19)
for a Kolmogorov energy spectrum, P( k) ∝ | k| 1-D-5/3 , where D is the space dimension, then P(k) ∝ k -5/3 in 1D and P( k) ∝ k -11/3 in 3D.

Compressible turbulence

While the Kolmogorov theory describes the incompressible turbulence, the velocity dispersions observed in the ISM shows that the turbulence is supersonic, and thus the medium must be highly compressible. As this problem is extremely complex, we can consider a simpler case given by a purely compressive

∂v x ∂t + v x ∂v x ∂x = ν ∂ 2 v x ∂x 2 (1.20)
where v x is the velocity, and ν is the viscosity. This equation can be solved exactly by using the Cole-Hopf transformation, which reads v x = -2ν∂φ/∂x, that allows to linearize the system. The energy spectrum for the Burgers equation is [START_REF] Girimaji | Spectrum and energy transfer in steady Burgers turbulence[END_REF], and corresponds to pure compressive turbulence.

E(k) ∝ k -2 (
Numerical simulations have shown that turbulence is composed by a mixture of solenoidal and compressive modes (Kritsuk et al. 2007).

MHD turbulence

As we recalled at the beginning of this section, turbulence in the ISM is magnetized, so if we want to shed some light on the physics behind the formation and evolution of molecular clouds, and on the star formation process, we need to deal with turbulence in a magnetized medium. The pioneering works of [START_REF] Iroshnikov | Turbulence of a Conducting Fluid in a Strong Magnetic Field[END_REF] and [START_REF] Kraichnan | Inertial-Range Spectrum of Hydromagnetic Turbulence[END_REF] were the first attemtps to characterize the MHD turbulence. Their main assumptions were to consider spatial isotropy, which lead to eddies with nearly the same dimensions in directions aligned and perpendicular to the magnetic field lines, and the existence of a local mean magnetic field. They derived a 1D energy spectrum of E(k) ∝ k -3/2 . [START_REF] Shebalin | Anisotropy in MHD turbulence due to a mean magnetic field[END_REF] showed that the turbulence was necessarily anisotropic. [START_REF] Goldreich | Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence[END_REF] performed a further development of the compressible MHD turbulence by allowing eddies to be distorted following field lines. This allows cascades to be also anisotropic; energy is transferred faster along k ⊥ than along k with respect to the field lines. For the perpendicular direction they found an energy spectrum E(k ⊥ ) ∝ k -5/3 ⊥ .

Magnetic Fields in the ISM

Magnetic fields participate in a wide variety of astrophysical processes, from sunspots to radiogalaxy jets. They play an extremely important role in the physics of interstellar gas. Different observations have shown that magnetic fields pervade the ISM. At the scale of the Galaxy the mean value of the magnetic field is ∼ 5µG. It is partially organised, but it presents an irregular component. In galaxies, and more precisely in molecular clouds, the magnetic fields contribute to increase the pressure in two ways: additional magnetic pressure, carried by magnetic field, and additional pressure coming from high energy charged particles trapped by the magnetic field.

Figure 1.8: The Galactic magnetic field as revealed by Planck. From Juan Soler webstite http://www.ias.u-psud.fr/soler/planckhighlights.html.

Interstellar magnetic fields play an important role in the structure, dynamics and evolution of molecular clouds. They can participate in the regulation of star formation by allowing the evacuation of prestellar angular momentum, and at the same time giving an additional support against gravity. They are central to the acceleration and propagation of cosmic rays, that are the main source of ionization, and they contribute to heat the inner part of the molecular clouds. Therefore their study is imperative to the understanding of several processes. Some reviews on this topic, highlighting the role of magnetic fields for the physics of the ISM, and describing different methods to detect, measure and map magnetic fields in the space, have been published in recent years. Among those that have inspired this section we can recall [START_REF] Ferrière | Interstellar magnetic fields[END_REF]; [START_REF] Zweibel | Magnetic fields in galaxies and beyond[END_REF]; [START_REF] Crutcher | Magnetic Fields in Molecular Clouds[END_REF][START_REF] Crutcher | Observations of Interstellar Magnetic Fields[END_REF].

In 1937 the swedish physicist Hannes Alfvén (Alfvén 1937a,b) argued that if the universe was permeated with plasma, charged particles would carry electric currents able to generate a galactic magnetic field, and pointed out that such magnetic fields would be able to confine cosmic rays. Based on this, [START_REF] Fermi | On the Origin of the Cosmic Radiation[END_REF] proposed that cosmic rays would be originated and accelerated in the ISM.

The first observational evidence of the presence of magnetic fields came with the observation of linear polarization of optical radiation from stars [START_REF] Hall | Observations of the Polarized Light from Stars[END_REF]Hiltner 1949a,b), in particular photometric observations of binary systems showed that the polarization likely originated in the interstellar gas rather than in the stars themselves. [START_REF] Davis | The Polarization of Starlight by Aligned Dust Grains[END_REF] showed that this polarization could arise as a consequence of anisotropic absorption and scattering by spinning alongated grains aligned with magnetic fields. Another important piece of observational evidence arose from observations of the Galactic radio emission. [START_REF] Kiepenheuer | Cosmic Rays as the Source of General Galactic Radio Emission[END_REF] argued that the Galactic radio emission could not be thermal, because it would require electronic temperatures of the order of 10 5 K (in other words ten times higher than the value obtained spectroscopically), and the nature of the radio emission would be the synchrotron radiation from spiraling charged particles confined by the magnetic fields. [START_REF] Shklovskii | Radioastronomiia[END_REF] showed that the optical emission from the Crab nebula would also be powered by synchrotron emission.

Another step in the search of magnetic fields was proposed by [START_REF] Bolton | On the Possibility of Measuring Interstellar Magnetic Fields by 21-CM Zeeman Splitting[END_REF], who suggested to search for the Zeeman splitting of the 21 -cm line of neutral hydrogen. In presence of weak magnetic field this line is split into three components due to the interaction of the magnetic moment of the atom and the magnetic field. Several attemps were made in the following years [START_REF] Galt | An attempt to detect the galactic magnetic field, using Zeeman splitting of the hydrogen line[END_REF][START_REF] Davies | A New Limit to the Galactic Magnetic Field Set by Measurements of the Zeeman Splitting of the Hydrogen Line[END_REF][START_REF] Weinreb | An Attempt to Measure Zeeman Splitting of the Galactic 21-CM Hydrogen Line[END_REF]), but results could only give upper limits to the strength of the Galactic magnetic field, due to the complexity of the measurement (small differences in gain for both polarizations, corrections that could introduce a spurious Zeeman effect). Since then, new instruments have been developed, able to probe interstellar magnetic fields, but observations are still challenging.

Observational Methods

Evidence for interstellar magnetic fields come from different observational methods.

Zeeman Splitting This method is used to obtain the strength and the sign of the component parallel to the line-of-sight of the magnetic field B .

The Zeeman splitting arises from the coupling between the magnetic moment of the atom, which is almost equal to the moment of the valence electrons, and the external magnetic field B. In the case of the hydrogen atom the splitting is about h∆ν ∼ µ B B, where µ B = e /2m e c is the Bohr magneton. As the splitting is proportional to the strenght of the field, it would be sufficient to measure the splitting to infer the magnetic field strength, but the level splitting caused by fields of strength 1 -100 µG is about ∆ν ∼ a few hertz, which is too small to be detected. The Zeeman splitting was first observed in the 21 -cm in absorption by [START_REF] Verschuur | Positive Determination of an Interstellar Magnetic Field by Measurement of the Zeeman Splitting of the 21-cm Hydrogen Line[END_REF][START_REF] Verschuur | Positive Determination of an Interstellar Magnetic Field by Measurement of the Zeeman Splitting of the 21-cm Hydrogen Line[END_REF] and later in emission by [START_REF] Heiles | Measurements of magnetic field strengths in the vicinity of Orion[END_REF].

The Zeeman splitting is very difficult to detect, nevertheless, the Zeeman effect produces a difference in the frequencies of both circular polariza-tions, thus the Stokes parameter V , which is proportional to B , is measured. The usual values obtained for the magnetic field in atomic clouds is B ∼ 6 µG, and for in molecular clouds is about B ∼ 10 -3000 µG [START_REF] Crutcher | Magnetic Fields in Molecular Clouds: Observations Confront Theory[END_REF][START_REF] Crutcher | Magnetic Fields in Molecular Clouds[END_REF].

Faraday Rotation The Faraday rotation is used in ionized regions, and just like the Zeeman splitting, it gives the strength and the sign of the component parallel to the line-of-sight of the magnetic field B .

Figure 1.9: Electromagnetic wave is rotated under the influence of the magnetic field.

The refraction index in a plasma is frequency-dependent, so the electromagnetic waves that propagates in the plasma will experience a delay that will depend on the frequency (lower frequencies will be the most delayed).

In the presence of a magnetic field, electrons will spiral around field lines in one specific direction, causing a different refraction index for right and left circularly polarized light. Linear polarization can be decomposed into two circular polarizations of opposite direction, thus as linear polarized light travels through the plasma, the phase velocity of both opposite circular polarizations will be different, producing a rotation of the plane of polarization of the wave. The rotation angle is given by ∆θ = RMλ 2 , where RM = C L 0 n e B ds is the rotation measure. C is a numerical constant, n e is the density of free electrons, and L is the length of the path that the radiation has cross through the plasma. Then B can be obtained by measuring the rotation angle for at least two frequencies. This technique has permitted the astronomers to characterize the Galactic magnetic field. We know that it has a regular and a dominating turbulent component, with values B reg 1.5µG and B turb 5µG in the vicinity of the Sun [START_REF] Rand | The local Galactic magnetic field[END_REF].

Synchrotron radiation Synchrotron radiation is produced by relativistic charged particles spiraling around magnetic field lines. This emission comes essentially from free electrons, and its intensity depends on both the magnetic field strength and the energy density of electrons. Assuming a power-law energy spectrum of relativistic electrons of f (E) = K e E -γ , the emissivity of synchrotron radiation is given by

E ν = F(γ)K e B (γ+1)/2 ⊥ ν -(γ-1)/2
, where F(γ) is a funtion of the spectral index of relativistic electrons, and B ⊥ is the strength of magnetic field projected on the plane of the sky. Furthermore, if the magnetic field is coherent, or ordered, the synchrotron emission will be linearly polarized and perpendicular to B ⊥ , and the polarized synchrotron emission will probe the strength and orientation of the ordered magnetic field B ord . The magnetic field in the ISM contains an ordered and a fluctuating components. Near the Sun, B ord /B 0.6 [START_REF] Beck | Galactic and Extragalactic Magnetic Fields[END_REF], which gives B ord 3µG.

Dust polarization This method is useful in dusty regions and gives only the orientation of the magnetic field projected onto the plane of the sky.

The mechanism is not fully understood, but observations seem to point out that irregular grains are subjected to radiative torques, causing them to spin around its shortest axis, and dust grains would gradually align their spinning axis with the magnetic field lines. Elongated dust grains block preferentially the component of starlight that is aligned with their long axis, which is perpendicular to the field line, producing a linear polarization. The dust thermal emission from the spinning grains is also anisotropic and contributes to visualize the magnetic field. Its maximum emission is aligned with the longest axis, in orther words, it is perpendicular to the magnetic field line.

The results of the combined observations show that the magnetic field is horizontal in the Galactic disk, and nearly azimuthal in the solar vicinity.

Origin and Maintenance of Galactic Magnetic Field

The origin and the actual mechanisms that maintain the Galactic magnetic field are still unclear [START_REF] Widrow | Origin of galactic and extragalactic magnetic fields[END_REF]), but it is widely accepted that magnetic fields were amplified by fluid motions from early cosmological, and extremely weak, fields.

The Primordial Theory postulates that galaxies are formed from galactic disks having large-scale magnetic fields and it supposes low magnetic diffusivity. It is hard to maintain a primordial magnetic field due to differential rotation, which would cause the field lines to wrap into a bisymmetrical spiral structure. In this configuration field lines with opposite polarity can reconnect and cancel. A second theory, that has been used before to explain different magnetic fields in other contexts, is the Dynamo Theory [START_REF] Beck | Magnetic spiral arms in the galaxy NGC6946[END_REF][START_REF] Brandenburg | Astrophysical magnetic fields and nonlinear dynamo theory[END_REF]. This theory supposes that the magnetic diffusivity is high, leading to magnetic flux leakage (Ferrière 1993a,b), and thus to a decay of the field, that is continuously regenerated through dynamo effect induced by large scale flows and small scale turbulent motions. Turbulent amplification is crucial for the efficiency of the dynamo effect. As we saw in the previous section, magnetic fields and turbulence are in close interaction.

Instruments

Recent and upcoming instruments represent huge leaps forward in how we understand cosmic magnetic fields, by offering more sensitivity, wider coverage of the sky and higher resolution. Among these "new generation" instruments we can cite the following projects: The most relevant studies for the physics of the ISM are "by-products" of the primary mission of measuring the CMB anisotropies. In particular, polarized dust emission maps have revealed an intrincate network of filaments in the ISM unveiling the projected structure of the Galactic magnetic field and its relative orientation (Planck Collaboration et al. 2014b,a, 2015b,a), Fig. 1.8 shows an all-sky map of the column density (in colors) and the Galactic magnetic field orientation obtained with Planck.

LOFAR (LOw
SKA (Square Kilometer Array) This is the largest array scheduled. It will be situated in the southern hemisphere in two hosting sites: Australia and South Africa. The final project (SKA2) will count with a total collecting area of 1km 2 , constituted by 250000 antennas and 2500 dishes. SKA will operate in the 50 MHz -3 GHz frequency range. There are several projects aiming at characterize magnetic fields in general via the Zeeman effect [START_REF] Robishaw | Measuring Magnetic Fields Near and Far with the SKA via the Zeeman Effect[END_REF], synchrotron emission from molecular clouds [START_REF] Dickinson | SKA studies of in situ synchrotron radiation from molecular clouds[END_REF], Faraday rotation [START_REF] Bonafede | Unravelling the origin of large-scale magnetic fields in galaxy clusters and beyond through Faraday Rotation Measures with the SKA[END_REF] among others. As LOFAR, it will be able to measure the Faraday rotation and the synchrotron emission, but it will be 100 times more sensitive than LOFAR.

PILOT (Polarized Instrument for Long wavelentgh Observation of the Tenuous ISM ) PILOT is a balloon-borne instrument currently developed at the IRAP (Institut de Recherche en Astrophysique et Planétologie) in Toulouse (France), specially conceived to measure polarized emission from dust grains present in the diffuse interstellar gas. It comprises a 1 m telescope and a photometer designed to operate in two photometric channels at wavelengths 240 and 550 µm. The instrument is still in the testing phase, but it will be completely dedicated to the study of the ISM from August 2015. Interstellar gas is strongly influenced by the interstellar radiation field (ISRF) in which it bathes. Radiation pressure can influence the dynamics of molecular clouds, gas and dust can be accelerated by anisotropic radiation, and radiation pressure can even disrupt molecular clouds [START_REF] Murray | The Disruption of Giant Molecular Clouds by Radiation Pressure & the Efficiency of Star Formation in Galaxies[END_REF]. The physical state of the gas is regulated by the radiation field. If photons are sufficiently energetic, they can either remove electrons from atoms, molecules, and grains, or dissociate molecules, regulating the chemical and ionization state of the gas. Temperature is highly dependent on the ISRF, in particular the heating in the ISM comes primarily from photoelectrons ejected from grains.

Interstellar Radiation Field

The distribution of the radiation is not constant throughout the Galaxy and it varies with the wavelength. Fig. 1.12 shows the dominant components of the ISRF in the solar neighborhood.

The Galactic synchrotron emission is produced by relativistic electrons spiraling around field lines (discussed in Section 1.5), which is dominant at frequencies ν 1 GHz. The all-sky survey at 408 MHz from [START_REF] Haslam | A 408 MHz all-sky continuum survey. II -The atlas of contour maps[END_REF] shows that the emissivity is higher toward the galactic plane and near supernova remnants. The second component of the ISRF is the fossil emission from the cosmic microwave background (CMB). It was discovered by [START_REF] Penzias | A Measurement of Excess Antenna Temperature at 4080 Mc/s[END_REF] as an excess, or residual noise, that was isotropic, constant in time and non-polarized. This emission is well described by a blackbody at a temperature of approximatly 2.73 K. Another component that contributes to the ISRF at frequencies ν 100 GHz is the free-free, free-bound and bound-bound transitions from a plasma at temperature T ∼ 10 4 . Of particular interest are the contribution from the infrared emission of dust, composed by the thermal emission of dust grains and molecular bands arising from polycyclic aromatic hydrocarbons (PAHs) [START_REF] Kennicutt | Star Formation in the Milky Way and Nearby Galaxies[END_REF], and the contribution from starlight that dominates the bulk of emission at high frequencies, specially in the UV. These last two contibutions are tightly coupled; dust grains absorbs preferently the UV photons and reradiate the energy at IR wavelengths. Near λ = 912 Å there is a sharp cut in the distribution of UV photons, known as the Lyman limit, as a result of photon absorption by neutral atomic hydrogen, which turns out to be the most abundant species in the ISM. Finally, a high-energy diffuse radiation, composed mainly of X-rays from hot plasmas produced by supernova explosions (T ∼ 10 5 -10 6 K), have been detected by the ROSAT satellite [START_REF] Snowden | Analysis procedures for ROSAT XRT/PSPC observations of extended objects and the diffuse background[END_REF][START_REF] Snowden | First Maps of the Soft X-Ray Diffuse Background from the ROSAT XRT/PSPC All-Sky Survey[END_REF].

The Far-Ultraviolet radiation

The far-ultraviolet (FUV) radiation is of great importance for the physics of the ISM. FUV photons paticipate actively to photoionization and to photodissociation processes, and drive the "warm chemistry" (Le [START_REF] Bourlot | Infrared and submillimetric emission lines from the envelopes of dark clouds[END_REF]Levrier et al. 2012).

It is thought that the ultraviolet background corresponds predominantly to the starlight from O and B stars scattered by the interstellar dust grains of high albedo. One of the first estimates of the intensity of the UV field was done by Habing (1968), who suggested that this radiation field would be approximately constant thoughout space. In the same work he proposed a typical value of u λ ≈ 40 × 10 -18 erg cm -3 Å -1 at wavelength λ = 1000 Å, giving a typical intensity of νu ν ≈ 4 × 10 -14 erg cm -3 . This value is used as a reference for the UV field, which is expresses as a function of this canonical value. Another Figure 1.13: Ultraviolet background in the solar neighborhood. Estimates from Habing (1968); Draine (1978) and [START_REF] Mathis | Interstellar radiation field and dust temperatures in the diffuse interstellar matter and in giant molecular clouds[END_REF] and observations from [START_REF] Henry | Far-ultraviolet studies. VII -The spectrum and latitude dependence of the local interstellar radiation field[END_REF] and [START_REF] Gondhalekar | Observations of the interstellar ultraviolet radiation field from the S2/68 sky-survey telescope[END_REF]. Figure from [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF].

parameter that is commonly used to describe the radiation field is G 0 G 0 = 13.6 eV 6 eV u ν dν 5.29 × 14 -14 erg cm -3 (1.21) which is the overal energy density integrated between 6 and 13.6 eV with respect to Habing's estimate integrated in the same region u Hab(6-13.5 eV) =

13.6 eV 6 eV u ν dν = 5.29 × 14 -14 erg cm -3 . Both values, (νu ν ) 1000A and u Hab (6 -13.5eV) are usually referred as Habing units.

A decade later Draine (1978) recalculated the intensity of the UV background using previous observations (Habing 1968;[START_REF] Henry | Far-ultraviolet studies. II -Galactic-latitude dependence of the 1530 A interstellar radiation field[END_REF][START_REF] Gondhalekar | The interstellar radiation field between 912 A and 2740 A[END_REF]Jura 1974;[START_REF] Hayakawa | Diffuse Component of the Cosmic Far UV Radiation and Interstellar Dust Grains[END_REF], obtaining a simple analytic expression for the angle-averaged photon flux F (E) = λ 3 u λ /(4πh 2 c) valid in the range between 5 and 13.6 eV, that reads

F (E) = [1.658×10 6 E-2.152×10 5 E 2 +6.919×10 3 E 3 ] photons cm -2 s -1 sr -1 eV -1
(1.22) where E is given in eV. This function is known as the "Standard UV field" or "Draine field", and defines a G 0 parameter of approximately 1.7 in Habing units. Fig. 1.13 shows energy density for the Draine field, among other estimates as [START_REF] Mathis | Interstellar radiation field and dust temperatures in the diffuse interstellar matter and in giant molecular clouds[END_REF], who gives a piecewise description of the UV background (MMP83), and lower limit estimates calculated by [START_REF] Gondhalekar | Observations of the interstellar ultraviolet radiation field from the S2/68 sky-survey telescope[END_REF] who measured direct starlight (GPW80).

Pressure equilibrium and equipartition of energy

Processes of different nature contribute to the global pressure in the ISM. The main contributions are

Thermal pressure In the ISM, the WNM phase and the CNM phases are roughly at thermal pressure equilibrium, with typical values (respectively for each phase) of n ∼ 10 cm -1 and T ∼ 100 K, and n ∼ 0.6 cm -1 and T ∼ 5000 K for the density and temperature, resulting in a typical thermal pressure of

P th = nkT ∼ 4.1 × 10 -13 erg cm -3 (1.23)
Magnetic pressure Magnetic fields pervade the space, and provide a magnetic support, or magnetic pressure. In the ISM the typical strength of the magnetic field is around 5 µG

P mag = B 2 8π B 5 µG 2 × 10 -12 erg cm -3 (1.24)
Turbulent pressure Macroscopic gas motions provide and additional pressure, also known as ram pressure or dynamic pressure.

P ram = ρv 2 n 1 cm -3 v 10 km s -1 2 × 10 -12 erg cm -3
(1.25)

Cosmic rays pressure relativistic particles, specially protons, with energies of a few MeV. The typical energy density is u CR = 0.8 eV cm -3 .

P CR = 1 3 u CR u CR 0.8 eV cm -3 × 4.3 × 10 -13 erg cm -3 (1.26)
Radiation pressure from stars and the cosmic microwave background radiation (CMB). The typical energy densities for these two sources are u stars = 0.5 eV cm -3 , and u CM B = 0.26 eV cm -3 (for T = 2.725 K), respectivey. The total radiation pressure is about

P rad ∼ 4 × 10 -13 erg cm -3 (1.27)
For the standard conditions that prevail in the ISM, all these pressures are comparable, ∼ 10 12 erg cm -3 , within an order of magnitude. This reflects a kind of equipartition of energy in the ISM, that is not yet well understood. In any case, what we could say is that none of the physical processes can be neglected, and that any description or model of the ISM must take them all into account.

2

The Molecular ISM In 1926, during Arthur Eddington's Bakerian Lecture, the question whether interstellar clouds would contain molecules opened a debate about nature of the composition of the ISM. Eddington himself stated that no processes were known that would form molecules efficiently in those extreme conditions, and that the presence of molecules in the interstellar gas would be very unlikely. During the late 1930s and early 1940s the first molecules in the ISM imposed a new paradigm. CH, CH + , and CN, detected respectively in 1937CN, detected respectively in , 1940CN, detected respectively in , and 1941, were observed as absorption bands from electronic transitions toward slightly reddened stars in optical wavelengths. But it was during the late 1960s, with the development of radioastronomy, and more recently with infrared telescopes, that the list of molecules found in the ISM increased considerably. Since then some 125 gas-phase molecules have been confirmed. Most of those molecules are organic, and their complexity ranges from diatomic molecules, as H 2 and CO, to very complex molecules containing more than 13 atoms, such as the extreme case of C 60 . Data bases, such as the Cologne Data base for Molecular Spectroscopy (CDMS http://www.astro.uni-koeln.de/cdms/molecules) or The Astrochymist ( http://astrochymist.org), show the variety of molecules that have been found in the space. This molecules exist in very different environments, so that they represent a probe for the physical conditons.

Among all these molecules molecular hydrogen is the simplest neutral molecule and the most abundant in the interstellar medium. Molecular hydrogen is found in a variety of astrophysical environments where the gas is well shielded against the surrounding UV field, typically where the visual extinction A V 0.01 -0.1 mag. H 2 is the first neutral molecule to be formed, and an important coolant for the formation of the first stars. H 2 can play a role in the thermal balance of the ISM and it intervenes in several reactions regulating the formation of other molecules.

Molecular hydrogen

In spite of being ubiquitous, molecular hydrogen is very hard to observe. H 2 is a symmetric and homonuclear diatomic molecule, linked by a covalent bound. As both atoms are identical, the center of mass coincides with the center of charge, producing a molecule that lacks of permanent dipole. Dipolar transitions are forbidden and only very weak electric-quadrupole transitions are allowed. Direct observation of H 2 is very challenging. Electronic transitions occur at energies between 11 and 14 eV, lying in the UV region, in the Lyman and Werner bands. Rotational-vibrational (rovibrational) and rotational transitions lie in the mid-IR region, for which the atmosphere is partially transparent, but due to their quadrupolar nature they are very hard to be detected. Furthermore, H 2 does not posses fine structrure nor hyperfine structure splitting in its ground state, so it cannot be detected by radio emision. Since the first laboratoy observation of H 2 quadrupole transitions made by [START_REF] Herzberg | Quadrupole Rotation-Vibration Spectrum of the Hydrogen Molecule[END_REF], it was necessary to wait until the advent of space telescopes to be able to observe H 2 in space. Several observations showed significant underabundances of neutral atomic hydrogen in dark interstellar clouds, which seemed to point out that the missing hydrogen would be in molecular form [START_REF] Solomon | Molecular and Solid Hydrogen in Dense Interstellar Clouds[END_REF]. [START_REF] Carruthers | Rocket Observation of Interstellar Molecular Hydrogen[END_REF], using rocket-based observations, detected H 2 in absorption toward ξ Persei, in the far-UV spectrum. This and other detections [START_REF] Smith | Interstellar Molecular Hydrogen Observed in the Ultraviolet Spectrum of Delta Scorpii[END_REF][START_REF] Spitzer | Spectrophotometric Results from the Copernicus Satellite.IV. Molecular Hydrogen in Interstellar Space[END_REF]) confirmed that large amounts of H 2 are present in the ISM, and were in good agreement with molecular fractions estimates (Hollenbach et al. 1971). Since then, several space-based missions have been launched in order to observe H 2 in the UV (Copernicus, ORPHEUS, FUSE), and in the mid-IR (ISO, Spitzer).

H 2 structure

Diatomic molecules can be modelled as two particles linked by a spring, that allows rotational and stretching motions. In addition to electronic levels, molecules possess energy levels related to the kinetic energy of nuclei. Diatomic molecules possess vibrational (stretching) modes, and rotational modes (Fig. 2.2). For a nonrotating diatomic molecule (J = 0), the vibrational energy, of vibrational quantum number v, can be approximated by an harmonic oscillator, of energy

E vib = ω e v + 1 2 (2.1
)

ω e = k µ (2.2) where µ = m 1 m 2 /(m 1 + m 2 )
is the reduced mass of the two nuclei, and the constant k = d 2 E potential /dR 2 depends on the form of the potential well.

Pure rotational modes are approximated by a rigid rotator, where the separation between the two atoms is constant and equal to the equilibrium distance R e . In this case the rotational energy, for a quantum rotational number J, is given by

E rot = B R J(J + 1) (2.3)
where I = µR e is the moment of inertia, and B R = 2 /(2µR 2 e ) is the rotational constant.

Ortho and Para H 2

Each hydrogen has a nuclear spin of I = 1/2. Both nuclei can combine in two different states when forming a molecule, generating a singlet (antisymmetric) and a triplet (symmetric) nuclear functions Ψ n . The total wave function Ψ = Ψ n Ψ rot , where Ψ n is the nuclear function, and Ψ rot is the rotational function, is constrained to be antisymmetrical. For example, if we analyze the ground electronic state of the hydrogen molecule, 1 Σ + g , where both electrons are in the lowest σ g orbital, which is a symmetric state and the ground vibrational state v = 0, wich is also a symmetric state, we can conclude that the rotational function must necessarily be antisymmetrical. In the same manner, as the total wave function must satisfy the antisymmetry, the rotational and the nuclear functions must fulfill the antisymmentry of Ψ n Ψ rot , leading to two different forms of H 2 , known as ortho and para H 2 Ortho-H 2 when the nuclear function is symmetric (parallel nuclear spins), the rotational function must be antisymmetric, and thus the quantum rotational number J must be odd. The statistical weight is

g J = 3(2J + 1)
Para-H 2 in this case the nuclear function is antisymmetric (antiparallel nuclear spins), and Ψ rot must be symmetric, forcing J to be even, with a statistical weight g J = 2J + 1.

As the nuclear spin does not change, rovibrational transitions must keep the same parity, an only transitions with ∆J = 0 and ∆J = ±2 are allowed, preventing the change between ortho and para states, which behave like different species. H 2 gas is formed in a mixture of both states, ortho and para hydrogen. Conversion between both forms can occur only in reactive collisions, and it happens at a very slow rate, so the ortho-to-para ratio (OPR) carries information about the physical conditions of H 2 formation.

At thermal equilibrium, the ortho-to-para ratio is a function of temperature. At low temperatures, the para-H 2 state dominates, concentrating more than 99% of molecules for temperatures lower than 20 K. For temperatures higher than room temperature (∼ 300 K), the ortho-to-para ratio reaches its maximum possible value of 3 (75% of ortho-H 2 , and 25% of para-H 2 ), which reflects the statistical weight of both species. The OPR for temperatures lower than room temperature can be estimated by supposing a thermalized Boltzmann population

n J = g J exp (-E vJ /kT ) = g J exp (-B R J(J + 1)) (2.4)
where the rotational constant can be calculated using

(E J=1 -E J=0 ) = 2B R = 2 /IR e ≈ 175 K.
As the most populated levels of H 2 are those of lowest energy, the population of ortho and para H2 can be approximated by the first lowest levels for each state

n O n P ≈ n(J = 1) + n(J = 3) n(J = 0) + n(J = 2) = 9 exp (-2B R /T ) + 21 exp (-12B R /T ) 1 + 5 exp (-6B R /T ) (2.5)
This expression is valid only for low temperature gas. Once the OPR has reached the value of 3, it must be kept at this value. The ortho-to-para ratio estimated in this way is depicted in Fig. 2.3.

Figure 2.3: Ortho-to-para ratio as a function of temperature for thermal equilibrium.

H 2 formation

The formation of molecular hydrogen remained a puzzling problem for many years. Nowadays we know that the most efficient H 2 formation process is through grain surface catalysis. In this section we will review the main formation mechanisms in gas phase as in grain surface.

Gas-phase Formation

The first mechanism that has been invoked is a two-body reaction called direct radiative attachment. This mechanism consists in a gas-phase reaction between two free particles, where the excess of energy is carried out by a photon.

H + H → H 2 + hν (2.6)
Figure 2.1 shows what happens when two hydrogen atoms, in their ground electronic state (1s), approaches each other to form a H 2 molecule. The molecule is formed in a repulsive state ( 3 Σ + u ). As the molecule lacks of electric dipole, it cannot remove the excess of energy by emitting a photon, thus the molecule is rapidly dissociated.

Another mechanism that could be invoked is the three-body association [START_REF] Palla | Primordial star formation -The role of molecular hydrogen[END_REF], where the excess of energy is carried out by a third atom. The main reactions are

H + H + H → H 2 + H + KE (2.7) H + H + H 2 → H 2 + H 2 + KE (2.8) H + H + He → H 2 + He + KE (2.9)
These reactions can occur at very high densities, when the gas reach number densities of 10 8 -10 9 cm -3 , the first reaction being dominant. But under the conditions prevailing in the ISM, a three body encounter is extremely rare, and a three-body reaction can be neglected.

A more plausible mechanism in gas-phase is the associative detachment [START_REF] Mcdowell | On the formation of H2 in H I regions[END_REF], that is to say a radiative association between H -and an electron,

H + e -→ H -+ hν (2.10)
with a reaction rate coefficient k = 3.37×10 -16 (T /300) 0.64 exp(-9.2/T ) cm 3 s -1 , followed by H 2 formation through associative detachment

H -+ H → H 2 (v, J) + e -+ KE (2.11)
with a reaction rate coefficient k = 4.82 × 10 -9 (T /300) 0.02 exp(-4.3/T ) cm 3 s -1 for temperatures T = 10 -100 K, and for higher temperatures T = 101 -3000 K, k = 4.32 × 10 -9 (T /300) -0.39 exp(-39.4/T ) cm 3 s -1 [START_REF] Mcelroy | The UMIST database for astrochemistry[END_REF]. But, independently of the rate coefficient, the formation rate of H 2 through this formation path is proportional to the number density of H -. In the ISM, specially in dense gas, the ionization fraction is very low, only about x H + = x e -= 10 -4 -10 -7 [START_REF] Tiné | Observational Indicators of Formation Excitation of H_2[END_REF], and in diffuse gas H -can be destroyed in several ways. It can be neutralized either by H + [START_REF] Moseley | H + + H -Mutual Neutralization Cross Section Obtained with Superimposed Beams[END_REF])

H -+ H + → H + H (2.12)
or by other positive ions [START_REF] Dalgarno | Heating and Ionization of HI Regions[END_REF])

H -+ M + → H + M (2.13)
where M stands for metals. In both cases the reaction rate is k = 7.51 × 10 -8 (T /300) -0.5 cm 3 s -1 . But in standard conditions for diffuse gas, H -is mainly photodestroyed by the interstellar ultraviolet radiation field (ISRF) at rate k = 1.43 × 10 -7 exp(-0.5A V ) s -1 , where A V is the local visual extinction

H -+ hν → H + e - (2.14)
Even though formation rates for this mechanism are about an order of magnitude faster than the direct radiative attachment mechanism, they are still very inefficient, specially at low temperatures, and it cannot explain the observed abundances of H 2 , then a different mechanism is needed.

Grain-Surface Catalysis of H 2

Nowadays it is widely accepted that molecular hydrogen is formed on dust grain surfaces. This mechanism was first suggested by van de Hulst, and first described by [START_REF] Mccrea | The formation of Population I stars, II. The formation of molecular hydrogen in interstellar matter[END_REF]. It was known that atoms were able to stick on the surface of solids exposed to atomic gas, and it was known that the adsorbed atoms could combine with atoms already present on the surface to form molecules.

The existence of dust were known since the works of Trumpler (1930a,b,c), who observed systematic variations in the brightness of stars clusters, that he interpreted as an extinction of starlight due to the presence of small solid particles, or fine cosmic dust. Now we know that dust represents only about 1% of the total mass of the ISM, but it plays a key role in astrophysical contexts. It shields the gas against UV radiation, and acts as a catalyst for several reactions, H 2 formation being the most relevant (Gould and Salpeter 1963;Hollenbach and Salpeter 1971).

H 2 can be formed on grain surfaces in two pathways: the Eley-Rideal (ER) mechanism and the Langmuir-Hinshelwood (LH) mechanism. Fig. 2.4 shows the two mechanisms. A detailed study of H 2 formation on dust grains, paying particular attention to these two pathways to form molecular hydrogen, has been recently preformed by [START_REF] Bron | Stochastic processes in the interstellar medium[END_REF] in his thesis work. In the case of ER mechanism, a chemisorbed atom reacts with a gas phase hydrogen atom. The adsorbed atom is chemically bonded to the grain surface, and it is not able to move. The second atom, still in gas-phase, hits the adsorbed atom to form an H 2 molecule, liberating some 4.5 eV. A part of the released energy is used to desorb the molecule.

The LH mechanism differs in the nature of the bond between hydrogen atoms and the grain surface. In this case the atom is physisorbed, and the bond is much weaker. This allows the molecule to wander over the grain surface until it encounters another adsorbed atom to form a new molecule. As in the ER mechanism, the molecule is desorbed by using part of the energy released by the reaction.

These mechanisms are essentially different from the mechanisms operating in gas-phase. In the case of H 2 formation on grain surfaces the excess energy is evacuated in a non-radiative manner.

Estimate of H 2 formation rate

The formation rate of H 2 can be estimated using the collision rate between hydrogen atoms and grains, considering that two collisions are needed to form one molecule. This can be written in the following manner:

dn(H 2 ) dt 1 2 n(HI)n gr Σ gr v th (2.15)
where n(HI) is the density of atomic hydrogen in gas phase, n gr is the density of grains, v th = (8kT /πm H ) 1/2 is the thermal velocity (where k and m H are the Boltzmann constant and the mass of the hydrogen atom) and Σ gr is the mean surface of grains. n gr Σ gr represents the total grain surface available to adsorb atoms, that can be inferred from the observed extinction, that relates the dust surface distribution to the number density of H nucleons (for a detailed treatment of the grain distribution see [START_REF] Weingartner | Dust Grain-Size Distributions and Extinction in the Milky Way, Large Magellanic Cloud, and Small Magellanic Cloud[END_REF])

n gr Σ gr = da dn da πa 2 = n tot Σ -21 × 10 -21 cm 2 (2.16
)

n tot = n HI + 2n H2
is the number density of total hydrogen. The value of Σ -21 will depend on the grain population, but is of the order of 1. Then the formation rate can be written:

dn(H 2 ) dt = n(HI)n tot k form (2.17)
where k form = 0.5 × 10 -21 (8kT /πm H ) 1/2 Σ -21 is the effective rate coefficient. The value of this expression can be estimated numerically as

k form = 7.25 × 10 -17 T 100 1/2 Σ -21 [cm 3 s -1 ] (2.18)
This estimate is in good agreement with the value found by Jura (1974[START_REF] Jura | Interstellar clouds containing optically thin H2[END_REF]) using Copernicus observations in the far UV. Assuming that formation and destruction rates of H 2 were equal, [START_REF] Jura | Interstellar clouds containing optically thin H2[END_REF] inferred a mean rate coefficient

k form,0 = 3 × 10 -17 cm 3 s -1 .
(2.19)

This value was confirmed by FUSE (Gry et al. 2002).

Sticking factor

Equation 2.17 corresponds to the case where all the hydrogen atoms that collides with a grain are adsorbed and lead to the formation of a molecule. Hollenbach and Salpeter (1971) included a sticking coefficient in order to take into account the fact that only a fraction of the incident atoms are adsorbed by the grain surface. Then a simple expression to describe the formation rate of H 2 is given by

k form = k form,0 T 100 K 1/2 × S(T ) (2.20)
This sticking coefficient S(T ) is estimated empirically, and it represents the probability of adsorption and mobility of the hydrogen atom to finally form a molecule. It is important to notice that this coefficient reflects the nature of the forming process. In our treatment we will consider a kind of mean effect, this is to say, any of the two mechanisms will be treated explicitly. Therefore, we would like to have a coefficient varying between 1 for low gas temperatures, where gas atoms can stick more easily, to 0 for very high gas temperature, where atoms in gas phase are more likely to bound from the surface instead of sticking.

A satisfactory expression to describe the sticking coefficient was derived by Le Bourlot et al. (2012) for the LH process, which is given by

S(T ) = 1 1 + T T2 β (2.21)
where S(T ) reflects the ability of the colliding H atom to dissipate the kinetic energy before bounding to form a H 2 molecule. Bron et al. ( 2014) have used T 2 = 464 K and β = 1.5 in order to match the expression of [START_REF] Sternberg | Chemistry in Dense Photon-dominated Regions[END_REF]. Figure 2.5 shows the effect of the dependence in temperature and the sticking factor. The dependence on T /100 comes from the thermal velocity in Eq. 2.15, manifesting the fact that more collisions happens as the temperature increases. The sticking coefficient, on the other hand, reflects the difficulty that experience incident atoms to be adsorbed. The global effect shows that the efficiency is maximum at some T = 300 K.

H 2 destruction

Molecular hydrogen can be destroyed collisionally at temperatures T > 4000 K [START_REF] Dalgarno | Collision-induced dissociation of interstellar molecules[END_REF], and by UV photons [START_REF] Black | Interstellar H2 -The population of excited rotational states and the infrared response to ultraviolet radiation[END_REF]). Collisional dissociation is not a major destruction process in the ISM due to the very low densities associated to hot gas. On the other hand, photodissociation is the principal process that destroys H 2 under the conditions that prevail in the ISM.

Direct Photodissociation

The simplest mechanism to destroy H 2 is the direct photodissociation. It is important to notice that in the same way as H 2 cannot be formed in gas phase by the emission of a photon, it cannot be destroyed by a transition between its ground state X1 Σ + g and a repulsive electronic state with energy lower than 13.6 eV. Thus the transition must be towards the next excited state able to dissociate the molecule, through the following transition

H 2 (X 1 Σ + g ) + hν → H 2 (B 1 Σ + g ) → H(1s) + H(2s, 2p) (2.22)
In this case a H 2 molecule in its electronic ground state X 1 Σ + g , where both electrons are in the lowest σ g orbital, absorbs an UV photon of energy higher than hν > 14.7 eV 1 , able to dissociate the molecule. Direct photodissociation is not expected to be the main photodestruction process, due to other processes with larger cross-sections that are in competition for these photons. Photons with energies hν > 13.6 eV more likely ionize HI, while photons with energies hν > 15.5 eV will be prone to photoionize H 2 .

UV Fluorescent Photodissociation

The main mechanism to destroy H 2 is a two step process called UV fluorescent photodissociation, that involve excited states of H 2 . The first step is the absorption of a UV photon in the Lyman or Werner bands, which excites the molecule from its ground electronic state X 1 Σ + g (Fig. 2.6). A photon in the Lyman band (∆E > 11.2 eV, λ < 1108Å) is absorbed preferently if its direction is parallel to the internuclear axis [START_REF] Malcolm | Polarisation of radiation from H2 excited by electron impact -Lyman and Werner emissions and Lymanalpha[END_REF]). This excites one of the electrons from the ground level σ g orbital to the lowest σ u orbital, that corresponds to the excited state

B 1 Σ + u of the H 2 molecule. H 2 (X 1 Σ + g , (v, J)) + hν → H 2 (B 1 Σ + u , (v , J )) (2.23) H 2 (B 1 Σ + u , (v , J ) → H 2 (X 1 Σ + g , (v , J )) + hν (2.24)
On the other hand, a photon in the Werner band (∆E > 12.3 eV, λ < 1008Å), absorbed preferently if its direction is parallel to the internuclear axis, will excite one of the electrons from the σ g orbital to the lowest Π u orbital, sending the molecule to the excited state

C 1 Π + u . H 2 (X 1 Σ + g (v, J)) + hν → H 2 (C 1 Π u , (v , J )) (2.25) H 2 (C 1 Π u , (v , J )) → H 2 (X 1 Σ + g , (v , J )) + hν (2.26)
In both cases, the electronically excited molecule will decay into the electronic ground state, in a different rovibrational state (v , J ), emitting an UV photon. In 90% of cases the molecule cascades toward bound rovibrational levels of the electronic ground state, populating the excited rovibrational levels, and only in 10% of the cases the absorptions are followed by transitions toward unbound vibrational levels (v ≥ 14) of the ground state, known as the vibrational continuum, dissociating the molecule.

The photodestruction rate of H 2 is

dn(H 2 ) dt = -n(H 2 )k ph (2.27)
where k ph is the rate coefficient for this process,

k ph = 4πJ ν u (σ ν,u p diss,ν ) dν, (2.28)
The rate coefficient depends on the intensity of the UV field, 4πJ ν , that sets the number of UV colliding photons, on the probability of photoexcitation toward all upper levels u, and on the probability of spontaneous decay, p diss,ν , for each excited level toward the vibrational continuum.

Photodissociation rate

The intensity of the UV field and the probability of photoexcitation and spontaneous decay can be approximated by an unshielded field attenuated by shielding factors. Therefore the mean rate coefficient can be described in a simpler form

k ph = e -τ d f shield k ph,0 , (2.29)
where e -τ d is a shielding factor due to dust, f shield is the self-shielding coefficient that takes into account the optical thickness of gas for the UV, and k ph,0 is the photodissociation rate coefficient for optically thin gas given by the following expression (Draine and Bertoldi 1996;Glover and Mac Low 2007)

k ph,0 = 3.3 × 10 -11 G 0 s -1 ,
(2.30) here G 0 is the intensity of the UV field in Habing units (Habing 1968). For the ISRF, or Draine field (Draine 1978), this value is G 0 = 1.7.

The shielding coefficients in Eq. 2.29 take into account the attenuation of radiation due to the presence of matter in the line of sight, but they operate in different ways.

Shielding by dust

The first shielding effect corresponds to the UV extinction due to dust absorption. Dust grains absorb and scatter photons at all wavelengths, but this absorption is more efficient at short wavelengths, particularly in the UV. Thus dust extinction removes a part of photodissociating photons, reducing the photodissociation rate. The optical depth for dust is given by

τ d = r 0 σ d n tot dr = σ d × N (2.31)
where σ d is the effective cross section, and N is the total column density. As dust extinction does not vary strongly in the range between 912 -1110 Å, Draine and Bertoldi (1996) approximate the effective cross section, by its value at λ = 1000 Å, σ d,1000 = 2 × 10 -21 cm -2 .

Self-shielding

When H 2 column densities become high enough, the photoexcitation transitions from the Lyman and Werner systems become optically thick, and H 2 molecules are shielded against photodissociation by other H 2 molecules, that interacts with radiation absorbing photons at specific wavelengths. The shielding function depends on the amount of molecular hydrogen that a UV photon must traverse before it can dissociate a molecule, and must be calculated in every direction.

For a single directions this can be approximated by the following expression (Draine and Bertoldi 1996)

f shield = 0.965 (1 + x/b 5 ) 2 + 0.035 (1 + x) 1/2 exp(-8.5 × 10 -4 (1 + x) 1/2 ), (2.32) x = N (H 2 ) 5 × 10 14 , (2.33) b 5 = b 10 5 cm s -1
(2.34)

here N (H 2 ) stands for the total column density of molecular hydrogen, and b = b 2 th + b 2 turb stands for the Doppler broadening parameter. This parameter represents the fact that absorptions do not occur at rest, due to thermal and turbulent gas motions.

Both shieldings operate in a line-by-line basis, so it is needed to compute them in all directions.

Thermal feedback from H 2

The presence of H 2 molecules, and its formation and destruction, modifies the thermal balance. H 2 offers new cooling channels, while H 2 formation and photodissociation processes heat the gas2 .

Heating

Heating by H 2 formation H 2 formation is a highly exothermic reaction. When a H 2 molecule forms, some ∆E = 4.5 eV are released, and a fraction of this energy contributes to heat the gas. The released energy is expected to be shared among all the degrees of freedom of the molecule and the grain. In other words, after using a part of this energy to desorb the molecule, the remnant energy is distributed into translational kinetic energy of the newly formed molecule, into H 2 internal energy, and into grain heating, but how this energy is distributed is not well known. Whether H 2 formation contributes to heat the grain, or which fraction of the binding energy actually heats the gas, or even in which state molecules are created are still open issues. Recent experimental works [START_REF] Islam | Formation Pumping of Molecular Hydrogen in Dark Clouds[END_REF][START_REF] Lemaire | Competing Mechanisms of Molecular Hydrogen Formation in Conditions Relevant to the Interstellar Medium[END_REF]) have shown that H 2 molecules leave the surface of the grain in rovibrational excited states.

The heating rate due to H 2 formation can be written in a general manner as

Γ H2form = form 7.2 × 10 -12 R H2 ergs s -1 cm -3 (2.35)
where the total available energy is ∆E = 4.5 eV = 7.2 × 10 -12 erg, and R H2 = k form n(HI)n tot is the H 2 formation rate. With respect to the fraction of energy that contributes to heating, form , there is no a global consensus. Glover and Mac Low (2007) have used a variable fraction,

form = 1 + n cr n -1
(2.36) that depends on the density n, and on the critical density n cr at which collisional de-excitation and radiative de-excitation happen at the same rate. In this case, for densities higher than the critical density most of the binding energy heats the gas. A more common approximation, that agrees better with the evidence of excited internal levels of nascent formed hydrogen molecules, supposes an equipartition of the binding energy. In this case form = 1/3

(2.37) leading to 1.5 eV that is converted into kinetic energy, and then heats the gas. A more detailed balance has been calculated for the ER mechanism by [START_REF] Sizun | Investigation of ZPE and temperature effects on the Eley-Rideal recombination of hydrogen atoms on graphene using a multidimensional graphene-H-H potential[END_REF]. Using two collision energies (E = 15 meV and E = 50 meV), for a binding potential of Bachellerie of the order of 0.1 eV, has found that most of the energy released by H 2 excites the internal rovibrational levels (E int ∼ 2.7 eV), about 1 eV heats the grain surface, and about 0.6 eV is converted into kinetic translational energy that contributes to heat the gas.

Heating by H 2 destruction

The photodissociation of H2 can also contribute to heat the gas. When a H 2 molecule is hit by a photodissociating UV photon, fluorescent photodissociation can follow in 10% of the cases. A fraction of the energy of the UV photon is converted into kinetic energy carried by the pair of dissociated hydrogen atoms. Black and Dalgarno (1977) have estimated that the amount of energy that contributes to heating is about 0.4 eV per photodissociation event, for a typical UV field. Following this result, Glover and Mac Low (2007) have estimated the heating rate by photodissociation of H 2 as

Γ H2dest = 6.4 × 10 -13 k ph n(H 2 ) ergs s -1 cm -3 (2.38)
where k ph is given in equation 2.29. Molecular hydrogen is the main coolant in primordial gas, when metal-line cooling is not possible, but when metallicities reach values higher than 0.1% of solar metallicity, metal-line cooling dominates the thermal balance [START_REF] Loeb | The First Galaxies in the Universe[END_REF][START_REF] Omukai | Thermal and Fragmentation Properties of Star-forming Clouds in Low-Metallicity Environments[END_REF][START_REF] Omukai | Low-metallicity Star Formation: Prestellar Collapse and Protostellar Accretion in the Spherical Symmetry[END_REF]. It has been shown that in some cases H 2 can play an important role as coolant. Gnedin and Kravtsov (2011) have shown that in galaxy formation simulations, cooling by H 2 can become dominant when temperatures are T 5000 K. Glover and Clark (2014), using a one-zone model, have shown how molecular hydrogen, under certain conditions, can contribute in great manner to the cooling in diffuse gas. In particular they have found that a necessary condition is that metallicities must be low enough, of about Z 0.1Z , and that the ratio of the strength of the UV radiation field to the local gas density, G 0 /n ∼ 0.01, must also be low, in such a way that H 2 is present in enough quantity to contribute to cooling. This last criterion implies that for the standard UV field, where G 0 = 1, cooling by H 2 is efficient for regions where the total gas density is about n = 100 cm -3 . These conditions are those of gas transiting from the WNM to the CNM, so H 2 could potentially play an important role in this transition for low metallicity gas.

Cooling

The cooling mechanism is through collisionally excited lines. As most of molecular hydrogen is found in the ground electronic level, the main contribution to cooling comes from rovibrational transitions within the electronic ground state. Figure 2.7 shows the internal energy levels of the ground electronic state of H 2 molecule. These excited rovibrational levels are populated collisionally, and de-excitation follows through quadrupole transitions. Pure rotation transitions of H 2 fall in the Mid-IR region, while rovibrational transitions occur in the Near-IR. In both cases photons can leave the gas, contributing to cool the gas.

Rovibrational transitions use to happen between close vibrational levels, and the probability of transitions for ∆v ≥ 2 is much weaker. Le Bourlot et al. (1999) have computed an accurate cooling function for H 2 , where they include transitions between all rovibrational levels of H 2 with energies E vJ ≤ 20000 K,

(v = 0, J ≤ 16; v = 1, J ≤ 13; v = 2, J ≤ 10; v = 3, J ≤ 8,
W (H 2 ) = 1 n(H 2 ) vJ,v J (E vJ -E v J ) n vJ A(vJ → v J ),
(2.39)

where W (H 2 ) is the amount of energy lost per H 2 molecule, E vJ and E v J are the energy levels, n vJ is the number of molecules in the rovibrational level (v, J), and A(vJ → v J ) is the probability of transition from the upper level of energy E v J toward the level of lower energy E v J . Figure 2.8 shows the cooling function computed for two extreme limits of the ortho-to-para ratio. This figure shows that the cooling is not sensitive to the OPR and it depends principally on the physical conditions of the cloud (total number density of hydrogen atoms, total number density of H 2 molecules, and temperature).

3

The Numerical Tools Numerical simulations are a modern tool that permits us to reproduce the behavior of different systems. They are based on models that contain the processes that are relevant to the problem. Numerical simulations are used in almost all areas of research and development, as fluid mechanics, meteorology, aeronautics, astrophysics, demography, sociology, and many others. They are particularly useful when the subjects are too complex to be solved analytically, or too complex to be recreated in laboratory conditions.

In the domain of astrophysics numerical simulations are of the greatest importance, due to complexity of processes, to the spatial and time scales involved, as well as the physical conditions, that cannot be reproduced in laboratory conditions.

Numerical methods

Astrophysical fluids, and fluids in general, can be described in two different, but equivalent ways.

The Lagrangian description keeps track of the location of individual fluid parcels. Each fluid parcel carries the description of its physical properties, density, velocity, pressure, magnetic field as a function of time:

• ρ i (t) = ρ( x i0 , t) • v i (t) = v( x i0 , t) • P i (t) = P ( x i0 , t) • B i (t) = B( x i0 , t)
where i stands for the i-th fluid parcel. Each variable depends on the initial position of the particle x i0 = x i (t = 0), and the entire flow is described by the aggregation of each individual parcel.

On the other hand, the Eulerian description, instead of follow fluid parcels, describes the variation of the physical properties of the flow in coordinates that are fixed in space, giving a field description of physical properties represented

• ρ = ρ( x, t) • v = v( x, t) • P = P ( x, t) • B = B( x i0 , t)
where each variable represents the respective field of the whole flow.

Numerical simulations are discretized representations of reality, and different numerical methods are adapted for each kind of description of the fluid. They can represent a physical fluid either by replacing the fluid by particles and follow these individual parcels of fluid or by defining a grid and see how the fluid moves through it.

As sand grains can collectivelly behave as a fluid, in the same way numerical simulations can treat fluids as a collection of tiny spherical particles. This first approximation fits better the Lagrangian approach. This method lies on a N-body description of the fluid, where nodes of the computational mesh follow the material particles of the fluid, giving their positions and velocities as a function of time. The main disadvantages of this kind of description are the difficulty to describe transport processes, the large number of particles needed to describe the fluid accurately, and the refinement strategy, that refines only on density. The smoothed particle hydrodynamics (SPH) technique [START_REF] Gingold | Smoothed particle hydrodynamics -Theory and application to non-spherical stars[END_REF] adds to the N-body description a smoothing length for their properties. Several of the widest used astrophysical codes are based on this technique as Seren (Hubber et al. 2011) or Gadget-2 [START_REF] Springel | The cosmological simulation code GADGET-2[END_REF].

The second formulation is the Eulerian approach. This is a patch-based family of codes, where a grid defines different volume controls, fixed to the space, that describe different properties of the fluid as fields, which means that each variable is expressed as a function of the space coordinates and time. The main advantages of this kind of codes are the refinement criteria, that allows a refinement of every field, and that they are well adapted for the description of turbulent motions allowing large deformations. The Sod shock tube test [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]) is used to check the ability of numerical codes to resolve sharp shock interface. Fig. 3.1 shows how AMR codes are better than SPH codes at capturing discontinuities. Some widely used astrophysical codes that uses this strategy are Castro (Almgren et al. 2010), Enzo (The Enzo Collaboration et al. 2013), Flash (Fryxell et al. 2000), Pluto (Mignone et al. 2012), and Ramses (Teyssier 2002;Fromang et al. 2006).

In any case, both prescriptions have to fulfill two criteria in order to give an accurate description of the fluid. Each parcel of fluid (particles or cells) must represent regions that are large in comparison to the mean spacing between molecules, and timesteps must be shorter than the characteristic timescale of the physical problem. The spatial discretization must be small with respect to the scale of the problem. Fig. 3.2 shows two simulations of a galaxy merging, performed with an SPH code and an AMR code, showing how both techniques give results in good agreement with observations. Grid and SPH techniques have been closely compared, finding a general good agreement between both techniques [START_REF] Commerçon | Which Numerical Method for Protostellar Collapse?[END_REF]), but stricking differences in the ability to treat dynamical instabilities and mixing [START_REF] Agertz | Fundamental differences between SPH and grid methods[END_REF]).

The Ramses code

All the simulations performed for this work have been done using the Ramses code. Ramses is an Eulerian code that solves the magnetohydrodynamics (MHD) equations in a grid-based scheme and permits us to follow the evolution of all the physical variables for each cell in the simulation and the exchanges at the interfaces of cells.

Ramses stands for Raffinement Adaptatif de Maillage Sans Effort Sur-humain, or adaptive mesh refinement without a superhuman effort, was created by Romain Teyssier to perform cosmological simulations to study the formation of large scale structures, but the code has evolved, and since then it has been modified in order to enhance its flexibility. The code can be "patched" to adapt it to different purposes and nowadays is used for a wide range of problems, including galaxy formation and interaction, dense core collapse, disk evolution, molecular cloud formation, and recently, chemical evolution.

This code uses the adaptive mesh refinement (AMR) strategy, that allows to refine certain parts of the grid, based on a cell-by-cell basis, using a given criterion. This strategy is particularly adapted for multiscale problems, where several scales are at play simultaneously, as in the case of the turbulent ISM. Furthermore, this code has been parallelized, and it is able to run over multiple CPUs.

The RAMSES code uses a fully threaded tree (FTT) data structure. The grids are divided recursively into octs (groups of eight cells). All the information is local, each cell knows all the information related to the parent cells recursively, making the code suitable for parallel computing by using the MPI library (Message Passing Interface). The hydrodynamical solver is a second order Godunov implementation of the constrained transport (CT) algorithm, able to capture discontinuities at a high precision level. The MHD version was developed by Fromang et al. (2006).

AMR and refinement strategy

The AMR method was introduced for the first time in Berger and Oliger (1984) as an adaptative finite difference method for solving partial differential equations using nested grids in a patch-based AMR, permitting to use cells with different sizes depending on the required resolution. The grid structure in Ramses is tree-based. In the fully threaded tree (FTT) (Khokhlov 1998) the tree is threaded at all levels and the cartesian mesh is refined recursively on a cellby-cell basis, using a refinement criterion, for each resolution level, defined by the user, as the density, mass contained in a cell, Jeans length, gradients of hydrodynamical variables, among the most used. The minimum and maximum levels of resolution, min and max , are also fixed by the user. Figure 3.3 shows how cells are refined. The refinement is done recursively over refinement levels, starting at the coarsest level = 0, that corresponds to the root of the tree and contains the whole volume. All the cells are refined at least until they reach the minimum level min . Cells will be refined recursively while they fulfill the refinement criterion, or until they reach the maximum refinement level, max .

When a cell at level reaches the refinement criterion, it is split into 2 n smaller identical children cells, where n is the dimension of the simulation. For example, in a 2D simulation cells will be split into 4 smaller children cells at level + 1, in this case the structure is called a "quadtree". For 3D simulation the structure is the "octtree", and cells are split into 8 children cells at level + 1, forming new octs. Cells stay refined as long as they satisfy the refinement criterion. Cells can be unrefined when such resolution is no longer needed.

A cell that do not posses internal substructures is called a leaf cell. As the resolution is fixed by the user, there will not be leaf cells with levels higher than max , nor lower than min . Leaf cells have the resolution required for the simulation, therefore only these cells are treated by the solver. Coarser levels contain mean values calculated using all the internal leaf cells.

The FTT assures a full connectivity. Each oct can access directly its parent cell and the 6 neighbors of the parent cell, and it can access its own children octs. Fig. 3.4 shows the full connectivity in the code.

Figure 3.4: Tree-data structure in Ramses. Each oct, at any resolution level points toward its parent cell, to the 6 cells neighbors of its "parent", and it points to its 8 "children" octs.

FTT data structures optimize the access to distant cells. Each cell can access other cells by walking the fully threaded tree recursively, and the numerical cost for this method is logarithmic (Barnes andHut 1986, 1989).

Finally, as the size of the timestep depends on the size of the cell, the time step is also adapted to the resolution level. Between two consecutive levels the length of the timestep doubles when passing to the next coarser level. Ramses solves the MHD equations iteratively, starting with the finest level of resolution max and finishing with the coarsest level demanded by the user, min .

Euler Equations and Godunov methods

Interstellar gas is a highly compressible fluid, and like fluids in general, its dynamics is governed by a set of equations, that reflects the main conservation laws. Fluids with high Reynolds numbers (Re 1) behaves as perfect gases. In the ISM, the Reynolds number is about Re ∼ 10 5 -10 7 , high enough to be considered as a perfect gas. Euler equations are a simplified version of Navier-Stokes equations when viscosity and thermal conduction can be neglected, and are well adapted to grid-based simulations. Euler equations consist in three conservation equations: the continuity equation, the momentum equation, and the energy equation, that express the conservation of mass, momentum and energy. The primitive variables are the gas density ρ, the fuid velocity v, and the pressure P (or equivalently the internal energy ρ ). Equations can be expressed in their conservative form by using combinations of primitive variables.

Hydro case

For an inviscid and non magnetized perfect fluid the conservative equations can be written as follows

∂ρ ∂t + ∇ • (ρ v) = 0, (3.1) ∂ρ v ∂t + ∇ • (ρ v ⊗ v) + ∇P = 0, (3.2) ∂E ∂t + ∇ • [(E + P ) v] = 0 (3.3)
The equation system is closed by the equation of state of ideal gas

P = ρkT µm H = (γ -1)ρ (3.4)
where γ = C p /C v is the ratio of specific heats, µ is the mean molecular weight, and m H is the mass of hydrogen. The total gas energy E in Eq. 3.3 contains the kinetic and the internal energy. Written explicitly this is

E = 1 2 ρ v 2 + ρ = 1 2 ρ v 2 + kT (γ -1)µm H (3.5)
Euler equations in Cartesian coordinates can be written in vector form in the following manner

∂ U ∂t + ∇ • F( U ) = ∂ U ∂t + ∂ F x ( U ) ∂x + ∂ F y ( U ) ∂y + ∂ F z ( U ) ∂z = 0 (3.6)
where U is the vector of conserved quantites, and F( U ) is a flux vector. Explicitly they are:

U =       ρ ρv x ρv y ρv z E       F x =       ρv x ρv 2 x + P ρv x v y ρv x v z (E + P )v x       F y =       ρv y ρv x v y ρv 2 y + P ρv y v z (E + P )v y       F z =       ρv z ρv x v z ρv y v z ρv 2 z + P (E + P )v z       (3.7)
This set of equations can be solved numerically in two different manners: finite difference method or finite volume method. Finite difference methods solves the evolution of the system in a grid-based scheme, where fluid variables and equations are discretized in the grid. On the other hand, volume finite methods also uses a grid, but it uses the volume surrounding the grid node. Volume integrals are converted into surface integrals, that are evaluated as fluxes through the cell walls. The flux at the surface between two adjacent cells must be the same, so the method is conservative.

The Ramses code uses the finite volume method. Physical quantites are volume averaged and have to obey the integral form of the conservation law. The discretization in time and in space is done by the AMR, that defines the sizes of the mesh and timesteps and the volume is that of the cell. For a cell of center (x i , y j , z k ) its boundaries will be at x i±1/2 , y j±1/2 , z k±1/2 . For simplicity, in the following we will analyse the 1D case. Other dimensions can be solved analogously. Then integrating the evolution between the current timestep t n and the next step at t n+1 , in the cell volume, Eq. 3.6 can be written in its integral form in 1D in the following manner

t n+1 t n x i+1/2 x i-1/2 dt dx ∂ U ∂t + ∂ F ( U ) ∂x = 0 (3.8) x i+1/2 x i-1/2 dx t n+1 t n dt ∂ U ∂t + t n+1 t n dt x i+1/2 x i-1/2 dx ∂ F ( U ) ∂x = 0 (3.9) x i+1/2 x i-1/2 dx U n+1 (x) -U n (x) + tn+1 t n dt F i+1/2 (t) -F i-1/2 (t) = 0 (3.10)
where U n (x) = U (t n , x), and F i±1/2 (t) = F (t, x i±1/2 ). In the finite volume approximation we use the volume averaged values of U and the time averaged flux function

U n i = 1 ∆x x i+1/2 x i-1/2 U n (x)dx (3.11) F n+1/2 i±1/2 = 1 ∆t tn+1 tn F i±1/2 (t)dt (3.12)
Then using the averaged values, Eq. 3.10 becomes

U n+1 i = U n i - F n+1/2 i+1/2 -F n+1/2 i-1/2 ∆x ∆t (3.13)
Eq. 3.13 is the exact solution of the problem, because this development supposes a full knowledge of the analytic expression of functions U (t, x) and F (t, x). The Godunov method consists in the interpretation of averages and the incorporation of approximations for the numerical computation of fluxes. Even if fluxes are not exact, the method is conservative by construction.

The original first order Godunov Method [START_REF] Godunov | A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations[END_REF]) was originally developed to solve compressible hydrodynamical equations without adding artificial viscosity to stabilize the scheme. Fluid elements in a simulation are represented by uniform cells at first order, using a piecewise representation of U , but the spatial accuracy can be increased by increasing the order of the description by using either a linear or a parabolic function.

The central problem in the Godunov method is to estimate accurately the Godunov fluxes F

The Riemann Problem

The Riemann problem consists in the determination of the flux exchange of physical quantities at the interface between two cells. The initial value is a piecewise data, with a single discontinuity between the two constant states, where each cell is described by a uniform state U L (left) and U R (right). Mathematically this is

∂ U ∂t + A ∂ U ∂x = 0; U (0, x) = U 0 (x) = U L , if x < 0, U R , if x > 0. (3.14)
where A is the Jacobian. Since no characteristic scale is involved, if U (t, x) is a solution, then U (αt, αx) for any constant α > 0 is also a solution. In other words, the problem is self-similar.

The simplest problem we can analyze is the simple linear advection

∂ U ∂t + a ∂ U ∂x = 0 (3.15)
in this expression a is the characteristic speed, and the solution of the Riemann problem is given by U (x + at), this is

U = U L , if x -at < 0, U R , if x -at > 0. (3.16)
Now if we consider a linear system of m equations and m variables

∂ U ∂t + A ∂ U ∂x = 0 (3.17)
where A is constant squared matrix of m × m elements. If A is diagonalizable, then it can be written as A = K -1 ΛK, where Λ is a diagonal matrix containing the eigenvalues λ i of A, and K is a matrix that contains the eigenvectors of A (AK i = λ i K i ). Defining W = K -1 U we can write the original system in the following manner

∂ W ∂t + Λ ∂ W ∂x = 0 (3.18)
where the solution for this system is given by W i = W i (x-λ i t). Each eigenvalue corresponds to an advection velocity, so the solution for the original variable U can be computed as a collection of advection equations.

The original system is given in Eq. 3.6, that in 1D can be written as

∂ U ∂t + ∂ F ∂x = ∂ U ∂t + ∂ F ∂ U ∂ U ∂x = 0 (3.19)
Introducing the Jacobian of the flux function J( U ) = ∂ F /∂ U , the equation can be written as

∂ U ∂t + J( U ) ∂ U ∂x = 0 (3.20)
To solve the simple hydrodynamical case, we define the independent variables in 1D: U = [ρ, ρv x , E] = [u 1 , u 2 , u 3 ], and we express the flux vector in terms of the independent variables, using Eq. 3.5 and the equation of state given by Eq. 3.4, as follows

F =     u 2 u 2 2 u1 + (γ -1) u 3 -1 2 u 2 2 u1 u 3 + (γ -1) u 3 -1 2 u 2 2 u1 u2 u1     (3.21)
The Jacobian matrix can be written

J =     0 1 0 -1 2 (γ -3) u2 u1 2 (3 -γ) u2 u1 (γ -1) -γ u2u3 u 2 1 + (γ -1) u2 u1 3 γ u3 u1 -3 2 (γ -1) u2 u1 2 γ u2 u1     (3.22)
Eq. 3.20 is the quasi-linear form of the system equation which will be hyperbolic if the eigenvalues of J( U ) are positive. In 1D, in a pure hydrodynamical case, the Jacobian matrix given by Eq. 3.22 possesses three eigenvalues, which are the speeds of waves that propagates in the fluid. The first wave is the rarefaction wave, or expansion wave, where quantities vary continously, the second wave is the contact discontinuity, located between the two cells, and the third one is the shock wave, that corresponds to a discontinuity in the fluid variables where the gas is compressed.

Magnetohydrodynamics

Interstellar gas is also magnetized, so in addition to the hydrodynamical description, it is needed to modify the standard Euler Equations in order to include the effect of magnetic fields. Therefore a full description of the dynamical evolution of astrophysical fluids is given by the equations of the ideal magnetohydrodynamics (MHD) in their eulerian form as follows

∂ρ ∂t + ∇ • (ρ v) = 0, (3.23) ∂ρ v ∂t + ∇ • (ρ v ⊗ v -B ⊗ B) + ∇P = 0, (3.24) ∂E ∂t + ∇ • [(E + P ) v -B( B v)] = 0, (3.25) ∂ B ∂t + ∇ • ( v B -B v) = 0, (3.26) ∇ • B = 0, (3.27)
where B is the magnetic field, and the rest of variables are the same than those of Eqs. 3.1 to 3.3. To the usual Euler equations, two additional equations have been included. Eq. 3.26 is the induction equation, and Eq. 3.27 reflects the absence of magnetic monopoles. This last constraint is known as the constrained transport condition (CT) or solenoidal constraint, and permits to avoid constraint-violating modes that would grow exponentially [START_REF] Meier | Constrained Transport Algorithms for Numerical Relativity. I. Development of a Finite-Difference Scheme[END_REF].

The total energy and total pressure must also be modified

E tot = 1 2 ρ v 2 + ρ + 1 2 B 2
(3.28)

P tot = P + 1 2 B 2 (3.29)
MHD equations can be written in the same manner as Eq. 3.6, where the vector of conserved quantities in 3D is

U = [ρ, ρv x , ρv y , ρv z , B x , B y , B z , E] (3.30)
If we analyze the 1D case, for only F x , we obtain that B x = cst, thus the problem is reduced to a system of 7 equations, where U = [ρ, ρv x , ρv y , ρv z , B y , B z , E], and F x is

F x =           ρv x ρv 2 x + P tot -B 2 x ρv x v y -B x B y ρv x v z -B x B z B y v x -B x v y B z v x -B x v z (E + P tot )v x -B x (v x B x + v y B y + v z B z )           (3.31)
This problem can be written as

∂ U ∂t + ∂ F ∂x = ∂ U ∂t + A ∂ U ∂x = 0 (3.32)
where A is the Jacobian matrix of the system. This matrix possesses 7 eigenvalues λ that corresponds to different wave speeds, namely

c sma = 1 2ρ γP + B 2 - γP + B 2 2 -4γP B 2 x 1/2
(3.34)

c fma = 1 2ρ γP + B 2 + γP + B 2 2 -4γP B 2 x 1/2 (3.35)
Each eigenvalue λ defines a wave that propagates from the discontinuity at speed dx/dt = λ. A usual way to represent these waves is the MHD wave fan shown in Fig. 3.6.

Riemann Solvers

There is no practical exact solver for the MHD Riemann problem, nevertheless there exist approximate solvers with different degrees of diffusivity. Different approaches are possible Exact Riemann solver (nonlinear) [START_REF] Giacomazzo | The exact solution of the Riemann problem in relativistic magnetohydrodynamics[END_REF] a full nonlinear solution of the Rieman problem is possible, but impracticable in numerical codes, due to the extremely expensive numerical cost of these kind of methods. Linearized Riemann solvers or Roe type solver. The Roe solver [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes[END_REF]) supposes a linear system by replacing the A matrix by an approached Roe matrix à that is constant between the two cells, which calculation it is not straightforward, used to calculate the eigenvector and eigenvalues for the seven waves. This method is good at capturing discontinuities and shocks. It has a low numerical diffusivity and viscosity. The main disadvantages are that it may create expansion shocks instead of smooth expansion waves, it can create non physical protrusions in some shocks, and it does not guarantee that fluxes are physical, manifested as negative energies or densities, leading to numerical problems. Other Roe-type solvers have been developed by [START_REF] Komissarov | A Godunov-type scheme for relativistic magnetohydrodynamics[END_REF]; [START_REF] Balsara | Total Variation Diminishing Scheme for Relativistic Magnetohydrodynamics[END_REF]; [START_REF] Koldoba | An approximate Riemann solver for relativistic magnetohydrodynamics[END_REF]; [START_REF] Antón | Relativistic Magnetohydrodynamics: Renormalized Eigenvectors and Full Wave Decomposition Riemann Solver[END_REF].

HLL-type Riemann solver This kind of solver is based on guessed signal speeds. HLL stands for Harten, Lax and van Leer. The first vesion of this method [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF]) uses only the two fastest waves. It assumes there is a uniform state U * , between the left and right states, limited by the two propagating waves λ 1,7 . The HLLC solver (HLL-Contact) adds the contact discontinuity, represented by the entropy wave at λ 4 , to increase the accuracy of the HLL solver. It assumes two intermediate states U * L between λ 1 and λ 4 , and U * R between λ 4 and λ 7 , using jump conditions. This solver is less diffusive than the HLL, but the generalization to the MHD case is not straightforward, due to the incompatibility of jump conditions across fast waves and across the contact discontinuity, due to the insufficient degrees of freedom. In the MHD case Alfvén waves, of speeds given by λ 2 and λ 6 , seem to be necessary to the accurate description of the physical processes at play. The HLLD Riemann solver [START_REF] Miyoshi | A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics[END_REF] adds these waves to the previous solver, using five waves and four intermediate states in total. This last method is less diffusive than the HLL solver. It is quite accurate and robust, being one of the most used solvers in astrophysical codes.

HLL as an exemple of Riemann solver

The Harten-Lax-van Leer Riemann solver [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF]) is based on a two wave approximation, using only the two fastest waves (the magneto-acoustic waves), which are the outermost waves, λ 1 and λ 7 , in the wave fan, as shown in Fig. 3.7.

U =      U L , if x/t < λ L , U * , if λ L < x/t < λ R , U R , if x/t > λ R . (3.36)
This intermediate state U * and the associated flux F * are calculated using the conservation laws. To build the HLL solver we define a control volume limited to the region between ±L, and the area in y -z of surface unity 1. Then integrating U in the control volume at time t * we obtain

L -L U (t * )dx = (L + λ L t * ) U L + (λ R -λ L )t * U * + (L -λ R t * ) U R (3.37)
The integral form in Eq. 3.8 can be integrated in the control volume between t = 0 and t * , with U (t = 0) given in Eq. 3.14

L -L U (t * )dx = L U L + U R -F R -F L t (3.38)
and equating these two last expressions, we obtain an expression for

U * U * = F L -F R -λ L U L + λ R U R λ R -λ L (3.39)
For obtaining an expression for the associated flux F * , we integrate U (t * ) and we use again the conservation law from Eq. 3.8, but this time we integrate in a control volume defined by x = -L and x = 0, and equating the expressions again we obtain

(L + λ L t * ) U L -λ L t * U * = L U L -F * -F L t * (3.40)
and replacing U * from Eq. 3.39 we obtain the following expression for the flux

F * = λ R F L + λ L F R + λ L λ R ( U R -U L ) λ R -λ L (3.41)
If in Eq. 3.40 we had used a control volume defined by x = 0 and x = L, we would have obtained the same expression for the flux F * .

Part II

Methods and Results

57

Tree-based method

Introduction

Radiative transfer plays a central role in the determination of the physical state of gas. The photoelectric heating, photoreactions, photoionization and photodissociation rates are highly dependent on the amount of radiation that reaches each parcel of fluid. Therefore an accurate description of the actual radiation field is of great importance for the correct modeling of the physics of molecular clouds.

In the following, we will refer to the standard UV radiation field, which is the most relevant radiation for the evolution of molecular clouds. In Section 1.6 we have described the standard radiation field, which is one of the most used in numerical studies. This radiation field is thought to be isotropic and constant (Habing 1968), which reduces the difficulty of the description of the field up to a certain extent. Radiation interacts with matter, which absorbs, re-emit, and diffuse light, varying the local properties of the radiation field. The local illumination will be the outcome of the interation of the incoming radiation field and the global distribution of matter. where I v (erg cm -2 s -1 Hz -1 sr -1 ) is the radiation specific intensity at a given frequency ν, and α ν is the absorption coefficient per volume (cm -1). From a microscopical point of view, α ν = nσ ν , where n is the number density of absorbers and σ ν is the effective cross section of absorbers at frequency ν. This equation reflects the loss of intensity of the incoming radiation due to the interaction with matter. The solution to this equation is

Radiative Transfer Equation in an absorbing medium

I ν (s) = I ν (s 0 )e -τν (4.2)
where I ν (s 0 ) is the unattenuated radiation intensity, and τ ν is the optical depth, a dimensionless parameter, given by

τ ν = s s0 α ν ds = s s0 n(s)σ ν ds (4.3)
and it is related to the mean free path of photons s, which is the mean distance that photons can travel before being absorbed, by

τ ν = 1 = s s0 α ν ds (4.4)

Optical Depth estimates in numerical simulations

The optical depth is the key parameter that describes the attenuation of light.

For a correct description of the total incoming radiation it is necessary to integrate optical depths, for every element in the simulation (cells or particles), in any possible direction in the whole espace. As flows evolve, the distribution of matter changes in time, and optical depths must be recalculated at each timestep for a proper description of the influence of the cloud structure on the gas dynamics and physical state.

In numerical simulations, even for simplified descriptions of the incoming radiation field, this is one of the most numerically demanding problems, with an associated computational cost of order N 5/3 , because it involves all the fluid elements in the calculation, and further approximations are needed. Among the most used approximations we can mention the following approximations Optically thin approximation As dust extinction becomes important for column densities above N 1.2 × 10 21 cm -2 , when the gas density is low enough, it can be considered as optically thin and the UV radiation penetrates the cloud without attenuation. As a consequence, temperature is generally overestimated.

Constant attenuation factor This approximation uses a constant flux scaling factor for the UV radiation field. It is in general used for high density gas. In the most extreme cases, as in Goldsmith (2001), when gas is well shielded, photoelectric heating can be completely ignored.

It is hard to tell whether the attenuation factor is representative for all the simulation box. As the flux scaling is the same everywhere, the radiation field outside clumps is also attenuated with the same optical depth used for dense gas.

Local shielding approximation This kind of approximation considers that the main contribution comes from local gas. Approximations can either use just the local contribution to attenuation, or estimate the column density by using the local number density of the fluid element, multiplied by the size of the cloud.

This kind of approximation can underestimate the amount of shielding for low density bubbles embedded in cold and dense gas, while the shielding for high density clumps can be overestimated.

The main advantage of these approximations is that they do not solve the radiative transfer, they use instead estimates of the optical depth and column densities, but as we know, molecular clouds present a wide range of densities, with fractal-like structures, so these approximations are not able to reflect the influence of the environment, thus structure-dependent approximations are needed.

N-rays approximation One of the ideas that comes intuitively is to calculate the exact ray tracing for a great number of directions, but this is numerically demanding. The choice is in general to use a limited number of directions, although this leads to an extremely coarse angular resolution. Inoue and Inutsuka (2012) have used a "two-ray" approximation by assuming that the gas was well shielded against ultraviolet radiation, except in the colliding direction, while other authors have used a "6-rays" approximation (Glover and Mac Low 2007;[START_REF] Glover | Modelling CO formation in the turbulent interstellar medium[END_REF], which is another commonly used choice. In these two cases the integration of column densities is straightforward, the contribution of each cell is equal to N i = ndx i , where n i is the number density of the i -th cell and dx its size. Nonetheless, when more directions are included geometrical corrections are needed, and the method turns out to be more expensive.

The main inconvenience is that dense clumps produce a non-physical shadowing effect on the gas, and the extinction can be overestimated for cells that are located along the rays. From the numerical point of view, this method is not suitable for parallel computing, because it would require the information stored in other CPUs.

Gather approach This method discretizes the space in the same way as in the N-rays approach, but instead of calculating the contibution along a ray, it defines solid angles where it gathers all the matter and calculate the mean contribution. This method requires an additional radial discretization as shown if Fig. 4.2. This method gives an accurate estimate of the attenuation factor, taking into account the contribution of all elements in the simulation without shadowing effects. The accuracy of this method depends on the discretization, increasing with the number of directions used. This method is prohibitely expensive in terms of computational effort, and not suitable for parallel computing, however it can be used for validation.

TreeCol method This method is also a gather approach. It was developed by Clark et al. (2012), and implemented for the SPH code Gadget2. TreeCol uses the tree data structure of the code to construct a map of column densities, so the numerical cost is of order N log N , for a simulation with N elements. It uses a criterium on the opening angle. When the angular size of a node is larger than the tolerance Θ tol , the node is "opened" at a higher resolution until the criterium is satisfied.

This method offers a fast and accurate alternative to radiative transfer. The importance of shielding and its influence on the chemistry, the structure and the evolution of molecular clouds have motivated our efforts toward a general scheme for AMR codes able to capture the influence of the structure of the cloud, in order to obtain accurate estimates of the fraction of radiation from the ISRF that actually reaches each parcel of fluid in the simulation, and in a more general manner give accurate estimates of column densities, not only for dust, but for any other chemical species.

The developpement of this work is presented in the next section, as the validation tests and the application to the UV screening problem.

4.2 PAPER: Tree-based method for estimating column densities [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF] In this work we presented a method to estimate column densities in AMR codes, with a numerical cost of order N log N , which is fast enough to be used in numerical simulations. We introduced an innovative optimization technique based on the self-similarity of the tree structure, and we applied the tree-based method to the UV extinction by dust and its consequences on the thermal balance.

Introduction

Radiation plays an important role in several astrophysical processes at different scales, and it is dynamically coupled to the behavior of the gas. In particular, the far ultraviolet (FUV) radiation influences the physical and chemical properties of molecular clouds, while the cosmic rays (CRs)

are the dominant ionizating source in both diffuse and dense media, controlling the chemistry. The UV radiation is the main heating source for the gas by regulating the molecule formation rates on the grain surface, but at the same time dust grains shield the inner regions of clouds favoring low temperatures.

The penetration of FUV radiation into molecular clouds has been studied theoretically and numerically (Flannery et al. 1980;Sandell & Mattila 1975;[START_REF] Whitworth | [END_REF]Goicoechea & Le Bourlot 2007) by modeling the properties of dust grains or by simplifying the geometry of the cloud, but crucial for describing the propagation of radiation is to estimate the column densities.

Another example is the cosmic-ray ionization rate, which also depends on the value of the column densities. Indeed, several works show that the CR ionization rate decreases with the increasing value of the column density (see Takayanagi 1973;Padovani et al. 2009;Indriolo & McCall 2012).

The correct treatment of the propagation of radiation and the estimate of ionization rates involve calculating column densities, but this calculation can be numerically challenging in high-resolution multidimensional simulations, hence the interest of developing fast and accurate methods.

Some common approximations include methods based on ray-tracing schemes, such as Razoumov & Cardall (2005), that define ray domains according to the photon travel direction, or stochastic integration methods, such as Cantalupo & Porciani (2011), that use a Monte Carlo combined with an adaptive mesh refinement (AMR) strategy for the ray casting. There is a much simpler approach by Inoue & Inutsuka (2012), who assume that the gas is shielded well from UV photons except in the colliding direction and have used a 'two-ray' approximation for dealing with the shielding for the UV radiation. But most of these approaches have several disadvantages, and in general they are computationally expensive or not accurate enough when the geometry of the problem is more complicated, as in a turbulent medium. In spite of being intuitive, ray-tracing methods are numerically demanding. For a simulation with N resolution elements (cells or particles), the number of operations is at least on the order of N 2 and it requires the exchange of large amounts of data between different CPUs when parallelized. In problems where the gas properties are dynamically affected by gravity and radiation, it is desirable to develop adapted numerical strategies that permit calculation of the radiative transfer on-the-fly. Recently some efforts have been made in this direction. Using the smoothed particle hydrodynamics code GADGET2, Clark et al. (2012) introduce a tree-based scheme, called TreeCol. This method uses the information already stored in the tree structure of the code to construct a full 4π sr map of column densities for each element with a gather approach. Since the column densities are calculated while the tree is walked, its computational cost is also on the order of NlogN.

In this work we present a simple scheme for estimate column densities that takes advantage of the tree data structure, implemented for AMR codes. Our tree-based method provides a fast and relatively accurate estimate of column densities that can be used in numerical simulations. We have implemented this method in the AMR code RAMSES (Teyssier 2002).

In the following section we present the radiative transfer problem. In section 3 we present an overview of the tree-data structure and introduce our tree-based grid method, as well as our implementation. We also describe a strategy for optimization consisting in precalculation of geometrical contributions, where we introduce further approximations. In section 4 we present two tests for validation. The first test consists in a uniform spherical cloud and the second one corresponds to a turbulent cloud. At the end of this section, we also include a test for validating the other approximations made during the precalculation of the geometrical contributions. In section 5 we present one application of this method to calculate the UV absorption by dust and its consequences on the dynamics of molecular clouds. Section 6 concludes the paper.

The radiative transfer problem

The full radiative transfer equation that describes the interplay of radiation and matter reads as

dI ν ds = -α ν I ν + j ν . (1)
where I ν is the specific intensity, α ν the absorption coefficient, and j ν the emission coefficient at frequency ν. Because we are interested in the early evolution of molecular clouds, when stars have not been formed yet, we can consider that there are no local sources. Then the solution of Eq. ( 1)

for a case where there is only absorption will be

I ν (s) = I ν (s 0 ) exp - s s 0 α ν (s )ds = I ν (s 0 ) e -τ ν (2)
This simplified version of the solution lets us solve the problem by just calculating the optical depth τ ν along several directions. But even with these simplifications for a simulation with N resolution elements, the cost of solving the radiative transfer problem is on the order of N 5/3 . This can be easily done in postprocessing, but if we are interested in the dynamical interaction between radiation and matter, it is necessary to calculate the optical depths at each time step. This is extremely expensive in terms of CPU time for relatively large simulations and not practical in parallel computing, so that approached methods must be used. In this paper we compare our results to those obtained with a ray-tracing method, which uses 50 rays, performed in a postprocessing step.

Tree-based method

Astrophysical fluids can be treated numerically by using either (i) a Lagrangian approach, with either N-body or SPH codes (Benz 1988;Springel et al. 2001;Hubber et al. 2011), where nodes follow the material particles, or (ii) an Eulerian approach with patch-based codes (Fryxell et al. 2000;Almgren et al. 2010;Mignone et al. 2012;[START_REF] Enzo Collaboration | An Adaptive Mesh Refinement Code for Astrophysics[END_REF], where each element is fixed and describes how material flows through the grids.

Several of these astrophysical codes are based on a tree-data structure that consists in a hierarchical structure where the simulation domain is recursively split into smaller units or nodes. These nodes can be cut recursively to four (2D) or eight (3D) smaller "daughter" nodes from the largest node, the root, which contains the whole simulation, to the leaves that do not contain any substucture. Each node knows its parent node and can access all the other nodes by walking the tree. The numerical cost for a simulation with N resolution elements is proportional to NlogN (Barnes & Hut 1986, 1989).

Since most astrophysical problems span wide ranges of spatial scales, to correcly describe them, it is necessary to use resolutions comparable to the smallest scales. Adaptive mesh refinement (AMR) techniques allow the resolution to be adapted in different regions. This method was introduced for the first time in Berger & Oliger (1984) as an adaptative finite difference method for solving partial differential equations using nested grids in a patch-based AMR. The fully threaded tree (FTT) was proposed by Khokhlov (1998), where the tree is threaded at all levels and the refinement is done on a cell-by-cell basis.

The RAMSES code (Teyssier 2002) is a grid-based solver with AMR. The refinement levels are labeled . The coarse level ( = 0) corresponds to the base of the tree data structure and contains the whole simulation box. The refinement is done recursively on a cell-by-cell basis by splitting the cell into 2 3 daughter cells that constitute an oct, the basic elements in the data structure. The cells in an oct are indexed by the oct index, and this index is distributed around the center of the oct given a specific distribution.

This code has been parallelized using the MPI library and it uses locally essential trees [START_REF] Warren | Proceedings of the 1993 ACM/IEEE conference on Supercomputing, Supercomputing '93[END_REF], thus all the information is local. Each CPU knows the full tree up to a given resolution, therefore each cell can recursively access the information from the coarser levels.

General Idea

The tree-based method for estimating column densities is based on the fact that any distant cell substends a small angle and therefore its contribution to the extinction along the line of sight will be diluted. Then it is suitable to approximate the distant structured cells with cells at lower resolution.

For each target cell, column densities can be estimated by summing up all the contributions of cells along each line of sight and decreasing resolution with distance. Because all the information can be accessed by walking the tree and knowing the density, the distance, and the size of a cell at a given resolution level, we can calculate the contribution to the column density as the product of the density and the distance covered through the crossed cell along the line of sight.

Implementing and calculating the extinction

Since molecular clouds are embedded in the interstellar radiation field (ISRF), we are particularly interested in the influence of the far-interstellar radiation field triggered by the radiation of UV photons from OB stars. Habing (1968) suggested that this radiation density could be constant throughout space. Here we consider an incident UV field, the Draine field (Draine 1978), which is supposed to be monochromatic, that corresponds to G 0 ≈ 1.7 in Habing units, attenuated by dust (Wolfire et al. 1995).

We have implemented the calculation of column densities in order to estimate the attenuation for the UV radiation from the ISRF, described by the parameter G 0 (Habing 1968). From Eq. ( 2)

we can define the attenuation factor χ, calculated as the mean value of the extinction seen by a given cell and calculated as

χ = 1 4π 4π e -τ(θ,φ) dΩ , (3)
τ(θ, φ) = r 0 K(θ, φ)dr = σN(θ, φ) (4)
where τ is the optical depth along the line of sight defined by θ and φ in spherical coordinates, K is the extinction coefficient, and σ the effective attenuation cross section for the dust grains at λ = 1000 Å. Here we use σ d,1000 = 2 × 10 -21 cm 2 (Draine & Bertoldi 1996), and N is the total column density of hydrogen.

To calculate the column density maps and build the extinction maps, there are two different possible approaches. The first one uses the ray-tracing approach, where the column densities are defined by a single 'ray'. The second one uses a gather approach, where all the matter that falls in a given solid angle is gathered and added to the column density. The approach that we have adopted is quite similar to the ray-tracing. However, it permits us to add part of the matter that belong to neighboring cells by diminishing the resolution for distant cells, in such a way that the density corresponds to a mean. This has the advantage of considering the contribution to the Article number, page 5 of 28page.28 screening from neighboring cells that could be missed by using a simple ray-tracing algorithm. This approximation represents a lower computational cost than for an exact gather approach, where we would be required to perform more complex calculations.

We call the target cell the cell for which we currently want to estimate column densities, in order to keep the same nomenclature as Clark et al. ( 2012), and we call treated cells those contributing to the column density seen by the target cell. For the resolution we use the notation 0 the resolution level of the target cell the resolution level required for a cell contributing to the column density ( < 0 ).

For this approach we define directions based on a spherical projection centered on the target cell. We discretize the azimuthal angle in N regular intervals (δφ = 2π/N) and the polar angle in M irregular intervals constrained to the δ cos θ = 2/M constant in order to cover equal solid angles.

(Fig. 1 shows how these directions are defined.) Then the directions are labeled by two indices m and n, and the representative angles for these directions are given by

θ m = arccos 1 + (1 -2m) M (5) φ n = 2π N (n -1) (6)
where m ∈ [1, M] and n ∈ [1, N]. Since solid angles are equal and cells have uniform densities, we can rewrite Eq. ( 3) as

χ = 1 M × N m n e -σN(m,n) (7) N(m, n) = i n i ∆x i (m, n) (8)
where the index i stands for the cell number, n i and ∆x i (m, n) correspond to the number density of the cell i and to the distance crossed through the cell i in direction (m, n) with respect to the target cell, respectively. For integrating column densities N(m, n) along the different directions, we define two main regions. The first region contains the cell itself and its siblings (the cells that belong to the same oct, sharing a common parent node), and the contribution to the column density is calculated at resolution level 0 . The outer region contains the rest of the cells in the simulation box, and the contribution to the column density is estimated by decreasing the resolution with the distance.

For the outer region, we construct cubic concentric domains, where cells are treated at the same resolution. For each level of resolution, starting from = 0 -1, we define a cubic shell as -We set the new center as the center of the grid that contains the target cell at level + 1 (or equivalently the center of the unrefined cell at level ); -We set the inner limit as the external limit of the previous shell;

-We define the outer limits for each cartesian direction by taking two or three neighboring cells into account at level in order to fully cover their parent grid.

The procedure is repeated recursively up to the border of the box. These cubic concentric shells define at which resolution the cells are taken into account for the calculation of the column density.

Figure 2 illustrates the construction of these concentric cubic domains.

To integrate the column densities in each direction we define three contributions: the internal contribution, corresponding to the target cell itself, the local contribution, given by the sibling cells, and the external contribution, given by the cells in the outer region.

For the internal and local contributions, directions are defined with respect to the center of the target cell. On the other hand, for the external contribution, directions are defined with respect to the center of the grid that contains the target cell at level 0 . This external contribution will be the same for all the sibling cells that belong to the same oct. Using Eq. ( 8) and rewriting the sum as a sum over concentric shells C , we can use the fact that all the cells that belong to the same shell have the same size dx to write the external contribution as

Article number, page 7 of 28page.28

N ext (m, n) = i∈C n i ∆x i (m, n, ) (9)
∆x i (m, n, ) = dx × K i (m, n, ) (10)
where C stands for the shell at level . The size of the cell at level is given by dx = 0.5 L, with L the size of the simulation box. The multiplicative factor K i contains the geometrical corrections to the distance crossed through the cell i with respect to the target cell in the direction (m, n) at level . Finally, the column densities along each line of sight can be estimated by adding all the contributions.

Optimization and precalculation module

The geometrical factor K i (m, n, ) in Eq. ( 10) is quite expensive to calculate; however, because the oct structure is self-similar this factor can in principle be calculated and stored at the beginning of the simulation. Then, using the fact that octs are self-similar at different resolution levels and knowing by construction that the outer limit for any cubic shell is at most five times the local dx in each cartesian direction (±x, ±y, ±z), we can define a standard cube of size 11 3 total cells (Fig. 3) of a one unit size. The cell in the center contains the target cell, so we can represent all the possible cases by defining its internal configuration, depending on the difference between the target cell level and the shell level, given by ∆ = 0 -. Therefore the total number of configurations is given by max -2 ∆ =0 (2 3 ) ∆ . For each one of these cases, we need to calculate the correction for any cell in the cubic shell in any direction. This implies M × N × 11 3 corrections for each case. This leads to a huge matrix that cannot be stored. We can, however, reduce the size of the matrix by restricting the number of cases considered. This is done by requiring that the target center is shifted for high values of ∆ .

As a first approximation for a shell at level , the center of the target cell can be approximated by the center of the grid that contains the target cell at level + 1. This is equivalent to replacing it by the target cell center at level . This is the exact configuration when 0 -= 1, but for other levels, the position of reference can be considerably drifted away, and then more precise corrections are needed. We introduce two correction levels in order to take the configuration of the target cell into account with respect to its coarser octs. The first correction level is used when 0 -= 2 and it considers the index of the target cell in its oct at level -1, so that there are eight possible positions.

The second correction level is used for the case 0 -≥ 3, where we consider two oct indices. The first one is the oct index of the cell that contains the target cell at level + 2, and the second one corresponds to the oct index at level + 3, generating 64 possible configurations. This leads to a total of 73 possible configurations. The effect of this approximation on the estimate of individual column densities is depicted in Fig. 4. Then for each configuration we calculate the distance crossed through every cell in the general cube for each one of the directions (m, n) in order to obtain the corrective factor K defined in Eq. ( 11). Finally the pre-calculation module generates a matrix that contains the corrective factors K for all the configurations considered and a boolean matrix that permits us to know if the geometrical correction is non-zero. These matrices are calculated just once at the beginning of the simulation.

The external contribution will be calculated as

N ext (m, n) = dx i ∈ C n i × K(r i -r target , m, n, ∆ ), (11)
then the corrective factors K can be found by knowing the indices of the required line of sight m and n, the configuration within the octs given by the ∆ = 0 -, and the relative position (r i -r target ) of the cells in the oct with respect to the target cell in dx units. As described later in Sect. 4.4, this speeds up the code by a fair amount without reducing too much the accuracy.

Verification of the method

In this section we analyze the reliability of our tree-based method by comparing our estimates to a reference. At the end of this section, we analyze the influence of the precalculation of geometrical terms and the approximation introduced in section 3.3 on the performance and accuracy.

To validate of our method we considered two test cases. The first one consisted of a uniform spherical cloud, and the second one corresponded to a turbulent cloud. For both cases we present the total column density maps integrated along each one of the main axes x, y and z with respect to the midplanes and the extinction maps as seen by the cells in the midplanes calculated with our treebased method for the following cases: 6, 12, 40, and 84 directions. To show the error dependence on the position angle, we present column density projections onto 4π sr maps as seen by a cell for both test cases.

The reference maps were calculated in postprocessing with two different methods. For the column density and extinction maps in a plane, we used a ray-tracing approach that uses 50 rays and takes all the cells of the simulation into account at the highest resolution. It calculates the exact contribution to the column density for each direction using all the cells that are intersected by the ray. For the 4π sr maps and for the mean column density maps, the reference was calculated with a gather approach. This method divides the computational domain in angular and radial bins. It finds the closest bin for each cell and calculates the fraction of the mass that falls in this bin and in angular neighboring bins. Because the geometrical corrections for cells are much more complicated than for spherical particles, the fraction was calculated assuming that the cell is seen perpendicular to one of its faces. The column density is calculated as the sum of the mean density in the bin multiplied by the radial thickness of the bin.

These two methods have several inconveniences. They require the information already stored in each CPU, so they are not suitable for parallel computing, and they are extremely expensive in terms of CPU time. In particular, the ray-tracing method has a numerical cost of about N 4/3 times the number of rays, while the gather method has a numerical cost of about N 2 , since for each cell all the cells are used for calculating the column densities. The reference maps for the extinction, column density, and mean column density are calculated at the same time. The cost of producing the maps for three slices using the ray-tracing method for a maximum resolution level = 9 is 90 hr of CPU-time, while for the maps done using the gather method, the cost for only one slice is more than 6800 hr of CPU-time. The enormous cost of the gather approach meant that for the mean column density, we performed only one map at the maximum resolution ( = 9) for the cut at z = 25 pc. For the other cuts, we used only = 8.

For all the comparisons, the fractional error is defined by Eq. ( 12). Because we were interested in how the UV field is shielded, the extinction maps were normalized to one, and the difference maps were calculated by Eq. ( 13):

∆N i = |N i -Ni | Ni ( 12
)
∆χ i = |χ i -χi | ( 13
)
where N and χ stand for the reference maps for the column density and extinction, respectively.

Spherical uniform cloud

To test our method we first considered the simple case of a spherical cloud of uniform density. The radius of the spherical cloud is 12.5 pc, and its number density 1000 cm -3 . The cloud is located at the center of a 50 pc cube with a number density of 10 cm -3 . We considered two AMR levels with a maximum resolution equivalent to 512 3 cells. The refinement criterion is the density set in such a way that the sphere is refined.

The column density maps integrated along the main axes (in Fig. 5) show that our approximation can reproduce the main features, and it is able to capture the discontinuity in density at the edge of the cloud. The error maps show that in general the error is lower than 10%, but at the edges, the error increases up to more than 100% in a narrow region. Since most of the volume around the sphere is filled with constant density gas, we calculated the mean error on the sphere alone. This mean error is about 15.2% when averaging for the three maps.

Since the error in the column density maps is higher at the edge of the sphere, we calculated the column densities as seen by a cell sitting at (12.25, 0.5, 0.5) pc. In Fig. 6 we show the column density projection in every direction as seen by a cell at the border of the sphere for 12, 40, 84, 144, and 220 directions, using the gather approach and the tree-based method. In the same figure, we present the relative error. For this cell, the center of the sphere is in direction (θ, φ) = (π/2, 0) so that one half of the sphere is seen at the beginning, and as the angle φ increases, the column density seen by the cell decreases. When φ = π, the direction points outward, where the density is weaker. As the angle φ keeps increasing, the direction approaches the sphere again. This figure

shows that the highest error happens at the edge of the sphere. A similar effect can be seen in Fig. 7 for TreeCol (Clark et al. 2012). In our case this error probably occurs because distant cells are not as well described by the tree-based method.

For the extinction, we present a comparison in Fig. 7 of the extinction maps for the cells in the midplanes, as defined above. The reference maps calculated with the ray-tracing method are shown in the top row. Then from top to bottom, we present the extinction maps calculated with the tree-based method for 6, 12, 40, and 84 directions. In Fig. 8 we present the respective difference ). The abscissa corresponds to the angle φ (from 0 to 2π), and the cos θ (from -1 to 1) is given in the ordinate. The reference maps (on the left) have been calculated using the gather approach. The panels at the center present the maps calculated using the tree-based method, while on the right we present the error maps calculated according to Eq. ( 12). From top to bottom: using 12, 40, 84, 144, and 220 directions. 13). For the case with six directions, the mean difference is about 0.036, while for the cases at 12, 40, and 84 directions the mean difference is less than 0.02, but the increasing number of directions does not improve the accuracy considerably (Table 1). This is probably because the resolution is getting coarser as we cross cells located farther away. While the number of directions is increasing, the angular resolution remains constant owing to limitations inherent in the method. In this manner, more intervals oversample the region without increasing the angular resolution.

Article number, page 12 of 28page.28 Fig. 7. Extinction maps as seen by the cells in the midplanes for x, y, and z. From top to bottom: the reference calculated with a ray-tracing method, then using the tree-based method for 6, 12, 40, and 84 directions.

Turbulent cloud test

The previous test presents a very simple geometry, but in real atrophysical problems the geometry of the medium can be quite complicated. With the purpose of studying the accuracy of the treebased method in a highly structured medium, we analyze the case of a turbulent molecular cloud.

This cloud has been formed through a simulation of colliding flows (Audit & Hennebelle 2005;Vázquez-Semadeni et al. 2007) with a turbulent velocity profile. This simulation produces a cloud that presents filaments and clumpy structures with densities ranging over more than six orders of magnitude. The size of the simulation box is 50 pc, and we used two AMR levels ( min = 7 and max = 9), with an equivalent maximal resolution of 512 3 or, equivalently, a spatial resolution of 0.1 pc. The criterion for refinement is again the density (n = 50 cm -3 for the first refinement and n = 100 cm -3 for the second one).

Article number, page 13 of 28page.28 Figure 9 presents the column densities integrated along the main cartesian directions as seen by the cells in the midplanes. At the top we show the reference maps, calculated with the ray-tracing method and, in the middle, calculated with our tree-based method. At the bottom we present the relative error maps calculated according to Eq. ( 12). The mean value of the relative error is about 50%, but for many cells it can reach up to more than 100%. The reference corresponds to an exact integration along the x, y, and z axes, so even far structures are taken into account at the highest resolution. With the tree-based method, these far structures are distributed over a larger region. This can be seen by comparing the column density maps obtained with the tree-based method to the reference. The missing structures can be interpreted as structures far from the midplane, while the clumps seen using the tree-based method are structures close to the midplane.

For the turbulent cloud we present a 4π sr map similar to Fig. 6. Figure 10 shows the column density projection in every direction as seen by a cell sitting at the center of the turbulent cloud for 12, 40, 84, 144, and 220 directions using the gather approach , the tree-based method, and the associated relative error. This map, unlike the case for the spherical cloud, does not show such high errors. Moreover, the errors seem to be distributed evenly throughout the map.

Figure 11 presents the extinction maps. At the top we present the reference maps calculated with the exact ray-tracing method, and the rest of the maps correspond to the tree-based method.

In Fig. 12 we present the difference maps calculated as defined before. These figures show how Article number, page 14 of 28page.28 the accuracy of the maps can be improved by including directions that are not aligned with the main cartesian directions, because more cells are taken into account. The mean and maximum differences are summarized in Table 1. When we use six directions, the mean difference is about 0.079, while for the rest it is about 0.05. Beyond the mean difference we see that the 'six-ray' method presents large and systematic variations that are absent with the multiray approach. As we have shown before for the spherical cloud test, the increasing number of directions does not improve the accuracy of the map considerably, particularly between 40 and 84 directions.

In Fig. 13 we present the mean column density seen by cells in the midplanes. At the top we present the mean colum density maps done with the ray-tracing method. This method present strong shadowing effects for all the cells that are aligned with dense clumps, while for directions that are not aligned, the entire contribution is missed. For this reason we have included reference maps done with the gather method using 220 angular bins (second row). The prohibitive cost of producing these maps led us to calculate the maps for the cuts at x = 25 pc and y = 25 pc at a lower resolution for a visual comparison. The map for the cut at z = 25 pc was calculated at level = 9, and it was used to calculate the relative error. This figure shows how the estimate of the value of the mean column density improves as we increase the number of directions used in the tree-based method. In particular the mean error decreases as the number of directions increases, as shown in Table 2.

General remarks

That the accuracy does not improve considerably is consistent with the results of Clark et al. (2012) for TreeCol, where they show that the efficiency depends on the relative sizes of the opening angle Article number, page 15 of 28page.28 Fig. 10. Column densities for the 4π sr as seen by a cell sitting at the center of the turbulent cloud. As in Fig. 6, the horizontal axis is the azimuthal angle φ (from 0 to 2π), and the vertical axis corresponds to cos θ. On the left we present the reference maps calculated using the gather approach. The panels in the center present the maps calculated using the tree-based method, and on the right we present the error maps calculated according to Eq. ( 12). From top to bottom: using 12, 40, 84, 144, and 220 directions. and the angular size of the node. This is true for a pure gather approach. It is important to notice that our implementation of the tree-based method is a hybrid method. The decreasing resolution for far cells mimics a gather approach, but calculation of the contributions to the column densities is done as in a ray-tracing method. With this in mind, we note that the increasing number of directions does not change the angular resolution for a given cell, but it does improve the quality of the description of the density field for calculating the extinction. This can be seen in Fig. 11, where the extinction maps are smoother and the maximum error is reduced as the number of directions is increased. Then the optimal number of directions for the tree-based method will correspond to the best compromise between resolution and numerical cost. Taking this into account and using as a criterium the variation of the error, a good compromise is found for 40 directions.

Article number, page 16 of 28page.28 In particular, the highest error happens toward the border of the computational domain for the spherical case, while the error in this region is much smaller for the turbulent cloud. This is probably due to the different distribution of sources. This effect is more evident for the spherical case than for the turbulent cloud owing to the anisotropy of the density field. For the turbulent cloud the distribution of sources is more isotropic, and for each cell it will be easier to cast a source, while for a compact central distribution, especially for cells not aligned with the central source, it will be more difficult to describe the density field. This same effect can be seen in Figs. 10 and6.

Validation of the precalculation approach

To quantify the influence on the accuracy of estimating the attenuation factor of the two-level approximation and the influence of precalculating geometrical terms on the speed of the code, we Article number, page 17 of 28page.28 and the two-level approximation. We present a comparison of the column density maps integrated along three different directions. For the extinction we present a comparison of the extinction maps and the difference on the estimated value of the extinction. We also present the calculation CPU times (t sim ) relative to the standard CPU time (t 0 ). The standard CPU time is the calculation time of a simulation that does not include the extinction. The two-level approximation accelerates the code, but it shifts the point of reference that defines the directions (Fig. 4) introducing errors on the calculation.

To estimate the gain in performance and the induced errors, we performed eight different calculations. The first one corresponds to the reference, calculated without including the screening for Article number, page 18 of 28page.28 Fig. 13. Mean column density as seen by cells in the midplanes. From top to bottom: using the ray-tracing method for 50 rays, using the gather approach for 220 directions, and using the tree-based method for 6, 12, 40, and 84 directions.

the UV. The second one uses a simple implementation of our method using just six directions, and it does not need any kind of geometrical correction. The six other simulations use two implementations of our method (calculation in situ and precalculation module) for 3 × 4 = 12, 5 × 8 = 40, and 7 × 12 = 84 directions. We use a turbulent cloud as initial condition, which is created by converging flows without any screening, and we compare the calculation time. For the comparison of the column density maps and the extinction maps we use the initial output, in order to use identical clouds. Table 3 presents the simulation times relative to t 0 , the calculation time without screening, and the mean value of the difference in the extinction maps. Calculating column densities for the Article number, page 19 of 28page.28 Fig. 14. Influence of the precalculation module. Column density maps for three directions not aligned with the cartesian directions for the case of calculation in situ (top), using the precalculation module (middle), and the relative difference map (bottom). The azimuthal angle is fixed at φ = 0, and from left to right the polar angle is θ = 0.841, 0.644, 0.541 rad. The mean fractional difference, defined as in Eq. ( 12), are ∆N = 6.3, 6.6, 6.8%, respectively. For estimating of the induced differences related to the two-level correction we calculated column density maps and extinction maps for both implementations. For the column density maps, we

Article number, page 20 of 28page.28 selected three directions that are not aligned with any of the main cartesian directions ±x, ±y, ±z, because these directions are not sensitive to this approximation. Figure 14 shows the column density maps for three directions where we fixed the azimuthal angle (φ = 0) and we varied the polar angle (θ = 0.841, 0.644, 0.541 rad). The relative difference map was calculated according to Eq. ( 12), and the mean value for all of them is less than 7%.

For the extinction we present a comparison between both implementations. We calculated the difference of the extinction maps as seen by the cells in the midplanes according to Eq. ( 13), for the case where we introduce the precalculation module with respect to the case where all the calculations are done in situ. This comparison was done using 12, 40, and 84 directions. Table 3 presents the mean values of the difference, which are better than 0.004 for all the cases, and it does not depend on the number of directions used. As an example, we present the extinction maps and difference maps for the cells in the midplanes in Fig. 15 as a comparison of both implementations for the case where we use 84 directions (7 intervals for the polar angle and 12 for the azimuthal angle).

Overall we conclude that the method with precalculation is sufficiently accurate and considerably faster than the in situ method.

Application

As an application we study the formation and the evolution of a molecular cloud formed from colliding streams of warm atomic gas (Audit & Hennebelle 2005;Vázquez-Semadeni et al. 2007; Article number, page 21 of 28page.28

V. Valdivia & P. Hennebelle: A fast tree-based method for estimating column densities in AMR codes Heitsch et al. 2005Heitsch et al. , 2006Heitsch et al. , 2008)). The 3D simulations were performed using the AMR code RAM-SES.

The set-up and the initial conditions are similar to those of Hennebelle et al. (2008) (see also Vázquez-Semadeni et al. 2007;Heitsch et al. 2005Heitsch et al. , 2006Heitsch et al. , 2008;;Inoue & Inutsuka 2012). We consider a cubic box of length L = 50 pc. We allow two AMR levels, with min = 8 and max = 10, reaching an effective numerical resolution of 1024 3 cells and a spatial resolution of about 0.05 pc.

The boundary conditions are imposed to mimic the large scale converging flows. The gas is injected from the left and right faces of the simulation box with a weakly turbulent velocity and a density of 2 cm -3 at a temperature of 8000 K. For the remaining faces, we use periodic boundary conditions.

The velocity field of the incoming gas V in depends on y, with an average velocity V 0 = 15 km s -1

and modulated by a function of amplitude = 0.5, as defined in Audit & Hennebelle (2005).

Initially the simulation box is filled with warm atomic gas with the same density and temperature as the inflowing gas and is uniformly magnetized with a magnetic field of strength 2.5 µG parallel to the x axis and therefore aligned with the incoming velocity field. The gas in the simulation box is heated by the background UV field, which corresponds to the Draine field G 0 = 1.7 in Habing units. Because this radiation field is assumed to be isotropic and constant, then we can define the effective UV field as

G 0 = χG 0 ( 14
)
where χ is the attenuation factor for the UV flux defined in Eq. ( 3) and calculated for each leaf cell in the simulation.

We use the same cooling and heating functions than Audit & Hennebelle (2005), but we have modified the heating by adding the screening for the UV and we included the heating by cosmic rays. The heating rates as implemented in the code are

Γ UV = 1.0 × 10 -24 ε G 0 n ergs cm -3 s -1 (15) Γ CR = 10 -27 n ergs cm -3 s -1 ( 16
)
where Eq. ( 15) is the heating rate due to the photoelectric effect on small grains and PAHs due to the FUV radiation (Bakes & Tielens 1994), n is the hydrogen density in cm -3 , ε is the heating efficiency, calculated as in Wolfire et al. (1995) :

ε = 4.9×10 -2 1.0+ G 0 T 1/2 /n e /1925 0.73 + 3.7×10 -2 (T/10 4 ) 0.7 1.0+ G 0 T 1/2 /n e /5000
, with n e , the electron density, given by the approximation proposed by Wolfire et al. (2003). Equation ( 16) is the heating rate due to cosmic rays according to the intermediate value given by Goldsmith (2001).

Our aim is to understand how the screening caused by the surrounding matter can affect the gas distribution and the distribution of structures formed by diminishing the amount of radiation Article number, page 22 of 28page.28 that arrives at the cells. Then we present two simulations in order to compare the influence of the screening for the UV. The first simulation does not include this effect, which is equivalent to considering the gas as optically thin and having an attenuation factor χ = 1. This means that the radiation reaches the cells unchanged. On the other hand, the second simulation takes the absorption due to the surrounding material into account. In this case the estimate is done using the treebased method, as described in section 3. Just for the sake of simplicity and time, we used 12 directions (M = 3 intervals for the polar angle and N = 4 for the azimuthal angle) for calculating the attenuation factor for the UV field (but see the Appendix A for a comparison between simulations at lower resolution using different number of directions to calculate the extinction for the UV field).

General structure of the cloud

To understand how the screening for the UV affects the structure of molecular clouds, we present a comparison between both simulations in Fig. 16 by presenting the local density in a slice cut through the middle of the simulation box at t = 10 Myr. This figure shows that the large scale structure of the molecular cloud seems to be barely affected by the UV screening, but the detail of the fragmentation of dense structures is substantially influenced by the extinction. In the case where the extinction is not taken into account, the structures formed tend to be bigger, while when the screening is included, the same region seems to be torn up, presenting smaller structures. This seems to indicate that the extinction mainly affects the dense parts of the gas. To go beyond this qualitative impression, we now turn to more quantitative studies.

Probability distribution function

The density probability distribution function (PDF) corresponds to the distribution of mass as a function of density and is one of the simplest statistical tools for understanding how the gas is distributed. In Fig. 17 we present a comparison of the PDFs and the mass-averaged temperature for each density bin. For both cases we can see that starting with only warm gas, which is caused by the cooling processes present in the simulation, the gas is able to transit from the WNM phase to the CNM phase (Audit & Hennebelle 2005). As the gas enters the box, the CNM phase develops. The PDF of the gas density and the temperature distribution show the bistable nature of the medium (Field et al. 1969;Wolfire et al. 1995Wolfire et al. , 2003)). The panel on the right shows how the screening affects mainly the dense medium, while the warm gas is almost not altered. Initially the distribution is not significantly affected. However the difference can be seen as dense gas forms. The bottom left hand panel shows that a part of the gas in the CNM phase presents higher density and lower temperature. Consequently, part of the gas has a lower density owing to the rarefaction caused by the development of more compact structures, and this explains why there is more gas with n ∼ 300 cm -3 in the case without screening.

Mass spectra

To investigate the influence of the UV screening on the statistical properties of the structures formed within the clouds, we analyzed the mass spectra of the clumps for both cases. The mass spectrum presents the number of structures per logarithmic mass interval for a given density threshold. For the clump extraction, we selected all the cells with a density higher than a given density threshold n th , and using a friend-of-friend algorithm, we identified the spatially connected regions that constitute a clump, rejecting isolated cells. In Fig. 18 we present the evolution of the clump distribution for t = 5, 7.5, and 10 Myr for different density thresholds. We show that the mass spectrum for low density thresholds does not vary significantly when we include the UV screening. On the other hand, the differences in the mass spectrum are more pronounced for more compact structures (n th ≥ 2500 cm -1 ). The number of compact clumps found is systematically larger in the simulation that includes the UV screening.

Conclusions

We have introduced a tree-based method for a fast estimate of column densities in astrophysical simulations. The general idea of this method can be implemented in any code with a tree-based data structure, and the implementation strategy is not unique. We presented a simple implementation on the AMR code RAMSES. In particular, we used a precalculation module that speeds up the calculation considerably without changing the accuracy of the estimate significantly (better than 7% for the column density maps and about 0.3% for the extinction maps).

Since the cells that contribute to the column density are taken into account at lower resolution as the distance increases, and if the tree is fully threaded, the tree-based method only needs the information stored locally in the essential tree and does not need any communication between dif-Article number, page 24 of 28page.28 ferent CPUs, making the method suitable for parallel computing. Finally, the errors on the column density maps are generally about 50%, while the extinction maps calculated from the estimated column densities have errors lower than 10%.

We found that the UV screening does not have a strong influence on the general structure of molecular clouds, but it has a significant impact on the details of fragments. Most notably, the extinction for the UV mainly affects the dense gas. The temperature for the WNM phase remains almost unchanged, while for the CNM it is lower by a factor up to 50% when the shielding effect of dust is included. Consequently, the local value of the Jeans mass is reduced, favoring the gravitational collapse of smaller structures. This is evident in the comparison of the mass spectra, where for low-density thresholds there are no noticeable differences; however, the number of compact structures for higher density thresholds is greater for the case that includes the screening.

Because an important part of the chemistry in the ISM depends on the UV radiation and on temperature, this method can be applied to the interstellar chemistry to more realistically estimate these parameters. Estimates of column densities can also be used to give a better value for the CR ionization rate. and that the differences remain within the limits of the statistics.

Article number, page 27 of 28page.28 H 2 formation and destruction

It has been shown that molecular gas, on galactic scales, is well correlated with star formation. Different works (Wong and Blitz 2002;Bigiel et al. 2008[START_REF] Bigiel | A Constant Molecular Gas Depletion Time in Nearby Disk Galaxies[END_REF][START_REF] Leroy | The Star Formation Efficiency in Nearby Galaxies: Measuring Where Gas Forms Stars Effectively[END_REF]) have found that the surface density of star formation, Σ SF R , is correlated with the molecular hydrogen surface density, Σ H2 , but, whether molecular gas is a necessary condition to star formation, or, on the other hand, it is a consequence, it is not clear. In seems that rather a causeconsequence relation, both quantities would be correlated to a third variable. [START_REF] Glover | Is molecular gas necessary for star formation[END_REF], performing numerical simulations, have concluded that H 2 and CO formation have a minor influence on the star formation rate, and the crucial parameter would be the ability of clouds to shield themselves against the interstellar radiation field, requiring large amounts of dense gas, a requisite for both star formation and molecule formation. What seems to be clear is that the H 2 formation and evolution of the molecular fraction is determined by the environment, as shown by Copernicus and FUSE observations of the HI-H 2 transition (Savage et al. 1977;Gillmon et al. 2006), that seems to point out the existence of a shielding threshold.

As we have seen in Chapter 2 H 2 is the first molecule to be formed, with a formation time of the order of t ∼ 10 9 /n yr for gas at density n in cm -3 . Molecular clouds, in their evolution, experiences fluctuations in density due to turbulent motion. Micic et al. (2012) have shown that overdensities created by compressional motions can enhance the formation of H 2 .

In order to understand all these effects dynamically, local density enhancement and structure able to shield the gas, it is needed to perform numerical simulations that take them into account. Previous works includes simplified descriptions of the self-shielding, which is the ability of molecular gas to absorb photodissociating photons (see Section 2.3.2), or do not include multi-phase gas, nor the full process of formation of a molecular cloud from WNM gas.

The self-shielding factor depends on the amount of H 2 gas that surrounds a parcel of fluid. In order to compute it properly it is needed to calculate H 2 column densities. We have adapted the Tree-based method to estimate H 2 column densities and self-shielding factors.

In the following section we present a study on the H 2 distribution and evolution during the formation of a multi-phase molecular cloud.

PAPER: H 2 distribution during formation of multi-phase molecular clouds

Realistic molecular clouds present clumpy and filamentary structures, where warm and cold phases are interwoven in a very complex structure. To understand the H 2 molecule formation process under the dynamical influence of a highly inhomogeneous structure, we have performed high-resolution MHD simulations of realistic molecular clouds formed through colliding streams of warm atomic gas. For this study we have used the Ramses code, where we have included the formation and destruction processes for H 2 , as well as the thermal feedback. The effects of dust shielding for the ultraviolet radiation, as the self-shielding due to H 2 molecules, have been included by using our tree-based method, detailed in [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF].

We find that H 2 molecular gas is formed faster than the usual estimates based on the mean density of the cloud. We interpret this as the result of the combined effect of the local density enhancement, that drives a faster H 2 formation, and the shielding provided by the global structure, that ensures the survival of molecular gas. Additionally, the mixing between different phases induce the presence of H 2 molecules in the warm phase, that under the shielding effects can survive in warm gas, explaining the presence of warm H 2 in the conditions that prevail in molecular clouds.

Our simulations reproduce total column densities as those of the translucent lines-of-sight observed by Copernicus and FUSE. We compared our results with these observations and we found a general agreement. In particular the molecular fraction as a funtion of the total column density present the same two regimes observed previously, as well as the HI-H 2 transition toward column densities N with log(N [cm -2 ]) 20.3, differences arising possibly from the different UV field strength conditions. We have also calculated the H 2 populations in the first levels at thermal equilibrium, as well as the Doppler broadening parameter b for each population. Excitation diagrams and b parameters are in good agreement with the values observed by Copernicus and FUSE. Altogether, these results suggest that excited populations might be the consequence of local temperature and density within molecular clouds, highlighting the multi-phase nature of molecular clouds.

The molecular richness observed in the ISM is not fully explained. The existence of warm H 2 can open new paths for further molecule formation and eventually could explain certain aspects of the warm chemistry. Further efforts must be done in this direction to understand whether our model can explain the observed amounts of other molecules in excited rotational states.

Introduction

It is well known that stars form in dense and self-gravitating molecular clouds and that the star formation rate per unit area (Σ SFR ) is on average relatively well correlated with the H 2 surface density (Σ H 2 ) (Lada et al. 2010(Lada et al. , 2012;;Wong & Blitz 2002;Bigiel et al. 2008).

Although H 2 is the simplest and the most abundant molecule in the interstellar medium (ISM), it is very hard to observe directly. Due to its homonuclear nature, H 2 lacks of permament dipole and only weak quadrupolar transitions are allowed. H 2 can be observed in emission by infrared (IR) rovibrational transitions (Burton et al. 1992;Santangelo et al. 2014;Habart et al. 2011), or in absorption at far ultraviolet (UV) wavelengths from the Lyman and Werner bands. These bands were first observed by [START_REF] Spitzer | Physical processes in the interstellar medium Spitzer[END_REF] using the Copernicus satellite. The Far Ultraviolet Spectroscopic Explorer (FUSE) (Moos et al. 2000) offered new perspectives on the study of H 2 in the ISM due to its sensibility 10 5 times higher than Copernicus in the far-UV part of the spectrum, providing measurements of H 2 column densities (total column density of H 2 and for several rotational levels J) along translucent lines of sight (Wakker 2006;Sheffer et al. 2008).

While it seems to be clear that molecular gas is well correlated with star formation at different scales, how the gas becomes molecular is still an open issue. Since understanding the atomic-tomolecular hydrogen (HI-to-H 2 ) transition is of the highest importance for understanding the star forming process, numerous models have been developed (e.g. Krumholz et al. 2008;Sternberg et al. 2014) and seem able to reproduce many of the observed constraints. While extremely useful, these models leave aside the dynamical aspects of H 2 formation. In particular, the question of the relatively long time scale that is needed to form the H 2 molecules has been for many years difficult to reconcile with short-living and quickly forming molecular clouds. Indeed the typical time scale for H 2 formation is on the order of 10 9 yr/n (Hollenbach et al. 1971), which would lead for molecular clouds of mean densities of the order of 10-100 cm -3 (Blitz & Shu 1980) to ages larger than 10 7 Myr. However, Glover & Mac Low (2007a) and Glover & Mac Low (2007b), who have been simulating the formation of molecular hydrogen in supersonic clouds, show that H 2 can form relatively quickly in the dense clumps induced by the shocks, leading to a time scale significantly shorter (typically of about 3 Myr). Recently, Micic et al. (2012) explored the influence Article number, page 2 of 31page.31 of the turbulent forcing and show that it has a significant influence on the timescale of H 2 formation. This is because compressible forcing leads to denser clumps than solenoidal forcing for which the motions are less compressible. This clearly shows that dynamics is playing an important role in the process of HI-to-H 2 transition, at least at the scale of a molecular cloud.

In the Milky-Way most of the molecular hydrogen is in the form of low temperature gas, nevertheless different observations have brought to light the existence of large amounts of fairly warm H 2 in the ISM (Valentijn & van der Werf 1999;Verstraete et al. 1999). Gry et al. (2002) and Lacour et al. (2005) have shown that the H 2 excitation resulting from UV pumping and from H 2 formation cannot account for the observed population of H 2 excited levels (J > 2). Excited H 2 is likely explained by the presence of a warm and turbulent layer associated to the molecular cloud. But in such warm gas, characterized by very low densities (n H ≈ 1 -10 cm -3 ) and very high temperatures (T ≈ 10 3-4 K), H 2 formation on grain surfaces becomes negligible therefore leading to an apparent contradiction. Recently, Godard et al. (2009) have been investigating the possibility that warm H 2 could form during intermittent high energy dissipation events such as vortex and shocks. They show in particular that under plausible assumptions regarding the distributions of these events, the observed abundance of excited H 2 molecules could be reproduced. Another possibility to explain the abundance of these excited H 2 molecules is that they could be formed in dense gas and then transported in more diffuse and warmer medium. This possibility has been investigated by Lesaffre et al. ( 2007) using a 1D model and a prescription to take into account the turbulent diffusion between the phases. In particular, they find that the abundance of H 2 molecules at low density and high temperature increases with the turbulent diffusion efficiency.

Therefore while different in nature, previous works (Glover & Mac Low 2007a,b;Lesaffre et al. 2007;Godard et al. 2009) consistently found that the formation of H 2 is significantly influenced by the dynamics of the flow in which it forms. Although the exact mechanism that leads to the formation of molecular clouds is still under investigation, it seems unavoidable to consider converging streams of diffuse gas (e.g. Hennebelle & Falgarone 2012;Dobbs et al. 2014). What exactly triggers these flows is not fully elucidated yet but is likely due a combination of gravity, large scale turbulence and large scale shell expansion. Another important issue is the thermodynamical state of the gas. The diffuse gas out of which molecular clouds form, is most likely a mixture of phases made of warm and cold atomic hydrogen (HI) and even possibly somewhat diffuse molecular gas. Such multi-phase medium presents large temperature variations (typically from 8000 K to less than 50 K). This makes the dynamics of the flow significantly different from isothermal or nearly isothermal flow (see Audit & Hennebelle 2010, for a comparison between barotropic and 2-phase flows).

Moreover the 2-phase nature of the diffuse atomic medium has been found to persist in molecular clouds (Hennebelle et al. 2008;Heitsch et al. 2008a;Banerjee et al. 2009;Vázquez-Semadeni et al. 2010;Inoue & Inutsuka 2012). This in particular produces a medium in which supersonic clumps of cold gas, are embedded into a diffuse phase of subsonic warm gas. Consequently, the question arises: what are the consequences of the multi-phase nature of a molecular cloud onto the forma-Article number, page 3 of 31page.31 tion of H 2 ? Clearly, various processes may happen. First of all, the dense clumps are continuously forming and accreting out of the diffuse gas because molecular clouds keep accreting. Secondly, some of the H 2 formed at high density is likely to be spread in low density, high temperature gas because of the phase exchanges.

In this work, we present numerical simulations using a simple model for H 2 formation within a dynamically evolving turbulent molecular cloud formed through colliding streams of atomic gas. Such flows, sometimes called colliding flows, have been widely used to study the formation of molecular clouds because they represent a good compromise between the need to assemble the diffuse material from the large scales and the need to describe accurately 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 section 2 we present the governing equations and the physical processes such as H 2 formation.

In section 3 we describe the numerical setup that we use to perform the simulations. In section 4 we present the results of our simulations, first describing the general structure of the cloud and then focussing on the H 2 distribution. We perform a few complementary calculations which aim at understanding better the physical mechanisms at play in the simulations. Then in section 5 we compare our results with various observations which have quantified the abundance of H 2 and find very reasonable agreements. Finally in section 6 we summarize and discuss the implications of our work.

Physical processes

Governing fluid equations

We consider the usual compressive magnetohydrodynamical equations that govern the behaviour of the gas. These equations written in their conservative form are:

∂ρ ∂t + ∇ • (ρv) = 0, (1)
∂ρv ∂t + ∇ • (ρvv -BB) + ∇P = -ρ∇φ, (2)
∂E ∂t + ∇ • [(E + P)v -B(Bv)] = -ρL, (3)
∂B ∂t + ∇ • (vB -Bv) = 0, (4)
∇ 2 φ = 4πGρ (5)
where ρ, v, B, P, E, φ are, respectively, the mass density, the velocity field, the magnetic field, the total energy and the gravitational potential of the gas, L is the net loss function which describes gas cooling and heating.

Chemistry of H 2

The fluid equations [START_REF] Audit | Thermal condensation in a turbulent atomic hydrogen flow[END_REF](2)(3)(4)(5) are complemented by an equation to describe the formation of the H 2 molecules

∂n H 2 ∂t + ∇ • (n H 2 v) = k form n(n -2n H 2 ) -k ph n H 2 , ( 6
)
where n is the total hydrogen density, n H 2 represents the density of H 2 , k form and k ph represent respectively the formation and the destruction rates of H 2 .

H 2 formation

When two hydrogen atoms encounter each other in gas phase, they cannot radiate the excess of energy due to the lack of electric dipole moment, and so forth H 2 formation in gas phase is negligible.

Nowadays it is widely accepted that H 2 is formed through grain catalysis (Hollenbach & Salpeter 1971;Gould & Salpeter 1963). Hydrogen atoms can adsorb on the grain surfaces and encounter another hydrogen atom to form a H 2 molecule through two mechanisms: Langmuir-Hinshelwood, where atoms are physisorbed (efficient in the shielded gas), and Eley-Rideal, or chemisorption (efficient on warm grains) (Le Bourlot et al. 2012;Bron et al. 2014). The detailed physical description of these two mechanisms is complex and the numerical treatment is computationally expensive.

For this reason we adopted a simpler description for H 2 formation on grain surfaces, given by the following mean formation rate:

k form,0 = 3 × 10 -17 cm 3 s -1 [START_REF] Miyoshi | A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics[END_REF] this rate was first derived by Jura (1974), using Copernicus observations and confirmed by Gry et al. (2002) using FUSE observations. The formation rate depends on the local gas temperature and on the adsorption properties of the grain, so this value is corrected by the dependence on temperature and by a sticking coefficient:

k form = k form,0 T 100 K × S (T ), (8)
S (T ) is the empirical expression for the sticking coefficient of hydrogen atoms on the grain surface as described in Le Bourlot et al. (2012):

S (T ) = 1 1 + T T 2 β , (9)
where we use the same fitting values as Bron et al. (2014), namely T 2 = 464 K, and β = 1.5.

H 2 destruction

The main mechanism that destroys the H 2 molecule is the photodissociation by absorption of UV photons in the 912 -1100 Å range (Lyman and Werner transitions). In optically thin gas, in presence of a UV field of strength G 0 , H 2 is photodissociated at a rate (Draine & Bertoldi 1996):

k ph,0 = 3.3 × 10 -11 G 0 s -1 , (10)
but the H 2 gas can protect itself against photodissociation by two shielding effects. The first shielding effect is the continuous dust absorption, while the second effect is the line absorption due to other H 2 molecules, called self-shielding. Draine & Bertoldi (1996) show that the photodissociation rate can be written:

k ph = e -τ d,1000 f shield (N H 2 )k ph,0 . (11)
The first term is the effect of the dust. Here τ d,1000 = σ d,1000 N tot is the optical depth along a line of sight (σ d,1000 = 2×10 -21 cm -2 is the effective attenuation cross section for dust grains at λ = 1000 Å and N tot is the total column density of hydrogen). The second term is the self-shielding factor and we use the same approximation:

f shield = 0.965 (1 + x/b 5 ) 2 + 0.035 (1 + x) 1/2 exp(-8.5 × 10 -4 (1 + x) 1/2 ), (12)
where x = N(H 2 )/5 × 10 14 cm -2 , b 5 = b/10 5 cm s -1 , where b is the Doppler broadening parameter.

We assume that the turbulent contribution dominates, then we use b = 2 km s -1 .

Thermal Processes

For the heating and cooling of the gas, we perform the same treatment as previous work [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF]) (see also Audit & Hennebelle 2005).

The dominant heating process in the gas is the due to the ultraviolet flux from the ISRF through the photoelectric effect on grains (Bakes & Tielens 1994;Wolfire et al. 1995), where we use the effective UV field strength, calculated using our attenuation factor χ for the UV field. We also include the heating by cosmic rays (Goldsmith 2001).

The primary coolant at low temperature is the 2 P 3/2 → 2 P 1/2 [CII] fine-structure transition at 158 µm (Launay & Roueff 1977;Hayes & Nussbaumer 1984;Wolfire et al. 1995). At higher temperature, fine-structure levels of OI can be excited and thus contribute to the cooling of the gas.

We include the [OI] 63 µm and 146 µm line emission (Flower et al. 1986;[START_REF] Tielens | The Physics and Chemistry of the Interstellar Medium Valdivia[END_REF]Wolfire et al. 2003). For the Lyman α emission, that becomes dominant at high temperatures, we use the classical expression of [START_REF] Spitzer | Physical processes in the interstellar medium[END_REF]. We also include the cooling by electron recombination onto dust grains using the prescription of Wolfire et al. (1995Wolfire et al. ( , 2003)), and Bakes & Tielens (1994). Note

that strictly speaking the present cooling function does not include molecular cooling in spite of the large densities reached in the simulation. It should however be noted that it leads to temperature which are entirely reasonable even at high densities and therefore appears to be sufficiently accurate in particular given that we are not explicitly solving for the formation of molecules apart from H 2 (see Levrier et al. 2012, for a quantitative estimate).

Thermal feedback from H 2

To the atomic cooling, described before, we added the molecular cooling by H 2 lines and the heating by H 2 formation and destruction.

During the H 2 formation process, about 4.5 eV are released. The distribution of this energy into translational energy, H 2 internal energy (rotational and vibrational excitation) and into the grain heating is not well known, and the fraction of this energy that actually heats the gas varies from one author to another (Le Bourlot et al. 2012;Glover & Mac Low 2007a). We consider an equipartition of the energy, so 1/3 of the energy released goes into heating the gas:

Γ form = 2.4 × 10 -12 k form n H n erg s -1 cm -3 (13)
H 2 photodissociation provides an additional heating source. Following Black & Dalgarno (1977) and Glover & Mac Low (2007a), we assume 0.4 eV released into the gas per photodissociation: H 2 contributes to the cooling of the gas through line emission and it can become important for the diffuse medium under certain conditions (Glover & Clark 2014;Gnedin & Kravtsov 2011).

Γ ph = 6.4 × 10 -13 k ph n H 2 erg s -1 cm -3 (14)
As H 2 is a homonuclear molecule, only quadrupolar transitions are allowed and thus, the ortho and para states can be treated independently. In our simulations we adopt the cooling function from Le Bourlot et al. (1999). This is a function of temperature, total density, n(HI)/n(H 2 ) relative abundance, and of the ortho-to-para-H 2 ratio (OPR), and it considers transitions between the first 20 levels. As the cooling function is quite insensitive to the OPR, we fixed it to 3 for simplicity.

Numerical setup and initial conditions

As already discussed, colliding flows are a practical ersatz to gather matter to form a molecular cloud (Inoue & Inutsuka 2012;Audit & Hennebelle 2005;Vázquez-Semadeni et al. 2007;Heitsch Article number, page 8 of 31page.31 et [START_REF] Springel | The cosmological simulation code GADGET-2[END_REF], 2006, 2008b;Micic et al. 2013) in order to mimic large scale converging flows, as in Galactic spiral arms or in super-bubble collisions.

Article number, page 9 of 31page.31

General setup

The setup is very similar to [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF]. The size of the simulation box is L = 50 pc and it is initially uniformly filled with atomic gas, with initial number density n tot = 1 cm -3

and initial temperature T = 8000 K. Atomic gas, with the same number density and temperature, is injected from the left and right faces of the box with an average velocity V 0 = 15 km s -1 , aligned with the x axis, that is modulated by a function of amplitude = 0.5, producing a slightly turbulent profile, as in Audit & Hennebelle (2005). We use periodic boundary conditions for the remaining faces. The gas is initially uniformly magnetized, with a moderate magnetic field (∼ 2.5 µG) aligned with the inflowing gas. This configuration avoids boundary issues and facilitates the building of the cloud. Introducing an angle between the magnetic and the velocity fields would be more realistic but in this very simplified framework, requires relatively low angles (Hennebelle & Pérault 2000;Körtgen & Banerjee 2015).

Radiative transfer and dust shielding

The whole simulation is embedded in an isotropic UV field (G 0 = 1.7 in Habing units), that heats the gas. The intensity of the UV field varies from one point to another in the simulation box according to the dust shielding factor χ, which is calculated at each point of the simulation box for each timestep using our tree-based method, fully described in [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF].

Hence, the strength of the local UV field, after dust shielding, can be written G 0 = χG 0 , where χ = e -2.5A V = e -1.3×10 -21 N tot , is the mean shielding factor due to dust. This mean value is calculated using a fixed number of directions. In our previous work [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF] we have shown that the number of directions is not crucial for the dynamical evolution of the cloud, thus, for a matter of simplicity and numerical efficiency, we used 12 directions (M = 3 intervals for the polar angle and N = 4 for the azimuthal angle) which has turned out to be an excellent compromise between accuracy and efficiency.

We have adapted the code to also compute the self-shielding (stated by Eq. 12). More precisely, along each direction we compute not only the total gas density but also the H 2 column density from which we compute the shielding coefficient, f shield , and the extinction due to dust, τ d,1000 . Since the UV flux is assumed to be isotropic, the final photodissociation rate is obtained by taking the mean value over all directions. We define χ shield as

χ shield = k ph k ph,0 =< e -τ d,1000 f shield (N H 2 ) > . (15)

Numerical resolution and runs performed

To perform our simulations we employ the adaptive mesh refinement code RAMSES (Teyssier 2002;Fromang et al. 2006). RAMSES solves the MHD equations using the HLLD Riemann solvers. It preserves the nullity of the divergence of the magnetic field thanks to the use of a staggered mesh. To solve Eq. ( 6), we use operator splitting, solving first for the advection (which is identical to the conservation equation) and then subcycling to solve for the right-hand side. In AMR codes the refinement is done on a cell-by-cell basis. In our simulations, the refinement criterion is the density. When a cell reaches a given density threshold, it is split into eight smaller cells, each one having the same mass and volume. The process is repeated recursively until the maximum resolution is reached. In our fiducial run, we allow two AMR levels min = 8 and max = 10, Article number, page 11 of 31page.31 leading to an equivalent numerical resolution of 1024 3 cells, and an effective spatial resolution of about 0.05 pc. For the first refinement level (from level = 9 to = 10) the density threshold is n thresh = 50 cm -3 and for the next refinement, the density threshold is n thresh = 100 cm -3 . To investigate numerical convergence, a particularly crucial issue when chemistry is considered, we also perform a high resolution run for which the resolution is doubled (that is to say, using levels from 9 to 11 with the same refinement criterion).

We run our simulations for about 20 Myr. Note that these calculations are a little demanding because of the short timesteps induced by high temperatures, therefore it was not possible to run the high resolution simulations up to this point and this run goes up to 15 Myr. To check further for convergence, we have also performed low resolution runs (128 3 and 256 3 see Appendix) that we compare with our highest resolution simulations. We also perform complementary runs in which we modify the physics of H 2 formation in order to understand better how H 2 is formed. More precisely, in these runs we suppress H 2 formation above a certain density threshold.

Results

We now present the results of our calculations. We start first by discussing the general properties of the clouds as it is essential to understand in which context the hydrogen molecules form. Then we discuss the abundance and distribution of H 2 itself. Finally, we perform various complementary calculations aiming at better understanding in which conditions H 2 forms.

Qualitative description of the cloud

As many colliding flow calculations have been performed during the last decade, we do not attempt here to display all the steps that the flow is experiencing (see references above for accurate descriptions) and quickly summarize the main steps.

First, the WNM flows collide and this triggers the transition for a fraction of the gas into moderately dense gas1 , that is to say densities on the order of 100-1000 cm -3 as a consequence of the thermal instability and ram pressure. Second, as enough gas has been accumulated, gravity becomes important and triggers infall, first at small scale and then at larger ones. This leads to a continuous increase of the dense gas fraction.

Left panels of Figs. 1 and2 show respectively the column density along the z-direction and a density cut in the xy plane at time 20 Myr. As will be seen below, this corresponds to a phase where gravity is already playing a significant role. As can be seen from the density cut, the flow is very fragmented. The dense gas is organised in dense clumps which are embedded in diffuse and warm gas. The column density spans typically 2-3 orders of magnitude from a few 10 21 cm -2 to at least 10 24 cm -2 in some very compact regions, the dense cores where gravity is triggering strong infall.

The column density map is also clumpy although less obviously structured than the density cut. Article number, page 13 of 31page.31

Density PDF and volume filling factor

In many respects the density distribution is an essential cloud property which reveals the dynamical state of the cloud and strongly influences the formation of H 2 .

Figure 3 shows the density probability distribution function (PDF) at time 5, 10, 15 and 20 Myr.

At time 5 Myr, the PDF presents two peaks, one at low density ( 1 cm -3 ) and one at a few 100 cm -3 . While the first one is simply a consequence of the initial and boundary conditions (which inject WNM from the boundaries), the second one represents the cold phase, whose formation has been triggered by the thermal instability. Because of the supersonic turbulence that develops in the cold gas, the density PDF is broad and presents a lognormal shape (Kritsuk et al. 2007;Federrath et al. 2008). At later times, more cold gas accumulates and the PDF broadens. It is seen that the high density tail tends to become less steep. At intermediate densities (ρ ∼ 10 3 -10 4 cm -3 ) the shape is compatible with a power-law with an exponent between -1 and -2. This is typical of what has been found in super-sonic turbulent simulations which include gravity (Kritsuk et al. 2011). At very high densities (> 10 4 -10 5 cm -3 ), the density field flattens. This is because, as we are not using sink particles, the very dense gas piles up and accumulate within few grid cells. This feature is thus a numerical artefact and represents a limit of the simulations.

Figure 4 shows the filling factor per density bins and for different computational box regions at time 20 Myr. The four lines shows the results for four different box regions (as shown in the label). The blue line shows the whole box while the cyan line is limited to the densest regions of the computational box. The volume filling factor is clearly dominated by the warm and diffuse gas.

Selecting gas of densities n < 3 cm -3 , we find that it typically occupies 70% for the region located between x = 27 and x = 37 (85% for the whole computational box). The dense gas, even in the densest region, occupies only a tiny fraction, for example we find that the gas denser than 100 cm -3 has a filling factor of about 3% (respectively 1.5% for the whole box). The remaining 26% (13%) are filled with gas of densities between 3 and 100 cm -3 . Thus we conclude that the interclump medium which occupies most of the volume, is itself made with two components, a warm gas that is similar to the standard WNM but can be a bit denser and a cold gas that is similar to the CNM but contains, as will be seen later, a significant fraction of H 2 . They typically occupy about 70% and 25% percent of the volume respectively.

Article number, page 15 of 31page.31 Another important diagnostic to characterise the dynamics of a medium are the different pressures.

Figure 5 shows the thermal, P th , dynamical, P ram , and magnetic, P mag , pressures equal to respectively P th = knT , P ram = ρV 2 and P mag = B 2 /8π.

As can be seen the dynamical pressure dominates the thermal one by typically one to two orders of magnitude in the dense gas (n > 10 cm -3 ) while they are more comparable in the diffuse gas.

Article number, page 16 of 31page.31 The magnetic pressure lies inbetween the two, with values a few times larger than the thermal one.

Also P ram and P mag increase with density and at densities on the order of n 10 3 cm -3 , they are about one order of magnitude larger than their mean values at n 10 cm -3 . Therefore while the low density gas which fills the volume provides some confining pressure, it has a limited influence on the clumps.

Altogether these results indicate that the molecular cloud produced in this calculation can be described by density fluctuations or clumps, induced by both ram pressure and gravity. The clumps occupies a tiny fraction of the volume which is filled by a mixture of warm diffuse and dense HI gas (occupying respectively a fraction >70% and >20% of the volume). This low density material is feeding the clumps, which grow in mass. This picture is in good agreement with the measurement performed by Williams et al. (1995), who found that the interclump medium has typically a density of a few particles per cm -3 and a large velocity dispersion of several km s -1 .

Statistics of clumps

As the cloud is organised in clumps, we further quantify the simulations by providing some statistics. To identify the clumps, we use a density threshold of 1000 cm -3 and a friend of friend algorithm. This structure is important for the chemistry evolution since, at the edge of clumps very significant density and temperature gradients arise. Fig. 6 shows the mass spectrum, the velocity dispersion and the kinetic β parameter (that is to say P ram /P mag ) of the clumps. As can be seen the mass spectrum for masses above a few 0.1 M presents a powerlaw with an index of about -1.7, which is very similar to what has been shown in related works (e.g. Audit & Hennebelle 2010;Heitsch et al. 2008b). The velocity dispersion within clumps as a function of their size is about σ 1km s -1 (R/1pc), which is close to the Larson relation (Larson 1981;Hennebelle & Falgarone Article number, page 17 of 31page.31 2012). From third panel, it is seen that the kinetic β parameter is typically of the order of, or slightly below 1, showing that the magnetic field plays an important role within the clumps.

Molecular hydrogen

We now turn to the H 2 abundance discussion. Middle and right panels of Fig. 1 show the H 2 column density and molecular abundance, f (H 2 ) = 2n(H 2 )/n. As expected the high column density regions are dominated by the H 2 molecules and values of f (H 2 ) close to 1 are obtained there. Obviously in the outer parts of the cloud, the values of f (H 2 ) are significantly lower.

The UV screening factor

Middle and right panels of Fig. 2 respectively shows a cut of n(H 2 ) and the value of the shielding parameter, χ shield = k ph /k ph,0 . This latter steeply decreases when entering inside the cloud where it takes values on the order of 10 -3 . A comparison between the density cut with k ph /k ph,0 reveals that there are regions of diffuse material in which the shielding paremeter is low. This is because these regions are surrounded by dense gas in which H 2 is abundant therefore providing an efficient self-shielding. As a consequence, there are low density regions which contain a relatively high abundance of H 2 (see left and middle panels of Fig. 2).

A more quantitative statement can be drawn from Fig. 7 that displays the distribution of the dust shielding, χ (top panel) and of the total shielding χ shield as a function of density (bottom panel). As expected most low density cells, are associated to χ values that are close to 1, that is to say are weakly shielded. There are however a fraction of them which present values of χ shield as small as 10 -5 . At higher densities, the mean value of χ decreases with a steep drop between n 10 3 and 10 4 cm -3 which correspond to the point where the dust significantly absorbed the UV external field. For the total shielding, a steep drop is observed at about 10 cm -3 . There is however a broad scatter, which increases from n = 1 to n = 10 4 cm -3 for the dust shielding and tends to decrease for the total shielding χ shield . This indicates that the mean value of χ shield is not a sufficient information to quantify the abundance of H 2 expected at a given density. This is a clear consequence of the complex cloud structure. Since χ shield plays a key role for the formation of H 2 , this complex distribution certainly makes molecular hydrogen formation in a multi-phase turbulent cloud, a complicated issue.

Global evolution of molecular hydrogen

Figure 8 shows f (H 2 ), the mean molecular fraction within the whole cloud, as a function of time.

During the first 5 Myr, f (H 2 ) exponentially increases from nearly 0 to about 0.1. After this phase, f (H 2 ) keeps increasing in a more steady almost linear way. Near 15 Myr, about 40% of the gas is in H 2 molecules. Since, as discussed before, the cloud is rather inhomogeneous in density and temperature, we have also plotted the time evolution per bins of density and temperature. The corresponding curves are displayed in Fig. 9. As expected, the total distribution closely follows Article number, page 18 of 31page.31 the values of the densest bins, that is to say corresponding to densities larger than 100-1000 cm -3 .

For these densities, the timescale for H 2 formation is thus of a few ( 6 -8) Myr. As recalled previously, this typical timescale is about 10 9 /n yr. Therefore the timescale we observe in the simulation agrees well with this scaling. Altogether, this is in good agreement with the results of Glover & Mac Low (2007a,b) which found that H 2 could be formed quickly in molecular clouds because it preferentially forms into clumps that are significantly denser than the rest of the cloud.

In the present case, the mean density of our cloud at time 5-10 Myr would be about 10-100 cm -3

(depending exactly which area is selected) and this would lead to timescale for H 2 formation on the order of 10 to 100 Myr.

What is somewhat more surprising is the relatively short amount of time ( 10 Myrs) after which there is a significant amount of H 2 in low density gas (n 1 -10 cm -3 ). The expected timescale for H 2 formation in this density range would be about 10 8 Myr. As can be seen from Fig. 9, the behaviour of f (H 2 ) is similar, although not identical, to the behaviour at higher densities.

There is a fast increase followed by a slowly increasing phase. There is however a change of slope at about 4 Myr that is different from the higher density evolution. The steady evolution starts at about 8 Myr instead of 5-6 Myr for the denser regions. This suggests that the presence of H 2 at low density is triggered by formation of H 2 in the dense gas. This can happen in various non-exclusive ways. First of all, some dense gas can expand back and mix with diffuse gas. This may happen either through a sonic expansion, through evaporation or through numerical diffusion. Second of all, some of the diffuse gas may be surrounded by dense gas in which f (H 2 ) is large and therefore has a low χ shield . Below we attempt to analyse further these mechanisms.

Finally, we note from the bottom panel of Fig. 9 that there is a small (1-10%) but nevertheless non negligible fraction of H 2 within relatively high temperature gas, i.e. T > 300 K. Since H 2 is the first molecule involved to produce other molecules and since some of them, such as CH + , are formed through reactions with high activation barriers (on the order of 4640 K for CH + ), the presence of H 2 could have consequences for the production of these species as proposed in Lesaffre et al. (2007). In the same way, as discussed below, this warm H 2 contains molecules in excited high J rotational levels, which can therefore contribute to the gas cooling.

Detailed distribution of molecular hydrogen

To further quantify the distribution of H 2 in the cloud, we display histograms of f (H 2 ) per bin of density (Fig. 10). We also draw in the same panel the distribution f (H 2 ) obtained at equilibrium.

Knowing the density, n, the temperature, T and the shielding factor, χ shield , we solve Eq. ( 6) at equilibrium. Note that this distribution is not fully self-consistent since the value of χ shield should in principle be recalculated to be consistent with this equilibrium distribution (this would imply performing several iterations). In practice, the goal here is simply to compare with the time-dependent distribution to gain insight on the H 2 production mechanisms, and it is therefore easier to have exactly the same formation and destruction rates. As time goes on, the differences between the equilibrium and time-dependent distributions become less important. For example the two distributions are obviously closer at time 20 Myr (right panels) than at time 5 Myr. However, they are not identical. Since all distributions at time 15 Myr and 20 Myr are similar, the persistence of the differences between the equilibrium and time-dependent distributions, indicates that this is likely the result of a stationnary situation. Most likely it reflects the accretion process, that is to say HI gas (possibly mixed with a fraction of H 2 ) is continuously being accreted within denser clumps and therefore the dense gas does not become fully molecular.

This interpretation is also consistent with the fact that the effect is more pronounced for the fourth density bin (line 4) than for the fifth and highest one (line 5). Looking now at the lowest density bin (first line), one finds the reverse effect at time 10 and 15 Myr, the time-dependent distribution dominates the equilibrium one. This is most likely an effect of the turbulent diffusion. From the low mass contained in this low density bin, this however remains a limited effect.

The second density bin (line 2), which corresponds to n between 1 and 10 cm -3 is a little puzzling. No significant differences are seen between the two distributions, which is surprising.

One possibility is that various effects compensate. That is to say the time delay to form H 2 (clearly visible for the 3 denser density bins) may be compensated by an enrichment from the denser gas. This is indeed possibly happening for the second density bin at time 20 Myr (fourth column, second line) where it is seen that there is overall a small excess of time-dependent H 2 but that for f (H 2 ) > 0.7 the equilibrium distribution dominates.

Spatial fluctuations

As a complement to the time evolution of f (H 2 ), it is also worth to investigate the spatial fluctuations at a given time. In this respect, the clumps constitute natural entities to study the spatial variations. The bottom panel of Fig. 6 shows f (H 2 ) within the clumps. For the most massive clumps, the value of f (H 2 ) is close to 1 and the dispersion remains weak. On the contrary, for the less Article number, page 20 of 31page.31 massive clumps, the dispersion becomes very significant and the value of f (H 2 ) can be in some circumstances rather low. This shows that clump histories, that is to say their ages and the local UV flux in which they grow, has a major influence on f (H 2 ).

Further analysis of H 2 formation

In order to better quantify the interdependence of the different density bins, we have conducted three complementary low resolution runs for which the formation of H 2 is suppressed (k form is set to zero) above a given density threshold. Figure 11 displays the result. In the top panel, no threshold is applied. In the middle and bottom panels a threshold of 1000 cm -3 and 100 cm -3 , respectively, is applied. As can be seen with the threshold 1000 cm -3 , the values of f (H 2 ) are not significantly changed except for the highest density range (for n > 1000 cm -3 ). On the other hand, with a threshold 100 cm -3 , the values of f (H 2 ) decrease by a factor of about 3 for all density bins. This clearly shows that most of the H 2 molecules in the low density gas form at a density of a few 100 cm -3 . While the shielding provided by molecules in gas denser than this value could contribute to enhance the more diffuse gas, the filling factor of the dense gas is too small to affect significantly the diffuse gas that has a much larger filling factor (see Fig. 4). Therefore we conclude that the H 2 abundance within the low density gas (<100 cm -3 ), is likely a consequence of turbulent mixing and gas exchange between diffuse and dense gas.

Thermal balance

To quantify the influence of the H 2 molecule on thermodynamics, we present the contributions of the various heating and cooling terms to the thermal balance as a function of density. Figure 12 shows the standard ISM cooling (blue curve) as well as the cooling by H 2 (green curve). As can be seen, the latter remains overall below the former except at density of about 3-5 cm -3 where they become comparable. While the heating due to H 2 destruction remains negligible at any density (light red curve), the heating due to its formation (red curve) turns out to be the dominant heating mechanism at high densities. It however does not affect very significantly the temperature which stays equal to about 10 K.

Finally, as can be seen for densities between 1 to 50 cm -3 , the cooling dominates the heating by a factor of a few. This indicates that there is another source of heating equal to the difference, which turns out to be due to the mechanical energy dissipation. This latter therefore appears to have a contribution comparable to the UV heating. This explains why warm gas is actually able to survive within molecular clouds. In the same way that density is much larger than in the rest of the ISM, kinetic energy is also larger and provides a significant heating (e.g. Hennebelle & Inutsuka 2006).

Article number, page 21 of 31page.31 (Savage et al. 1977), and FUSE (Gillmon et al. 2006;Rachford et al. 2002Rachford et al. , 2009)).

Confrontation to observations

We now compare our results with various observations which fall in two categories. The first ones have attempted to measure the molecular fraction, i.e. the total column density of H 2 with respect to the total column density. The second ones have been measuring the excited levels of H 2 , therefore presumably tracing high temperature gas. These two sets of observations are therefore complementary and very informative.

Molecular fraction vs column density

The H 2 column density has been estimated along several lines of sight mainly (though not exclusively) with the Copernicus (Savage et al. 1977) and FUSE observatory (Gillmon et al. 2006;Rachford et al. 2002Rachford et al. , 2009)). These observations are probing lines of sight spanning a wide range of column densities and therefore constitute a good test for our simulations although a difference to keep in mind is that the lines of sight extracted from our simulations are all taken from a 50 pc region and therefore are more homogeneous than the observed lines of sight.

Figure 13 shows the molecular fraction f (H 2 ) as a function of total column density for all lines of sight of the simulations (taken along the z-direction) and the available lines of sight reported in (Gillmon et al. 2006;Rachford et al. 2002Rachford et al. , 2009)). As can be seen, the simulation results and the observations agree rather well, many observational data directly fall into the same regions than simulated data. In particular, the two regimes, i.e. the vertical branch at column density between 10 20 and 3 × 10 20 cm -3 as well as the denser column density area, are clearly seen both in observations and in the simulation.

This said, there are also data points which are not reproduced by any of the lines of sight from the simulations. This is particularly the case for column densities larger than 3 × 10 21 cm -2 and for the Copernicus survey at column densities around 10 21 cm -2 . The most likely explanation is that Article number, page 22 of 31page.31 This is beyond the scope of the present paper.

Excited levels of H 2

We now turn to the comparison with the observations of rotational levels of the H 2 molecule, which require high temperatures to be excited.

As investigated by previous authors (Lacour et al. 2005;Godard et al. 2009), these data cannot be explained for example by UV excitation. While Lacour et al. (2005) concluded that to explain the observations, a warm layer associated to the cold gas should be present, Godard et al. (2009) perform detailed models entailling shocks or vortices. The general idea is that in such dissipative Article number, page 23 of 31page.31 small scale structures, the temperature is reaching high values because of the intense mechanical heating. Interestingly, Godard et al. (2009) obtained nice agreement between their model and the data. Let us stress that the small dissipative scales which are determinant in these models, cannot be described in our simulations which have a limited resolution on the order of 0.1 pc. On the otherhand as our simulations contained large quantities of warm H 2 , it is worth investigating whether the inferred populations of excited H 2 can be quantitatively reproduced.

Calculation of population levels and column densities

We calculate the population of the first six rotational levels (J) of H 2 as in Flower et al. (1986), based on the Elitzur & Watson (1978) approach. This calculation is done in post-processing for all grid cells. It needs the total hydrogen density, the H 2 density, the He density (we assume 10% of total H), the temperature, and the ortho-to-para ratio. It supposes thermal equilibrium, and it includes spontaneous de-excitation, and collisional excitation by HI, by He, and by other H 2 molecules.

Another important parameter is the ortho-to-para ratio (OPR), which is not very well known in the ISM (e.g. Le Bourlot 2000) and may need a specific time-dependent treatment, which is beyond the scope of the paper. For this reason, we will consider two different simple assumptions. First we will assume an OPR equal to 3, which simply corresponds to the statistical weights, and second we will assume that the OPR is given by thermal equilibrium. These two assumptions likely represents two extreme regimes and the exact value likely lies inbetween.

Once we obtain the density of the H 2 rotational levels, column density maps are produced by integration along the x-direction.

Comparison with observed column densities

Figure 14 shows the distribution of the column densities corresponding to the various rotational levels in the simulation. The top panel shows results for OPR equal to 3 while the bottom one shows results for the thermal equilibrium assumption. The color coding shows the total column Article number, page 24 of 31page.31 The simulation results split into 2 groups of points, very low column densities coming from the WNM that is located outside the molecular cloud (blue points) and higher column densities

Article number, page 25 of 31page.31

(green and yellow points), which come from the gas belonging to the molecular cloud. For the J = 0 and 1 levels, the column densities are more or less proportional to the total column density. This is indeed expected since most of H 2 molecules are in one of these two levels. Their respective abundance significantly depends on the assumption used for the OPR. While there is about a factor 10 between the first two rotational levels for N J /g J when the OPR is equal to 3, it is almost a factor thousand when the OPR is assumed to be at thermal equilibrium.

The situation is different for the higher levels. The largest column densities of the excited rotational levels do not correspond to the largest total hydrogen column densities and generally there is no obvious correlation between the two. This is due to the fact that the high J levels are coming from the warm gas (with temperatures between few hundreds and few thousands Kelvin), which itself has a low column density and is largely independent of the cold gas.

The comparison with the observational data is very enlighting. Overall the agreement with the simulation results is very good. The observational data points have column densities which are very similar to the simulation data except for J = 1 for which the simulation data points are above the observational ones for an OPR equal to 3 and below for an OPR at thermal equilibrium. This probably suggests that the actual OPR is lying inbetween. There is a possible tension with some of the data of Gry et al. (2002) and Lacour et al. (2005) for J = 2, by a factor on the order of 3-10. If confirmed this would indicate that another mechanism, such as the one proposed by Godard et al. (2009) could operate and contribute to excite H 2 .

In order to understand what is the origin of the excited rotational levels in the simulation, Fig. 15 displays the mean density of the H 2 excited levels for five density bins as well as the complete distribution. As expected, for J = 0 the main contributor is the dense gas. For the excited level the highest contributions are found in gas of total density between 1 and 10 cm -3 . Then, in roughly equal proportions, very diffuse gas with n < 1 cm -3 and moderately dense one (with n between 10 and 100 cm -3 ). It is worth recalling that the peak of the cooling contribution of H 2 , precisely lies in this density domain (see Fig. 12).

Velocity dispersion of individual rotational levels

So far the conclusion we can draw is that the population of the excited levels of H 2 (J ≥ 2) is dominated by the warm H 2 that is inspersed between the cold and dense molecular clumps.

This "layer" of warm gas is directly associated with the molecular cloud and agrees well with the proposition made by Lacour et al. (2005). Interestingly these authors made an other interesting observation. They found that the width of the lines associated to the excited levels increases with J. In order to investigate whether our simulation exhibits the same trends, we have calculated the mean velocity dispersion along several (125 along each axis) lines of sight for each excited level. Figure 16 shows the mean dispersion parameter, b linked to the velocity dispersion σ by Article number, page 26 of 31page.31

b 2 = 2 1/2 σ 2 , and where

σ 2 = ρ (v i -< v i >) 2 + C 2 s dx ρdx . ( 16
)
Figure 16 displays the results. Since the colliding flow configuration is highly non-isotropic, we have calculated these velocity dispersion along x-axis (top panel), y-axis (middle panel) and z-axis (bottom panel). While a more accurate treatment would consist in simulating the line profiles and then applying the same algorithm that the authors used, this simple estimate is already illustrative.

From σ 2 , we can infer the width of the line and the parameter, b.

As can be seen the trends in the simulation and in the observations are similar. Higher J levels tends to be associated with larger b. Quantitatively the agreement is satisfying for the y and zdirections. As the observations, the simulated lines of sight present slightly larger velocity dipersion for higher J. The agreement is poor for the lines of sight along the x-direction as they present a dispersion that is too high with respect to the observations. This is most likely an artefact of the colliding flow configuration.

The trend of larger b for higher J stems for the fact that higher levels need warmer gas to be excited. Therefore the fluid elements, which are enriched in high rotational levels have larger temperatures and therefore larger velocity dispersions (since both are usually correlated).

Altogether these results suggest that the excited rotational level abundances reveal the complex structure of molecular clouds that is 2-phase in nature and entails warm gas deeply interspersed with the cold gas due to the mixing induced by turbulence.

Conclusions

We have performed high resolution magneto-hydrodynamical simulations to describe the formation of molecular clouds out of diffuse atomic hydrogen streams. We particularly focus on the formation of the H 2 molecule itself using a tree-based approach to evaluate UV shielding (see [START_REF] Valdivia | A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds[END_REF]).

In accordance with previous works (Glover & Mac Low 2007a), we find that H 2 is able to form much faster that simple estimates, based on cloud mean density, would predict. This is because, since the clouds are supersonic and have a 2-phase structure, H 2 is produced in clumps which are much denser than the clouds on average.

Due to a combination of phase exchanges and high UV screening deep inside the multiphase molecular clouds (numerical convergence tests suggest that numerical diffusion is not responsible of this process), a significant fraction of H 2 develops even in the low density and warm interclump medium. This warm H 2 contributes to the thermal balance of the gas and in the range of density 3-10 cm -3 its cooling rate is comparable to the standard cooling of the ISM.

Article number, page 27 of 31page.31

Detailed comparisons with Copernicus and FUSE observations show good agreement overall.

In particular the fraction of H 2 varies with the total gas column density in a very similar way showing a steep increase between column densities 10 20 and 3×10 20 cm -3 and a slow increase at larger column densities. There is however a trend for the high column densities regions to present H 2 fractions significantly below the values inferred from the simulations. This is a possible consequence of the constant UV flux assumed in this work. Interestingly, the column densities of the excited rotational levels obtained at thermal equilibrium reproduce fairly well the observations. This is a direct consequence of the presence of H 2 molecules in the warm interclump medium and suggests that the H 2 populations in excited levels reveal the 2-phase structure of molecular clouds. this suggests that numerical diffusion is not primary responsible of the fraction of warm H 2 that we observe at low density (as discussed in Section 4.2). To understand the effect of the total shielding, dust shielding for the UV and self-shielding due to other H 2 molecules, as the influence of the magnetic field and the influence of diffusivity of the Riemann solver we have performed additional simulations.

5.A Influence of the total shielding

We have performed two identical MHD simulations to understand the influence of the total shielding for the photodissociation χ shield on the formation of molecular hydrogen. The set-up is the same as the one used in the previous section, but for this test we use a lower resolution ( 0.2 pc). Both simulations include the shielding for the UV radiation in the thermal balance, but only the first one includes the total shielding factor χ shield for the photodissociation of H 2 .

Figure 5.1 shows the evolution of the molecular fraction during simulations. When χ shield is not included the presence of molecular gas is globally suppressed, and in a period of 15 Myrs the simulation converts less than 10% of total gas into molecular gas.

We find that when the photodissociation rate includes the combined shielding effects, the amount of gas in H 2 form is about 5 times higher than when this effect is not included. An interesting result comes from the most diffuse gas, for which the molecular fraction when the shielding is included is about 2 orders of magnitude higher than the case without shielding. This indicates that the shielding effect of dust and H 2 self-shielding is crucial to explain the warm H 2 . It is well known that MHD simulations produce clouds with a different morphology than pure hydro simulations (Hennebelle et al. 2008). We have performed two simulations using the standard set up described in Section 5.1, at the same resolution. The first one is a pure hydro simulation ( B = 0), and the second one is a MHD simulation, with a moderate magnetic field aligned with the gas inflow (B ∼ 2.5 µG).

5.B Influence of magnetic field

Figure 5.3 shows the general evolution of the cloud in terms of total column density. The hydro case produces higher column densities at earlier times, as a consequence of the lack of magnetic support. In the MHD simulation, magnetic fields offer an additional support against gravity, delaying the collapse, as shown by Hennebelle et al. (2008). The cloud formed in the hydro simulations displays a more compact structure, reaching higher densities that favor an earlier forma-tion of molecular gas (Fig. 5.4). Higher column densities favor higher molecular fractions, as shown in Fig. 5.5. An interesting remark is that MHD simulations produce clumpy and filamentary structures, allowing a mixture of warm and cold gas that produces an interwoven structure with different amounts of H 2 .

The evolution of the amount of gas in molecular form is shown in Fig. 5.2. The general evolution of the total amount of H2 (top panel) is very similar in both cases, but the hydro case reaches 10% of molecular hydrogen about 1 Myr before than the MHD case. The bottom panel shows the molecular fraction evolution for different density bins. The difference between both cases is more pronounced at densities of few ∼ 1 -10 cm -3 . As the amount of molecular gas depends on the ability of gas to protect itself against photodestruction, we analize the distribution of the total shielding factor during the simulation. Fig. 5.7 shows that in the hydro case more gas is found at higher densities, in regions that are well shielded, which is consistent with an earlier formation of molecular gas for the hydro case.

As a complement to the column density plots we present results for a slice of gas in Fig. 5.6. This figure shows the total number density map, temperature map, and total shielding factor χ shield map for a slice cut through the midplane of the simulation box at 15 Myrs. These three maps give a detailed view of the local structure and the effect of shielding. The hydro case shows clumpy structures, for which the shielding factor can reach values as low as χ shield ∼ 10 -5 , but gas is not as intermixed as the MHD structure. Another noteworthy feature in the plots for the MHD case is a hot and low density bubble near x = 30 pc, y = 40 pc, for which the photodissociation rate is reduced to 1% of the value for similar regions, in terms of density and temperature conditions, that are more exposed to the UV field, revealing the importance of the global structure on the survival of molecular gas. This can be seen in Fig. 5.7, that displays the total shielding factor χ shield distribution. This figure shows how more cells in the hydro simulation reach higher shieldings (lower χ shield ) at earlier times, and in general more cells are found at high densities and low χ shield .

Figure 5.8 shows the distribution of the molecular fraction as a function of the total density. Both cases present a vertical transition separating both distributions at low molecular fraction and at high molecular fraction, but their shapes is very different. The hydro case present a compact distribution and a much narrower transition that the MHD case. We interpret this as the outcome of a much uniform physical conditions, specially a much uniform mass distribution. As the cloud evolves toward a more complicate structure, the distribution becomes wider. Finally, Figs. 5.9 and 5.10 show the populations of the first 6 rotational levels, and the Doppler broadening parameter b (in the z-axis). Both simulations are in good agreement with observations, but the hydro case seems to reproduce better the populations of levels J = 2 and J = 3. This difference could be caused by the additional constrain imposed by the magnetic field. Fig. 5.10 shows that in general the b parameter is higher in the MHD case. This can be explained as a consequence of the difference in the global structure and velocity field. It is important to note that as the lines-of-sight are defined by a fixed grid, not all the lines of sight are represented. As the MHD structure is more clumpy, it is easier to miss lines-of-sigh with high densities, on the other hand, the hydro structure is more uniform and denser, so high column densities are more easily represented.

5.C Influence of the solver

To understand the influence of the solver we compare the HLL [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF], and the HLLD [START_REF] Miyoshi | A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics[END_REF] Riemann solvers (see Section 3.3.4). We use the standard set up described in Section 5.1, and we use the output at 10 Myr as the starting point. We perform two identical simulations except for the solver. We let the simulations evolve during 5 Myr, until they reach 15 Myr.

Figure 5.11 shows the total number density, the number density of molecular hydrogen, and the temperature for a slice cut through the midplane of the simulation box (z = 25 pc). Even though the general structure is very similar, the higher diffusivity of the HLL solver produces noticeable differences between both simulations. Total number density plots show structures that are more diffuse for the HLL simulation than for the HLLD simulation. In the same manner, the H 2 number density plots show that molecular hydrogen is concentrated in dense clumps in the HLLD simulation, while the H 2 distribution is much smoother in the HLL case (this can be seen specially in the region toward x = 30 pc, y = 20 pc. The temperature distribution is also modified, in particular hot/warm bubbles seen in the HLLD case seem to be diluted in the HLL simulation, having consequences on the physical state of gas, and thus on the determination of the determination of excited populations of H 2 . Figure 5.12 mass in H 2 form per density bin shows a deviation between the HLL simulation and HLLD simulation. HLL simulation show less gas in H 2 form at higher densities, and systematically for all density bins it presents less gas at higher molecular fraction and an excess at lower molecular fractions. This can be interpreted as a consequence of the higher diffusivity of the HLL Riemann solver.

128

6

Conclusions and Perspectives

The work presented in this thesis focused on the influence of the ultraviolet radiation and the chemistry on the formation and evolution of molecular clouds, but also, and I would say more importantly, on the influence of the structure of the molecular cloud and its evolution on the chemistry and on the effective radiation field. At the begining of this thesis we presented a global view of the different physical processes that influence the evolution of molecular clouds, emphasizing their interdependence. In adition to the thermal balance, we paid particular attention to the role of turbulence, magnetic field, and radiation field, which are major players in the physics of the ISM. In order to understand the relative importance of different processes it has been mandatory to count with accurate models for the thermal balance in the interstellar medium, and for the formation of molecular gas, as well as numerical tools adapted to represent all the relevant scales at play in molecular clouds.

One of the first challenges was the correct estimate of the amout of radiation that reached the gas at any point. As dust absorbs ultraviolet photons, and reemits this energy at lower frequencies (in the infrared region) that can escape the gas easily, less UV photons are available to heat the gas through photoelectric heating. Then, in order to include this effect it is needed to estimate column densities in any possible direction for any element in the simulation box. As the cloud evolves in time, this value changes with it, so it is necessary to recalculate column densities at each timestep. The main complexity of this problem lies on the onerous cost of this calculation. Each cell in the simulation box has its own "illumination", that varies from point to point, and in time. One of the first results of this work was the development of the Tree-based method for AMR codes, detailed in Chapter 4, which is able to estimate column densities and screening factors on-the-fly in order to count with the right amount of UV photons able to heat the gas. The optimization of the method was based on the self-similarity of octs at different resolution levels and on the pre-calculation of geometrical corrections, making the calculation up to three times faster than the case where the calculation was done in situ, and keeping a good accuracy. accurate estimates for the self-shielding. As the Tree-based method is a general method that permits us to obtain column densities, it was modified in order to estimate H 2 column densities, that used along with total column densities permitted the estimate of photodissociation rates including the effect of the UV shielding by dust, and the self-shielding by other H 2 molecules. We found that this shielding effect is of key importance for the survival of gas. In this part of the work we concluded that H 2 formation occurs faster than predicted by simple estimates based on mean density of gas. Turbulence drives local enhancements in density, that accelerates the formation of molecular gas, while the global structure of the cloud protects it. We have shown the presence of fairly warm molecular gas in the relatively diffuse and warm interclump medium, that was likely formed at higher densities and migrated toward less dense gas through the phase mixing operating in turbulent clouds.

Excited H 2 has been observed by Copernicus and FUSE satellites in several lines of sight, but its presence remains unexplained. We analyzed whether excited H 2 could be explained by the local physical conditions that reign in molecular clouds. For this reason we have calculated the population distribution of H 2 for the first 20 levels, and we compared the first 6 levels to the populations observed with FUSE. We used two prescriptions for the ortho-to-para ratio that represents two cases: the standard ratio OP R = 3, and at thermal equilibrium OP R(T ). Both prescriptions are equivalent at high temperatures, so the first two levels are the most sensitives to variations in the OP R. We have found that the observed OP Rs lie in between of both prescriptions, indicating that they are rather a mixture of local values of this ratio. The analysis of excited levels turned out to be more interesting. We found a good agreement between simulations and observations for different lines of sight. In particular excitation diagrams are well reproduced by simulations, as well as the Doppler broadening parameter b. The analysis of the mean density of H 2 excited levels showed that populations arise from regions with different physical conditions. These results suggest that the excited populations observed toward many lines of sight would arise as a consequence of the multi-phasic nature of molecular clouds, where warm and cold phases are mixed. The observation of these levels would prove the multi-phase nature of interstellar gas.

These results open new perspectives on the interpretation of observations. As H 2 formation is an important prerequisite for more complex molecule formation, the presence of fairly warm H 2 , explained through a mixture of local density enhancement, shielding effects and turbulent mixing, with excited populations arising as a direct consequence of the phase mixture of turbulent molecular clouds, would offer a new simple explanation to the presence of other excited molecules. To this purpose it is needed to develop more detailed chemical networks that includes the formation of molecules such as CH + and H + 3 . Two aproaches are possible: the first one is the inclusion of a modest, but more detailed network in Ramses, in order to solve dynamically the evolution of different species, while, a second approach would be the calculation in postprocessing of the abundances at equilibrium. Both approaches have been udertaken with different degrees of progress. For the first approach, we have included a reduced version of the chemical network used in Chemses, which is the merging of the code Dumses (Ramses for the Dummies, without the AMR) and the shock code Paris-Durham. We have been able to extract the solver and merge it to the modified Ramses code version. An exploratory work that does not uses the shielding estimates, have been performed using the Chemses network. We found that simulation times are slower by two orders of magnitude. The next steps will be the coupling between the chemistry and the estimates for shielding, and the obtainment of a minimal network. To accomplish this last task two reduced networks have been proposed, the first one containing only hydrogen chemistry, and the second one extending the network to OI, OII, C, CII, SI, SII. The first tests have not been conclusive, and further tests are needed. The second approach is not totally separated from the first one, but for some species it is independent. In the current state of this approach, the equations for the calculation of the abundances of CH + and H +

1. 1

 1 Flow of baryons in the Milky Way, from Draine (2011). . . . . . 1.2 Schematic energy level structure for CII and OI, showing the transition probabilities A ul for the main cooling transitions. . . . 1.3 Standard equilibrium curve: Thermal pressure P/k vs total hydrogen density (left panel), and cooling (solid lines) and heating (dashed lines) processes in the ISM (right panel), from Wolfire et al. (1995). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 The interstellar cooling function for different ionization fraction

  the same work, but among the most interesting regions we can recall the Taurus Molecular Cloud Complex (T), the ρ Ophiuchus Complex (ρ), and the Orion Complex (O). . . . . . . . . . . . . . . . . . . . . . 1.6 Eddies breaking into smaller eddies. Sketch of obstacles in a water flow (on the top) and a free water jet falling issuing from a square hole into a pool. Drawing by Leonardo da Vinci, circa 1500. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Kolmogorov cascade showing the inertial range and the k -5/3 law. 1.8 The Galactic magnetic field as revealed by Planck. From Juan Soler webstite http://www.ias.u-psud.fr/soler/planckhighlights. html. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Electromagnetic wave is rotated under the influence of the magnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Polarized starlight Heiles and Crutcher (2005) and dust thermal emission (total intensity and vectors) from Planck. J-Ph. Bernard, Planck collaboration (ESLAB 2013 ). . . . . . . . . . . III 1.11 Left panel: LOFAR stations in Europe. Right panel: Total magnetic field of the Whirpool galaxy M 51. Image from Mulcahy et al. (2014) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Interstellar radiation field showing the main components. Figure from Draine (2011). . . . . . . . . . . . . . . . . . . . . . . . . .

  Figure from Tasker et al. (2008) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The top panel shows a photograph of the Antennae galaxies, one of the most outstanding example of a galaxy collision, obtained with one of the NOAO telescopes from the ground. Credit: NOAO/AURA/NSF, B. Twardy, B. Twardy, and A. Block (NOAO)).

IV 3 . 6

 36 Wave fan for MHD waves. . . . . . . . . . . . . . . . . . . . . . . 53 3.7 Three states used in the HLL Riemann solver. . . . . . . . . . . . 54 4.1 Loss of intensity in the beam due to the absorption. . . . . . . . 57 4.2 Angular and radial discretization of space in the gather approach. 60 4.3 Treecol method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.1 Evolution of the total molecular fraction and the molecular fraction per density, at t = 5, 10, and 15 Myr: standard MHD simulation (left), and MHD simulation without shielding (right). . . . 123 5.2 Evolution of the total molecular fraction and the molecular fraction per density, at t = 5, 10, and 15 Myr: Pure Hydro case (left), and the standard MHD simulation (right). . . . . . . . . . . . . . 124 5.3 Total column density, at t = 5, 10, and 15 Myr. . . . . . . . . . . 128 5.4 H 2 column density, at t = 5, 10, and 15 Myr. . . . . . . . . . . . 129 5.5 Total molecular fraction f (H 2 ), at t = 5, 10, and 15 Myr. . . . . 130 5.6 Slices cut through the midplane at t = 15 Myr for a pure hydrodynamical simulation (left panel) and a mhd simulation (right panel). From top to bottom: total number density, number density of H 2 , and gas temperature . . . . . . . . . . . . . . . . . . . 131 5.7 Distribution of the total shielding factor (dust and self-shielding) for the photodissociation as a function of density, at t = 5, 10, and 15 Myr. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.8 Molecular fraction as a function of total column density, at t = 5, 10, and 15 Myr. . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.9 Populations of the first 6 rotational levels calculated at thermal equilibrium, using OP R = 3, at t = 5, 10, and 15 Myr. . . . . . . 134 5.10 Total Doppler broadening parameter b in the z-axis, at t = 5, 10, and 15 Myr. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5.11 Slices cut through the midplane at t = 15 Myr for HLL solver (left panel) and HLLD solver(right panel). From top to bottom: total number density, number density of H 2 , and temperature. . . 136 5.12 Comparison of the mass in H 2 form per density bin. From left to right and from top to bottom: n ∈ (0.1, 1), (1, 10), (10, 100), (100, 1000), (1000, 10000) [cm -3 ].137 V List of Tables 1.1 Physical conditions of the different phases of the ISM . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Flow of baryons in the Milky Way, from Draine (2011).

Figure 1 . 2 :

 12 Figure 1.2: Schematic energy level structure for CII and OI, showing the transition probabilities A ul for the main cooling transitions.

Figure 1 .

 1 Figure 1.3: Standard equilibrium curve: Thermal pressure P/k vs total hydrogen density (left panel), and cooling (solid lines) and heating (dashed lines) processes in the ISM (right panel), from Wolfire et al. (1995).

Figure 1 .

 1 Figure 1.4: The interstellar cooling function for different ionization fraction x = n e /n H , from[START_REF] Dalgarno | Heating and Ionization of HI Regions[END_REF]. The discontinuity toward 10 4 K is due to the cooling by Lyman α emission, while at lower temperatures the cooling is due to CII and OI.

Figure 1 .

 1 Figure 1.5: 3D internal velocity dispersion of molecular clouds and condensations as a function of its maximum linear dimension L. The fit corresponds to σ v [km s -1 ] = 1.1L [pc]0.38 . Figure fromLarson (1981). Symbols are described in Table1in the same work, but among the most interesting regions we can recall the Taurus Molecular Cloud Complex (T), the ρ Ophiuchus Complex (ρ), and the Orion Complex (O).

Figure 1 .

 1 Figure 1.6: Eddies breaking into smaller eddies. Sketch of obstacles in a water flow (on the top) and a free water jet falling issuing from a square hole into a pool. Drawing by Leonardo da Vinci, circa 1500.

Figure 1 . 7 :

 17 Figure 1.7: Kolmogorov cascade showing the inertial range and the k -5/3 law.

Figure 1 .

 1 Figure 1.10: Polarized starlight Heiles and Crutcher (2005) and dust thermal emission (total intensity and vectors) from Planck. J-Ph. Bernard, Planck collaboration (ESLAB 2013 ).

  Frequency ARray) It is the largest existing network of radiotelescopes composed by a total of 47 stations spread in Europe (and a total of ∼ 25000 omni-directional antennas), centered in the Netherlands. It observes in the 10 -250 MHz frequency range, and it can measure the Faraday rotation and the synchrotron emission.

Figure 1 .

 1 Figure 1.11: Left panel: LOFAR stations in Europe. Right panel: Total magnetic field of the Whirpool galaxy M 51. Image from Mulcahy et al. (2014)

Figure 1 .

 1 Figure 1.12: Interstellar radiation field showing the main components. Figure from Draine (2011).

Figure 2 . 1 :

 21 Figure 2.1: Electronic states of H 2 . Potential energy curves as a function of the distance between both atoms, from Field et al. (1966).

Figure 2 . 2 :

 22 Figure 2.2: Rotations and vibrations of a diatomic molecule. Vibration corresponds to stretching modes. Rotations around the nuclear axis (x) carry little kinetic energy. Rotations around y and z axis are indistinguishable.

Figure 2 . 4 :

 24 Figure 2.4: Two main paths for the formation of molecular hydrogen on grain surface: (a) shows the Eley-Rideal mechanism, while (b) shows the Langmuir-Hinshelwood mechanism.

Figure 2

 2 Figure 2.5: Temperature dependence of the H 2 formation rate.

Figure 2 . 6 :

 26 Figure 2.6: First electronic levels of H 2 showing the Lyman and Werner transitions. Credits: Franck Le Petit

Figure 2 . 7 :

 27 Figure 2.7: Schematic energy level structure for H 2 , showing the first vibrational levels v = 0 to 3, and the rotational levels J for ortho and para-H 2 .

Figure 2

 2 Figure 2.8: Cooling function per hydrogen molecule, for n tot = 10 4 cm -3 and n(H)/n(H 2 ) = 10 -2 . Cooling functions have been calculated for two different OP R = 10 -4 and 3. From Le Bourlot et al. (1999).

Figure 3 . 1 :

 31 Figure 3.1: Comparison of Sod test for two AMR codes (Flash and Enzo), and two SPH codes (Gadget-2 and Hydra).Figure from Tasker et al. (2008)

  Figure 3.1: Comparison of Sod test for two AMR codes (Flash and Enzo), and two SPH codes (Gadget-2 and Hydra).Figure from Tasker et al. (2008)

Figure 3

 3 Figure 3.2: The top panel shows a photograph of the Antennae galaxies, one of the most outstanding example of a galaxy collision, obtained with one of the NOAO telescopes from the ground. Credit: NOAO/AURA/NSF, B. Twardy, B. Twardy, and A. Block (NOAO)). Bottom pannels shows two simulations of a galaxy merger performed with two of the most popular simulation techniques. The panel on the left: SPH simulation performed with the code Gadget-2 (image from Gadget-2 website http:// www.mpa-garching.mpg.de/gadget). The panel on the right: AMR simulation using the Ramses code (Teyssier et al. 2010).

Figure 3

 3 Figure 3.3: Refinement strategy. The figure on the left shows the tree at different levels, while the figure on the right shows the leaf cells.

Figure 3

 3 Figure 3.5: Mesh in space and time showing cells used in the Godunov scheme.

  at the boudaries between adjacent cells, each one represented by uniform states 1 U is the flux between the cell centered at x = x i and the cell centered at x = x i±1 , as shown in Fig 3.5. Time averaged flux functions are calculated by solving the Rieman problem at the interface of cells.

Two

  Alfvén waves (λ = v x ± c a ): or transverse modes, analogous to the vibration in a string, Their propagation does not change density nor pressure, Two slow magneto-acoustic waves (λ = v x ± c sma ): that rise from the coupling between the Lorentz force and the thermal pressure, when the magnetic field B and the density ρ are anticorrelated Two fast magneto-acoustic waves (λ = v x ± c fma ): also rise from the coupling between Lorentz force and thermal pressure, but B and ρ are correlated, One entropy wave (λ = v x ): this is a contact discontinuity, and it changes only density, where c a is the Alfvén speed, c sma and c fma are respectively the slow and the fast magneto-acoustic speeds. Both magneto-acoustic (alson known as magnetosonic waves), are compresional modes, and propagate changing density and pressure. The Alfvén speed and slow and fast magneto-acoustic waves are

Figure 3

 3 Figure 3.6: Wave fan for MHD waves.

Figure 3 . 7 :

 37 Figure 3.7: Three states used in the HLL Riemann solver.

Figure 4 . 1 :

 41 Figure 4.1: Loss of intensity in the beam due to the absorption.

Figure 4 . 2 :

 42 Figure 4.2: Angular and radial discretization of space in the gather approach.

Figure 4 . 3 :

 43 Figure 4.3: Treecol method.

Fig. 1 .

 1 Fig. 1. Example of the discretization for N = 4 and M = 3. The red lines show the bounds of the solid angles, while the blue arrows show the representative directions (θ, φ) for the column densities.

Fig. 2 .

 2 Fig.2. Central region containing the target cell and the cells that belong to its oct (sibling cells) at the resolution level of the target cell 0 . The subsequent domains are constructed around the grid decreasing the resolution consecutively. In this manner the i-th shell will contain cells at level 0 -i.

Fig. 3 .

 3 Fig. 3. Standard cube showing all the possible configurations for the two-level approximation. The size of each cell in the standard cube is considered to be one unit in terms of the local dx.

Fig. 4 .

 4 Fig. 4. Effect of the two level approximation on the estimate of single column densities. The parent grid that contains the target cell is marked in red.

Fig. 5 .

 5 Fig.5. Column density maps integrated along the x, y and z axis for the reference (on the top) and for the tree-based method (center) calculated with respect to the mid planes. The map at the bottom corresponds to the fractional error calculated according to Eq. (12).

Fig. 6 .

 6 Fig.6. Column densities for the 4π sr as seen by a cell sitting on the edge of the sphere(12.5, 25, 25) [pc]). The abscissa corresponds to the angle φ (from 0 to 2π), and the cos θ (from -1 to 1) is given in the ordinate. The reference maps (on the left) have been calculated using the gather approach. The panels at the center present the maps calculated using the tree-based method, while on the right we present the error maps calculated according to Eq. (12). From top to bottom: using 12, 40, 84, 144, and 220 directions.

Fig. 8 .

 8 Fig.8. Difference maps for the extinction as seen by cells in the midplanes for x, y, and z. From top to bottom: using the tree-based method for 6, 12, 40, and 84 directions.

Fig. 9 .

 9 Fig.9. Column density maps for a highly structured cloud, integrated along the x, y, and z axes (from left to right): using the ray-tracing method (on the top), our tree-based method (middle), and the relative error (on the bottom).

Table 2 .

 2 Mean error for the mean column density maps.

Fig. 11 .

 11 Fig. 11. Extinction maps as seen by the cells in the midplanes. On the top we show the reference, calculated with the ray-tracing method. The other maps (from top to bottom) have been calculated with the tree-based method using 6, 12, 40, and 84 directions.

  Fig. 12. Relative difference maps for the tree-based method. From top to bottom: 6, 12, 40, and 84 directions.

Table 3 .

 3 Calculation time relative to the case without screening and the average difference of extinction maps.Directionst sim /t 0 |χ precχ in situ | M of the tree-based method. The first implementation calculates each contribution to the column density for each cell without doing any approximation. This implementation is called in situ. The second implementation uses the precalculation of geometrical terms

Fig. 15 .

 15 Fig. 15. Influence of the precalculation module. Extinction maps calculated using 84 directions (M = 7, N = 12). At the top we present the extinction map for the case where all the calculations are done in situ, and in the middle using the precalculation module. At the bottom we present difference maps.

Fig. 16 .

 16 Fig. 16. Slices cut through the midplane at t = 10 Myr showing the detail of difference in fragmentation. The left panel shows the case without screening, and the right panel shows the case with screening. The fragments zoomed in have a size of 5 pc and are centered at x = 15, y = 2.75 pc with respect to the center of the box.

Fig. 17 .

 17 Fig. 17. Comparison of the probability distribution function (PDF) of gas (left panel) and the corresponding temperature per density bin (right panel) at t = 5 Myr (top) and at t = 10 Myr (bottom) for the cases with and without screening.

Fig. 18 .

 18 Fig. 18. Comparison of the mass spectra. From top to bottom at t = 5, 7.5, 10 Myr. From left to right for the density thresholds: n th = 500, 1000, 2500, 5000 cm -3 .

Fig. A. 1 .

 1 Fig. A.1. Comparison of the PDF (top) and the respective temperature per density bin (bottom) at t = 10 Myr, obtained including the screening calculated with the tree-based method using 6, 12, 40, and 84 directions.

Fig. A. 2 .

 2 Fig. A.2. Comparison of the mass spectrum of clumps at t = 10 Myr for different estimates of the column density at low resolution. From top to bottom: n th = 500, 1000, 2500, 5000 cm -3 .

Fig. 1 .

 1 Fig. 1. Column density maps. Total hydrogen column density map (left), H 2 column density map (center), and molecular fraction (right). The gas enters the box from the left and right sides (x-axis), and the conditions for the other faces are periodic. The gas is compressed along the x direction and the cloud forms toward the middle of the box.

Fig. 2 .

 2 Fig. 2. Slices cut through the mid-plane. The panel on the left shows the total number density of hydrogen nucleons, the panel on the center shows the number density of molecular hydrogen (on the same scale), and the panel on the right shows the mean shielding factor for the H 2 photodissociation rate coefficient, calculated as exp(-τ d,1000 ) f shield .

Fig. 3 .

 3 Fig. 3. Density PDF evolution. Color lines show the density PDF at t = 5 (blue), 10 (green), 15 (red), and 20 Myr (light blue).

Fig. 4 .

 4 Fig. 4. Volume filling factor as a function of density in different regions. V tot is the volume fraction of each region.

Fig. 5 .

 5 Fig. 5. The top panel shows the thermal pressure, P th = knT , the middle panel shows the ram pressure, P ram = ρV 2 , while the bottom panel shows the magnetic pressure, P mag = B 2 /8π in the simulation box. The red line shows the mean value per density bin. The color scale indicate the number of points.

Fig. 6 .

 6 Fig. 6. Statistics of clumps. The top panel shows the mass spectrum, the second panel their internal velocity dispersion as a function of size, the third panel the kinetic β = P ram /P mag and the bottom panel displays the H 2 fraction as a function of mass.

Fig. 7 .

 7 Fig. 7. UV screening factor vs density (top panel) and total shielding coefficient for H 2 vs density (bottom panel) at t = 20 Myr. The blue line shows the mean value per density bin, and the color scale indicates the number of points.

Fig. 9 .

 9 Fig. 9. Molecular fraction evolution. On the top: per density bin (from purple to red the density increases), on the bottom: per temperature bin (the temperature increases from purple to salmon).

Fig. 10 .

 10 Fig. 10. Mass in H 2 form per density bin. Comparison between the simulation (blue) and the calculation at equilibrium (gray) for different timesteps. From left to right: t = 5, 10, 15, and 20 Myr. From top to bottom: density bins n ∈ (0.1, 1), (1, 10), (10, 100), (100, 1000), (1000, 10000) [cm -3 ]. Note the change of vertical scale.

Fig. 11 .

 11 Fig. 11. Evolution of the molecular fraction, where the formation of H 2 has been suppresed for gas denser than a fixed threshold. The top panel corresponds to the standard case, with no suppresion of H 2 formation. The middle panel shows the evolution for a density threshold of n = 1000 [cm -3 ], while the panel on the bottom one corresponds to a density threshold of n = 100 [cm -3 ]

Fig. 12 .

 12 Fig.12. Various cooling and heating contributions as a function of density. As can be seen H 2 does not have a dominant influence except at densities of about 4 cm -3 where H 2 cooling (green line) becomes comparable to the standard ISM cooling (blue line) and at high densities where the heating by H 2 formation (red line) is significant (but without modifying substantially the temperature).

Fig. 13 .

 13 Fig. 13. Molecular fraction as a function of total hydrogen column density. Comparison with Copernicus survey(Savage et al. 1977), and FUSE(Gillmon et al. 2006;Rachford et al. 2002Rachford et al. , 2009)).

Fig. 14 .

 14 Fig. 14. H 2 population distribution along different lines of sight, using an ortho-to-para ratio of 3 (top), and an ortho-to-para ratio at thermal equilibrium (bottom). The symbols correspond to the observational data of Wakker (2006) (cross), Rachford et al. (2002) (circle), Gry et al. (2002) (square) and Lacour et al. (2005) (diamond). Colors correspond to the total column density range.

Fig. 15 .

 15 Fig. 15. H 2 level population distribution for individual cells in the simulation calculated using an ortho-to-para ratio at thermal equilibrium.

Fig. 16 .

 16 Fig. 16. Mean velocity dispersion along a line of sight for excited levels of H 2 and comparison with the data from Lacour et al. (2005). Top panel is for lines of sight along the x-axis, middle one is along y-axis and bottom panel along z-axis.

Fig. A. 1 .

 1 Fig. A.1.Density PDF for a series of numerical resolutions. While the discrepancy between between low and high resolution runs is very significant. The difference between the resolution of the standard and high resolution runs is not very high except for the highest density bins, which contain little mass.

Figure 5 . 1 :

 51 Figure 5.1: Evolution of the total molecular fraction and the molecular fraction per density, at t = 5, 10, and 15 Myr: standard MHD simulation (left), and MHD simulation without shielding (right).

  Figure 5.2: Evolution of the total molecular fraction and the molecular fraction per density, at t = 5, 10, and 15 Myr: Pure Hydro case (left), and the standard MHD simulation (right).

  Figure 5.3: Total column density, at t = 5, 10, and 15 Myr.

Figure 5 . 6 :

 56 Figure 5.6: Slices cut through the midplane at t = 15 Myr for a pure hydrodynamical simulation (left panel) and a mhd simulation (right panel). From top to bottom: total number density, number density of H 2 , and gas temperature

Figure 5 . 8 :

 58 Figure 5.8: Molecular fraction as a function of total column density, at t = 5, 10, and 15 Myr.

Figure 5 . 9 :

 59 Figure 5.9: Populations of the first 6 rotational levels calculated at thermal equilibrium, using OP R = 3, at t = 5, 10, and 15 Myr.

Figure 5 .

 5 Figure 5.10: Total Doppler broadening parameter b in the z-axis, at t = 5, 10, and 15 Myr.

Figure 5 .

 5 Figure 5.11: Slices cut through the midplane at t = 15 Myr for HLL solver (left panel) and HLLD solver(right panel). From top to bottom: total number density, number density of H 2 , and temperature.

Figure

  Figure 5.12: Comparison of the mass in H 2 form per density bin. From left to right and from top to bottom: n ∈ (0.1, 1), (1, 10), (10, 100), (100, 1000), (1000, 10000) [cm -3 ].

  . . . . . . . . . . . . . . . . . . 2.1 Electronic states of H 2 . Potential energy curves as a function of the distance between both atoms, from Field et al. (1966). . . . . 2.2 Rotations and vibrations of a diatomic molecule. Vibration corresponds to stretching modes. Rotations around the nuclear axis (x) carry little kinetic energy. Rotations around y and z axis are indistinguishable. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Ortho-to-para ratio as a function of temperature for thermal equilibrium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Two main paths for the formation of molecular hydrogen on grain surface: (a) shows the Eley-Rideal mechanism, while (b) shows the Langmuir-Hinshelwood mechanism. . . . . . . . . . . . . . . 2.5 Temperature dependence of the H 2 formation rate. . . . . . . . . 2.6 First electronic levels of H 2 showing the Lyman and Werner transitions. Credits: Franck Le Petit . . . . . . . . . . . . . . . . . . 2.7 Schematic energy level structure for H 2 , showing the first vibra-

Table 1 .

 1 1: Physical conditions of the different phases of the ISM

	Phase	n H (cm -3 )	T (K)
	Dense H 2	10 3 -10 6	10 -50
	Diffuse H 2	∼ 100	∼ 50
	CNM (cool HI)	30	∼ 100
	WNM (warm HI)	0.6	∼ 5000
	HII gas	0.2 -10 4	10 4
	HIM (coronal gas)	∼ 0.004	10 5.5

  51 rovibrational levels in total). The cooling function per H 2 molecule is, thus, given by

Table 1 .

 1 Mean difference relative to the reference maps for the extinction.

	Directions M × N	∆χ spherical cloud turbulent cloud mean max mean max
	6 3 × 4 = 12 5 × 8 = 40 7 × 12 = 84 0.018 0.090 0.036 0.203 0.017 0.094 0.014 0.068	0.079 0.280 0.054 0.204 0.048 0.164 0.047 0.145
	maps calculated according to Eq. (

  Article number, page 19 of 31page.31In particular the equilibrium distribution is expected to be different from the time-dependent one for at least two main reasons. First, since the H 2 formation time-scale is somewhat long, if the time-dependent f (H 2 ) lies below the equilibrium one, it will indicate that what has limited the H 2 abundance is the time to form the molecule. Second, if the time-dependent f (H 2 ) is above the equilibrium value, it will indicate that f (H 2 ) has increased because of an enrichment coming from denser gas. As discussed above this could operate through the expansion or the evaporation of cold clumps, a process that we call turbulent diffusion, or through numerical diffusion. The latter is quantified in the appendix by performing convergence studies which suggest that it remains limited. Another possibility is that the H 2 excess has been produced during a phase where the fluid particle was more shielded by the surrounding material and therefore the UV field was lower. The two effects are probably acting simultaneously and are difficult to disentangle.At time 5 Myr (left panels), all density bins show an excess of the equilibrium distribution with respect to the time-dependent one. This is clearly due to the long time-scale needed to form H 2 .

Hydrogen is ionized with 13.6 eV, photodissociation of H

begins at 14.7 eV, and H 2 photoionization of H 2 begins at 15.4 eV

When H 2 is destroyed by collisions, it contributes to cool the gas

High order schemes uses a distribution around the mean value. For example, in the second order scheme the states of each cell are described by the mean value and a slope.

This approximation is one of the most used due to its simplicity(Heitsch et al. 2005;Micic et al. 2013; Vázquez-Semadeni et al. 2007;Banerjee et al. 2009; Audit and Hennebelle 2005), but is valid only for very low density gas, with very low column densities in all directions.

if the incoming flow is fast enough a collision is not even necessary as it is the case for example in the study presented byKoyama & Inutsuka (2002)

Article number, page 31 of 31page.31

in post-processing have been provided by Maryvonne Gerin.

Acknowledgements

Acknowledgements. V.V. acknowledges support from a CNRS-CONICYT scholarship. This research has been partially funded by CONICYT and CNRS, according to the December 11, 2007 agreement. P.H. acknowledges the finantial support of the Agence National pour la Recherche through the COSMIS project. This research has received funding from the European Research Council under the European Community's Seventh Framework Program (FP7/2007-2013 Grant Agreement No. 306483). Acknowledgements. V.V. acknowledges support from a CNRS-CONICYT scholarship. This research has been partially funded by CONICYT and CNRS, according to the December 11, 2007 agreement. P.H. acknowledges the finantial support of the Agence National pour la Recherche through the COSMIS project. This research has received funding from the European Research Council under the European Community's Seventh Framework Program (FP7/2007-2013 Grant Agreement No. 306483). M.G. and P.L. thank the French Program "Physique Chimie du Milieu Interstellaire" (PCMI) This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche

List of Objects Appendix A: Low resolution study

In order to ensure that the difference on the statistical properties depends on the inclusion of the screening but not on the method used for estimating it, we present a low resolution study of the influence of the method. We have performed four identical simulations, as decribed in section 5, where the screening has been calculated using 6, 12, 40, and 84 directions. We allow two AMR levels ( min = 7 and max = 9), with an effective numerical resolution of 512 3 cells and spatial

List of Objects Appendix A: The issue of numerical convergence: resolution study

Numerical resolution is a crucial issue particularly in the context of chemistry. Here we present a comparison between runs of various resolutions. and 2048 3 (for the two last and higher resolutions this represents the highest effective resolution since AMR is used as described above). The density PDF of low resolution runs is quite different from the higher ones showing the importance of numerical resolution. The two highest resolution runs present similar though not identical PDF. This suggests that for the highest resolution runs, numerical convergence for the density PDF is nearly reached although strictly speaking, it would be necessary to perform even higher resolution runs. This conclusion agrees well with the results of Federrath (2013) who find that for isothermal self-gravitating simulations, the density PDF converges for a resolution inbetween 1024 3 and 2048 3 .