Dr Rosine Lallement

Dr Chiara Ferrari

Dr Cécile Gry

Dr Pierre Hily-Blant

Dr Marc-Antoine Miville-Deschenes

Michel Perault

Grazie Ad Alessia

E Bellomi
email: elena.bellomi@obspm.fr

B Godard

P Hennebelle

V Valdivia

G Pineau Des Forêts

P Lesaffre

M Pérault

3D chemical structure of diffuse turbulent ISM I. Statistics of the HI-to-H 2 transition

Keywords: che sono stati fantastici come sempre, lontani o vicini sempre presenti, sempre positivi ISM: structure -ISM: molecules -ISM: kinematics and dynamics -ISM: clouds -methods: numericalmethods: statistical

for taking the time to read my manuscript and for the interesting discussion we had around it. I am deeply grateful to all the people that contributed to improving this manuscript with their corrections and words of advice. In particular, my Ph

IV Conclusion and

We look up at the night sky and, if we are lucky enough to be far from the cities' light pollution, we can observe a myriad of bright dots lighting up the sky's darkness. Those white, blue, red luminous points can be stars in our Galaxy or even entire galaxies, so far away from us that they appear like dots, or again, as the Greeks observed, if the position of these dots moved from one night to another they can be planets of our solar system (Planetes in ancient Greek meaning wanderer). Stars have played a fundamental part in human history. For instance, they helped the first explorers to navigate the seas and oceans, they ignited passionate debates between philosophers throughout the ages, and always made humankind dream. After all, hasn't the first scientific experiment been to observe the sky? Allowing us to grasp this mystery that contains us and inhabits us. As Carl Sagan said, "We are made of star-stuff. We are a way for the Cosmos to know itself". Indeed, in their life cycles, stars produce heavy nuclei in their cores thanks to multiple fusion processes and during supernovae explosions. These cataclysmic cosmic events are the "death" of certain stars. Last but not least, they produce light! Stars usually live in galaxies, where they are not quiet guests, they continuously change their surrounding environment, with their continuous feedback processes, radiation, and winds.

Contrary to human conception, the space separating the stars and filling the Galaxy is diluted but far from empty. This intermediate space called the interstellar medium (ISM) is composed of gas and dust grains (small solid particles mixed with the interstellar gas). The ISM has an average density of about one particle per cubic centimeter. The ISM is a very rich, turbulent, and magnetized medium. Gas and dust are in constant interactions with magnetic fields, cosmic rays, i.e. high energy particles, gravitational perturbations, turbulence, and electromagnetic radiations (photons): together, they govern the dynamics of the ISM. All these constituents have a great impact on the dynamic of the gas. For instance, magnetic pressure works against gravitational collapse in molecular clouds, while magnetic tension helps to confine filaments of the ISM increasing their lifetimes. Turbulence mixes the gas, transporting the energy to smaller and smaller scales, hence, connecting the different physical scales. Photons interact with matter thus heating the medium and changing its chemical content. The different constituents listed above are in almost equipartition of energy density, meaning that, generally speaking, one does not prevail on the others, and none is negligible compared to the others. Moreover, these constituents are all coupled with one another, influenced by each other, making the ISM a complex and highly dynamic medium.

The ISM is the birthplace of stars, their nursery, their home, and eventually their death bed and grave. Molecular clouds, formed from diffuse atomic gas, become increasingly denser subsequently condensing into sheets and filament structures and cores that in time collapse, under the effect of their own gravitational pull, into stars when gravity overcomes thermal and magnetic pressure. Stars will burn their fuel, usually hydrogen and helium at first through thermonuclear fusion. This nuclear process will synthesize heavier elements, and modify the surrounding environment through the injection of energy, momentum, and mass in the ISM. All these injections are called feedback processes, and examples are stellar winds, the ionizing radiation from massive stars, the HII regions surrounding such Figure 1: Picture taken from https://chandra.harvard.edu/ "NGC 6357, a region where radiation from hot, young stars is energizing the cooler gas in the cloud that surrounds them. This composite image contains X-ray data from NASA's Chandra X-ray Observatory and the ROSAT telescope (purple), infrared data from NASA's Spitzer Space Telescope (orange), and optical data from the SuperCosmos Sky Survey (blue) made by the United Kingdom Infrared Telescope. Located in our galaxy about 5,500 light years from Earth, NGC 6357 is actually a "cluster of clusters," containing at least three clusters of young stars, including many hot, massive, luminous stars. The X-rays from Chandra and ROSAT reveal hundreds of point sources, which are the young stars in NGC 6357, as well as diffuse X-ray emission from hot gas. There are bubbles, or cavities, that have been created by radiation and material blowing away from the surfaces of massive stars, plus supernova explosions." Image credit: X-ray: NASA/CXC/PSU/L. Townsley et al; Optical: UKIRT; Infrared: NASA/JPL-Caltech Figure 2: This image is taken with the NASA/ESA Hubble Space Telescope. "When a massive new star starts to shine while still within the cool molecular gas cloud from which it formed, its energetic radiation can ionize the cloud's hydrogen and create a large, hot bubble of ionized gas. Amazingly, located within this bubble of hot gas around a nearby massive star are dark compact globules of dust and gas, some of which are giving birth to low-mass stars." stars, and UV radiation. The end of life of a star will differ depending on its mass. When massive enough they will burst in a supernova explosion thus injecting a large quantity of mechanical energy and heavy nuclei into the surrounding medium.

The interstellar medium interacts greatly with the stars, which are photon emitters. As mentioned before, photons and matter interact in a complicated dance. It is mainly thanks to photon-matter interactions that we get information on the medium in which they travel. Photons are the prime messengers of the Universe surrounding us. Everything we see is light or the absence of light. Our understanding of the ISM, the Galaxy, and even the Cosmos essentially come from photon measurements. However, photons are not the only messengers of our Universe. Neutrinos and more recently gravitational waves take now part in our new era of multi-messenger astrophysics. Indeed neutrinos which are fundamental particles that only interact through weak interaction are the second most abundant particles in the Universe and can, thanks to their properties, give us information on the very early Universe even before recombination times when the Universe stopped being opaque to photons. The same goes for gravitational waves which propagate at the speed of light in the space-time fabric and give us early warning of huge astrophysical phenom-ena such as black-hole merging. It is important to notice that since almost all information comes from one messenger, and is two-dimensional, the interpretation of observations can be challenging. Many quantities observed are integrated quantities along the line of sight, e.g. column densities, and it is sometimes difficult to link the 3D physics with the 2D data coming from the observations. A natural way to study the ISM on its whole, is to access to its properties though the observation of the local interstellar medium. Ongoing efforts to observe the local diffuse ISM unveiled a richness of features that characterize the space around us. Observations show a complex structure of matter characterized by dense, and cold, and almost fractal regions surrounded by hot and diffuse gas. These phases, defined by very different densities and temperatures, exchange matter and energy, and coexist because they are approximately in thermal pressure equilibrium. A detailed observation of the ISM, through spectroscopy also revealed a richness in molecules. Molecules not only exist in a medium that is very hostile to their survival because rich in UV photons, but they are also abundant enough to be detected. As we go around the cycle of matter, more and more complex chemical species are found, showing the importance of chemistry as an observable tracer of the evolution of the ISM. All these characteristics seem to be global properties of the local ISM because they appear to be observed in all directions.

These observations of the diffuse local ISM reveal a complex medium and require a strong theoretical counterpart to be explained and interpreted. Usually, the theoretical efforts are divided into two categories: on one side numerical simulations to explain the formation of the multi-scale structures; and on the other chemical models to explain the degree of ionization and the abundances of the different chemical species.

Ordinarily, the chemical solvers are one-dimensional models that solve a large and complex chemical network at steady-state, because it is virtually impossible to compute 3D time-dependent chemistry. An example is PDR codes (photons-dominated regions) which model regions in which the chemistry is driven by the UV photons. These codes usually explain well the abundances of many chemical species, however fail to explain the abundances of some species which are formed through highly endoenergetic reactions. Moreover, there is no coupling between the chemistry and the dynamics, therefore the structure of the medium and its subsequent permeability to photons are not considered. Other models, such as TDR (turbulent dominated regions) rose as attempts to couple chemistry and dynamics at small scales and explain the abundances of species with endothermic formation reactions. TDR codes give good results on the abundances of species whose formation is endothermic, nevertheless, the dynamics and the chemistry are not coupled at large scales, and the distribution of the positions and the velocities of the turbulent dissipation are introduced as ad hoc parameters.

To describe the evolution of physical quantities like the density, the temperature, the velocity, and the magnetic field, extensive numerical simulations are required. On one side, the magnetohydrodynamical equations describing the dynamics of the ISM are highly non-linear and therefore no analytical solution exists. On the other side, the physical con-Figure 3: Picture taken from https://science.nasa.gov/. "The Great Nebula in Orion, an immense, nearby starbirth region, is probably the most famous of all astronomical nebulas. Here, glowing gas surrounds hot young stars at the edge of an immense interstellar molecular cloud only 1500 light-years away. In the featured deep image in assigned colors highlighted by emission in oxygen and hydrogen, wisps and sheets of dust and gas are particularly evident. The Great Nebula in Orion can be found with the unaided eye near the easily identifiable belt of three stars in the popular constellation Orion. In addition to housing a bright open cluster of stars known as the Trapezium, the Orion Nebula contains many stellar nurseries. These nurseries contain much hydrogen gas, hot young stars, and stellar jets spewing material at high speeds. Also known as M42, the Orion Nebula spans about 40 light years and is located in the same spiral arm of our Galaxy as the Sun." Image Credit & Copyright: Josep M. Drudis & Don Goldman ditions of density and temperature are too extreme (e.g. n H ∼ 0.01 cm -3 , T ∼ 10 5 K) to be able to conduct laboratory experiments. Even with powerful calculators, the challenge remains difficult, because on the one hand simulations need to cover a huge range of scales (from stars formation to the galaxy rotation), and on the other, they need to treat several competing and entangled physical processes.

These two facetschemistry and dynamical simulationswere treated in parallel and independently from each other until recently. Recent theoretical efforts have started to couple the two approaches by implementing the time-dependent evolution of a simplified chemical network in the simulations and by recomputing the abundances of other chemical species at each point in post-treatment. This leads to a more dynamic and coherent treatment of the problem. In particular, it is important to include the processes of formation and destruction of molecular hydrogen in the MHD simulations of ISM mainly for three reasons. Firstly, H 2 is the most abundant molecule in the ISM, secondly, it has long timescales compared to other species, i.e. it is more likely to be out-of-equilibrium, and finally, it is the precursor to the rest of the chemistry.

In this Ph.D. thesis, we continue the effort of combining chemistry and dynamics in a more coherent picture, where at least part of the chemical and dynamical evolutions are coupled. Our goal is to understand the physical processes governing the evolution of the local diffuse ISM and the effect of the multiphase structure on its 3D chemical composition. Numerical simulations allow to conduct a series of experiments for different initial conditions and parameters. Even though these numerical experiments provide important information on their own, the real power resides in their connection with observations. Indeed, the combinations of observations and simulations is a powerful tool to constrain the physics of the simulations, to highlight the influence of different parameters, and to interpret the observational data. In this context, the rising amount of data collected by spectrometers from radio to ultraviolet (UV) wavelengths opens a new era where the statistical and chemical information contained in the observations can be used concomitantly. The originality of this work stands in developing new tools and defining a framework to effectively compare observations and simulations, taking advantage of the union of multiphase simulations and a chemical solver on one side, and the increasing number of observations on the other side. The tools developed in this thesis are a first step toward an N-dimensional investigation in which, multiple observables would be studied simultaneously maximizing the exploitation of the statistical power contained in the observations. Structure of the Ph.D. Manuscript In the first chapters, i.e. Part I, we present the diffuse interstellar medium in general, its components ( §1), the thermodynamics of the multiphase diffuse gas ( §2), the interaction of matter and photons ( §3), and the three observables we concentrate on in this work ( §4): the hydrogen and its transition from its atomic to molecular form ( §4.1), the neutral carbon ( §4.2) and the elusive CH + ( §4.3). In Part II, we present different aspects of the RAMSES simulations ( §5), the turbulent forcing, the out-of-equilibrium molecular hydrogen, and the computation of the excitation of neutral carbon ( §6) and molecular hydrogen ( §7). In Part III, we present the tools to compare ob-servations and simulations: the reconstruction of lines of sight ( §8), which allows to build the observed distribution of lines of sight from the simulations, and our modifications of the , including some tests and application. In the last chapters of Part III, we present the results obtained on the HI-to-H 2 transition ( §10), CH + ( §11), and CI excitation ( §12).

Chapter 1

ISM Composition

The interstellar medium is a dilute mixture of gas, high energy particles, and dust grains. The ISM consists mainly of hydrogen which carries around 70% of the mass, followed by helium which accounts for 28% of the mass. The remaining 2% is carried by heavier elements that are commonly called metals (we will adopt this terminology).

Since the medium is rich in atomic hydrogen, which has an energy potential of 13.6 eV, UV photons with hν > 13.6 eV are effectively absorbed by the interstellar matter. This leads to the existence of a neutral phase in which species with an ionization potential smaller than 13.6 eV are in their singly ionized form while species with a potential larger than that of hydrogen are primarily neutral.

The state and evolution of the diffuse neutral ISM are at the root of the present thesis. In this chapter, we make a brief presentation of its composition, its structure related to the other phases of the ISM, and the energy sources that need to be considered.

Dust

Even though dust grain carry only 1% of the mass of the ISM, they play an important role in the physics and chemistry of the medium by locking elements in their mantles and cores, by affecting the thermodynamics through heating and cooling, by absorbing and scattering the UV photons coming from the stars and acting as catalysts for chemical reactions.

The existence of dust grains was first indicated by the dimming of starlight. In the 30s Trumpler tried to determine how much light was attenuated by the "absorbing medium" between us and the stars. Using a pair method (stars with same spectral type), the extinction due to the interstellar medium as a function of the wavelength could be traced. This curve is called the "extinction curve" and presents a bump at 217.5 nm that gives away the presence of dust grains: since this feature is too broad to be associated with atoms and small molecules, it traces the existence of large carboneus-rich complexes. Moreover, since different wavelengths are absorbed and scattered differently, and dust grains can only absorb photons with wavelengths smaller than their size, the observation of the extinction curve provides constraints on the size of grains. Mathis et al. (1977) built a model of a grain size distribution that could reproduce well the extinction curve, assuming graphite and silicate spheres. This study led to the wellknown MRN model in which: n(a) ∝ a -3.5

(1.1)

where n(a) is the density of grain of size a per unit size between ad a/2 and a + d a/2. This distribution has the interesting property that the total surface of grains

a max a mi n a 2 n(a)d a ∝ a -1/2 mi n
is mostly given by the small grains, while their total mass a max a mi n a 3 n(a)d a ∝ a1/2 max is controlled by the large grains.

Nowadays the theory of dust leans on many other observational pieces of evidence, such as the abundance of elements in the gas phase compared to the solar abundances, the polarization of light, the infrared emission, and the abundances of several chemical species, like molecular hydrogen, which requires grains to be formed. Among those studies, polarization is particularly interesting. At the end of the 40s, a correlation between the polarization of the light and the reddening was pointed out: the stronger the reddening the more the light is polarized; in addition, stars close to one another show similar polarizations. These observations led to the conclusion that dust polarizes the light coming from stars. This is presently understood as a consequence of an alignment between dust and the magnetic field, showing that dust grains are not spherical and are coupled dynamically to the magnetic field lines.

Phases

Observations of the ISM performed at different wavelengths reveal a medium with different thermal and chemical states, characterized by very distinctive physical conditions. For instance, X-ray and UV reveal the presence of gas at temperatures larger than 10 5 K, while observations at millimeter wavelengths highlight the presence of Giant Molecular Clouds at temperatures below 100 K. All these observations show that the ISM is highly inhomogeneous and can be classified in phases with different temperatures, densities, and ionization fractions (Lequeux 2005;Draine 2011).

HIM The Hot Ionized Medium, also called coronal gas is an ionized phase with a typical temperature T ∼ 10 5 -10 6 K, and low density n H ∼ 10 -2 cm -3 . This very diffuse and hot medium is heated and ionized by supernovae explosions and stellar winds and can take up to ∼ 50% of the volume of the disk in the galaxy. It is observed through UV absorptions by highly ionized species (e.g. CIV, SIV, NV, OVI), continuum emission generated by free electrons accelerated in the Coulomb field of positive ions or the radiative recombination of hydrogen and helium ions with free electrons, soft X-rays 1 , and spectral emission lines (OVII, OVIII) excited by collisions followed by radiative de-excitation.

WIM The Warm Ionized Medium is an ionized phase, with characteristic temperature and density of T ∼ 8000 K, and n H ∼ (0.2 -) 0.5 cm -3 . Its important ionization fraction is mainly due to photons (stellar radiation) and partially to shocks and collisions (ejection of matter from HII regions). This medium can take up to ∼ 10 -25% of the volume of the ISM and takes up ∼ 90% of the ionized hydrogen (the rest is in HII regions). This medium is traced by the recombination of forbidden and optical line emissions of atoms and ions: Hα (at 6563 Å), SII, NaII, OI, OIII, HeI (Sembach et al. 2000), and by ionic absorption in the visible and the ultraviolet.

WNM

The Warm Neutral Medium is the atomic neutral phase characterized by the same temperature as the WIM T ∼ 8000 K, and a density n H ∼ 0.5 (-0.6) cm -3 . The WNM is believes to occupy ∼ 30 -40% of the volume of the Galactic disk and to carry a significant fraction of the total mass of the neutral hydrogen. Its main observational tracer is the emission line of HI at 21 cm.

CNM

The Cold Neutral Medium is the atomic neutral phase characterized by T ∼ 70 K (30-100K) and n H ∼ 50 -100 cm -3 . It occupies a very small percentage of the volume of the galactic disk, only a few percent, but since it is quite dense, it takes up 30% -60% of the mass of the ISM. The CNM is observable through absorption in UV and visible, and line emission of C + which is an important coolant for this phase. Dense molecular clouds are colder (T ∼ 10-20 K) and denser (n H ∼ 10 3 -10 6 cm -3 ) environment within the CNM bounded by self-gravitation. This fully molecular medium is the residence of complex molecules, the playground of dust chemistry, and the site where star formation occurs. Their main tracers are molecular absorption and emission lines observed in the millimeter domain.

The phases described above are not a universal classification and can be referred to with different names. Several authors prefer to use different criteria, such as the extinction or the chemical state of hydrogen. For instance, Snow & McCall (2006) classify different CNM and molecular clouds based on the value of the visual extinction A V . Their classification is more refined and leads to the definition of the diffuse atomic, diffuse, translucent and dense molecular cloud.

Reservoirs of Energy

Non-isolated, the ISM is constantly pervaded by several energy sources. A major ingredient is the stars themselves. Stars not only irradiate the ISM with photons at different wavelengths but regularly inject mass and mechanical energy in the medium through stellar winds and supernovae explosions. A second important energy source is the cosmic rays. These relativistic particles accelerated in supernovae explosions propagate through the ISM with energy sufficient to heat and ionize the gas. At last, the galactic dynamics itself and the sustained large-scale magnetic field provide another important energy input. The magnetic field not only provides dynamical support to the gas through magnetic tension and magnetic pressure but also affects the propagation of cosmic rays.

Mechanical Energy and Turbulence

The ISM is a very dynamic medium. This dynamical state is maintained by the gravitational well of the Galaxy, its differential rotation, the stellar winds, the supernovae explosions of massive stars, and the protostellar outflows which constantly inject massive amounts of mechanical energy in the gas. As a result, the line profiles observed through spectroscopy exhibit supra-thermal broadenings, indicating that the medium is highly turbulent. Indeed, the Reynold number estimated from the velocity dispersion at a scale of 10 pc is of the order of 10 7 : the mechanical energy injected at this scale or above, necessarily induces a turbulent cascade that transport the kinetic energy to the viscous scales where it is dissipated. Because of this long dissipation process, mechanical energy is a major energy reservoir of the ISM. Quantitavely, the density of mechanical energy, or ram pressure, can be written as function of the mass density ρ and the speed v, as

P dyn = ρv 2 ∼
n H 10 cm -3 v 3.5 kms -1 2 × 1.7 × 10 -12 erg cm -3 .

(1.2)

Magnetic Field

Magnetic fields are ubiquitous in the ISM, in the Galaxy, and in the Universe. Even though the origin and the evolution of the magnetic field are still open questions, its importance on the ISM is well established. Because the ISM is at least partly ionized, the magnetic field interacts with gas and dust through magnetic tension and magnetic pressure. This not only affects the entire turbulent cascade and the transport of mechanical energy between scales (Elmegreen & Scalo 2004), but also the structures and evolutions of the different phases of the ISM. Moreover, the magnetic field participates in the regulation of star formation: not only magnetic pressure opposes gravity but also the structure of the magnetic field is key for the evacuation of prestellar angular momentum during the phase of collapse.

Several indirect methods exist to measure the different components of the magnetic fields (parallel or perpendicular to the line of sight): the Zeeman splitting, the synchrotron emission, the Faraday rotation, and the dust polarization. With a typical magnetic field B ∼ 5 µm the magnetic pressure in the local ISM writes as

P mag = B 2 8π ∼ B 5 µm 2 × 10 -12 erg cm -3 (1.3)
and is therefore an important energy reservoir of the ISM.

Radiation Field

The photons constitute another important energy reservoir. The photons pervading the ISM originate from several sources including the Cosmic Microwave Background, the stars, and the dust emission, which all contribute to the radiation energy density at different wavelengths. Separately, these sources amount to typical energy densities of u stars ∼ 10 -12 erg cm -3 , u CMB ∼ 4.2 × 10 -13 erg cm -3 , and u IR,dust ∼ 5 × 10 -13 erg cm -3 (see Table 12.1 of Draine 2011). The corresponding total energy density contained in photons therefore sums up to P rad ∼ 2.2 × 10 -12 erg cm -3 .

(1.4) which is comparable to the magnetic and mechanical energy densities. Among all photons, the UV radiation field is of particular importance. Indeed, the interaction of the UV photons with interstellar matter plays a key role in the thermodynamics and chemistry of the ISM (see §2 and §3).

Cosmic Rays

The final main energy reservoir of the ISM is the cosmic rays. Cosmic rays (CR) are energetic (up to 10 20 eV) particles that propagate in the ISM at relativistic speed. These particles are protons for the most part (H + ), electrons (1%), and some heavy nuclei (10% alpha particles and 1% heavy nuclei). Since CR are very energetic particles, they are an important source of ionization, they can penetrate deeper than UV photons in clouds, and therefore dissociate molecules and heat the medium. The energy spectrum of CR is very sharp: the very energetic CR are rare and have little interaction with the medium, while the low-energy CR are numerous and have large interaction cross-sections with the gas. CR with galactic origins can be generated in supernovae and can reach energies up to 10 15 eV through diffusive shock acceleration mechanism, i.e. the Fermi mechanism. OB-type stars may be the major source of the Galactic CR at low energy. These CR have small gyroradii and are easily deflected by the magnetic field, thus their propagation is approximated to a diffusive process (Strong et al. 2007). For these reasons, the low-energy CR distribution is often supposed to be isotropic (Draine 2011).

Overall, the contribution of this isotropic flux of CR to the energy density is about (Webber 1998) P CR ∼ 1.8 eV cm -3 ∼ 2.9 × 10 -12 erg cm -3 .

(1.5)

Equipartition

As shown in the previous sections, all the energy reservoirs presented above have comparable energy densities, which are also comparable to the thermal pressure of the gas

P th ∼ nk B T ∼ n H 30 cm -3 T 100 K × 10 -12 erg cm -3 . (1.6)
The physics of complex systems is often described by considering only the dominant processes and neglecting, at least in first approximation, all the other processes. The fact that all the energy densities derived above are comparable indicates that they could all play an important role in the structure and evolution of the interstellar matter, and should therefore not be neglected. Considerations on the energy densities alone are, however, not sufficient to draw a definitive conclusion. The importance of any energy reservoir depends on its capacity to interact with the matter and the corresponding interaction rates. These aspects will be described in more detail in the following chapters.

Chapter 2

Thermodynamics

The neutral diffuse ISM is not an isolated system but a system that exchanges energy with its surrounding environment, principally through the absorption and emission of photons. The interaction of this radiation and the gas is responsible for several physical processes taking place in the interstellar medium. These microscopic processes influence the dynamics of the medium and its macroscopic structure. Indeed, observations show that the diffuse ISM is not homogeneous in density and temperature but presents different phases (see §1.2) which are a consequence of thermodynamical processes. In the next sections, the different heating and cooling mechanisms are presented and modeled with formulae taken from Wolfire et al. (1995) and Wolfire et al. (2003), and shown in Fig. 2.3. The thermal processes indicated with filled triangles ( ) are those taken into account in the RAMSES code (see §5) used to run all the simulations used in this work.

Heating Processes

The diffuse ISM is rich in photons that can photoionize atoms, molecules and grains, and produce electrons with an excess kinetic energy that heats the medium through collisions.

This heating process is mainly due to UV photons coming from OB stars. Since atomic hydrogen is abundant in the ISM and has an ionization potential of 13.6 eV, a neutral phase can develop. This neutral ISM is opaque to photons with energy greater than 13.6 eV. Photons less energetic can interact and ionized heavier species (ex: C, Mg, Si, Fe..) and dust that have lower ionization thresholds.

Photoelectric Effect on Grains

The process of ejection of an electron due to an incident radiation field is called the photoelectric effect. The UV photon, propagating through the diffuse ISM, frees an electron (often called photoelectron) from the dust grain surface. This electron has an excess energy that is the difference between the photon energy and the ionization potential (which is the energy barrier for the detachment of an electron). The surplus of energy is given to the surrounding environment via thermalization through collisions.

The photoelectric effect is a quantum process. Its probability to happen depends on the charge and the surface of the grains. The more positively charged are the grains the largest is the ionization potential that needs to be overcome by the incident photon. This has two negative effects on the heating: (1) a smaller fraction of photons have sufficient energy to ionize the grain, and (2) the ejected electrons have a smaller excess energy.

The heating rate of this process can be written as follow:

Γ ph = 1.3 × 10 -24 ǫ G eff erg s -1 (2.1)
where the heating efficiency is

ǫ = 4.9 × 10 -2 1.0 + [κ/1925] 0.73 + 3.7 × 10 -2 (T /10 4 K) 0.7 1.0 + [κ/5000] (2.2) and κ = G eff T 1/2 n e φ PAH .
(2.

3)

The efficiency of this heating process depends on the incident UV radiation field, G eff (seen in more details in §3), the temperature T and the recombination parameter of electrons on PAHs, φ PAH , discussed in Wolfire et al. (2003) and set to 0.5. The electronic density, n e in the above formula, can be written at ionization equilibrium as (Wolfire et al. 2003):

n e = 2.4 × 10 -3 ζ CR 10 -16 s -1 1/2 T 100 K 0.25 G 1/2 eff φ PAH + n H x C + cm -3 (2.4)
where ζ CR is the total ionization rate including primary and secondary ionization of H by EUV photons, soft X-ray photons, and cosmic rays (set to ζ CR = 1.3 × 10 -16 s -1 ), x C + is the abundance of C + relative to n H , x C + = n(C + )/n H , which is set to 1.4×10 -4 , i.e. its Solar Neighborhood value assuming a 40% depletion of the carbon onto grains and that the carbon remaining in the gas phase is singly ionized. This equation for the density of electrons differs from Wolfire et al. (2003) by the addition of C + which is the most abundant ion in the diffuse and transluscent cold neutral medium (Snow & McCall 2006). Since the size distribution of the grains goes as n(a) ∝ a -3.5 (see §1.1), the majority of the surface that undergoes the photoelectric effect is provided by small grains and PAHs.

The efficiency of the photoelectric effect is often shown as a function of the parameter κ for different temperatures (see Fig. 2.1). But since κ and T are not two independent variables, and since we described the electronic density also as a function of the temperature, we also plot in Fig. 2.1 the photoelectric efficiency ǫ as a function of κ taking into account the dependence of this parameter on the temperature κ(T ), for two different density and a G eff = 1 in Habing units.

Cosmic Rays and X-rays ionization Cosmic rays and X-ray photons are important contributors to the heating. By ionizing the main constituents of the gas (i.e., H, H 2 , and He which has a large cross-section with X-rays), CR and X-rays free electrons. These electrons have an excess energy which has a wide distribution with an average of 30-35 eV. The freed electron can either produce secondary ionizationsimportant in chemistryor thermalize, heating the gas mainly through electron-electron collisions. X-rays decrease as the column density increases, so the photoionization by X-rays becomes less important in shielded regions. In contrast, CR can penetrate and ionize the gas effectively in those regions. The resulting heating can be written as

Γ CR ∼ 10 -27 ζ CR 10 -16 s -1 erg s -1 .
(2.5)

△ Mechanical Energy Macroscopical movements of the gas inject mechanical energy at large scales that cascades down to small scales where it is dissipated (see §5.3). This dissipation contributes to the heating of the gas, but since the medium is magnetized, nonadiabatic, and highly compressive, mechanical heating is difficult to evaluate. In the past, several theoretical efforts have been made to try to model the turbulent heating rate. On average, the turbulent heating rate can be deduced from the Larson's law (Pan & Padoan 2009) and modeled as

Γ turb = 3 × 10 -27 L 1pc 0.2 n ergs -1 cm -3 (2.6)
which gives a function of the scale considered. However, such a prescription for turbulent heating is a great misconception. Indeed, the turbulent heating, as turbulence itself, is highly intermittent. At a certain point in space and time, a lot of the energy can dissipate and a lot of heating can be produced. This dissipation by bursts can heat the gas through viscous dissipation and have strong consequences on the chemistry (Falgarone & Puget 1995;Joulain et al. 1998). This type of heating is not directly included in RAMSES because the physical scale required to described dissipation regions is not resolved numerically.

△ H 2 Formation and Destruction H 2 formation and destruction also release energy in the medium. The chemical reaction of formation of H 2 onto grains is exothermic and releases 4.5 eV that can heat up the surrounding gas, be partially stored in H 2 as internal energy and/or heat the grain on which H 2 forms. Since there is no consensus on the details and the fraction of energy going into these different processes, equipartition of energy is usually assumed. As a consequence, when H 2 is formed onto a grain, 1.5 eV is given to the thermal energy of the gas.

Destruction of H 2 usually happens after a radiative excitation in Lyman or Werner bands. This excitation leads 10% of the times to a radiative dissociation that releases on average 0.4 eV (see §5.4). In the other 90% of the cases H 2 de-excites to rovibrational excited states of the ground electronic state. This process is called fluorescence.

To correctly take into account these processes into the heating of the gas, it is not only necessary to compute the formation and destruction of molecular hydrogen on-the-fly in a simulation, but also to compute its excitation in order to obtain the dissociation rate. The implementation of these processes in RAMSES is an ongoing work.

Cooling Processes

Cooling processes correspond to the energy lost by the system, i.e. the amount of thermal energy converted into radiation that escapes the cloud. This conversion mainly happens through collisions, which internally excite some species followed by the reemission of photons. The main cooling processes are, therefore, line emissions (radiative cooling) from the most abundant species.

Lyman α

The most abundant species in the diffuse ISM is atomic hydrogen which has permitted dipole transitions. At high temperatures, the first electronic level of neutral hydrogen can be collisionally excited by an electron. This process is followed by the emission of a Lyman α photon that escapes the medium. The associated cooling rate reads (Spitzer 1978):

Λ HI = 7.3 × 10 -19
x e e -118400/T erg cm3 s -1 , (2.7) where electrons are the main collisional parters of H and their abundance relative to the proton density, x e , has been previously defined in eq. 2.4. The temperature for this transition is T = 118400 K which correspond to a photon of 10.2 eV and a wavelength of 1216 Å. This cooling is important at high temperatures and becomes dominant for T 10 4 K (cf Fig. 2.2).

When the temperatures become higher, the abundance of atomic hydrogen decreases, and Lyman α emission becomes less effective in cooling the gas. Even though temperature regimes over 10 5 K are not relevant in our study, we show in Fig. 2.2 a sketch of the total cooling rates for T > 10 5 K taken from Gnat & Ferland (2012).

Fine-Structure Lines Permitted electronic transitions are difficult to excite at temperatures below 10 4 K because the number of electrons with sufficient energy exponentially decreases. In this regime of temperatures, forbidden electrical-dipole transitions become the dominant coolants. Some species have a fine-structure splitting in their ground electronic state. This splitting is due to the coupling between the spin, i.e. the intrinsic magnetic momentum, and the orbital electronic momenta, i.e. the magnetic dipole due to the orbital motion of the electrons bound to the atom. The separation between the fine-structure energy levels is of the order of 10 -2 eV, which corresponds to some tens or hundreds of K, depending on the species. This small separation between the fine-structure levels makes it possible to be excited at low temperatures. Species with such fine-structure splitting are, for instance, CI, OI, CII, and SiII. The associated transitions are forbidden dielectric dipole and as a consequence, the rates are orders of magnitude lower than the electric dipole permitted transitions, e.g. Lyman α transitions 1 . Both oxygen and carbon ion, C + , which are the most abundant species after hydrogen and helium, present fine-structure levels.

Since neutral carbon has an ionization potential of 11.26 eV, most of the carbon is in the form of C + in the neutral diffuse ISM. The fine-structure of the ground state of C + includes two energy levels separated by ∆E /k B = 92 K. It follows that the radiative de-excitation of the excited fine-structure level can be an effective coolant of the medium at temperatures of a few tens of kelvins (cf Fig. 2.2). Including collisions by atomic hydrogen and electrons, the cooling rate reads:

Λ CII = 2.25 × 10 -23 + 10 -20 T 100 K -0.5 n e n H e -92/T x C + [erg cm 3 s -1 ].
(2.8)

where x C + = 1.4×10 -4 , i.e. the abundance of C + is set to its Solar Neighborhood abundance and assuming a depletion of 40% in grains (Sofia et al. 1997;Sembach et al. 2000).

Because oxygen has an energy potential higher that the Lyman limit (13.6 eV), oxygen exists mainly in its neutral form in the neutral diffuse ISM. Neutral oxygen presents three fine-structure energy levels at 0 K, 228 K and 326 K. Taking the solar abundance of oxygen and a depletion of 37% we obtain x O = 3.2 × 10 -4 (Meyer et al. 1998;Sembach et al. 2000). The cooling rate due to oxygen line emission can be written:

Λ OI = 7.81 × 10 -24 T 100 K 0.4 e -228/T x O [erg cm 3 s -1 ],
(2.9)

where collisions with the primary collider H are considered and the two terms arise from the two fine-structure excited levels of atomic oxygen. Wolfire et al. (2003) included only the transition from the first excited level at 228 K. However, oxygen has also a second excited level at 326 K. In §5.2.2, we discuss the effect of the inclusion of this level (see Fig. 5.2).

Oxygen and carbon are the dominant coolants in the CNM phase. Even though oxygen is more abundant, its fine-structure levels have a separation in energy larger than C + , making this last species the most effective coolant at low temperatures (see Fig. 5.1).

▽ Metastable Fine-Structure Lines Both CII and OI also have transitions between metastable states. The cooling rates induced by these transitions can be expressed as (Hollenbach & McKee 1989):

Λ CII,meta = 6.2 × 10 4 k B × 9/5(2.3 × 10 -8 x e × T -0.50 4 + 10 -12 )e -6.2/T 4 x C + (2.10) Λ OI,meta = 2.3 × 10 4 k B × 5/7(5.1 × 10 -9 T 0.57 4
x e + 10 -12 )e -2.3/T 4 x O (2.11) + 4.9 × 10 4 k B /7(2.5 × 10 -9 T 0.57 4

x e + 10 -12 )e -4.9/T 4 x O (2.12)

+ 2.6 × 10 4 k B /5( 5.2 × 10 -9 T 0.57 4

x e + 10 -12 )e -2.6/T 4 x O (2.13)

with T 4 = T /10 4 K. The main collisional partners are hydrogen atoms, for which we used a constant rate coefficient of 10 -12 cm 3 s -1 given by Federman (1983), and electrons (Hollenbach & McKee 1989). Because the levels are separated by several 10 4 K, these metastable transitions play a role only at high temperature. As shown in §5.2.2, these cooling rates are important only if Lyman α is not included. We, therefore, neglect these terms in our simulations.

Electron Recombination onto Charged Grains Grains and PAHs can be positively charged. When electrons recombine on a dust grain, the gas loses kinetic energy. As shown by Bakes & Tielens (1994), this cooling can be expressed as Λ rec = 4.65 × 10 -30 T 0.94 κ β x e φ PAH erg cm 3 s -1

(2.14) with β = 0.74/T 0.068 and φ PAH the recombination parameter.

As the photoelectric effect, the electron recombination depends on the parameter κ =

G eff T 1/2
n e φ PAH which takes into account the increase of the charge of the grain with the temperature, thus, the increase of Coulomb interaction. This cooling is important in all temperature regimes in the diffuse neutral ISM (see Fig. 2.2) and needs to be included in the computation of the thermal processes.

Molecular Lines Molecular species can also contribute to the cooling. Compared to atoms and ions, they have more degrees of freedom, thus they not only present electronic levels but also vibrational and rotational levels that can be excited. Similarly to the fine-structure splitting, the separation in energy of these levels is orders of magnitude smaller than for the electronic states.

▽ Molecular hydrogen does not have a permanent dipole (see §5.4 for more details). Consequently, the transitions between neighboring rotational levels are strongly forbidden (∆J must be even). The first allowed rotational transition is J = 2 → 0, which is separated by ∆E /k B = 510 K. Because of this large energy separation, this transition is extremely difficult to excite at low temperatures (T < 100 K), i.e. in regions where molecular hydrogen is the most abundant. Glover & Clark (2014) showed that the cooling is dominated by CII and OI and that molecular hydrogen radiative cooling has an impact on the thermodynamics of the medium only for the two following scenarios: (1) if the metallicity of the local ISM is 90% lower than its solar value, i.e. metals like C and O are extremely rare compared to H 2 , or (2) in molecular shocks, where the temperatures can be important enough to excite rovibrational levels of H 2 (Hollenbach & McKee 1979, 1989). H 2 cooling is very dependent on H 2 fraction, not only because the radiative cooling depends on the abundance of the species under consideration, but also because H 2 -H 2 collisions are 10 times more effective than H-H 2 collisions in exciting molecular hydrogen, as long as T < 500 K.

▽ Hydrogen deuteride, HD, has a permanent dipole moment, therefore J = 1 → 0 transition is possible. Since the corresponding separation energy, ∆E /k B , is lower the J = 1 level can be collisionally excited at lower temperatures. One could think that this molecule could be effective at cooling the gas at T < 100 K, since its formation rate is exothermic, i.e. energetically favorable, and its destruction endothermic. However, its abundance is limited by deuterium abundance, which is smaller than 2.5×10 -5 (Cooke et al. 2014). Consequently, this cooling process is less important than CII and OI even when considering a fully molecular medium, and will be ignored in our study. At high density, it could potentially be an important source of cooling, but those densities are beyond the scope of this Ph.D. thesis.

▽ Another molecule that can cool the medium is carbon monoxide, CO, which is the second most abundant molecule in the ISM. It is heavier than the two molecules presented before. Its energy levels are closer to each other, which implies that CO can be an effective coolant at low temperatures. In unshielded regions, the carbon is in its ionized form. For CO to be abundant, C + must have transitioned to C and finally to CO. This requires As shown in Fig. 2.4, the thermal equilibrium curve in the P /k B -n H space is not monotonic and therefore reveals a thermal instability process. Branches where the pressure increases with the density correspond to stable equilibrium states, while branches where the pressure decreases when the density increases correspond to unstable equilibrium states.

Let's consider an isobaric line with pressure P crit > P > P min , indicated with a green dotted line in Fig. 2.4, and let's call A, B, C the intersections of this line and the thermal equilibrium curve. If we consider an isobaric perturbation of a parcel of gas at density n B , e.g. an increase of the density n H > n B , the parcel of gas enters the cooling domain. Its temperature decreases while its density increase until it reaches the stable equilibrium state C. If the density is perturbed again, n H > n C , the parcel enters the heating zone of the graph. Its temperature increases and its density decreases back to the state C. If the initial isobaric perturbation of the parcel of gas in B corresponds to a decrease in density, i.e. n H < n B , the parcel enters the heating zone. Its temperature increases and its density decreases until it reaches the stable equilibrium state A. This cartoon vision of the condensation and evaporation of a parcel of gas, starting from the B state, allows to understand the three regimes of the equilibrium curve which correspond to the different phases of the ISM. The two stable branches at high and low temperatures correspond to the WNM and CNM respectively, while the unstable branch and the surrounding states correspond to the Lukewarm Neutral Medium (LNM), i.e. gas in transition between the two stable states.

Chapter 3

Radiation -photons

The ISM is a complex medium, in which the different components and energy sources are coupled and entangled. A striking example is the coupling of radiation and matter that influence the dynamical evolution of the medium.

The physical state of the interstellar gas is determined in large part by the interaction of the gas and dust with the photons. As seen before, the thermal state also depends on this radiation field through different processes of heating and cooling which have dynamical consequences on the structure of the medium. Moreover, photons govern the chemical and ionization states of the gas, which greatly depends on the rates of photoionization and photodissociation.

In the following sections, we are presenting the radiation field at different wavelengths, focusing on the waveband of interest in our study, the extinction of this radiation field due mainly to dust grains, and finally, we cast an eye on the sources responsible for the UV part of the spectrum, OB stars.

Interstellar Radiation Field

There are different sources of radiation in the ISM, dominating the ISRF at different wavelengths in the Solar Neighborhood, which will be presented in the next paragraph by increasing frequency. All energy densities are taken from Draine (2011) book and the global spectrum is shown in Fig. 3.1. Particular attention will be given to the far-ultraviolet (FUV) photons which are major actors in the chemistry of the medium, being responsible for the photodissociation of most of the molecular species.

Galactic Synchrotron Emission

Relativistic electrons interact with the magnetic field emitting synchrotron radiation, which is the dominating radiation in the radio waveband (ν 1 GHz). Relativistic electrons are not characterized by a Maxwellian velocity distribution, therefore the synchrotron emission is non-thermal. This emission dominates the radiation field at frequencies ν 1 GHz. Since the relativistic electrons are produced by supernovae and are accelerated through the Fermi process in shocks, the distribution of this emission can vary in space. A spatial mean of the synchrotron emission over the whole sky gives a total energy density integrated overall frequencies of 2.7 ×10 -18 erg cm -3 .

Cosmic Microwave Background

The CMB has nearly a constant temperature of 2.73 K (2.7255 ± 0.0006 K) and it is very close to emitting as a perfect blackbody. The radiation is basically isotropic. Even though it dominates the energy budget of the ISRF between λ = 600 µm and λ = 30 cm, with an energy density of ∼ 4 × 10 -13 erg cm -3 , it does not play an important role in the overall energy balance of the ISM, because the gas and the dust are typically warmer. The CMB interacts little with the ISM.

Infrared emission

Dust grains absorb UV radiation coming from young stars and re-emit in IR. This radiation dominates at wavelengths between λ ∼ 600 µm and λ ∼ 5 µm. This process is called extinction and depopulates the medium of UV photons, enriching it in infrared photons. Around 2/3 of the power of this radiation is emitted thermally by dust at λ ∼ 50 µm, and 1/3 by PAH vibrational bands. The total energy density of dust radiation is 5 × 10 -13 erg cm -3 . Old stars as well contribute to the IR radiation.

Cold stars Cold stars are responsible for radiation in the visible and near-infrared part of the spectrum. Mathis et al. (1983) modeled this emission as a combination of three diluted blackbody spectra at 6184, 6123, and 2539 K. An update used by the Meudon PDR code (Le Petit et al. 2006) reads:

I λ = 4π(2.07 × 10 -14 B λ (6184K ) + 1.09 × 10 -14 B λ (6123K ) + 1.52 × 10 -12 B λ (2539K )) (3.1)
where

B λ (T ) = 2hc 2 /λ 5 e hc/λk B T -1 (3.2) u λ d λ [erg cm -3 ] u λ cd λ [erg cm -2 s -1 ] u λ /E λ cd λ [photon cm -2 s -1 ]
scaling factor Habing (1968), Mathis et al. (1983), andDraine &Salpeter (1978).

u ν = 4πI ν c with λ in Å.
Another prescription of the UV radiation field was proposed by Habing (1968). The spectrum of this radiation field has been fitted by Draine & Bertoldi (1996) In Table 3.1 and Fig. 3.2, we compare the properties of the UV ISRF prescribed by Mathis, Draine, and Habing in a range of wavelength going from 910 Å to 1110 Å. In the RAMSES simulations (see §5) used in this work, the external radiation field illuminating the simulated gas at the integral scale is set to the Habing radiation field scaled with the parameter G 0 .

Plasma emission

The hot ionized gas is characterized by free-free (bremsstrahlung), freebound (recombination), bound-bound (heavy ions not completely ionized) transitions that emit EUV and soft X-rays (13.6 eV < hν < 10 3 eV). It is important to notice that the ISM is mostly opaque to EUV and transparent to high energy X-rays (> 500 eV). The energy density of X-rays between 1 -2 keV is of 10 -17 erg cm -3 , i.e. very small compared to the energy densities of the radiation field in other bands, nevertheless, these photons can be an important source of ionization and electrons (see §2).

Extinction

The medium is not transparent to photons. The UV photons propagating through the ISM are absorbed by dust and re-emitted at longer wavelengths (IR), dimming the light coming from the stars. This extinction of the radiation depends on many factors as the wavelength of the incident radiation, the size distribution of the dust grains, and their spatial distribution, i.e. the structure of the matter that dictates the permeability to photons. The attenuation of light at different wavelengths is characterized by the extinction curve. The extinction is a function that increases from longer to shorter wavelength. Because the attenuation of the radiation due to the scattering and absorption of photons by dust grains is greater for shorter wavelengths, the propagating light appears redder than the emitted light: this phenomenon is called reddening.

The shielding effect due to grains is crucial in the physics and chemistry of the ISM. It affects the thermodynamics, through the photoelectric effect, and the photodissociation of the chemical species in the medium. Since the permeability of photons depends on the structure and chemical composition of the ISM, radiation, chemistry, and dynamics are once again coupled.

The intensity of the radiation field at a wavelength λ, in the UV visible and near IR part of the spectrum, goes as:

I (λ) = I 0 (λ)e -τ(λ) (3.7)
where I 0 (λ) is the emitted intensity and τ(λ) is the optical depth that depend on the wavelength, λ. The extinction at a wavelength λ is defined as:

A(λ) = -2.5 log I (λ) I 0 (λ) ≈ 1.086τ(λ) (3.8)
where τ is related to the size of the absorbing particles via their cross sections.

The extinction in the visible, A V , can be parametrized by the ratio of extinction in the V band1 to the reddening:

R V ≡ A V A B -A V = A V E (B -V ) = 1 τ B /τ V -1 . (3.9)
R V depends on the dust properties and varies considerably from one line of sight to another. In the Solar Neighborhood the typical value is R V = 3.1 (Schultz & Wiemer 1975;Fitzpatrick 1999). E (B -V ) is the color excess is defined as:

E (B -V ) = (B -V ) obseved -(B -V ) intrinsic (3.10)
Since the color excess depends on the amount of dust on the line of sight, it correlates with the proton column density (Bohlin et al. 1978)

N H E (B -V )
∼ 5.8 × 10 21 atoms cm -2 mag -1 .

(3.11)

OB Stars

As mentioned before, FUV photons correspond to the most relevant part of the spectrum for our study. UV radiation mainly comes from young stars. These stars belong to two spectral types called O and B. They are massive stars belonging to the main sequence characterized by high temperatures and short lifetimes. Their ionizing radiation changes the chemistry of the surrounding environment, injecting energy, and ionizing and heating the entire diffuse interstellar gas.

OB stars' spatial distribution is known to be in clusters, called OB associations, that disperse in time. Ambartsumian (1947) defined OB associations as regions in which the density of this type of stars is higher than the average field. Not gravitationally bound, these associations are expanding and dispersing in time.

As shown by the recent 3D studies of the distributions of stars based on the Hipparcos and GAIA Catalogs (e.g. Bouy &Alves 2015 andZari et al. 2018), the typical distances separating two associations in the Solar Neighborhood range between 50 pc and a few hundreds of pc, i.e. several times the mean distance deduced from the integrated surface densities of OB stars (∼ 1.6 × 10 -3 pc-2, Maíz-Apellániz 2001). Interestingly, such distances are not only comparable to the heights of the molecular gas (∼ 75 pc) and the cold HI gas (∼ 150 pc) above the galactic plane deduced from CO and HI all-sky surveys (e.g., Dame et al. 2001, Dickey & Lockman 1990, and Kalberla & Kerp 2009), but they also correspond to the typical size of HI superclouds (Elmegreen & Elmegreen 1987).

Chapter 4 Observations

Thanks to the development of a great variety of telescopes during the past century by the international community, the diffuse ISM has been observed at all wavelengths, not only in our Galaxy but also in extragalactic environments. This accumulation of data has led, over time, to a more complete view of the structure, composition, and evolution of the diffuse gas. In this context, spectroscopic observations are of particular importance. By giving access to atomic and molecular spectra, spectroscopic observations provide individual measurements, and sometimes even cartographies, of the chemical composition of the ISM. In turn, and because the chemical and excitation processes are highly dependent on local physical conditions, each observation provides formidable constrains on the thermodynamical properties of the gas.

As shown by Gudennavar et al. (2012), the local diffuse ISM has been observed in absorption or emission over thousands of lines of sight, at all Galactic longitudes and latitudes, with more than 20 different atomic and molecular tracers. The complete description of all observational data, chemical correlations, and theoretical deductions are far beyond the scope of the present thesis. In this chapter, we, therefore, focus on three different datasets which are at the center of our entire work: the observations of atomic and molecular hydrogen, the observations of the excited levels of the neutral carbon, and the observations of the methylidyne cation CH + .

Atomic and Molecular Hydrogen

Since hydrogen is the most abundant element in the ISM, the observations of HI and H 2 are of great value to characterize the structure and morphology of the diffuse matter. Moreover, because the conversion of the gas from atomic to molecular is a fundamental step to trigger a more complex chemistry, direct measurements of the molecular fraction is a primary goal for any coherent study of the diffuse gas.

In the past, extensive observations of molecular hydrogen have been done in absorption, starting from its first detection (Carruthers 1970), followed by short-term missions, such as the HUT (the Hopkins Ultraviolet Telescope, e.g., Blair et al. 1996), ORFEUS (the Orbiting and Retrievable Far and Extreme Ultraviolet Spectrometer, e.g., Barnstedt et al. 2000), and IMAPS (e.g., Jenkins & Peimbert 1997). At the beginning of the 70s, a large survey on molecular hydrogen was performed with the Copernicus space mission (1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981), (e.g. Savage et al. 1977;Bohlin et al. 1983), followed by the FUSE mission (Far Ultraviolet Spectroscopic Explorer, e.g. Rachford et al. 2002b;Gillmon et al. 2006;Rachford et al. 2009), which was capable to observe translucent clouds thanks to its high spectral resolution.

In this section, we present the set of lines of sight where both the total hydrogen column density N H and the molecular hydrogen column density N (H 2 ) have been observed. The goal is to study the combined distribution of these two quantities and derive observational conclusions regarding the molecular fraction. If HI and H 2 are observed along the same line of sight, N H is simply computed as

N H = N (H) + 2N (H 2 ); if not, N H is derived from the reddening E (B -V ) as N H = 5.8 × 10 21 E (B -V ) cm -2 ,
assuming a standard Galactic extinction (Bohlin et al. 1978).

Observational Results

Column Densities All the column densities of H 2 presented below are obtained through direct observations of its FUV absorption lines. The sample results from a combination of the Copernicus survey of nearby stars (e.g., Savage et al. 1977;Bohlin et al. 1983) and the FUSE survey of the Galactic disk and the Galactic halo (e.g., Rachford et al. 2002a;Lehner et al. 2003;Cartledge et al. 2004;Pan et al. 2004;Gillmon et al. 2006;Jensen & Snow 2007a,b;Rachford et al. 2009) toward distant stars or even extragalactic sources.

In contrast, the HI column densities are derived from both direct methods, which include Lyα absorption studies, EUV observations of stellar spectra, and observations of the 21cm line, and indirect methods, which include the curve of growth of neutral and singly ionized atoms and optical interstellar absorption lines of NaI which are both found to correlate with N (H) (de Boer et al. 1986;Ferlet et al. 1985). The column densities of HI, therefore, result from a combination of FUV and optical studies for nearby stars (e.g., Bohlin et al. 1978Bohlin et al. , 1983;;Diplas & Savage 1994;Fitzpatrick & Massa 1990;Jensen & Snow 2007a,b ) and radio studies of the 21cm line for extragalactic sources at high Galactic latitude (Wakker et al. 2003;Gillmon et al. 2006).

Finally, the color excess E (B -V ) results from a compilation which includes direct measurements of the star reddening compared to its intrinsic (B -V ) 0 color (e.g., Savage et al. 1977;Fitzpatrick & Massa 1990;Diplas & Savage 1994;Rachford et al. 2002a) and dust emission maps from the IRAS telescope1 (Schlegel et al. 1998). magnitude over the entire sample. (3) The column density of H 2 measured in the Solar Neighbourhood appears to be highly dependent on the total column density: while N H varies over three orders of magnitude, N (H 2 ) varies over eight orders of magnitude. Finally, (4) the molecular fraction, computed as

f H 2 = 2N (H 2 ) N H (4.1)
reveals a sharp evolution as a function of N H : starting from a value of ∼ 10 -5 for N H < 10 20 cm -2 , the molecular fraction suddenly "jump" by three to four orders of magnitude at Gillmon et al. 2006). Because of this sharp and bimodal evolution of the molecular fraction, Fig. 4.2 is often referred to as "the HI-to-H 2 " transitions and is interpreted as a signature of a transition from atomic to molecular hydrogen, and more generally from atomic to molecular gas.

N H ∼ 3 × 10 20 cm -2 (
In addition to this chemical information, Fig. 4.2 also contains several statistical information. To highlight these statistical features, we divide the observational sample into 5 subsamples A, B, C, D, and E, as follows A) 5 × 10 18 N H 5 × 10 19 , 5.6 × 10 -8 f H 2 10 -3 which encompass almost all the observational points and whose statistical properties are summarized in Table 4.1. 48 lines of sight out of 360 are found to be not detected in H 2 . Most of these upper limits are obtained for a total column density smaller than 10 21 cm -2 and about half of them provide strong constraints on the molecular fraction with f H 2 10 -5 . 3% of the 312 detections belong to subsample A, 13% to subsample B, 16% to subsample C, and 65% to subsample D. Interestingly, the subsample E is empty and no line of sight is observed with N H > 10 22 cm -2 . While the mean value of the logarithm of the molecular fraction strongly increases from subsamples A to D, the dispersion simultaneously decreases by about a factor of three, probably revealing an effect of average over long distances.

B) 5 × 10 19 N H 7 × 10 20 , 5.6 × 10 -8 f H 2 10 -3 C) 1.6 × 10 20

Theoretical Interpretation and Open Question

Over the last decades, great efforts have been devoted to propose analytical descriptions of the HI-to-H 2 transition in homogeneous clouds with plane-parallel or spherical geometries (e.g., Sternberg 1988;Krumholz et al. 2008;McKee & Krumholz 2010;Sternberg et al. 2014 and references therein), compute this transition in detailed 1D chemical models assuming chemical equilibrium (e.g., van Dishoeck & Black 1986;Abgrall et al. 1992;Le Bourlot et al. 2012;Bron et al. 2014) or not (e.g., Lee et al. 1996;Goldsmith et al. 2007; Lesaffre f H 2 = 1 (Eq. 4.1). The red dashed-dotted line indicates the theoretical molecular fraction derived in an unshielded WNM-type environment with a density of 0.5 cm -3 and a temperature of 8000 K, illuminated by a UV photon flux of 10 8 cm -2 s -1 (see Eqs. 5.14 & 5.18). The regions A, B, C, D, and E defined in the main text correspond to an arbitrary separation of the observational sample.

et al. 2007), treat the chemistry of H and H 2 in subgrid models applied to simulations of galaxy formations (e.g., Gnedin et al. 2009;Christensen et al. 2012;Thompson et al. 2014;Diemer et al. 2018), or solve it in 3D isothermal or multiphase simulations of the diffuse ISM using various treatments of the radiative transfer (e.g., Glover et al. 2010;Valdivia et al. 2016b;Hu et al. 2016;Bialy et al. 2017;Nickerson et al. 2018). These studies led to the following conclusions. At the scale of a homogeneous cloud, the molecular content, the sharpness of the HI-to-H 2 transition, and the asymptotic column density of HI are controlled by the ratio of the intensity of the UV field to the gas density and the dust-to-gas ratio, or equivalently, the metallicity (Sternberg et al. 2014). At larger scales, the integrated column densities of HI and H 2 also depend on the distribution of clouds of various densities along the line of sight and on the porosity to the UV radiation field. Because of these effects the total column densities of H and H 2 necessarily depend on the strength, the scale, and the compressibility of the turbulent forcing which control the clumpiness of the CNM clouds (Gnedin et al. 2009;Christensen et al. 2012). Despite these achievements, very few works have been dedicated so far to the analysis of the HI-to-H 2 transition in a turbulent multiphase medium at a scale sufficient to resolve the formation of CNM structures. In addition and while the predictions of analytical models (e.g., Krumholz et al. 2008) and simulations (Gnedin et al. 2009;Valdivia et al. 2016b) were able to reproduce the trend of the HI-to-H 2 transition observed by Copernicus and FUSE in the local ISM (e.g., Savage et al. 1977;Gillmon et al. 2006;Rachford et al. 2009), the LMC and the SMC (e.g., Browning et al. 2003;Gillmon et al. 2006;Leroy et al. 2007), no detailed comparison with the statistical properties of these observations have been proposed. As a result, the occurrence of lines of sight with large molecular fractions predicted by numerical simulations often exceeds what is deduced from the observations (Valdivia et al. 2016b). Finally, and while statistical studies of 1D probability distribution functions (PDF) have become a common tool to understand the formation and the dynamics of molecular clouds (Körtgen et al. 2019), few statistical studies have been performed to date on 2D probability distribution functions using combined observations of different molecular tracers. In this thesis, we revisit the interpretation of the HI-to-H 2 transition in terms of probabilistic ordering of CNM clouds along lines of sight across a turbulent multiphase environment.

Neutral Carbon

Carbon is one of the most abundant elements in the ISM. Neutral carbon presents a finestructure splitting of the electronic ground level into three levels. The fine-structure excited levels are at energies of E /k B = 23.6 and 62.4 K and can be excited and de-excited by collisions with the other species at the typical densities and temperatures of the diffuse ISM. Since the balance between the collisional excitations and radiative decays determines the fine-structure populations, the population ratios of CI are a diagnostic of the local physical conditions, i.e. the density and the temperature of the medium containing neutral carbon. In a multiphase gas, where the cold and warm phases are believed to be in thermal pressure equilibrium, it is capital to measure the pressure. CI is an excellent probe for that.

CI can be observed in absorption in the UV part of the spectrum. The first large-survey on CI absorption lines in the UV was conducted with the Copernicus satellite by Jenkins & Shaya (1979) using hot stars as background sources. This early observational study was extented by later surveys (e.g. Jenkins et al. 1983;Smith et al. 1991;Jenkins & Wallerstein 1995;Jenkins et al. 1998;Jenkins & Tripp 2001), the largest of which was done by Jenkins & Tripp (2011). These last authors used a broad wavelength coverage in the UV spectrum of STISthe Space Telescope Imaging Spectrograph on Hubble Space Telescopewhich records many wavelengths elements at once, and is able of measuring 10 different CI multiplets simultaneously. The CI multiplets consist of complex blends of individual lines from different fine-structure levels. An example is given in Fig. 4.3, where both the original data from Jenkins & Tripp (2011) (black curves) and their reconstructions (red) are shown. Fig. 4.3 shows almost direct measurements of the apparent optical depth as a function of the radial velocity. Their spectra have a resolving power of 2.6 km s -1 (or 1.5 km s -1 for the sources of Jenkins & Tripp 2001) in radial velocity. To be able to differentiate the stellar features, from those originating from the diffuse gas, they chose background sources with large radial velocities and they rejected background sources with small radial velocities, ending up with 460 independent velocity samples.

Observational Results

Sources

The positions and distances of the 89 background sources reported by Jenkins & Tripp (2011) are shown in Fig. 4.4. Oppositely to the HI and H 2 survey presented in the previous section, most of the background sources are located toward the Galactic disk: quantitatively, 90% of the sources have a latitude smaller than 20 • . The target stars of the 21-stars survey carried out by Jenkins & Tripp (2001) are clustered around Galactic longitudes of 120 • and 300 • because they were chosen in the continuous viewing zone2 of HST, while the larger survey of Jenkins & Tripp (2011) is more evenly spread at all longitudes. The distances of the background sources range between 120 pc and about 3 kpc.

Theoretical Interpretation and Open Questions

The striking result of Fig. 4.5 is that none or few of the observational points follow the theoretical predictions, regardless of the density, temperature, or pressure of the gas. Not only most of the mass of CI is found to have an f 2 ratio far larger than the predicted value but large variations of this ratio are observed with no apparent correlation with the f 1 ratio. To solve this issue, Jenkins & Tripp (2011) interpreted this result as the presence along any line of sight of two kinds of regions: a low-pressure component (P ∼ 10 3 -10 4 K cm -3 ) contributing to most of the mass and volume of CI and a large-pressure component (P ∼ 10 5 -10 7 K cm -3 ) contributing to a small fraction of the mass and volume of CI. By analyzing each datapoint as a mixture of these two components, Jenkins & Tripp (2011) inferred the distribution of pressure in the low-pressure component rich in neutral carbon. They found that the pressure of this component weighted by the mass density, follows a lognormal distribution, with an average of log 10 P /k B = 3.58 K cm -3 . As shown by Saury et al. (2014), this distribution of pressure can be reproduced by numerical simulation of the multiphase turbulent ISM assuming a turbulent forcing with a corresponding large-scale velocity dispersion in agreement with the 1D velocity dispersion derived from HI profiles seen in emission.

These remarkable findings raise, however, a major issue: what is the origin of the highpressure component? So far, several scenarios have been proposed (e.g. Jenkins & Tripp 2011): the high-pressure component could originate from regions of dissipation of interstellar turbulence or from interstellar shocks at high velocities propagating in the neutral diffuse gas; it could be linked to the background star itself; it could be due to the recoil of H atoms following the photodissociation of H 2 at the border of molecular clouds; or, finally, it could be the trace of gravitational collapse of molecular clouds along the line of sight. In addition to this astrophysical riddle, a major question has been left out from theoretical interpretation: can the statistical distribution of the abundances of CI and of its excited states displayed in Fig. 4.5 be explained by the current astrochemical models? In this thesis, we explore these questions by computing the chemistry and excitation of CI in turbulent simulations of the multiphase ISM.

CH +

Diatomic molecules were detected for the first time in the diffuse medium at the end of the 1930s. Interestingly, the methylidyne cation, CH + , was one of the first molecules to be observed (Douglas & Herzberg 1941). Oppositely to HI, H 2 , and CI, CH + is neither a tracer of the mass, the distribution, nor the pressure of the entire volume of the multiphase ISM. At first sight, the study of such species could therefore appear to be secondary. Such considerations would be a mistake. Indeed, because of its chemical properties, the presence of CH + in the ISM is a deep mystery rooted in the evolution of the interstellar matter. In the diffuse gas, the destruction of CH + occurs through hydrogenation or dehydrogenation, two exothermic processes with very short timescales ∼1 yr (50 cm -3 /n H ). The only formation pathway efficient enough to balance these fast destruction mechanisms is the hydrogena-tion of C + , C + + H 2 → CH + + H, which is a highly endothermic reaction (∆E/k B ∼ 4620 K). This simple argument sets the mystery of CH + : its formation requires H 2 which is mostly formed in the cold neutral medium, i.e. in regions where the temperature is far too cold to activate the above endothermic reaction. The very presence of CH + in the diffuse gas challenges the validity of the predictions of astrochemical models at equilibrium. Even if it does not provide information on the global properties of the diffuse ISM, CH + can therefore be a formidable tracer of out-of-equilibrium processes induced by turbulent advection or the dissipation of suprathermal energy sources. In this section, we present the set of lines of sight where the total column density and the column density of CH + have been observed.

Observational Results

Sources

The column density of CH + across the local diffuse gas have been measured over the past 80 years through its UV absorption lines detected towards nearby OB stars (e.g. Federman 1982;Allen 1994;Crane et al. 1995;Gredel 1997;Gredel et al. 2002;Sheffer et al. 2008;Weselak et al. 2008). The distance and the position of the 190 background sources observed so far are shown in Fig. 4.6. Similarly to the observational sample of CI, about 90% of the lines of sight are directed towards the Galactic disk at latitudes smaller than 15 • . In contrast, these lines of sight are well distributed in Galactic longitudes covering homogeneously all the Galactic quadrants. The distance of the background sources ranges between 50 pc and 3.2 kpc. It follows that the resulting set of lines of sight sample the local diffuse gas located in a radius of about 3 kpc around the sun with very few data at high Galactic latitude.

CH + Column Densities

The resulting sample is shown in Fig. 4.7 which displays the column density of CH + as a function of the total column density of the gas N H . As for HI, H 2 , and CI, these observations reveal several chemical and statistical information. The column densities of CH + are found to be proportional to N H and to display a very large dispersion. The relative abundance of CH + , N (CH + )/N H , ranges over more than two orders of magnitude from 2 × 10 -10 to 8 × 10 -8 with an average of 8.9 × 10 -9 . Interestingly, the relative abundance of CH + are not stacked around the mean but appear to be homogeneously distributed in logarithmic space. Analysis of the fraction of non-detections also provides invaluable information. Indeed four regimes can be identified as a function of N H : (1) CH + is systematically not detected for N H < 1.7×10 20 cm -2 ; (2) for N H between 1.7×10 20 cm -2 and 9 × 10 20 cm -2 , there is an equipartition between non-detections and detections of CH + ; (3) for N H between 9 × 10 20 cm -2 and 2.1 × 10 21 cm -2 , the percentage of non-detections drops to around 3%; and finally, (4) CH + is systematically detected for N H 2.1 × 10 21 cm -2 . This very peculiar distribution has strong similarities with the HI-to-H 2 transition presented in §10. Because the formation of CH + requires molecular hydrogen, the range of total column density at which CH + starts to be detected corresponds to the range of column density at which the HI-to-H 2 transition occurs.

Theoretical Interpretation and Open Questions

Because CH + cannot be formed in cold environments, the predictions of static astrochemical models taking into account heating and chemical processes induced by UV photons only have been shown to systematically underestimate the abundance of CH + by two to three orders of magnitude (see Fig. 4.8). To reconcile theoretical predictions and observations several scenarios have been invoked. One possible scenario is that CH + traces the dissipation of a suprathermal energy source such as magnetized turbulence (Falgarone et al. 1995). In this scenario, the large abundances of CH + result from the out-ofequilibrium chemistry triggered in coherent and intermittent dissipative structures such as magnetized shocks (Draine & Katz 1986;Pineau des Forêts et al. 1986;Lesaffre et al. 2013;Godard et al. 2019), magnetized vortices (Joulain et al. 1998;Godard et al. 2009) or magnetic reconnection (Xie et al. 1995). As shown by Godard et al. (2014), such a scenario could explain the observed abundances of CH + providing that all the energy contained in the turbulent cascade is dissipated by ambipolar diffusion in small-scale dissipative structures (∼ 100 AU, see Fig. 4.8). Another viable scenario is that CH + could trace the turbulent mixing or the turbulent transport of warm and cold gas at the interfaces between the WNM and the CNM (Lesaffre et al. 2007;Valdivia et al. 2017). In this case, the formation of CH + would not be activated by the dissipation of a suprathermal energy source but by the outof-equilibrium H 2 transported by turbulence from the CNM to the warm LNM (Lesaffre et al. 2007;Valdivia et al. 2017). Although this scenario was shown to underestimate by a factor of five the observed abundances of CH + (Valdivia et al. 2017), no parametric study was performed on the subject.

Despite these pioneer works, the origin of CH + is still unknown. Indeed, while all the scenarios given above provide viable solutions, none really reach an unshakable consensus. On the one hand, the models describing the dissipation of suprathermal energy are limited by the fact that the dynamical structures responsible for the dissipation of interstellar turbulence are still unknown. On the other hand, the scenario advocating for turbulent transport and turbulent mixing is based on numerical simulations whose parameter domains are difficult to explore and which are ultimately limited by the spatial resolution.

In this thesis, we revisit the interpretation of the formation of CH + through turbulent transport of out-of-equilibrium H 2 from the CNM to the WNM by performing a parametric study of numerical simulations and by using simultaneously the chemical and statistical information contained in the observations of CH + . Warm H 2 A possible related issue is the presence of warm molecular hydrogen in the diffuse ISM. Indeed, since the 70s, observations of molecular hydrogen have shown that the population of H 2 in its excited rotational levels does not follow a Boltzmann distribution at the kinetic temperature T of the gas (Spitzer et al. 1973(Spitzer et al. , 1974;;Spitzer & Cochran 1973), indicating the presence of what is now called warm H 2 . This warm H 2 was later observed in cold molecular gas (e.g., Gry et al. 2002;Falgarone et al. 2005). Intermediate and high-J levels of H 2 can only be effectively populated through UV pumping or in gas with T 300 K. Since the formation of CH + is closely related to the abundance of molecular hydrogen, and requires high temperatures, the problems of warm H 2 and of the abundance of CH + Gredel (1997); Weselak et al. (2008). The CH + and total hydrogen column densities inferred from the far-infrared rotational transitions of CH + observed with Herschel/HIFI instrument and the 21 cm emission lines of HI are shown in blue and are from Falgarone et al. (2010), andGodard et al. (2012). The PDR model predictions (filled symbols, top panel) are computed for several densities: 10 (crosses), 30 (triangles), 50 (squares), and 100 cm -3 (circles). The TDR models (filled symbols, bottom panel) are computed for two different dissipation rates ǫ = 10 -24 (green) and 10 -23 (magenta) erg cm -3 s -1 , see Godard et al. 2014). could be linked.

Summary

The observables presented in the previous section, trace different physical and chemical aspects of the diffuse ISM. The hydrogen traces the total amount of matter, its molecular fraction, and its division into phases. The populations of the neutral carbon trace the pressure of the cold phases. Finally, CH + traces dynamical events taking place in small volumes and warmer phases.

In this work, we explore the effects of local physical conditions (turbulence, magnetic field, thermal instability, and gravity) on the observational tracers presented in this section. Our goal is to address the following questions. What physical processes control the HI-to-H 2 transition, the abundance of CH + , and the abundance and excitation of neutral carbon? What are the origins of the variations from one line of sight to another of the different observables? Does the distribution of lengths of the lines of sight influence the observed statistical properties? Can the observed heterogeneity of lines of sight be accounted for by the joint action of turbulence and the thermal instability?

To answer these questions we perform a parametric exploration of numerical simulations (see §5) of the multiphase ISM and compare the results to observations with methods detailed in §8.1 and §9.2.

Part II

Modeling

Chapter 5

Simulations: RAMSES

In this Ph.D. thesis, we explore the physical processes at play in the diffuse neutral ISM using state-of-the-art numerical simulations. The goal is to perform, for the first time, a large parametric study of the evolution and steady-state of a large patch of the multiphase diffuse gas as a function of its mean density, its size, the nature and the strength of the turbulent forcing, the strength of the magnetic field, and the strength of the UV radiation field.

All the simulations explored below and obtained with RAMSES (Teyssier 2002;Fromang et al. 2006), a well-known grid-based solver with Adaptative Mesh Refinement (AMR, Berger & Oliger 1984), and a fully-treated tree (Khokhlov 1998) data structure. The tree data structure is at the root of RAMSES efficiency and flexibility. In a nutshell, the grid is divided recursively into groups of 8 cells (octs) sharing the same parent cell. Each one of these octs has access to its parent cell address in memory, but also to the six neighbors of the parent cell and to its own children octs. The fully-treated tree data structure, not only makes the code suitable for parallelization with the Message Passing Interface (MPI) library, but it also allows access to distant cells by walking the tree.

In the following sections, we will address different aspects of the interstellar medium and how they are integrated into the RAMSES code. These include the heating and cooling processes (see §5.2), the turbulence (see §5.3), and the chemistry of molecular hydrogen (see §5.4). In the last section of this chapter, we will present the range of parameters studied in this work and the steady-state nature of the biphasic diffuse gas.

The diffuse matter in the Solar Neighborhood of our Galaxy is simulated over a box of size L with periodic boundary conditions. The matter, defined by a mean proton density n H , is assumed to be illuminated on all sides by an isotropic spectrum of UV photons set to the standard interstellar radiation field (Habing 1968) and scaled with a factor G 0 , and turbulence is induced through an artificial forcing.

MHD: Conservation Laws

Because the diffuse gas is a low density, the collision rates are low compared to many spontaneous rates. If we consider scales much larger than the mean free path of the species composing the medium, we can describe the medium as a fluid. The mean free path is

λ = 1 σn ∼ 10 15 n H 1 cm -3 -1 cm, (5.1)
For a cross section of an atom or a small molecule of σ = π(2 • a 0 )2 ∼ 10 -15 cm 2 with a 0 = 5.29 × 10 -9 cm, the Bohr radius of the hydrogen atom, we obtain λ ∼ 1 AU in dense clouds, and 100 AU in the WNM1 . These lengths are much shorter than the typical scales of interest in our study and thus the continuum hypothesis holds. Since the condition of the fluid approximation is verified (λ ≪ L), we can define macroscopic quantities like the mass density (ρ), velocity ( v), temperature (T ) and pressure (P) and use fluid equations (Navier-Stokes equations) and magnetohydrodynamics equations to describe the ISM dynamics.

An important dimensionless number in hydrodynamics is the Reynolds number (Reynolds 1895). This parameter is the ratio between the advection and the viscosity term in the moment conservation equation

Re = V L/ν (5.2)
where V is the characteristic velocity of the flow, L the largest scale dynamically relevant, and ν the kinematic viscosity. When the Reynolds number is small, Re ≪ 1, the fluid is dominated by viscous stresses. As the Reynolds number increases, the flows remains laminar for a wide range 2 of value Re. When Re ≫ 1, the advection term dominates over the viscous term and the flow is turbulent. It is important to notice that the advection term is highly non-linear and that the bigger it gets the more chaotic is the flow pattern. This case corresponds to a turbulent flow, characterized by strong velocity gradients and a partition of the kinetic energy on different spatial scales.

Since the interstellar medium is magnetized, the dynamical problem is described by magnetohydrodynamical (MHD) equations. In all its complexity, MHD should be treated with a multifluid approach, taking into account possible decoupling between charged and neutral particles and resistive processes. Unfortunately, such an approach would require a range of spatial and temporal scales far beyond the capabilities of modern computers.

To simplify the problem, RAMSES solves the equations3 of ideal MHD, an approximation which holds as long as the scales considered are far larger than the decoupling and resistive scales. In their conservative forms, those equations include the mass conservation,

∂ρ ∂t + ∇ • (ρ v) = 0, (5.3) the momentum conservation, ∂ρ v ∂t + ∇ • (ρ v v -B B ) + ∇P = -ρ∇Φ + ρ f , (5.4) the energy conservation, ∂E ∂t + ∇ • [(E + P ) v -B ( B • v)] = -ρ v • ∇Φ + ρ f • v -ρL , (5.5)
and the induction equation,

∂ B ∂t + ∇ • ( v B -B v) = 0, (5.6) ∇ • B = 0. (5.7)
Where ρ, v, B , P and E are the mass density, the velocity field, the magnetic field, the total pressure, and the total energy density, respectively. L is the net cooling function per unit mass (see §2), and f is the acceleration due to the turbulent driving (see §5.3).

The axis x, y, and z are chosen so that z corresponds to the direction perpendicular to the Galactic disk, and x corresponds to the direction of the mean magnetic field initially parameterized by a constant value of B x .

Gravity To take all gravitational forces into account, including self-gravity and the action of stars and dark matter, the gravitational potential Φ is divided into two terms:

Φ = φ gas + φ gal .
(5.8)

The self-gravity potential, φ gas , is deduced from the Poisson's equation:

∇ 2 φ = 4πGρ.
(5.9)

Following Kuijken & Gilmore (1989) and Joung & Mac Low (2006), we assume that the Galactic potential along the direction z perpendicular to the Galactic disk can be written as (5.10) where the first term is the contribution of the stellar disk parametrized by z 0 = 0.18 kpc and a 1 = 1.42 × 10 -3 kpc Myr -2 and the second term is the contribution of the spherical dark halo parametrized by a 2 = 5.49 × 10 -4 Myr -2 .

φ gal (z) = a 1 z 2 + z 2 0 -z 0 + 2a 2 z 2 ,

Diffusion

In the above equations of conservation, there are no terms of diffusion included, either thermic, i.e. conduction, or of momentum, i.e. viscosity. This implies that the diffusion in these simulations is entirely due to the numerical scheme. The scales at which the diffusion takes place are much smaller than the ones resolved by the problem, and therefore the numerical diffusion takes place at scales bigger than the physical scale at which viscosity and conduction act.

Thermodynamics

As shown by Field (1965), the ISM has a multiphasic nature. This results from the thermal balance of the gas and thus from its net cooling function L defined by

ρL = n 2 H Λ -n H Γ (5.11)
where n 2 H Λ and n H Γ are the cooling and heating rates of the medium (in erg cm -3 s -1 ) and n H the proton density (in cm -3 ).

To correctly describe the thermal state of the diffuse ISM, we include the heating induced by the photoelectric effect and the decay of cosmic ray particles and the cooling induced by Lyman α photons, the recombination of electrons onto grains, and the finestructure lines of OI and CII. All these processes, described in §2, are modeled with the analytical formulae given by Wolfire et al. (2003). A few tests to validate the model are presented in the next section, and a comparison with (Wolfire et al. 2003) is shown in Fig. 5.3.

Radiative Transfer

The extinction of UV photons has a direct impact on the dynamics of the medium because it affects the heating. To describe the interaction of the radiation and the matter, a treatment of radiative transfer is, therefore, necessary. However, a full radiative transfer would be too expensive numerically. Several methods have been proposed to tackle the problem (e.g. Glover et al. 2010;Valdivia et al. 2016b;Hu et al. 2016;Bialy et al. 2017;Nickerson et al. 2018). Here we present the tree-based method used in RAMSES and implemented by Valdivia & Hennebelle (2014).

To access the effective radiation field G eff in each cell of our simulation we need to compute the absorption of UV photons by dust. With the tree-based method, we can compute the column densities around each cell of the simulation, using the information already stored in the tree structure of RAMSES code. This approach has the advantage of being fast numerically because the column densities are calculated while the tree is walked. This approach is also key for the calculation of the formation and destruction of molecular hydrogen in the simulation, in which the molecular hydrogen is subject to self-shielding and, thus, the column densities of H 2 around a point affect its destruction rate.

The column densities N H (and N (H 2 ) for H 2 self-shielding presented in §5.4.3) are estimated by adding all the contributions of the cells along each line of sight with decreasing resolution as a function of the distance with the general assumption that distant cells have a smaller contribution to the local extinction.

At each point the effective radiation field G eff (in Habing units) is computed as

G eff = G 0 〈e -2.5A V 〉, (5.12)
where A V is the visual extinction along a given ray, deduced from the integrated proton column density4 computed from the border of the box to the current point (5.13) and 〈e -2.5A V 〉 is an average performed over 12 directions, treated as solid angles evenly spread in polar coordinates.

A V = 5.34 × 10 -22 N H cm -2 ,
It is also important to understand how certain assumptions made on the parameters at play in the thermodynamics, like φ PAH , ζ CR , x C + , x O , affect the thermal processes considered (see Fig. 5.2).

Since oxygen line emission is the main cooling at temperatures of a few hundred kelvins, its abundance influences the thermal equilibrium curve in the intermediate unstable branch, without affecting the stable branches. Oxygen presents two excited levels in the finestructure, one at 228 K and the other at 326 K. In the second panel of Fig. 5.2, we show how the inclusion of the second level increases the total cooling rate. The thermal equilibrium curve lowers slightly towards lower pressures.

C + line emission is an important coolant at intermediate temperatures and becomes the main coolant at lower temperatures because the excited level of its fine structure is at 92 K. Therefore, its abundance affects the cold stable branch. It is interesting to notice that, as for the oxygen, an increase in the abundance of this coolant lowers the equilibrium curve, but oppositely to the oxygen, this corresponds to an increase in the slope of the unstable branch, and consequently an increase of the pressure range over which the gas is bistable. For the oxygen the opposite happens, because the oxygen affects more the low-density end of the unstable branch lowers, while the high-density end is less affected.

The inclusion of metastable-states of O and C + has a slight impact on the high temperature range. Both metastable states and Lyman α lines have the effect of lowering the low density (high temperature) stable branch of the equilibrium curve, but since Lyman α dominates the cooling in this range of parameters (see Fig. 5.1), the impact of the metastable state on the total net cooling function is negligible and thus ignored in our study.

Wolfire takes values of the PAH recombination parameter φ PAH between 0.25, which is considered low PAH recombination efficiency, and 1.0, which is considered high PAH recombination efficiency. In Fig. 5.2 we show the effect of φ PAH = 0.25, 0.5 and 1 on the equilibrium curve. The impact of this parameter is to slightly reduce the pressure range (see Fig. 2.4) in which there is a coexistence of different stable phases, and raising this range towards larger pressures. As a recombination parameter, φ PAH plays a role in the efficiency of the photoelectric effect. Higher φ PAH implies less positively charged PAH, and in such a case the photoelectric efficiency increases, increasing the heating of the medium.

In Wolfire et al. (2003) (dashed blue line), the electronic density does not include the electrons coming from the ionization of the carbon atom. In contrast, the prescription we included in RAMSES code does (solid black line). This slightly affects the curve mainly at high densities, where a higher electronic density increases the efficiency of the photoelectric effect (see eq. 2.2).

As the reference value for the cosmic ray ionization rate, we take a fiducial value of

ζ CR = 10 -16 s -1 (see §1.3.4
). An increase in the CR ionization rate results in an increase of the heating rates. ζ CR not only affects CR heating rate, but also the efficiency of the photoelectric effect through its effect on electronic densities. As shown by the last two panels of Fig. 5.2, the thermal equilibrium curve mainly depends on the local illumination of the gas by the UV radiation field. When the UV radiation field is weaker, i.e. G 0 is smaller, or the extinction is stronger, i.e. A V is higher, the photoelectric effect onto grains decreases and the equilibrium curve moves towards lower

Turbulence

"According to an apocryphal story, Werner Heisenberg was asked what he would ask God, given the opportunity. His reply was: "When I meet God, I am going to ask him two questions: Why relativity? And why turbulence? I really believe he will have an answer for the first."

Overview on Theory

Molecular clouds show supra-thermal line widths. This non-thermal broadening is attributed to the turbulent motions of the gas. Turbulence is one of the greatest challenges in physics, that has been studied for centuries (e.g. Leonardo da Vinci). In 1922, Richardson proposed the concept of turbulent cascade: "We realize that big whirls have little whirls that feed on their velocity, and little whirls have lesser whirls and so on to viscosity." This theory was later developed by Kolmogorov (1941) who formalized it for an incompressible flow, under the hypotheses that the small-scale turbulent motions are statistically isotropic and that the statistics of those small scales are entirely determined by the kinematic viscosity and the rate of energy dissipation.

Turbulence has a heuristic definition which defines it as an instability of the laminar flow. This instability develops when the inertial forces, i.e. advection, become larger than the viscous ones. The ratio bewteen the advection and the viscous term is given by a dimensionless number, the Reynolds number (see eq. 5.2). The ISM has very high Reynolds numbers (Re ∼ 10 5 -10 7 ), meaning that the advection term in the momentum conservation is very large, i.e. the medium is highly turbulent, and that the viscosity is negligible at large scales.

Turbulence is a multiscale phenomenon that transfers energy from large scales (big eddies) to small scales (very small eddies) and it plays a central role in the evolution of the ISM and the structure formation. In the ISM, there are many sources of mechanical energy to which the medium is subject. As seen in §1.3, the mechanical energy is an important energy reservoir of the medium, with roughly the same energy density as the other energy reservoirs. Several sources can inject mechanical energy in the ISM: supernovae (Padoan et al. 2016a,b;Pan et al. 2016;Körtgen et al. 2016), stellar winds, superbubbles, expanding HII regions (Matzner 2002;Krumholz et al. 2006;Goldbaum et al. 2011), the differential rotation of the galaxy, gravitational instabilities (Federrath et al. 2011;Sur et al. 2012), jets and outflows (Li & Nakamura 2006;Wang et al. 2010;Federrath et al. 2014), and accretion (Klessen & Hennebelle 2010;Vázquez-Semadeni et al. 2010;Lee & Hennebelle 2016).

The mechanical energy is injected at large scales, and different energy sources can have different energy spectra, i.e. inject energy at different scales (Norman & Ferrara 1996). This injected energy is transferred to lower scales through the turbulent cascade until it reaches scales small enough so that the energy dissipates through ambipolar friction, viscosity, or resistivity. Many theoretical efforts have been done to characterize and describe turbulence and the turbulent cascade, and both are still debated.

In numerical simulations, there are different methods to inject turbulence. Some in-clude the source of the mechanical energy at large scales, e.g. SN feedback. Others mimic the turbulence through a mathematical prescription of a forcing term in the momentum equation.

In this work, we adopt a pseudo-spectral approach to mimic the turbulence. Oppositely to converging flows studied by Audit & Hennebelle (2005) and Valdivia & Hennebelle (2014), this prescription allows to conserve the mass in the simulations to directly compare our results to other works, such as Seifried et al. (2011) and Saury et al. (2014), and study the influence of the injection scales, frequencies, and the relative power in solenoidal and compressive modes on the structure of the matter and in the chemical evolution.

Implementation in RAMSES: Turbulent Forcing

To mimic the injection of mechanical energy in the diffuse ISM, a large-scale turbulent forcing is applied. Following Schmidt et al. (2009) and Federrath et al. (2010), this forcing, modeled by an acceleration f in the momentum conservation equation, is driven through an Ornstein-Uhlenbeck process using a pseudo-spectral method. At regular time intervals ∆τ, random fluctuations of the forcing term are generated and applied over an autocorrelation timescale τ. To excite only large-scale modes, the forcing is modeled as a paraboloid in Fourier space covering a small interval of wavenumbers 1 k 3 and centered on k = 2. Using the notations of Seifried et al. (2011) and Saury et al. (2014), the total magnitude of these perturbations is set with either an acceleration parameter F or, equivalently, a velocity parameter V related by F = V 2 /L drive , where L drive is the main driving scale, L drive = L/2. A Helmholtz decomposition is finally applied, in order to control the powers injected in compressive and solenoidal modes. Using the classical notation, these powers are set with a parameter5 ζ ranging from a pure solenoidal field (ζ = 1) to a pure compressive field (ζ = 0).

Molecular Hydrogen

In addition to the thermodynamical properties described above, RAMSES also computes the chemistry of molecular hydrogen. The timescale required for the abundance of molecular hydrogen to reach its equilibrium value is known to range over several orders of magnitude, depending on the physical conditions of the ISM (e.g., Goldsmith et al. 2007;Tabone 2018). In the diffuse gas, this timescale varies typically between a few 10 3 yr and a few 10 7 yr (Valdivia et al. 2016b), hence over a range of values that often exceeds the dynamical timescales. To take this important aspect of the diffuse interstellar chemistry into account, the out-of-equilibrium abundance of H 2 is computed self-consistently in the simulation, using the formalism introduced in RAMSES by Valdivia et al. (2016b).

In this section, we will introduce the main formation and destruction processes, and the self-shielding of H 2 against photodissociation.

Formation of H 2

In this section we will first address the different processes that can lead to the formation of molecular hydrogen, starting from the gas-phase reactions.

Radiative Association

The formation of H 2 by radiative association,

H + H → H 2 + hν,
has an extremely small rate coefficient. This is due to the fact that the excess of energy of the two initial atomic hydrogens can only be radiated away through a quadripole transition, because molecular hydrogen is homonuclear, hence has no permanent dipole.

H -Route One catalytic way to form H 2 is the H -route that happens through a radiative attachement, H + e -→ H -+hν, followed by an associative detachment, H -+ H → H 2 + e -+ KE, which is an exothermic reactions. Since ion-neutral reactions (e.g. H -associative detachment) have bigger rates than neutral-neutral reactions (e.g. H radiative association), this route is more efficient in forming molecular hydrogen than the previous one. Nevertheless, H -formed by the radiative attachment can be destroyed by a photodetachment, H -+ hν → H + e -, which has a rate proportional to the external radiation field, and therefore predominant in dense and dust-poor PDRs (Tabone et al. 2020). This route is then negligible in the diffuse medium under study in this work.

Three-body Association Three-body association (Palla et al. 1983) is another mechanism, H + H + C → H 2 + C, where C is a third body that can be hydrogen, molecular hydrogen, helium, or another abundant species of the ISM. In this formation process, the excess of energy is carried out by a third atom. The density at which three-body collisions can be important is much higher than the densities of the diffuse ISM we are modeling here, and this formation route will be neglected in our studies.

Grain Surface Catalysis

The dominant process of H 2 formation in our Galaxy is through dust-grain catalysis. The dust-grain surface acts like a catalyser for the formation of molecular hydrogen. There are two-pathways of this process: the physisorption (or the Langmuir-Hinshelwood mechanism) in which a hydrogen atom is weakly bonded and able to wander on the surface of the grain until it encounters another physisorbed hydrogen atom, and the chemisorption (the Eley-Rideal mechanism), in which the atomic hydrogen is chemically bonded to the surface of the grain and unable to move, but can react with a hydrogen atom in the gas phase that collides with the grain.

The formation rate of these mechanisms depends on the rate of collisions between atomic hydrogen and dust grains, and on the probability of an atomic hydrogen to stick to the surface of the dust grain. The formation of H 2 onto grains in physisorption and chemisorption sites is the main formation processes of H 2 in the diffuse ISM and it is modeled in RAMSES with the simplified rate of Le Bourlot et al. (2012) (5.14) where n H and n(H) are the local proton and atomic hydrogen densities, and

k f = 3 × 10 -17 n H n(H) T 100 K S(T ) cm -3 s -1 ,
S(T ) = 1 1 + ( T T 2 ) β (5.15)
is the sticking coefficient of H onto grain, parametrized by T 2 = 464 K and β = 1.5 (Bron et al. 2014).

Destruction of H 2

Photodissociation of H 2 In the diffuse medium and PDRs, the principal destruction process is photodissociation,

H 2 + hν → H + H.
The photodissociation of H 2 can be direct, i.e. one molecular hydrogen is hit by a photon and dissociates in two atomic hydrogens, or it can be a two-step process, where the molecular hydrogen is excited by a photon and then it decays to the vibrational continuum.

The direct photodissociation from X ground state to b unbound state (H 2 → H(1s)+H(1s)) is highly forbidden by selection rules (∆S = 0), while the direct dissociation from X electronic ground state to B or C vibrational continua requires a minimal energy of 14.7 eV. Since the diffuse ISM is rich in atomic hydrogen, whose ionization potential 13.6 eV, this photodissociation is negligible because the medium is not rich in photons with hν > 14.7 eV capable of directly photodissociating H 2 .

The two-step photodissociation process, happens through a first step in which the molecular hydrogen absorbs a UV photon and goes from its ground electronic state to the Lyman band (if hν > 11.2 eV, i.e. photons with λ < 1108 Å, see Fig. 7.2) or the Werner band (if hν > 12.3 eV, i.e. photons with λ < 1008 Å, see Fig. 7.2), and a second step that can lead to a transition/decay to a vibrational excited level (∼ 85-90% of the cases) of H 2 ground state or to the vibrational continuum of the ground state (∼ 10-15% of the cases), leading to the dissociation of H 2 . Both these 2 nd steps imply an emission. Since the diffuse ISM is rich in UV photons, this photodissociation process is the main one in our domain of study.

Photoionization Another process that destroys hydrogen molecules is the photoionization,

H 2 + hν → H + 2 + e -.
Since the photoionization threshold for molecular hydrogen is 15.4 eV (see Fig. 7.2), photons with hν > 15.4 eV can ionize H 2 , destroying it. The diffuse ISM is opaque to such photons, and this destruction route of H 2 is negligeble in our study.

Collisional Dissociation of H 2 H 2 can also be destroyed by collisions with other species. Dalgarno & Roberge (1979) showed that the collisional dissociation can happen at both high temperatures and high densities (e.g. T > 4000 K and n H > 10 4 cm -3 ). Thus, this type of dissociation is negligible in the diffuse interstellar medium.

H 2 Shielding

Since the main process of destruction of H 2 is photodissociation by UV photons, the absorption and scattering of UV photons is a key ingredient in the computation of molecular hydrogen. The extinction of UV photons can lower the efficiency of the photodissociation of H 2 in certain regions, leading to an increase of molecular hydrogen. The UV radiation field is attenuated by dust, as explained in §3.2, but also by molecular hydrogen itself. H 2 molecules absorb UV photons in lines. These molecules are, thus, sources of opacity for photons at those wavelengths. In the inner regions of a cloud, the external radiation field is attenuated by all of the H 2 molecules in the outer regions of the cloud, leading to the so-called 'self-shielding' of H 2 . Different methods exist to compute the shielding of H 2 due to dust and to molecular hydrogen itself. While an exact radiative transfer can be computed in PDR codes, this approach is very expensive numerically. As a consequence other approximative approaches are often used, such as the FGK transfer (Federman et al. 1979) , which considers the absorption at the line center to compute the radiative transfer. In this work we use an analytical prescription to compute the self-shielding of H 2 (Draine & Bertoldi 1996):

f shield (x) = 0.965 (1 + x/b D ) 2 + 0.035 e -8.5×10 -4 1+x 1 + x , (5.16)
where x = N (H 2 )/5 × 10 14 cm -2 (5.17) and b D is the Doppler broadening parameter expressed in km s -1 . This prescription needs the amount of H 2 traversed by UV photons to compute an estimate of the self-shielding.

In RAMSES, this quantity is computed with a tree-based method developed by Valdivia & Hennebelle (2014) to provide a fast estimate of H 2 column densities along 12 different directions. These amounts are then averaged to obtain f shield (x) and compute the photodissociation rate of Eq. 5.18. Valdivia & Hennebelle (2014) method uses solid angles to compute the intensity of the incident external radiation field in each point in the simulation, taking into account the attenuation due to dust, e -σ d N H , and due to the self-shielding coefficient, f shield . For the shielding induced by dust, we adopt here an effective dust attenuation cross section at λ = 1000 Å, σ d = 2 × 10 -21 cm 2 (Sternberg et al. 2014). Following all these considerations, the photodestruction rate of H 2 is modeled as

k d = k 0 d n(H 2 )G 0 〈e -σ d N H 〉 f shield (x) cm -3 s -1 , (5.18)
where k d ,0 = 3.3×10 -11 s -1 is the inverse free-space dissociation timescale of H 2 in an isotropic Habing field (Glover & Mac Low 2007), n(H 2 ) is the local density of the molecular hydrogen, e -σ d N H is the shielding induced by dust and f shield the self-shielding function.

Note that the prescription of H 2 self-shielding used in this work (Eq. 5.16) is taken from Draine & Bertoldi (1996). As discussed by Wolcott-Green et al. (2011), such a prescription is reliable for diffuse gas at low temperature but becomes less and less reliable for high-temperature environments (T > 500 K) where efficient collisional excitation of H 2 in its rovibrational levels reduces the self-shielding. In Fig. 5.4, we compare the analytical self-shielding at the scale of a CNM cloud. The impact of this broadening parameter on the comparison of the HI-to-H 2 transition with the observations will be discussed in §10.

Grid of Simulations

Steady-State

All simulations start at an homogeneous density and temperature of T = 8000 K. The gas evolves under the joint actions of turbulence, gravity, magnetic field, radiation field, and thermal instability, following the MHD equations given in §5 over several tens of megayears. The outputs are saved periodically to be able to access the physical quantities in every cell of the cube at different times. The steady-state is characterized by a plateau reached by all mean physical quantities. Because of the turbulence and the thermal instability, the steady-state has a statistical nature: while the maps of the column densities can appear very different from one time to another, the PDFs remain constant.

During its evolution, the gas splits up in three different phases at thermal pressure equilibrium, the WNM, the CNM, and the LNM (see Fig. 5.5). The formation of dense environments well shielded from the destructive UV radiation field triggers the formation of H 2 and the medium jointly evolves from a purely atomic state to a partly molecular state (see the cyan line in Fig. 5.6). If the mass of the gas is conserved, as imposed by the periodic boundary conditions, the medium progressively tends towards a steady-state where the mean pressure (see the red curve in Fig. 5.6), 〈P /k B 〉, the mass fractions of the different phases, f WNM and f CNM , their velocity dispersion, σ turb , and their mean molecular fractions are roughly constant. This steady-state is typically reached after a few turnover timescales, provided that the corresponding time is longer than the damping time, i.e. the typical timescale of decay of the turbulence.

The turbulent forcing and the subsequent turbulent cascade induce pressure variations and shear motions at all scales which trigger mass exchanges between the different phases. Any pressure or density structure is therefore a transient system which is usually described by its contribution to probability distribution functions. This steady yet everchanging environment is the reason for the sustained presence of a substantial amount of gas in the LNM at densities and temperatures out of thermal equilibrium (e.g. Marchal et al. 2019). All the results shown throughout this work are taken at steady-state, i.e. at times ranging from a few tens of Myr up to 100 Myr depending on the strength of the turbulent forcing.

Parametric Study

A parametric study, covering a wide range of physical conditions, is conducted to explore the impact of different physical parameters and unveil the joint effects of turbulence, magnetic field, radiation field, and gravity on the structure and the evolution of the ISM and the degeneracies between all these components, with the help of comparisons with different observables (e.g. HI-to-H 2 transitions, pressure distribution and correlations between different chemical species, see §4).

The framework, described in §5, lays on several independent parameters, which are all related to key physical ingredients of the ISM. The influence of each ingredient is studied here through several grids of simulationsincluding a total of 305 runscovering a broad range of physical conditions and centered around a fiducial setup6 . The reference value adopted for each parameter, and the range of values explored in this work, are summarized in Table 5.1. Among all parameters, L, n H and G 0 are of particular importance.

With our assumptions, L simultaneously corresponds to the scale of illumination of the gas by UV photons and twice the integral scale of turbulence, that is twice the scale of injection of mechanical energy L drive . Taking into account the distances between of OB stars associations in the Solar Neighborhood, the height of molecular gas above the Galactic plane, and the typical size of HI superclouds discussed in §3.3, we adopt a fiducial simulation with L = 200 pc and explore values down to a few tens of pc.

The mean density of the gas, n H , represents the mass of the diffuse neutral ISM contained in a volume L 3 , and also controls the porosity of the matter to the impinging radiation field. In this work, we adopt a fiducial value n H = 2 cm -3 , a value slightly larger than the standard Galactic midplane density of HI at a galactocentric distance of 8.5 kpc (Kalberla & Kerp 2009). With our fiducial value of L, the mass of gas contained in the box accounts for all the mass surface density of HI in the Solar Neighborhood (Σ HI ∼ 10 M ⊙ pc -2 , Nakanishi & Sofue 2016; Miville-Deschênes et al. 2017).

G 0 controls the intensity of the radiation field, hence the thermodynamical state of the gas. Since G 0 is normalized to the Habing field, we choose a fiducial value of G 0 = 1. The corresponding UV energy density of the fiducial setup is therefore slightly smaller than that contained in the standard UV radiation fields given by Draine & Salpeter (1978) and Mathis et al. (1983).

The fiducial values of the two parameters F and B x are set to 9 × 10 -4 kpc Myr -2 (see vertical dotted line in Fig. 5.11) and 3.8 µG, respectively. Those values are chosen so that the velocity dispersion of the gas and the strength of the magnetic field are close to the values observed in the diffuse ISM.

Physical State of the Gas

Mean Pressure As shown in Fig. 5.7, the mean pressure of the gas is primarily regulated by the thermal equilibrium curve (see §2) which depends on the local illumination of the gas by the UV radiation field G eff . Since G eff results from the absorption of the external UV field by the surrounding environments, 〈P /k B 〉 is not only sensitive to G 0 but also to the total mass of the simulation set by n H and L. Larger values of G 0 or smaller values of L or n H leads to larger 〈P /k B 〉.

The turbulent forcing induces pressure fluctuations and shearing motions at all scales.

the turbulent forcing, we introduce the velocity dispersion. Scaling laws link the velocity dispersion to the size of molecular clouds (Larson 1981;Hennebelle & Falgarone 2012).

The velocity dispersion of the WNM is mostly given by the strength of the turbulent forcing and the driving scale (L drive ∼ L/2), with a slight dependence on n H and ζ. As proposed by Saury et al. (2014), a realistic value for the turbulent velocity dispersion of the WNM can be estimated by looking at the HI 21 cm emission spectra with the fewest components observed at high Galactic latitude (Kalberla et al. 2005). Towards these directions, Haud & Kalberla (2007) derive a total velocity dispersion σ tot = (σ 2 turb + σ 2 th ) 1/2 ∼ 10 km s -1 , where σ th is the 1D thermal velocity dispersion (∼ 8.2 km s -1 for the WNM). In the present work, the turbulent forcing applied to the fiducial simulation (see Table . 5.1 and Fig. 5.11) is chosen so that σ turb ∼ 4 -6 km s -1 (see dotted lines in Fig. 5.11), in fair agreement with the observations at high Galactic latitude. While this value is chosen as a reference, the velocity dispersions obtained in all the simulations explored in this work range between 1 and 15 km s -1 .

There are many ways to define the turbulent velocity dispersion σ turb . It could be defined as the volume-weighted or mass-weighted velocity dispersion computed over the entire volume (Audit & Hennebelle 2005;Federrath et al. 2010), the average of the massweighted velocity dispersion computed along individual lines of sight (Miville-Deschenes & Martin 2007; Saury et al. 2014), or the dispersion of the mass-weighted velocity centroids computed along independent directions (Henshaw et al. 2019). All these definitions give velocity dispersions that differ from one another by factors of a few and are not equally relevant for the comparison with observed quantities. To relate the velocity dispersion to a kinetic energy and provide values that could be compared to the observations of broad HI emission profiles at high Galactic latitude, we choose to compute σ turb in the WNM only as

σ 2 turb = 1 3 ρ|| v -v|| 2 ρ (5.21) with v = ρ v ρ . (5.22)
Correlation Time and Damping Time In Fig. 5.10, we show the impact of different correlation times on the 1D turbulent velocity dispersion for different box sizes for F = 1.5 × 10 -4 and 1.5×10 -3 kpc Myr -2 . The curves are approximately flat as a function of the correlation time, and they depend more strongly on the driving strength and the size of the box. Theoretical studies and numerical simulations (McKee & Ostriker 1977;Korpi et al. 1999;Mac Low & Klessen 2004) show that the main source of turbulence is the mechanical energy injected by supernovae explosions. To follow this trail, we set the time scale to ∆τ ∼ 0.4 Myr, which is the typical interval separating two SN explosions in a box of (200 pc) 3 . The characteristic damping time of the turbulence τ is approximately set to the turnover timescale of the diffuse ISM τ ∼ 33(L/200pc) 0.6 Myr (Larson 1981;Hennebelle & Falgarone 2012).

Driving Strength In Fig. 5.11, we show how different box sizes L and different driving strengths, F , affect the turbulent velocity dispersion. These 15 simulations are run with- out magnetic field and with a mean density n H = 1 cm -3 and an external radiation field of G 0 = 1. The turbulent velocity dispersion is computed following eq. 5.21.

Another parameter left free in this work is ζ, that imposes the fraction in solenoidal and compressible modes. In Fig. 5.12, we show the dependence of the turbulent velocity dispersion in 1D to the driving strength F for three different compression parameters, ζ = 0.1, 0.5, 0.9. In this case, instead of using the acceleration parameter, F , we trace the velocity dispersion as function of the velocity parameter V = F L drive , to show how the velocity dispersion scales roughly linearly with this parameter. These 15 simulations are run on a box of L = 200 pc, i.e. the main driving scale is L drive = 100 pc, and all the rest of the parameters (apart from F and ζ) are set to their reference value in Table 5.1. Compare to the driving strength and the size of the box, the effect of ζ is small on the velocity dispersion.

Chapter 6

CI Excitation

Carbon is an abundant element in the ISM, and one of the most abundant metals1 with oxygen. In the diffuse neutral ISM, carbon is mainly carried by three different species: C, C + , CO. Since the populations of the fine-structure levels of CI are often used as a measure of pressure (see §4), we describe in this section its energy structure and transitions along with the code we built to compute the excitation of neutral carbon in post-processing.

The neutral carbon atom has an electronic ground state configuration of 1s 2 2s 2 2p 2 : the two outer electrons are in p orbitals, meaning that they have an orbital momentum l = 1. These two orbital momenta build up to a total orbital momentum, L, that can have values 0, 1, or 2. Both these electrons have a spin of 1/2 which can build up a total spin momentum S of 0 or 1. From Pauli's exclusion principle the total wavefunction needs to be asymmetric with respect to the exchange of electrons, which is equivalent to say that either both S and L are odd or even. Looking at all possible combinations gives the electronic levels 1 S 0 (S = 0, L = 0), 1 D 2 (S = 0, L = 2) and 3 P 0,1,2 (S = 1, L = 1). The superscript indicates the degeneracy of the level due to the total spin momentum (subscript = 2S + 1). Both 1 S and 1 D are singlets states, while 3 P is a triplet that shows a fine-structure splitting (see Fig. 6.1).

Electronic Ground State: Fine-Structure Levels

A good quantum number to differentiate the fine-structure levels is J = |L +S|, which in this case has 3 possible values: J = 0, 1, 2, as shown in Fig. 6.1. In the following, we call the first excited level CI * ( 3 P 1 ) and the highest level CI * * ( 3 P 3 ). Because they lie at energies 23 K and 62 K, these levels can be populated collisionally at low temperatures.

In §2, we have already discussed the fine-structure transitions (the one responsible for cooling at T < 10 4 K) and pointed out that these transitions are highly forbidden. Selection rules 2 in atomic physics dictate that allowed electronic dipole transitions are given by: ∆S = 0 (6.1)

∆L = ±1 (6.2) ∆J = 0, ±1 (6.3)
All CI fine-structure levels have the same orbital momentum L = 1: the transitions are therefore forbidden. It follows that absorption of radiation and induced emission are probably negligible processes in the excitation and de-excitation of the fine-structure levels.

CI Excitation

Radiative and Collisional (De-)Excitation

Collisions with other species present in the medium can excite and de-excite the atom and therefore change the population in the fine-structure levels of the ground electronic state of CI. The main collisional partners of CI are hydrogen (whatever its forms: H + , H, H 2 ), helium and electrons. The collision rate for each partner 3 can be written as n C k C i j (in s -1 ), where n C is the density of the collisional partner and k C i j is the collision rate coefficient (in cm 3 s -1 ) and depends on the cross-section, the relative velocity of the two species and their reduced mass.

Spontaneous emission can also depopulate the higher levels. An atom can emit a photon and radiatively de-excite with a rate which is given by the Einstein coefficient of spontaneous emission A uℓ (in s -1 ).

2 Selection rules, or transition rules constrain the possible transitions from one quantum state to another of a system 3 indicated with the letter C for collisional partner 70

At steady-state d n i d t = 0 , we can write

n i j =i n C k C i j = j =i n j n C k C j i + j >i n j A j i . (6.4)
Dividing by the abundance of the element and putting all terms depending on the abundances on one side we obtain

AX = B

with:

A = M rad + M coll ; X =    x 0 x 1 x 2    ; B =    0 0 0    ; M coll =    -C 01 -C 02 C 10 C 20 C 01 -C 10 -C 12 C 21 C 02 C 12 -C 20 -C 21    ; M rad =    0 A 10 A 20 0 -A 10 A 21 0 0 -A 20 -A 21    .
The collisional matrix is composed of the sum of the collisions for each collision partner

C i k = c n c k c i j
where the colliders are H, H + , H 2 , He, and e -. The collision excitation rates are linked to the collision de-excitation rates by the detailed balance principle:

C ℓu = C uℓ g u g ℓ e -∆E uℓ k B T (6.5)
where

∆E uℓ = E u -E ℓ .
This system of equations is not closed, i.e. the lines forming the A matrix are not linearly independent. Therefore the system needs to be closed by a conservation relation which is

i x i = 1,
or equivalently n(CI) = i n i . This equation can replace any equation of the system. For numerical reasons, we choose to introduce this closure equation instead of the steadystate equation of the level i with the largest abundance x i .

We built the excitation matrix taking into account radiative decay (Froese Fischer & Saha 1985) and collisions with atomic hydrogen (Launay & Roueff 1977), molecular hydrogen (Schroder et al. 1991), helium (Staemmler & Flower 1991), ionized hydrogen (Roueff & Le Bourlot 1990) and electrons (Johnson et al. 1987). creases until it becomes larger than f 1: this behavior makes the curve look like an upward and backward bent arch. Interestingly, in the low-pressure regime, one given pressure can result in markedly different values of f 1 and very close values of f 2. This behavior is inversed at high pressure. All these results are in perfect agreement with theoretical prediction of Draine (2011) and Jenkins & Tripp (2011). This validates the post-treatment excitation code of neutral carbon. In the following, we estimate the impact of the chemical pumping and optical pumping on the populations of the fine-structure levels of CI.

Chemical Pumping

Until now we only considered excitation from collisions and de-excitation due to collisions and spontaneous emission. However, excitation and de-excitations can also occur during chemical reactions, i.e. at the formation and destruction of neutral carbon. In the diffuse ISM, neutral carbon is mainly formed through the recombination of C + with electrons and PAHs, and mainly destroyed by photoionization.

These processes add a formation and destruction term for every level. At steady-state, the equation system is modified to

n i C j =i n C k C i j + k D = C j =i n j n C k C j i + j >i n j A j i + K i F . (6.6)
Given the large exothermicity of recombination reactions, we assume that the formation rate of a level, K i F , is proportional to the degeneracy of the level i . Assuming that all levels are destroyed at a single rate k D , the system of equations to solve reads AX = B , where

A = M rad + M coll + M D , M D =    -k D 0 0 0 -k D 0 0 0 -k D    and B =    -p 0 k D -p 1 k D -p 2 k D    .

Optical Pumping

Through the absorption of optical photons neutral carbon, can be excited to higher electronic states and then radiate to end up in one of the three fine-structure levels of the ground state. This process, referred to here as optical pumping, can be easily modeled by adding an optical pumping matrix to A which reads (Jenkins & Tripp 2011). Since we express the radiation field in Habing units, and the tabulated optical pumping values are given for a Mathis et al. (1983) radiation field, we use the conversion factor (1.27) calculated in Table 3.1 between the energy densities of these two radiation fields. Interestingly, the inclusion of the chemical and optical pumping has almost no impact on the excitation and de-excitation of the fine-structure levels of CI. This result is in line with the conclusion of Jenkins & Tripp (2011). Because of the relatively large density of the CNM and the weak energy density of optical photons in the ISRF, both processes are negligible compared to collisional excitations.

M OptPump =    -P

CI radiative transfer

Computing the excitation of CI is not enough to allow a comparison with the observations described in §4.2. Indeed, comparing simulated data with observations requires to access the relative population of CI in velocity channels. To do so, it is necessary to implement the radiative transfer and create line profiles of the excited states of CI.

The evolution of the specific intensity I ν (in erg cm -2 s -1 sr -1 Hz -1 ) can be written as

d I ν d z = -k ν I ν + η ν (6.7)
where k ν is the atomic absorption coefficient, (in cm -1 ), and η ν is the emissivity (in erg cm -3 s -1 sr -1 Hz -1 ). Assuming that there is no scattering and that the optical emission of CI is neglible, the transfer equation simplifies into

d I ν d z = -k ν I ν (6.8)
where the absorption coefficient reads

k ν = λ 2 i j 8π g i A i j n j g j φ ν i j . (6.9)
with a Gaussian line profile,

φ v i j = 1 2πσ exp - (v -v z ) 2 2σ 2 (6.10)
where σ is the total velocity dispersion, which is a quadratic sum of the 1D thermal velocity dispersion and the micro-turbulent velocity dispersion, b:

σ 2 = k B T m + b 2 .
(6.11)

Resulting Line Profiles

Our procedure to compute the line profiles of CI across simulations of the diffuse multiphase ISM is based on the chemical post-process code developed by Valdivia et al. (2017) (see §13). This code provides the atomic and molecular abundances of 151 chemical species including neutral carbon and all its collisional partners. Implementing the excitation of CI in this code along with the simple radiative transfer described above, not only gives access to the densities n i of the three fine-structure levels of CI but also their column densities in any velocity channel. In practice, the dynamical information from RAMSES output, e.g. the velocity in each direction and the temperature, are used to determine the line profile per velocity, φ v i j , defined above. The quantity n i φ v i j is then integrated along the line of sight in velocity intervals of 0.5 km s -1 (see Jenkins & Tripp 2011), to find the column density of the level i per velocity channel

N v i = v+0.25 v-0.25 l los 0 n i φ v i j d ℓd v. (6.12)
Each cell of the simulation, therefore, contributes as a gaussian to N v i with a width that depends on the temperature T of the cell, a height proportional to the density n i , and a position that depends on the velocity along the line of sight v z of the cell. In Fig. 6.4, we show an example of those variables along one line of sight of our fiducial simulation. We note that, on that particular line of sight, both CI * and CI * * exhibit two density peaks, which correspond to two different velocities at around 2 km s -1 and -13 km s -1 , and to two different phases, the CNM at low temperature and WNM with T ∼ 8000 K.

The corresponding column density per velocity channel N v CI * is shown in the left panel of Fig. 6.5 assuming a microturbulent velocity dispersion of 1 km s -1 . The range of velocity, from -30 to 30 km s -1 , is determined by the minimal and maximal velocity in the simulation and the maximal velocity dispersion, min(v z ) -max(σ) < v < max(v z ) + max(σ).

N v

CI * exhibits two peaks centered at the velocities given above: the peak at high temperature appears broader, because of the thermal broadening. The thinner colored lines are the contributions of each cell integrated in velocity bins of 0.5 km s -1 , n i φ v i j d v, and the thicker line is the integration over the entire line of sight. In the right panels of Fig. 6.5, we show other examples of profiles integrated over other lines of sight. The top panel displays a line profile with two separated peaks in velocity because the corresponding peaks of n(CI * ) along the line of sight are at velocities that differ by 10 km s -1 and belong to the CNM. In contrast, the bottom panel shows a line profile that appears as a mixture of a large component at high temperature and a narrow component at a lower temperature, with a separation in velocity similar to the previous case. All these results show that the combination of the chemical solver, the CI excitation code, and the radiative transfer code provide a powerful tool to study the column densities of neutral carbon integrated in separate velocity channels. Such a tool allow a direct comparison of the simulations with observations of neutral carbon across the diffuse neutral ISM collected by Jenkins & Tripp (2011) (see §4), but also allows to study the kinematic signature of CI. The comparison of the simulated data and the observations will be shown in §12.

Chapter 7 H 2 Excitation

Molecular hydrogen, H 2 , is a homonuclear molecule formed by two atoms of hydrogen linked by a covalent bond, with a binding energy of 4.476 eV which corresponds to a temperature of 51942 K. Since this molecule is symmetric, it lacks of a permanent dipole moment. The excitation of H 2 can be an important diagnostic of the turbulent multiphase diffuse ISM. Although the goal of computing absorption and emission profiles of H 2 was not completely achieved in this thesis, we present in the next sections the levels and transitions of H 2 , and the code we built to compute H 2 excitation in post-processing.

H 2 Levels

Rotational Levels In first approximation, H 2 can be treated as a linear rigid rotor with a moment of inertia I m that takes into account the reduced mass µ of the molecule and the nuclear distance through I m = µr 2 0 , where r 0 is the inter-nuclei distance at rest. The energy levels are quantified and the rotational energy of a level i can be written as function of the moment of inertia of the molecule

E (J ) = J (J + 1)ħ 2 2I m = B J (J + 1) (7.1)
where J is the rotational quantum number, ħ is Planck constant devided by 2π, and B = ħ 2 2I m the rotational constant, which has a value of 85.25 K for H 2 . The distance between two neighboring rotational levels can be written as

∆E = 2B J (7.2)
with J the rotational quantum number of the upper level.

Ortho and Para Molecular hydrogen has two isomeric forms. Each hydrogen forming the molecule has a nuclear spin of 1/2. The two proton nuclear spins can be parallel forming ortho-hydrogen with total nuclear spin I = 1, or the two proton spins can be antiparallel forming para-hydrogen. Para-H 2 has an anti-symmetric nuclear wavefunction and since the total wavefunction must be anti-symmetric with respect to exchange of the two protons (from Pauli's principle), the rotational wavefuction needs to be symmetric with respect to permutation of the two protons, hence J must be even. Ortho-H 2 , on the other hand, has a symmetric nuclear wavefunction therefore J needs to be odd. The transitions ortho-H 2 ↔ para-H 2 are highly forbidden, because they need a change in the nuclear spins of the molecule. For these reasons, ortho-H 2 and para-H 2 are often considered as two independent species.

Consequently, only rotational transitions with even ∆J are possible and therefore the first rotational accessible excited level from the ground state is J = 2 at 510 K. The transitions between ortho-H 2 and para-H 2 can only occur through reactive collisions with H, H + , H 2 and H + 3 that are able to change the orientation of nuclear spins.

Vibrational Levels Vibrational levels are the result of the quantification of the different vibrational modes and their quantum number is v which corresponds to the number of nodes in the vibrational wavefunction. Diatomic and linear molecules can vibrate along their internuclear axis and the quantified energy of the vibrational levels can be written as

E (v) = hν 0 (v + 1 2 ) (7.3)
where ν 0 is the fundamental vibrational frequency.

Rovibrational Transition

One ro-vibrational level (v, J ) can be written as the sum of its rotational and vibrational levels:

E (v, J ) = hν(v + 1 2 ) + B J (J + 1) (7.4)
The transition between different levels are usually indicated with letters to indicate the jump in rotational levels as indicated in the table 7.1. Electronic Levels On the top of ro-vibrational states, there are the electronic states, that are more separated in energy.

Each electronic state of a molecule, is actually a superposition of a multitude of rovibrational levels. The different electronic transitions of H 2 are observed as broad bands which are the results of the superposition of the ro-vibrational lines.

The selection rules for molecular electronic transitions are ∆S = 0, ∆Λ = ±1, Σ + ↔ Σ + , Σ -↔ Σ -and g ⇄ u for homonuclear molecules. Therefore, the transition from the bound electronic state X to the unbound electronic state b is highly forbidden, while electronic transitions to Lyman and Werner bands are allowed. Photons with energy hν > 11.2eV, or equivalently wavelength of λ < 1108 Å, can excite molecular hydrogen from the ground state to the Lyman Band, while photons with ∆E > 12.3eV, or equivantely λ < 1008Å can excite the Werner band.

Rovibrational Excitation of H 2 in its Electronic Ground State

Where is the molecular hydrogen? Looking at simulations, with the help of the 2D probability histograms in the space of P /k B vs n H weighted by the mass of H 2 , we can characterize the physical conditions in which we find molecular hydrogen (see right panl Fig.

cient A i j in s -1 , which is proportional to the absorption oscillator strength f j i A i j = 8π 2 ν 2 i j e 2 m e c 3 g j g i f j i . (7.5)

As shown in 7.4, the Einstein coefficient of spontaneous emission for a pure rotational de-excitation with ∆v = 0 and ∆J = -2 increases with J i (e.g. for v = 0, A i j ∈ [2.9 × 10 -11 , 5.5 × 10 -6 ] s -1 ). The transition can also link rotational states with different vibrational values. The rovibrational transitions (∆v = 0) considered are quadripole transitions (O,Q,S). We do not consider induced emission and absorption (B i j and B j i ), because the occupation number of the photons in the wavelength of interest is low. As seen in the previous section, the rovibrational transitions are in the infrared waveband, and the ISM medium is not rich enough in such photons.

Collisional Transition

Even in the diffuse ISM, collisions can be an important source of excitation and de-excitation, especially in warm gas.

Since the collisional rates are proportional to the density of the partner and the crosssection is inversely proportional to the square root of the reduced mass, the most important collisional partners are the ones that are both abundant and light. We will, thus, consider collisions between H 2 and atomic hydrogen (Wrathmall & Flower 2007), helium (Flower et al. 1998), and molecular hydrogen (Flower 1998).

The rate of de-excitation per unit volume of a level due to a collision with the collisional parter C is written as

n i n C k C i j [cm -3 s -1 ] (7.6)
where n i is the density of the level i considered, n C is the density of the collisional partner,

k C
i j is the rate coefficient of the collision that depends on the collisional partner and the starting and ending level of H 2 . k C i j = 〈σv〉 i → j is the integrated cross-section, σ the crosssection, and v the relative velocity distribution.

The rate coefficients used in this work for collisional de-excitation result from a fit using three parameters a, b, c:

q(T ) = a + b t + c t 2 (7.7) k H i j (T ) = 10 q(T ) (7.8)
where t = T (K)/10 3 + δt , and δt = 1.0 to prevent the divergence of the rate coefficient when t goes to zero. i j indicate the transition form the state i = (v, J ) to a state at lower energy j = (v ′ , J ′ ). The excitation rate coefficients k j i is given by the detailed balance principle:

k j i = k i j g i g j e (-∆E i j /k B T ) cm 3 s -1 (7.9)
H 2 + H Rate coefficients for collisional de-excitations in reactive scattering due to collisions with H are taken from Le Bourlot et al. (1999):

k H r v J → v ′ J ′ = k H nr v J → v ′ J ′ × E J -J ′ even (7.11) k H r v J → v ′ J ′ = 1 2 k H nr v J → v ′ J ′ -1 + k H nr v J → v ′ J ′ + 1 × E J -J ′ odd and J even (7.12) k H r v J → v ′ J ′ = 1 6 k H nr v J → v ′ J ′ -1 + k H nr v J → v ′ J ′ + 1 × E J -J ′ odd and J odd (7.13)
where

E = exp -max 0, 3900 -E v J -E v ′ J ′ T (7.14)
The factor of 3900 K accounts for the activation energy barrier, k H r indicates the reactive collision rate with H, and k H nr non-reactive collision rate with H (see previous section).

Chemical Pumping

Until now we have considered excitation from collisions, both reactive or not, and deexcitation due to collisions and spontaneous emission. However, de-excitation and excitation can also occur during chemical reactions, i.e. at chemical formation and destruction of molecular hydrogen.

Assumptions

The formation of molecular hydrogen has an exothermicity of 4.4781 eV, i.e. the binding energy of H 2 . A common assumption is that one third of this energy is tranformed into internal energy of H 2 , one third into the grain on which H 2 forms and one third heats the medium (see §2). Under these assumption, the energy that internally excites H 2 is usually assumed to be spread among the levels following a Boltzmann distribution with a temperature T F . The probability of forming H 2 in a certain level i therefore reads:

p i = N g i exp(-E i /k B T F ), (7.15)
where N , is the normalization factor

N = 1/ i g i e (-E i /k B T F ) (7.16)
which basically corresponds to the inverse of the partition function.

If the total energy distributed in the internal rovibrational levels of H 2 is equal to 1/3 of the formation energy,

N lvl i =1 g i E i exp(-E i /k B T F ) N lvl i =1 g i e (-E i /k B T F )
= 4.4781/3 (7.17)

where E i is in eV and the Boltzmann constant k B = 8.617 × 10 -5 eV K -1 . We can easily see that T F depends on the number of level included N lvl . To avoid this problem, we choose a constant T F = 17323 K which correspond to one third of the formation energy in kelvins.

The downside of this approach is that including only a small number of levels in the computation, implies that the energy recovered is far smaller than 4.4781 eV/3 (e.g. if N lvl = 10, the recovered energy is 0.75 eV). Respectively, if the number of levels included is large, the energy distributed in the internal levels of H 2 is larger than 1/3 of 4.4781 eV.

Population as function of physical conditions

Determining the population of the different levels, under the steady-state assumption, means solving the equation:

n i C i = j n C k C i j (T ) + i < j A i j + k D = C i = j n C n j k C i j (T ) + i > j A j i n j + K i F (7.18)
for each level i of H 2 . n i and n j are the densities of the levels i and j , n C the density of a collisional parner and k C i j the rate coefficient to go from the level i to level j due to a collision with a partner C , A i j the spontaneous Einstein coefficient, k D the destruction rate of the level i (in s -1 ), K i F is the formation rate of the level i (in cm -3 s -1 ). On the left side of the equation we have all the processes de-populating the level i and on the right side the processes populating i . To solve this system of equations, we can write it in matrix form AX = B and look for the solution X = A -1 B . A is the matrix including all excitation and deexcitation rates

A = M rad + M col + M D ,
where M D is the destruction matrix, which contains the k D on the diagonal and zeros elsewhere. M rad is the radiation matrix containing the Einstein coefficients, and M col the collisional matrix containing C n C k C i j at the position i j . X is the population vector composed of the abundances of the levels x i , or the density vector composed of n i , and B is a vector composed of -p i k D . For the matrix to be inversible and the solution unique, this system needs to be closed by the conservation equation x i = 1, or n i = n(H 2 ). The solution X is then found with a Newton-Raphson algorithm.

Validation of the H 2 Excitation Code

In this section, we perform standard tests on the population of H 2 computed with the code.

At Local Thermal Equilibrium (LTE), the population of the levels should follow a

Boltzmann distribution

n i n(H 2 ) = g i e -E i /k B T i g i e -E i /k B T (7.19)
where g i and E i are the degeneracy and the energy of the level i and i n i = n(H 2 ).

2. In the ISM, the densities are generally lower than the critical densities, and the levels can be far from LTE, depending on the species considered and the physical conditions. We expect the levels with low critical densities, n i ,crit to be thermalized.

part of the chemical solver and be applied in post-processing to the ensemble of the simulations. This chapter presented the first step toward this more ambitious goal.

The main goal of this work is to compare the outcomes of the simulations with the observed samples introduced in Part I ( §4). In order to achieve this goal, we need to create a simulated sample of lines of sight statistically comparable to the observed sample. This is achieved by creating a set of synthetic lines of sight, patching together extractions from our MHD simulations ( §8). The comparison is optimally achieved with a modified KS test that optimizes the comparison of our observed and simulated samples, on a multivariate statistical basis ( §9). These comparisons are then analyzed in the subsequent chapters of this last part of the manuscript.

Chapter 8

Reconstruction Algorithm

This chapter presents and discusses how we have constructed sets of simulated lines of sight through neutral diffuse ISM, using the outputs of the parametric study presented in §5, with the goal to perform reasonably accurate comparisons with the observed samples presented in §4. While we had adopted a fiducial size of 200 pc for the MHD simulated lumps of diffuse neutral gas, we now have to take two important facts into account: the distribution of the lengths of the lines of sight through the ISM in the observed sample, and the heterogeneity of the medium along these lines of sight, as only a fraction of the lines of sight is expected to be neutral medium.

We first describe the method used to create our sets of simulated lines of sight (reconstruction). We then apply this method to two of our observed column density samples: the atomic vs molecular hydrogen sample (the HI-to-H 2 transition, see §8.2), and the CH + molecular ion vs total hydrogen sample (see §8.3). In the last part of this chapter (see §8.4), we discuss some of the limitations of the method adopted. The scientific analysis and the conclusions will be drawn later in §10, §11, and §12.

Method

Let's consider an observational sample. The first step is to extract the distribution of lengths of the line of sight, l los . Since the observational sample may contain extragalactic sources or galactic sources at high latitude, and since the amount of molecular gas in the Milky Way decreases exponentially as a function of the distance from the midplane, the length l los occupied by the observed diffuse gas cannot always be identified to the distance of the background source. For the sake of simplicity, adopting a scale height of 100 pc above the galactic midplane for the cold ISM component, we estimate the length of one line of sight through the ISM as

l los = min 1" p , 100 sin(|b|) pc, (8.1)
where b is the Galactic latitude of the background source and p is its parallax.

The distribution of these estimated lengths for a given sample is built as a histogram regularly sampled on a logarithmic length scale. A simulated sample is then built with

Tests and Caveats

The construction of synthetic lines of sight, explained in §8.1, is based on the concatenation of randomly selected cuts through a single simulation output cube, assuming a uniform filling factor ϕ for the ionized gas in the local interstellar medium. The consequences of these simplifying assumptions are successively examined in the present section through the effects of changing the random concatenation process, patching cuts from different simulations or from independent outputs of a single simulation, and varying the value of the assumed filling factor.

The Random Drawing Process

Reconstruction by Translation of the Box We introduce here an alternative method to generate lines of sight with the same distribution of lengths as the observational sample. When the column densities are generated as explained previously, the random process involved decorrelates the spatial structure (see the left panel in Fig. 8.7).

Instead of drawing independent segments for each of the reconstructed lines of sight, a long patch can be created once with the concatenation of random translations of the same cube. All lines of sight can then be drawn from this long patch. Indeed, since the boundary conditions of the simulation are periodic, the simulated box repeats itself every period L in the x, y, z directions and each point of the box (X,Y) can be seen as the origin of a new translated box. Fig. 8.6 shows a schematic view on how the box-translation is built. A point (X, Y) is drawn randomly in the simulation box and divides the cube into 4 parallelepipeds (A,B,C,D). These parallelepipeds are rearranged to create a new simulation box (C,D,A,B) which is a simple translation of the original box. In the box-translation reconstruction, the column densities of the generated lines of sight are computed as an integration over boxes with different origin points. Following the explanation in the precedent section, if (1-ϕ)l los > L, we create L/[(1-ϕ)l los ] translated boxes by drawing random origin points (see Fig. 8.6) and then we add the respective individual column densities to create the column densities of lengths l los , therefore, the building blocks of the generated lines of sight are boxes translated by a random distance the one compared to the other.

In the right panel of Fig. 8.7, we can see how the reconstruction with the box-translation algorithm keeps information on the structure of the matter. This alternative method should be considered for studies involving the spatial structure.

Comparison Translation and Random Reconstruction

In Fig. 8.7, we show the reconstructed total column density of length l los = 800 pc for the fiducial simulation. On the left with the random method explained in the previous section, and on the right with the boxtranslation method. We notice that they both extend on the same range of column densities, but the random one does not contain any information on the original structure. In Fig. 8.8, we show the isocontours of the 2D probability distributions found with both methods. The two methods give the same statistical results on the column densities, allowing us to compute the probability histograms, necessary to our study, with the random method, which has the advantage of being simpler. Despite the fact that this reconstruc-tion method has no use in cases where the information on the structure is key, it is a very useful method in our work, because compared to the box-translation method it has another perk. The different lines of sight to the observed sources may start from the same local ISM patch, but further from the observer, the patches traversed become less and less correlated. Hence drawing the successive cuts independently from each other is better than drawing all lines of sight from a single concatenation/realization. For these reasons, the random reconstruction seems a better method to model the observed lines of sight in which no correlation between the different neutral clouds is assumed/observed.

The Choice of Concatenated Patches

Different Simulations Yet a better way of generating lines of sight of length l los would be to take different simulation runs with the same set of parameters and stock them one behind the other to reconstruct our lines of sight. This procedure would allow to have a coherent structure in each piece of neutral medium and conserve structure information in the column densities map, without assuming any spatial correlation between contiguous pieces of the neutral medium along each line of sight. This method has a high computational cost because the number of simulations required (that need to be run and post-processed) to perform the parametric study would be at least 8 times larger than the current one. To address this issue, but still be able to analyze a possible result, we used different outputs of the same simulation. We chose two outputs separated by one correlation time because after a correlation time the initial conditions are erased. Since the correlation time of the turbulent forcing was set to the turbulent turnover time in our simulations, the time difference between the two outputs needs to be larger than τ l = l /σ l where l can be any scale. Since we want the turnover time of the driving scale, l = L drive , the velocity dispersion at L drive = 100 pc is of σ 100pc ∼ 6 km s -1 , giving a correlation time of around 16 Myr.

Effect of Different Outputs on the Random Reconstruction

To see the effect of this approach, we analyze the random reconstruction done on two different outputs. In the reconstruction explained in the previous paragraph, each output would model a piece of neutral interstellar medium along lines of sight, therefore there would be a contribution of different outputs along one line of sight. Here we chose to apply our original method to the two different, assumingly independent, outputs, in order to assess the sensitivity of the simulated histograms to the choice of the simulation output.

In the left panel of Fig. 8.9, we can see that the total column density maps integrated over L = 200 pc essentially show no correlation in the structure of the matter. The right panel shows both histograms for the reconstructed lines of sight using the later output as a building block for the colored probability histogram and the solid lines, and the early output for the dashed lines. We notice that the difference between the dashed and solid curves is insignificant, validating the random reconstruction approach. We remind the reader that the more complex method would consist in taking not one building block (as done here), but L/[(1 -ϕ)l los ] different building blocks, each one corresponding to a different output, thus any difference visible on the right panel of Fig. 8.9 would be diminished. Hill et al. 2018) on the volume of the HIM. These latest studies estimate that 20% ϕ 90% in the Galactic midplane.

Highly uncertain, the volume filling factor of the HIM can also vary from one line of sight to another (Fig. 1 of Hill et al. 2018). In Fig. 8.10, we tested the results of the reconstruction algorithm assuming different values of ϕ. Increasing ϕ has the effect of stretching the predicted 2D probability histogram toward low values of N H along the x-axis, and gives a wider HI-to-H 2 transition, that can extend on almost 3 orders of magnitude from N H ∼ 10 18 cm -2 to N H ∼ 10 21 cm -2 in the extreme case in which most of the line of sight is composed of the hot and warm ionized medium. The balance of probabilities of occurrence of lines of sight at high and low molecular fraction is also modified, and the occurrence of low column densities increases as a function of ϕ.

Because of the lack of additional information and since no major difference is seen in the probability histogram distributions for value of ϕ given in the litterature, i.e. 0.25 < ϕ < 0.75, we assume here a conservative value of ϕ = 0.5. Taking into account a realistic distribution of ϕ would require to simulate the Galactic disk and halo over a scale of several kiloparsecs and to properly model and follow the impact of supernovae explosions. This is far beyond the scope of the present work.

One Realization We note that, while spatially uncorrelated, the pieces of diffuse neutral gas used in the reconstruction algorithm correspond to random realizations obtained with a single simulation. The approach of assuming that the local ISM can be built out of a single simulation was chosen in order to highlight the natural variations induced by turbulence and thermal instability alone in a diffuse neutral gas with a known averaged density and UV illumination factor. However, because the medium probed by the observations extends in all directions around the sun over a maximum distance of 3 kpc (see Fig. 8.2), it stands to reason that potential variations of all parameters should be taken into account, not only from one line of sight to another but also along a single and outstretched line of sight. We ignore this aspect here because it would drastically increase the parameter space to explore. However, such consideration needs to be kept in mind for the comparison with the observations. Chapter 9

Kolmogorov-Smirnov Statistics

As shown in the previous chapter, the originality of the work presented here is to compare simultaneously the chemical and statistical information of observed and simulated samples by the mean of probability histograms (PHs), or equivalently, probability distribution functions (PDF). In many ways, the direct comparison of probability histograms, such as those shown in Figs. 8.8 and 8.10, provides valuable information regarding the statistical nature of the correlation between different species. Unfortunately, such an approach is impractical to use for the comparison of an observational sample to a large set of simulations (obtained with a parametric study) or if the chemical correlations are studied in dimensions larger than two.

For all these reasons, we developed a procedure to measure with a single value, the degree of similarity between two probability distribution functions. This procedure, based on the well-known Kolmogorov-Smirnov (KS) test, can be used in any dimension and can take into account all observational information including actual detections and lower or upper limits. In the following sections we will introduce the classic Kolmogorov-Smirnov test ( §9.1), its extension to N-dimensional samples, the modification, and definition of the KS distance, R KS , and its application on two observational samples.

Classic KS test

The KS test is a cumulative test originally built to compare two PDFs in one dimension. In a nutshell, the KS test consists in computing the KS difference which measures the maximum value of the absolute distance between two distribution functions (two-sample KS test), (9.1) or between one sample and a model (one-sample KS test),

D = max -∞<x<∞ S N 1 (x) -S N 2 (x)
D = max -∞<x<∞ |S N (x) -P (x)| . (9.2)
S N (x), S N 1 (x) and S N 2 (x) are cumulative distribution functions, i.e. functions giving the fraction of datapoints smaller than x, and P (x) is a known cumulative probability distribution.

quadrant. The observational dataset O used to compute the KS distance is chosen as the subsample such that all quadrants scanned contain at least n thr observed lines of sight. With this assumption, the error on the merit function is calculated by taking into account only the statistical errors on the number of simulated lines of sight. For each quadrant, we assume that the "true" merit function ranges between

M minus = log 10 f sim -3 ( f sim /S)
f obs and (9.5)

M plus = log 10 f sim + 3 ( f sim /S) f obs , (9.6)
where S is the total number of simulated lines of sight. Because the errors are asymmetric,

M plus or M minus can tend toward infinity. If so, the KS distance is a lower limit, even if the infinite error bar is obtained for another quadrant than the one that maximizes M . In short, for each observational datapoint, we compute M minus , M , M plus , and R KS . If one of the M minus goes to infinity, R KS is considered as a lower limit. The lower end of the errorbar is given by the minimum value between M minus and M plus evaluated in the quadrant that maximizes the KS difference, R KS .

M min = min max(M ) (M minus , M plus ), (9.7
)

M max =    max O (M minus , M plus ) if either → +∞ max max(M ) (M minus , M plus ) if both finite (9.8)
The definition adopted above, not only allows to compute the merit function value but also to compute the errors on the merit function. We should bear in mind that the number of scanned points during the KS procedure, and the errors on the merit function depend on the minimal number of observations per quadrant, n thr . We will show the impact of n thr in the following section.

Influence of the Distribution of Datapoints

Another important aspect for the stability of the KS test is linked to the shape of the observational dataset. Since the cumulative probability is defined as the fraction of points per quadrant and since the division of quadrants is performed on a cartesian grid, the result of the KS test, explained in the previous section, depends on the choice of the axes.

To explain this concept, let's imagine we observe two quantities X and Y , which are strongly correlated, as in the left panel of Fig. 9.3. We can see a real correlation between X and Y and if we impose n thr = 1 datapoint, we scan only 15% of the observed datapoints. If we perform a rotation on the observational sample (see right panel of Fig. 9.3), and thus we plot Y -aX , where a is the slope of the correlation between X and Y , we see that the distribution is now flat and the percentage of the scanned point increases greatly (to 90% in the example) and therefore the subsample over which the KS distance is computed is more populated.

For any observed quantities, X and Y , the best mathematical option would be to identify the principal components of the observational sample using proper orthogonal decomposition or singular value decomposition algorithms. Such a method could even be to observationsor try to extract some information from the non-detections. To do so, and integrate the non-detections in the KS distance, we assume that a non-detection is equivalent to an upper limit, Y upl , i.e. it signals a detection at an unknown column density Y < Y upl and X = X upl .

The subsample O of the scanned datapoints assuring a minimal fraction of observations per quadrant does not change from the previous version, i.e. only detections give receivable couples of values (X obs , Y obs ). The difference between the previous version of the KS test resides in the fraction of observed datapoints per quadrant, which now can include also upper-limits. The contribution to the fraction of observations, f obs , of each non-detection depends on its value Y upl . For all non-detections with Y upl < Y obs , we can assume that the "real" detection belongs to the corresponding quadrant and, therefore, contribute to the fraction of observations in a quadrant for 1/n(O ). For non-detections in the upper quadrants (B and D), i.e. Y upl > Y obs , nothing sure can be deduced to such a value and it is therefore excluded, i.e. it does not count in any quadrant. In Fig. 9.10, we give an example of the method. The red points are the detections that are scanned imposing a threshold on the number of observationsboth detections and non-detectionsin each quadrant. When a detection point is scanned (X obs , Y obs ), it divides the space into four quadrants (dark red). Non-detections are indicated with triangles, the upper limits with Y upl < Y obs are counted in quadrants A or D as detections (light green triangles), the upper limits with Y upl > Y obs are excluded (gray triangles). This is one way to tackle the problem, in §9.5.2 we discuss different assumptions and subsequent methods to count for the contribution of an upper-limit in a quadrant.

Application on CH +

We applied the modified KS test, including the non-detections as explained in the previous section and we show the result of this modified KS test in Fig. 9.11. The orientation of the axis is N (CH + ) vs N H and n thr = 6, as before. Taking into account also information from the upper-limits the KS distance has improved (cf Fig. 9.9). While from the previous KS test on detections, we could only conclude that the main difference between the observational dataset and the simulated one was due to the non-detections (see Fig. 9.9), in the new KS test case we can extract more valuable information. The upper-limits are located in the same range of N H as the low-N (CH + ) bump, so the A quadrant is not the quadrant that maximizes the merit function. The quadrant that maximizes the merit function is now quadrant D (indicated with a gray rectangle in 9.11), where the fraction of simulated lines of sight is around three times larger than the fraction of observed lines of sight. This indicates that either the dispersion of the simulated lines of sight in the high-N (CH + ) bump is larger than the observed one, or the simulated low-N (CH + ) bump extends toward unrealistically large N H , which could suggest that the transition between the two simulated bumps happens at sightly too larger total column densities.

In general, the KS distance is smaller when the upper limits are taken into account. Not all distances computed including the upper limits show as much variation as the fiducial case. In §11 we show the KS differences computed with and without the non-detections.

When including the upper-limits in the KS test applied to CH + we obtain a value of 0.47 for the KS distance. This means that in the worst quadrant the fractions of simulations and observations differ by a factor of less than three.

Different Assumptions to count the contribution in a quadrant of a non-detection

The KS test including the non-detections is a great improvement on the method. Here we proposed two additional more convoluted methods.

Half Method An attempt to take into account all upper-limits, also the one with Y upl > Y obs , i.e. belonging to B or C quadrants (see 9.12), would be to assume that the corresponding detection has the same probability of being larger or smaller column densities than Y obs . This translates into giving a weight of one half in both A and B (or C and D) quadrants to all detections that belong to B (or C). Taking the same example as before, we show the concept in Fig. 9.12. Detections (red points) in each quadrants count as 1/n(O ) in f obs , upper-limits with Y upl < Y obs counts as detections, i.e. 1/n(O ) (light green triangles), and upperlimits with Y upl > Y obs counts as half in both the upper and lower quadrant (A-B or C-D depending if X upl < X obs or X upl > X obs ), i.e. a contribution of 0.5/n(O ) in each (dark green points).

Distance Method

We could imagine that a non-detected species has the same probability, in logarithmic space, to have whatever value between Y min up to Y upl , where Y min is an arbitrary column density and can be set to whatever value between zero and the minimal column density of the observed species.

With these assumptions, we can define a contribution to the fraction of observations due to a non-detection in each quadrant, defined by the datapoint (X obs , Y obs ). Each Y upl counts as a fraction of 1/n(O ) in the quadrant A and B (or C and D) if X upl ≤ X obs (or X upl > X obs ) that is normalized to 1/n(O ). The contribution to the fraction of observations, f obs , of each non-detection, Y upl , depends on the distance of the non-detection in question to the origin of the quadrant, Y upl -Y obs . Fig. 9.13 shows that if Y upl > Y obs the non-detection contributes to the fraction of observations of both the upper quadrant (B or C) and the lower quadrant (A or D), and this contribution is proportional to Y upl -Y obs and Y obs -Y min , respectively. Normalizing, we obtain the contribution of each upper-limit in the case

Y upl > Y obs 1 n(O ) Y upl -Y obs
Y upl -Y min for B and C quadrants, (9.9)

1 n(O ) 1 - Y upl -Y obs
Y upl -Y min for A and D quadrants, (9.10) while if Y upl < Y obs the contribution of that upper-limit to the fraction of observations in a quadrant counts as 1 n(O ) in the quadrants A or D, depending if X upl < X obs or not.

Comparisons on CH + In Fig. 9.14, we apply the "half" (left panel) and the "distance" (right panel) method to the 2D probability histogram of N (CH + ) vs N H for the fiducial sim-

Chapter 10

Atomic-to-Molecular Transition

The new and original analysis tools presented in the previous chapters have been applied to several sets of observations (see §4) in combination with a parametric study of the diffuse, neutral, and multiphasic ISM. The first application dedicated to the analysis of the HI-to-H 2 transition has been published in Bellomi et al. (2020). In a nutshell, the grid of simulations run with RAMSES (see §5) were used to compute the distributions of the total column densities, N H and molecular hydrogen column densities, N (H 2 ) taking into account the distribution of the observational sample (see §4.1) through the reconstruction method (see §8.1). The Kolmogorov-Smirnov difference (see §9.2) was then computed to facilitate the comparison of the simulations with the observations. The results are presented and discussed in the paper which is included below. In the next section, we simply present two additional results. The first one, discussed but not shown in the paper, focuses on the comparison of the predictions of simulations assuming two different prescriptions for the self-shielding of H 2 . The second section shows the impact of including the upperlimits in the statistical analysis of the HI-to-H 2 transitions.

H 2 Self-shielding at High Temperature

The prescription of H 2 self-shielding used in the paper (Eq. 5.16) is taken from Draine & Bertoldi (1996). As discussed by Wolcott-Green et al. (2011), such a prescription is reliable for diffuse gas at low temperature but becomes less and less reliable for high-temperature environments (T > 500 K) where efficient collisional excitation of H 2 in its rovibrational levels reduces the self-shielding. To estimate the impact of this process, we ran the fiducial simulation with the alternative self-shielding function proposed by Wolcott-Green et al.

(2011):

f shield (x) = 0.965 (1 + x/b D ) α + 0.035 e -8.5×10 -4 1+x 1 + x , ( 10.1)
where α = 1.1 and x = N (H 2 )/5 × 10 14 cm -2 (10.2) and b D is the Doppler broadening parameter expressed in km s -1 . In Fig. 10.1, we compare the 2D probability histograms of the HI-to-H 2 transition with the self-shielding computed

Introduction

The multiphase nature of the interstellar medium (ISM) is at the root of the regulation of star formation in galaxies (e.g., Hill et al. 2018). As shown by the emission profiles of the HI 21 cm line (Heiles & Troland 2003a,b;Murray et al. 2015Murray et al. , 2018)), the diffuse neutral ISM is composed of two stable thermal states at thermal pressure equilibrium (Jenkins & Tripp 2011), the warm neutral medium (WNM, T ∼ 7000 K) and the cold neutral medium (CNM, T ∼ 70 K), coexisting with a third unstable state, the lukewarm neutral medium (LNM), whose temperature is comprised between those of the CNM and the WNM (e.g., Marchal et al. 2019). Through condensation and evaporation processes, turbulent transport, and turbulent mixing, the diffuse matter flows from one stable state to the other eventually leading to the formation of dense and cold clouds massive enough to trigger gravitational collapse (e.g., Ostriker et al. 2010). While this picture is widely accepted, the intricated effects of turbulence, gravity, radiation field, and magnetic field on the exchange of mass and energy between the different phases and on the formation of structures at all scales has yet to be unveiled. Following the illustrious analytical descriptions of the thermal instability process (Field 1965;Wolfire et al. 1995Wolfire et al. , 2003;;Bialy & Sternberg 2019), several analytical and numerical studies have been dedicated to understand the dynamical evolution of the gas, focusing on the formation of CNM structures, molecular clouds, and collapsing cores (e.g., Hennebelle To extend the predictions of simulations to a larger set of observational diagnostics, recent numerical studies have undertaken the challenging task of solving the chemical evolution of turbulent and/or multiphase environments. Originally dedicated to the formation of CO in molecular clouds and to the analysis of the CO-to-H 2 conversion factor in galaxies and the CO dark-gas (e.g., Glover et al. 2010;Smith et al. 2014;Richings & Schaye 2016a;Seifried et al. 2017;Gong et al. 2018), numerical simulations are now used to study a variety of atomic and molecular tracers, including CII, CI, CH + , OH + , H 2 O + , and ArH + (e.g., Richings & Schaye 2016b;Valdivia et al. 2017;Clark et al. 2019;Bialy et al. 2019). All these works demonstrate the predictive power of astrochemistry. The column density distribution of each atom and molecule has a unique signature that provides detailed information on the thermodynamical state of the diffuse matter (Clark et al. 2019). In turn, the confrontation with the predictions of numerical simulations can be used to estimate the scale and strength of the injection of mechanical energy by stellar feedback (Bialy et al. 2019), the large-scale turbulent transport and the interfaces between CNM and WNM (Valdivia et al. 2017), and the nature of the turbulent dissipation processes (Lesaffre et al. 2020).

In this context, understanding the formation and survival of molecular hydrogen has long been recognized as a major topic of investigation. As the most abundant molecule in space, H 2 is at the root of interstellar chemistry and the growth of molecular complexity. In addition, and because its formation preferentially occurs in dense environments, H 2 naturally correlates with the star formation rate of galaxies (e.g., Lupi et al. 2017) and therefore offers a valuable proxy to understand the limit in the Kennicutt-Schmidt relation above which star formation occurs (e.g., Bigiel et al. 2008Bigiel et al. , 2011;;Schruba et al. 2011;Leroy et al. 2013).

Over the last decades, great efforts have thus been devoted to propose analytical descriptions of the HI-to-H 2 transition in homogeneous clouds with plane-parallel or spherical geometries (e.g., Sternberg 1988;Krumholz et al. 2008;McKee & Krumholz 2010;Sternberg et al. 2014 and references therein), compute this transition in detailed 1D chemical models assuming chemical equilibrium (e.g., van Dishoeck & Black 1986;Abgrall et al. 1992;Le Bourlot et al. 2012;Bron et al. 2014) or not (e.g., Lee et al. 1996;Goldsmith et al. 2007;Lesaffre et al. 2007), treat the chemistry of H and H 2 in subgrid models applied to simulations of galaxy formations (e.g., Gnedin et al. 2009;Christensen et al. 2012;Thompson et al. 2014;Diemer et al. 2018), or solve it in 3D isothermal or multiphase simulations of the diffuse ISM using various treatments of the radiative transfer (e.g., Glover et al. 2010;Valdivia et al. 2016;Hu et al. 2016;Bialy et al. 2017;Nickerson et al. 2018).

Thanks to all these works, a global picture of the formation of H 2 in galaxies is now emerging. At the scale of a homogeneous cloud, the molecular content, the sharpness of the HI-to-H 2 transition, and the asymptotic column density of HI are controlled by the ratio of the intensity of the ultraviolet (UV) field to the gas density and the dust-to-gas ratio, or equivalently, the metallicity (Sternberg et al. 2014). At larger scales, the integrated column densities of HI and H 2 also depend on the distribution of clouds of various densities along the line of sight and on the porosity to the UV radiation field. Because of these effects, the statistical properties of the total column density are found to depend on the strength, the scale, and the compressibility of the turbulent forcing in simulations of CNM gas (Micic et al. 2012;Bialy et al. 2017). The amount of molecular gas depends on the "clumpiness factor" used for the subgrid models in simulations of galaxy formation (Gnedin et al. 2009;Christensen et al. 2012).

Despite these achievements, very few works have been dedicated so far to the analysis of the HI-to-H 2 transition in a turbulent multiphase medium at a scale sufficient to resolve the formation of CNM structures. In addition and while the predictions of analytical models (e.g., Krumholz et al. 2008) and simulations (Gnedin et al. 2009;Valdivia et al. 2016) were able to reproduce the trend of the HI-to-H 2 transition observed by Copernicus and FUSE in the local ISM (e.g., Savage et al. 1977;Gillmon et al. 2006;Rachford et al. 2009), the LMC and the SMC (e.g., Browning et al. 2003;Gillmon et al. 2006;Leroy et al. 2007), no detailed comparison with the statistical properties of these observations have been proposed. As a result, the occurrence of lines of sight with large molecular fractions predicted by numerical simulations often exceed what is deduced from the observations (Valdivia et al. 2016). Finally, and while statistical studies of 1D probability distribution functions (PDF) have become a common tool to understand the formation and the dynamics of molecular clouds (Körtgen et al. 2019), few statistical studies have been performed to date on 2D probability distribution functions using combined observations of different molecular tracers. In that perspective, the recent work of Bialy et al. (2019) opens new horizons for the analysis of chemistry in the diffuse matter.

In the first paper of this series, we extend these pioneer statistical studies to the measurements of the atomic-to-molecular transition observed in the diffuse and translucent ISM located in a radius of ∼3 kpc around the sun. We perform a parametric exploration of numerical simulations of the multiphase ISM and compare the results with the observed 2D probability histogram (PH) of total and molecular hydrogen column densities in order to identify the physical processes that control the molecular content of CNM clouds and the probability of occurrence of lines of sight. The observational dataset and the distribution of sizes of the sampled medium are presented in Sect. 2. The different setups of the simulations and the method used to reconstruct the 2D PH are described in Sect. 3. The comparisons with the observations are shown in Sect. 4 which also highlights the influences of the different parameters. The paper finally ends with Sects. 5 and 6 where we discuss the validity of our approach and summarize our main conclusions.

Observations of the HI-to-H 2 transition

Observational sample and distances

The observational sample studied in this work is built from the database of Gudennavar et al. (2012) A.1). The blue dashed line indicates the maximum value of N(H 2 ) derived from a purely molecular medium with an integrated molecular fraction f H 2 = 1 (Eq. ( 2)). The red dashed-dotted line indicates the theoretical molecular fraction derived in an unshielded WNM-type environment with a density of 0.5 cm -3 and a temperature of 8000 K, illuminated by a UV photon flux of 10 8 cm -2 s -1 (see Eqs. ( 13) and ( 15)). The regions A, B, C, D, and E defined in Table 1 correspond to an arbitrary separation of the observational sample used for quantitative comparisons with the results of simulations (see Sect. 4). 65% to subsample D. Interestingly, the subsample E is empty and no line of sight is observed with N H > 10 22 cm -2 . While the mean value of the logarithm of the molecular fraction strongly increases from subsamples A to D, the dispersion simultaneously decreases by about a factor of three, probably revealing an effect of average over long distances.

All these statistical properties, and more precisely the probability of occurrence of a given line of sight, are the subject of this paper. What physical processes control the HI-to-H 2 transition? How does the distribution of lengths of the lines of sight influence its observed statistical properties? What are the origins of the variations of the molecular fraction from one line of sight to another?

Physics and numerical method

To study the physical processes at play in the HI-to-H 2 transition, we performed numerical simulations of the multiphase diffuse ISM, using the RAMSES code (Teyssier 2002;Fromang et al. 2006), a grid-based solver with adaptative mesh refinement (Berger & Oliger 1984). The methodology applied in this paper follows the works of Seifried et al. (2011), Saury et al. (2014), and Valdivia et al. (2015, 2016).

The diffuse matter in the Solar Neighborhood of our galaxy is simulated over a box of size L with periodic boundary conditions. The matter, defined by a mean proton density n H , is assumed to be illuminated on all sides by an isotropic spectrum of UV photons set to the standard interstellar radiation field (Habing 1968) and scaled with a factor G 0 .

Fluid equations

Within this framework, RAMSES computes the evolution of the gas solving the classic equations of ideal magnetohydrodynamics (MHD):

∂ρ ∂t + ∇ • (ρu) = 0, ( 3
)
∂ρu ∂t + ∇ • (ρuu -BB) + ∇P = -ρ∇Φ + ρ f , (4)
∂E ∂t + ∇ • [(E + P)u -B(Bu)] = -ρu • ∇Φ + ρ f • u -ρL, and (5)
∂B ∂t + ∇ • (uB -Bu) = 0, (6)
where ρ, u, B, P and E are the mass density, the velocity field, the magnetic field, the total pressure, and the total energy density, respectively. The net cooling function per unit mass, L, and the acceleration due to the turbulent driving, f , are described in Sects. 3.2 and 3.3. The axis x, y, and z are chosen so that z corresponds to the direction perpendicular to the Galactic disk, and x corresponds to the direction of the mean magnetic field initially parametrized by a constant value B x .

To take into account all gravitational forces, including selfgravity and the action of stars and dark matter, the gravitational potential Φ is divided into two terms:

Φ = φ gas + φ gal . ( 7
)
The self-gravity potential, φ gas , is deduced from the Poisson's equation:

∇ 2 φ = 4πGρ. ( 8
)
Following Kuijken & Gilmore (1989) and Joung & Mac Low (2006), we assume that the Galactic potential along the direction z perpendicular to the Galactic disk can be written as

φ gal (z) = a 1 z 2 + z 2 0 -z 0 + 2a 2 z 2 , ( 9
)
where the first term is the contribution of the stellar disk parametrized by z 0 = 0.18 kpc and a 1 = 1.42 × 10 -3 kpc Myr -2 and the second term is the contribution of the spherical dark halo parametrized by a 2 = 5.49 × 10 -4 Myr -2 .

A36, page 4 of 33 E. Bellomi et al.: 3D chemical structure of diffuse turbulent ISM. I.

Thermal processes and radiative transfer

As shown by Field (1965), the multiphase nature of the ISM results from the thermal balance of the gas and thus from its net cooling function L defined by

ρL = n 2 H Λ -n H Γ, (10)
where n 2 H Λ and n H Γ are the cooling and heating rates of the medium (in erg cm -3 s -1 ) and n H is the proton density. To correctly describe the thermal state of the diffuse ISM, we include in this work the heating induced by the photoelectric effect and the decay of cosmic ray particles and the cooling induced by Lyman α photons, the recombination of electrons onto grains, and the fine structure lines of OI and CII. All these processes, described in Appendix B, are modeled with the analytical formulae given by Wolfire et al. (2003).

The absorption of UV photons by dust, and its subsequent impact on the photoelectric effect, is treated with the treebased method proposed by Valdivia & Hennebelle (2014). At each point the effective radiation field G eff (in Habing units) is computed as

G eff = G 0 e -2.5A V , ( 11
)
where A V is the visual extinction along a given ray, deduced from the integrated proton column density 1 computed from the border of the box to the current point

A V = 5.34 × 10 -22 N H cm -2 , ( 12
)
and e -2.5A V is an average performed over 12 directions, treated as solid angles evenly spread in polar coordinates.

Turbulence forcing

To mimic the injection of mechanical energy in the diffuse ISM, a large scale turbulent forcing is applied. Following Schmidt et al. (2009) and Federrath et al. (2010), this forcing, modeled by an acceleration f in the momentum conservation equation, is driven through an Ornstein-Uhlenbeck process using a pseudo-spectral method. At regular time intervals ∆τ, random fluctuations of the forcing term are generated and applied over an autocorrelation timescale τ. To excite only large scale modes, the forcing is modeled as a paraboloid in Fourier space covering a small interval of wavenumbers 1 k 3 and centered on k = 2. Using the notations of Seifried et al. (2011) and Saury et al. (2014), the total magnitude of these perturbations is set with either an acceleration parameter F or, equivalently, a velocity parameter V related by F = V 2 /L drive , where L drive is the main driving scale, L drive = L/2. A Helmholtz decomposition is finally applied, in order to control the powers injected in compressive and solenoidal modes. Using the classical notation, these powers are set with a parameter 2 ζ ranging from a pure solenoidal field (ζ = 1) to a pure compressive field (ζ = 0). Throughout this work, we adopt ∆τ ∼ 0.4 Myr which roughly corresponds to the time interval separating two supernova events occurring in a volume of (200 pc) 3 . The characteristic damping time of the turbulence τ is approximately set to the turnover 1 In this work we use the relation between N H and A V deduced from the observations of the mean Galactic extinction curve (Fitzpatrick & Massa 1986). 2 We note that for ζ = 0.5, the power of the compressive forcing corresponds to 1/3 of the total power. timescale of the diffuse ISM τ ∼ 33(L/200pc) 0.6 Myr (Larson 1981;Hennebelle & Falgarone 2012). F (or V) and ζ are left as free parameters.

H 2 chemistry

The timescale required for the abundance of molecular hydrogen to reach its equilibrium value is known to range over several orders of magnitude, depending on the physical conditions of the ISM (e.g., Goldsmith et al. 2007;Tabone 2018). In the diffuse gas, this timescale varies typically between a few 10 3 yr and a few 10 7 yr (Valdivia et al. 2016), hence over a range of values that often exceeds the dynamical timescales. To take into account this important aspect of the diffuse interstellar chemistry, the outof-equilibrium abundance of H 2 is computed self-consistently in the simulation, using the formalism introduced in RAMSES by Valdivia et al. (2015Valdivia et al. ( , 2016)).

The formation of H 2 onto grains in physisorption and chemisorption sites is modeled with the simplified rate of Le Bourlot et al. (2012)

k f = 3 × 10 -17 n H n(H) T 100 K S (T ) cm -3 s -1 , ( 13
)
where n H and n(H) are the local proton and atomic hydrogen densities, and

S (T ) = 1 1 + T T 2 β (14)
is the sticking coefficient of H onto grain, parametrized by T 2 = 464 K and β = 1.5.

The destruction of H 2 by UV photons is computed using the formalism, described by Draine & Bertoldi (1996) and Sternberg et al. (2014), which is classically introduced in many astrochemical models (e.g., Lesaffre et al. 2013;[START_REF] Bellomi | 3D chemical structure of diffuse turbulent ISM[END_REF]Bialy et al. 2017). In each cell, the photodestruction rate of H 2 is modeled as

k d = k 0 d n(H 2 )G 0 e -σ d N H f shield (x) cm -3 s -1 , (15)
where x = N(H 2 )/5 × 10 14 cm -2 , k d,0 = 3.3 × 10 -11 s -1 is the inverse freespace dissociation timescale of H 2 in an isotropic Habing field, n(H 2 ) is the local density of the molecular hydrogen, e -σ d N H is the shielding induced by dust and f shield the selfshielding function. We adopt here an effective dust attenuation cross section at λ = 1000 Å, σ d = 2 × 10 -21 cm 2 (Sternberg et al. 2014). Following Draine & Bertoldi (1996), the self-shielding function is computed as

f shield (x) = 0.965 (1 + x/b D ) 2 + 0.035 e -8.5×10 -4 √ 1+x √ 1 + x , ( 16
)
where b D is the Doppler broadening parameter expressed in km s -1 . As done for the photoelectric heating rate (see Sect. 3.2), both the shielding by dust and the self-shielding are calculated along 12 different directions and then averaged to obtain the photodissociation rate of Eq. ( 15).

Fiducial model and grids of parameters

The framework described above lays on several independent parameters, which are all related to key physical ingredients of the ISM. The influence of each ingredient on the HI-to-H 2 A36, page 5 of 33

A&A 643, A36 (2020) With our assumptions, L simultaneously corresponds to the scale of illumination of the gas by UV photons and twice the integral scale of turbulence, that is twice the scale of injection of mechanical energy L drive . OB stars, which are the dominant sources of the interstellar UV field, are not uniformly distributed in the sky but are known to be clustered in associations (Ambartsumian 1947). As shown by the recent 3D studies of the distributions of stars based on the HIPPARCOS and Gaia Catalogs (e.g., Bouy & Alves 2015, Zari et al. 2018), the typical distances separating two associations in the Solar Neighborhood range between 50 pc and a few hundreds of pc, that is several times the mean distance deduced from the integrated surface densities of OB stars (∼1.6 × 10 -3 pc -2 , Maíz-Apellániz 2001). Interestingly, such distances are not only comparable to the heights of the molecular gas (∼75 pc) and the cold HI gas (∼150 pc) above the Galactic plane deduced from CO and HI all-sky surveys (e.g., Dame et al. 2001, Dickey & Lockman 1990, Kalberla & Kerp 2009), but they also correspond to the typical size of HI superclouds (Elmegreen & Elmegreen 1987). For all these reasons, we, therefore, adopt a fiducial simulation with L = 200 pc and explore values down to a few tens of pc.

The mean density of the gas, n H , represents the mass of the diffuse neutral ISM contained in a volume L 3 , and also controls the porosity of the matter to the impinging radiation field. In this work, we adopt a fiducial value n H = 2 cm -3 , a value slightly larger than the standard Galactic midplane density of HI at a galactocentric distance of 8.5 kpc (Kalberla & Kerp 2009). With our fiducial value of L, the mass of gas contained in the box accounts for all the mass surface density of HI in the Solar Neighborhood (Σ HI ∼ 10 M ⊙ pc -2 , Nakanishi & Sofue 2016, Miville-Deschênes et al. 2017).

G 0 controls the intensity of the radiation field, hence the thermodynamical state of the gas. Since G 0 is normalized to the Habing field, we choose a fiducial value of G 0 = 1. The corresponding UV energy density of the standard setup is therefore slightly smaller than that contained in the standard UV radiation fields given by Draine (1978) and Mathis et al. (1983).

The standard values of the two parameters F and B x are set to 9 × 10 -4 kpc Myr -2 and 3.8 µG, respectively. As we will show later, those values are chosen so that the velocity dispersion of the gas and the strength of the magnetic field are close to the values observed in the diffuse ISM.

Steady-state

The evolution of the multiphase environments simulated here is identical to the description already given in many papers (e.g., Seifried et al. 2011;Saury et al. 2014;Valdivia et al. 2015Valdivia et al. , 2016)). Starting from an homogeneous density n H = n H , the gas evolves under the joint actions of turbulence, gravity, and thermal instability, and splits up in three different phases at thermal pressure equilibrium, the WNM, the CNM, and the LNM. The formation of dense environments well shielded from the destructive UV radiation field triggers the formation of H 2 and the medium jointly evolves from a purely atomic state to a partly molecular state. If the mass of the gas is conserved, as imposed by the periodic boundary conditions, the medium progressively tends toward a steady-state where the mean pressure, the volume filling factors of the different phases, their velocity dispersion, and their mean molecular fractions are roughly constant. This steady-state is typically reached after a few turnover timescales, providing that the corresponding time is longer than the damping time (see Sect. 3.3).

Because of turbulence and thermal instability, the steadystate has a statistical nature. The turbulent forcing and the subsequent turbulent cascade induce pressure variations and shear motions at all scales which trigger mass exchanges between the different phases. Any pressure or density structure is therefore a transient system which is usually described by its contribution to probability distribution functions. At steady-state, only PDFs remain constant. This steady yet ever changing environment is the reason for the sustained presence of a substantial amount of gas in the LNM at densities and temperatures out of thermal equilibrium (e.g., Marchal et al. 2019). All the results shown throughout this paper are taken at steady-state, at times ranging from a few tens of Myr up to 100 Myr depending on the strength of the turbulent forcing.

Properties of the multiphase medium

The steady-state values of the mean pressure P/k , the turbulent velocity dispersion σ tur , and the fractions of mass of the WNM and the CNM, f WNM and f CNM , are shown in Fig. 4 for a set of 60 different simulations. For the sake of simplicity, we assume A36, page 6 of 33 2).

that the WNM is composed of all cells with a temperature T 3000 K, the CNM of all cells with T 300 K, and the LNM of all cells with 300 K < T < 3000 K, hence

f WNM = ρ |T 3000 K ρ ( 17
)
and

f CNM = ρ |T 300 K ρ . ( 18
)
The mean pressure P/k (with k the Boltzmann constant) is classically computed as an average over the entire volume. While the above conventions are well established, there are many ways to define the turbulent velocity dispersion σ tur . It could be defined as the volume-weighted or mass-weighted velocity dispersion computed over the entire volume (Audit & Hennebelle 2005;Federrath et al. 2010), the average of the mass-weighted velocity dispersion computed along individual lines of sight (Miville-Deschênes & Martin 2007;Saury et al. 2014), or the dispersion of the mass-weighted velocity centroids computed along independent directions (Henshaw et al. 2019). All these definitions give velocity dispersions that differ from one another and are not equally relevant for the comparison with observed quantities.

To relate the velocity dispersion to a kinetic energy and provide values that could be compared to the observations of broad HI emission profiles at high Galactic latitude (see below), we chose to compute σ tur in the WNM only as 20) and where the sums are performed over all cells with T 3000 K. The results displayed in Fig. 4 are very similar to those obtained in the parametric studies of Seifried et al. (2011) and Saury et al. (2014) and are in line with the expectations of models of turbulent multiphase environments (Wolfire et al. 2003;Ostriker et al. 2010). The mean pressure of the gas is primarily regulated by the thermal equilibrium curve (see Appendix B) which depends on the local illumination of the gas by the UV radiation field G eff . Since G eff results from the absorption of the external UV field by the surrounding environments, P/k is not only sensitive to G 0 but also to the total mass of the simulation set by n H and L. Larger values of G 0 or smaller values of L or n H leads to larger P/k . In turn, the fractions of mass contained in the WNM and the CNM are controlled by the mean pressure and the total mass of the gas. Larger pressure implies larger densities of the WNM (and the CNM). The fraction of mass of the Fig. 9. Top frame: comparisons of the observational dataset (black points) with the 2D probability histograms of N H and N(H 2 ) computed from the reconstruction algorithm (see Sect. 3.8) applied to the simulations. Bottom frame: fraction of lines of sight (%), and mean value µ and dispersion σ of the logarithm of the molecular fraction computed from the simulated histograms in the regions A, B, C, D, and E defined in Table 1. Numbers correspond to the values of %, µ, and σ. The color code indicates a measure of distance (in arbitrary units) between the observed and simulated values in order to guide the eye. These comparisons are shown in each frame for 15 different simulations with G 0 varying from 0.5 (left panels) to 4 (right panels) and n H varying from 0.5 cm -3 (top panels) to 4 cm -3 (bottom panels). All other parameters are set to their fiducial values (see Table 2).

σ 2 tur = 1 3 ρ||u -u|| 2 ρ (19) with u = ρu ρ (
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4
A36, page 11 of 33 A&A 643, A36 (2020) to upper limits on N(H 2 ). Finally, the simulated sample clearly shows that a substantial fraction of the lines of sight are in the region E, with an integrated probability of 4%. If the observational sample of 360 targets is unbiased, between 10 and 18 lines of sight should have been observed in this region, which is not the case. These discrepancy will be discussed in more details in Sect. 5.

Impact of the resolution

The comparison of the two panels of Fig. 8 shows that the kingfisher diagram is independent of the resolution over about one order of magnitude of scales, from 64 3 to 512 3 . Even if not systematically shown, this unusual result is not limited to the fiducial setup but is a general feature of all the simulations explored in this work (see Fig. D.2 for instance). Our interpretation is based on the analytical model presented in Appendix C and summarized in the previous section. Evidently, high resolution simulations are important to accurately describe small scale structures. In particular, large resolution are required for modeling the formation of dense and gravitationally bound environments and follow their collapse. The fact that the kingfisher diagrams are independent of resolution for R 64 3 therefore suggests that the cold and dense structures between 0.4 and 3 pc have no influence on the distributions of N H and N(H 2 ) for the fiducial simulations and are insignificant in the mass and the volume budgets of H and H 2 . This is in line with the conclusions deduced from the analytical model. Indeed, as explained in the previous section, the turbulent forcing at large scale induces density fluctuations in the diffuse gas that extend over ∼20 pc and a distribution of dense structures with sizes smaller than ∼10 pc. While the total quantity of matter is inferred to be contained in the diffuse and the largest dense components, H 2 is exclusively built up in dense components smaller than ∼10 pc. The results of Fig. 8 combined with conclusions deduced from the analytical model therefore imply that the structures contributing the most to the mass and volume of H 2 are above 3 pc and below 10 pc for the fiducial simulation.

Beside the physical insights on the typical scales participating to the build-up of column densities, this result also provides a strong justification for using simulations with moderate numerical resolution for the study of the HI-to-H 2 transition. In all this work we therefore adopt a standard resolution of R = 256 3 unless indicated otherwise.

Impact of G 0 and n H

Figures 9 and 10 show comparisons between the observational dataset and the 2D PHs extracted from the simulations for 0.5 cm -3 n H 4 cm -3 , 0.5 G 0 4, and two sets of values of the box size and the turbulent forcing, L = 200 pc and F = 9 × 10 -4 kpc Myr -2 (Fig. 9) and L = 50 pc and F = 3.6 × 10 -3 kpc Myr -2 (Fig. 10). While the trend and the statistics of the HI-to-H 2 transition weakly depend on the resolution, they strongly depend on the total mass of the gas parametrized by n H and on the UV illumination factor. As G 0 increases or n H decreases, (1) the fraction of lines of sight with large f (H 2 ) drops to the benefit of lines of sight with low f (H 2 ), (2) the transition is shifted toward larger total column density and its width increases, and (3) the molecular fraction globally decreases over all lines of sight while its dispersion increases. Interestingly, and because n H and G 0 have opposite effects, the simulations that reproduce the most accurately the observed statistics of the HI-to-H 2 transition follow a trend with G 0 /n H ∼ 0.5-1. While similar, the effect of these two parameters are, however, not symmetrical. In particular, n H has an obvious and strong impact on the fraction of lines of sight in region E, regardless of G 0 . Likewise, the fraction of lines of sight at low column densities (regions A and B) and the mean molecular fraction at large column densities (regions C and D) are not constant for a given G 0 /n H ratio but respectively decrease and increase with n H . All these properties effectively break the degeneracies between the two parameters. All things considered, the tightest concordance between observed and simulated data is obtained for n H ∼ 1-2 cm -3 , in agreement with the Galactic midplane density deduced from HI surveys in the Solar Neighborhood (Kalberla & Kerp 2009).

At first sight, the results described above seem obvious as they somehow mimic the dependencies of the HI-to-H 2 transition found with detailed models of photodissociation regions (e.g., Krumholz et al. 2009;Sternberg et al. 2014). Such an interpretation is, however, a dangerous misconception. Indeed, while G 0 is tightly linked to the effective radiation field G eff that locally illuminates the gas, the mean density n H should not be mistaken for the local density. Moreover, the results displayed in Figs. 9 and 10 cannot be compared to PDR models because they are statistical in nature. For instance, simulations at high G 0 do not preclude the existence of clouds with high molecular fractions. Indeed, increasing G 0 may lead to denser local environments with larger molecular fractions whose probability of occurrence along a line of sight is simply reduced. It implies that the results shown here are very different from PDR model predictions. They rather reflect the complex link between the global properties of the simulation (mass, illumination, driving scale) on the one side, and the local conditions and their probability distribution functions on the other side.

Because the local abundance of H 2 is inversely proportional to G 0 , increasing G 0 naturally reduces the local self-shielding. Similarly, increasing G 0 or decreasing n H reduce the large-scale self-shielding. As a result, the density threshold n tr H required to produce highly molecular environments (see item 3. of Sect. 3.9) rises by a factor of three when either G 0 is multiplied or n H is divided by a factor of eight. While significant, such an effect on the local conditions of the gas is, however, too shallow to fully explain the variations observed in Figs. 9 and 10.

Indeed, regardless of local conditions, G 0 and n H have a major impact on the probabilistic reconstruction of lines of sight. As shown in Sect. 3.7, increasing G 0 or decreasing n H strongly reduce the fractions of mass and volume occupied by the dense and cold gas. This not only reduces the molecular fraction averaged over the entire simulation (red star in Fig. 7) but also favors the occurrence of lines of sight with low or intermediate f H 2 (cases (a) and (b) in Fig. 7). When combined with the distribution of sizes l los , the HI-to-H 2 transition is naturally shifted toward larger N H , and the asymptotic molecular fraction at high N H drops. Moreover, because the central limit theorem requires larger lines of sight to apply, the HI-to-H 2 transition is naturally wider and the dispersion of lines of sight at large molecular fraction increases. This final property is particularly well seen in the kingfisher diagram obtained for the simulation at G 0 = 0.5 and n H = 4 cm -3 where most of the lines of sight follow an homothetic transformation of the mean normalized column densities N H /l los and N(H 2 )/l los (red star in Fig. 7).

Impact of the box size L

The impacts of the box size revealed by comparing Figs. 9 and 10 partly follow the results of the previous section. Reducing L by a A36, page 12 of 33 E. Bellomi et al.: 3D chemical structure of diffuse turbulent ISM. I. factor of four drastically reduces the total amount of matter in the simulation, hence the absorption of the impinging UV radiation field and the large-scale self-shielding. As shown in Sect. 3.7, the mean pressure of the gas rises while the mass and volume occupied by the CNM decreases. All the local and statistical effects described in the previous section therefore apply and modify the kingfisher diagrams accordingly. Changing L has, however, two additional and specific consequences.

G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm
Because L drive is four times smaller in the simulations displayed in Fig. 10 than in those displayed in Fig. 9, the decorrelation scales y diff dec and y dens dec are correspondingly smaller. According to the interpretation given in Sect. 3.9, all the reconstructed lines of sight are therefore considerably larger than individual density fluctuations. This not only favors the occurrence of lines of sight at intermediate and high molecular fractions (cases (b) and (c) of Fig. 7) but also magnify the impact of the central limit A36, page 13 of 33 A&A 643, A36 (2020) assess whether this result reveals an actual and important statistical discrepancy or a simple observational limit. It implies that the likelihood of a simulation to be representative of the local diffuse gas cannot be estimated from this criterion alone. It must involve other observational signatures such as the average and the dispersion of the molecular fraction in regions A, B, C, and D, or even the chemical and statistical signatures of other atomic and molecular lines. This latter aspect is currently under development and will be the subject of the next paper of this series.

The fact that the limitations differ between FUSE and Copernicus surveys finally raises the question of the validity of studying the two samples simultaneously. We find, however, that performing comparisons with simulations on the two samples separately gives very similar results and does not impact any of our conclusions. It is so because the observational bias discussed above occurs at an extinction which increases as the natural logarithm of the square root of the instrument sensitivity. The largest total column density probed by FUSE is thus only three times larger than that observed by Copernicus (Gillmon et al. 2006;Rachford et al. 2009). Moreover, the number of lines of sight observed by FUSE that are above the maximum extinction seen by Copernicus corresponds to a small fraction of the entire sample (Gillmon et al. 2006).

Variations of physical conditions in the local ISM

The reconstruction of the simulated sample of lines of sight and the subsequent comparisons with the observations are done assuming that the local ISM can be built out of a single simulation (see Sect. 3.8). This approach was chosen in order to highlight the natural variations induced by turbulence and thermal instability alone in a diffuse neutral gas with a known averaged density and UV illumination factor. However, because the medium probed by the observations extends in all directions around the sun over a maximum distance of 3 kpc (see Fig. 2), it stands to reason that potential variations of all parameters should be taken into account, not only from one line of sight to another but also along a single and outstretched line of sight. Indeed, such considerations would offer a natural explanation for observations whose existence is in contradiction with the corresponding simulated probability of occurrence (see Sect. 4.1).

The total proton mass surface density deduced from HI and CO all sky surveys appears to be rather constant in the Galactic layer located between 5 and 12 kpc from the Galactic center (Fig. 9 of Miville-Deschênes et al. 2017). Taking into account variations of the ISM scale height above the Galaxy (Kalberla & Kerp 2009), the midplane proton density is expected to vary by less than a factor of two over the corresponding volume.

Using the Galactic star distribution of Wainscoat et al. (1992) and the grains composition and size distribution of Weingartner & Draine (2001), Porter & Strong (2005) and Moskalenko et al. (2006) estimated the energy density of the radiation field across the Galaxy. Similarly to the midplane density, the mean UV radiation field is estimated to vary by a factor two over the volume considered in this paper (Fig. 2 of Porter & Strong 2005). Interestingly, this estimation is far smaller than the variations derived by Jenkins & Tripp (2011) (Fig. 8 in their paper) from the observation of the fine structure line of CI in the local gas. Such discrepancies could be explained by the fact that Jenkins & Tripp (2011) perform local measurements: the observed gas could thus be located close to or far from an irradiating star. Alternatively, we note that the results obtained by Jenkins & Tripp (2011) are derived from models at equilibrium which do not take into account the uncorrelated perturbations of pressure and density in a turbulent multiphase medium: this naturally favors large fluctuations of the UV radiation field.

The cosmic ray ionization rate inferred from submillimeter observations of several molecular ions, including OH + , H 2 O + , ArH + , and H + 3 , shows a wide distribution across the Galactic disk (Indriolo et al. 2015). A recent estimation performed by Neufeld & Wolfire (2017) suggests that this rate could vary by a factor of five in the gas located between 5 and 12 kpc from the Galactic center.

Finally, potential variations of the composition and the size distribution of grains in the local diffuse gas should also be considered. While the extinction curve is found to be surprisingly uniform in the Milky Way (Schlafly et al. 2016), the local distribution of grains could change, not only along the line of sight but between the atomic and ionized phase and the molecular clouds. This would modify the efficiency of the photoelectric effect and the equilibrium between the two stable states of the neutral ISM, and would also have an impact the H 2 formation rate.

Interestingly, if uncorrelated, the expected variations of n H and G 0 alone would help to explain the slight discrepancies described in Sect. 4.3 but it would also lead to a dispersion of lines of sight far larger than those observed (see Figs. C.2 and C.3). The fact that the predicted statistics of the HI-to-H 2 transition is close to the observed sample therefore suggests that the variations of all the parameters described above cannot be considered independently but must follow strong correlations which apply locally (as discussed, for instance, by Bialy et al. 2019) but also across the Galactic disk.

Fraction of ionized gas

The fraction of volume ϕ occupied by the ionized phases of the ISM, the warm ionized medium (WIM) and the hot ionized medium (HIM), plays an important role in our reconstruction algorithm. As illustrated in Fig. 5, this parameter controls the filling factor of the neutral medium along a line of sight of length l los . Unfortunately, its value in the Milky Way is highly uncertain and still debated.

The consensus is that the volume of the HIM far exceeds that of the WIM and results from an interplay between supernovae explosions, which regularly produce hot gas in the Galactic disk, and buoyancy, which lifts this gas into the halo, releasing the pressure in the midplane. Early analytical studies neglecting buoyancy (McKee & Ostriker 1977) predicted a large fraction of HIM in the midplane (ϕ ∼ 95%). In contrast, early numerical simulations, including the cycle of matter and energy between the disk and the halo, led to considerably smaller predictions with ϕ ∼ 25% (de Avillez & Breitschwerdt 2004). This value is now considered as a lower limit by the most recent numerical simulations which reveal the importance of the driving mode of supernovae explosions (Walch et al. 2015) and of the photoelectric heating (Li et al. 2015;Hill et al. 2018) on the volume of the HIM. These latest studies estimate that 20% ϕ 90% in the Galactic midplane.

Highly uncertain, the volume filling factor of the HIM can also vary from one line of sight to another (Fig. 1 of Hill et al. 2018). In this paper, we adopt a constant and conservative value ϕ = 0.5 for every line of sight. Changing ϕ would have the effect of either squeeze or stretch the predicted 2D PHs displayed in Figs. 9 and 10 along the x axis, and to modify the balance of probabilities of occurrence of lines of sight at high and low molecular fraction. Taking into account a realistic distribution of ϕ would require to simulate the Galactic disk and halo over a scale of several kiloparsecs and to properly model and follow the impact of supernovae explosions. This is far beyond the scope of the present paper.

Galactic vertical structure

The simulated 2D PH of the HI-to-H 2 transition are found to slightly depend on the Galactic gravitational potential. This odd result is nothing but an artifact of the physical scales considered here. Since the sizes of all simulations are below 200 pc, the gas expands, at the most, over 100 pc above the Galactic plane, a distance far smaller than the characteristic scale of variation of the thermal pressure expected for a turbulent gas in hydrostatic equilibrium. It follows that the column densities show no significant variation as a function of the position of the line of sight or its angle with the Galactic plane.

This setup, initially chosen to favor the physical resolution, is a strong shortcoming which prevents us from using and studying the information carried by the Galactic latitude of each observation. Indeed, the comparison between the FUSE halo survey and the data collected by FUSE and Copernicus in the Galactic disk indicates that the HI-to-H 2 transition at high latitude occurs at a total hydrogen column density ∼2 times smaller than that in the disk (Gillmon et al. 2006). Similarly to the previous section, studying this effect would requires to compute the local vertical structure of our Galaxy over several hundreds of parsecs, taking into account the hot and ionized component of the ISM.

Doppler broadening parameter

The self-shielding of H 2 depends on the dispersion of the Lyman and Werner lines which is usually modeled with a turbulent Doppler broadening parameter b D (Draine & Bertoldi 1996 and Eq. ( 16)). In single cloud models, this parameter has the effect of shifting the HI-to-H 2 transition toward larger total column density without modifying any of the asymptotic states (e.g., Bialy et al. 2017). In this paper, we identify this parameter with the turbulent velocity dispersion of the gas at large scales and therefore adopt b D = 8 km s -1 for the fiducial simulation. This is done to prevent an overestimation of the H 2 self-shielding at large scales, at the cost of underestimating the self-shielding at the scale of a CNM clouds.

To estimate the effect of this parameter, all the grids presented in this paper were also run assuming b D = 2 km s -1 , which roughly corresponds to the velocity dispersion expected at a scale of 10 pc for the fiducial simulation. While locally important, b D is found to have a relatively small impact on the 2D PHs of N H and N(H 2 ): increasing b D by a factor of four slightly increases the width of the HI-to-H 2 transition and the fraction of lines of sight in region B. We interpret this limited effect as a consequence of the fact that the asymptotic molecular states of any cloud are independent of b D .

Even so, it should be stressed that b D has a strong impact on the local molecular fraction in transition regions. Therefore, and while inconsequential for the global statistics of H and H 2 , the value of the Doppler broadening parameter might be paramount for any chemical species preferentially formed at the border of molecular clouds. As proposed by Bialy et al. (2019), the H 2 selfshielding could be computed self-consistently using the velocity and density fields of the simulation and a cost effective radiative transfer method. This would prevent the dilemma of favoring large-scale or small-scale self-shielding.

H 2 self-shielding at high temperature

The prescription of H 2 self-shielding used in this paper (Eq. ( 16)) is taken from Draine & Bertoldi (1996). As discussed by Wolcott-Green et al. (2011), such a prescription is reliable for diffuse gas at low temperature but becomes less and less reliable for high temperature environments (T > 500 K) where efficient collisional excitation of H 2 in its rovibrational levels reduces the self-shielding. To estimate the impact of this process, we ran the fiducial simulation with the alternative self-shielding function proposed by Wolcott-Green et al. (2011) (Eq. (12) in their paper). This prescription leads to a similar probability histogram and therefore does not influence the global analysis of the kingfisher diagram presented in this paper. However, we note that adopting this alternative prescription slightly increases the width of the HI-to-H 2 transition, and induces more lines of sight at intermediate molecular fraction (region B, see Fig. 3), in better agreement with the observations.

Summary and conclusions

This paper presents an exhaustive parametric study of the transition from atomic to molecular hydrogen in the local diffuse ISM. Using state-of-the-art MHD simulations, and an ensemble of 305 runs, we quantify separately the impacts of the mean density, the UV radiation field, the integral scale, the resolution, the turbulent forcing, the magnetic field, and the gravity on the molecular content of multiphase environments. The original feature of this work is to not only focus on the production of individual column densities but also on their statistics, meaning the probabilities of occurrence of these column densities along random lines of sight. For the first time, both the chemical and statistical information are used concomitantly, through the socalled kingfisher diagrams, to interpret the distribution of H and H 2 observed toward 360 lines of sight across the local interstellar matter.

The results of the simulations are interpreted with a semianalytical model which attempts to separate the effects of local conditions from those induced by the probabilistic reconstruction of individual lines of sight. To compare the kingfisher diagrams to the observational sample, we propose a new version of the Kolmogorov-Smirnov test which can be generalized and used for the comparison of two probability histograms or distribution functions in any dimension larger than one.

Taking into account the distance of each background source and simulating random lines of sight over the same distribution of distances is paramount to explain the range of observed column densities and their corresponding statistics. Once this aspect is included, the joint actions of thermal instability and large-scale turbulence in the standard simulation are found to produce a wealth of lines of sight which reproduce the observed position and width of the HI-to-H 2 transition, and whose probabilities of occurrence match those derived from the observations. The agreement is so remarkable that it is almost unnecessary to invoke variations of physical conditions along lines of sight or from one line of sight to another.

The minimal KS distance obtained over the entire grid is ∼0.5. Such a value implies that there exists a small group of lines of sight in the observational sample whose probability of occurrence is under-or over-predicted by about a factor of three. However, it also implies that the probability of occurrence of any other group of observed lines of sight, small or large, is reproduced to a better level.

A36, page 17 of 33 A&A 643, A36 (2020) The distribution of column densities computed from the simulations strongly depends on the Galactic midplane density parametrized by the mean density n H , the density of OB stars parametrized by the UV scaling factor G 0 , and the scale of neutral diffuse clouds parametrized by the box size L. It is so because these three parameters not only regulate the mean pressure of the gas, hence the fractions of mass and volume occupied by the CNM and WNM, but also control the typical scale of density fluctuations in the WNM and the distribution of sizes of the CNM structures where H 2 forms. The tightest concordance between the observed and simulated samples is obtained for a mean density n H = 1-2 cm -3 and a UV radiation field scaling factor G 0 = 1 (in Habing units), in good agreement with the values deduced from HI and CO all sky surveys and from direct observations of the UV radiation field in the Solar Neighborhood. The range of observed column densities of H and H 2 requires a box size L = 200 pc which corresponds to the estimated scale of HI superclouds.

Within this setup, the column densities of HI are inferred to be built up in large-scale WNM and CNM structures correlated in densities over ∼20 and ∼10 pc, respectively. In contrast H 2 is inferred to be built up at smaller scales. However, the fact that the kingfisher diagram is independent from the resolution of the simulation suggests that most of the mass and volume of H 2 is contained in CNM structures between ∼3 and ∼10 pc. All these values are given for the standard simulation (L = 200 pc and L drive ∼ 100 pc) but naturally depend on the size of the box and the mean density of the gas.

In spite of the strong influences of n H , G 0 , and L, the statistical properties of the HI-to-H 2 transition are otherwise remarkably stable. Admittedly, the kingfisher diagram depends on the strength of the turbulence if most of the forcing is injected in solenoidal modes; however such a configuration prevents to reproduce the observational sample unless the large-scale velocity dispersion of the gas is unrealistically small. In contrast, if most of the kinetic energy is injected in compressive modes, the kingfisher diagram is found to weakly depend on the strength of the forcing. Similarly, the HI-to-H 2 transition is almost not affected by gravity and is found to weakly depend on the Doppler broadening parameter and the strength of the magnetic field, as long as B x 4 µG. The 2D PH of the column densities of H and H 2 is therefore a valuable tool to constrain the nature of the turbulent forcing at large scale; however, it provides few or no information regarding the velocity dispersion of the gas, the amount of gravitationally bound environments and the strength of the magnetic field. Other observational tracers are required.

All these results open new perspectives for the study of the chemical state of the ISM in which any observation must be understood through the combination of local physical conditions and the probabilistic ordering of these conditions along the line of sight. In particular, similar studies should be applied to all atomic and molecular species with observational samples large enough to conduct statistical analysis. It also invites to expend the study of PHs to higher dimensions, taking into account simultaneously the joint information contained in the column densities of several species. In this context, the generalization of the Kolmogorov-Smirnov test proposed in this paper will be very valuable. All these aspects are currently under development and will be the subjects of the following papers of this series. the emission profiles of the 21cm line can also be highly uncertain because the measurements are done over a beam far larger than the pencil-beam sampled by H 2 data but also because it requires to identify in the HI profiles the components associated to the molecular gas. It should thus be kept in mind that while the errors on the column densities of H 2 are somehow limited, those on HI can be sometimes larger than a factor of five, in particular for the lines of sight at high Galactic latitude (Gillmon et al. 2006). The color excess E(B-V) given in Table A.1 finally results from a compilation which includes direct measurements of the star reddening compared to its intrinsic (B-V) 0 color (e.g., Savage et al. 1977;Fitzpatrick & Massa 1990;Diplas & Savage 1994;Rachford et al. 2002) and dust emission maps from the IRAS telescope4 (Schlegel et al. 1998).

With all these data at hand, we adopt the following methodology to derive the total proton column densities N H : if the column densities of HI and H 2 are available, then N H is computed as N(H) + 2N(H 2 ); if not, N H is derived from the reddening E(B-V) as N H = 5.8 × 10 21 E(B-V) cm -2 assuming a standard Galactic extinction curve and the average interstellar ratio R V = A V /E(B-V) = 3.1 (Fitzpatrick & Massa 1986;Fitzpatrick 1999). It should be noted that the N H /E(B-V) ratio observed at low column density is larger than the standard value used in this work (Liszt 2014;Lenz et al. 2017). Because of this and the uncertainties on the HI column densities discussed above, the values of N H derived here should be considered as estimates.

Examples of the underlying uncertainties on N H can be seen in Table A.1 where N(H) sometimes exceeds slightly the total column density derived from E(B-V). (32) 21.18 (32) 20.32 (43) 21.29 HD 133.960 -4.990 1.692 0.411 (21) 21.34 (18) 20.51 (43) 21.45 HD 134.580 -4.960 2.268 0.460 (18) 21.26 (18) 20.67 (43) 21.44 HD 135.080 -3.820 2.558 0.480 (18) 21.45 (18) 20.43 (43) 21 (23) 21.01 (42) 20.48 (42) 21 (5) 21.18 (42) 20.38 (42) 21 (18) 21.16 (32) 20.85 (31) 21.46 HD 207.070 -16.260 0.339 0.717 (21) 20.91 (41) 20.64 (46) 21.23 HD 237.290 -27.100 0.466 0.046 (35) 19.75 (3) 15.51 (2) 19.75 HD 214.510 -18.500 0.520 0.070 (3) 20.52 (3) 15.68 (2) 20.52 HD 194.070 -5.880 3.378 0.300 (45) 21.30 (18) 19.53 (38) 21.31 HD 183.970 0.840 0.663 0.150 (3) 21.08 (8) 19.73 (2) 21.12 HD 180.090 4.340 4.000 0.460 (46) 21.50 (41) 20.58 (41) 21.59 HD 041117 189.650 -0.860 1.000 0.450 (8) 21.40 (8) 20.69 (41) 21.54 HD 187.750 1.770 1.400 0.290 (45) 21.40 (18) 20.52 (41) 21.50 (46) 21.27 (46) 20.87 (46) 21.52 HD 241.630 -20.780 0.616 0.020 (45) 20.30 (7) <14.85 (7) (45) 20.26 (7) 14.95 (7) 20.26 HD 262.060 -10.420 0.363 0.050 (45) 20.26 (7) <14.60 (7) 20.46 HD 266.680 -12.320 0.139 0.020 (7) <20.74 (7) <14.78 (7) 20.06 HD 263.480 -10.280 1.117 0.080 (7) 20.54 (7) 15.08 (7) 20.54 HD 245.430 2.050 5.747 0.200 (41) 21.23 (41) 19.72 (41) 21.26 HD 255.980 -4.710 0.668 0.040 (3) 19.95 (18) 14.45 (2) 19.95 HD 262.800 -7.690 0.479 0.040 (10) 19.78 (10) 14.23 (2) 19.78 HD 254.520 -1.330 1.292 0.180 (26) 21.06 (42) 19.73 (41) 21.10 HD 266.830 -5.820 1.718 0.360 (8) 21.18 (32) (45) 21.15 (18) 20.26 (38) 21.25 HD 265.860 -1.900 2.041 0.270 (32) 21.08 (32) 20.20 (32) 21.18 HD 267.360 2.250 1.299 0.300 (32) 21.20 (32) 20.72 (32) 21.42 HD 277.690 -7.370 0.151 0.040 (7) 20.78 (7) <17.90 (7) (45) 21.00 (18) 19.30 (29) 21.02 HD 285.160 -2.130 3.205 0.280 (26) 21.15 (18) 19.92 (41) 21. (45) 21.34 (8) 19.70 (41) 21.36 HD 091651 286.550 -1.720 1.934 0.300 (18) 21.15 (18) 19.07 (41) 21.16 HD 285.700 0.070 2.331 0.270 (32) 21.12 (32) 19.85 (32) 21.16 HD 285.880 0.050 4.255 0.260 (32) 21.17 (32) 20.14 (32) 21.24 HD 287.600 -2.020 4.587 0.340 (45) 21.28 (18) 18.93 (18) 21.55 (18) 19.83 (38) 21.57 HD 287.570 -0.710 2.227 0.420 (8) 21.40 (8) 19.77 (44) 21.42 (18) 21.34 (8) 19.83 (43) 21.37 HD 287.670 -0.940 1.101 0.330 (45) 21.34 (18) 19.52 (44) 21.35 HD 287.740 -1.020 2.941 0.370 (26) 21.40 (29) 19.84 (43) 21. (5) 21.04 (18) 19.28 (43) 21.05 HD 288.240 -0.900 2.625 0.340 (8) 21.33 (18) 19.61 (41) 21.35 HD 295.690 -14.730 0.267 0.180 (28) -20.70 (42) 21.02 HD 272.830 29.170 0.387 0.140 (25) 20.90 (25) 19.06 (25) (46) 20.54 (27) 20.50 (27) 20.99 HD 296.850 -1.020 2.381 0.210 (18) 21.16 (18) 19.82 (41) 21.20 HD 297.450 -0.340 2.315 0.220 (18) 21.11 (18) 19.98 (39) 21.17 HD 298.230 3.790 0.086 0.020 (45) 20.04 (7) <14.08 (7) (45) 21.35 (42) 19.95 (42) 21.38 HD 301.920 -15.360 0.341 0.220 (46) 20.86 (27) 20.49 (27) 21.13 HD 301.710 -9.880 2.755 0.260 (18) 21.04 (19) 20.04 (41) 21.12 HD 301.960 -0.200 0.420 0.520 (12) 20.85 (27) 20.64 (27) 21. (9) 21.34 (42) 20.34 (43) 21.42 HD 305.760 0.150 2.101 0.490 (18) 21.36 (42) 20.63 (42) 21.50 HD 306.060 0.220 2.268 0.400 (45) 21.41 (18) 20.58 (43) 21.52 HD 308.230 10.680 1.675 0.130 (45) 21.04 (18) 19.63 (38) 21.07 HD 316.000 51.000 0.084 0.030 (3) 18.83 (13) 12.95 (2) 18.83 HD 116781 307.050 -0.070 2.045 0.340 (41) 21.18 (41) 20.08 (41) 21.24 HD 304.880 -16.130 22.727 0.210 (31) 20.96 (32) 19.83 (43) 21.02 HD 118716 310.190 8.720 0.168 0.040 (10) 19.90 (10) <14.08 (7) 20.37 HD 100.700 65.320 0.030 0.080 (45) <20.90 (8) 13.38 (2) 20.67 HD 314.240 19.120 0.119 0.100 (7) 20.40 (7) <14.78 (7) 20.76 HD 314.070 14.190 0.120 0.020 (45) 19.28 (19) 12.80 (2) 19.28 HD 333.970 55.840 3.425 0.090 (26) 20.71 (18) 18.70 (39) 20.72 HD 311.770 1.250 0.160 0.060 (45) 19.52 (3) 12.80 (2) 19.52 HD 312.260 1.790 2.387 0.298 (35) 21.26 (32) 20.31 (42) 21.35 HD 312.670 -0.420 1.808 0.530 (18) 21.39 (42) 20.52 (43) 21.49 HD 322.770 16.670 0.095 0.050 (7) 20.11 (7) <14.18 (7) 20.46 HD 320.130 -2.640 0.835 0.220 (3) 21.08 (3) 19.77 (2) 21.12 HD 352.020 39.230 0.093 -19.38 (19) 14.34 (6) 19.38 HD 331.320 13.820 0.115 0.020 (10) 20.18 (7) <14.26 (7) 20.06 HD 336.720 18.860 0.822 0.250 (5) 21.00 (42) 20.56 (42) 21.24 HD 333.190 11.890 0.137 0.030 (7) 20.23 (7) <14.26 (7) 190 (8) 21.15 (10) 19.41 (2) 21.17 HD 353.190 23.600 0.161 0.210 (35) 21.09 (10) 19.83 (2) 21.13 HD 352.750 22.770 0.142 0.220 (10) 21.18 (10) 20.04 (2) 21.24 HD 339.040 8.420 0.266 0.350 (28) 21.07 (42) 20.77 (42) 21.37 HD 354.610 22.700 0.135 0.270 (3) 21.15 (3) 19.89 (2) 21.20 HD 351.310 17.000 0.100 0.400 (8) 21.38 (18) 19.79 (2) 21.40 HD 352.450 17.630 0.181 0.640 (16) 21.43 (16) 20.78 (16) 21.59 HD 344.860 10.090 0.295 0.390 (5) 21.41 (42) 20.68 (42) 21.55 HD 352.250 16.850 0.140 0.740 (16) 21.50 (16) 20.90 (16) 21.68 HD 353.650 17.710 0.092 0.520 (32) 21.71 (32) 20.57 (32) 21.77 HD 352.860 17.040 0.139 1.090 (16) 21.46 (16) 21.37 (16) (32) 21.80 (32) 19.88 (32) 21.81 HD 353.100 15.800 0.117 0.100 (10) 20.95 (10) 18.74 (2) 20.96 HD 336.370 -0.220 1.135 0.660 (18) 21.60 (18) 20.71 (40) 21.70 HD 339.380 2.510 0.842 0.370 (3) 21.12 (18) 20.44 (2) 21.27 HD 340.540 3.010 1.316 0.680 (8) 21.40 (8) 20.79 (46) 21 190 (45) 21.32 (42) 20.36 (42) 21.41 HD 346.120 3.910 0.268 0.050 (7) 20.40 (7) <14.26 (7) 20.46 HD 343.480 1.220 2.300 0.400 (45) 21.35 (8) 20.29 (44) 21.42 HD 343.030 0.870 1.403 0.680 (8) 21.77 (18) 20.73 (46) 21 (32) 21.37 (32) 20.47 (32) 21.47 HD 344.620 1.610 1.500 0.330 (45) 21.28 (18) 20.21 (38) 21.35 HD 344.810 1.610 16.667 0.460 (18) 21.43 (18) 20.33 (39) 21.49 HD 349.970 3.220 1.217 0.820 (16) 21.00 (27) 21.16 (27) 21.59 HD 352.590 2.870 0.994 0.230 (45) 21.08 (10) 19.92 (45) 19.23 (20) 12.70 (2) 19.23 HD 301.040 -4.720 0.202 0.083 (35) 20.19 (18) <14.23 (7) (8) 21.23 (18) 19.85 (39) 21.26 HD 30.990 13.370 0.143 0.190 (7) 20.82 (7) 19.85 (7) (46) 21.95 (46) 20.19 (41) 21.96 HD 6.060 -1.200 1.185 0.310 (26) 21.18 (18) 20.00 (38) 21.23 HD 6.050 -1.330 1.235 0.380 (45) 21.20 (26) 20.22 (38) 21.28 HD 343.330 -13.820 0.279 0.060 (45) 20.85 (3) 18.95 (2) 20.86 HD 6.120 -1.480 1.276 0.360 (45) 21.36 (18) 20.20 (38) 21.42 HD 6.400 -1.560 1.996 -21.41 (42) 20.15 (42) 21. 1.205 0.120 (45) 21.11 (32) 16.53 (32) 21.11 (18) 21.15 (10) 20.28 (2) 21.25 HD 18.250 1.680 2.033 1.040 (27) 21.60 (27) 20.85 (27) 21.73 HD 16.940 0.840 2.100 0.790 (27) 21.65 (18) 20.68 (27) 21.73 HD 5.820 -6.310 2.488 0.240 (45) 21.11 (18) 20.10 (41) 21.19 HD 17.540 -0.670 2.128 1.120 (16) >19.95 (16) 21.16 (16) 21.81 HD 21.060 -0.530 0.231 0.480 (46) 21.15 (18) 20.86 (27) 21. 191 (21) 20.87 (18) 19.43 (39) 20.90 HD 88.220 -2.080 0.931 0.170 (9) 20.99 (41) 20.00 (38) 21.07 HD 80.980 -10.050 0.187 0.130 (35) 20.68 (7) 19.15 (7) 20.70 HD 87.610 -3.840 0.587 0.320 (8) 21.00 (3) 20.29 (2) 21.14 HD 100.510 8.620 2.611 0.600 (18) 21.20 (42) 20.60 (42) 21.38 HD 309.460 -31.740 0.292 0.280 (47) 21.27 (32) 20.64 (32) 21.44 HD 90.560 -2.330 0.223 0.720 (27) 21.48 (27) 21.00 (27) (16) 21.30 (27) 20.86 (27) 21.54 HD 99.800 3.620 0.958 0.440 (32) 21.09 (32) 20.44 (32) 21.25 HD 103.140 6.990 1.025 0.590 (8) 21.34 (8) 20.83 (27) 21.55 HD 103.110 6.820 1.026 0.520 (15) 21.20 (41) 20.76 (41) 21.44 HD 101.600 4.670 0.838 0.640 (46) 21.34 (18) 20.91 (27) 21.58 HD 102.710 4.980 0.911 0.520 (40) -20.87 (34) 21.48 HD 104.030 6.440 0.829 0.290 (31) 21.23 (18) 20.34 (34) 21. (32) 21.21 (32) 19.42 (32) 21.22

where R is the resolution of the box of size L (see Table 2), i 0 is a random starting index for the integration of the column densities, and i l is the final index deduced from l. Because of spatial correlations of the density n H , this computation is not equivalent to the sum of random realizations of n H drawn out of the 1D probability distribution (Fig. C.1). It depends, instead, on how and over which distance the values of n H are correlated. As proposed by Bialy et al. (2017Bialy et al. ( , 2019)), and since the correlations of density in a turbulent medium decrease over long distances, we assume that these correlations can be modeled with a parameter y dec , called the decorrelation scale. The density is supposed to be constant over distances smaller than y dec , and uncorrelated over larger distances. In this framework, if y dec is the same for all densities, the total column density integrated over a distance l becomes equivalent to the sum of 1 + l/y dec random realizations of n H . In isothermal simulations, the 1D probability distribution of the gas density is found to follow a lognormal distribution with a dispersion proportional to the Mach number, dependent on the nature of the turbulent driving, and independent on the mean density of the gas (e.g. Padoan et al. 1997;Federrath et al. 2008). Because of this property, Bialy et al. (2017) were able to establish a relation between the dispersion σ N H /l of the distribution of the averaged densities and the dispersion σ n H of the distribution of the proton density:

σ N H /l σ n H = 1 + L drive y dec l L drive -1/2 , (C.2)
where L drive is the turbulence driving scale. Fitting the σ N H /l /σ n H ratio as a function of l/L drive , Bialy et al. (2017) estimated y dec = 0.2 × L drive . 2). The black histogram correspond to the extracted data. The red dashed curve shows an example of the sum of two log-normal components and a power-law tail at high density, for comparison. The blue line indicates the inflection point of the PH between the diffuse and the dense components.

Unfortunately, this method cannot be applied here. Indeed, oppositely to isothermal simulations and as illustrated in Fig. C.1, the PH of n H derived from simulations of the multiphase ISM is usually described by the sum two log-normal distributions plus a power-law tail at high density that could be a signature of the CNM which is known to behave like a polytropic gas with an exponent γ < 1 (Passot & Vázquez-Semadeni 1998) or a signature of gravity (Federrath 2013;Girichidis et al. 2014). To overcome this issue, and for the sake of simplicity, we assume that a two phase medium is described by two decorrelation scales. All densities below a limit n lim H are supposed to belong to the diffuse log-normal component and to be correlated over a scale y diff dec . Similarly, all densities above n lim H are supposed to belong to the dense log-normal component and to be correlated over a smaller scale y dens dec . We identified here the limit n lim H between the two log-normal distributions with the inflection point of the 1D probability histogram of the gas density (blue line in Fig. C.1). Since the diffuse component behave like an isothermal gas at the temperature of the WNM, we assume that y diff dec /L drive is constant for all simulations and adopt the value given by Bialy et al. (2017), y diff dec = 0.2 × L drive , where L drive = L/2 is the main driving scale used for the turbulent forcing (see Sect. 3.3). In contrast, and because the diffuse component occupies most of the volume, we state that y dens dec depends on the total mass of the gas or equivalently its mean density n H . To simplify, we propose that

y dens dec ∝ n H 1/3 , (C.3)
which means that the impact of changing the total mass of the gas is to change the typical volume of the dense structures by the same factor. Within this framework, the semi-analytical model proposed here therefore depends on a single parameter: the decorrelation scale of the dense component for the fiducial simulation. In the following, we assume y dens dec = 10 pc for n H = 2 cm -3 which implies a decorrelation length of 6.3, 7.9, and 12.6 pc for n H = 0.5, 1, and 4 cm -3 , respectively.

It is quite optimistic to believe that the dense and cold component of the ISM can be modeled by a single decorrelation scale y dens dec . This component is indeed likely to follow a distribution 10 -2 10 -1 10 0 10 1 10 2 PH N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 PH N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 10 2 N H / l (cm -3 ) 10 -6 10 -5 10 -4 10 -3 10 -2 of sizes which decrease with the local pressure, hence the local density n H . However, and as we show below, such an assumption allows us to highlight important features of the simulations. The essential impact of the distribution of sizes is proven and discussed in Appendix C.4.

G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm -3 l=50

C.2. Comparison with simulations

The goal of the model is to offer an explanation on how the PHs of local densities translate into PHs of column densities in a simulation of the multiphase ISM. To test its validity, we generate a series of N fictitious lines of sight of size l and compare the PHs to those obtained with an equivalent sample of lines of sight extracted from numerical simulations. For each line of sight, we draw a sequence of random realizations of n H out of its known 1D PH using the rejection method. For each draw, the density is supposed to be constant over a distance y diff dec if n H < n lim H , and over a distance y dens dec otherwise 5 . The contribution of this piece of 5 To ensure that the resulting sample matches the original PH, the probability associated to each density is weighted by the inverse of the component size 1/y diff dec or 1/y dens dec .

gas to the total column density is therefore computed as n H y diff dec or n H y dens dec . In contrast, the contribution of this piece to the column density of H 2 is inferred from the expected density profile of H 2 over a 1D slab of size y diff dec or y dens dec . Throughout the slab, the density of H 2 is calculated at equilibrium taking into account local and large-scale extinction and self-shielding as

e -σ d (N loc H +N ext H ) f shield N loc (H 2 ) + N ext (H 2 ) 5 × 10 14 cm -2 , (C.4)
where N loc H and N loc (H 2 ) are the local column densities of protons and H 2 computed from the border of the slab. The mean values in this expression are calculated from six random realizations of the large-scale column densities N ext H and N ext (H 2 ) drawn out of the sample a lines of sight under construction. The construction of each line of sight ends when its size reaches the integration scale l. The entire sample of lines of sight is finally reconstructed until convergence of the large-scale shielding processes described above.

The comparisons between the 1D and 2D PHs reconstructed from the semi-analytical model and extracted from the simulations are shown in ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) lines of sight, 15 different simulations, and different values of the integration scale l, from L down to y diff dec . Unexpectedly, setting the decorrelation length of the dense component to ∼10 pc for the fiducial simulation (n H = 2 cm -3 ) leads to a remarkable agreement between all the fictitious and actuals PHs. Once this parameter is set, the model not only reproduces surprisingly well the shapes of the 1D PHs and of the HI-to-H 2 transition, but also their global trends depending on n H and G 0 and their deformations depending on the chosen integration scale l. This result is not straightforward and highly depends on the decorrelation lengths y diff dec and y dens dec . The agreement observed in Figs. C.2 and C.3 therefore suggests that the model somehow captures an essential property of the simulations, namely some characteristic lengths of the diffuse and dense components of a multiphase gas with a turbulence driven at a scale L drive .

G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm -3 l=100
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm -3 l=50
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm -3 l=25

C.3. Interpretation of the model

The statistics of the HI-to-H 2 transition derived from the model result from a combination of effects. Locally, the fraction of H 2 of a given slab depends on the density and the size of the slab y diff dec and y dens dec , which set the local self-shielding, and on the surrounding environment, which sets the large-scale self-shielding. How these local properties contribute to the integrated quantities N H and N(H 2 ) depends, in turn, on the sizes of the slabs and on the PH of the density which both control the reconstruction of the line of sight.

1. Since y dens dec is fixed, the HI-to-H 2 transition induced by the local self-shielding alone occurs in any slab with a density larger than

n tr H ∝ G 1/2 0 (y dens dec ) -1/2 ∝ G 1/2 0 (n H ) -1/6 . (C.5)
This equation can be obtained from the expression of the column densities of HI envelopes (Sternberg et al. 2014, Eq. ( 40)) in the weak field limit. ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 10 -2 10 -1 10 0 10 1 N H / l (cm -3 ) 2. Lines of sight with very low molecular fraction necessarily result from the combination of slabs with n H < n tr H . The occurrence of such events depends on the volume filling factor of the diffuse gas, hence on the 1D PH of low density material, and on the integration length l: as l increases, their likelihood decreases. Such a scenario occurs for a maximum normalized column density of

G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm -3 l=100
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm
G 0 = 0.5 G 0 = 1.0 G 0 = 2.0 G 0 = 4.0 - n H =0.5 cm -3 - n H =1.0 cm -3 - n H =2.0 cm -3 - n H =4.0 cm
N H /l = n tr H . (C.6)
For l ∼ y diff dec , such a high normalized column density of atomic gas is a likely event. As l increases, it becomes, however, unlikely to throw a line of sight composed of components with identical densities n tr H . Therefore, while the above limit is still valid, the maximum normalized column density of weakly molecular gas appears to decrease.

3. Lines of sight with large molecular fraction necessarily contain at least one slab with n H > n tr H . Oppositely to the previous case, the occurrence of such events depends on the volume filling factor of the dense gas, hence on the 1D PH of large density material, and on the integration length l: as l increases, their likelihood increases. Such a scenario occurs for a minimum normalized column density

N H /l = n tr H y dens dec /l. (C.7)
4. These two limits for lines of sight with low and large f H 2 (items 2. and 3.), set the width of the HI-to-H 2 transition seen in column densities. As shown in Fig. C.3, this width is somehow smaller than that obtained from the numerical simulations. We will discuss this point in the next section.

5. At last, lines of sight with intermediate H 2 fraction (10 -4 f H 2 10 -2 ) mostly result from a combination of slabs of low and moderate densities (n H n tr H ). Because such events are unlikely, the model predicts a small fraction of lines of sight at intermediate f H 2 , in contradiction with results extracted from the simulations. This point will also be discussed further in the next section.

All these properties fully explain the behaviors of the analytical model observed in Fig. C.3. The transition density for the fiducial simulation is n tr H = 4 cm -3 . As expected, the corresponding lower and upper limits of N H /l required to activate the HI-to-H 2 transition are in agreement with the limits found for the lowest integration scale l = 25 pc (bottom right panels of Fig. C.3). Increasing the integration length has three effects: (a) to squeeze the 2D PH along the x and y axis as both N H and N(H 2 ) progressively tend toward Gaussian distributions centered on the means, (b) to shift the HI-to-H 2 transition to lower N H /l according to the limits derived above, and (c) to increase the occurrence of lines of sight with large f H 2 to the detriment of lines of sight with low f H 2 as the probability of intercepting a slab with n H > n tr H rises. The dependence of the distributions on G 0 and n H are also straightforward. As G 0 increases or n H decreases, the HI-to-H 2 transition is shifted according to the dependence of simulated datasets. To compare these values, we define a merit function

M = log 10 f sim f obs , (D.1)
and the KS distance between the two distributions as the maximum value of M computed over all quadrants of an observational dataset O,

R KS = max O (M). (D.2)
This procedure, not only provides a measurement of the difference between the two distributions, but also the datapoint and the quadrant that maximize the merit function (see Fig. D.1). The interpretation is also straightforward. For instance, a KS distance of 1 implies that one of all the quadrants scanned contains 10 times fewer or 10 times more observations than simulated lines of sight, and that all the other quadrants have smaller distances.

The stability of the procedure depends on the errors made on the merit function and therefore on the absolute numbers of observed and simulated lines of sight contained in each quadrant. The observational dataset O used to compute the KS distance (red points in Fig. D.1) is chosen as the subsample such that all quadrants scanned contain at least 10 observed lines of sight. With this assumption, the error on the merit function is calculated by taking into account only the statistical errors on the number of simulated lines of sight. For each quadrant, we assume that the "true" merit function ranges between

M min = log 10        f sim -3 ( f sim /S ) f obs        and (D.3) M max = log 10        f sim + 3 ( f sim /S ) f obs        , (D.4)
where S is the total number of simulated lines of sight. Because the errors are asymmetric, M min can tend toward infinity. If so, the KS distance is a lower limit, even if the infinite error bar is obtained for another quadrant than the one that maximizes M.

In short, for each observational datapoint, we compute M min , M, M max , and R KS . If one of the M min goes to infinity, R KS is considered as a lower limit. Because the division in quadrants is performed on a Cartesian grid, the procedure also depends on the choice of the axes. Mathematically, the best option would be to identify the principal components of the observational sample using proper orthogonal decomposition or singular value decomposition algorithms. Such a method could even be applied to subsamples in order to adaptively rotate the system of axes and follow more precisely the distribution of observations. In any case, the resulting system would be a linear combination of N H and N(H 2 ) which could be difficult to relate to the underlying physical properties. Because the molecular fraction is bimodal as a function of N H , we choose here N H and f H 2 as primary variables. This choice facilitates the physical interpretation of the KS test while ensuring some homogeneity of the spread of the observational data in all quadrants (see Fig. As a proof of concept, we display in Fig. D.2 the results of the KS test applied to a grid of simulations obtained for different values of n H and G 0 and two different resolutions. Despite its simplicity, this test appears to capture the main behaviors described in Sects. 4.2 and 4.3. In particular, the KS distance between the simulations and the observations is found to strongly depend on both n H and G 0 , and more loosely on the resolution. Moreover, the simulations that minimize R KS are found to be identical to those identified in Sect. 4.3 to be in closest agreement with the observations. Interestingly, the value of R KS obtained for n H = 4 cm -3 and G 0 = 4 is relatively small, in apparent contradiction with the conclusions of Sect. 4.3. This is due to the limit imposed on the minimum number of observations contained in each quadrant. Because of this limit, several observations at large N H are not included in the analysis which reduces the merit function at large column density (region E, see Table 1). Keeping in mind these border effects, the KS test turns out to be a valuable tool for estimating the distance between two distributions without performing a detailed comparison of the samples. In this paper, we apply this method to display our results in a synthetic manner (see Sects. 4.5-4.7) and only give additional details when necessary.

A36, page 33 of 33

Chapter 11

The Elusive CH +

In this chapter, we are addressing the problem of the abundance of CH + in the diffuse ISM. As explained in §4, the very presence of CH + is a long-standing conundrum of the physics of the ISM. What are the physical processes that control its abundance? Can the dispersion observed be explained by the joint action of the thermal instability and turbulence? Without taking into account the details of turbulent dissipation, can part of the observed column densities of CH + be reproduced? To try and answer those questions, we perform a parametric study ( §5), that we post-process chemically ( §13). The lines of sight are reconstructed with the algorithm presented in §8.3, so that the simulated lines of sight have the same distribution in lenghts as the observed sample. The distributions are finally compared through a modified version of the KS test ( §9).

Chemical Consideration

CH + formation The formation of CH + is a slow process because the routes of its formation have either very low rates, e.g. the radiative association Barinovs & Hemert (2006) C + + H → CH + +hν k C + ∼ 2 × 10 -17 T 300 K -0.42 cm 3 s -1 (11.1)

or they are very endothermic, e.g. the ion-neutral reaction (Hierl et al. 1997)

C + + H 2 → CH + + H k H 2 ∼ 7
.4 × 10 -10 exp(-4537/T ) cm 3 s -1 . (11.2) where 7.4 × 10 -10 cm 3 s -1 and E A /k B = 4537 K are the best fit parameters found by Hierl et al. (1997), and E A is the activation energy of this reaction. There are other formation routes for CH + , involving species that are orders of magnitude less abundant than the ones considered above. For this reason, they are usually secondary in the standard conditions of diffuse ISM, and therefore we do not discuss them here. Nevertheless, they are taken into account in the chemical solver. (11.5) and molecular hydrogen: (11.6) Because of its rapid destruction, CH + is a short-lived species, i.e. a species that exists exclusively where it forms because it cannot be transported.

CH

CH + + H 2 → CH + 2 + H k d,H 2 ∼ 1.2 × 10 -9 cm 3 s -1 ,

Impact of Warm-H 2

In each simulation, the abundance of CH + can be computed with the chemical solver presented in Appendix §13

As explained in §13, the chemical composition of the gas can be computed with a known density of molecular hydrogen, calculated out-of-equilibrium with RAMSES, or not. In this case, the density of H 2 is then calculated at equilibrium with the rest of the chemistry. In Fig. 11.1, we show the reconstructed column density of CH + vs the total column density N H (see §8.3), for those two cases applied to the fiducial simulation and we compare them to the observational dataset (black circles indicate detections and pink triangles indicate upper limits). When H 2 is computed out of equilibrium, the 2D probability histogram shows a bimodality, i.e. two peaks of probability, one at low column densities, and one at high column densities of CH + . The high-N (CH + ) regime extends over 1.5 orders of magnitude in N H and 2 orders of magnitude in fraction of CH + (dotted lines), and matches the bulk of detections (black circles). The bulk of non-detections, on the other hand, is located at the same total column densities as the peak at low-N (CH + ), between N H ∼ 2 × 10 19 and N H ∼ 9 × 10 20 cm -2 (see 4.7). While this low-N (CH + ) regime is at an almost constant fraction of CH + , f (CH + ) = N (CH + )/N H ∼ 10 -11 and is similar in the right and left panels, the difference between the location of the high-N (CH + ) regime is striking.

If the abundance of molecular hydrogen is post-processed at equilibrium (Fig. 11.1, right panel), the simulated lines of sight fail to reproduce the column densities of CH + by at least two orders of magnitude, and the bimodality in the distribution of N (CH + ) disappears. Since the output of the RAMSES simulation is the same in both panels, the structure of the matter is the same and, therefore, the PDF of the total column is the same. The difference on the distribution of N (CH + ) shown in Fig. 11.1 indicates that the observed column densities of CH + highly depend on the presence of out-of-equilibrium molecular hydrogen. Where does the CH + come from? In which physical condition is it produced and what chemical species are involved? the B x does not necessarily need to be lower than the equipartition to reproduce the observed 2D probability histograms of N (CH + ) vs N H .

Conclusion

The transport induced by large-scale turbulence in the multi-phase ISM is a viable scenario that explains simultaneously the large column densities of CH + detected in the diffuse medium, its trend as a function of N H , and its observed probability distribution function. The variations of the integrated fraction of CH + from one line of sight to another can be obtained with one set of parameters, i.e. without taking into account variations of large-scale physical conditions along one line of sight or from one line of sight to another. The tightest agreement is obtained for a neutral diffuse gas modeled over ∼ 200 pc, with a mean density n H = 2 cm -3 , illuminated by the standard interstellar UV radiation field G 0 = 1, and stirred up by a large-scale compressive turbulent forcing. This was also found to be the best set of parameters required to explain the properties of the HI-to-H 2 transition observed in the local ISM.

Combining the chemical and statistical information contained in the observations of CH + shows that out-of-equilibrium H 2 is capital in our study. CH + proves to be a tracer of the warm and out-of-equilibrium phase of the gas, warm-H 2 , and turbulence. A further improvement would be to take into account the excitation of the warm H 2 . Indeed, the excitation of H 2 would reduce the energy barrier of the formation of CH + , making it more effectively produced in warm gas. This is especially true in regions where the radiation field is high, i.e. where UV pumping effectively populates the high energy levels of molecular hydrogen.

Chapter 12

Excitation of Neutral Carbon

As explained in §4.2, Jenkins & Tripp (2011) derived the pressure through the observation of the levels populations of CI, by dividing the spectrum in velocity bins and computing the abundance of each level relative to the total amount of CI in each velocity channel. In this chapter, we explore the production and excitation of neutral carbon in the diffuse ISM in our simulations. Our goal is to analyze if the variation of pressure, density, and temperature due to the joint action of thermal equilibrium and turbulence are sufficient to reproduce the dispersions of f 1 and f 2 defined in eqs. 4.2. To achieve such a goal, we compute the abundance of CI with the chemical solver (see §13) and the population of the fine-structure levels, including all the excitation and de-excitation processes presented in §6. We reconstruct the lines of sight (see §8) with the same distribution of l los observed (see Fig. 12.1) and we then evaluate the integrated line profile (see §6.4) to have the population fraction of the two excited levels of the fine structure per velocity channel.

Distribution of lengths of LOS From the positions (see Fig. 4.4) and the distances (see Fig. 15.2) of the background sources observed by Jenkins & Tripp (2011), we compute the length of lines of sight with eq. 8.1. The normalized histogram (see Fig. 12.1) give us the weights for each l los to be used in the recontruction algorithm §8.

Reconstruction

We create a dataset of f 1 and f 2 which contain the N (CI * )/N (CI tot ) ratio and the N (CI * * )/N (CI tot ) ratio for each velocity channel of width 0.5 km s -1 . In Fig. 12.2, we show a comparison of the results with the raw simulation of l los = L = 200 pc (on the left panel) and taking into account the distribution of lengths of the lines of sight of the observational sample (right panel). The entire probability histogram appear to belong to the low-pressure regime of the f 2 vs f 1 graph (see Fig. 6.3). Moreover, the bulk of the probability histogram is at low values of f 1 and f 2. The peak of the reconstructed probability histogram (right panel) is found to be shifted and spread toward larger f 1 and f 2, compared to those of the initial distribution (left panel). This is particularly true for models with a higher Doppler parameter b (see Fig. 12.3). All these results were expected. Indeed, long lines of sight or large Doppler parameter favor the mixing in any velocity channel of gas at different physical conditions. This naturally favors larger f 1 and f 2 values compared to the original distribution.

transition, the C + -to-CI transition requires either larger sizes of CNM clouds (for a given density) or larger densities (for a given size). The size distribution of CNM clouds set by the combined action of turbulence and thermal instability implies that H 2 is built up in CNM clouds of typical sizes comprised between 3 and 10 pc (Bellomi et al. 2020). Because the size distribution is fixed, the mass of CI is necessarily built up in denser clouds which are considerably smaller than these scales. Increasing the resolution has two combined effects. It allows the gas to form denser entities which are favorable for the production of CI. Because these denser entities occupy a small fraction of the volume but also have a high surface filling factor, increasing the resolution increases the probability of intercepting them along random lines of sight. This is why the probability distributions shown in Fig. 12.5 shift towards larger densities as the resolution increases and why the 2D distribution of the f 1 and f 2 ratio sample larger values of f 1 and f 2 at larger resolution (Fig.

12.2).

Turbulence In Fig. 12.6, we show the effect of turbulence on the populations of CI, for turbulent strengths, F = 4.6 × 10 -4 , 1.5 × 10 -3 , and 4.6 × 10 -3 kpc Myr -2 , and compressibility factors, ζ = 0.1, 0.5, and 0.9. The corresponding distributions of n(CI * ) in a pressure-density diagram are displayed in Fig. 12.7. As shown by Seifried et al. (2011), Saury et al. (2014), andBellomi et al. (2020), increasing the turbulent forcing quickens the cycle between the CNM, LNM, and the WNM and favors the exchange of mass between these different phases. The mean pressure therefore increases, the 1D PDF of the density broadens and its bimodal nature progressively disappears (Piontek & Ostriker 2005;Walch et al. 2011). These broader variations of the pressure and the density naturally imply larger variations of the f 1 and f 2 ratios displayed in Fig. 12.6. Because high turbulent forcing drastically increases the mass and volume occupied by the LNM phase, the ratio f 2/ f 1 slightly increases compared to low and moderate turbulent forcing. The effect on the f 2/ f 1 ratio is only slight because, while the mass of CI produced in the LNM increases with the turbulent forcing, it is still low compared to the mass of CI built up in the CNM phase (see Fig. 12.7). Interestingly, the compressibility factor ζ has only a slight effect on the combined distribution of f 1 and f 2. This is surprising because ζ was found to have a large impact on the HI-to-H 2 transition (Bellomi et al. 2020). The main reason is that f 1 and f 2 respectively correspond to column density ratios N (CI * )/N (CI tot ) and N (CI * * )/N (CI tot ) per velocity channel and are therefore less sensitive to the total mass of CI.

Comparison with Observations

In Fig. 12.8, we compare the f 1 and f 2 ratios observed by Jenkins & Tripp (2011) to the results of the fiducial simulation at a resolution of 512 3 and all the other parameters set to their fiducial values. All the values of f 1 and f 2 computed in the simulation are found between the theoretical curves at low and high temperature T = 15 K and 10 4 K (dotted and dashed lines) and to be confined in the low-density regime, hence the low-pressure regime (see Fig. 6.3). Despite the variations of physical conditions induced by turbulence and the thermal instability, the predictions of simulations clearly fail to reproduce the distribution et al. 2014). This being said, the results of this chapter are unambiguous: the origin of the high-pressure component seen in the observations of CI is a deep and fundamental mystery that cannot be explained by thermal instability, turbulent transport between phases, or gravitational collapse. Several alternative scenarios need to be explored.

• The high-pressure component could originate from high-velocity J-type shocks propagating in the WNM. High-velocity J-type shocks (V S > 40 km s-1) propagating in a low-magnetized environment could bring the gas to pressure and density regimes high enough to significantly populate the second level of CI. In principle, the probability of occurrence of these events and their effect on the local chemistry and excitation processes should be captured by numerical simulations of the diffuse multiphase ISM. In practice, because the scales considered (∼ 0.4 -200 pc) are much larger than the viscous scale (∼ 0.3 mpc for the WNM) and because the numerical resolution is finite, such structures are smoothed out by numerical simulations. There is almost no way around that because simulations performed at a much smaller scale will fail to capture the thermal instability process. The only way to study these events would be to run 1D state-of-the-art models of interstellar shocks and ascribe a probability distribution function to such events based on the analysis of turbulence at larger scales.

• As proposed by Jenkins & Tripp (2011), the high pressure component could also be due to the intermittent dissipation of interstellar turbulence in CNM clouds (Falgarone et al. 1995;Joulain et al. 1998;Godard et al. 2009;Lesaffre et al. 2013;Godard et al. 2014;Lesaffre et al. 2020). Dissipative coherent structures in the CNM could locally increase the temperature (hence the pressure) of the gas to a level sufficient to excite CI. In particular, dissipation through ambipolar diffusion, i.e. the differential velocities between neutral and charged species, would have the effect of enhancing the collisional excitation rate of CI in collisions with charged species such as H + or e -.

In principle, such a scenario could explain the observations providing that a CNM cloud of density of ∼ 100 cm -3 is heated up to temperatures of about 1000 K or if the ambipolar diffusion drift is larger than ∼ 4 km s -1 . Again, the ideal MHD simulations presented in this work do not resolve dissipation processes and dissipative scales. Moreover, the simulations only consider a single fluid and therefore do not follow the decoupling between ions and neutrals. One way to study this effect would be to run state-of-the-art numerical simulations of the CNM designed to capture dissipative processes (Lesaffre et al. 2020, Richard et al., in prep.) including the decoupling between ions and neutrals.

• At last, and as proposed by Jenkins & Tripp (2011), the high-pressure component could be uncorrelated with the diffuse interstellar gas and rather originate from the gas associated with the background sources. If so, this component can not be captured by the present simulations which neglect the presence of stars and their surrounding environments.

The study of the diffuse interstellar medium is rich in terms of spectroscopic observations. In this Ph.D. work, we attempt not only to use the chemical information contained in the observations but also to take advantage of their statistical properties.

The local ISM is an open and complex system. One of the main challenges to be faced is to understand the joint action of the thermal instability and the turbulent motions, and the roles of different components of the ISM on observed quantities. To quantify separately the impacts of the mean density, the UV radiation field, the integral scale, the resolution, the turbulent forcing, the magnetic field, and the gravity on the observables, we used stateof-the-art MHD simulations, and an ensemble of 305 runs (see §5). Taking into account the distance of each background source and simulating random lines of sight over the same distribution of distances is paramount to explain the range of observed column densities and their corresponding statistics (see §8). To compare the results of the simulations to the observational sample, we propose a new version of the Kolmogorov-Smirnov test (see §9) which is a valuable tool for estimating the distance between two distributions and the errors on the estimated distance. Moreover, it provides information on where the maximal difference between the distributions occurs, bringing key insights into the origin of the discrepancies between observations and simulations. The original feature of this work is to not only focus on the production of individual column densities but also on their statistics, meaning the probabilities of occurrence of these column densities along random lines of sight. This exhaustive parametric study, combined with a chemical solver and the tools mentioned above, allows the investigation of the thermodynamical state and the evolution of the local diffuse interstellar medium. In this work, we applied the study to three chemical observables: H 2 , CH + , and CI.

The HI-to-H 2 transition The joint actions of thermal instability and large-scale turbulence in the fiducial simulation are found to produce a wealth of lines of sight which reproduce the observed position and width of the HI-to-H 2 transition, and whose probabilities of occurrence match those derived from the observations. The agreement is so remarkable that it is almost unnecessary to invoke variations of physical conditions along lines of sight or from one line of sight to another.

The distribution of column densities computed from the simulations strongly depends on the Galactic midplane density parametrized by the mean density n H , the density of OB stars parametrized by the UV scaling factor G 0 , and the scale of neutral diffuse clouds parametrized by the box size L. It is so because these three parameters not only regulate the mean pressure of the gas, hence the fractions of mass and volume occupied by the CNM and WNM, but also control the typical scale of density fluctuations in the WNM and the distribution of sizes of the CNM structures where H 2 forms. The tightest concordance between the observed and simulated samples is obtained for a mean density n H = 1 -2 cm -3 and a UV radiation field scaling factor G 0 = 1 (in Habing units), in good agreement with the values deduced from HI and CO all-sky surveys and from direct observations of the UV radiation field in the Solar Neighborhood. The range of observed column densities of H and H 2 requires a box size L = 200 pc which corresponds to the estimated scale of HI superclouds.

Within this setup, the column densities of HI are inferred to be built up in large-scale WNM and CNM structures correlated in densities over ∼ 20 pc and ∼ 10 pc, respectively. In contrast, H 2 is inferred to be built up at smaller scales. However, the fact that the kingfisher diagram is independent of the resolution of the simulation suggests that most of the mass and volume of H 2 is contained in CNM structures between ∼ 3 and ∼ 10 pc. All these values are given for the standard simulation (L = 200 pc and L drive ∼ 100 pc) but naturally depends on the size of the box and the mean density of the gas.

In spite of the strong influences of n H , G 0 , and L, the statistical properties of the HI-to-H 2 transition are otherwise remarkably stable. Admittedly, the kingfisher diagram depends on the strength of the turbulence if most of the forcing is injected in solenoidal modes; however, such a configuration prevents to reproduce the observational sample unless the large-scale velocity dispersion of the gas is unrealistically small. In contrast, if most of the kinetic energy is injected in compressive modes, the kingfisher diagram is found to weakly depend on the strength of the forcing. Similarly, the HI-to-H 2 transition is almost not affected by gravity and is found to weakly depend on the Doppler broadening parameter and the strength of the magnetic field, as long as B x 4 µG. The 2D probability histogram of the column densities of H and H 2 is, therefore, a valuable tool to constrain the nature of the turbulent forcing at large scales; however, it provides few or no information regarding the velocity dispersion of the gas, the amount of gravitationally bound environments and the strength of the magnetic field. Other observational tracers are required.

The elusive CH + The distribution of column densities N (CH + ) vs N H computed from the simulations not only strongly depends on the mean density n H , the UV scaling factor G 0 , and L, but also on the turbulent forcing, which greatly affects the width of the PDF of N (CH + ). As for the HI-to-H 2 transition, the observed variations of the integrated fraction of CH + from one line of sight to another can be obtained with a single simulation computed with one set of parameters, i.e. without taking into account variations of large-scale physical conditions along one line of sight or from one line of sight to another. The best agreement is obtained for a neutral diffuse gas modeled over ∼ 200 pc, with a mean density n H = 2 cm -3 , illuminated by the standard interstellar UV radiation field G 0 = 1, and stirred up by a large-scale compressive turbulent forcing. The study we conducted indicates that out-of-equilibrium H 2 is capital for reproducing the observed column densities of CH + .

The Neutral Carbon With the CI excitation code (see §6) we computed the populations of the fine-structure levels of neutral carbon and the fractions f 1 and f 2 per velocity channel, which is a measure of the pressure of the diffuse ISM. The probability histograms of f 2 vs f 1 are greatly affected by the resolution, the turbulent forcing strength, and the doppler parameter. Contrary to the case of CH + and molecular hydrogen, the distribution of pressure derived by Jenkins & Tripp (2011) cannot be fully reproduced by the combination of the thermal instability and the turbulent transport between the phases, or gravitational collapse. While the distribution of pressure in the low-pressure regime is in agreement with that deduced from the observed populations of the first two excited fine-structure levels of CI, the high-pressure component is not reproduced. These discrepancies require the exploration of other scenarios such as high-velocity J-type shocks, or intermittent dis- Table 13.1: The elements included in the chemical code, and their initial state (neutral or ionic). The elemental abundances take into account the solar abundances and a depletion factor. On the third column the first ionization potential of the neutral species in eV. Species with ionization potential lower than the Lyman limit (13.6 eV) are in initially their ionic state.

Reactions

The chemical solver takes as inputs not only the physical condition of each cell of the simulation but also an input file describing a chemical network, i.e. the γ i , α i , β i , E i for each of the reactions i included in the code. This network is used to create the system of equations which is then solved using the Newton-Raphson algorithm to compute the densities of all 151 species.

γ i , α i , β i , E i are parameters given by quantum mechanical calculations or experiments that allow to compute the reaction rates k i . Chemical reaction rates, k i , for two-body collisions (e.g. ion-neutral, dissociative recombination, radiative association) are modeled with an Arrhenius-Kooij's law:

k i = γ i T 300 α i exp (-E i /k B T ) (13.1)
an expressed in cm 3 s -1 , where T is the effective temperature of the collision between the two reactants. In the above formuma, E i is the activation energy for endothermic reactions, while γ i and α i are the coefficients setting the amplitude and temperature dependence of the reaction rate. The total reaction rate between two species j and k with densities n j and n k in cm -3 , reads k i n j n k , and is expressed in cm -3 s -1 .

Cosmic-ray processes, both direct or that include secondary photons are computed as

k CR = γ i ζ H 2 (13.2)
where k CR is the reaction rate in s -1 and ζ H 2 is the cosmic ray dissociation rate of molecular hydrogen, and γ i is the scaling parameter for each atom and molecule compared to H 2 .

Photoreaction rates are finally computed as

k γ = γ i G 0 exp -β i A V (13.3)
where k γ is the rate in s -1 , G 0 the scaling factor of the standard ISRF and β i describes the dependence of the rate on the visual extinction A V .

relation seems to appear between ζ H 2 and N H , but a drop happens once the cloud goes from diffuse to dense and the low-energy CR are not energetic enough to penetrate into the cloud and "stop" in the outer region of the cloud.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Interstellar continuum radiation field in an HI cloud in the Solar Neighborhood as presented in Draine's book. Spectral lines are not included. Solid line is the sum of all components for hν ≤ 13.6 eV.

Figure 4 . 2 :

 42 Figure 4.2: H 2 column density as a function of the total column density of protons N H . Open circles correspond to detections of H 2 while arrows correspond to upper limits. The blue dashed line indicates the maximum value of N(H 2 ) derived from a purely molecular medium with an integrated molecular fraction

Figure 4 . 3 :

 43 Figure 4.3: Normalized intensity observed by Jenkins & Tripp (2001) exhibiting multiple lines from different fine-structure levels of CI. The ground state is indicated with CI, and CI * and CI * * are the first and second excited levels, respectively. The original data normalized to a continuum is shown in black for the multiplets in the spectrum of HD 210839. The best fit model is shown in red.

Figure 4 . 5 :

 45 Figure 4.5: From Jenkins & Tripp (2011): Measurements of f 1 and f 2 for all velocity bins (each having a width of 0.5 km s -1 ) for the 89 background stars in the survey. The area of each dot is proportional to the value of N (CI tot ), with a normalization in size indicated on the top portion of the plot. The white X located at f 1 = 0.209, f 2 = 0.068 represents the "center of mass" of all of the dots. The curves indicate the expected level populations for three different temperatures, for a CNM cloud with a molecular fraction f (H 2 ) = 0.6 and a fraction of helium x(He) = 0.094 (more detail on the computation of the theoretical curve is given in §6). The large open circles on the curves indicate integer values of log(P /k B ) with accompanying numbers to indicate their values.

Figure 4 . 8 :

 48 Figure 4.8: From Godard et al. (2014): Observations (blue and red open circles) are compared to the predictions of PDR (top panel) and TDR (bottom panel) models. The CH + and total hydrogen column densities inferred from visible absorption lines and extinction measurements are shown in red and are fromCrane et al. (1995);Gredel (1997);Weselak et al. (2008). The CH + and total hydrogen column densities inferred from the far-infrared rotational transitions of CH + observed with Herschel/HIFI instrument and the 21 cm emission lines of HI are shown in blue and are fromFalgarone et al. (2010), andGodard et al. (2012). The PDR model predictions (filled symbols, top panel) are computed for several densities: 10 (crosses), 30 (triangles), 50 (squares), and 100 cm -3 (circles). The TDR models (filled symbols, bottom panel) are computed for two different dissipation rates ǫ = 10 -24 (green) and 10 -23 (magenta) erg cm -3 s -1 , seeGodard et al. 2014).

Figure 5 . 10 :

 510 Figure 5.10: Turbulent velocity dispersion along one axis as a function of the correlation time, for different box physical size and driving strength. The two horizontal lines indicate the range of the observed turbulent velocity dispersion. The driving strength in this graph is in machine units, that can be converted into physical units by the factor 1.5 × 10 -7 kpc Myr -2 .

Figure 6 . 1 :

 61 Figure 6.1: CI fine-structure energy levels in ground state. This splitting is due to LS coupling.

Figure 7 . 1 :

 71 Figure 7.1: Vibration-rotation energy levels of the ground electronic state of molecular hydrogen with J ≤ 29 taken fromDraine (2011) with indication of the 1 -0S(1) transition, taken as an example. 1 -0S(1) transition refers to v = 1 → 0 and ∆J = 2 ending in J = 1 so the rotational transition J = 3 → 1. This transition occurs at 2.1218 µm

  Fig.3. H 2 column density as a function of the total column density of protons N H . Open circles correspond to detections of H 2 while arrows correspond to upper limits (see TableA.1). The blue dashed line indicates the maximum value of N(H 2 ) derived from a purely molecular medium with an integrated molecular fraction f H 2 = 1 (Eq. (2)). The red dashed-dotted line indicates the theoretical molecular fraction derived in an unshielded WNM-type environment with a density of 0.5 cm -3 and a temperature of 8000 K, illuminated by a UV photon flux of 10 8 cm -2 s -1 (see Eqs. (13) and (15)). The regions A, B, C, D, and E defined in Table1correspond to an arbitrary separation of the observational sample used for quantitative comparisons with the results of simulations (see Sect. 4).

  E.Bellomi et al.: 3D chemical structure of diffuse turbulent ISM. I.

Fig. 4 .

 4 Fig. 4. Colored tables of the mean pressure expressed in K cm -3 (first line), the turbulent velocity dispersion σ tur (second line), and the fractions of mass f WNM and f CNM contained in the WNM phase (third line) and the CNM phase (fourth line). First and second columns: these quantities as functions of n H and G 0 , for L = 50 pc and F = 3.6 × 10 -3 kpc Myr -2 (first column) and for L = 200 pc and F = 9 × 10 -4 kpc Myr -2 (second column). Third and fourth columns: these quantities as functions of the acceleration parameter F and the compressive ratio ζ, for L = 50 pc (third column) and L = 200 pc (fourth column). All other parameters are set to their standard values (see Table2).

Fig. 10 .

 10 Fig. 10. Same as in Fig. 9 for simulations with a box size of L = 50 pc instead of 200 pc. The turbulent forcing is adjusted as in Fig. 4 (F = 3.6 × 10 -3 kpc Myr -2 ) to obtain similar velocity dispersions for L = 50 pc and L = 200 pc.

  Fig. C.1. Probability histogram of the proton density n H extracted from the fiducial simulation (see Table2). The black histogram correspond to the extracted data. The red dashed curve shows an example of the sum of two log-normal components and a power-law tail at high density, for comparison. The blue line indicates the inflection point of the PH between the diffuse and the dense components.

  Fig. C.2. Comparisons of the 1D probability histograms of the total column density normalized to the integration scale, N H /l, extracted from the simulations (black histograms) and constructed with the semi-analytical model described in the main text (green histograms) for l = 200 (top left), 100 (top right), 50 (bottom left), and 25 (bottom right) pc. Each of these four main panels displays the comparisons performed for 15 simulations with different n H and G 0 around the fiducial setup defined by n H = 2 cm -3 and G 0 = 1. All probability histograms inferred from the semi-analytical model are obtained assuming a fixed correlation length of the diffuse component y diff dec = 0.2 × L drive = 20 pc and a correlation length of the dense component y dens dec = 10 × (n H /2cm -3 ) 1/3 pc (see main text).

A36

  Figs. C.2 and C.3 for samples of N = 64 000

Fig. C. 3 .

 3 Fig. C.3.Comparisons of the 2D probability histograms of the total column density normalized to the integration scale, N H /l, and the column density of H 2 normalized to integration scale, N(H 2 )/l extracted from the simulations and constructed with the semi-analytical model described in the main text for l = 200 (top left), 100 (top right), 50 (bottom left), and 25 (bottom right) pc. The results of the simulations are indicated with contour plots with isoprobabilities of 10 -4 , 10 -3 , and 10 -2 . The colored histograms correspond to the results obtained with the semi-analytical model. As in Fig.C.2, each of the four main panels displays the comparisons performed for 15 simulations with different n H and G 0 around the fiducial setup defined by n H = 2 cm -3 and G 0 = 1. All probability histograms inferred from the semi-analytical model are obtained assuming a fixed correlation length of the diffuse component y diff dec = 0.2 × L drive = 20 pc and a correlation length of the dense component y dens dec = 10 × (n H /2 cm -3 ) 1/3 pc (see main text).

Fig. C. 4 .

 4 Fig. C.4. Same as Fig. C.3, assuming that y dens dec = 10 × (n H /2 cm -3 ) 1/3 (n H /n lim H ) -0.35 pc. The exponent used here corresponds to a simple test to show the effect of a dense decorrelation scale y dens dec varying with density.

  D.1).

  + destructions In the diffuse ISM, the principal destruction routes for CH + are the photodissociation by UV photons (van Dishoeck & Black 1988) CH + + hν → C + H + k γ ∼ 2.5 × 10 -10 G eff s -1 , (11.3) the dissociative recombination with electrons (Florescu-Mitchell & Mitchell 2006) CH + + e -→ C + H k d,e -∼ 1.5 × ×10 --neutral reaction collision of CH + with atomic hydrogen, CH + + H → C + + H 2 k d,H ∼ 7.5 × 10 -10 cm 3 s -1 ,

Table 3 . 1 :

 31 Comparison of the UV radiation field integrated between 910 and 1110 Å, as described by

Table 7 . 1 :

 71 Notation for the rotational transitions Electronic potential well as function of the internuclei distance (in Å) with energies written in eV on the extreme right, in cm -1 on the left, where also the transition wavelength in Å is indicated. Beside every well there is the term symbol indicating the electronic quantum state. Ground state is indicated by the notation X 1 Σ + g , the first electronic excited state, i.e. the Lyman band by B 1 Σ + u , and the second excited electronic state, i.e. Werner band by C 1 Π u . Between the ground state and the first electronic excited state there is an unbound state b 3 Σ + u .

	Figure 7.2:

Table 1 .

 1 Statistical properties of H 2 observations in the subsamples A, B, C, D, and E defined in footnote and shown in Fig.3.Notes. Only the lines of sight where H 2 has been detected (312 sources out of 360, see TableA.1) are considered. The mean and dispersion values, µ and σ are computed on the logarithm of the molecular fraction f H 2 observed in the corresponding subsample. Definitions of subsamples: region A: 5 × 10 18 N H 5 × 10 19 , 5.6 × 10 -8 f H 2 10 -3 ; region B: 5 × 10 19 N H 7 × 10 20 , 5.6 × 10 -8 f H 2 10 -3 ; region C: 1.6 × 10 20 N H 10 21 , 10 -3 N H

	Region Number % µ log( f H 2 )	σ log( f H 2 )
	A B C D E	10 41 13 3 50 16 204 65 0 0	-5.02 -4.62 -1.10 -0.70 -	0.99 0.70 0.49 0.37 -
	10 21 N H 10 22 , 10 -3 N H 10 20 10 -3 N H 10 20 0.9 f H 2 1.	0.9	0.9 f H 2 1; region E: 10 22 N H 10 23 , f H 2 1; region D: 10 20

Table 2 .

 2 Fiducial model and range of parameters explored in this work. here through several grids of simulationsincluding a total of 305 runscovering a broad range of physical conditions and centered around a fiducial setup 3 . The reference value adopted for each parameter, and the range of values explored in this work, are summarized in Table2. Among all parameters, L, n H and G 0 are of particular importance.

	Parameter	Notation	Ref	Range	Units
	Box size	L	200	20-200	pc
	Mean density UV radiation field Resolution	n H G 0 R	2 1 256 3	0.5-4 0.5-4 64 3 -512 3	cm -3 Habing field -
	Turbulent forcing Compressive ratio	F ζ	9 × 10 -4 10 -5 -10 -2 0.1 0.1-0.9	kpc Myr -2 -
	Doppler broadening Initial magnetic field Self-gravity Galactic well	b D B x --	8 3.8 on on	1-8 0-40 on-off on-off	km s -1 µG --
	transition is studied

Table A .

 A 1. continued.

	Source ID	Longitude [ • ] Latitude [ • ] Distance [kpc] E(B-V) log 10 (N(H)) log 10 (N(H 2 )) log 10 (N H )
	HD	132.910	-5.870	2.809	0.410

  A36, page 22 of 33 E. Bellomi et al.: 3D chemical structure of diffuse turbulent ISM. I.

	Table A.1. continued.
	Source ID	Longitude [ • ] Latitude [ • ] Distance [kpc] E(B-V) log 10 (N(H)) log 10 (N(H 2 )) log 10 (N H )
	HD	287.570	-0.710	2.688	0.370

  A36, page 24 of 33 E. Bellomi et al.: 3D chemical structure of diffuse turbulent ISM. I.

	Table A.1. continued.
	Source ID	Longitude [ • ] Latitude [ • ] Distance [kpc] E(B-V) log 10 (N(H)) log 10 (N(H 2 )) log 10 (N H )
	HD HD	10.760 10.460	-1.580 -1.740	2.079 0.861	0.310 (3) 0.330	21.08 (3)	20.18 (2)	21.18

  The large-scale self-shielding not only increases the fraction of H 2 in atomic slabs (i.e., for n H < n tr H ) but also shifts n tr H toward lower values, two effects which depend on the size of the box L. For L = 200 pc, including the large-scale selfshielding is found to reduce n tr H by a factor of two.

A36, page 30 of 33 E. Bellomi et al.: 3D chemical structure of diffuse turbulent ISM. I.

Soft X-ray have typical energies below 1 keV

Fine-structure transitions are 10 14 times weaker than Lyman α because (1) fine-structure transitions are magnetic dipole transitions, therefore their transition matrix elements are of the order of α

F S ∼ 5 × 10 -5 , where α F S = 1/137 is the fine structure constant, thus 10 5 times smaller than a permitted transition; (2) the splitting between finestructure levels is smaller and since the A i j ∝

ν 3 i j , the Einstein coefficients of the fine-structure transitions are 10 9 smaller than Lyman α Einstein coefficients.

V it the visible passband at λ ∼ 550 nm, and B the blue passband at λ ∼ 440 nm

Available from the NED (http://nedwww.ipac.caltech.edu) and on the more recent plateform https://irsa.ipac.caltech.edu.

The continuous viewing zone, CVZ, exclude the zone where Earth occultations happen

As a comparison, the meam free path in air is λ air ∼ 1 µm 1 AU ∼ 1, 496 × 10 13 cm

The range depends on the geometry of the system

These MHD equations are written in the code in a condensed form, using the definition of the tensor product as A = B v, where A is a second rank tensor, also called dyatic tensor, A i j = B i v j . Its flux is ∇ • A = ∂ i A i j .

In this work we use the relation between N H and A V deduced from the observations of the mean Galactic extinction curve(Fitzpatrick & Massa 1986).

We note that for ζ = 0.5, the power of the compressive forcing corresponds to 1/3 of the total power.

The grids have been run on the computing cluster Totoro funded by the ERC Advanced Grant MIST. The computational time of the fiducial simulation is ∼ 6 days using 40 processors.

We remind that in astronomy all elements heavier than helium are called metals

The grids have been run on the computing cluster Totoro funded by the ERC Advanced Grant MIST. The computational time of the standard simulation is ∼6 days using

processors.

Notes. The distance of each source is computed from the parallax measured by Gaia if the data is given in the DR2 catalog (Gaia Collaboration 2018); otherwise, the distance of the source is taken fromGudennavar et al. (2012). Column densities are expressed in cm -2 . The total proton column densities N H are computed as N(H) + N(H 2 ) if the column densities of HI and H 2 are available, or derived from the reddening E(B -V) as N H = 5.8 × 10 21 E(B -V) cm -2 assuming a standard Galactic extinction curve and the average interstellar ratioR V = A V /E(B -V) = 3.1(Fitzpatrick & Massa 1986; Fitzpatrick 1999). References.(1) York (1976);(2) Savage et al. (1977);(3) Bohlin et al. (1978);(4) Hobbs (1978);(5) Neckel & Klare (1980);(6) Federman (1982);(7) Bohlin et al. (1983);(8) Shull & van Steenberg (1985);(9) Savage et al. (1985);(10) Jenkins et al. (1986);(11) van Steenberg & Shull (1988);(12) van Dishoeck & Black (1989);(13) Welsh et al. (1990);(14) Welsh et al. (1991);(15) Garmany & Stencel (1992);(16) Palazzi et al. (1992);(17) Federman et al. (1994);(18) Diplas & Savage (1994);(19) Fruscione et al. (1994);(20) Roth & Blades (1995);(21) Schlegel et al. (1998);(22) Wolff et al. (1999);(23) Le Coupanec et al. (1999);(24) Linsky et al. (2000);(25) Ryu et al. (2000);(26) Savage et al. (2001);(27) Rachford et al. (2002);(28) [START_REF] Ambartsumian | The evolution of stars and astrophysics[END_REF];(29) André et al. (2003);(30) Lehner et al. (2003);(31) Cartledge et al. (2003);(32) Cartledge et al. (2004);(33) Oegerle et al. (2005);(34) Pan et al. (2005);(35) Gnacinski & Krogulec (2006);(36) Gillmon et al. (2006);(37) Ritchey et al. (2006);(38) Snow et al. (2007);(39) Burgh et al. (2007);(40) Sheffer et al. (2007);(41) Jensen & Snow (2007a);(42) Cartledge et al. (2008);(43) Sheffer et al. (2008);(44) Snow et al. (2008);(45) Bowen et al. (2008);(46) Rachford et al. (2009);(47) Burgh et al. (2010);(48) Gudennavar et al. (2012).

Available from the NED (http://nedwww.ipac.caltech.edu) and on the more recent plateform https://irsa.ipac.caltech. edu.

Acknowledgements

Acknowledgements. We thank the referee for his careful reading and his valuable comments and suggestions. The research leading to these results has received fundings from the European Research Council, under the European Community's Seventh framework Programme, through the Advanced Grant MIST (FP7/2017-2022, No 742719). We would also like to acknowledge the support from the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. We would finally like to thank E. Falgarone and M. Gerin for the stimulating discussions we had and their precious comments regarding this work.

A36, page 19 of 33 A&A 643, A36 (2020)

Part III

Comparison with Observations

Appendix A: Observations of HI and H 2 in the local diffuse ISM

The complete observational dataset used in this work is given in Table A.1, including the sources identifiers, the coordinates, and the column densities of HI and H 2 . The column densities of H 2 are obtained through direct observations of its FUV absorption lines. The sample presented in Table A.1 therefore results from a combination of the Copernicus survey of nearby stars (e.g., Savage et al. 1977;Bohlin et al. 1983) and the FUSE surveys of the Galactic disk and the Galactic halo (e.g., Rachford et al. 2002Rachford et al. , 2009;;Lehner et al. 2003;Cartledge et al. 2004;Pan et al. 2004;Gillmon et al. 2006;Jensen & Snow 2007a,b). In contrast, and as reviewed by Fruscione et al. (1994), the HI column densities are derived from both direct methods, which include Lyα absorption studies, EUV observations of stellar spectra, and observations of the 21cm line, and indirect methods, which include curve of growth of neutral and singly ionized atoms and optical interstellar absorption lines of NaI which are both found to correlate with N(H) (de Boer et al. 1986;Ferlet et al. 1985).

The column densities of HI given in Table A.1 therefore result from a combination of FUV and optical studies for nearby stars (e.g., Bohlin et al. 1978Bohlin et al. , 1983;;Diplas & Savage 1994;Fitzpatrick & Massa 1990;Jensen & Snow 2007b,a) and radio studies of the 21cm line for extragalactic sources at high Galactic latitude (Wakker et al. 2003;Gillmon et al. 2006).

As reported by all these authors, the indirect methods relying on the observations of metals can be subjects to uncertainties mostly due to the assumptions made regarding the elemental abundances. In addition, the column densities of HI derived from Finally, the cooling rate due to electron recombination onto charged grains and PAHs is set to Λ rec = 4.65 × 10 -30 T 0.94 κ β x e φ PAH erg cm 3 s -1 , (B.9)

with β = 0.74/T 0.068 . To validate the calculations of the heating and cooling terms, we compare in Fig. B.1 the thermal equilibrium curve obtained with RAMSES to the predictions of Wolfire et al. (2003).

Appendix C: Analytical description of 1D and 2D probability histograms

In order to interpret the results found in Sect. 4, we propose a semi-analytical prescription to predict the 1D and 2D PHs of the total column density and the column density of H 2 obtained with numerical simulations. This prescription is based on the work of Vázquez-Semadeni & García ( 2001) and Bialy et al. (2017Bialy et al. ( , 2019) ) who showed that lines of sight across isothermal simulations of turbulence can be accurately modeled as a series of random density fluctuations.

C.1. Decorrelation scales

In Fig. C.1, we display an example of the volume-weighted distribution of the proton density computed for the fiducial simulation (see Table 2). The total column density integrated along the x direction over a random line of sight of size l L is

Part IV

Conclusion and Perspectives

sipation of interstellar turbulence.

Perspectives

All these results open new perspectives for the study of the chemical state of the ISM in which any observation must be understood through the combination of local physical conditions and the probabilistic ordering of these conditions along the line of sight. In particular, similar studies should be applied to all atomic and molecular species with observational samples large enough to conduct statistical analysis. Since the modifications of the Kolmogorov-Smirnov test make it applicable to N dimensions, the study can be expended to probability histograms in higher dimensions, taking into account simultaneously the joint information contained in the column densities of several species.

Part V Appendix

189

Chapter 13

Postprocessing: Chemistry

To explore the effects of the local physical conditions (density, temperature, and local radiation field) and the global properties of the ISM (turbulence, magnetic field, mean radiation field, thermal instability, and gravity) on observational tracers such as CH + and the level population of CI, we need to compute the chemical composition of the multiphase medium. To achieve such a goal, we post-process the simulation output with a state-ofthe-art chemical solver that computes the density of 151 species in each cell of the simulation, by solving a chemical network, at equilibrium, that includes around 2700 reactions.

The chemical network includes gas-phase reactions (e.g., neutral-neutral, ion-neutral), reactions onto grains (e.g., the formation of H 2 ), photoreactions due to UV photons (photodissociation and photoionization), and reaction with cosmic rays. The code takes as input the local conditions of a fluid cell, i.e. the temperature T , the proton density n H , the radiation field G 0 , the extinction by dust grains A V , the self-shielding of H 2 , f shield , and optionally the abundance of H 2 . The gas-phase chemical abundances are taken to be solar abundance with a depletion factor that depends on the element (see Table 13.1).

The chemical solver, fully described by Valdivia et al. (2017) (see Appendix A) computes the abundances of all chemical species at equilibrium assuming two different prescriptions depending on the choice of the user. If the out-of-equilibrium of H 2 computed in RAMSES is given in input, the code computes the chemistry at equilibrium except for H and H 2 . If not, the code computes the entire chemistry at equilibrium, including the abundances of H and H 2 . Such a method allows to compute the chemistry partly outof-equilibrium and is justified by consideration on chemical timescales. Indeed as found by Valdivia et al. (2017) the equilibrium timescales of H 2 are larger than that of any other species over a wide range of physical parameters, as long as n H > 3 cm -3 or f shield < 10 -2 . Although these conditions do not cover the entire physical conditions of the simulations, it covers most of the physical conditions over which the chemistry is activated by the formation of H 2 . The chemical solver combined with the separated prescriptions for the treatment of H 2 not only provides the entire chemical composition of the diffuse neutral matter but also allows to estimate the impact of out-of-equilibrium molecular hydrogen on the chemistry of the ISM.

Chapter 14

Cosmic-Rays Ionization Rate

Cosmic rays are high-energy particles presented in §1.3.4. Since low-energy cosmic rays are not observed directly, their ionization rate needs to be deduced indirectly, constrained by other chemical observers, e.g. H + 3 , which is the most direct tracer of CR ionization, due to its rather simple chemistry. The formation of H + 3 happens in two steps: CR photodissociate H 2 creating a H + 2 ion that through a neutral-ion reaction with another H 2 creates H + 3 . The H + 2 formed by the first reaction can also be destroyed by dissociative recombination with electrons or charge transfer to atomic hydrogen, but both these processes are slow in comparison with the reaction that forms H + 3 . The first step of the formation of H + 3 is called the rate-limiting step because it is very slow compared to the second, and so its reaction rate can be taken as the formation rate of H + 3 . For the destruction of H + 3 the predominant process in the diffuse medium is the dissociative recombination with an electron, with a rate of k e . If we assume steady state (0), the formation of H + 3 it is equal to its destruction:

where k(H + 3 |e -) ∼ 2 × 10 -7 cm 3 s -1 is the dissociative recombination rate which measured experimentally by McCall et al. (2004), N (H + 3 ) and N (H 2 ) are the column densities of H 2 and H + 3 which can be observed. Equation (1) assumes that the ratio between the density of H + 3 and H 2 does not change much on a line of sight and equation (2) takes into account that C + is the main source of electrons because it has an ionization potential lower than Lyman limit and it is the most abundant element after hydrogen (Sofia et al. 2004). Indriolo & McCall (2012) computed this calculation also for a more complete and complex chemical reaction networks taking into account the destruction of H + 3 by CO and O, showing that in diffuse clouds the previous prescription gives good results. McCall et al. (2003) found a CR primary ionization ratio for atomic hydrogen ζ H ∼ 5 × 10 -16 along sightline to ζ Persei. Indriolo & McCall (2012) found that the mean value of the ionization rate is ζ H 2 ∼ 3.5 × 10 -16 s -1 with big variations around this mean value, values can vary greatly from one line to another from 10 -17 to 10 -15 s -1 , apparently due to small-scale and local effects, e.g. the distance between the acceleration the ionization site. ζ H 2 was expected to decrease with increasing column density N H (Padovani et al. 2009) because CR would lose their energy entering the cloud and ionizing more and more matter. However, no cor-