Aknowledgements

First of all, I would like to thank Pierre Lesaffre, Edith Falgarone and Guillaume Pineau des Forêts for their invaluable help and support during the writing of this thesis. I am also grateful to Maryvone Gerin, director of the Laboratoire de Radioastronomie and Michel Pérault, director of LERMA and my colleagues in the laboratory for providing an excellent working environment. I would also like to thank my family: my father Kostas and my sisters Amalia, Myrto and Eleni for their moral support. Next, I would like to thank Nicolas and Nadia Behr for their invaluable support during the most difficult parts of the writing of the thesis. Last but not least, I would like to thank my good friends Tarina Laspopoulou and Cédric Goubert for being there for me. 7

Abstract

Turbulent energy dissipation is a key process in the cold interstellar medium (ISM), not only on the road to star formation but also as a source of suprathermal energy able to open new chemical routes, otherwise inactive at the low gas temperature. Such routes are required, though, to explain the high abundance of species such as CH + and SH + observed in the ISM. In this context, the space-time intermittency of energy dissipation is particularly relevant because it drives injection of suprathermal energy in the ISM locally far above the average level. The detailed characteristics of the spatial distribution and the geometrical properties of the energy dissipation rate can provide valuable inputs to chemical models.

We study them here with the aid of direct numerical simulations with unprecedented dedication to resolve the dissipation processes numerically. As the sonic Mach number in the diffuse interstellar medium can take values in a wide range, we bracket the possible physics by considering two categories of decaying turbulence models. On the incompressible extreme, we perform pseudo-spectral simulations of viscous and resistive magnetohydrodynamics, with a particular emphasis on ambipolar diffusion due the ion-neutral drift. On the compressible extreme we consider grid-based (Godunov) simulations of isothermal resisitive and viscous magnetohydrodynamics where our focus is on numerical dissipation.

Our incompressible simulations show that ambipolar diffusion leads to force-free magnetic fields at small scales. As a result, the typical scale of ion-neutral friction heating is displaced to large scales in the inertial range, much greater than dimensional analysis would predict. The structures of high dissipation are spatially coherent sheets, each with a single nature of dissipation (viscous, ohmic or ambipolar). We reveal their statistical scaling laws and compute their intermittency exponents.

We show that compressible simulations are subject to a lot of numerical dissipation: in our set up, less than half of the total dissipation is accounted for by the physical terms, the rest is produced by the numerical scheme. We design a method to recover locally the energy lost in the scheme and we use it to examine the sheet-like structure of the dissipation field as in our incompressible simulations. We show that numerical dissipation prevents us to assess the nature of dissipative structures. For instance, although we confirm previous results that physical dissipation in shearing sheets rapidly dominates over shocks, the balance could be reversed if numerical dissipation were shock dominated.

Finally, we examine the efficiency of various observational tracers to characterize the structures of high dissipation. In particular, we find that increments of molecular line centroid velocity or of polarization Stokes parameters correlate very well on the plane of the sky with specific structures of high dissipation. We also compute the intermittency exponents measured for these tracers and find they span a broad range of possible values. At last, we mix the Fourier phases to demonstrate the crucial role of coherence in producing the filamentary structure of observable maps of increments such as recently produced by the Planck collaboration. A.1 Adjusted viscosity of shocks with constant ν, constant µ or no viscosity at all. In the case without viscosity, it is the adjusted constant ν f which is given, as it fits the shock profile much better. The percentages in parenthesis give the difference between the total irreversible heating as recovered from the method 1 and the theoretical one. . . . . . . . . . . . . . . . . isotropic turbulence generated behind a grid [START_REF] Corke | A new view on origin, role and manipulation of large scales in turbulent boundary layers[END_REF]. Top right: Phytoplankton bloom reveals turbulent ocean currents, an example of geophysical turbulence [START_REF] Signorini | Environmental factors controlling the barents sea spring-summer phytoplankton blooms[END_REF]. Bottom left: Turbulence in the interstellar medium as revealed by the gradient of linear polarization [START_REF] Gaensler | Low-mach-number turbulence in interstellar gas revealed by radio polarization gradients[END_REF]. Bottom right: Current density magnitude in a high resolution numerical experiment of homogeneous isotropic magnetohydrodynamic turbulence [G. Momferratos, unpublished]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Streamwise velocity component recorded by a hot-wire anemometer in a classical wind tunnel experiment of homogeneous isotropic turbulence behind a grid [START_REF] Sreenivasan | Fractals and multifractals in fluid turbulence[END_REF]. . . . . . . 2.3 Schematic representation of the energy cascade process (from [START_REF] Ecke | The turbulence problem -an experimentalists perspective[END_REF]). . . . . . . . . . . 2.4 Viscous dissipation rate along a fixed line from a 1024 3 incompressible magnetohydrodynamics simulation with initial condition based on the ABC flow. The snapshot shown is at the temporal peak of total (ohmic + viscous) dissipation [G. Momferratos, unpublished]. . . . . . 2.5 A realization of one-dimensional Brownian motion. "Wiener process zoom". Licensed under Creative Commons Attribution-Share Alike 3.0 via Wikimedia Commons. . . . . . . . . . . . 2.6 Pdf of the longitudinal velocity increment from a high resolution incompressible hydrodynamical simulation for various separations r n = 2 n ∆x, ∆x being the distance between grid points in the simulations. (from [START_REF] Ishihara | Study of high-reynolds number isotropic turbulence by direct numerical simulation[END_REF]). . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Comparison of experimental data with the predictions of the structure function scaling exponents of various intermittency models (from [START_REF] Anselmet | Turbulent flows and intermittency in laboratory experiments[END_REF]). . . . . . . . . . . . . . . B y . Red : v x . Cyan: v y . To be compared with [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF], figure 1a 4.1 Compensated kinetic (top) and magnetic (bottom) energy spectra for OT runs 10, 12 and 14 at the temporal peak of dissipation. The assumed limits of the inertial and dissipation ranges are shown in dashed vertical lines. The vertical blue and red dashed-dotted lines with square symbols correspond to the expected AD critical wavenumber k a (see equation (4.5)) in runs 12 and 14, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Time evolution of volume integrated dissipation rates for the OT runs 10, 12 and 14. Solid lines show the Ohmic plus viscous dissipation which we use to define peak dissipation. . . . 4.3 Log-normal core and power-law tail fit for the pdf of the total dissipation (Run 12 AD -OT).

List of Tables

List of Figures

The core follows a log-normal distribution with mean µ l -4.27 and standard deviation σ l 1.03 while the tail follows a power-law with exponent -τ -2.61. Vertical lines show the mean value (red) and thresholds located at 1 (black), 2 (green) and 3 (blue) standard deviations above the mean value. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Cumulative probability density function of the total dissipation for run 12 (AD -OT). . . . . 4.5 Compensated dissipation spectra for run 12 (AD-OT, Re a = 100). Blue solid line: AD dissipation, red dashed line: Ohmic dissipation, green dotted line: viscous dissipation. We plot ke (k) in a log-lin plot, so that the area under the curve over any interval shows directly the amount of power inside this interval. We mark the position of the maximum value of of Run 12 (AD -OT). They are defined as connected sets of points having values of the total dissipation three standard deviations above the mean value. Each little sphere has a 2-pixels diameter, ie: about the size of the viscous (or equivalently Ohmic) dissipation length. . . . . 4.16 One of the largest structures extracted from run 12 . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Scaling relations P i ∝ L D P i from Run 12 (AD -OT), at the peak of dissipation, with a threshold of two standard deviations above mean value. The dotted line shows the effect of adopting the slope found by [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF] 4.19 We look at the subset of pixels above a given threshold of total dissipation in run 12 (AD-OT, Re a = 100) at the dissipation peak. For each value of the threshold, we plot the fraction of the total energy dissipation on this subset versus the volume of this subset (black curve). We also give the fraction of the total dissipation on this subset for each nature of dissipation (red: Ohmic, green: viscous, blue: AD). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Comparison of scaling exponents with three-sigma error bars between pure MHD (red circles) and AD MHD (green squares Re a = 100 and blue triangles Re a = 10) -ABC Runs 9,11,13 . 4.21 Same as Figure 4.20 but for the OT Runs 10,12,14. . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Scatter plots of the ratio ε o /ε t versus the ratio ε a /ε t for AD-OT runs 12 and 14. . . . . . . . 5.1 Time evolution of the ratio of kinetic to magnetic energy E k /E m . Time is in units of l 0 /u r.m.s. . Blue: OT initial condition M s = 0 (incompressible). Green: OT initial condition M s = 1. Red: OT initial condition M s = 4. Dashed lines: same colors for ABC initial condition. . . . 5.2 Time evolution of total physical and numerical heating rates. Blue: M s = 0 (incompressible).

Green: M s = 1. Red: M s = 4. Solid lines: total heating rate -Ėtot (see text). Dashed lines: physical heating rate ε phys . Diamonds: total heating rate ε corr as recovered by method 1 of the appendix (see text). Left: ABC initial condition. Right: OT initial condition. . . . . . 5.3 Time evolution of the ratio of viscous to total physical dissipation rates ε v /( ε v + ε o ).

Blue: M s = 0 (incompressible). Green: M s = 1. Red: M s = 4. Solid lines: OT initial conditions. Dashed lines: ABC initial conditions. . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Time evolution of the ratio of dilatational to total viscous dissipation rates (vortical plus dilatational) ε vc /( ε vc + ε vω ). Green: M s = 1. Red: M s = 4. Solid lines: OT initial conditions. Dashed lines: ABC initial conditions. . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Density cuts for runs 1-4. Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT).

Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT). . . . . . . . . . . . 5.6 Column density (i.e. density field integrated over the whole cube) maps for runs 1-4. Top left:

Run 1 (M s = 1, ABC). The diffuse interstellar medium

The interstellar medium (ISM), defined as the gas and dust that lies between the stars in a galaxy, is a complex multi-component system, and a place where a wide range of physical processes take place. In this introductory chapter we provide a very brief description of the ISM, focusing mainly on its diffuse component.

For a comprehensive treatment of general ISM physics, the reader is referred to the classic book of [START_REF] Spitzer | Physical Processes in the Interstellar Medium[END_REF], the recent book of [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF], as well as the review articles of [START_REF] Ferrière | The interstellar environment of our galaxy[END_REF] and [START_REF] Cox | The Three-Phase Interstellar Medium Revisited[END_REF].

Components

The ISM is a very important component of every galactic ecosystem. Stars, the other most important component, are formed by gravitational collapse of the densest regions of the ISM. Their stellar winds during their lifetime or their supernovae explosions and planetary nebulae at the end of their life-cycle are able to enrich the ISM with elements heavier than helium (metals), thus building up a complex chemical environment in which a rich variety of molecules is present. As a galaxy evolves, the ISM is converted to stars which in turn enrich it with metals. Part of the ISM may also leave the galaxy in the form of galactic winds while conversely in-falling gas from the intergalactic medium may add to the ISM. A broad definition of the ISM as everything in the galaxy that lies between the stars would include several constituents. Following the beautiful and inspiring text of [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF]:

• Interstellar gas: Gas-phase ions, atoms or molecules with a nearly thermal velocity distribution.

• Interstellar dust: Small solid particles with diameter in general less than ∼ 1 µm, or very large molecules.

• Cosmic rays: Electrons, protons or heavier nuclei with kinetic energies much larger than the thermal energy, often in the relativistic regime.

• Electromagnetic radiation: Photons from various sources such as the cosmic microwave background, stellar photo spheres, radiation emitted by interstellar ions, atoms or molecules, thermal emission by dust particles that have been heated by starlight, free-free emission from interstellar plasma, synchrotron radiation from relativistic electrons, gamma rays emitted in nuclear transitions and π 0 decays.

• Interstellar magnetic field: The magnetic field produced by electric currents in the ISM. This field guides the cosmic rays and in some parts of the ISM it is strong enough to be dynamically important.

• The gravitational field: The field due to the ISM gas, stars, stellar remnants and dark matter. In the densest regions of the ISM, the gas-phase contribution leads to self-gravitating clouds.

• Dark matter: The extent to which this component interacts with baryons in a non-gravitational way is currently unknown.

Although the ISM accounts for most of the volume of a galaxy, its estimated mass fraction does not exceed 15 % of the total mass of the galactic disk in the case of the Milky Way. If dark matter is included, the contribution of the ISM is as low as 0.5 %. This mass fraction is distributed 99 % in the gas phase and 1 % in dust. The density distribution of the ISM is highly inhomogeneous, with half of the mass contained in high-density regions occupying only 1-2 % of the total volume.

The chemical composition of the ISM in the Milky Way is estimated to be 90.8 % by number (70.4 % by mass) of hydrogen, 9.1 % (28.1 %) helium and 0.12 % (1.5 %) of metals [START_REF] Spitzer | Physical Processes in the Interstellar Medium[END_REF]. About 60 % of the interstellar hydrogen is in the form of neutral H atoms, 20 % is in the form of H 2 molecules and the remaining 20 % is ionized. A significant fraction of the heavier elements is locked inside dust grains, but the depletion factor may vary significantly between different elements as well as between different regions of the ISM. On average, the most abundant metals such as C, N, and O are depleted by a factor of ∼ 1.2 -3, while elements such as Mg, Si or Fe can have depletion factors as large as ∼ 10 -100.

Phases

The ISM is not a homogeneous medium, on the contrary it is organized in several phases with distinct physical characteristics. Again following [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF] and many others, the interstellar medium can be divided into the following phases:

• The hot ionized medium: This phase, also called the coronal gas, consists of gas that has been heated to temperatures of the order of 10 6 K by shock waves originating from supernova explosions. The primary ionization mechanism is collisional ionization, with ions such as OIV present. The hot ionized medium fills a significant fraction, almost 50 %, of the total ISM volume. Typical number densities in this phase are of the order of 4 × 10 -3 cm -3 .

• HII gas: In this phase, hydrogen is photoionized by ultraviolet radiation from hot O-type stars. The photoionized gas may be dense material from a nearby cloud (HII region) or lower density intercloud medium (diffuse HII or warm ionized medium). The temperature is of the order of 10 4 K while the number density can vary greatly from 0.3 to 10 4 cm -3 .

• Warm HI: Gas heated to temperatures of the order of 10 4 K, but remaining in atomic form. A characteristic scale for the number density of this phase is 0.6cm -3 . This phase fills a significant fraction of the total volume, perhaps as high as 40 %.

• Cool HI: Often referred to as the cold neutral medium (CNM), this phase consists of predominantly atomic gas at temperatures T ∼ 10 2 K and number densities n H ∼ 30cm -3 . The volume fraction of this phase is of the order of 1 %.

• Diffuse molecular gas: In this phase, the physical conditions are similar to those of cold HI clouds (T ∼ 50K, n H ∼ 10 2 cm -3 ), but the column densities are high enough for hydrogen to be in predominantly molecular form.

• Molecular clouds: These clouds are gravitationally bound and have temperatures in the range of 10 to 50 K and number densities in the range of 10 3 to 10 6 cm -3 .

Early models of the interstellar medium [START_REF] Mckee | A theory of the interstellar medium -Three components regulated by supernova explosions in an inhomogeneous substrate[END_REF] classified interstellar matter into three phases: the cold neutral medium (CNM), the warm phase which comprises of the warm ionized medium (WIM) and the warm neutral medium (WNM), and the hot ionized medium (HIM). In these models, the cold neutral medium is considered as a collection of clouds, while the hot ionized medium is considered as an intercloud medium. The warm phase is considered as a boundary layer between the clouds and the intercloud medium. These three phases were thought to be in approximate pressure equilibrium with each other.

The cold neutral medium comprises a variety of clouds characterized by a wide range of physical conditions. Observationally, the densest clouds which are most protected from ultraviolet radiation are referred to as dense clouds, dark clouds or molecular clouds. The most tenuous clouds which are fully exposed to starlight are usually called diffuse clouds. Clouds between these two extremes have been called translucent clouds. The above classification has not been fully consistent in the literature though.

Cloud type

Diffuse atomic Diffuse Molecular

Translucent Dense Molecular Definition f n H2 < 0.1 f n H2 > 0.1, f n C + > 0.5 f n C + < 0.5, f n CO < 0.9 f n CO > 0.9 A V (minimum) 0 ∼ 0.2 ∼ 1-2 ∼ 5-10 Typ. n H (cm

Classification of clouds

Snow [START_REF] Snow | Diffuse Atomic and Molecular Clouds[END_REF] have proposed a systematic classification of interstellar clouds based on predominant form of hydrogen and carbon. Let the number density of a particular element be denoted by n, for example n H for hydrogen. This number density corresponds to the total number of nuclei in a cubic centimeter, whatever their molecular form. For hydrogen, n H n(H) + 2n(H 2 ). For the case of carbon, most nuclei are found in the form of C + ,C or CO, so we can approximately write n C n(C + ) + n(C) + n(CO). The fraction of an atom's nuclei that is found in a particular form is referred to as f n , for example f n C for neutral carbon. The classification proposed by [START_REF] Snow | Diffuse Atomic and Molecular Clouds[END_REF] is shown in table 1.1.

Diffuse atomic clouds are the most tenuous regions of the neutral interstellar medium. The low column densities of these clouds make them essentially fully exposed to starlight. The predominant form of hydrogen is neutral atomic, since most of molecular hydrogen is photo-dissociated by starlight. Atoms with ionization potentials lower than hydrogen, of which carbon is a notable case, are almost fully ionized, providing abundant electrons. The paucity of molecules due to photo-dissociation implies that very little chemistry happens in these clouds. Diffuse atomic clouds typically have a fairly low total number density n ∼ 10 -100 cm -3 and temperatures in the range 30 -100K.

In diffuse molecular clouds, the interstellar radiation field is sufficiently attenuated, at least in the individual wavelengths that dissociate H 2 , so that the local fraction of hydrogen in molecular form f n H2 becomes substantial (> 0.1). However, the interstellar radiation field is still strong enough to photo-ionize any atomic carbon present and dissociate any CO, so that carbon remains predominantly in the form of C + (f n C + > 0.5). In steady state, diffuse molecular clouds must necessarily be surrounded by layers of diffuse atomic gas, which provide the shielding of radiation. Thus, most sightlines that cross a diffuse molecular cloud will also cross diffuse atomic gas. The presence of abundant H 2 in diffuse molecular clouds activates certain chemical routes that lead to formation of molecules such as CO, CH, CN, C 2 and C 3 which are visible in UV/visible absorption, CO and H + 3 which are visible in the infrared, and HCO + , OH and C 2 H which are visible at millimeter wavelengths. Diffuse molecular clouds typically have densities that range from 100 -500 cm -3 and temperatures that range from 30 -100 K.

In a translucent cloud, the increased attenuation of the interstellar radiation field allows carbon to begin the transition from ionized atomic (C + ) to neutral atomic (C) and molecular (CO) form. The chemistry is qualitatively different than in a diffuse molecular cloud, both because of the decreased electron fraction and the presence of the highly reactive C atoms. Translucent clouds are defined as clouds where f n C + < 0.5 and f n CO < 0.9. This definition reflects the fact that C + is no longer the predominant form of carbon, but also differentiates a translucent cloud from a dense molecular cloud, where the carbon is almost exclusively in CO form. In steady state, a translucent cloud is surrounded by diffuse molecular cloud material, which is in turn surrounded by diffuse atomic cloud material in an onion-like structure.

If the visual extinction is further increased, carbon becomes almost completely molecular (f n CO ∼ 1), defining the regime of dense molecular clouds. There is again a qualitative difference in chemistry, as the electron abundance is very low, with cosmic-ray ionization being the dominant source, and the reactive C is replaced by the very stable CO. Again, dense molecular cloud material is found surrounded by translucent material in an onion-like structure. These clouds are often self-gravitating, with densities n at least 10 4 cm -3 and kinetic temperatures in the range of 10 -50 K. Chemistry in dense molecular clouds is very rich, as most of the 130 currently known interstellar molecules were found through observations in such clouds.

Energy content

The following section is, again, closely inspired by [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF]. In the interstellar medium, energy is found in a number of forms: turbulent kinetic energy of the gas motions 1 2 ρu 2 , magnetic energy 1 8π B 2 , cosmic ray kinetic energy u CR , thermal energy 3 2 nkT , and energy of the interstellar radiation field, which can be further subdivided in the cosmic microwave background, far-infrared emission form dust and starlight. It is a remarkable fact that in the local interstellar medium, all these forms of energies have comparable energy densities in the range of 0.2 -2 eV cm -3 , as can be seen in table 1.2.

Although the equality of the energy density of the cosmic microwave background with the other energy densities could be coincidental, the other energy densities are coupled by dynamical processes which lead to equipartition, as the following examples show. First, the turbulent energy density of the gas motion 1/2ρu 2 is coupled to the magnetic energy density B 2 /8π through the equations of magnetohydrodynamics describing the turbulent motion of an ionized gas threaded by a dynamically important magnetic field. More specifically, the magnetic field feeds back on the velocity field through the Lorentz force in the momentum equation, whereas in turn the turbulent velocity field continuously stretches and distorts magnetic field lines (see section 2.2). Next, if the cosmic ray energy density were much higher than the magnetic energy density, the galactic magnetic field would not be able to confine the cosmic rays inside the galactic disk; they would freely escape into intergalactic space. This negative feedback limits cosmic ray energy density to a level comparable to the energy density of the galactic magnetic field.

Component

Energy density u (eV cm [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF]).

The interstellar medium is constantly in a state far from thermodynamic equilibrium, and what allows it to remain in such a state is the constant input of free energy in the form of starlight, kinetic energy from supernovae ejecta and infall from the intergalactic medium. Ultimately, all of this free energy escapes the galaxy in the form of photons propagating in the cold extragalactic sky.

Chemistry in the diffuse ISM

Despite the very low number densities in the diffuse ISM, since the 1930s and the advent of high-resolution optical spectroscopy there has been strong evidence of a rich molecular chemistry [START_REF] Eddington | Diffuse matter in interstellar space[END_REF][START_REF] Swings | Les constituants physiques de l'espace interstellaire[END_REF]. These then-newly discovered spectral lines where attributed to molecules such as CN and CH + . Today, a rich variety of molecules have been detected in interstellar space [START_REF] Tielens | The molecular universe[END_REF]. In this section, we provide a brief review of the paths leading to the formation of molecules in the diffuse ISM.

Molecules in the diffuse ISM

Gas in the diffuse ISM is in a very dilute and cold state (table 1.1). Collisions are thus very infrequent, with inelastic collisions being even less frequent. Prior to the discovery of interstellar molecular lines mentioned above, chemistry was thought to take little part in the physics of the diffuse ISM. However, even if collisions are very infrequent, they can still lead to the formation of molecules if the time scales involved are large enough.

One important factor in interstellar molecular chemistry is dust, since H 2 is known to form in the surface of dust grains, which act as reaction catalysts [START_REF] Draine | Physics of the Interstellar and Intergalactic Medium[END_REF]. The presence of H 2 opens the possibility of a variety of reactions to take place. In these reactions, a proton is removed from the H 2 molecule by an inelastic collision with an atom or ion X and it is replaced by it in the product: H 2 + X → HX + H.

(1.1)

These reactions however often have activation barriers of several thousand K, depending on the reactant X [START_REF] Lesaffre | Low-velocity shocks: signatures of turbulent dissipation in diffuse irradiated gas[END_REF]. This raises the question of how these molecules are formed in such a dilute and cold medium.

One possible formation route is the sequence of reactions

H 2 + H + 2 → H + 3 + H H + 3 + X → XH + + H 2 .
(1.2)

In this sequence, the formation of H + 3 requires a strong ionizing source such as cosmic rays or irradiation. Another possible route is to have a heating source which is able to raise the temperature, at least locally, to values above the activation barrier of the reaction 1.1. This heating source can be provided by hydrodynamic processes such as shocks and turbulent dissipation.

Energy dissipation and molecular chemistry

Different natures of dissipation

The irreversible processes which convert ordered energy (magnetic or kinetic) into heat can be of different nature. For instance, in the J-type shock above, viscous friction is the main actor. But magnetized shocks come in an other flavor: in C-type shocks, charges tied to the magnetic field and the neutrals species have different velocities and it is the ambipolar friction between the ions and neutrals which is the main heating process.

The heating rate is usually much smaller (cf. figure 1.1), so it takes a longer time and length scale to dissipate the energy flow which enters the shock: viscous friction makes much smaller structures than ambipolar heating for the same energy dissipated (both shocks on the figure have same velocity and magnetic field). As a result, the temperature inside a C-type shock is much smaller. But the velocity difference between the two fluids provides extra energy in the frame in which neutral-ion reactions take place: the effective temperature for neutral-ions reactions is raised by the velocity drift.

For instance, the effective temperature, which is defined as

T ef f [H 2 , C + ] = T r + T s (1.7)
where

T r = µ n T n + µ i T i µ n + µ i (1.8)
is the ion-neutral kinetic temperature and is the contribution of the relative macroscopic kinetic energy of the ionized and neutral particles, is plotted in figure 1.1. µ n is the atomic weight of the neutrals, µ i the atomic weight of the ions and µ in = µ i µ n /(µ i +µ n ).

It reaches a maximum of 4550 K and allows the reaction

C + + H 2 → CH + + H (1.10)
which has an endothermicity of 4640 K to proceed. This opens the subsequent chain of

CH + + H 2 → CH + 2 + H CH + 2 + H 2 → CH 3 + H 2 (1.

Models of turbulent dissipation

We have focused here on the example of shocks, but many other ways of dissipating energy exist, such as current sheets or vortices, for example. The turbulent dissipation region (TDR) model of [START_REF] Godard | Models of turbulent dissipation regions in the diffuse interstellar medium[END_REF] is an attempt to model the effects of turbulent energy dissipation on ISM chemistry using the magnetized Burgers vortex as a building block. This model yields a different chemical signature compared with the chemical signature of models based on J and C-type shocks.

The current thesis aims at understanding better the local structure of dissipation in magnetohydrodynamic turbulence, through the statistical study of the geometrical and dynamical characteristics of structures of high dissipation, as they are observed in numerical simulations.

Chapter 2

Turbulence

Turbulence1 is an ubiquitous state of fluid flow in which the flow variables, such as velocity or temperature, exhibit complex variation in both space and time. It is a phenomenon widely observed in both terrestrial environments such as the atmosphere and the oceans and extraterrestrial environments such as stellar interiors, the ISM and the intra-cluster medium in galaxy clusters (figure 2.1). An important difference between terrestrial and astrophysical turbulence is that while terrestrial turbulence is always a flow of a neutral fluid, such as water or air, astrophysical turbulence almost always concerns the motion of an ionized fluid threaded by a dynamically important magnetic field. Turbulence in the ISM is an important example of astrophysical turbulence, and has significant consequences on ISM physics and chemistry.

In this chapter we will attempt a self-contained presentation of the vast turbulence problem, with particular emphasis on issues relevant to energy dissipation and intermittency. After a brief introduction to the general problem, we will present the governing equations before proceeding to a summary of the statistical theory of turbulence followed by a presentation of the phenomenological theories of the inertial range. After a section on intermittency, we conclude this chapter with a section on interstellar turbulence.

Introduction

Turbulent flows lie at the interface between order and chaos: although individual realizations of an experiment flow variables can exhibit chaotic variation, there often exists predictable behavior in the statistical sense, if a sufficiently large number of physical experiments is considered. This important property has led to the development of a statistical theory, which is fundamentally based on the notion of an ensemble average: an experiment is repeated a large number of times in such a way that although the initial and boundary conditions are macroscopically identical, small perturbations in them lead to a sample of turbulent flow fields from which an average can be extracted. In practice, the ensemble average is often replaced with an average with respect to space or time, assuming a corresponding ergodic hypothesis.

This variability between different experiments with macroscopically identical initial conditions is known as sensitive dependence on initial conditions and is a signature of chaotic dynamics of nonlinear systems. This phenomenon has been studied extensively in the context of low-dimensional dynamical systems typically described by ordinary differential equations. Turbulence, being a phenomenon extended in both space and time, has in principle an infinite number of degrees of freedom and is the prototypical example of spatiotemporal chaos, where a very large number of degrees of freedom interact non-linearly.

But how much are these degrees of freedom excited, as a function of their characteristic scale? Contrary to other complex phenomena studied in statistical physics, in turbulence there is no clear-cut separation of scales. Energy as a function of scale is rather characterized by a continuous spectrum, which can be accounted for by simple phenomenological assumptions. Turbulence is thus a multiscale phenomenon, where different degrees of freedom, each having a specific scale, interact with each other to produce a continuous energy spectrum. [START_REF] Corke | A new view on origin, role and manipulation of large scales in turbulent boundary layers[END_REF]. Top right: Phytoplankton bloom reveals turbulent ocean currents, an example of geophysical turbulence [START_REF] Signorini | Environmental factors controlling the barents sea spring-summer phytoplankton blooms[END_REF]. Bottom left: Turbulence in the interstellar medium as revealed by the gradient of linear polarization [START_REF] Gaensler | Low-mach-number turbulence in interstellar gas revealed by radio polarization gradients[END_REF]. Bottom right: Current density magnitude in a high resolution numerical experiment of homogeneous isotropic magnetohydrodynamic turbulence [G. Momferratos, unpublished].

This interaction between structures at different scales leads to the phenomenological picture of the energy cascade. Energy is injected at large scales, either by a driving force or by appropriate boundary conditions. At scales much smaller than the energy containing scale, but still large enough so that viscosity is unimportant, nonlinear interaction between modes transfer energy from large to small scales in a conservative way; this range of scales is called the inertial range. Finally, energy reaches scales small enough such that viscosity cannot be neglected, and is dissipated into heat there; this range of scales is known as the dissipative range.

Energy dissipation is however not uniformly distributed in space and time, but is characterized by strong localized bursts alternating between larger quiescent regions. This property of fluid turbulence, known as intermittency, has important theoretical as well as practical consequences. On the theoretical side, it is capable of modifying the phenomenological picture based on the energy cascade in a nontrivial way. On the practical side, the intermittency of the energy dissipation rate is capable of heating the interstellar medium above the threshold of certain important chemical reactions, thus activating new chemical routes.

Governing equations 2.2.1 Incompressible models

Hydrodynamics

Incompressible turbulence in a neutral fluid is governed by the Navier-Stokes equation [START_REF] Landau | Fluid Mechanics[END_REF] ∂

t u + (u • ∇) u = -∇p + ν∇ 2 u (2.1)
together with the incompressibility condition

∇ • u = 0. (2.2)
In the above, u is the velocity field and ν = µ/ρ is the kinematic viscosity, with µ the dynamic viscosity and ρ the constant mass density. p = p/ρ is the actual pressure p divided by the density.

If the velocity is rescaled by a reference velocity U , spatial coordinates are rescaled by a reference length L and the time coordinate is rescaled by L/U , the Navier-Stokes equation remains unchanged except from the substitution of the viscosity ν with the inverse of the Reynolds number

ν -1 → Re = U L ν . (2.3)
Two flows in domains which are geometrically similar are identical if their Reynolds number are equal. From a physical point of view, the Reynolds number measures the relative importance of the nonlinear advection term (u • ∇) u with respect to the viscous term ν∇ 2 u. A flow with low Reynolds number is a smooth laminar flow dominated by viscosity. As the Reynolds number is increased, successive instabilities take place and the laminar flow gradually transitions to turbulence. A flow with very high Reynolds number, presumably an asymptotic state, is known as fully developed turbulence.

The quadratic nonlinear term, the second term on the left hand side of equation (2.1), can be identified as the source of the complexity of turbulence. As we will see below, this term is responsible for the interactions that take place between different scales of the flow.

In incompressible flows, the pressure p is not a true dynamical variable, but simply an agent enforcing the incompressibility condition. This can be seen by taking the divergence of equation (2.1), keeping in mind that the velocity field is solenoidal. One gets the following Poisson equation for the pressure

∇ 2 p = -∂ i u j ∂ j u i = -∂ i ∂ j (u i u j ) (2.4)
where the equivalence of the two forms of the source term on the right hand side follows from the incompressibility condition (2.2). The Poisson equation (2.4) is an elliptic equation, and this has the important consequence that a local change in the velocity field is propagated instantaneously across the whole flow field through the pressure term. In physical terms, this corresponds to an infinite sound velocity.

Using standard vector identities, one can rewrite the Navier-Stokes equation in rotational form

∂ t u + ω × u = -∇ p + 1 2 u 2 + ν∇ 2 u (2.5)
where ω = ∇ × u is the vorticity, a measure of the local rotation rate of the fluid. Taking the curl of this equation, one obtains an evolution equation for the vorticity ω

∂ t ω = ∇ × (u × ω) + ν∇ 2 ω (2.6)
which is equivalent to

∂ t ω + (u • ∇)ω = Dω Dt = (ω • ∇)u + ν∇ 2 ω.
(2.7)

Thus, as the vorticity vector is transported, it is stretched by velocity gradients and diffused by viscosity. Indeed, let us consider a line L moving with the fluid and let δl be an element on this line. If v is the velocity of the fluid at the one end of the line, say x, the velocity at

x + δl is v + (δl • ∇)u.
Thus in a time interval dt, the line element δl changes by dt(δl • ∇)u. Therefore

Dδl dt = (δl • ∇)u (2.8)
and, if we neglect viscous effects, the evolution of the vectors δl and ω are given by identical formulae. It follows that if these two vectors are initially in the same direction, they will remain parallel, and the ratio of their moduli will remain constant. If we define the vortex lines as the integral curves of the vorticity field, ie. the solutions of the system of ordinary differential equations

dr dt = ω(r, t) (2.9)
two fluid particles will always be on the same vortex line, and the value of |ω| will be proportional to the distance between the particles. This vortex-stretching process is the fundamental mechanism for the transfer of energy from large to small scales in three-dimensional incompressible hydrodynamic turbulence. It is also the primary source of difference between three-dimensional and two-dimensional turbulence, as the term (ω • ∇)u is identically zero in the latter case, the vector ω being by definition perpendicular to the plane where the vector u lies.

It is possible to investigate the role of the nonlinear term on scale interactions by considering the Navier-Stokes equation (2.1) in Fourier space:

u(x) = u(k) exp(ik • x) dk.
(2.10) By this transformation, the velocity field is decomposed into plane waves of wavevector k, which are completely localized in wavenumber space. It is the most simple way to separate the different contributions to the velocity field coming from different scales and study their interactions. The Fourier coefficients are given by

u(k) = 1 (2π) 3 u(x) exp(-ik • x) dx.
(2.11)

By applying the Fourier transform (2.11) to the Navier-Stokes equation (2.1) and using equation (2.4) to eliminate the pressure term, one can derive the following spectral form

∂ t u i (k) + νk 2 u i (k) = k i P ij (k) k=p+q u i (q)u j (p) dpdq (2.12)
where the projection operator

P ij (k) = δ ij - k i k j k 2 (2.13)
appears when removing the pressure by use of the Poisson equation (2.4). The presence of this projection operator ensures the conservation of incompressibility. The scale interaction mediated by the nonlinear term is now evident in the right hand side of equation (2.12). The nonlinear term, which is local in physical space, is non-local in Fourier space, effectively a convolution integral. Two perturbations with wavevectors p and q interact to excite a perturbation with wavevector k, but only if these three wavevectors form a triangle: k = p + q. Although this form of the Navier-Stokes equation clearly shows the existence of scale interactions, it provides no information on their locality or non-locality in scale. It also provides no information on the direction of such interactions, ie. if there is, on average, a flow of energy from large to small scales or the converse.

In the ideal limit where viscosity tends to zero, the dynamics of the fluid is governed by the Euler equation

∂ t u + (u • ∇) u = -∇p, ∇ • u = 0. (2.14)
The incompressible Euler equation (2.14) is compatible with the following symmetries [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]:

• Space translations: g space ρ : t, r, u → t, r + ρ, u, ρ ∈ R 3 .
• Time translations: g time τ : t, r, u → t + τ, r, u, τ ∈ R.

• Galilean transformations:

g Gal U : t, r, u → t, r + Ut, u + U, U ∈ R 3 . • Parity P : t, r, u → t, -r, -u. • Rotations g rot A : t, r, u → t, A • r, A • u, A ∈ SO(R 3 ). • Scaling g scal λ : t, r, u → λ 1-h t, λr, λ h u, λ ∈ R, h ∈ R + .
where SO(R 3 ) is the groups of rotations in R 3 . When the viscosity is equal to zero, the Euler equations possess an infinite number of scaling groups labeled by a scaling exponent h. When the viscosity is finite, the Navier-Stokes equation (2.1) allows only h = -1. This value corresponds to the similarity principle of fluid dynamics: in this case, the scaling transformations g scal λ leave the Reynolds number unchanged. The Euler equation (2.14) possesses two invariant quantities: the first is kinetic helicity

H k = V u • ω dV (2.15)
and the second is kinetic energy

E k = V 1 2 u 2 dV (2.16)
which is quadratic and positive-definite. Quadratic invariants are particularly important because they remain invariants if the fields are developed in Fourier series and then truncated at a finite wavenumber.

In the Navier-Stokes case where viscosity is finite, kinetic energy is dissipated locally at the rate

ε = ν 2 3 i,j=1 (∂ i u j + ∂ j u i ) 2 .
(2.17)

The mean value of the energy dissipation rate is equal to

-Ėk = 1 V V ε(r) d 3 r = νΩ (2.18)
where Ω is the enstrophy

Ω = V ω 2 dV. (2.19)
The positive definite quadratic invariant E k , concentrated at the energy containing scale L i , where it is ejected by either a driving force or boundary conditions, is dissipated at the rate ε . By examining the form of the viscous term in Fourier space

ν∇ 2 u(x) → -ν|k| 2 u(k) (2.20)
it is easy to see that for ν small enough, dissipation occurs at small scales l d L i . Thus for ν small enough, there exists a range of scales l d l L i where both viscosity and boundary condition or injection effects are absent, and energy is conserved by nonlinear interactions; this range of scale is called the inertial range. The range of scales where viscosity dissipates kinetic energy is called the dissipative range.

The kinetic helicity H k is a measure of the topological complexity or the knottedness of vortex lines [START_REF] Moffatt | Helicity in laminar and turbulent flow[END_REF]]. In the ideal limit (zero viscosity), it is exactly conserved. It is not a positive definite quantity and, for a non-zero viscosity it is dissipated at the rate

-Ḣk = ν V ω • (∇ × ω) dV.
(2.21)

There is also an important local conservation law in ideal incompressible hydrodynamics, the law of conservation of circulation [START_REF] Landau | Fluid Mechanics[END_REF]]

dΓ dt = 0 (2.22)
where Γ is the velocity circulation

Γ = C u • dl (2.23)
considered around a closed contour C moving with the fluid, whose element is dl. The law of conservation of circulation can be recast in terms of the vorticity

Γ = A ω • dS (2.24)
where A is a surface enclosed by the contour C and dS is an element of this surface. This form of the law of conservation of circulation is closely related to vortex stretching: if the flow is such that the area of the surface A is greatly reduced, the modulus of the vorticity |ω| can reach very high values.

Magnetohydrodynamics

The interstellar medium is an ionized gas threaded by the galactic magnetic field, whose energy density is comparable with the energy density corresponding to the motion of the gas itself. Thus, the magnetic field is dynamically important and a model describing its interaction with the gas motion is called for. Incompressible magnetohydrodynamics describes the motion of a completely ionized incompressible fluid threaded by a magnetic field. The equations of incompressible magnetohydrodynamics can be derived from the Maxwell and Navier-Stokes equations under the assumptions of incompressibility and non-relativistic fluid velocities [START_REF] Landau | Electrodynamics of continuous media[END_REF].

∂ t u + (u • ∇) u = -∇p + j × b + ν∇ 2 u ∂ t b = ∇ × (u × b) + η∇ 2 b (2.25)
together with the two solenoidal conditions

∇ • u = 0 ∇ • b = 0 (2.26)
where b = B/ √ 4πρ is the Alfvén velocity with B the magnetic field, j = ∇ × b is the current density, j × b is the Lorentz force and η is the resistivity. Although the induction equation (the second of equations (2.25)) is linear in b the system is strongly non-linear because of the coupling of the two equations through the Lorentz force j × b on the right hand side of the momentum equation. In analogy with the term (u • ∇)u in the momentum equation, the advection term ∇ × (u × b) is responsible for the stretching of magnetic field lines which leads to interaction between structures at different scales.

If the equations of magnetohydrodynamics are written in non-dimensional form, then in addition to the kinetic Reynolds number U L/ν there also appears the magnetic Reynolds number

Re M = U L η (2.27)
replacing the inverse of the resistivity η in the induction equation. Here, as previously, U is a characteristic velocity scale and L a characteristic length scale. The magnetic Reynolds number measures the relative importance of advection versus diffusion in the dynamics of the magnetic field. In a flow with low magnetic Reynolds number, magnetic field dynamics is dominated by Ohmic diffusion whereas in the case of a high magnetic Reynolds number, it is dominated by advection and stretching.

The presence of both a viscosity and a resistivity in the equations of magnetohydordynamics, two quantities which have the same dimensions, introduces one more non-dimensional number, the magnetic Prandtl number P r = ν η .

(2.28)

This number measures the relative importance of viscosity with respect to resistivity. In a fluid with large magnetic Prandtl number, viscous diffusion of the velocity field is much stronger than Ohmic diffusion of the magnetic field. Just as in the case of the vorticity equation of incompressible hydrodynamics (2.6), the induction equation can be developed in the form

∂ t b + (u • ∇)b = Db Dt = (b • ∇)u + η∇ 2 b (2.29)
which shows that the magnetic field is stretched by velocity gradients as it is transported. An essential difference in the dynamics of the vectors b and ω is that the magnetic field is an independent vector, whereas the vorticity is a function of the velocity field: ω = ∇ × u.

Using standard vector identities, the Lorentz force on the right hand side of the momentum equation can be developed into

j × b = (b • ∇)b -∇ b 2 2 (2.30)
where the first term corresponds to the contribution of magnetic field line tension and the second term corresponds to an isotropic magnetic pressure.

In the ideal limit where both viscosity and resistivity are equal to zero, the equations of ideal magnetohydrodynamics

∂ t u + (u • ∇) u = -∇p + j × b ∂ t b = ∇ × (u × b) (2.31)
possess three invariant quantities. The first is total energy, which is the sum of kinetic and magnetic energy

E t = E k + E m , with E k = V 1 2 u 2 dV and E m = V 1 2 b 2 dV.
(2.32)

Total energy is a positive definite quantity, and is dissipated at the rate

-Ėt = ν ω 2 + η j 2 (2.33)
The first part is viscous dissipation, while the second part is Ohmic dissipation, essentially Joule heating. An important difference between viscous and Ohmic dissipation is that in the case of Ohmic dissipation, both global and local values are given by the same formula η j 2 , whereas in the case of viscous dissipation ε = ν ω 2 but ε(x) = ν(ω(x)) 2 locally. This is because the velocity field is a true vector, whereas the magnetic field is a pseudo-vector defined as the curl of the vector potential [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF].

The second conserved quantity is magnetic helicity

H m = V a • b dV (2.34)
where a is the vector potential defined by b = ∇ × a.

(2.35)

In analogy with kinetic helicity, magnetic helicity is a measure of the topological complexity ("knottedness") of magnetic field lines [START_REF] Moffatt | Field Generation in Electrically Conducting Fluids[END_REF]. It is not a positive definite quantity, and it is dissipated at the rate

-Ḣm = η V j • b dV. (2.36)
The last conserved quantity is the cross-helicity

H c = V u • b dV.
(2.37)

The cross-helicity is a measure of the correlation (alignment in the mean sense) between the velocity and the magnetic field. As will be shown below, if this correlation is perfect the non-linearity of the equations is completely switched off. It is not a positive definite quantity and it is dissipated at the rate

-Ḣc = ν + η 2 V j • ω dV. (2.38)
In analogy with the conservation of velocity circulation in ideal incompressible hydrodynamics, ideal incompressible magnetohydrodynamics possess the local conservation law of magnetic flux

DΦ Dt = 0 (2.39)
where the magnetic flux Φ is defined as

Φ = A B • dS (2.40)
where A is a surface moving with the fluid and dS is an element of this surface. The equations of magnetohydrodynamics (2.25) can be written in a more symmetric form by introducing the Elsässer variables z ± = u ± b.

(2.41)

It is then a matter of simple algebra to transform equations (2.25) into

∂ t z ± + z ∓ • ∇ z ± = -∇ P + ν + ∇ 2 z ± + ν -∇ 2 z ∓ (2.42)
where P = p + 1/2b 2 is the total pressure and ν ± = (ν ± η)/2. In the presence of a mean magnetic field the nonlinear term (z ∓ • ∇) z ± represents counter-propagating Alfvén wave-packets, which play an important role in the physics of magnetohydrodynamic turbulence. These are transverse waves propagating in both directions along magnetic field lines with the velocity B 0 / √ 4πρ, where B 0 can be either a global mean field or a local average field.

The Elsässer form (2.42) of the equations of incompressible magnetohydrodynamics reveals two important properties. The first is that the nonlinearity is completely switched off in the case where the velocity is aligned (or anti-aligned) with the magnetic field. Indeed, if u = b then z -= 0, whereas if u = -b then z + = 0 and the equations (2.42) become linear. The second is that if a mean magnetic field B 0 is present, it appears with a plus sign in the variable z + and with a minus sign in the variable z -, and thus cannot be removed by a Galilean transformation. This is in contrast with incompressible hydrodynamics, where the form of the nonlinear term allows a mean velocity field to be removed.

Ambipolar diffusion magnetohydrodynamics

Many diffuse clouds are characterized by extremely low ionization fractions. Although the magnetic field remains dynamically important in these clouds, the classical model of incompressible magnetohydrodynamics is not sufficient to describe the evolution of their gaseous component. The most complete description of the incompressible motion of a partially ionized gas is the two-fluid model, which comprises of one momentum and continuity equation for the neutral component, another pair of equations for the ion component and the induction equation for the magnetic field

∂ t ρ n + ∇ • (ρ n u n ) = 0 ρ n (∂ t u n + (u n • ∇) u n ) = -∇p n + µ n ∇ 2 u n + γρ i ρ n (u n -u i ) ∂ t ρ i + ∇ • (ρ i u i ) = 0 ρ i (∂ t u i + (u i • ∇) u i ) = -∇p i + ρ i (j × b) + µ i ∇ 2 u i + γρ i ρ n (u i -u n ) ∂ t b = ∇ × (u i × b) + η∇ 2 b.
(2.43)

Here, u n is the neutral velocity, p n is the neutral pressure, µ n is the neutral dynamic viscosity, ρ n is the neutral density and γ is the coefficient of ion-neutral friction. The corresponding quantities of the ionized component appear with the index i , and b = B/ √ 4πρ i is the Alfvén velocity, with B the magnetic field. Because the neutrals do not interact directly with the magnetic field, the Lorentz force only appears in the ion momentum equation and in the induction equation, the magnetic field is transported by the ion velocity.

For low ionization fractions ρ i /ρ n , one can neglect ion inertia and viscosity, and the ion momentum equation becomes

ρ i (j × b) -∇p i = γρ i ρ n (u n -u i ).
(2.44) Substituting this expression in the neutral momentum equation and in the induction equation one gets

∂ t u + (u • ∇) u = -∇p + j × b + µ∇ 2 u ∂ t b = ∇ × (u × b) + α∇ × ((j × b) × b) + η∇ 2 b, α = (γρ) -1 .
(2.45)

The above approximation is known as the strong coupling approximation or the single fluid approximation.

The new velocity u is the center-of-mass velocity, which is approximately equal to the neutral velocity for low ionization fractions

u = ρ n u n + ρ i u i ρ n + ρ i u n (2.46)
while the density ρ is equal to the sum of the two components

ρ = ρ n + ρ i ρ n (2.47)
The nonlinear diffusion term α∇ × ((j × b) × b) that appears in the induction equation is the ambipolar diffusion term. It is third order in the magnetic field and second order in spatial derivatives. Physically, it expresses the diffusion of magnetic field lines due to the ion-neutral drift.

If the equations of ambipolar diffusion magnetohydrodynamics are written in non-dimensional form, then in addition to the kinetic and magnetic Reynolds numbers, there also appears, in front of the ambipolar diffusion term, the ambipolar diffusion Reynolds number

Re a = Lγρ i U = L l a = T t a (2.48)
where U is a characteristic velocity scale, L is a characteristic length scale, T = L/U while l a and t a are the characteristic length and time scales of ambipolar diffusion. The ambipolar diffusion Reynolds number measures the relative importance of ambipolar diffusion with respect to advection. In a flow with a low ambipolar diffusion Reynolds number the magnetic field dynamics is dominated by the ion-neutral drift, whereas the limit Re a → ∞ corresponds to ordinary magnetohydrodynamics. Note that [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF] constructed the ambipolar diffusion Reynolds number of eddies of length scale and velocity U based on the AD diffusion coefficient: R AD ( ) = U/λ AD = / AD where AD = U t a /M 2 a with M a = U/c a . [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF] then argue that only eddies of length scales below AD should be affected by AD. We prefer to get a similar estimate by comparing the Fourier amplitudes of the AD e.m.f. (Re -1 a (j × b) × b → Re -1 a kb 3 ) and the inertial e.m.f. (u × b → ub) in the induction equation (2.25). Wave numbers above the critical wave number

k a = Re a u 2 / b 2
(2.49) should be AD dominated. We hence define a = 2π/k a accordingly, as the length scale below which AD should be effective. Similarly, we can estimate the length scale below which the strong coupling approximation breaks down by comparing the magnitude of the neglected inertial term D t ρ i u i to the coupling term ρ n (u nu i )/t a . Assuming that u n , u i and u nu i all share the typical magnitude u, we then get a critical wavenumber k two-fluids ρ n /ρ i /l a above which the strong coupling approximation fails and the two-fluids approximation is needed. Provided ρ i /ρ n is small (it is typically lower than 10 -3 in the interstellar medium, for example if the main charges are C + ions), the strong coupling approximation breaks down at scales much smaller than the typical AD diffusion scale. Other authors [Oishi andMac Low, 2006b, Padoan et al., 2000] have claimed that the strong coupling approximation breaks down as soon as < AD or R AD < 1, where the ions inertia does not appear explicitly. But Fig. 1 of Li et al. [2006] shows a C-shock computed with the two-fluid approximation (solid) compared to an analytical solution (dashed) using the strong coupling approximation, and the agreement is perfect. We hence believe that the strong coupling approximation is a very good one in the low ionized ISM where ρ i /ρ n 1, even in cases where R AD > 1. The introduction of the ambipolar diffusion term in the equations of ideal incompressible magnetohydrodynamics breaks the conservation of total energy and cross-helicity, but not the conservation of magnetic helicity. This is because this term is written in the form

AD = ∇ × (u AD × b) .
(2.50)

The above form also guarantees the conservation of magnetic flux: although ambipolar diffusion is a dissipative process, it is unable to reconnect field lines.

Compressible models

The wide range of sonic Mach numbers observed in the interstellar medium implies that there exist areas where the incompressible models are not applicable. In a weakly compressible flow, density perturbations grow as the square of the sonic Mach number: ∆ρ/ρ ∝ M2 s 2 . Flows with sonic Mach numbers up to about 0.3 can perhaps be approximated as incompressible, but flows with M s ∼ 0.5 and above certainly cannot. Given that in the interstellar medium, the range of observed sonic Mach numbers is roughly 0.1 to 10 [ [START_REF] Elmegreen | Interstellar turbulence i: Observations and processes[END_REF], there certainly exist regions for which the incompressible models are completely inappropriate. In this subsection, we will present the compressible models which are capable of describing high sonic Mach number flows where shock waves are present.

Hydrodynamics

In compressible hydrodynamics, the law of conservation of mass is expressed by the continuity equation [START_REF] Landau | Fluid Mechanics[END_REF] ∂ t ρ + ∇ • (ρu) = 0 (2.51) which, by expanding the divergence, can be written in the equivalent form

∂ t ρ + (u • ∇)ρ + ρ∇ • u = 0 (2.52)
or, in terms of the material derivative of the density

D ln ρ Dt = -∇ • u. (2.53)
Written in terms of the specific volume

v = 1 ρ (2.54)
the above equations reveals the physical significance of the divergence of the velocity:

1 v Dv Dt = ∇ • u (2.55)
Positive divergence corresponds to an increase of the specific volume (dilatation) while negative divergence corresponds to a decrease (compression).

The conservation of momentum is expressed by the Navier-Stokes equation

ρ(∂ t u + (u • ∇)u) = -∇p + µ∇ 2 u + ζ + 1 3 µ ∇∇ • u (2.56)
where µ is the coefficient of shear dynamic viscosity and ζ is the coefficient of bulk viscosity, which models the viscous processes that take place because of changes in the specific volume of the gas. The bulk viscosity coefficient expresses the inability of the gas to adjust its density to rapid pressure changes as a result of a transfer of kinetic energy to internal degrees of freedom, such as rotational of vibrational molecular degrees of freedom [START_REF] Zeldovich | Physics of shock waves and high-temperature hydrodynamic phenomena[END_REF]. Thus, ζ is identically zero for a monoatomic gas which has no internal degrees of freedom. It is also zero for an isothermal gas, since in this case the density is a one-to-one function of the pressure. Below, we will take ζ equal to zero unless otherwise noted. Our starting point will thus be the following momentum equation

ρ(∂ t u + (u • ∇)u) = -∇p + µ∇ 2 u + µ 3 ∇∇ • u (2.57)
which is valid for a monoatomic gas with an arbitrary equation of state, or a (polyatomic or monoatomic) isothermal gas.

The system of equations of continuity (2.51) and momentum (2.57), which comprises four scalar equations, contains five unknowns, namely the three components of the velocity u, the pressure p and the density ρ. Thus, one more relation is needed to close the system. Such a relation is provided by the ideal gas law

p = (γ -1)ρe (2.58)
where γ = c p /c v is the ratio of specific heats and e is the internal energy per unit mass. c p = dQ/dT | p is the specific heat at constant pressure and c v = dQ/dT | v is the specific heat at constant volume. The ideal gas law calls for an additional partial differential equation for the internal energy

∂ t 1 2 ρu 2 + ρe = -∇ • ρu 1 2 ρu 2 + h -u • σ -κ∇T (2.

59)

where h = e + p/ρ is the enthalpy,

σ ij = µ ∂ i u j + ∂ j u i - 2 3 ∂ k u k δ ij + ζ∂ k u k δ ij (2.60)
is the viscous stress tensor, T is the temperature and κ is the coefficient of thermal diffusion. For a calorically perfect gas the temperature is related to the energy by

e = c v T (2.61)
and to the enthalpy by h = c p T.

(2.62)

In compressible hydrodynamics, density, velocity and pressure perturbations travel at the finite isentropic speed of sound c = γRT (2.63)

where R = c pc v is the universal gas constant.

The local heating rate due to viscosity is given by

ε = σ ij ∂ i u j .
(2.64)

If isothermality (T = const) is assumed, pressure is proportional to density:

p = c 2 i ρ (2.65)
where

c i = √ RT (2.66)
is the constant isothermal sound speed, at which perturbations travel. In this case, the energy equation is redundant, since the pressure is a function of the density alone.

If viscous processes are neglected, the dynamics of the ideal gas is governed by the Euler equations

∂ t u + (u • ∇)u = - ∇p ρ .
(2.67)

The continuity equation (2.51) remains unchanged, and energy conservation is best expressed though the entropy equation

∂ t S + (u • ∇)S = 0 (2.68)
where the entropy of an ideal gas is defined as

S = c v ln T -R ln ρ. (2.69)
Taking the curl of the Euler equation, we have

∂ t ω = ∇ × (u × ω) + 1 ρ 2 ∇ρ × ∇p (2.70)
The last term on the right hand side is the baroclinic term which represents generation of vorticity by density gradients. For a barotropic gas, of which an isothermal gas is a special case, p = f (ρ) and ∇ρ is parallel to ∇p and the baroclinic term is zero. Thus for an isothermal gas the vorticity transport equation is identical to the one describing an incompressible fluid.

For the compressible Euler equation (2.67), the kinetic energy

E = V 1 2 ρu 2 dV (2.71)
is an invariant, but contrary to incompressible ideal hydrodynamics, the kinetic helicity is not conserved in the general case, except in the case of a barotropic gas. The same applies for the law of conservation of circulation.3 An isothermal gas possesses the same invariants as an ideal incompressible fluid: kinetic energy and kinetic helicity, but in this case the kinetic energy is not a quadratic invariant, because of the presence of the density.

For the compressible Navier-Stokes equation the conserved quantity is the total energy e + 1 2 u 2 , which involves the internal energy of the gas.

Magnetohydrodynamics

Magnetic fields remain dynamically important in the regions of the interstellar medium where the gas flow is highly compressible, and the equations of compressible hydrodynamics are not sufficient to describe the physics of such regions. The equations of compressible magnetohydrodynamics describe the motion of a fully ionized magnetized compressible gas, again under the assumption of non-relativistic fluid velocities. The continuity equation is not modified by the presence of the magnetic field

∂ t ρ + ∇ • (ρu) = 0.
(2.72)

In the momentum equation, magnetic field feedback on the velocity field is represented by the addition of the Lorentz force J × B on the right hand side of the momentum equation

ρ(∂ t u + (u • ∇)u) = -∇p + J × B + µ∇ 2 u + ζ + 1 3 µ ∇∇ • u. (2.73)
The induction equation is not modified with respect to incompressible magnetohydrodynamics

∂ t B = ∇ × (u × B) + η∇ 2 B (2.74)
whereas the energy equation takes the form

∂ t 1 2 ρu 2 + B 2 8π + ρe = -∇ • u 1 2 ρu 2 + h + u • σ -κ∇T + 1 4π B × (u × B) - η 4π B × (∇ × B) .
(2.75) Again, this system of equations must be supplemented by an equation of state, which can be the ideal gas law or the assumption of isothermality.

If the effects of viscosity, resistivity and heat conduction are neglected, the system of ideal compressible magnetohydrodynamics reads

∂ t ρ + ∇ • (ρu) = 0 ρ(∂ t u + (u • ∇)u) = -∇p + J × B ∂ t B = ∇ × (u × B) (2.76)
In the absence of dissipative processes and assuming isothermality, the total energy

E = V 1 2 ρu 2 + B 2 8π (2.77)
is conserved. Magnetic helicity is also conserved because the induction equation (2.74) is identical to that of incompressible magnetohydrodynamics and it is easy to check that the cross-helicity is also conserved. The dissipation rates for the two helicities are identical to the ones corresponding to incompressible magnetohydrodynamics, with J replacing j, while the dissipation rate for the total energy is given by the sum of viscous dissipation (2.64) and Ohmic dissipation η J 2 .

Statistical theory

As experimental flow visualizations indicate, turbulence is a phenomenon that lies at the interface between order and chaos. The deterministic approach to the turbulence problem has certainly a lot to offer at the level of the mathematical study of solutions to equations of hydrodynamics and magnetohydrodynamics [START_REF] Gallavotti | Some rigorous results about 3d navier-stokes[END_REF][START_REF] Foias | Systems of conservation equations with a convex extension[END_REF], or the modeling of transition to turbulence by low-dimensional dynamical systems [START_REF] Manneville | Instabilités, chaos et turbulence[END_REF]. The statistical approach, starting with the pioneering work of [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF], [START_REF] Keller | Differentialgleichung für die turbulente bewegung einer kompressiblen flüssigkeit[END_REF] and [START_REF] Taylor | Statistical theory of turbulence[END_REF] has led to important predictions on the properties of fully developed turbulence. In this section, we will attempt a concise presentation of the statistical theory of turbulence. For a more complete exposition, the reader is referred to the classic books of [START_REF] Batchelor | The theory of homogeneous turbulence[END_REF], Monin and Yaglom [2007a] and Monin and Yaglom [2007b], as well as the modern books by [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF] and [START_REF] Lesieur | Turbulence in fluids[END_REF].

Statistical formulation of the turbulence problem

In a turbulent flow, the flow variables such as the velocity field4 exhibit erratic, complex behavior in both space and time. As an example, figure 2.2 shows the streamwise velocity component as recorded by a hot-wire anemometer in a classical wind tunnel experiment of homogeneous isotropic turbulence behind a grid. Two realizations of a turbulent flow under macroscopically identical initial conditions will produce two velocity fields that are completely different in their details. The turbulent velocity field is modeled as a stochastic vector field

u = u(x, t, ) (2.78)
where the stochastic argument represents the identity of a particular experiment. At least formally, the turbulence problem corresponds to the following: Given the statistical properties of the stochastic vector field u(x, t, ) at t = 0, predict its statistical properties for t > 0. [Batchelor, 1953, Monin andYaglom, 2007a]. A rigorous formulation of this problem is possible [START_REF] Vishik | Mathematical problems of statistical hydromechanics[END_REF], but rather intricate and to remote from the physics of turbulence. In the context of this thesis, we will not work with the general stochastic vector field u(x, t, ) but with its moments. [START_REF] Sreenivasan | Fractals and multifractals in fluid turbulence[END_REF].

Ensemble averaging

Macroscopically identical realizations are found in practice to possess many reproducible properties. These properties are revealed if one considers an ensemble of realizations of the turbulent field under the same initial and boundary conditions. Then, instead of focusing on the details of individual realizations one considers properties averaged over the ensemble. This process is known as ensemble averaging [Monin and Yaglom, 2007a], and it is symbolized by angle brackets within this section. For example, the ensemble-averaged energy dissipation rate ε is symbolized by ε . Mathematically, the ensemble average corresponds to the expectation value of the stochastic variable with respect to the stochastic argument

X = E(X( )).
(2.79)

Moments

The expectation value of a random variable, say the value of the x-component of the velocity at a particular point x 1 and time t 1 , can be expressed as an integral with the aid of the probability density function

E(u x (x 1 , t 1 )) = up (x1,t1) 1 (u) du. (2.80) Intuitively, p (x1,t1) 1
(u) du is the probability of the random variable u 1 (x 1 , t 1 ) to have a value in an interval of length du around u. One can define similarly the n-dimensional joint probability density function

p (x1,t1)(x2,t2)...(xn,tn) i1i2...in (y 1 , y 2 , ..., y n ) (2.81)
which expresses the joint probability density of the values of the random variables

u i1 (x 1 , t 1 ), u i2 (x 2 , t 2 ), ..., u in (x n , t n ).
(2.82)

In principle, one has complete information on the statistics of a turbulent flow only if all possible probability distributions for any value of n and any configuration of space-time points (x n , t n ) and indices i n are known. In practice however, such detailed knowledge is not possible, and one has to concentrate on lowerorder information. Such information is provided by the moments of the stochastic vector field u(x, t, ), which are defined as

M i1i2...in p1p2...pn (x 1 , t 1 , x 2 , t 2 , ..., x n , t n ) = u p1 i1 (x 1 , t 1 )u p2 i2 (x 2 , t 2 )...u pn in (x n , t n ) (2.83)
where the sum m = p 1 + p 2 + ... + p n is called the order of the moment.

Statistical symmetries

In practice, it is not always possible to build an ensemble of identical experiments, physical or numerical.

If the turbulent flow possesses a statistical symmetry, one can replace ensemble averaging by time or space averaging invoking a corresponding ergodic hypothesis. In addition, the presence of the statistical symmetry is able to simplify moment relations considerably.

A turbulent flow is said to be spatially homogeneous, or simply homogeneous if its statistical properties do not depend on the position vector x. Let us use the correlation tensor of the velocity field as an example,

B ij (x, t, x , t ) = u i (x 1 , t 1 )u j (x 2 , t 2 ) (2.84)
keeping in mind that what holds for the correlation tensor holds for any moment. It is obvious that if the statistical properties of the turbulence do not depend on position, then the correlation tensor cannot depend on x and x separately, but only on their difference r = xx:

B ij (x, t, x , t ) = B ij (r, t, t ), r = x -x.
(2.85)

Similarly, if the statistical properties of the turbulence do not depend on time, the flow is said to be statistically stationary, or simply stationary. This implies that the correlation tensor depends only on the difference between the two time arguments τ = tt.

B ij (x, t, x , t ) = B ij (x, x , τ ), τ = t -t.
(2.86)

Last, if the statistical properties of the turbulence are invariant under rotations, the flow is said to be statistically isotropic or simply isotropic. In this case, the correlation tensor cannot depend on the direction of r, but only on its magnitude r.

B ij (x, t, x , t ) = B ij (r, t, t ), r = |x -x|.
(2.87)

Isotropy implies homogeneity but not conversely. An example is axisymmetric turbulence, where all statistical properties are invariant under rotations about a specific axis. This flow is necessarily homogeneous but clearly anisotropic. In many cases, global homogeneity and isotropy is an assumption that does not correspond to the actual physics of the flow. A very useful relaxed assumption is that of local homogeneity and isotropy. This assumption is based on the concept of a random field u(x, t) that is not necessarily homogeneous but has homogeneous increments. This means that the statistical properties of the random field

u(x + r, t) -u(x, t) (2.88)
do not depend on the specific point x but only on the vector r, and perhaps the time t. If, in addition, all statistical properties depend only on the magnitude of r and not its direction, the field is called locally isotropic. It is important to note that the statistical properties of a locally homogeneous random field u(x) do depend on x in general; the field is not necessarily homogeneous. Only its increments are guaranteed to have statistics independent of position.

Structure functions

There is no obvious simplification for the moments of a locally homogeneous field, thus it is natural to consider the structure functions, which are moments of the increments

S i1i2...in p1p2...pn (r 1 , r 2 , ...r n ) = (u i1 (x 1 + r 1 ) -u i1 (x 1 )) p1 (u i2 (x 2 + r 2 ) -u i2 (x 2 )) p2 ...(u in (x n + r n ) -u in (x n )) pn .
(2.89) Because of the homogeneity of the increments, the structure functions depend only on the vectors r n and not on the x n . If in addition the field is locally isotropic, the structure function depends only on the magnitudes r n .

The ergodic hypothesis

One can use the assumption of stationarity to replace ensemble averages with time averages, or the assumption of statistical homogeneity to replace ensemble averages with space averages. By performing this replacement, one invokes an ergodic hypothesis which for the case of stationarity is

X = lim T →∞ 1 T T 0 X(t) dT.
(2.90)

Isotropy cannot be used in this sense, because it is not possible perform a rotation by an infinite amount [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]. In the case of turbulence, no proof of the ergodic hypothesis has been achieved, although it is strongly supported by experiment [Monin and Yaglom, 2007a].

Kinematic relations

Under the assumptions of isotropy and incompressibility, it is possible to obtain certain simplified expressions for the moments as well as the structure functions [START_REF] Landau | Fluid Mechanics[END_REF], Monin and Yaglom [2007b].

In particular, the correlation tensor can be expressed as

B ij (r) = (B LL (r) -B N N (r)) r i r j r 2 + B N N (r)δ ij (2.91)
where the longitudinal correlation function is defined by

B LL (r) = u L (x + r)u L (x) (2.92) with u L (x) = u(x)
• r the component of the velocity along the vector r,a signed quantity. Similarly, the transverse correlation function is defined by

B N N (r) = u N (x + r)u N (x) (2.93)
where u N (x) is the component of the velocity component in an arbitrary direction transverse to r. In addition, incompressibility implies

B N N (r) = B LL (r) + r 2 dB LL (r) dr (2.94)
The third order moment5

B ij,l (r) = u i (x)u j (x)u l (x + r) (2.95)
is particularly important because, as it will be seen below, it represents nonlinear energy transfer. In isotropic turbulence, it can be written as

B ij,l (r) = B LL,L (r) -B N N,L (r) -2B LN,N (r) r 3 r i r j r l + B LN,N (r) r (δ jl r i + δ il r j ) + B N N,L (r) r δ ij r l (2.96)
where

B LL,L (r) = u L (x) 2 u L (x + r) , B N N,L (r) = u N (x) 2 u L (x + r) and B LN,N (r) = u L (x)u N (x)u N (x + r) .
(2.97)

The third order moments are related due to the constraint of incompressibility:

B N N,L (r) = - 1 2 B LL,L (r), B LN,N (r) = 1 2 B LL,L (r) + r 4 dB LL,L (r) dr .
(2.98)

The above relations can be obtained under the sole assumptions of isotropy and incompressibility. They are kinematic relations and the Navier-Stokes equation is not used in their derivation [Monin and Yaglom, 2007b]. Two important length scales can be defined in terms of the longitudinal correlation function B LL (r). The first is the integral length scale defined by

L i = ∞ 0 B LL (r) dr B LL (0) . (2.99)
The integral length scale provides an estimate of the maximum spatial extent over which correlations exist in a turbulent flow. A turbulent flow with a large integral length scale exhibits larger-scale correlations in comparison with a flow having a smaller integral length scale. The second length scale is the Taylor microscale defined by

λ = - B LL (0) 2(B LL ) (0) 1 2 (2.100)
It expresses the rate of decay of the longitudinal correlation function near the origin. The smaller the Taylor microscale, the more rapid is the decrease of the correlation function near the origin. The Taylor microscale can also be expressed in terms of the kinetic energy and the mean energy dissipation [Monin and Yaglom, 2007b]

λ = 15νE ε 1 2
.

(2.101)

For a locally isotropic field, the simplified relations for the structure functions are completely analogous to the relations (2.91,2.96) for the moments of an isotropic field. In particular the second order structure function tensor can be written as

S ij (r) = (S L 2 (r) -S N 2 (r)) r i r j r 2 + S N 2 (r)δ ij (2.102)
where

S L 2 (r) = δu L (x, r) = (u L (x + r) -u L (x)) 2 and S N 2 (r) = δu N (x, r) = (u N (x + r) -u N (x)) 2 (2.103)
are the longitudinal and transverse second order structure functions respectively. Due to the incompressibility constraint, they are related by

S N 2 (r) = S L 2 (r) + r 2 dS L 2 (r) dr .
(2.104)

If the turbulence is globally isotropic, the second order structure functions are obviously related to the correlations by

S L 2 (r) = 2(B LL (0) -B LL (r)), S N 2 (r) = 2(B N N (0) -B N N (r)). (2.105)
where

B(0) = B LL (0) = B N N (0) = u 2 (x, t) 3 (2.106)
The third order structure tensor

S ijl (r) = (u i (x + r) -u i (x))(u j (x + r) -u j (x))(u l (x + r) -u l (x)) (2.107)
can be written as

S ijl (r) = (S L 3 (r) -3S LN N 3 (r)) r i r j r l r 3 + S LN N 3 (r) r i r δ jl + r j r δ il + r l r δ ij (2.108)
where

S L 3 (r) = (u L (x + r) -u L (x)) 3 , S LN N 3 (r) (u L (x + r) -u L (x))(u N (x + r) -u N (x)) 2 .
(2.109)

For a solenoidal velocity field,

S LN N 3 (r) = 1 6 S L 3 (r) + r dS L 3 (r) dr (2.110)
Again, the above relations follow directly from the assumptions of incompressibility and isotropy, with no input from the dynamical equations.

Spectral quantities

The definition of the Fourier representation of a homogeneous random field

u(x) = u(k)e ik•x dk (2.111)
presents a particular difficulty: because the definition of a homogeneous field implies that it cannot decay at infinity, the Fourier coefficients are not ordinary random functions but random distributions [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]. An alternative formulation in terms of Fourier-Stieltjes integrals is possible [Batchelor, 1953, Monin andYaglom, 2007b], involving only ordinary random functions, but this formulation is more cumbersome in actual calculations. Here, we will adopt the first approach. The Fourier coefficients are defined by

u(k) = 1 (2π) 3 u(x)e -ik•x dx.
(2.112)

The Fourier transform of the correlation tensor with respect to the separation r defines the spectral tensor

F ij (k) B ij (x) = F ij (k)e ik•r dr. (2.113) Now, B ii (0) = u(x) 2 thus 1 2 u(x) 2 = 1 2 F ii (k) dk (2.114)
The above relation shows that the spectral tensor provides a measure of the energy associated with a particular wavevector k. It is often useful to consider the energy associated with a particular wavenumber k = |k|:

1 2 u(x) 2 = ∞ 0 E(k) dk (2.115)
where the energy spectrum E(k) is defined by integrating the spectral tensor over all angles in wavevector space

E(k) = |k|=k 1 2 F ii (k) dk = 2πk 2 F ii (|k| = k).
(2.116)

In isotropic turbulence, the energy spectrum completely characterizes the spectral energy distribution, as the spectral tensor can be written as

F ij (k) = E(k) 4πk 2 δ ij - k i k j k 2 .
(2.117)

The energy spectrum E(k) can be expressed in terms of the longitudinal correlation function, and vice-versa Monin and Yaglom [2007b]:

E(k) = 1 π ∞ 0 (kr sin kr -k 2 r 2 cos kr)B LL (r) dr, B LL (r) = 2 ∞ 0 - cos kr (kr) 2 + sin kr (kr) 3 E(k) dk.
(2.118)

Similarly, for locally isotropic turbulence, the longitudinal second order structure function can be expressed in terms of the energy spectrum Monin and Yaglom [2007b]

S L 2 (r) = 4 ∞ 0 1 3 + cos kr (kr) 2 - sin kr (kr) 3 E(k) dk.
(2.119)

The above relation implies that if the energy spectrum is a power-law of the form

E(k) = Ck -1-γ (2.120)
then the longitudinal second order structure function is also a power-law:

S L 2 (r) = C r γ . (2.121)
Using the constrains of isotropy and incompressibility, it can be shown that the spectral tensor corresponding the third order moment

F ij,l (k) = 1 (2π) 3 B ij,l (r)e -ik•r dr (2.122)
can be expressed in terms of a single scalar function F 3 (k):

F ij,l (k) = iF 3 (k) δ jl k i k + δ il k j k - 2k i k j k l k 3 . (2.123)
In terms of the energy spectrum, the integral length scale is

L i = π 2 ∞ 0 k -1 E(k) dk ∞ 0 E(k) dk (2.124)
while the Taylor microscale is given by

λ = 3 ∞ 0 E(k) dk 2 ∞ 0 k 2 E(k) dk 1 2
.

(2.125)

Dynamical relations

Starting from the incompressible Navier-Stokes equation (2.1), it is possible to derive equations for the moments such as the correlation tensor [START_REF] Landau | Fluid Mechanics[END_REF], Monin and Yaglom [2007b]. However, because of the quadratic nonlinearity of the governing equations, such an equation for a m-th order moment will contain unknown moments of order m + 1. Thus the system of moment equations always involves more unknown quantities than equations. To find a way to close this system and obtain a solution is known as the closure problem.

For the more general case of homogeneous turbulence, the transport equation for the correlation tensor reads

∂ t B ij (r, t) + ∂ k (B ik,j (r, t) -B i,jk (r, t)) = 1 ρ (∂ i B pj (r, t) -∂ j B ip (r, t)) + 2ν∇ 2 B ij (r, t) (2.126) where B pj (r, t) = p(x + r, t)u j (x, t) (2.127)
is the pressure-velocity correlation. This equation is the general Kármán-Howarth equation for homogeneous turbulence. It describes the evolution of the second order moment B ij (r, t) in terms of the third order moment B ij,k (r, t). The terms involving B ij,k (r, t) represent dissipationless nonlinear energy transfer, whereas the term 2ν∇ 2 B ij (r, t) represents viscous dissipation. In Fourier space, the Kármán-Howarth equation becomes

∂ t F ij (k, t) = Γ ij (k, t) + Π ij (k, t) -2νk 2 F ij (k, t) (2.128) where Π ij (k, t) = 1 (2π) 3 (∂ i B pj (r, t) -∂ j B pi (r, t))e -ik•r dr (2.129)
is the Fourier representation of the terms corresponding to the pressure-velocity correlation. In this case, nonlinear energy transfer is expressed by the tensor Γ ij (k, t) while viscous dissipation is represented by the term 2νk 2 F ij (k, t).

In the case of isotropic turbulence, the general Kármán-Howarth equation can be considerably simplified to the following form

∂ t B LL (r, t) = ∂ r + 4 r B LL,L (r, t) + 2ν∂ r B LL (r, t) . (2.130)
This equation is known as the isotropic Kármán-Howarth equation. Here, nonlinear energy transfer is represented by the third order moment B LL,L (r, t) while viscous dissipation is represented by the term 2ν∂ r B LL (r, t). In spectral space, the isotropic Kárman-Howarth equation assumes the form

∂ t E(k, t) = T (k, t) -2νk 2 E(k, t), T (k, t) = -8πk 2 F 3 (k, t).
(2.131)

In this case, nonlinear energy transfer is expressed in terms of the single scalar function F 3 (k, t). In the above equation, the term -2νk 2 E(k, t) represents viscous dissipation, whereas the term T (k, t) defined in terms of F 3 (k, t) represents nonlinear energy transfer. An important property of the nonlinear energy transfer term is that its integral across the spectrum vanishes [Batchelor, 1953, Monin andYaglom, 2007b]:

∞ 0 T (k, t) dk = 0. (2.132)
Thus the nonlinear energy transfer term only redistributes energy across different Fourier modes without modifying the global energy budget.

Observing that both ∞ 0 2νk 2 E(k) dk and ν ω 2 are equal to the mean dissipation rate ε , we see after multiplying equation (2.131) with k 2 that the term k 2 T (k, t) represents the rate of change of the enstrophy due to nonlinear vortex-stretching interactions. It can be shown [START_REF] Batchelor | The nature of turbulent motion at large wave-numbers[END_REF]] that this term is related to the skewness of the velocity derivative

Σ 3 = (∂ x u) 3 (∂ x u) 2 3/2 = lim r→0 S L 3 (r) (S L 2 (r)) 3/2
(2.133) by the formula

Σ 3 (t) = - 135 98 1/2 ∞ 0 k 2 T (k, t) dk ∞ 0 k 2 E(k, t) dk . (2.134)
This formula reveals the connection between the third order structure function S 3 L (r) and nonlinear energy transfer explicitly. If vortex stretching is present, the right hand side is positive and the skewness of the velocity derivative is negative.

For locally isotropic turbulence, it is possible to express the Kármán-Howarth relation in terms of structure functions [START_REF] Landau | Fluid Mechanics[END_REF]]

- 2 3 ε - 1 2 ∂ t S L 2 (r, t) = 1 6r 4 ∂ r (r 4 S L 3 (r, t)) - ν r 4 ∂ r (r 4 ∂ r S L 2 (r, t)) (2.135) where ε = -∂ t 1 2 u 2 (2.136)
is the mean energy dissipation rate. If the turbulence is stationary, ∂ t S L 2 (r, t) = 0 and the equation can be integrated to give the following relation

S L 3 (r) = - 4 5 ε r + 6ν∂ r S L 2 (r). (2.137)
This relation is known as the 4/5-ths law and was first derived by [START_REF] Kolmogorov | Dissipation of energy in locally isotropic turbulence[END_REF]. It is an exact result that follows from the Navier-Stokes equations under the assumptions of local isotropy and stationarity. If, in addition, the limit of infinite Reynolds number is taken, the term 6ν∂ r S L 2 (r) can be neglected and the law assumes the simple form

S L 3 (r) = - 4 5 ε r (2.138)
The above form of the 4/5-ths law contains the additional assumption of the dissipative anomaly: the mean dissipation rate ε remains finite in the limit of infinite Reynolds number. This important assumption remains unproven, but is supported by both experimental [START_REF] Sreenivasan | The phenomenology of small-scale turbulence[END_REF] and numerical [START_REF] Kaneda | High-resolution direct numerical simulation of turbulence[END_REF] data. An important implication of the 4/5ths law is that, under the assumptions that lead to its derivation, the velocity increment is expected to scale as the third power of the distance

u L (x + r) -u L (x) ∝ |r| 1/3 . (2.139)

Phenomenology

The nonlinear nature of the governing equations makes it difficult to derive exact results about turbulent flows. It is however possible to derive several useful results starting from plausible assumptions without direct use of the governing equations. This phenomenological approach was pioneered by Kolmogorov [1941a] and Obukhov [1941] following the intuitive ideas of [START_REF] Richardson | {Weather Prediction by Numerical Process}[END_REF]. It is based on the concept of the energy cascade: energy is fed in the large scales, and a hierarchy of eddies of ever smaller scale is created by successive instabilities. As large eddies break into smaller eddies, the influence of the energy input and the boundary conditions is forgotten and the turbulent flow tends towards local isotropy. At some point, eddies become small enough for dissipation to be important and the cascade terminates. Thus, energy flows continuously from large scales to small scales, where it is dissipated into heat by viscosity in the dissipative range of scales. For length scales much smaller than the scale where energy is injected, but also much larger than the dissipative range, the influence of both energy injection and dissipation can be neglected. This range is governed by dissipationless nonlinear energy transfer mediated by the inertial terms of the equations of motion, and is known as the inertial range.

Incompressible hydrodynamic turbulence: the Kolmogorov [1941a]-Obukhov [1941] theory (K41)

Let us concentrate on the properties of a particular eddy lying in the inertial range. The notion of an eddy is not a precisely defined mathematical concept but a loosely defined physical entity: it represents a structure of the turbulent flow localized in space, time and scale. Thus, an eddy is characterized by a length scale λ and a velocity scale δu λ , which is a typical velocity difference across the distance λ. From these two quantities we can form a timescale τ λ = λ/δu λ . During the cascade process, an eddy of scale λ in the inertial range breaks into several smaller eddies of size λ < λ under the action of nonlinear instability. These "daughter" eddies go unstable in their turn and break up into the next generation of still smaller eddies. This dissipationless cascade process is repeated until the eddies produced are small enough for the viscosity to damp the instability that induces the breakdown. The cascade process is shown schematically in figure 2.3.

For sufficiently high Reynolds numbers, the separation of scales between large scale energy input and small scale dissipation is so large that it is possible to define the inertial range such that the eddies lying there do not feel the effects of the forcing or the dissipation. In addition, if the inertial range is long enough the cascade process contains a sufficient number of eddy generations for the large-scale details to be forgotten, leading the small scales towards local isotropy. For the inertial range of scales, it is natural to assume that the characteristics of each eddy depend only on the mean energy dissipation rate ε and the length scale λ, and not on the scale of the energy input L or the viscosity ν. This assumption is the content of the first hypothesis of similarity of Kolmogorov [1941a]. The energy associated with an eddy of size λ is (δu λ ) 2 . Assuming that this energy cascades towards smaller scales in the timescale τ λ , we have ε = (δu λ ) 3 /τ λ . Because the energy dissipation rate is, by assumption, a constant quantity independent of scale, we have

δu λ ∝ ε 1/3 λ 1/3 ∝ λ 1/3 . (2.140)
This relation is supported by the 4/5-ths law (2.138), which implies that the velocity increment should scale as the 1/3-power of the distance. In the inertial range, the structure functions scale as

S L p (r) = (δu L (r)) p ∝ ( ε r) p/3 (2.141)
with a similar expression for the transverse structure functions (δu N (r)) p . Because of relations (2.120),(2.121), the energy spectrum scales as

E(k) ∝ ε 2/3 k -5/3 (2.142)
This expression for the energy spectrum was first derived by Kolmogorov [1941a] and Obukhov [1941]. It was later re-derived indepedently by [START_REF] Onsager | The distribution of energy in turbulence[END_REF] and [START_REF] Heisenberg | On the theory of statistical and isotropic turbulence[END_REF]. One can use the prediction (2.140) to define an eddy Reynolds number

Re(λ) = δu λ λ ν ∝ λ 4/3 (2.143)
which, as expected, decreases as the scale λ decreases: smaller eddies are more and more influenced by viscosity. At the scale l d where Re(l d ) ∼ 1, viscosity is expected to dominate. This scale can be taken as a representative scale of the dissipative range

l d = ν 3 ε 1 4
.

(2.144)

The above estimate can provide information on the number of active degrees of freedom in a turbulent flow.

Assuming that the largest scale is the integral length scale L i , we have

N = L i l d 3 ∝ Re 3/4 . (2.145)
Relation (2.144) for the dissipative scale can equivalently be derived under the assumption that in the dissipative range of scales, the properties of the turbulence should depend only on ε and ν. This assumption forms the content of the second similarity hypothesis of Kolmogorov [1941a].

Incompressible magnetohydrodynamic turbulence

The Iroshnikov-Kraichnan theory

In incompressible magnetohydrodynamic turbulence, the eddy concept is replaced by the concept of an Alfvén wave. Two colliding Alfvén waves are deformed by the nonlinear term in such a way that total energy and cross-helicity remain constant, since these quantities are only dissipated by the diffusive terms. This property allowed [START_REF] Iroshnikov | Turbulence of a conducting fluid in a strong magnetic field[END_REF] and [START_REF] Kraichnan | Inertial-range spectrum of hydromagnetic turbulence[END_REF] (IK) to construct a phenomenological theory of the inertial range. They considered the Elsässer equation (2.42) in the presence of a mean magnetic field B 0 :

∂ t z ± ∓ v A • z ± + z ∓ • ∇ z ± = -∇ P + ν + ∇ 2 z ± + ν -∇ 2 z ∓ (2.146)
where v A = B 0 / √ 4πρ is the Alfvén velocity defined in terms of the mean field. The mean field B 0 can be either a global imposed field or a local average.

If the Alfvén wave amplitude at a scale λ is δu λ , during a collision of two waves it changes by the magnitude of the nonlinear term δu 2 λ /λ multiplied by the interaction time λ/V A : ∆δu λ = δu 2 λ /V A . If one further assumes that collisions change the energy of the waves as a random walk, the number of collisions required to deform each wave considerably is

N = (δu λ /∆δu λ ) 2 = (V A /δu λ ) 2 . The timescale of the cascade is τ IK = N λ/V A = (λV A )/δu 2 λ .
Assuming a constant energy flux through the inertial range ε = δu 2 λ /τ IK , one obtains the scaling δu λ ∝ ( ε B 0 λ) 1/4 . For the energy spectrum, IK predicts

E(k) ∝ C IK ( ε B 0 ) 1/2 k -3/2 while the dissipative scale is l IK d = (B 0 η 2 / ε ) 1/3 .

The Goldreich-Sridhar theory

In the IK theory, the energy spectrum is isotropic. [START_REF] Goldreich | Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence[END_REF] (GS95) introduced the assumption of critical balance: the Alfvén term ∓(v A • ∇)z ± is of the same order as the nonlinear term (z ∓ • ∇)z ± . This implies V A /l ∼ δb λ /λ with l the longitudinal scale of the wave and λ its transverse scale. l is equal to the product of the Alfvén velocity V A and the eddy turnover time τ GS = λ/δu λ , implying δu λ ∼ δb λ . τ GS is chosen as the timescale of the cascade, and the constancy of the energy flux δu 2 λ /τ GS gives δu λ ∼ δb λ ∝ λ 1/3 in the inertial range. This scaling is identical to K41, except that the scale λ is perpendicular to the magnetic field, resulting in an anisotropic energy spectrum

E(k ⊥ ) ∝ ε 2/3 k -5/3 ⊥ . In fact, critical balance implies l ∝ ε -1/3 B 0 λ 2/3
, so that the cascade becomes increasingly anisotropic at small scales.

Compressible turbulence

From the point of view of phenomenology, an important difference between incompressible and compressible turbulence is the absence, in the compressible case, of a quadratic conserved quantity which would allow the construction of a cascade argument. In incompressible turbulence this role is played by the kinetic energy 1/2 u 2 which is exactly conserved by nonlinear interactions involving triads of wavenumbers. In compressible turbulence the conserved quantity is the total energy e + 1/2ρu 2 which involves the internal energy of the gas and cannot be obviously related to the statistics of the velocity field. Another complicating factor is that nonlinear interactions in compressible turbulence are fundamentally more complex due to the presence of the density field. Even in the case of isothermal turbulence, where the conserved quantity is 1/2ρu 2 , the cubic nonlinearity due to the presence of the density field complicates things considerably. Theoretical results in compressible turbulence are even scarcer than incompressible turbulence. [START_REF] Galtier | Exact relation for correlation functions in compressible isothermal turbulence[END_REF] have derived an exact expression, analogous to the 4/5ths law, for the third order structure function in isothermal turbulence. [START_REF] Banerjee | Exact relation with two-point correlation functions and phenomenological approach for compressible magnetohydrodynamic turbulence[END_REF] have generalized this result for isothermal magnetohydrodynamic turbulence.

Intermittency

All the phenomenological theories presented in section 2.4 are built on the assumption that the energy dissipation rate ε(x) is uniformly distributed in space. Both experiments [START_REF] Batchelor | The nature of turbulent motion at large wave-numbers[END_REF][START_REF] Grant | Turbulence spectra from a tidal channel[END_REF][START_REF] Anselmet | High-order velocity structure functions in turbulent shear flows[END_REF] and numerical simulations [START_REF] Siggia | Numerical study of small-scale intermittency in three-dimensional turbulence[END_REF][START_REF] Kerr | Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF] have shown that this is not the case: in fact, the energy dissipation rate is concentrated on a few intense structures separated by large quiescent regions, a property known as intermittency. Figure 2.4 shows an example of an intermittent function: the viscous energy dissipation rate along a fixed line, taken from a high Reynolds number incompressible magnetohydrodynamics simulation. The dissipation rate is clearly concentrated in a small subset of the line, assuming high values there. The property of intermittency has obvious important practical consequences for the heating in a turbulent flow, but it also has important theoretical consequences for the general scaling properties of the turbulence.

Self-similarity and intermittency

Let us focus on the longitudinal component of a turbulent velocity field u L (x). The property of intermittency can be quantified if one considers the random variable

θ(r) = δu L (x, r) ( ε r) 1/3 = u L (x + r) -u L (x) ( ε r) 1/3 (2.147)
where homogeneity and isotropy of the velocity field have been assumed. Under the assumptions leading to K41, the probability density function of the random variable θ(r) should be independent of the scale r: the velocity field is self-similar. Perhaps the most well-studied example of a self-similar random function is Brownian motion, a realization of which is presented in figure 2.5. In the figure, it is intuitively clear that the statistical properties of the curve do not depend on scale. For an intermittent velocity field however, this is no longer the case: as smaller separations r are considered, the tails of the distributions corresponding to extreme velocity bursts become increasingly important. For example in many cases the distribution p(δu L (x, r)) is very close to Gaussian for large r but assumes a very different form, such as a stretched exponential for small r. This can be seen clearly in figure 2.6, where the pdf of the longitudinal velocity increment from a high resolution incompressible hydrodynamical simulation is shown for various separations. For a function that exhibits no intermittency, the non-dimensional ratios6

F (r) = S 2p L (r) (S p L (r)) 2
(2.148) are independent of scale: there is no statistical difference between large and small scales. For an intermittent function, by contrast, there is by definition much more extreme bursts in the small scales than in the large scales, and thus the ratios F (r) are expected to grow with decreasing r:

F (r) ∝ L i r δ . (2.149)
This is because the structure functions S 2p (r) in the nominator are of higher order than those in the denominator, and thus sample the tails of the probability distribution of the velocity increment corresponding to extreme events. For an intermittent function, these tails become stronger and stronger at small increments leading to a faster increase of the quantity S 2p (r) in comparison with (S p (r)) 2 . Thus the introduction of the property of intermittency breaks the self-similarity of the velocity field by introducing the integral length scale L i in the scaling of the statistics.

The phenomenological theories presented in section 2.4 are clearly self-similar. In K41 for example, the structure functions scale as S L p (r) ∝ r p/3 (2.150) which implies that F (r) is constant. These phenomenological theories are therefore incompatible with intermittency.

The refined similarity hypothesis

The measures of intermittency of the velocity field, expressed by the scaling of the structure functions S L p (r), can be related to the intermittency characteristics of the dissipation field ε(x). Obukhov [1962] suggested that ε , the mean energy dissipation rate could be replaced by the quantity

ε r = 1 4/3πr 3 B ε(x) dx (2.151)
which is the energy dissipation rate averaged over a sphere B of radius r. One can use this quantity to define the time scale T r = r 2/3 ε -1/3 r , the velocity scale U r = (rε r ) 1/3 , the scale-dependent Reynolds number Re r = U r/ν = r 4/3 ε 1/3 r /ν and the dissipation scale l r d = ν 3/4 ε -1/4 r

. [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF] went further by introducing a refined similarity hypothesis which states the statistical properties of the non-dimensional velocities

w(ξ k , τ k ) = u(ξ k r, τ k T r ) U k (2.152)
where |ξ k | and |τ k | are of order unity, depend only on Re r . In additon, if Re r → ∞ then the aforementioned statistical properties become universal. It follows from this hypothesis that for l r d r L i , the structure functions can be expressed as

S L p (r) = C p ε p/3 r r p/3 (2.153)
where the C p are constants of order unity.

The refined similarity hypothesis allows to bridge the scaling exponents of the dissipation field and the scaling exponents of the structure functions.

S L p (r) ∝ ε p/3 r ζ L p , ζ L p = p/3 + τ p/3 (2.154)
where τ p is the scaling exponent of the dissipation field defined by

ε p r ∝ ε p r L i τp .
(2.155)

It must be stressed that the introduction of relation (2.155) and the scaling exponent τ p is a second assumption of extended self-similarity, independent from assumption (2.152).

Scaling properties of the dissipation field

The refined similarity hypothesis suggests that the K41 theory is only a first approximation because it ignores fluctuations in the energy dissipation rate. Thus K41 is only a mean field theory of turbulence and fluctuations in the dissipation rate will introduce at least one more universal scaling exponent. In this context, the most simple relevant correlation is

ε(x + r)ε(x) = 1 2 d 2 (r 2 ε r ) dr 2 (2.156)
where the right hand side follows from homogeneity and the definition of ε r . For r in the inertial range of scales, it is natural to assume that this correlation function scales as a power of r. Assuming further that the viscosity plays no role at these scales, we expect, because of assumption (2.155)

ε(x + r)ε(x) ∝ ε 2 L i r µ .
(2.157)

Thus, the fluctuations of the dissipation field introduce a new scaling exponent µ. For incompressible hydrodynamic turbulence, the value of µ found in experiments and numerical simulations is 0.20 ± 0.05 [START_REF] Anselmet | Turbulent flows and intermittency in laboratory experiments[END_REF].

As an example of how the exponent µ can enter the scaling of the velocity field, let us suppose that indeed in the inertial range ε r scales as (2.155) and the structure functions S L p (r) scales as (2.154). Then, the refined similarity hypothesis implies that ζ L 6 = 3µ. The scaling properties of the energy dissipation rate can be characterized in more details by introducing the notion of multifractality [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]Parisi, 1985, Paladin and[START_REF] Paladin | Anomalous scaling laws in multifractal objects[END_REF]. The energy dissipation field is said to be multifractal if it exhibits the scaling

ε r (x) ε ∝ r L i α-1
(2.158) for all points x lying in a set D α which has fractal dimension dimD α = F (α). Thus, a multifractal dissipation field has singularities of exponent α-1 on sets of dimensions F (α). The function F (α) is called the spectrum of singularities of the multifractal field. For a three-dimensional turbulent flow, the probability of a randomly chosen sphere centered at the point x with a radius equal to r to lie in a given set D α scales as [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]

P(x ∈ D α ) ∝ r L i 3-F (α)
.

(2.159)

The quantity 3 -F (α) is called the co-dimension of the set D α . This implies that the moment of the dissipation is given by the integral

ε q r = ε q αmax αmin r L i q(α-1)+3-F (α)
dα.

(2.160) By a steepest-descent argument [START_REF] Bender | Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory[END_REF], the above integral yelds the dissipation scaling exponent τ q : τ q = min α {q(α -1) + 3 -F (α)}.

(2.161) An important consequence of multifractality which is apparent in the above relation is that the scaling exponent τ q becomes in general a nonlinear function of the order q. Through the bridge relation (2.154), this implies that the longitudinal structure function scaling exponents ζ L p are also a nonlinear function of p. This is to be contrasted with the phenomenological theories of section 2.4 which all predict a linear relation between p and ζ L p , but also with the unifractal beta model of [START_REF] Frisch | A simple dynamical model of intermittent fully developed turbulence[END_REF]. The above definition of multifractality assumes that singularities are indeed present in the dissipation field, a fact that is far from obvious based on the current understanding of the Navier-Stokes equation [START_REF] Gallavotti | Some rigorous results about 3d navier-stokes[END_REF][START_REF] Foias | Systems of conservation equations with a convex extension[END_REF]. One can circumvent this difficulty by adopting the following probabilistic definition of multifractality [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]: the dissipation field is said to be multifractal if there is a function F (α) which maps real scaling exponents α to scaling dimensions F ≤ 3 such that for any α we have

lim r→0 ln C r (r α-1 ) ln r = 3 -F (α) (2.162)
where C r (ε) is the cumulative distribution function of the energy dissipation rate ε r . In addition to avoiding the unnecessary assumption of the existence of singularities in the dissipation field, the above probabilistic definition leads to the interpretation of intermittency in terms of large deviations theory [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF], which will not be presented here.

Models of intermittency

In this subsection, we will present specific intermittency models which have been useful either conceptually or in analyzing and interpreting data from experiments, observations and numerical simulations. All these models are purely phenomenological, as they do not deduce intermittency properties directly from the governing equations but use instead plausible physical assumptions for the cascade process. For a more complete review of intermittency models the reader is referred to the second volume of Monin and Yaglom [2007b] and to the book of [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF].

Random cascade models

Let us consider a turbulent flow inside a cube of side l 0 in which the energy dissipation rate is uniform, equal to ε . In a random cascade model, this uniform average value represents the "mother eddy" of generation n = 0. The eddies of the next generation n = 1 are obtained by subdividing the cube into eight equal sub-cubes of side l 1 = l 0 /2. In each new sub-cube, the value of the energy dissipation rate is defined as the product of the value corresponding to the previous generation with a chosen random variable W , which has the following properties:

W ≥ 0, W = 1, W q < ∞, q > 0.
(2.163)

The nth generations consists of 2 n eddies of size

l = l 0 2 -n . (2.164)
Each of these eddies is characterized by a uniform value of the energy dissipation rate which is equal to

ε l = ε n i=1 W i (2.165)
where the random variables W i are assumed independently and identically distributed.

The property W i = 1 implies that the mean value of each of the ε l s is equal to ε . Nevertheless, the sum of the dissipations of the eight daughter eddies ε l/2 is not directly related to the value of the dissipation of the mother eddy ε l : the cascade is nonconservative. In addition, multiplication of a great number of random variables can lead to very large fluctuations in the small scales.

Using relations (2.164) and (2.165), we can calculate the moments of the random variable ε l :

ε l = ε q l l 0 τq (2.166)
where the exponent τ q is given by

τ q = -log 2 W q .
(2.167)

A particularly simple case of a random cascade model is the black and white model of [START_REF] Novikov | The intermittency of turbulence and the spectrum of energy dissipation fluctuations[END_REF]. In this model, only two values are permitted for the random variable W :

W = 1/β with probability β 0 with probability 1 -β (2.168)
This implies

τ q = -(1 -q) log 2 β, ζ p = p 3 -1 - p 3 log 2 β.
(2.169)

The [START_REF] Novikov | The intermittency of turbulence and the spectrum of energy dissipation fluctuations[END_REF] model is precisely equivalent to the unifractal β model of [START_REF] Frisch | A simple dynamical model of intermittent fully developed turbulence[END_REF], predicting the same linear relation between the scaling exponents and the order of the moment. However, if the random variable W takes more than one non-vanishing value, one can show using the Schwartz inequality that the model produces a nontrivial concave relation between τ p and p. Hence, random cascade models exhibit multifractal scaling properties in general.

The log-normal model

The first intermittency model considered by [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF], is a random cascade model for which the probability density function of the random variable W is assumed to be log-normal. This assumption follows from considering the logarithm of the dissipation ε l in the context of a random cascade model:

ln ε l = ln( ε n i=1 W i ) = ln ε + n i=1 ln W i .
(2.170)

Hence for small scales l the logarithm of ε l is equal to the sum of a large number of independent identically distributed random variables ln W i . By the central limit theorem, the random variable ln ε l will be normally distributed and therefore ε l will have a log-normal probability density function.

To illustrate the predictions of the log-normal model, let us assume that

W = 2 -m (2.171)
where m is a Gaussian random variable with mean m = m (2.172) and variance

σ 2 = (m -m) 2 .
(2.173)

The relation

W = 1 implies 2 m = σ 2 ln 2 (2.174)
Observing that the sum

n i=1 m i (2.175)
is also a Gaussian random variable, one can use the formalism of large deviations theory [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF] to compute the singularity spectrum of the dissipation field:

F (α) = 3 - (α -1 -µ/2) 2 2µ (2.176)
where the parameter µ defined by µ = 2 m (2.177)

is the same quantity that appears in the scaling of local energy dissipation rate ε r in (2.157). From relation (2.161) and the bridge relation (2.154), the values of the scaling exponents predicted by the log-normal model are

τ q = µ 2 (q -q 2 ), ζ L p = p 3 + µ 18 (3p -p 2 ).
(2.178)

From the above relations we see that the second order moment of the dissipation scales as l -µ , and thus indeed the quantity µ appearing in equation (2.157) is equal to 2 m in the context of the log-normal model. We see that the log-normal model predicts a nonlinear relation between the scaling exponents and the order of the moment, and thus exhibits multifractal scaling. However, for

p > p * = 3 2 + 3 µ (2.179)
ζ L p becomes a decreasing function of the order p. It can be shown that this is a serious inconsistency which invalidates the incompressibility assumption [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]. In addition, the log-normal model violates an inequality first proven by [START_REF] Novikov | Intermittency and scale similarity in the structure of a turbulent plow[END_REF] which states that τ q + 3q ≥ 0 q ≥ 0,

τ q + 3q ≤ 0 q ≤ 0.
(2.180)

The physical reason for which the log-normal model violates the Novikov inequality is the nonconservative nature of the cascade [START_REF] Mandelbrot | Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence[END_REF]. Another property of the log-normal model is the fact that the dimension F (α) becomes negative for α > 1 + µ 2 + 6µ.

(2.181) However, this is not an inconsistency since it can be interpreted as simply meaning that the co-dimension 3 -F (α) is greater than 3, implying that the probability p l of encountering a scaling exponent α -1 goes to zero with l faster than l 3 [START_REF] Mandelbrot | Negative fractal dimensions and multifractals[END_REF][START_REF] Mandelbrot | Random multifractals: negative dimensions and the resulting limitations of the thermodynamic formalism[END_REF]. According to [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF], the inconsistencies of the log-normal model stem from the invocation of the central limit theorem, which completely ignores large deviations.

The log-Poisson model

A very successful intermittency model was introduced by She and Leveque [1994b] for the case of incompressible hydrodynamic turbulence. The model was later generalized by [START_REF] Grauer | Scaling of high-order structure functions in magnetohydrodynamic turbulence[END_REF] and [START_REF] Politano | Model of intermittency in magnetohydrodynamic turbulence[END_REF] for incompressible magnetohydrodynamic turbulence and by [START_REF] Boldyrev | Scaling relations of supersonic turbulence in star-forming molecular clouds[END_REF] for compressible hydrodynamic turbulence. [START_REF] Dubrulle | Intermittency in fully developed turbulence: Log-poisson statistics and generalized scale covariance[END_REF], [START_REF] She | Quantized energy cascade and log-poisson statistics in fully developed turbulence[END_REF] reinterpreted the model in terms of log-Poisson statistics. Apart from being successful in reproducing experimental, observational and numerical results, this model is unique in incorporating geometrical information about dissipative structures.

In the original formulation of the log-Poisson model by She and Leveque [1994b], the energy dissipation field averaged over a ball of radius l ε l is characterized by a hierarchy of fluctuation structures ε (p) l , p = 0, 1, 2, ... which are defined by the ratio of successive dissipation moments

ε (p) l = ε p+1 l ε p l .
(2.182)

The first structure of the hierarchy ε (0) l is corresponds to the mean energy dissipation rate ε and the last structure ε (∞) l corresponds to the small scale dissipative structures, assumed to be vortex filaments. The intensity of the pth order structures ε (p) l is given by

ε (p) l = ε p+1 l ε p l = ε p+1 l P (ε l ) dε l ε p l P (ε l ) dε l = ε l Q p (ε l ) dε l (2.183) so that ε (p) l → ε (∞) l
as p → ∞. The A p are scale-independent non-universal constants. Equations (2.183) and (2.189) imply

τ p+2 -(1 -β)τ p+1 + βτ p + 2 3 (1 -β) = 0 (2.190)
a relation which is consistent with the asymptotic relation (2.188). Setting τ p = -2/3p + 2 + f (p) with lim p→∞ f (p) = 0 we obtain a homogeneous difference equation with solution

f (p) = αβ p .
(2.191)

The solution must satisfy the boundary conditions τ 0 = 0, which follows from the assumption of finite dissipation in the limit of zero viscosity, and τ 1 = 0 which follows from the 4/5ths law (2.138) and the bridge relation (2.154). These two boundary conditions uniquely determine α and β, leading to the final expression

τ p = - 2 3 p + 2 1 - 2 3 p (2.192)
for the scaling exponents of the dissipation moments and

ζ L p = p 9 + 2 1 - 2 3 p/3 . (2.193)
For the energy spectrum, the model predicts

E(k) ∝ k 29/9+2(2/3) 2/3 k -5/3-0.03 (2.194)
while for the intermittency exponent µ the predicted value is µ = 2/9.

(2.195)

The log-Poisson model has met considerable success in predicting experimental, observational and numerical results [START_REF] Anselmet | Turbulent flows and intermittency in laboratory experiments[END_REF]. Three important ingredients that render the model unique are the postulated scaling relation (2.189) which attempts to capture the nonlinear dynamics of the cascade, the consideration of the filamentary character of the most intense dissipative structures and the assumption of a divergent scaling law for the energy flux. This success of the log-Poisson model is illutrated in figure 2.7, which shows a comparison of experimental data with the predictions of the structure function scaling exponents of various intermittency models.

Interpretation in terms of log-Poisson statistics

The above model admits an alternative interpretation in terms of log-Poisson statistics [Dubrulle, 1994, She and[START_REF] She | Quantized energy cascade and log-poisson statistics in fully developed turbulence[END_REF]. This interpretation is formulated in terms of the non-dimensional energy dissipation

π l = ε l ε (∞) l .
(2.196)

The scaling relation (2.189) implies that

π l = B p B (1-β p )/(1-β) 1 π l (1-β p )/(1-β) (2.197)
where the B p are non-universal constants. This is a relation between moments that is satisfied by a whole family of probability distribution functions. Introducing the variable Y = log π l / log β, these functions can be expressed as and Leveque [1994b]. The Poisson distribution is known to be associated with the general concept of intermittency, as it appears in the context of discontinuous random processes [Skorokhod and [START_REF] Gikhman | Introduction to the theory of random processes[END_REF], and is also the limit of a wide class of statistical distributions involving rare events [START_REF] Arratia | Poisson approximation and the chen-stein method[END_REF].

P λ(l) (π l ) dπ i = F λ(l) dY, F λ(l) (Y ) = P λ(l) G(Y -Z) dZ (2.198)

Extension to incompressible magnetohydrodynamics

Grauer et al. [1994] and [START_REF] Politano | Model of intermittency in magnetohydrodynamic turbulence[END_REF] have extended the log-Poisson model for the case of incompressible magnetohydrodynamic turbulence. Taking into account the inertial interactions of Alfvén waves, the refined similarity hypothesis (2.152) assumes the form

ε l ∝ (δz(l)) g lB g+3 0 (2.200)
where δz(l) is the Elsässer increment, g is the inverse of the scaling exponent in the inertial range defined by δz(l) ∝ l 1/g and B 0 is the quasi-uniform field that guides the Alfén waves. For IK phenomenology, g = 4 while for K41 g = 3. The bridge relation (2.154) remains identical.

Assuming the scaling laws

ε (p+1) l = A B p ((ε l ) (p) ) β B ((ε l ) (∞) ) 1-β B , 0 < β B < 1 (2.201)
for the dissipation hierarchy and

t l ∝ l 1-1/g (2.202)
for the time scale of the cascade, we obtain

τ p = - p g + (C B 0 ) 1-p 1 - 1 g p (2.203)
for the dissipation exponents and

ζ B p = p g 1 - 2 g + C B 0 1 -1 - 2 gC B 0 p g (2.204)
where the co-dimension C B 0 is given by

C B 0 = 1 -g 1 -β B .
(2.205)

Thus the model contains two free parameters: the inertial range scaling exponent g and the co-dimension of the dissipative structures C B 0 . In magnetohydrodynamicturbulence, the co-dimension C B 0 = 1 as most of the dissipation occurs in quasi two-dimensional sheet-like structures. The exponent parameter g can be either 3 for K41 scaling or 4 for IK scaling. This generalization of the log-Poisson model has met considerable success in predicting results of numerical simulations [START_REF] Müller | Scaling properties of three-dimensional magnetohydrodynamic turbulence[END_REF] and solar wind observations [START_REF] Politano | Model of intermittency in magnetohydrodynamic turbulence[END_REF].

Extension to compressible turbulence [START_REF] Boldyrev | Scaling relations of supersonic turbulence in star-forming molecular clouds[END_REF] assumed that the inertial range of compressible magnetohydrodynamic turbulence is governed by K41 scaling leading to the choice g = 3. The co-dimension is taken equal to 1 because the dissipation is assumed to take place in shocks, which are two dimensional structures. This "Kolmogorov-Burgers" model, which is identical with the model used in the numerical study of [START_REF] Müller | Scaling properties of three-dimensional magnetohydrodynamic turbulence[END_REF] in the context of incompressible magnetohydrodynamics, has been successfully compared against numerical simulations of isothermal hydrodynamic and magnetohydrodynamic compressible turbulence [START_REF] Padoan | Structure function scaling in compressible superalfvénic mhd turbulence[END_REF][START_REF] Kritsuk | The Statistics of Supersonic Isothermal Turbulence[END_REF].

Chapter 3

Numerical Methods

Turbulence in the ISM can exhibit a wide range of sonic and Alfvénic Mach numbers [START_REF] Elmegreen | Interstellar turbulence i: Observations and processes[END_REF]. The high variability of the sonic Mach number, which in general lies in the interval [0.1, 10], implies that there are regions where the flow is highly compressible and shock-dominated (high Mach number), and regions where the flow is essentially incompressible (low Mach number). This dichotomy calls for different numerical methods in order to capture the different physics of incompressible and compressible flows.

In the case of incompressible flow, spectral methods introduced by Orszag and Patterson [1972] (see also [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF]) are the standard because of their high accuracy for the representation of the small scales. However, their oscillatory character near discontinuities prohibits their extension to compressible flow. An efficient class of methods for strongly compressible, shock-dominated flows are the Godunov methods introduced by [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] (see also [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]).

In this chapter, we will present the principles of spectral and Godunov methods, which were used in the numerical simulations of compressible and incompressible ISM turbulence.

Spectral Methods

Principles

In a spectral method, the velocity and magnetic fields are developed in a truncated Fourier series (trigonometric polynomial) using the Fast Fourier Transform (FFT) algorithm. In the following, we will use the velocity field as an example, but the case of the magnetic field is completely analogous.

The velocity field is defined on a regular Cartesian grid of equidistant points with coordinates

C x = 2π N (n 1 , n 2 , n 3 ), 0 ≤ n j < N, j = 1, 2, 3 . (3.1)
Boundary conditions are, by definition, periodic in all directions while the length of the computational domain is L = 2π. Because of periodicity, one can develop the velocity field in a truncated Fourier series with coefficients

u(k) = 1 N 2 k∈C k u(x) exp(-ik • x) (3.2)
where k is the discrete wave-vector belonging to the set

C k = {(k 1 , k 2 , k 3 ), -K ≤ n j < K, j = 1, 2, 3} (3.3)
with K = N/2. The velocity field can be reconstructed from the Fourier coefficients u(k) as

u(x) = k∈C k u(k) exp(ik • x). (3.4)
Spatial derivatives are mapped in Fourier space into algebraic operations (multiplication by the wave-vector

k) ∇φ(x) → ikφ(k), ∇ • u(x) → ik • u(k), ∇ 2 u(x) → -|k| 2 u(k). (3.5) A quadratic nonlinear term U ij (x) = u i (x)u j (x) is transformed in Fourier space into a convolution sum U ij (k) = 1 N 3 k=p+q p,q∈C k u i (p)u j (q), k ∈ C k . (3.6)
Evaluation of this sum in Fourier space requires O(N 6 ) arithmetic operations, a prohibitive cost. One can use the convolution theorem to compute this sum in physical space. Practically, an inverse Fourier transform of u i (k) and u j (k) is performed, U ij (x) = u i (x)u j (x) is calculated in physical space and finally a forward Fourier transform gives the required U ij (k). The cost of this procedure is O(N 3 log N ) operations, essentially the cost of the three-dimensional FFT.

In the discrete Fourier space, the Navier-Stokes equation becomes

(∂ t + νk 2 )u(k) = NL(k). (3.7)
The nonlinear term NL(u(k)), which includes the pressure, is calculated by returning to physical space. In fact, it is the only point in the algorithm where it is necessary to return to physical space, as the above equation can be advanced in Fourier space, using a standard integration algorithm such as the Runge-Kutta method.

Dealiasing

However, straightforward application of the fast Fourier transform introduces "aliasing" errors. The fast Fourier transform, being a discrete transform, is unable to distinguish between wave-vectors components modulo 2K, where K is the maximum wave-vector component considered in the calculation. Therefore, what is actually obtained by applying the inverse transform to the nonlinear term U ij (x) = u i (x)u j (x) is not U ij (k) but rather the aliased sum

A ij (k) = 1 N 3 k=p+q p,a(q)∈C k u i (p)u j (a(q)), k ∈ C k (3.8)
where a i (q) = q i -2K if q i ≥ K,a i (q) = q i + 2K if q i < -K and a i (q) = q i otherwise.

The terms involving a(q) = q are the ones who introduce the aliasing error and are not desirable. They can be avoided altogether by restricting the wave-vector domain C k . The most popular dealiasing scheme is the two-thirds rule, which consists of a truncation of all wave-vectors that have a component greater that

N/3 C 2/3 k = {k ∈ C k , k i < N/3, i = 1, 2, 3}.
(3.9)

Using this scheme, the quadratic interaction of the two largest possible wave-vector components, equal to N/3, excites a wave-vector component equal to 2N/3, which gets aliased at -N/3-1, just outside the allowed range {-N/3, N/3}. This dealiasing scheme is very efficient computationally, as it does not demand extra operations, but has a low efficiency in representing the small scales of the flow due to the severe truncation: in a three-dimensional simulation, the set C 2/3 k represents only 31% of the Fourier modes in the set C k . Thus, less than one third of the Fourier modes remain active, the rest being kept equal to zero throughout the computation. Patterson and Orszag [1971b] proposed an alternative method for dealiasing, which involves the wavevector set

C P O k = {k ∈ C k , |k| < √ 2 3 K}. (3.10)
The above truncation is however not sufficient to remove the aliasing errors completely, and an additional evaluation of the nonlinear term in a grid shifted by ∆x/2 is required. In practice, this shift is efficiently implemented as a multiplication by a phase factor exp(i∆x/2) in Fourier space followed by an inverse Fourier transform. The method of Patterson and Orszag [1971b] introduces the extra cost of an additional evaluation of the nonlinear term, but the set C P O k corresponds to 45% of the modes contained in the set C k . In [START_REF] Orszag | Numerical simulation of incompressible flows within simple boundaries: accuracy[END_REF], it is pointed out that if one gives up the isotropy feature of the set C P O k , it is possible to extend it by considering the set

C O k = {k ∈ C k , |k i ± k j | ≤ 2N 3 , i, j = 1, 2, 3, i = j}. (3.11)
This polyhedral truncation is sufficient to remove all aliasing errors if it is complemented by the above phase-shift method. It is not isotropic in k-space, but the set C O k corresponds to 52% of the modes in the set C k , giving a better representation of the small scales.

The case of a cubic term

In ambipolar diffusion magnetohydrodynamics, the induction equation includes the cubic ion-neutral drift term ∇ × ((J × B) × B).

(3.12)

In [START_REF] Hossain | On computing high order galerkin products[END_REF], it is shown that computing (J × B) first, dealiasing it using a standard second order method, and then multiplying by B before dealiasing again is not sufficient to remove aliasing errors.

If one chooses to extend the standard two-thirds rule, it is easy to see that in order to remove all aliasing errors from the cubic term, one has to consider the set

C * k = {k ∈ C k , k i < N/2, i = 1, 2, 3}. (3.13)
This is a very severe truncation, as in three dimensions the above set contains only 14% of the modes in the set C k .

As is shown in [START_REF] Hossain | On computing high order galerkin products[END_REF] however, one can use either the isotropic truncation of Patterson and Orszag [1971b] or the polyhedral truncation of [START_REF] Orszag | Numerical simulation of incompressible flows within simple boundaries: accuracy[END_REF], complemented by the phase-shift method, to completely dealias a cubic term without further additions or modifications in the scheme. Thus ambipolar diffusion magnetohydrodynamics can be treated with the exact same numerical scheme as ordinary magnetohydrodynamics, although there is a higher order nonlinearity in the induction equation.

Time integration

The application of the spectral methods results in a set of coupled ordinary differential equations for the Fourier coefficients u(k)

du(k) dt = NL(u(k)) -ν|k| 2 u(k) (3.14)
The linear terms (that is, viscous and Ohmic terms) are advanced semi-implicitly using the Crank-Nicholson scheme [START_REF] Hirsch | Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics[END_REF]

du(k) dt = u i+1/2 (k) -ν|k| 2 ( 1 2 u i+1 (k) + 1 2 u i (k)) (3.15)
where the index i represents the value at timestep i. The term u i+1/2 (k) is given by explicit advancement of the nonlinear term using a 4th order Runge-Kutta method [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF]

u i+1 (k) = u i (k) + ∆t 6 (k 1 + 2k 2 + 2k 3 + k 4 ). (3.16)
The coefficients k i are given by the expressions

k 1 = NL(u i (k)), k 2 = NL(u i (k)) + ∆t 2 k 1 ), k 3 = NL(u i (k)) + ∆t 2 k 2 ), k 4 = NL(u i (k)) + ∆tk 3 ).
(3.17)

We have found that the introduction of the Crank-Nicholson scheme for the linear terms increases overall stability in comparison with a pure Runge-Kutta scheme.

The timestep is determined by the Courant-Friedrichs-Lewy (CFL) condition [START_REF] Hirsch | Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics[END_REF] ∆t

≤ CFL ∆x C (3.18)
where C is the maximum velocity of propagation of information. In the case of incompressible hydrodynamics C is given by the spatial maximum of the velocity

C HD = max x |u(x)| (3.19)
whereas is the case of incompressible magnetohydrodynamics, it is given by the spatial maximum of the Alfvén velocity, in the frame where the fluid is locally at rest

C M HD = max x |u(x) ± b(x)|. (3.20)
For a stable time integration scheme, CFL must be a number < 1. In practice, we have used a value of 0.8 for magnetohydrodynamic simulations, whereas a value of 0.1 was necessary for ambipolar diffusion simulations.

The code ANK

The spectral method described above was implemented in the ANK1 code. I wrote the ANK code in the Fortran 95 programming language, employing a hybrid approach for parallelization, using a combination of distributed and shared memory parallelism which allowed me to efficiently utilize cluster architectures of multi-core processors. Distributed parallelism is achieved using the Message Passing Interface (MPI), while shared memory parallelism is achieved using OpenMP. The code relies on the FFTW2 free software library for the computation of Fast Fourier Transforms. The scaling of the code has been tested for up to 256 processor cores with excellent results. ANK is capable of simulating incompressible hydrodynamics, incompressible magnetohydrodynamics and incompressible ambipolar diffusion magnetohydrodynamics (using the single-fluid approximation). Any of these models can be coupled with passive scalar transport, inertial particle transport as well as Lagrangian particle transport.

Validation

It is a well-known fact that incompressible ideal hydrodynamics and magnetohydrodynamics possess several invariant quantities (see section 2.2). The exact conservation of these integral quantities provides a convenient benchmark for the nonlinear part of the code ANK.

Let us begin with incompressible ideal magnetohydrodynamics. In this case, the conserved quantities are the total energy

E = 1 2 u 2 (3.21)
and the kinetic helicity

H k = u • ω . (3.22)
In order to test the conservation of these quantities, we evolved the velocity field setting the viscosity equal to zero, using a grid with a linear resolution of N = 64 elements and single precision arithmetic. The results, shown in figure 3.1, clearly show that non-conservation of the invariants is at the level of round-off errors.

Similarly, for the case of incompressible ideal magnetohydrodynamics, we advanced the equations setting both the viscosity and the resistivity equal to zero. Again, the grid has a linear resolution of N = 64 elements and the arithmetic is done in single precision. Figure 3.2 shows the absolute value of the error as a function of time. It is again clear that the non-conservation of the invariants is at the level of round-off error. We must note here that in practice, numerical conservation of the invariants cannot hold for arbitrary time intervals. Because of the extremely low numerical dissipation of the spectral method, and the complete absence of physical dissipation, there is no mechanism to smooth-out discontinuities. Yet both incompressible ideal hydrodynamics and magnetohydrodynamics have the tendency to build very strong gradients characterized by ever-decreasing length scales. Thus, it is inevitable that at some point the characteristic scale of these gradients will become comparable to the mesh size. Then, the solution will develop spurious oscillations due to the Gibbs phenomenon and the numerical conservation of the invariants will be broken.

Resolution and convergence

Spectral methods are known to have excellent convergence properties: for linear problems, one can prove that the error drops exponentially with the number of grid points [START_REF] Gottlieb | Numerical analysis of spectral methods[END_REF]. This is to be contrasted with finite difference or finite volume methods where the error always drops algebraically (with a power equal to the order of the method). For non-linear problems, no such proof exists but spectral methods still exhibit excellent convergence properties.

In hydrodynamics and magnetohydrodynamics, the a priori requirement for convergence is that the product of the dissipative scale l d with the maximum representable wave-number k max is a given number greater than unity:

l d k max = C d > 1. (3.23)
A common choice that we adopted in all our simulations is C d = 2, which corresponds to two dissipative scales per grid interval.

One can check convergence a posteriori by examining the behavior of the energy spectra in the range of high wavenumbers. In a well-resolved simulation, there should be a range of exponential decrease in the spectrum following the power-law inertial range of scales. In a badly-resolved simulation, there is a tendency for structures smaller than the distance ∆x between two grid points to form. If such a structure forms, then from the point of view of the Fourier transform it is equivalent to a discontinuity in the field, and the Figure 3.2: Absolute value of the error in the conservation of total energy, cross-helicity and magnetic helicity in incompressible ideal magnetohydrodynamics. Gibbs phenomenon will lead to spurious oscillations which will be manifested as a pile-up of energy in the high-wavenumber range of the energy spectrum. Because the energy spectrum of a discontinuous field is in general a power-law, in under-resolved simulation one expects power-law energy spectra up to the maximum wavenumber k max .

Figure 3.3 compares energy spectra from two identical hydrodynamic simulations, one well-resolved with k max l d = 2.0 and another under-resolved with k max l d = 0.5. In the case of the under-resolved simulation, there is spurious pile-up of energy in the high-wavenumber range which leads to a much shallower, power-law spectrum.

Godunov Methods

Because of the use of a truncated Fourier series as the fundamental representation of the solution, spectral methods are well-suited for incompressible flow but perform poorly in compressible, shock-dominated flow. This poor performance is the consequence of the Gibbs phenomenon near discontinuities of the flow variables such as shocks. Although there has been attempts to extend spectral methods in compressible flows [START_REF] Passot | Numerical simulation of compressible homogeneous flows in the turbulent regime[END_REF], Vazquez-Semadeni et al., 1995, Passot et al., 1995] Godunov methods [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF] present a much more viable choice for the study of shock-dominated flows. Indeed, Godunov methods have many desirable properties such as the exact conservation of mass, linear momentum and energy. They are also relatively robust and can handle easily shocks and other discontinuities without the need of adding artificial viscosity.

In this section, we will present the basic principles of Godunov methods for solving the equations of ideal compressible magnetohydrodynamics

Conservative form of the governing equations

In the context of Godunov methods, we will focus on the equations of ideal compressible magnetohydrodynamics ignoring the effects of viscosity and resistivity. This is justified by the fact that the Godunov method focuses on the computation of the inviscid flux, whereas viscous and resistive fluxes are in our case computed by central finite differences. The most relative form of the governing equations is the conservative one:

∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρuu -BB) = -∇P, ∂ t E + ∇ • ((E + P )u -B(B • u)) = 0, ∂ t B + ∇ • (uB -Bu) = 0, ∇ • B = 0. (3.24)
In the above equations, ρ is the density, u is the velocity, B is the magnetic field, P is the total pressure defined by

P = p + 1 2 B 2 (3.25)
and E is the total energy defined by

E = e + 1 2 ρu 2 + 1 2 B 2 (3.26)
with e the internal energy of the fluid. We will assume a perfect gas equation of state:

p = (γ -1)e, γ = 5 3 . (3.27)
It is useful to write equations 3.24 in the following form:

∂ t U + ∂ x F + ∂ y G + ∂ z H = 0 (3.28)
where U = (ρ, ρu x , ρu y , ρu z , B x , B y , B z , E) T (3.29) and the flux vector F is given by 

           ρu x ρu 2 x + P -B 2 x ρu x u y -B x B y ρu x u z -B x B z 0 B y u x -B x u y B z u x -B x u z (E + P )u x -B x (B • u)             (3.30)
with corresponding expressions for the fluxes G and H being completely symmetric.

Principles of the method

Godunov methods are well suited to handle shocks and discontinuities. When using Godunov methods, there is no need to introduce artificial viscosity or resistivity as discontinuities are naturally resolved within a few computational cells. Godunov methods are an example of a finite volume method, which means that each computational cell represents an elementary fluid volume element. For a first order method, this element is characterized by a uniform density, velocity, pressure and magnetic field which represents an average value. Second and third order methods consider respectively linear and parabolic profiles for the fundamental variables within the computational cell. At the interface between neighboring cells, fluxes of matter, momentum and energy are exchanged in order to advance the fundamental quantities forward in time.

There are several ways of recovering the slope of the profiles in a cell in order to reproduce better the monotonicity of the hydrodynamic variables. Methods of slope computation that preserve the monotonicity are called TVD (total variation diminishing) methods, and are known to exhibit superior stability and convergence properties. In our numerical experiments, we use exclusively the MINMOD slope computation method [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF].

Let the center of each cell be located at the point (x i , y j , z k ). For each given cell, faces normal to the x-axis have coordinates x i±1/2 and cover a surface element defined by the inequalities y j-1/2 < y < y j+1/2 and z k-1/2 < z < z k+1/2 . The coordinates of the other faces can be defined similarly.

The fundamental variables (density, momentum, energy, magnetic field) are volume-averaged over a cell and the discrete values are defined at the cell center. For example:

ρ i,j,k = 1 ∆x∆y∆z x i+1/2 x i-1/2 y j+1/2 y j-1/2 z k+1/2 z k-1/2
ρ(x , y , z ) dx dy dz .

(3.31)

After spatial integration over a cell and time integration from t n to t n+1 , equation 3.28 becomes

U n+1 i,j,k -U n i,j,k ∆t + F n+1/2 i+1/2,j,k -F n+1/2 i-1/2,j,k ∆x + G n+1/2 i,j+1/2,k -G n+1/2 i,j-1/2,k ∆y + H n+1/2 i,j,k+1/2 -F n+1/2 i,j,k-1/2 ∆z = 0 (3.32)
where superscripts n and n + 1 refer to quantities considered at the time instants t n and t n+1 respectively. U n i,j,k and U n+1 i,j,k are the cell-averaged variables at time t n and t n+1 , while the time-and face-averaged flux is defined by i,j,k±1/2 are approximately calculated. However, even if the fluxes are not completely accurate, the method conserves mass, momentum and energy by construction, since the amount retrieved from one cell is exactly given to its neighbors. The expression (3.32) to update the variables does not contain derivatives but only flux differences and this is why discontinuities are well-resolved.

F n+1/2 i+1/2,j,k = 1 ∆t∆x∆z t n+1 t n y j+1/2 y j-1/2 z k+1/2 z k-1/2 F (x i+1/2 ,

The Riemann problem

In a Godunov method, the most important point is to calculate the flux exchange between two neighboring states U L and U R . This problem is known as the Riemann problem, and since no characteristic scale is involved in its definition, it is self-similar. This statement means that the solution at (x, t) can be deduced from the solution at (x , t ): U (x/t) = U (x /t ). An important implication of this is that the flux exchanged between the two states is constant in time. For a Godunov method, the importance of an accurate Riemann solver which is able to provide the fluxes in equation (3.32) cannot be overstated.

For a general hyperbolic system

∂ t U + ∂ x F (U ) = 0 (3.34)
we have the important result that when the Jacobian of the above system is a constant matrix, the solution of the Riemann problem can be found explicitly [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]. Let us consider the simple linear scalar advection equation

∂ t u + a∂ x u = 0. (3.35)
The solution of this equation is given by f (xat), where f is a function depending on the initial conditions. The rate of change of u along the characteristic curve defined by dx/dt = a is evidently zero (a is the characteristic speed). In this case, the solution of the Riemann problem is given by

u = u L , x -at < 0, u = u R , x -at > 0. (3.36)
Let us generalize to a linear system of m equations:

∂ t U + A∂ x U = 0 (3.37)
where A is the constant Jacobian matrix. This matrix can be put in diagonal form, A = K -1 ΛK, where Λ is a diagonal matrix whose elements are the eigenvalues of A, denoted by λ i . The matrix K = [K 1 , K 2 , ..., K m ] contains the eigenvectors of A defined by the relation AK i = λ i K i . Defining the matrix W by W = K -1 U and multiplying the advection equation by K -1 , we have the equation

∂ t W + Λ∂ x W = 0. (3.38)
Since the matrix Λ is diagonal, the components of the vector W are completely decoupled. The general solution of the system is thus W i (xλ i t), or in terms of the U variable,

U (x, t) = m i=1 W i (x -λ i t)K i . (3.39)
Let us now consider the Riemann problem defined by U (x, 0) = U L for x < 0 and U (x, 0) = U R for x > 0. The constants U L and U R can be written as

U L = m i=1 α i K i , U R = m i=1 β i K i . (3.40) Thus, W i = α i if x < 0 and W i = β i if x > 0.
For each given position x and time t there is an eigenvalue λ i such that λ i < x/t < λ i+1 . This implies that the solution of the multidimensional linear Riemann problem can be expressed in the form

U (x, t) = m i=L+1 α i K i + L i=1 β i K i .
(3.41) Physically, the above form means that the two states have been decomposed into m eigen-modes, each with associated wave speed λ i . It is further possible, for each location, to determine whether the corresponding It is now possible to estimate the flux exchanged between two neighboring cells. The value of U at the interface, U (0), can be written as:

U (0) = m i=L+1 α i K i + L i=1 β i K i , = m i=1 α i K i + L i=1 (β i -α i )K i , = m i=L+1 (α i -β i )K i + m i=1 β i K i , = 1 2 (U L + U R ) + 1 2 m i=1 sign(λ i )(β i -α i )K i (3.42)
where the last expression was obtained by summing the two previous expressions. Since the problem is linear (the Jacobian A is a constant matrix), the flux F is equal to A(U (0)):

F (U (0)) = 1 2 (A(U L ) + A(U R )) + 1 2 m i=1 sign(λ i )(β i -α i )AK i , = 1 2 (F L + F R ) + 1 2 m i=1 |λ i |(β i -α i )K i .
(3.43)

In our numerical studies, we used exclusively the Lax-Friedrichs scheme for treating the nonlinear terms [START_REF] Foias | Systems of conservation equations with a convex extension[END_REF]Lax, 1971, LeVeque, 1992]. In this scheme, the flux is calculated according to the formula

F (U ) = 1 2 (F L + F R ) + |λ max | 1 2 (U R -U L ) (3.44)
where λ max is the maximum wave propagation speed. This scheme is the most stable among Riemann solvers, but it is also the most diffusive.

Riemann solvers for magnetohydrodynamics

In the case of compressible magnetohydrodynamics, the wave fan entails 7 waves, namely: 2 fast magnetosonic waves, 2 transverse Alfvén waves, 2 slow magneto-sonic waves and one entropy wave which does not propagate and moves with the fluid. Since the fast waves propagate faster than the others, they always are at the boundaries of the wave fan. In this section, we describe various Riemann solvers which have been developed to solve the MHD equations. It is worth stressing that no exact solver is known in the MHD case except in restricted cases. Even if it would be the case, such a solver would probably, as it is the case for the exact hydrodynamic solver, require an iterative procedure which would require a large amount of computing power. In Godunov schemes, there are various methods to meet the requirement of zero magnetic field divergence. In our case, the magnetic field vector B is located on the faces of the computational cell whereas all other variables are located at the center [START_REF] Fromang | A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics[END_REF].

The code DUMSES

The code DUMSES is a simplified version of the RAMSES code [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES[END_REF][START_REF] Fromang | A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics[END_REF] where the capability of adaptive mesh refinement is removed. Thus, DUMSES is essentially a code implementing the Godunov method on a regular grid. This results in a much simpler code in terms of development, but also has as a result a significant speed improvement over RAMSES.

Validation

In the context of this thesis, we modified DUMSES, which was previously only able to treat a polytropic gas with γ > 1, rendering it capable of treating the case of an isothermal gas. In order to validate these changes, we implemented a two-dimensional tilted shock test following the test proposed by [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF], which is in turn based on the proposal of [START_REF] Brio | An upwind differencing scheme for the equations of ideal magnetohydrodynamics[END_REF].

As a first test, we simulated a one-dimensional shock with initial conditions

ρ L = 1, ρ R = 0.1, B L y = 5 √ 4π , B R y = 2 √ 4π , B L z = 0, B R z =, v L x = 0, v R x = 0, v L y = 0, v R y = 0, v L z = 0, v R z = 0. (3.45)
where the superscript L corresponds to values left (upstream) of the shock and the superscript R to values right (downstream) of the shock. The subscript x corresponds to a component perpendicular to the shock, while the subscripts y , z correspond to components parallel to the shock surface.

Our results are shown in figure 3.5, to be compared with figure 1a of [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF], which is reproduced here in figure 3.6. We observe excellent agreement between the results of DUMSES and those of [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF].

In order to validate the ability of DUMSES to properly handle multi-dimensional shocks, we implemented a tilted two-dimensional shock test proposed by [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF]. In this test, which again corresponds to figure 1a in [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF] (reproduced as figure 3.6 here), the shock is tilted 45 degrees with respect to the two-dimensional computational domain. Thus, the surface of the shock runs from the upper left corner to the lower right corner of the domain. The initial conditions are the same as in the one-dimensional test, and are given by relations (3.45).

The results of DUMSES again show very good agreement with those of [START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF], although in two dimensions there is some additional smoothing of the surfaces of discontinuity. Energy dissipation in ambipolar diffusion magnetohydrodynamics 1

Abstract

The dissipation of kinetic and magnetic energy in the interstellar medium (ISM) can proceed through viscous, Ohmic or ion-neutral friction, also referred to as ambipolar diffusion (AD). It occurs at very small scales compared to the scales at which energy is presumed to be injected. This localized heating may impact the ISM evolution but also its chemistry, thus providing observable features.

Here, we perform 3D spectral simulations of decaying magnetohydrodynamic (MHD) turbulence. We include viscous and Ohmic dissipation with or without AD and study the properties of dissipative structures at the time of the dissipation peak. We find that the AD heating power-spectrum peaks at scales in the inertial range. This is due to a strong alignment of the magnetic and current vectors in the dissipative range, in the presence of AD. AD also mildly affects the magnetic and kinetic power spectrum, with a shallower slope over the inertial and dissipative range of scales. AD affects much greater scales than the expected AD scale (defined up to a dimensionless factor as the product of the ions-neutrals coupling time scale with the root mean square velocity of the flow), which appears to play no specific role in our simulations. We conjecture this is due to a large coherence time for the AD effects.

We find that energy dissipation is concentrated on thin sheets, each sheet being predominantly either Ohmic, viscous or AD. Dissipation is highly concentrated in space: 30 percent of it takes place in less than 3 percent of the volume. Its probability density function follows a log-normal law with a power law tail at high values which hints at intermittency. We quantify the intermittency by use of anomalous structure function exponents for the velocity and magnetic fields, which we match with generalized She & Levêque models. Initial conditions (Orzsag-Tang or ABC-flow) and the inclusion of AD or not all affect the intermittency. Finally, we extract structures of high dissipation, defined as connected sets of points where the total dissipation is most intense and we measure the scaling exponents of their geometric and dynamical characteristics: the inclusion of AD favors small sizes in the dissipative range.

Introduction

In partly ionized astrophysical fluids, magnetic fields remain attached to the charged particles. The neutral fluid does not feel the Lorentz force and usually drifts with respect to the charges. Only the ion-neutral drag can remind the neutrals about the existence of magnetic fields. The fields are then able to slip through the bulk of the fluid: this process is called ambipolar diffusion (AD). [START_REF] Mestel | Star formation in magnetic dust clouds[END_REF] were the first to realize its importance in the context of star formation, where it would help magnetic fields to diffuse out of a contracting dense core and allow it to form a star. [START_REF] Mullan | The structure of transverse hydromagnetic shocks in regions oflow ionization[END_REF] discovered how ambipolar diffusion could influence the dynamics of shocks which then yielded a wealth of papers on the chemical signatures of C-type shocks, starting with [START_REF] Draine | Magnetohydrodynamic shock waves in molecular clouds[END_REF] and [START_REF] Flower | Theoretical studies of interstellar molecular shocks. I -General formulation and effects of the ion-molecule chemistry[END_REF]. [START_REF] Tóth | Simulations of the Wardle Instability of C-Type Shock Waves[END_REF] produced the first multi-dimensional simulations with AD for the stability of such shocks, and this opened the way to a collection of analytical and numerical studies in various astrophysical contexts. [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF] and [START_REF] Brandenburg | Effects of Pressure and Resistivity on the Ambipolar Diffusion Singularity: Too Little, Too Late[END_REF] envisaged that ambipolar diffusion could form very sharp structures, which would then induce the Ohmic resistivity to reconnect the magnetic field in the interstellar medium (ISM): they argued this could be a key element in solving observational problems with the galactic dynamo [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF]. In discs [START_REF] Blaes | Local shear instabilities in weakly ionized, weakly magnetized disks[END_REF], Brandenburg et al. [1995] and Mac [START_REF] Low | Incorporation of ambipolar diffusion into the ZEUS magetohydrodynamics code[END_REF] showed AD was able to modulate the magneto-rotational instability. In the context of clouds and stars formation [START_REF] Nakamura | Quiescent Cores and the Efficiency of Turbulence-accelerated, Magnetically Regulated Star Formation[END_REF], [START_REF] Nakamura | A Numerical Model of Hercules A by Magnetic Tower: Jet/Lobe Transition, Wiggling, and the Magnetic Field Distribution[END_REF], [START_REF] Kudoh | Three-dimensional Simulation of Magnetized Cloud Fragmentation Induced by Nonlinear Flows and Ambipolar Diffusion[END_REF] and [START_REF] Kudoh | Formation of Collapsing Cores in Subcritical Magnetic Clouds: Three-dimensional Magnetohydrodynamic Simulations with Ambipolar Diffusion[END_REF] computed how the magnetic support of clouds can leak out to let the gas condense and form dense cores and stars. Finally, in the context of ISM turbulence [START_REF] Padoan | Ambipolar Drift Heating in Turbulent Molecular Clouds[END_REF], [START_REF] Zweibel | Ambipolar drift in a turbulent medium[END_REF], [START_REF] Oishi | The Inability of Ambipolar Diffusion to Set a Characteristic Mass Scale in Molecular Clouds[END_REF], [START_REF] Li | Sub-Alfvénic Nonideal MHD Turbulence Simulations with Ambipolar Diffusion. I. Turbulence Statistics[END_REF], McKee et al. [2010], Li et al. [2012] and Li et al. [2012] focused on how magnetic fields decouple from the neutrals velocity or density and estimated the heating resulting from the ion-neutral drift.

In the diffuse interstellar medium (ISM), turbulent energy dissipation can be an important source of suprathermal energy driving hot chemistry [START_REF] Falgarone | The intermittency of turbulence in interstellar clouds: implications for the gas kinetic temperature and decoupling of heavy particles from the gas motions[END_REF]. This may be evidenced by the observed high values of the column density of species such as CH + and SH + . The formation of such species requires energy barriers of the order of at least 2000 K to be overcome, in clouds where the average temperature is known to be a few tens of K. One possible explanation is that these cold clouds contain pockets of hot gas heated by intermittent turbulent dissipation. Hot chemistry is activated there, and it is possible to construct models of turbulent dissipation that account for the high column densities of CH + [START_REF] Godard | Models of turbulent dissipation regions in the diffuse interstellar medium[END_REF].

In a turbulent magnetized fluid, dissipation can of course be due to viscosity or resistivity, but when the fluid is partially ionized and AD is at play, there can also be a significant contribution from the heat released by ion-neutral friction, as demonstrated by several authors in the context of the ISM [START_REF] Scalo | Heating of dense interstellar clouds by magnetic ion slip -A constraint on cloud field strengths[END_REF][START_REF] Zweibel | Hydromagnetic wave dissipation in molecular clouds[END_REF][START_REF] Elmegreen | Energy dissipation in clumpy magnetic clouds[END_REF][START_REF] Padoan | Ambipolar Drift Heating in Turbulent Molecular Clouds[END_REF], Li et al., 2012]. Not only does the heating help to raise the temperature which increases the rate of some chemical reactions, but the ion-neutral drift velocity provides additional energy in the reaction frame for ion-neutral reactions. In some instances, this can open new chemical routes which would otherwise be blocked by reaction barriers. The places of strong AD heating are thus expected to bear specific chemical signatures such as the ones encountered in magnetized vortices [START_REF] Godard | Models of turbulent dissipation regions in the diffuse interstellar medium[END_REF] or C-shocks [START_REF] Lesaffre | Low-velocity shocks: signatures of turbulent dissipation in diffuse irradiated gas[END_REF]. We would hence like to characterize the geometry and statistical properties of the regions of strong turbulent dissipation in the ISM.

Before us, [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF] conducted a thorough study of the statistics of strong dissipation in the context of incompressible pure magnetohydrodynamic (MHD) turbulence. In this paper we make some progress towards the physics of the ISM and we work with incompressible MHD turbulence with or without AD. We stay within the model of incompressible MHD in a first step to link our work with [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF] and to allow the use of spectral methods which are well suited for the study of small-scale dissipative structures because of their very low level of numerical dissipation. In section 2 and 3 we briefly describe the equations and the numerical method used and we present the simulations that were performed. In section 4 we present an overall picture of the dissipation fields through the time evolution of their average values, their pdfs and their spectra. We also provide a qualitative view of the dissipation field in physical space through color maps. Section 5 deals with the extreme dissipative events, it begins with a discussion of the structure functions of the velocity and magnetic fields and concludes with results from the statistical analysis of the geometrical and dynamical properties of structures of high dissipation. We discuss and conclude our results in section 6.

The equations 4.3.1 Ambipolar drift

In a partly ionized fluid, the time-dependent evolution of both the neutral and the ionized fluids should in principle be followed. However, in circumstances that we will make explicit below (see subsection 4.3.4), we can adopt the strong coupling approximation. In this approximation, we neglect the inertia, pressure and viscosity of the ions in the ion momentum evolution and we are left with the balance between the ion-neutral drag and the Lorentz force

γρ i ρ n (U i -U n ) = J × B (4.1)
where the current density

J = 1 4π (∇ × B),
ρ i and ρ n are the ion and neutral mass density (ρ i ρ n ), U i and U n are the ions and neutrals respective velocities, and where γ is the coefficient of ion-neutral drag

γ = σv in m i + µ
with m i and µ the ions and neutrals mass per particle and σv in the ion-neutrals collision rate. Assuming that n He = 0.2n H2 , we find µ = 2.33m p for molecular gas, where m p is the mass of the proton. In diffuse clouds, the average mass per ion is m i = 12m p as the dominant ion is C + . Following [START_REF] Draine | Magnetohydrodynamic shock waves in molecular clouds[END_REF], we take σv in = 1.9 × 10 -9 cm 3 s -1 and we arrive at γ = 6.7 × 10 13 cm 3 s -1 g -1 . Within these approximations, the above balance (4.1) expresses the ion-neutral drift velocity as a function of the magnetic field. When this is plugged into the induction equation

∂ t B = ∇ × (U i × B) + η∇ 2 B
one recovers the ambipolar diffusion term: [START_REF] Balbus | Linear Analysis of the Hall Effect in Protostellar Disks[END_REF], which can be developed into

∇ × 1 γρ i ρ n (J × B) × B [cf.
∇ × - B 2 γρ i ρ n J + J.B γρ i ρ n B
from which properly speaking only the first term takes the form of a diffusion, with diffusion coefficient λ AD = B 2 γρiρn [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF], but the qualificative is usually retained for the whole term. In particular, [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF] recognized that the second term steepens the magnetic field profile near magnetic nulls. It should finally be noted that AD itself is not per se able to reconnect the field lines: this requires Ohmic diffusion.

Incompressible MHD

We now take u 0 a typical velocity and l 0 a typical length scale as unit velocity and unit length to normalize our equations. We also define t 0 = l 0 /u 0 as the unit of time. We write the non-dimensional velocity u = U cdm /u 0 where U cdm is the center of mass velocity

U cdm = ρ i U i + ρ n U n ρ
with ρ ∼ ρ n the total density of the gas, ρ = ρ i + ρ n . We write the non-dimensional Alfvén velocity b = B/ √ 4πρ/u 0 . The non-dimensional current is simply j = ∇×b where ∇ is now understood as derivatives in coordinates in units of l 0 : ∇ → l 0 ∇.

In the diffuse ISM, the sonic Mach number M s = u 0 /c s (where c s is the speed of sound) as well as the Alfvén Mach number M a = 1/|b| take values in the range 10 -1 -10 [ [START_REF] Elmegreen | Interstellar turbulence i: Observations and processes[END_REF]. This wide range of values suggests that although most of ISM turbulence is highly compressible, incompressible turbulence is not irrelevant since in a weakly compressible flow the density fluctuations ∆ρ/ρ ∼ M 2 s are of the order of the square of the sonic Mach number. Hence a turbulent flow with M s < 0.3 can be adequately described by the incompressible equations. For example, the studies of [START_REF] Brandenburg | Effects of Pressure and Resistivity on the Ambipolar Diffusion Singularity: Too Little, Too Late[END_REF], [START_REF] Zweibel | Ambipolar drift in a turbulent medium[END_REF] or [START_REF] Godard | Models of turbulent dissipation regions in the diffuse interstellar medium[END_REF] on turbulent dissipation with AD were all based on the incompressible equations of motion.

We use the equations in [START_REF] Balbus | Linear Analysis of the Hall Effect in Protostellar Disks[END_REF] and the above notations to derive the equations of incompressible, viscous, resistive, AD MHD :

∂ t u + (u • ∇) u = -∇ p + j × b + Re -1 ∇ 2 u ∂ t b = ∇ × (u × b) + Re -1 a ∇ × (j × b) × b + Re -1 m ∇ 2 b (4.2)
where u and b satisfy ∇ • u = 0 and ∇ • b = 0 and p = P/(u 2 0 ρ) is the non-dimensional pressure, with P the actual thermal pressure. Equations (4.2) are parameterized by three non-dimensional numbers Re, Re a and Re m , for which we now give estimates.

Reynolds numbers

Firstly, the kinetic Reynolds number Re = u 0 l 0 /ν, where ν is the molecular viscosity of the fluid, expresses the relative importance of inertial terms in comparison to the viscous term. In the neutral ISM, assuming that the most significant contribution to viscosity is given by H 2 collisions, we have ν ∼ 1 3 λ H2 c s where λ H2 is the mean free path of the H 2 molecule and c s = (Γk B T /µ) 1/2 is the isentropic sound speed, with Γ 5/3 the ratio of specific heats, and k B the Boltzmann constant. The mean free path is given by λ H2 ∼ (n H2 σ H2 ) -1 where n H2 is the number density of H 2 and σ H2 = 3 × 10 -15 cm 2 is an estimate of the cross section of H 2 collisions [START_REF] Monchick | Theoretical studies of H2-H2 collisions. II. Scattering and transport cross sections of hydrogen at low energies: Tests of a new ab initio vibrotor potential[END_REF]. For molecular gas n H2 = 0.5n H where n H is the hydrogen nuclei density. Under these assumptions, the kinetic Reynolds number is of the order

Re ∼ 1.8 × 10 7 n H 100 cm -3 u 0 1 km s -1 l 0 10 pc T 100 K -1 2 .
Elmegreen and [START_REF] Elmegreen | Interstellar turbulence i: Observations and processes[END_REF] quote typical values of the kinetic Reynolds number in the cold ISM ranging from 10 5 to 10 7 .

Secondly, Re m = u 0 l 0 /η is the magnetic Reynolds number, where η is the resistivity. The magnetic Reynolds number expresses the relative importance of advection in comparison to Ohmic diffusion in the dynamics of the magnetic field. In a system with a large value of the magnetic Reynolds number, the dynamics of the magnetic field is dominated by advection and stretching. The value of the resistivity is given by

η = 234 n n e T 1/2 cm 2 s -1
[ [START_REF] Balbus | Linear Analysis of the Hall Effect in Protostellar Disks[END_REF], where n is the total number density and n e is the electron density. The order of magnitude of the magnetic Reynolds number is

Re m = 2.2 × 10 17 l 0 10 pc u 0 1 km s -1 × × n e 10 -4 n H T 100 K -1/2
where we assumed n = 0.6n H for molecular gas.

Lastly, the AD Reynolds number Re a helps to measure the ratio of the ambipolar to advective electromotive forces in the induction equation:

Re a = t 0 t a , t a = 1 γρ i (4.3)
where t a can be recognized in equation (4.1) as the ion-neutral friction time scale. The quantities Re a and l a should not be confused with their more usual definitions R AD and AD as introduced by [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF], for example. For instance, the usual values depend on the r.m.s. velocity and magnetic field, whereas our definition encompasses only the ion-drift time: see the subsection 4.3.4 for more details.

We can also write it as a ratio of length scales

Re a = l 0 l a
where we define l a = t a u 0 (4.4) which gives a typical length scale for ion-neutral decoupling. The AD Reynolds number is an increasing function of the ionization fraction x = ρ i /ρ. Using C + as the dominant ion of the ISM, we find Re a = 4.9 × 10 3 l 0 10 pc

n C + 10 -4 n H n H 100 cm -3 u 0 1 km s -1 -1 .
In the ISM, we have Re a Re Re m which suggests the ordering l a l ν l η for the ambipolar, viscous and resistive dissipation scales. However, our finite computing power does not allow much dynamics of scales and here we can only afford l a > l D where l D = l ν ∼ l η is a single dissipative scale. In this study the magnetic Prandtl number Pr m = ν/η is therefore taken equal to unity, so that the hydrodynamic and magnetic Reynolds numbers are equal. This choice results in a single dissipative range of scales for both the velocity and the magnetic fields, a fact which simplifies the analysis considerably.

Although this numerical study is confined to Reynolds numbers that are several orders of magnitude lower than those found in the interstellar medium, it is relevant to ISM physics in the sense that it allows a detailed quantitative study of the dissipation field as well as the relative importance of the different types of turbulent dissipation, which can all be important as a heating source for ISM chemistry.

Lengths scales associated with AD

We outline here various scales introduced by AD physics. The case of C-shocks allows to clearly separate them. We have already introduced the length scale l a = u 0 t a which corresponds to the re-coupling length scale between ions and neutrals. It is the scale of variation of the neutral's velocity in a C-shock with entrance velocity u 0 , which therefore has a length on the order of l a [START_REF] Flower | Non-thermal sputtering of interstellar grains in magnetohydrodynamic shocks[END_REF].

Under the strong coupling approximation, the scale of variation of the ions velocity in the same C-shock is l ai = u 0 t a /M2 a where M a is the transverse Alfvénic Mach number of the shock. For shocks with a large Alfvénic Mach number, this length-scale is significantly smaller than l a and the structure of the shock consists of a front in the ions velocity followed by a smoother transition for the neutral velocity [see Fig. 1 of Li et al., 2006, for example].

In the case of typical ISM turbulence, though, M a is of order one, and both length scales do not differ significantly. In fact, [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF] constructed the AD diffusion Reynolds number of eddies of length scale and velocity U based on the AD diffusion coefficient: [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF] then argue that only eddies of length scales below AD should be affected by AD. We prefer to get a similar estimate by comparing the Fourier amplitudes of the AD e.m.f. (Re -1 a (j × b) × b → Re -1 a kb 3 ) and the inertial e.m.f. (u × b → ub) in the the induction equation (4.2). Wave numbers above the critical wave number

R AD ( ) = U/λ AD = / AD where AD = U t a /M 2 a with M a = U/c a . 2 Zweibel
k a = Re a u 2 / b 2 (4.5)
should be AD dominated. We hence define a = 2π/k a accordingly, as the length scale below which AD should be effective. Note that AD and l 0 a differ only by a factor of 2π.

Similarly, we can estimate the length scale below which the strong coupling approximation breaks down by comparing the magnitude of the neglected inertial term D t ρ i U i to the coupling term ρ n (U n -U i )/t a .

Assuming that U n , U i and U n -U i all share the typical magnitude u 0 , we then get a critical wavenumber k two-fluids ρ n /ρ i /l a above which the strong coupling approximation fails and the two-fluids approximation is needed. Provided ρ i /ρ n is small (it is typically lower than 10 -3 if the main charges are C + ions), the strong coupling approximation breaks down at scales much smaller than the typical AD diffusion scale. Other authors [Oishi andMac Low, 2006b, Padoan et al., 2000] have claimed that the strong coupling approximation breaks down as soon as < AD or R AD < 1, where the ions inertia does not appear explicitly. But Fig. 1 of Li et al. [2006] shows a C-shock computed with the two-fluid approximation (solid) compared to an analytical solution (dashed) using the strong coupling approximation, and the agreement is perfect. We hence believe that the strong coupling approximation is a very good one in the low ionized ISM where ρ i /ρ n 1, even in cases where R AD > 1. In particular, k two-fluids is at least a few ten times above the largest wave number in all our AD simulations, which amply justifies our use of the strong coupling approximation.

Finally, the observed emission of the ISM depends on the chemical and thermal state of the gas, which are strongly linked to the heating. The scale at which the heating takes place may not necessarily be directly connected to the scale a where AD undergoes a change of dynamical regime. In fact we will see that it is not the case in the present paper, and we will be forced to introduce yet an other length scale * a for the typical thickness of sheets of strong AD heating.

The simulations

Method

Many compressible methods have been devised to treat AD in the strong coupling or in the two-fluids approximations: see [START_REF] Masson | Incorporating Ambipolar and Ohmic Diffusion in the AMR MHD Code RAMSES[END_REF] and references therein. Here, we solve the strong-coupling incompressible equations (4.2) in 3D using a spectral method for various values of the parameters and various initial conditions. Our spectral code ANK 3 is fully de-aliased by use of the phase-shift method of Patterson and Orszag [1971a] and uses polyhedral truncation [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF]. Polyhedral truncation considers only these wave-vectors for which the sum of any two of their components does not exceed 2N/3. Similarly with the widely used two-thirds rule, but contrary to an isotropic spherical truncation, this truncation scheme does not possess a sharp limit in wavenumber space. Polyhedral truncation allows us to keep 55 percent of the modes active, in comparison to 33 percent for the standard two-thirds rule, resulting in a more accurate description of the small scales. The Fourier transforms are computed with FFTW with single precision accuracy, and a standard fourth-order Runge-Kutta method is used for time integration. We checked that our code gives the correct solution for Alfvén waves damped by AD. As a resolution check for all simulations, we also looked for bumps in the kinetic and magnetic energy spectra near the truncation limit. In the case of the kinetic energy spectra we found bumps no larger than 15 %, whereas no bumps were present in the magnetic energy spectra. Note that we do not include a driving force in equations (4.2): our simulations are freely decaying. The MHD simulations with 512 3 resolution take about 5000 CPU hours until the peak of dissipation but the equivalent AD-MHD simulations require ten times more CPU time due to the more stringent time-step requirement.

In the spectral simulations performed, the velocity and magnetic fields are defined on a regular Cartesian grid of points, while boundary conditions are periodic in all directions. Note the total length of the computational domain is 2π and the smallest non-zero wave-vector has a norm of one.

Initial conditions

We use two types of initial conditions, corresponding to two different situations for the magnetic and crosshelicities.

In the first case, the three lowest non-zero wave numbers of both the velocity and the magnetic field are initially loaded with a superposition of different Arnol'd-Beltrami-Childress [ABC, see [START_REF] Dombre | Chaotic streamlines in the ABC flows[END_REF] flows (u x , u y , u z ) = (A sin(kz) + C cos(ky), B sin(kx)

+ A cos(kz), C sin(ky) + B cos(kx)). (4.6)
Different values of the coefficients A, B, C are chosen for the first three wave numbers from a uniform random number generator. In higher wave numbers a random field with energy spectrum

E(k) = C E k -3 exp -2 (k/k c ) 2 , k c = 3 (4.7)
is superposed. The phases are chosen from a uniform random number generator with the same seed for all simulations.

In the second case, the large scale initial flow is the Orszag-Tang (OT) vortex and in higher wave numbers a random velocity field with the same properties as above is superposed. In order to keep the initial value of magnetic helicity close to zero, no random magnetic field is added to the OT initial condition, in contrast to the ABC initial condition.

The compressive components of the initial velocity and magnetic fields are subtracted so that the initial condition is purely solenoidal. In all cases, the constant C E in equation (4.7) is chosen such that u 2 = b 2 = 1, so that we start from equipartition between kinetic and magnetic energy. The energy of the initial condition fields is concentrated on large scales k < k c due to the exponential cutoff in (4.7).

In the case of the ABC initial condition, the non-dimensional cross-helicity

H c = 2 u • b u 2 b 2 is ∼ 2 × 10 -3
, corresponding to a low initial correlation between the velocity field and the magnetic field.

The mean magnetic helicity

H m = a • b
where a is the vector potential with b = ∇ × a, is considerable, ∼ 0.2. In the case of the OT initial condition the non-dimensional cross-helicity is ∼ 0.1 while the mean magnetic helicity is almost zero, ∼ 1 × 10 -9 . Thus these two different initial conditions represent evolution under different constraints: in the ABC case, low cross-helicity and sizable magnetic helicity whereas in the OT case sizable cross-helicity but low magnetic helicity. This fact is important because in the ideal MHD limit (inviscid and non-resistive) the energy, cross-helicity and magnetic helicity are all conserved during the evolution. If AD is included in the ideal MHD equations, the conservation of magnetic helicity remains while energy and cross-helicity conservation are broken. This is a consequence of the form of the AD term in the induction equation, which takes the form of an advection term ∇ × (u d × b) , with u d = Re -1 a (j × b) the non-dimensional ion-neutral drift velocity. This form also implies that although AD is a dissipative process, it conserves magnetic flux and is thus unable to reconnect field lines.

Parameters

The parameters of the simulations performed are shown in table 4.1. Throughout this paper, we focus mainly on the analysis of the OT initial condition with Re a = 10, 100, ∞, and we discuss the differences with respect to the ABC initial conditions only when they arise. The Taylor microscale Reynolds number Re λ is defined as Re λ = U λRe where U = u 2 is the r.m.s. velocity and

λ = 2π ∞ 0 e(k) dk ∞ 0 k 2 e(k) dk 1 2
is the Taylor microscale, with ∞ 0 e(k) dk = 1 2 u 2 + b 2 the total energy and e(k) the total energy spectrum. The value of Re λ is given at the peak of viscous plus Ohmic dissipation where

ε = ε o + ε v ε o = Re -1 m j 2 ε v = Re -1 2 3 i,j=1 (∂ i u j + ∂ j u i ) 2 .
The time when the peak of dissipation occurs is appropriate for analysis because for a given integral length scale

L = 2π ∞ 0 k -1 e(k) dk ∞ 0 e(k) dk (4.9)
and dissipative scale (assuming Kolmogorov scaling)

l d = Re -3 ε 1 4 (4.10)
the scale separation L/l d between the energy-containing scales and the dissipative scales is maximum. Another desirable property at the peak of dissipation is quasi-stationarity, due to the time derivative of the dissipation rate which cancels at the peak, by definition. In all the MHD simulations (with Re -1 a = 0), two snapshots of the fields were recorded for analysis: one at the peak of dissipation and a second one one eddy turnover time later. We define the macroscopic eddy turnover time as

T = √ 3 L U (4.11)
where the one-directional r.m.s. velocity U/ √ 3 and the integral length scale L were both computed at dissipation peak to estimate when to output the next snapshot. For the AD simulations (with Re -1 a > 0), we could not afford to compute beyond the dissipation peak.

Power-spectra

The kinetic and magnetic energy spectra of the high resolution OT runs 10-12-14 are shown in Figure 4.1, at the temporal peak of dissipation. The spectra are compensated by the Kolmogorov law k -5/3 and normalized by U 2 . The extent of the inertial range, as defined by the portion of the spectra that has slope -5/3 is very limited, especially in the cases with AD. In the same Figure we show the limits of the inertial and dissipation ranges assumed for the analysis of section 4.6.2: they are taken from [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF] as [0.21, 1.3] and [0.025, 0.18] respectively (in units of l 0 ).

The kinetic and the magnetic energy appear to remain in approximate equipartition across all scales except for the smallest scales in the AD runs where magnetic energy dominates the kinetic energy (the tick-marks on the vertical lines can guide the eye to estimate the relative position of the curves between the upper and the bottom panel).

The vertical dash-dotted lines with square symbols correspond to the wave-number k a defined in equation (4.5), for the cases of runs 12 and 14. Surprisingly, departures from MHD spectra start at about the same wave number for all AD runs (in the range k ∼ 5 -8 for both the kinetic energy spectrum and the magnetic one): this hints at the fact that a = 2π/k a is not the proper scale to assess the dynamical importance of AD in our simulations. In particular, dynamics can be affected at scales much larger than that in run 12. Although the difference between the pure MHD and AD MHD spectra is modest, there is a clear tendency for AD to flatten the energy spectra, especially for the magnetic energy. This is in line with the idea by [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF] that AD diffuses magnetic fields on the one hand, but on the other hand helps to build sharper magnetic structures in specific places.

The dissipation field 4.5.1 Total dissipation

We present in Figure 4.2 the time evolution of volume integrated dissipation rates. The viscous and Ohmic dissipation rates follow each other closely and we don't separate their respective contributions in this figure. Most of the time, the total dissipation rates due to viscosity, resistivity and AD are of comparable magnitude. But the total AD dissipation rate is seen to peak before the Ohmic plus viscous dissipation rate, especially at low values of Re a . This makes our choice of the temporal peak of Ohmic plus viscous dissipation more appropriate to avoid the initial transient spike of AD dissipation (this spike is even more pronounced in the ABC case run 13).

We note that Ohmic dissipation is not enhanced by the presence of AD. On the contrary, the peak value of the Ohmic plus viscous dissipation decreases as Re a is decreased. Even though [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF]'s idea that AD sharpens magnetic structures at small scales is valid in our simulations, AD also smooths the fields at intermediate scales and the net effect on the global Ohmic heating is to decrease it. However, this may be due to the finite dynamical range in our simulations. Higher Re m Reynolds number simulations, if they yield enhanced magnetic power in a more extended range at small scales, could result in a globally enhanced rate of reconnection, in agreement with [START_REF] Zweibel | Current sheet formation in the interstellar medium[END_REF]'s model.

Probability distribution function

The probability density function (pdf) of the total dissipation rate 4.3 for the high-resolution run 12. The core of the pdf is very close to the log-normal distribution

ε t = ε o + ε v + ε a where ε a = Re -1 a (j × b) 2 is shown in Figure
P c (ε t ) ∝ exp - (ln ε t -µ l ) 2 σ 2 l
with mean µ l -4.27 and standard deviation σ l 1.03, while the tail of the distribution can be fitted by a power-law

P t (ε t ) ∝ ε -τ t 76 10 -3 10 -2 10 -1 10 0 10 1 ε t 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 P(ε t ) µ µ + σ µ + 2σ µ + 3σ
Total dissipation PDF Core lognormal fit Tail power-law fit The core follows a log-normal distribution with mean µ l -4.27 and standard deviation σ l 1.03 while the tail follows a power-law with exponent -τ -2.61. Vertical lines show the mean value (red) and thresholds located at 1 (black), 2 (green) and 3 (blue) standard deviations above the mean value.

with exponent τ 2.61. This power-law is one of the signatures of intermittency of dissipation [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF]. For still higher values of the total dissipation the pdf has an exponential cut-off, although these high dissipation values are close to the sampling limit. In the analysis of the next section (extraction of structures of high dissipation) effectively only the power-law range of the distribution is sampled.

In Figure 4.4 we present the cumulative probability density function of the total dissipation for run 12. It flattens out before reaching unity, which shows that high values of the dissipation are concentrated in a small volume subset of the spatial domain. The analysis of the next section concerns events that take place in this high-dissipation plateau.

Power-spectrum

We now investigate the distribution of the energy dissipation with respect to spatial scale. This can be discerned thanks to the power spectra of the dissipation fields, displayed in figure 4.5. The AD heating peaks at larger scales in comparison to Ohmic and viscous dissipation. The scale of this peak * a is actually only a few times smaller than the integral length scale. At this range of scales AD dissipation is much more important than Ohmic and viscous dissipation. This suggests that the heating due to AD has a characteristic length scale * a , which can be much larger than the dimensional estimate a = 2π/k a (4.5) for run 11. This length scale is relevant to the heating, and hence may manifest itself in structures revealed by chemical tracers of the warm chemistry of the ISM.

The AD dissipation term, Re -1 a (j × b) 2 , is proportional to the square of the Lorentz force. It is hence interesting to look at the influence of AD on the characteristics of the Lorentz force, in comparison with the pure MHD case. Figure 4.6 shows the spectrum of the Lorentz force for OT runs 10, 12 and 14. The inclusion of AD has a significant effect on the total power of the Lorentz force, reducing it importantly especially in the dissipative range. By contrast, as seen on figure 4.1, AD results in a deficit in the magnetic energy spectrum (and hence the spectrum of the current vector) which is much smaller in comparison to the deficit of the cross-product of these two vectors (the Lorentz force), and is only present on intermediate scales.

This is explained if AD has the effect of aligning the vectors j and b [as was also found in the simulations of Brandenburg et al., 1995], with a stronger effect at small scales. To put it in an other way, AD leads the magnetic field at small scales closer to a Lorentz force free configuration, where the feedback of the magnetic field evolution on the velocity field dynamics is weaker than for MHD. This was also found in the simulations by [START_REF] Brandenburg | Effects of Pressure and Resistivity on the Ambipolar Diffusion Singularity: Too Little, Too Late[END_REF]. Figure 4.5: Compensated dissipation spectra for run 12 (AD-OT, Re a = 100). Blue solid line: AD dissipation, red dashed line: Ohmic dissipation, green dotted line: viscous dissipation. We plot ke (k) in a log-lin plot, so that the area under the curve over any interval shows directly the amount of power inside this interval. We mark the position of the maximum value of ke a , at k = k * a = 2π/ * a , and the position of k = k a .

Figure 4.6: Power spectra of j × b for high resolution runs 10,12 and 14. The field becomes force-free at small scales when the strength of the AD is increased.

Here we attempt to trace this tendency back to the evolution equations. We write j ⊥ the component of the current vector j: the double cross-product (j × b) × b can then be simply written b 2 j ⊥ . This allows to write

(∂ t j) withAD = (∂ t j) withoutAD + ∇ × [∇ × (Re -1 a b 2 j ⊥ )] (4.12)
which shows that in regions where b 2 is smooth enough, the effect of AD is to diffuse out the component of the current perpendicular to b, and it does so faster at small scales like any diffusion process. Hence AD brings the field closer to a force-free state, more efficiently at small scales, except perhaps at the smallest scales where b 2 varies and the behavior of equation (4.12) is less easy to predict.

Spatial structure

Next, we consider qualitatively the different contributions to the bulk of the dissipation field. For this purpose, each different mechanism of dissipation (Ohmic, viscous and ambipolar diffusion) is assigned to a color channel: Ohmic dissipation is assigned to the red channel, viscous dissipation to the green channel and AD dissipation to the blue channel. To emphasize the structures in the bulk of the dissipation we first compute the total dissipation value ε l below which 10% of the heating occurs and the value ε u below which 90% of the dissipation occurs. We discard the pixels with total dissipation ε t < ε l , we saturate the intensity of the pixels with ε t > ε u (while keeping their intrinsic color) and we apply a logarithmic scaling for the intensity in between these two thresholds. The color of each channel is hence given by the ratio of each type of dissipation to the total dissipation

Red = ε o ε t I Green = ε v ε t I Blue = ε a ε t I
with the intensity I given by

I =    0 if ε t < ε l log(εt)-log(ε l ) log(εu)-log(ε l ) if ε l ≤ ε t ≤ ε u 1 if ε u < ε t
The color maps of a slice through the dissipation fields are shown in figures 4.7-4.9 for all high-resolution runs with the ABC initial conditions. In the pure MHD case (Figure 4.7), viscous and Ohmic dissipation are in general concentrated on thin sheets: a slice by slice inspection of the full cube reveals continuously evolving filaments at the intersection between the sheets and the plane of the slice. The length of the sheets is comparable to the integral length scale (4.9) while their thickness is comparable to the dissipation scale (4.10).

The case of AD MHD (figures 4.8,4.9) is similar in the sense that viscous and Ohmic dissipation are again concentrated on thin current sheets, although these sheets are fewer in number (more connected), and the voids of low dissipation between them are smaller. AD dissipation is significantly more diffuse than both Ohmic and viscous dissipation, and is concentrated on much thicker structures. The thickness of the AD dissipation structures seems to coincide with the AD heating length * a measured on the power spectra. In some cases the structures of strong AD dissipation are seen to surround the sheets of Ohmic or viscous dissipation: AD sandwiches Ohmic dissipation, much like in the [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF] picture. This is probably also seen on the left panel of Figure 7 of [START_REF] Padoan | Ambipolar Drift Heating in Turbulent Molecular Clouds[END_REF] where the structures of AD dissipation often go in pairs, except Ohmic dissipation is absent from their simulations, so reconnection proceeds only through numerical truncation errors. Between figure 4.7 and Figure 4.8, there is little difference in the spatial structure of the dissipation field, with similar size, shape and position for each sheets of dissipation.

The above color maps reveal only a small fraction of the pixels have mixed colors (such as cyan, yellow or purple). This suggests that the dissipative structures of different nature (Ohmic, viscous or AD) are well separated.

4.6 Intermittency and structures of high dissipation 4.6.1 Structure functions Figures 4.7-4.9 suggest that the dissipation field is not smoothly distributed in space, but is characterized by a high degree of intermittency with regions of extreme dissipation values alternating with relatively quiescent regions. The intermittent distribution of the energy dissipation rate is expected to have an impact on the structure functions of the velocity and the magnetic field. A longitudinal velocity structure function of order p is the p-th moment of the longitudinal increment of the velocity field

S u p (r) = (δu (r)) p , δu (r) = (u(x + r) -u(x)) • r
where r is the unit vector in the direction of r. The definition for the structure function of the magnetic field S b p (r) is completely analogous. According to standard Kolmogorov [1941b] phenomenology (hereafter K41), in the inertial range the structure functions exhibit power-law scaling

S u p (r) ∝ r ζ u p , l d r L
where the exponents ζ u p vary linearly with p, ζ u p = p/3. Intermittency considerations [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF] lead to deviations from the linear K41 prediction. A particularly successful model of intermittency that was introduced for hydrodynamic turbulence by She and Leveque [1994a] and generalized for MHD turbulence by [START_REF] Politano | Model of intermittency in magnetohydrodynamic turbulence[END_REF] predicts

ζ GSL p (g, C) = p g 1 - 2 g + C 1 -1 - 2 gC p g
where g is the inverse of the inertial range scaling exponent of the velocity increment

δu (r) ∝ r 1/g
and C is the co-dimension of the dissipative structures, C = 2 for filaments and C = 1 for sheets in three space dimensions. K41 phenomenology predicts g = 3, while Iroshnikov-Kraichnan (IK) MHD phenomenology [START_REF] Iroshnikov | Turbulence of a conducting fluid in a strong magnetic field[END_REF][START_REF] Kraichnan | Inertial-range spectrum of hydromagnetic turbulence[END_REF] predicts g = 4.

In order to estimate the structure function exponents from numerical data, we used the extended selfsimilarity (ESS) method introduced by [START_REF] Benzi | Extended self-similarity in turbulent flows[END_REF]. The exponents were estimated by computing the logarithmic slope

ζ u p = d log S u p (r) d log S * (r)
where S * (r) is a structure function whose scaling behavior is known from theory. In the case of MHD turbulence, S * (r) is given by the two functions

S ± z (r) = δz ∓ (δz ± ) 2
where δz ± is the increment of the Elsässer fields z ± = u±b. Starting from the MHD equations and assuming statistical homogeneity, isotropy and stationarity, one can derive analytically

S ± z (r) = - 4 3 ε ± r l d r L (4.13)
where ε ± are the dissipation rates of (z ± ) 2 [START_REF] Politano | Dynamical length scales for turbulent magnetized flows[END_REF]]. The linear scaling of S ± z (r) in the inertial range is confirmed approximately by the flattening seen in Figure 4.10 which shows compensated plots. Although the derivation of the law (4.13) is not proven in the case of AD MHD, the linear scaling of S ± z (r) with r is not further from linearity in comparison to pure MHD. The velocity and magnetic field structure function exponents up to order 8, calculated for high resolution runs 9-14 using extended self-similarity are shown in figures 4.11-4.14. We see a departure from the linear prediction of K41, a sign of intermittency, especially in the case of the magnetic field. The structure function exponents follow closely, but not exactly, the predictions of the generalized She & Levêque model with a scaling parameter g = 3 and a co-dimension for the structures of high dissipation around C = 1 (with exceptions at C = 2 and some below C = 1). This fact suggests that in the inertial range δu(r) and δb(r) are proportional to r 1/3 , which is the prediction of K41 phenomenology, rather than r 1/4 as predicted by IK phenomenology. The value C = 1 of the co-dimension suggests that the structures of high dissipation in the inertial range are sheet-like, in accordance with figures 4.7-4.9 and the analysis of section 4.6.2. The velocity exponents for runs 9, 10 and 12 hint towards C = 2 (filaments) and the magnetic exponents for the OT runs exhibit a greater degree of intermittency. This indicates that even though the exponents appear to follow generalized She & Levêque models, the model's phenomenology probably does not subtend the dissipation in our simulations.

The effect of AD on intermittency is not easy to discern from these scaling exponents. In the case of the ABC initial condition, the deviation from the K41 values is larger for both the velocity and the magnetic field. The situation is reversed in the case of the OT initial condition, where the pure MHD fields appear to be more intermittent.

Extraction of structures

Following [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF] (hereafter UR10) and Moisy and Jiménez [2004], we implemented an algorithm for the extraction of structures of high dissipation. In the problem considered, there are three types of local dissipation rates: the viscous dissipation rate,

ε v = Re -1 2 3 i,j=1 (∂ i u j + ∂ j u i ) 2 (4.14)
the Ohmic dissipation rate

ε o = Re -1 m j 2 = Re -1 m (∇ × b) 2 (4.15)
and the AD dissipation rate

ε a = Re -1 a (j × b) 2 . (4.16)
A structure of high dissipation is defined as a connected set of points x where ε(x) > ε + jσ ε , j =1,2 or 3. (4.17)

In the above relation, ε can be any of the three dissipation rates or the total dissipation rate

ε t = ε v + ε o + ε a (or ε t = ε v + ε o in the case AD is absent)
. ε is the spatial mean value of the dissipation rate and σ ε its standard deviation. For example, the three thresholds we use on total dissipation for run 12 are displayed on figures 4.3 and 4.4 over the PDF and the CDF of the total dissipation. The algorithm is capable of isolating the structures of high dissipation so that a statistical analysis of their geometric and dynamical characteristics can be performed. The extracted structures are generally sheet-like, as can be seen in figure 4.15, where all structures extracted from run 12 (AD -OT) whose characteristic linear sizes L i (see below) lie in the inertial range are shown. An example of a more complex structure can be seen in Figure 4.16. This structure is one of the largest extracted from this dataset. It is sheet-like, with its length larger than the integral length scale (4.9) and its thickness is comparable to the dissipation scale (4.10). It is overall characterized by a high degree of geometrical complexity.

After the extraction of the structures, the following quantities were computed for statistical analysis

L i = δ max m,l∈Λi |r m -r l | (4.18a) L {x,y,z}i = δ max m,l∈Λi |{x, y, z} m -{x, y, z} l | (4.18b) R i = L 2 xi + L 2 yi + L 2 zi (4.18c) V i = δ 3 m∈Λi 1 (4.18d) A i = δ 2 m∈Λi M (m) ⊂Λi 1 (4.18e) P i = δ 3 m∈Λi ε(r m ) (4.18f)
In the above definitions, Λ i is the i-th structure, m its m-th point, M (k) the set of 26 neighbors of the point m and δ = 2π/N is the grid spacing. L i is the characteristic linear scale of the structure, L {x,y,z}i is the linear scale of its projection on the three axes, R i is the characteristic linear scale of the smallest volume embedding the whole structure, V i is its volume, A i its surface and P i the volume-integrated dissipation rate.

In the next section, all statistical quantities are computed on the sample of different extracted structures. Note that the definition of R i makes it dependent of the orientation of the structure (its value can change by tilting slightly the axis of the domain).

In order to validate our extraction procedure, we implemented it in two different ways. First, we implemented the same algorithm as UR10 (ie: recursive, using breadth-first search as explained in UR10). Second, we implemented a non-recursive algorithm: we parse the whole cube to find a pixel above threshold which is not yet included into a structure ; we tag it and we scan the whole cube several times to tag the neighboring pixels of this growing seed which happen to be above threshold, until we find no new pixel to attach to this structure ; finally we reiterate to find another pixel not yet included in a structure until we fail to find any new pixel above threshold. The second algorithm is much more CPU time consuming, but keeps the memory usage constant and is easier to implement. We checked that both algorithms identify strictly the same structures for the two implementations in the low resolution cases.

Like in UR10, all the above quantities are found to exhibit power-law scaling with respect to structure linear size L i , with different scaling behavior in the inertial and dissipative ranges. The quantity X i , which could be any of L i , R i , V i , A i , P i scales as

X i ∝ L D X i
with different scaling exponents D X in the inertial and dissipative ranges, while the pdf of X i scales as

P(X i ) ∝ X -τ X i
with different scaling exponents τ X in the inertial and dissipative ranges. As an example, the scaling relations P i ∝ L D P i and P(P i ) ∝ P -τ P i are shown in figures 4.17 for the structures extracted from Run 12 (AD -OT), at the peak of dissipation, with a threshold of two standard deviations above mean value. The limits of the inertial and dissipation ranges are also shown (as used by UR10, see section 4.4.4). The upper limit of the dissipation range is just below the lower limit of the inertial range, while the lower limit of the dissipation range is ∼ 2 times the numerical resolution.

Comparison with UR10

In this section we compare the results of the statistical analysis of structures of high dissipation with those of UR10. These authors consider pure MHD and study the structures of high Ohmic dissipation or high enstrophy from Run 12 (AD -OT), at the peak of dissipation, with a threshold of two standard deviations above mean value. The dotted line shows the effect of adopting the slope found by [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF] instead of our own slope.

ε ω = Re -1 ω 2 , ω = ∇ × u -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 L i -3 -2 -1 0 1 2 3 4 P i
The results for the scaling exponents are shown in figure 4.18, for the case of the OT initial condition. The case of the ABC initial condition (not shown) is similar. The exponents are calculated from run 10 which in terms of initial condition and Reynolds number is similar with run III of UR10. As in the present paper, the snapshot analyzed is on the peak of dissipation. The structures of high dissipation are defined as connected sets of points having values of the Ohmic dissipation two standard deviations above the mean value. We keep the same definitions as UR10 for the inertial and dissipative ranges (cf. section 4.4.4).

Figure 4.18 shows that the agreement between our results and those of UR10 is not completely satisfactory. Although most three-sigma error bars are compatible and the errors are on the same order, there remains systematic differences, especially for our pdf exponents which appear to yield shallower pdfs than UR10. The least square method used by UR10 to estimate the slope of power-law pdfs is known to introduce some bias and the maximum-likelihood estimate method (MLE) should be used instead [see [START_REF] Clauset | Power-Law Distributions in Empirical Data[END_REF]. We computed the exponents with the MLE, and found them to be very close to our least-square values. We turned to explore the effects of some systematics due to the uncertainty on the boundaries of the inertial and dissipative range and the size of the bins used to produce the pdfs. We varied randomly these parameters with factors in an octave centered on their initial chosen values. The excursion of the resulting three-sigma error bars over a thousand of such realizations are plotted as thin red error bars on Figure 4.18. A slight displacement of the inertial or dissipative range boundaries incorporates (or leaves out) new data near the edge of the fitting intervals, where their leverage on the fitted slope is quite important. The resulting figure shows that such systematics can account for nearly all discrepancies with respect to UR10, except for the correlation between the total dissipation of structures and their linear size. We illustrate the corresponding discrepancy on Figure 4.17 where the dotted red line shows the slope followed by UR10 data: this line clearly falls below our data at small scales. This suggests that our smallest structure have a higher value of dissipation. This may perhaps be traced back to the slightly more refined dealiasing rule which we use. The picture is the same for the ABC runs, except for the error bar 0.02 on D A in Run I of UR10 (see their Table II) which is probably a typo as the value they quote does not correspond to the scatter displayed in their Figure 3.

To validate further our results, we compared the computed scaling exponents for the dataset corresponding to the temporal peak of total dissipation of the pure MHD run with ABC or OT initial condition (Run 9 and 10) with those computed from a snapshot taken one macroscopic eddy turnover time later, in the decay period of the turbulence. In agreement with the results of UR10, we find no statistical difference between the two snapshots (one-sigma error bars are compatible). Similarly, we compared the exponent values computed based on j 2 with those computed based on ω 2 for the dataset corresponding to the peak of total dissipation

D L D R D V D A D P τ L τ R τ V τ A τ P -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Inertial range, OT Uritsky et al. Present work D L D R D V D A D P τ L τ R τ V τ A τ P -1 0 1 2 3 4
Dissipative range, OT of the same runs (Runs 9 and 10). Again the exponents were seen to be compatible within one-sigma error bars, as in UR10.

Structures based on total dissipation

In this section we focus on the statistical analysis of structures extracted based on the total dissipation ε t = ε o + ε v + ε a for both pure MHD and AD MHD. The analysis based on the total dissipation is more relevant to the heating of the ISM because all three different types of dissipation can be important heating agents. AD has an additional specificity because the ion-neutral drift increases the effective temperature of the chemical reactions, but we do not consider this yet in the present work. In the following, we will note D X for the linear size exponents and τ X for the probability exponents of a characteristic X. All the structures discussed in this section are defined as connected sets of points having values of the total dissipation two standard deviations above the mean value.

In Table 4.2 we present the results of the structure extraction algorithm. The relative amount of dissipation contained by all the detected structures does not depend much on the Reynolds number, but the volume filling factor of the structures decreases as the Reynolds number increases. The presence of AD results in fewer detected structures than without AD (except for run 13). However, they fill roughly the same volume fraction, thus AD structures tend to be larger. This difference between the number of different structures in pure MHD and AD increases with the Reynolds number. Figure 4.19 gives a more detailed view of the fraction of total dissipation contained in structures with a value above a given threshold as a function of the volume fraction occupied by these structures. The curve rises steeply near the origin, so that 30 percent of the dissipation in contained in less than 3 percent of the total volume. The steepness at the origin is seen to be mainly due to the Ohmic heating: this is in line with the original picture of [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF] where AD forms sharp features in which Ohmic dissipation is favored.

Figures 4.20 and 4.21 summarize all results on the exponents for the total dissipation two standard deviations above the mean. Although most three-sigma error bars are compatible, in particular for the D X exponents which are unchanged with AD, some systematic differences exist for the pdfs exponents τ X . The pdfs exponents are in general steeper in the dissipation range: AD seems to favor more fragmented structures in the dissipative range. This confirms the tendency for more intermittency with AD that was suggested by the structure functions analysis. However, we see no clear cut tendency in the inertial range. If we discard the strong AD results (Re a = 10), the inertial range shows a behavior opposite from the dissipative range (shallower pdfs slope, ie: larger structures are favored when AD is present). However if we now look only at the MHD runs and the Re a = 10 runs, the inertial range sees no change in the τ X exponents for the OT Figure 4.19: We look at the subset of pixels above a given threshold of total dissipation in run 12 (AD-OT, Re a = 100) at the dissipation peak. For each value of the threshold, we plot the fraction of the total energy dissipation on this subset versus the volume of this subset (black curve). We also give the fraction of the total dissipation on this subset for each nature of dissipation (red: Ohmic, green: viscous, blue: AD).

runs, but bigger exponents for the ABC runs, in agreement with the dissipative range and contrary to the Re a = 100 runs... This complicated picture might perhaps not be genuine, as the huge systematics experienced for the comparison with UR10 show. However, the dependence of the intermittency statistics on the initial conditions (seen both in the structure pdfs slopes and in the structure functions exponents) points at their difference in magnetic helicity content. In the OT initial condition, the initial value of magnetic helicity is almost zero, and the equation of magnetic helicity evolution is unchanged by the inclusion of the AD term. As far as the effects of viscosity and resistivity are neglected, both the pure MHD and AD MHD solution evolve under the same constraint of zero magnetic helicity. A zero value for the magnetic helicity is an important constraint because it implies statistical reflection invariance, a property which the ABC runs will not share. In the ABC initial condition, the constraint of very low cross-helicity is broken by AD, which provides a source term in the cross-helicity equation.

The above analysis does not give any information on the relative amount of Ohmic, viscous and AD dissipation contained within each structure. To answer this question, we show in figure 4.22 scatter plots of the ratio of total Ohmic dissipation to total dissipation ε o /ε t versus the ratio of total AD dissipation to total dissipation ε a /ε t in the two OT runs 12 and 14. In all AD cases there is a tendency for the structures to cluster close to the line of zero viscous dissipation (dashed line in figure 4.22). This shows that within the structures of high total dissipation, viscous dissipation is relatively less important. The relative values of Ohmic and AD dissipation span however the whole spectrum, in contrast to the impression given by our RGB slices (see figures 4.7 to 4.9): this could be a genuine difference between the extreme dissipation events and the bulk of the dissipation shown on the RGB figures, or the extraction algorithm of connected structures could merge nearby sheets with different dissipation natures. We checked on a pixel by pixel scatter plot similar to Figure 4.22 that it is indeed a genuine difference. For the strongest AD runs, though, the intense dissipation structures tend to be predominantly due to AD heating (top left corner in Figure 4.22).

Concluding remarks

We performed a 3D numerical study of the structures of high dissipation in MHD turbulence, with the inclusion of ambipolar diffusion. At the Reynolds numbers studied, the total dissipation due to viscosity, resistivity and AD are of comparable magnitude.

Kinetic and magnetic energy spectra show that ambipolar diffusion enhances the turbulent energy to small scales at the expense of intermediate scales. This agrees with the idea of [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF] that AD can sharpen the magnetic gradients, but the effect is not strong enough to increase the total Ohmic dissipation rate in our simulations. Previous authors [START_REF] Li | Sub-Alfvénic Nonideal MHD Turbulence Simulations with Ambipolar Diffusion. I. Turbulence Statistics[END_REF] and [START_REF] Oishi | The Inability of Ambipolar Diffusion to Set a Characteristic Mass Scale in Molecular Clouds[END_REF] have examined the case of driven two-fluids compressible turbulence with a mean magnetic field and they find different results: [START_REF] Oishi | The Inability of Ambipolar Diffusion to Set a Characteristic Mass Scale in Molecular Clouds[END_REF] find no effect on the slope of the spectrum (for Re a = 2.5 and Re a = 5) while [START_REF] Li | Sub-Alfvénic Nonideal MHD Turbulence Simulations with Ambipolar Diffusion. I. Turbulence Statistics[END_REF] find that AD steepens it, although they don't display the results for their runs with Re a = 12 and Re a = 120, the only ones which have a in the computed range of scales as in our study. It should be noted that both these papers neglect Ohmic diffusion, and rely on truncation errors only to reconnect the field: the present work is the first to account for AD in the presence of controlled Ohmic and viscous dissipation. This is important because the dissipation physics can be quite different from the numerical dissipation as was demonstrated by Fromang and Papaloizou [2007], Fromang et al. [2007] in magneto-rotational turbulence.

As in [START_REF] Oishi | The Inability of Ambipolar Diffusion to Set a Characteristic Mass Scale in Molecular Clouds[END_REF], we fail to detect a significant change of regime in the spectra at the expected AD length scale a , but in our simulations it happens at a greater length scale. This length scale a is predicted from the balance between the moduli of the Fourier coefficients for the inertial e.m.f. and the AD e.m.f.. We would underestimate a if AD was more coherent in time than the advection of the field, and so we conjecture that AD terms have a greater coherence time than the inertial terms.

In our simulations, we observe that AD shuts off the Lorentz force at small scales: this shifts the peak of the AD heating power spectrum to larger scales. The position of that peak defines a scale which seems to match the characteristic thickness of the sheets where AD heating is strong. This scale * a might be revealed by the characteristic chemistry of AD heating where neutral-ion endothermic reactions are favored (e.g.: CH + or SH + formation, see [START_REF] Godard | Models of turbulent dissipation regions in the diffuse interstellar medium[END_REF]. Table 4.3 sums up the characteristics of various environments of the ISM. We identify the integral length scale and the r.m.s. velocity in our simulations to their corresponding physical values in each considered ISM components to apply our results. Although we find that Re -1 a in the ISM varies from 10 -3 to 10 -2 , we use * a = 0.6 as measured in our AD simulations with Re -1 a = 0.01 to estimate the scale l 0 * a of the AD diffusion heating. The AD dissipation heating is always about a fourth of the typical scale whereas l 0 . a can be much smaller.

The qualitative picture of the dissipation field suggests that it is dominated by intermittent sheet-like structures which alternate with large voids of low dissipation. The sheets of various dissipative nature (Ohmic, viscous or AD) appear to be clearly separated, except for the highest dissipation rates, where viscosity fades out and Ohmic and AD heating blend. AD heating sheets are often seen to sandwich much thinner regions of strong Ohmic or viscous dissipation, as in the simple case studied by [START_REF] Brandenburg | The formation of sharp structures by ambipolar diffusion[END_REF]. The high degree of intermittency is confirmed by the computation of the structure function exponents for the velocity and the magnetic field, as well as by the pdf of the total dissipation, which exhibits a log-normal core and a strong power-law tail for high values of the dissipation. We compared the statistics of the structures of strong dissipation with those of UR10 and we obtain only marginal agreement, probably because of the systematics linked with the definition of the inertial and dissipative ranges. The statistical analysis of structures of high total dissipation reveals the highly intermittent nature of the dissipation field, as more than 30% of the dissipation takes place in less than 3% of the volume. No significant difference in the scaling laws between pure MHD and AD MHD was found, but the slope of the power-law pdfs is affected in the dissipative range, with a statistical preference towards more fragmented structures.

In future work, we intend to make progress towards a more realistic picture of the ISM, relaxing the incompressible hypothesis, with the final aim of including realistic cooling. We also hope our statistical results will provide new ways to approach the observed characteristics of the ISM, intermediate between direct post-processing of 3D numerical simulations and the building of line of sights with elementary models such as shocks or vortices [START_REF] Godard | Models of turbulent dissipation regions in the diffuse interstellar medium[END_REF].

95

Chapter 5

Energy dissipation in compressible turbulence

Turbulent flows in diffuse clouds are characterized by a wide range of sonic Mach numbers, from 0.1 to 10 or even more [START_REF] Elmegreen | Interstellar turbulence i: Observations and processes[END_REF]. This means that at least for some regions of diffuse clouds, models based on the incompressible equations of motion are insufficient and one should use compressible models instead. In this chapter, we investigate the properties of homogeneous and isotropic compressible magnetohydrodynamic turbulence as a model of flow inside a diffuse cloud.

Introduction

In a compressible flow, the sonic Mach number is defined as

M s = u 0 c s (5.1)
where u 0 is a characteristic velocity scale, for example the root mean square velocity, and c s is the speed of sound. This number expresses the relative importance of compressibility effects. As mentioned above, the root mean square sonic Mach number of the flow in the ISM can be anywhere from 0.1 to 10 or even more [START_REF] Elmegreen | Interstellar turbulence i: Observations and processes[END_REF]. In a weakly compressible flow, the relative density perturbation grows as the square of the sonic Mach number [START_REF] Landau | Fluid Mechanics[END_REF]:

∆ρ ρ ∝ M 2 s .
(5.2)

As a consequence, flows with a sonic Mach number M s 0.1 can be regarded as practically incompressible, as the relative density variations ∆ρ/ρ are of the order or 1%. Thus, some regions of the diffuse interstellar medium are governed by incompressible dynamics. However, incompressible models are certainly insufficient for flows with M s 0.5 and above, because relative density variations can be as high as 25%, and compressible models should be used in this case. The Alfvénic Mach number which is defined as

M a = u 0 v a , v a = |B 0 | √ 4πρ (5.3)
where v a is the Alfvén velocity, describes the relative importance of the magnetic field. Here, B 0 can be taken as either the global mean field (if present), or a local average field. The relative importance of the magnetic field can be also characterized by the beta parameter, defined as

β = p E M (5.4)
where p is the thermal pressure and E M the magnetic energy density. Due to the approximate equipartition between kinetic and magnetic energy, the Alfvénic Mach number M a and the β parameter are of order unity in the diffuse interstellar medium. There have been many numerical studies of compressible turbulence in the literature, but the overwhelming majority do not explicitly consider physical viscosity and resistivity but instead rely on numerical dissipation. Studies of forced isosthermal MHD turbulence include [START_REF] Vestuto | Spectral Properties of Compressible Magnetohydrodynamic Turbulence from Numerical Simulations[END_REF], [START_REF] Kritsuk | Simulations of Supersonic Turbulence in Molecular Clouds: Evidence for a New Universality[END_REF] and [START_REF] Federrath | Comparing the statistics of interstellar turbulence in simulations and observations. Solenoidal versus compressive turbulence forcing[END_REF]. [START_REF] Kritsuk | Comparing Numerical Methods for Isothermal Magnetized Supersonic Turbulence[END_REF] performed a detailed study of decaying compressible MHD in the context of comparison of different numerical codes, while [START_REF] Hennebelle | On the structure of the turbulent interstellar atomic hydrogen. I. Physical characteristics. Influence and nature of turbulence in a thermally bistable flow[END_REF] performed a study of two-phase compressible MHD turbulence in the context of the dynamics of the cold neutral medium (CNM). Studies of compressible ambipolar diffusion magnetohydrodynamics using the two-fluid model under the heavy-ion approximation include [Oishi and Mac Low, 2006a, Li et al., 2006, McKee et al., 2010] and [Li et al., 2012]. [START_REF] Federrath | The Turbulent Dynamo in Highly Compressible Supersonic Plasmas[END_REF] performed a study of magnetic Prandtl number dependence of the compressible turbulent dynamo by explicitly taking into account physical viscosity and resistivity.

In our simulations, we chose two values of the rms sonic Mach number: M s = 1 and 4, while the Alfvénic Mach number is always unity (cf. table 5.1). These choices correspond to reasonable assumptions for trans/supersonic, trans-alfvénic flow in diffuse clouds (see table 4.3 for typical values of the ISM). Our computational resources allowed us to explore two different kinetic Reynolds numbers: Re = 6.2 × 10 3 and Re = 8.9 × 10 3 . The magnetic Prandtl number is set equal to unity so the magnetic Reynolds number Re m is equal to the kinetic Reynolds numbers Re. We present here only the low resolution runs (Re = 6.2 × 10 3 , at 512 3 grid-based resolution) which are compared to the incompressible runs of the previous chapter with similar Reynolds numbers (at 256 3 spectral resolution).

The initial conditions of the compressible simulations, just as the incompressible simulations of chapter 4 were based on macroscopic fields derived from the ABC flow and the Orszang-Tang vortex, as defined by relations (4.6) or (4.8).

Id.

N are computed using centered differences. The implementation of physical viscosity and resistivity has been validated by tests of the decay of Alfvén waves [START_REF] Lesaffre | Exact shearing box solutions of magnetohydrodynamic flows with resistivity, viscosity and cooling[END_REF]. The minimum diffusion time-scale

τ d = min ∆x 2 ν , ∆x 2 η (5.5)
where ∆x is the distance between two grid points, is calculated and used in addition to the inviscid, nonresistive CFL relation to constrain the timestep. In all simulations a CFL number equal to 0.7 was used.

In our compressible simulations, the kinematic viscosity ν = µ/ρ is a fixed given quantity, which implies that the dynamic viscosity µ varies in space and time. This choice is motivated by the fact that, in a shock, the width of the viscous layer is proportional to a viscous length given by λ v = ν/u = µ/(ρu) where u is the shock speed. For fixed µ, this width can reach very low values in high density regions, and there is no a priori guarantee that our grid will be able to resolve it. Thus, with a fixed value of the kinematic viscosity ν we impose a single length scale for viscous dissipation in shocks, and are therefore able to resolve it more uniformly. For the same reason, we set our Mach number M s = 4 simulations by decreasing the sound speed c s at fixed initial r.m.s. velocity u r.m.s. rather than by increasing u r.m.s. at fixed c s . In other words, we set simulations at higher Mach number while keeping the Reynolds number fixed.

Note that in the simulations with linear resolution N = 1020, we use a value of the viscosity which is almost equal to the value used for the N = 512 incompressible spectral simulations of chapter 4. This is because a spectral method is able to achieve a spatial resolution two times higher than the one achieved by a finite volume method using the same number of grid points [START_REF] Gottlieb | Numerical analysis of spectral methods[END_REF].

Time evolution

In this section we discuss the time evolution of average quantities such as the mean dissipation rate. Our simulations are freely decaying and start from smooth initial conditions, so that variables that depend on the velocity and magnetic field derivatives have small values at the beginning, and are expected to increase as the turbulent structure of the field grows.

E tot = 1 2 ρu 2 + 1 8π B 2 + p log ρ (5.6)
at each snapshot of the simulation. We then compute its time derivative -Ėtot (solid lines) which provides the total actual irreversible heating rate during each interval between two snapshots. The physical heating rate ε phys (dashed lines) is computed from

ε phys = ε v + ε o (5.7)
where ε o and ε v are the Ohmic and viscous heating rates (the viscous heating rate now includes the full compressible viscous tensor). The physical heating rate is indistinguishable from -Ėtot in the case M s = 0 (incompressible runs), which confirms that spectral methods are very performant at resolving the dissipation. But in the compressible runs, (M s = 1 and 4), the physical heating rate is less than half the actual energy loss. This is a sign that a lot of numerical dissipation due to the scheme is present. In order to study dissipation in these simulations, we hence need a means to estimate the dissipation produced by the scheme and we devise new methods to do so in the appendix. We denote ε corr the total energy dissipation rate as recovered energy from method 1 in the appendix. Figure 5.2 shows that this quantity slightly overestimates the total actual energy losses (the worst case is OT initial condition, M s = 1 where the overestimation is up to 40%). Figure 5.2 also illustrates that in runs with higher Mach number, the dissipation peak appears earlier. Finally, it is worth noting that the temporal peak of total dissipation is not far from the corresponding peak of physical dissipation, which we used to estimate the actual dissipation peak. The physical dissipation sometimes undergoes secondary local peaks, and we in fact settle on the last of these peaks to record our dissipation peak snapshot.

In figure 5.3 we plot the time evolution of the ratio of viscous to total physical dissipation rates ε v /( ε v + ε o ). In general, the ratio remains close to 1/2 throughout the time evolution, with an equal amount of viscous and Ohmic dissipation. There is a slight relative increase of viscous dissipation as the Mach number increases, but the difference is not very pronounced.

In the case where the dynamic viscosity µ is constant, the mean viscous dissipation rate ε v can be decomposed into a compressible and a vortical part [START_REF] Kritsuk | The Statistics of Supersonic Isothermal Turbulence[END_REF]]

ε vc = 4 3 µ (∇ • u) 2 , ε vω = µ ω 2 .
(5.8)

Note that this decomposition is only valid for the mean viscous dissipation and not for the pointwise viscous dissipation. The decomposition is not strictly possible in the case where the kinematic viscosity ν is constant, because there is an additional term proportional to ∇ρ. However, we have found that this term is small and the analogous decomposition

ε vc = 4 3 νρ(∇ • u) 2 , ε vω = νρω 2 .
(5.9) serves as a proxy of the relative importance of the vortical and dilatational components of the mean viscous dissipation.

In figure 5.4 we present the time evolution of the ratio of dilatational to total viscous dissipation rates ε vc /( ε vc + ε vω ). The ratio peaks early (before the temporal peak of total physical dissipation rate) at 0.7, which may correspond to a transfer of energy from rotational to dilatational degrees of freedom (the initial conditions do not contain any energy in the dilatational modes). The ratio then decreases to 0.2 or even less, a value which is quite low. Again, this may correspond to a preferential decay of the dilatational modes. There also appears to be surprisingly little difference between the various Mach numbers and initial conditions. For most of the time evolution, the physical heating is dominated by the rotational component, not the dilatational one. In other words, the physical dissipation is dominated by the contribution from shearing sheets and not from shocks. However, we recall that a great deal of the dissipation in our simulations is due to the numerical scheme (more than half of the total, ie: more than twice the total viscous dissipation ε vc + ε vω ). If most of the numerical heating comes from shock contributions as in our tests from the appendix, then we should conclude on the contrary that shocks dominate the dissipation in our simulations... To summarize, our method for recovering energy losses from the scheme is currently unable to discriminate between vortical and dilatational losses. Similarly, it is unable to detect how much magnetic energy is degraded vs. how much kinetic energy is degraded. Hence more work needs to be done before we can safely assess which is the main actor of dissipation in the ISM: shearing sheets, current sheets or shocks ?

The density field

In compressible simulations, the absence of the solenoidal constraint ∇ • u = 0 allows the density field ρ(x) to vary in space and time. In this section, we investigate the spatial variation of the density field.

In figure 5.5, we present density cut color maps from runs 1-4. In the top row, which corresponds to a Mach number M s = 1, the appearance of the density field is qualitatively similar to that of a passive scalar in incompressible hydrodynamic turbulence. In the bottom row, which corresponds to an increased Mach number M s = 4, the density contrast is much greater and sharp density ridges strike out.

We define the column-density in the direction of observation z as:

N (x, y) = 1 L box L box 0 ρ(x, y, z) dz.
(5.10)

Note that since the average density ρ is always equal to ρ = 1, we similarly have N = 1. The column density maps for runs 1-4 are shown in figure 5.6. In the top row which corresponds to a Mach number M s = 1 there is relatively less sheet-like structure in comparison to the bottom row which corresponds to a Mach number M s = 4. In the bottom row, higher compression yields higher column-density contrasts and column-density peaks appear at the crossing between density ridges. The pdfs of the density are presented in the left panel of figure 5.7. In the cases corresponding to M s = 1, the pdfs can be adequately fitted by a log-normal distribution where the standard deviation σ s depends on the Mach number according to the relation

P(s) = 1 2πσ 2 s exp (s -s ) 2 2σ 2 s , s = ln(ρ/ ρ ) (5.11)
σ 2 s (M s ) = ln(1 + b 2 M 2 s ) (5.12)
with b a free parameter. In the cases corresponding to M s = 4 the fit is much worse because the distribution is skewed towards lower densities. Federrath et al. [2010a] report values of b from 1 to 0.35 in forced compressible magnetohydrodynamic turbulence, depending on the nature of the forcing (solenoidal versus compressive). In our case of freely decaying turbulence, we find somewhat lower values of b, in the range 0.29 -0.34 depending on the initial conditions. In contrast with the pdfs of the density, the pdfs of the column density are close to a lognormal distribution even for the cases corresponding to M s = 4. The integration process from densities to column-densities smears out the extreme values, and we naturally find much lower values for b, in the range 0.06 -0.12 depending on the initial condition.

The dissipation field

Here, we investigate the spatial structure of the dissipation field qualitatively starting with the same method as that used in chapter 4 to produce figures 4.7-4.9, that is we only show the "bulk" of the dissipation: the bottom 10%, calculated based on the cumulative distribution function, is truncated to black, and the top 10% are saturated.

There is not much qualitative difference between figure 5.8 and figures 4.7-4.9 which correspond to the M s = 0 (incompressible) case. Dissipation is concentrated on thin, sheet like structures whose projections appear as filaments in figure 5.8. Ohmic and viscous dissipation appear colocated on the same structures, and as the Mach number is increased from M s = 1 in the top row to M s = 4 in the bottom row, the structures become fewer and scarcer and the black areas between them larger. Although the decomposition of the total physical dissipation in a compressive and a solenoidal part is strictly valid only in the case of global mean values and constant dynamic viscosity, the corresponding local quantities

ε vc (x) = 4 3 νρ(∇ • u) 2 , ε vω (x) = νρω 2
(5.13) can serve as a useful proxy for the relative importance of dilatational and rotational effects in the dissipation. Indeed, we have observed that ε vc + ε vω was highly correlated to ε v by looking at the two corresponding maps. In figure 5.9, Ohmic dissipation ηJ 2 is assigned to the red channel, viscous rotational dissipation ε vω is assigned to the green channel and viscous dilatational dissipation ε vc is assigned to the blue channel. Again, only the "bulk" of the dissipation is shown: the bottom 3% based on the cumulative distribution function is truncated to black and the top 3% is saturated.

Comparing with figure 5.8, we observe that the different types of dissipation are now less colocated on the same structures. The sheets of intense blue colour (strong dilatational viscous dissipation), whose integrations in thin slices appear as blue strokes in the color-maps of figure 5.9, are presumably correlated with the position of shock waves. In figure 5.2, we observed that throughout the time evolution of the simulation, a significant amount of the total dissipation is due to numerical dissipation inherent to the numerical scheme. Here, we investigate the spatial distribution of numerical dissipation by applying method 1 of the appendix. In figure 5.10, Ohmic dissipation ηJ 2 is assigned to the red channel, viscous dissipation is assigned to the green channel and numerical dissipation ε

corr. v + ε corr. o -(ε v + ε o
) is assigned to the blue channel. The color truncation is the same as in figure 5.9 (top and bottom 3%). The dissipative structures still assume sheet-like shapes, as the incompressible case.

The pdfs of the total dissipation plotted in figure 5.11 show that there is a general tendency for the pdf to broaden with increasing Mach number. This means that there is a higher probability for both low and high values of the local dissipation. In particular, the tail at large dissipation rates is more pronounced as the Mach number is increased: the intermittency of the dissipation field is increased at higher Mach numbers. The actual picture is probably even more clear cut: although our method to recover energy losses from the scheme performs well to recover the spatially integrated irreversible heating, we note in the appendix that we smear out the highest dissipative structures. This means that very likely we underestimate the statistical weight of the highest dissipation rates and the compressible runs should have even shallower tails on the high dissipation side.

with the lag l expressed in pixels or resolution elements.

In the following, the centroid velocity increment is computed by weighting the cells by the local dissipation rate in the integral along the line of sight. This is for the following reason. The lines that have been used to build the statistics of the CVI are 12 CO lines in diffuse molecular gas. In this gas component, the production of CO is likely to be dominated by non-equilibrium chemistry triggered by turbulent dissipation in intense velocity-shears and vortices as in [START_REF] Godard | Chemical probes of turbulence in the diffuse medium: the TDR model[END_REF], or in low-velocity shocks [START_REF] Lesaffre | Low-velocity shocks: signatures of turbulent dissipation in diffuse irradiated gas[END_REF] because the large observed abundances and isotopic ratios cannot be understood in the framework of UV-driven chemistry [START_REF] Hily-Blant | Small-Scale Dissipative Structures of Diffuse ISM Turbulence: I -CO Diagnostics[END_REF][START_REF] Hily-Blant | Dissipative structures of diffuse molecular gas. III. Small-scale intermittency of intense velocity-shears[END_REF]. We therefore anticipate that this molecular species in that medium be a good tracer of the regions of intense turbulent dissipation. We have compared the results obtained with different weightings or no weighting, and the results are only slightly different. The energy-dissipation averaged centroid velocity increments are defined as:

CV ε (x, y) = 1 L d ε LoS (x, y) L d 0 ε(x, y, z) u z (x, y, z) dz and CVI ε (x, l) = |CV ε (x + r) -CV ε (x)| C l |r| = l (6.4)
where ε LoS (x, y) is the LoS integrated energy dissipation rate:

ε LoS (x, y) = 1 L d L d 0 ε(x, y, z) dz (6.5)
The above weighted quantities are expected to correspond better to the quantities which are actually observed in emission.

The magnetic field (intensity and orientation) can be measured by different techniques depending on the ISM phase studied. A comprehensive review of all the methods is given in [START_REF] Crutcher | Magnetic Fields in Molecular Clouds[END_REF]. In the following, we discuss one observable: the polarization angle of the thermal dust emission. It is a tracer of the orientation of the magnetic field projected on the plane-of-the-sky (PoS) in the cold ISM. The dust thermal emission is linearly polarized because dust grains are elongated. They spin at supra-thermal velocities about an axis perpendicular to their long-axis, and have their angular momentum precessing about the local magnetic field. They eventually get their angular momentum aligned with the field. As a result, the dust thermal emission is polarized along a direction perpendicular to the local magnetic field. The alignment mechanisms are still debated, but the most efficient is due to radiative torques exerted on dust grains by stellar photons ( [START_REF] Hoang | Radiative torque alignment: essential physical processes[END_REF], see the review of [START_REF] Andersson | Interstellar grain alignment[END_REF]).

The observables are the three Stokes parameters I, Q and U where I is the intensity of the dust thermal emission and Q and U involve the positional angles of the magnetic field with respect to the plane of sky. In a numerical simulation, in order to obtain information on the local structure of the magnetic field, one defines the local polarization angles γ and φ as follows:

γ(x, y, z) = tan -1 (B z (x, y, z), B x (x, y, z) 2 + B y (x, y, z) 2 ) φ(x, y, z) = tan -1 (B y (x, y, z), B x (x, y, z)) (6.6)
The quantities Q and U are then defined in terms of integrals of these angles along the line of sight:

Q(x, y) = 1 L d L d 0 cos(2φ(x, y, z)) cos 2 (γ(x, y, z)) dz and U (x, y) = 1 L d L d 0 sin(2φ(x, y, z)) cos 2 (γ(x, y, z)) dz.
(6.7)

Finally, the integrated polarisation angle ψ is defined as:

ψ(x, y) = 1 2 tan -1 (U, Q). (6.8)
The increments of the polarization quantities are defined in a similar way as for the centroid velocity

δQ 2 (x, l) = |Q(x + r) -Q(x)| 2 C l , δU 2 (x, l) = |U(x + r) -U(x)| 2 C l and δψ 2 (x, l) = |ψ(x + r) -ψ(x)| 2 C l .
(6.9)

All the above increments are expected to show increased values in regions where the original LoS integrated variables exhibit strong gradients, such as vortices and current sheets, especially for small lags. Figure 6.1 shows an example of an observational CVI map in the Polaris flare computed in a map of CO(J=2-1) emission. Figure 6.2 displays the polarization angle dispersion function S (here noted δψ) as measured by the Planck satellite over a large portion of the sky.

Qualitative link with structures of high dissipation

The increments of the diagnostic quantities introduced above are expected to reach high values in regions where there are rapid variations of the LoS integrated velocity and the magnetic field orientations. If a structure of high dissipation is present in a part of the computational box, the strong gradients at its edges could be detected by high values of the increments of the diagnostic quantities if not erased but the LoS integration. In this section we present a qualitative link between structures of high dissipation and increments of projected diagnostic quantities.

In figures 6.3 to 6.5, integrated dissipation field and increment maps of the diagnostic quantities CVI ε , δQ and δψ are displayed for a lag l = 3 pixels and integration depths of 1/64th, 1/8th and the whole box size. In the dissipation color-map, we use the same display mode as in chapter 4: only the bulk of the dissipation is shown: the top and bottom 10% of extreme values, based on the cumulative distribution function of the total dissipation are truncated. Increment maps present worm-like structure which, at least qualitatively, appears to be correlated to the integrated dissipation.

Figure 6.6 shows color maps of the integrated dissipation field overlaid with contours of the diagnostics CVI ε , δψ, δU and δQ, for a lag l = 3. The contours of the diagnostics closely follow the integrated structures of high dissipation. We also note that the velocity and magnetic field diagnostics rarely coincide spatially but they seem to be complementary in the sense that they trace different parts of the same structures of enhanced dissipation. Still, not all structures of high dissipation are captured. This is because of projection effects which cause some variations of velocity or magnetic field components to remain undetected. The CVI diagnostic traces only the sharp variations of v z in the plane of the sky, and is unable to capture any variation of the other velocity components nor variations along the line of sight. Similarly, the δψ diagnostic captures variations of the orientation of the projection of the magnetic field in the plane of the sky. δU , δQ and δψ all trace variations of the orientation of the magnetic field, but with different nuances as explicit in equations (6.7) to (6.9).

In a qualitative sense, the correlation of the contours with the integrated dissipation is better for smaller integration depths. Integration over larger volumes blurs out the correlation, although not entirely. The subset of the dissipation field that is traced by the increment diagnostics is smaller and smaller as the depth of the integration increases. Qualitatively also, the contours of the diagnostics for the different integration depths are clearly more coherent spatially at small depth than at larger depths.

Even though the structures in the contours of δψ at large integration depth do not always correspond to projected structures of high dissipation, they still testify for coherent structures in the dissipation field. Indeed, it is shown below that if the phases of the velocity and the magnetic field are mixed randomly, the contours of δψ and CVI lose most of their spatial coherence.

Quantitative link with structures of high dissipation

In this section we attempt a quantitative analysis of the correspondence between structures of high dissipation and increments of projected diagnostic quantities. One can introduce a quantitative criterion of how well a given diagnostic captures the dissipation by defining its efficiency as the ratio of the surface-averaged integrated dissipation contained within a given contour C of the diagnostic to the volume-averaged dissipation over the whole plane of sky P : If the diagnostic captures the integrated dissipation efficiently, E eff is greater than unity, whereas for a randomly placed contour with no correlation with the integrated dissipation it is close to unity. If the contour C avoids the regions of high dissipation, E eff is less than unity.

E eff = ε LoS C ε LoS P (6.10)
In the following, increment maps are normalized by their maximum value and contours are labeled by their level. Thus, contour level 0.5 is the contour at half the maximum value of the increment on the whole increment map. For a level close to one, the contour focuses on the strongest increment structure in the box: the efficiency is zero if the structure is off-center or very large if it corresponds to an actual structure of high dissipation. For a level close to zero, the contour encompasses almost the whole box and the efficiency is one. In figure 6.7, the efficiency as a function of the contour level is increasing for the lowest levels, becomes staggered due to the low statistics of the highest contour levels. At large values of the contour level, only the most intense structures remain. If the peaks of these structures are located at an (x, y) position with a high value of integrated dissipation, a high efficiency is obtained. If the structures are slightly off the projected structures of high dissipation, the efficiency obtained is low. Hence we consider only contour levels ≤ 0.5. Figure 6.8 shows that, as expected, E eff peaks at small lags and for a small integration depth. At larger depths, the diagnostic power is blurred by the projection effects.

The efficiency depends on the nature of the dissipation. For instance, figure 6.9 illustrates the efficiency of CVI ε at tracing various types of dissipation. As expected, it best traces viscous dissipation for small integration depths. For large integration depths it is a weak and unselective tracer of dissipation. Similarly (not shown) δQ, δU and δψ are good tracers of Ohmic dissipation at small integration depths. At large depths, δU and δQ become as weak and unselective tracers as CVI ε . However, δψ at large depths has an efficiency of one (i.e.: it is not efficient at detecting dissipation) for all types of dissipations.

Efficiency for a given nature of dissipation can also depend on the tracer. In figure 6.10, we examine the efficiency of our four diagnostics to trace the integrated Ohmic dissipation for two projection depths. As expected, for small integration depths, δQ, δU and δψ are significantly better tracers than CVI ε . At large integration depths, δQ, δU become comparable to CVI ε and as mentioned above δψ is inefficient. For viscous dissipation, the quantity CVI ε is the best diagnostic and for A.D. dissipation, there is no good diagnostic amongst those we have studied (not shown here).

Intermittency exponents: importance of initial conditions

One classical signature of the intermittent structure of a field is the presence of non-Gaussian, often exponential, wings in the pdf of its increments. Figure 6.11 presents the pdfs of the variables δψ, δQ, δU and CVI for various lags l, integrated over the whole box. As commonly found for increments of intermittent fields, the pdfs are Gaussian at large lags and develop fat wings at small lags. In other words, extreme events at small scale are not so rare.

In this section we discuss the intermittency exponents of the structure functions of the projected diagnostic quantities. These exponents provide a means to quantify the spatial intermittency properties of the various projected fields.

The structure function of a projected diagnostic quantity, for example the polarization angle ψ, is defined as the moment of the increment with respect to the corresponding probability distribution 6.11) where P l (δψ) is the pdf of the increment δψ for a specific lag l. The structure functions of the other quantities are defined similarly

S ψ p (l) = δψ p (x, l) = ∞ -∞ |δψ| p P l (δψ) dδψ (
S CVε p (l) = |CVI ε | p (x, l) , S U p (l) = |δU | p (x, l) , S Q p (l) = |δQ| p (x, l) (6.12)
with an analogous definition for the weighted quantities. In practice, the pdfs P l (δX) for a variable X are computed as pdfs of signed values of all the increments with lags between l and l + 1 pixels. Figure 6.11 shows the pdfs of P l (CVI ε ), P l (δU ), P l (δQ) and P l (δψ) integrated over the whole box for various lags l. As the lag l is decreased, the tails of the pdfs become more pronounced. This is identified as a signature of small-scale intermittency, because pronounced pdf tails are associated with an increased probability of extreme events. For large lags, the pdfs become increasingly Gaussian, with the exception of P l (δψ) which becomes uniform (polarisation angles are independent from each other at large distances).

One way to quantify the intermittency properties of the above pdfs is to measure their moments of order p, which are none other than the structure functions introduced above. Because high order moments sample the tails of the pdf preferentially, high order structure functions are expected to characterize the small scale intermittency present in the projected fields. In the inertial range, structure functions of a variable X are expected to scale as

S X p (l) ∝ l ζ X p (6.13)
where the exponent ζ X p is in general a non-linear function of the order p (cf. section 2.5). These exponents can be measured as logarithmic slopes:

ζ X p = d log S X p (l) d log l . (6.14)
In numerical simulations or experiments where the extent of the inertial range is limited, one can try to extend the range of scales where a power-law relation holds by applying the method of extended self-similarity (ESS) introduced by [START_REF] Benzi | Extended self-similarity in turbulent flows[END_REF]. This method consists in measuring the relative exponent

ζ X * p = d log S X p (l) dS X * (l) (6.15)
where S X * (l) is a reference structure function . In our case, we have adopted S X * (l) = S X 3 (l), inspired by the case of incompressible hydrodynamics, where this choice is supported by the 4/5-ths law (cf. chapter 4). Using these results, we compute the structure function exponents as the average slope of these lines and estimate the error using standard linear regression. Figure 6.13 presents the ESS intermittency exponents for all diagnostics and all our simulations (see chapters 4 and 5). It is customary to characterize the level of intermittency by the departure from the straight line p/3 (commonly referred to as K41). Both ABC and OT simulations span a large range of possible intermittency exponents between K41 and the generalised She-Lévêque model with g = 3 and C = 1, just as we previously found for the 3D intermittency exponents in the corresponding simulations.

There is however no one to one correspondence between the 3D exponents and their projected counterparts. Nevertheless, several trends emerge. The intermittency curves almost never cross each other. All diagnostics display more intermittency for high A.D. Reynolds number (Re AD = 100) than for low A.D. Reynolds number (Re AD = 10). In all cases, CVI ε displays more intermittency than CVI. δQ and δU always follow each other except in two cases (ABC for incompressible MHD, and OT at Mach 4). For the incompressible runs with OT initial conditions, CVI also follow δQ and δU , while for ABC initial conditions, there are strong differences in the values of the exponents. Higher Mach numbers naturally display more intermittency than lower Mach numbers. To summarise, these figures illustrate the high sensitivity of observable intermittency exponents to all parameters of the simulations, and in particular to the initial conditions.

The detailed comparison of these intermittency exponents to those observed by [START_REF] Hily-Blant | Dissipative structures of diffuse molecular gas. III. Small-scale intermittency of intense velocity-shears[END_REF] is made even more difficult since both our computed diagnostics CVI and CVI ε are greatly simplified with respect to the actual line formation. The CO lines result from a combination of non-linear thermal and chemical processes at the sites of dissipation, followed by complex radiative transfer and excitation of the molecule, not to mention instrumental noise. Nevertheless, our computations confirm that the CVI as measured are indeed a signature of intermittency, but it is too early to use them to discriminate between models.

Importance of coherence

In this section we investigate whether the diagnostics quantities are faithful tracers of the coherence of the structures of high dissipation present in the flow. A structure of high dissipation is considered coherent if it is localized in space and survives for several eddy turnover times. Coherence implies correlations between the phases of the Fourier coefficients of the velocity and magnetic fields, and it is also linked to intermittency [START_REF] Frisch | Turbulence: the legacy of A. N. Kolmogorov[END_REF].

For the purpose of investigating the role of coherence, we introduce the weakly phase-mixed fields u , b and the strongly phase-mixed fields u , b . The weakly phase-mixed fields are defined as

f (k) = f (k) exp [2πiφ(k)]
(6.16)

Here, f represents the Fourier transform of any component of the velocity or the magnetic field and φ(k) is a random phase of uniform distribution in [0, 2π] that depends on the wave-vector. In weak phase mixing, the random phase function φ(k) is the same for all components of the velocity and the magnetic field. Strong phase mixing is defined similarly, but now the random phase function is different for each component. We have found that weak phase mixing is sufficient to remove coherence, and for this reason we only discuss weak phase mixing results in this subsection. Figure 6.14 shows a color-map of the integrated dissipation field in the case of weak phase mixing. The dissipation field has clearly lost the sheet-like structure which appears in figures 6.6. The alignment of the vectors b and j = ∇ × b is destroyed, and this results in enhanced ambipolar diffusion dissipation, which is proportional to (j × b) 2 . This is why AD dissipation dominates the two other dissipation processes in this case. We have found that the above features do not change with the integration depth. Similarly, the viscous and ohmic dissipation fields have lost their sheet-like structure (figures 6.15 and 6.16). Figure 6.17 shows a color-map of the weighted centroid velocity increment CVI ε (x, l) for a lag l = 3, integrated over the whole box. Weak phase mixing has removed all filamentary structure from the CVI, except for that present at very small scales, and has reduced its intensity considerably. In addition, correlation with any of the dissipation fields is completely lost. Figure 6.18 shows a color-map of the weighted increment δψ(x, l) for a lag l = 3, integrated over the whole box. There is still some worm-like structure visible, but the structures are much smaller, more disrupted and all correlation to the integrated dissipation fields is lost.

The pdfs of the weighted centroid velocity increments (figure 6.19) still show exponential wings at small lags, but these wings are much steeper/narrower that those corresponding to the unmixed field. This suggests that an important part of the intermittent structure of the velocity field is lost after the application of weak phase mixing. This is further supported by figure 6.20, where we present intermittency exponents for all variables in the case of weak phase mixing. Weak phase mixing is sufficient to remove all structure from the fields and bring the exponents very close to K41 values, which correspond to a non-intermittent field.

The above results indicate that the long worm-like structures observed in the diagnostics maps are the signature of coherent structures of high dissipation. In addition, the weakening of the tails of the pdfs of the centroid velocity increments after phase mixing implies that these structures are linked to intermittency. Thus, fields with a power-law energy spectrum and random phases, of which fractionated Brownian motion is an example, are not a good representation of the structures of high dissipation.

Conclusions and perspectives

Turbulence in the interstellar medium is a complex, multi-scale, multi-physics phenomenon where a wide range of processes take place: chemistry, radiation, hydrodynamics, magnetic field... In the context of this thesis, we focused on the magnetohydrodynamics of the turbulent motion, neglecting additional phenomena such as chemistry and radiation. We hope that even while limiting our study to the hydrodynamic process of turbulent energy dissipation, we are still able to make valuable contributions to the study of the complex ISM life-cycle.

Turbulence is a most challenging theoretical problem of nonlinear statistical physics, and analytical difficulties due to the interplay of nonlinearity and stochasticity leave its theory in an incomplete state, with scarce exact results. Classical turbulence theory is not able to contribute much by itself in the study of turbulent energy dissipation in the ISM, but it leads up naturally to the formulation of the much more useful theory of intermittency. In the context of this theory, a wide range of models are able to predict local bursts in the dissipation and their impact on the statistical properties of the velocity and magnetic fields, such as structure function exponents. Even though they focus on the properties of the inertial range and not on the properties of the dissipative range, these models are able to contribute to the study of ISM chemistry. This is because quantities such as the mean time between two bursts of dissipation along a Lagrangian trajectory of a fluid parcel, which are very important for chemical models, are properties of the inertial range and can be described in terms of standard intermittency models.

In the future, we plan to use the statistical analysis of the dynamical and geometrical properties of the structures of high dissipation to construct synthetic dissipation fields or time series along a Lagrangian trajectory. These can then be used as input for chemical models of the ISM, that simulate a detailed chemical network but are unable to accurately model turbulence because of severe time-scale constraints. Another important prospect is the study of the properties of the structures of high dissipation as a function of the magnetic Prandtl number. This parameter study will allow us to investigate whether the approximate equality of different kinds of dissipation continues to hold if the viscosity is different from the resistivity. Let us note here that the ISM is characterized by a huge value of the magnetic Prandtl number, which brings forward serious theoretical and numerical challenges.

Incompressible models

Incompressible numerical models based on spectral methods are able to capture the intermittency properties of the velocity and magnetic fields. They have the interesting property to introduce very little numerical dissipation in comparison to grid-based codes. In the future, we plan to release the code ANK online. We believe that making this spectral code available will be a valuable contribution to the community. Ambipolar diffusion greatly reduces the magnitude of the Lorentz force at small scales: because of the increase of magnetic and kinetic energy at small scales, we conclude that this is due to an efficient alignment of the vectors j and b. This reduction shifts the peak of ambipolar diffusion heating at larger scales. The position of the peak is close to the thickness of the sheets of strong dissipation. This shift of the peak presumably sets the characteristic scale of the regions of molecule formation amplified by ambipolar diffusion heating, and its observational signature may be traced by molecular emission maps of specific species in diffuse gas: HCO + , CH + and SH + .

Qualitatively, dissipation seems to be concentrated on quasi two-dimensional sheet-like structures. Quantitatively, the high degree of intermittency is confirmed by the presence of a power-law tail in the pdf of the dissipation. This fact brings forward the importance of extreme heating events which are crucial for molecule formation in the ISM. Accordingly with the qualitative results, structure function exponents can be adequately fitted by a log-normal model which assumes Kolmogorov scaling in the inertial range and dissipative structures of dimension two. This rules out Iroshnikov-Kraichnan scaling based on turbulent dynamics dominated by Alfvén wave interaction in favor of dynamics based on interactions between eddies of neighboring size and vortex stretching. It also sets the geometry of the regions heated by turbulent dissipation in the ISM: they are sheets of dimension two. These results open interesting possibilities for qualitative and quantitative comparisons of numerical maps to observational maps of molecular emission.

Compressible models

Compressible numerical models based on Godunov methods are able to account for an additional, independent type of energy dissipation: energy dissipation in shocks. This leads to an increased dissipation of kinetic energy in comparison to magnetic energy. Indeed, in the majority of the simulations performed in the context of this thesis, magnetic energy was found to dominate over kinetic energy throughout the evolution.

In the future, we would like to investigate numerically a wider range of Mach numbers and initial conditions: rms sonic Mach numbers up to 10, or even more, are relevant to the ISM. As far as initial conditions are concerned, we would like to investigate the influence on the turbulent dynamics of the ratio of dilatational to rotational energy present in the initial velocity field. An important, although preliminary, finding is that only subsets of the regions of high dissipation follow the density peaks: dissipation takes place in areas of relatively low density.

For most of the evolution, the heating is dominated by the rotational component, and not the dilatational one. This hints to a reduced relative importance of shock waves. However, the relative contribution of numerical heating is large; we have made dedicated efforts to estimate the level of numerical dissipation and have found that it can be locally as high as the physical heating. We plan to use the methods presented in the appendix to compute the ratio of physical (i.e. corrected for numerical dissipation) viscous dilatational heating to viscous rotational heating. This should provide a much more reliable diagnostic for the relative importance of shock waves, and compressibility effects in general, in comparison with the vortical dynamics.

Another prospect is the use of the first method described in the appendix in the following way: after extracting ridges of high viscous dissipation from a simulation data cube, we could fit models of isothermal shocks in order to better evaluate the shock contribution to the total dissipation. We have to note here that in general, efficient and accurate extraction of shocks as well as shock fitting in a three-dimensional context presents serious algorithmic challenges. This line of work will address the important question: what fraction of total dissipation is due to shocks?

Comparison to observations

In chapter 6, we attempted a comparison of numerical results with observations using increments as diagnostics of various projected quantities, namely the increments of the centroid velocity, polarization direction and the Stokes parameters. We focused here on using results from incompressible simulations in this observational context. In the future, we would like to use more results from compressible simulations in our comparisons to observations, in order to investigate the importance of compressibility effects.

The contours of the diagnostic quantities were found to follow the structures in dissipation field integrated along the line of sight. However, because of projection effects, not all structures of high dissipation are captured by the diagnostic quantities. In general, the correlation is better for smaller integration lengths. This correlation is of particular importance to observers, because it provides evidence for a link between contours of observed diagnostic quantities and regions of elevated heating due to turbulent dissipation in the ISM.

The pdfs of the increments of the diagnostic quantities exhibit growing exponential tails at small scales, a signature of small-scale intermittency. At larger scales, the pdfs become increasingly Gaussian. Intermittency in the pdfs is reflected on the dependence of the structure function exponents on their order. These facts open up interesting possibilities for comparison between simulations and observations, where the same properties have been found. As an example of an interesting prospect, we mention the possibility of comparison with statistics of magnetic field diagnostics Q and U originating from observations of the Planck satellite.

As a future observational direction, we plan to use data cubes from our simulations to produce synthetic molecular spectra by simulation of the relevant radiative transfer processes. We can then compare these synthetic spectra to observational results. One important advantage here is that, contrary to most previous studies, we have access to the detailed spatial and velocity structure of the dissipation field which is crucial in the computation of radiative transfer in lines and in molecular line synthesis.

Coherent structures

The property of coherence is very important for the above results: indeed, if the Fourier phases of the velocity and magnetic fields are randomly mixed, the dissipation field, as well as the contours of the diagnostics, lose all of their structure and mutual correlation. This suggests that the worm-like structures in the projected maps are the signature of coherent structures of high dissipation, and these structures are indeed linked to intermittency. Fields with a power-law energy spectrum and random phases, which are candidates for providing input to chemical models, do not allow a good representation of the structures of high dissipation. Whether they are spatially coherent over a finite timescale that can be measured or spatially coherent as a propagating wave remains to be established.

As mentioned in [START_REF] Uritsky | Structures in magnetohydrodynamic turbulence: Detection and scaling[END_REF], the statistical analysis of the geometric and dynamical characteristics of high dissipation in the dissipative range admits an interpretation in terms of self-organized criticality, a fact which puts our study in a wider context. In this theory, a complex spatio-temporal phenomenon is interpreted in terms of avalanches; there is thus an incremental build-up of small-scale activity until the avalanche threshold is reached and a dissipative burst takes place.

In any case, our work shows, for the first time, that this level of spatio-temporal coherence corresponds well to that revealed by the statistical signatures of our kinematic and magnetic observables. It can be taken as a strong encouragement to further elaborate on the impact of these structures of intense dissipation on chemistry that proceeds swiftly, but also on reconnection of the magnetic field lines.

Our detailed knowledge of the geometry of the current sheets, which are identical to the structures of high Ohmic dissipation, opens up interesting possibilities for studies of magnetic reconnection. These studies could proceed by analyzing the topology of magnetic field lines in the vicinity of a current sheet. This investigation involves the challenge of extending our analysis of structures in the temporal domain, since magnetic reconnection is an inherently dynamical process. ). This method is impractical for any other dissipative structure than a planar shock, but in the present case it allows us to benchmark the other methods. In particular, the integral under both black lines is the same and is determined solely by the entrance conditions of the shock: the total irreversible heating through an isothermal shock does not depend on the value of ν. We would like our method to at least recover this integral, if not the exact width of the shock.

ν f or µ f /10 -3 n =128 n =256 n =512 n =1024
n =2048 cst. ν 0 = 10 -3 , u = 1 3.9 (58%) 2.3 (54%) 1.5 (41%) 1.2 (25%) 1.1 (13%) same, for u = 4 9.6 (-13%) 5.1 (-1%) 2.9 (1%) 1.9 (3%) 1.4 (4%) cst. µ 0 = 10 -3 , u = 1 5.8 (60%)

3.2 (57%) 2.0 (45%) 1.4 (35%) 1.1 (20%) same, for u = 4 26.7 (-6%) 13.4 (-2%) 7.0 (-3%) 3.8 (-2%) 2.1 (-2%) No visc. u = 1 3.3 (63%) 1.6 (63%) 0.8 (54%) 0.4 (60%) 0.2 (54%) No visc. u = 4 8.8 (-21%) 4.3 (-3%) 2.1 (-4%) 1.1 (-4%) 0.5 (-6%)

Table A.1: Adjusted viscosity of shocks with constant ν, constant µ or no viscosity at all. In the case without viscosity, it is the adjusted constant ν f which is given, as it fits the shock profile much better. The percentages in parenthesis give the difference between the total irreversible heating as recovered from the method 1 and the theoretical one.

Method 3

∂ t ( 1 2 ρu 2 + 1 2 B 2 ) + ∇.F 3 = p∇.uq where F 3 = F 1 + up. We evaluate q as in the previous two methods. The shape of the heating is positive but its shape is asymmetric and the total heating is grossly overestimated (by 180% for n = 128 and u = 1, not shown on figure A.1).

Conclusion

We adopt method 1 as the best compromise between methods 1, 2 and 3. With this method, we locally recover rather well the amount of heating, but we should be aware that the strongest bursts of dissipation are in reality about 4 times thinner than our simulations show. On the other hand, the slow convergence to the shock solution justifies our use of a moderate resolution associated to this recovering scheme: we would not gain much by running our simulations at twice the resolution. Finally, note that such a slow convergence could lure an unaware numericist into thinking his simulations are converged... Although thanks to this method we gained access to the total numerical dissipation, we could not find an accurate way to separate the numerical dissipation of magnetic fields from the numerical dissipation of kinetic energy. In order to compute corrected values for the viscous heating and the ohmic heating, we simply shared between each of them the total numerical heating in proportion to their relative physical values, namely: ε corr. . Whenever the numerical dissipation is negative (i.e: q < ε v + ε o ), we simply set ε corr.

Additional notes

Method 0 can be used in various circumstances to estimate the effective viscosity of the scheme. We computed various shocks (u = 1 or u = 4) at various resolutions and with a variety of modes of dissipation: a physical viscosity with constantν 0 or with constantµ 0 , or no viscosity at all. To each of these shocks, we fitted steadystate shocks model parametrised by either a constant dynamic viscosity µ f or a constant viscosity coefficient ν f . Table A.1 summarises the results of this short study.

The estimated ν f (or µ f ) gives an idea of the effective width of the shock compared to the physical width it should have.

The adjusted viscosity depends on the shock velocity u, with larger spread for higher velocities. Constant µ 0 viscosity, as should be used in principle for an isothermal gas, is much harder to resolve due to the density dependence of the viscous length. The resulting fits are usually much worse, because the heating shape is asymmetric, with a sharper edge at higher densities on the post-shock side. This is felt especially at large velocity u.

Shocks without physical viscosity are much better adjusted by a constant ν f than by a constant µ f : it is tempting to associate the constant size of the zones to an effective mean free path. However, the best fit value depends on the velocity of the shock: the length of the effective mean free path depends on the shock velocity... This non-linear effect of numerical dissipation is also seen in [START_REF] Lesaffre | Exact shearing box solutions of magnetohydrodynamic flows with resistivity, viscosity and cooling[END_REF] in the case of the ZEUS3D code.

Interestingly, in the case of shocks at our fiducial resolution (n = 128), there is no obvious gain between adding a physical viscous term or letting the scheme take care of viscosity. However, we postulate that it makes a difference in the subsonic regions of the flow, where the scheme dissipation is in principle very small.

Finally, due to the form of λ v = ν/ max(c, |u|), it is not equivalent to perform Mach 4 simulations by raising the r.m.s. velocity or by diminishing c, for a given vicous coefficient ν. We decreased c from 1 to 1/4, everything else being kept fixed. That way, the resolution power of our scheme remains constant with respect to the shocks.

1. 1

 1 Classification of interstellar cloud types (from Snow and McCall [2006]), see text for notation. 1.2 Energy densities in the local interstellar medium (from Draine [2011]). . . . . . . . . . . . . . 1.3 Comparison of molecular column densities for different types of shocks (J-type and C-type) at the same level of energy dissipation (units: cm -2 ). . . . . . . . . . . . . . . . . . . . . . . . 4.1 Parameters of the simulations. N : linear resolution, L: integral length scale at the peak of dissipation (pd) , λ: Taylor microscale (pd), l d : dissipative scale (pd), Re λ : Taylor microscale Reynolds number U λRe (pd), Re: kinetic Reynolds number, Re m : magnetic Reynolds number, Re a : AD Reynolds number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Results of the structure extraction algorithm for all runs and structures defined as connected sets of points having a value of the total dissipation two standard deviations above mean value. 4.3 Characteristics of various components of the ISM. Dimensions are recovered from our simulations by assuming L 2a where U a is the line-of-sight r.m.s. Alfvén velocity. The quantity γρ i is computed by assuming the ions are essentially C + ions with a number density n i = 10 -4 n H . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Parameters of compressible numerical simulations. N : linear resolution. M s : rms sonic Mach number. M a : rms Alfvénic Mach number. ν: viscosity coefficient (in units l 0 u r.m.s. with L box = 2πl 0 and u r.m.s. is the initial r.m.s. velocity). η: physical resistivity. Re: kinetic Reynolds number. Re m : magnetic Reynolds number. . . . . . . . . . . . . . . . . . . . . . . .

1. 1

 1 Comparison of heating rates and temperature in two different types of shock: J-shock and C-shock. The x axis is the distance from the first integration point. . . . . . . . . . . . . . . . 1.2 Abundances of OH, H 2 O, HCO + and CO inside a J-shock. . . . . . . . . . . . . . . . . . . . . 2.1 Examples of turbulent flows. Top left: Classical laboratory experiment of homogeneous

3. 1

 1 Absolute value of the error in the conservation of kinetic energy and kinetic helicity in incompressible ideal hydrodynamics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Absolute value of the error in the conservation of total energy, cross-helicity and magnetic helicity in incompressible ideal magnetohydrodynamics. . . . . . . . . . . . . . . . . . . . . . 3.3 Energy spectra from two identical hydrodynamic simulations, one well-resolved with k max l d = 2.0 (black line) and another under-resolved with k max l d = 0.5 (red dashed line). . . . . . . . . 3.4 The wave fan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 One-dimensional shock test results: DUMSES with Lax-Friedrichs solver. Blue: ρ. Green:

  (reproduced as figure 3.6 here). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Tilted shock test results: figure 1a from Kim et al. [1999]. Upper left: ρ. Upper middle: B y . Upper right: B z . Lower left: v x . Lower middle: v y . Lower right: v z . . . . . . . . . . . . . . . 3.7 Tilted shock test results: DUMSES with Lax-Friedrichs solver. Upper left: ρ. Upper middle: B y . Upper right: B z . Lower left: v x . Lower middle: v y . Lower right: v z . . . . . . . . . . . .

  ke a , at k = k * a = 2π/ * a , and the position of k = k a . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Power spectra of j × b for high resolution runs 10,12 and 14. The field becomes force-free at small scales when the strength of the AD is increased. . . . . . . . . . . . . . . . . . . . . . . 4.7 Color maps of a slice of the dissipation fields, run 10 (MHD -OT). Red: Ohmic dissipation, green: viscous dissipation. All snapshots are taken at the peak of dissipation. . . . . . . . . . 4.8 Same as Figure 4.7, Run 12, (AD -OT, Re a = 100) with blue: ambipolar diffusion heating. . 4.9 Same as Figure 4.7, Run 14, (AD -OT, Re a = 10) . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Compensated plot S ± z (r)/r for the OT runs 10,12 and 14. . . . . . . . . . . . . . . . . . . . . 4.11 ESS velocity field structure function exponents for ABC runs 9-11-13. . . . . . . . . . . . . . 4.12 ESS velocity field structure function exponents for OT runs 10-12-14. . . . . . . . . . . . . . 4.13 ESS magnetic field structure functions for ABC runs 9-11-13. . . . . . . . . . . . . . . . . . 4.14 ESS magnetic field structure functions exponents for OT runs 10-12-14. . . . . . . . . . . . . 4.15 Inertial range structures extracted from the dataset corresponding to the peak of dissipation

  instead of our own slope. . . . . . . . . . . . 4.18 Scaling exponents for structures extracted based on the Ohmic dissipation (red circles), comparison with results of UR10 (blue squares). Upper panel: Inertial range exponents for our Run 10 (MHD-OT) compared with the corresponding Run III of UR10. Lower panel: Dissipative range exponents for the same runs. Thick error bars are three-sigma error brackets. Thin red error bars estimate the systematics related to the choice of the boundaries of the inertial and dissipative ranges as well as the bin size for the definition of the pdfs: see text for details. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Pdfs of the density and column density. Green: OT initial condition, M s = 4. Red: ABC initial condition, M s = 1. Cyan: ABC initial condition, M s = 4. Left: Pdf of the density. Right: Pdf of the column density. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Dissipation color-maps for runs 1-4, integrated over 1/64th of the box. Red: Ohmic dissipation. Green: viscous dissipation. Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT). . . . 5.9 Color map of Ohmic, viscous dilatational and viscous rotational dissipation integrated over 1/64th of the box. Red: Ohmic dissipation. Green: viscous rotational dissipation ε vω . Blue: viscous dilatational dissipation ε vc . Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT). . . . 5.10 Color map of Ohmic, viscous and numerical dissipation. Red: Ohmic physical dissipation. Green: viscous physical dissipation ε vc . Blue: viscous numerical dissipation ε corr. v + ε corr. o -(ε v + ε o ). Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT). . . . . . . . . . . . . . . . . . . . 5.11 Pdf P(log 10 ε corr. ) of the total dissipation. Blue: M s = 0 (incompressible). Green: M s = 1. Red: M s = 4. Left: ABC initial condition. Right: OT initial condition. . . . . . . . . . . . . 6.1 CVIs based on the 12 CO(J = 2 -1) line observed in the Polaris Flare. Left: Normalized pdfs measured over variable lags expressed in units of 15 (upper right corner of each panel). Dashed red lines are Gaussian pdfs of same dispersion (given in km.s -1 at the bottom of each panel). Right: Spatial distribution of the azimuthally averaged CVI for a lag l = 4 × 15 . The color scale gives the values of the velocity increments in km.s -1 . The largest values (blue lanes) of CVI are associated to the tails of its pdf. From Hily-Blant et al. [2008a]. . . . . . . 6.2 Map of the polarization angle dispersion function S at 353 GHz with 1 deg resolution and for lag δ = 30 . The map is shown in log 10 scale over the range 0.1 deg < S < 70 deg. Only sky regions where the S/N on S is larger than 3 are shown (see Planck Collaboration Int. XIX, 2015, in press, A & A) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 For an integration depth of 1/64th of the box. Top left: Dissipation field. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. Top right: CVI ε . Bottom left: δQ. Bottom right: δψ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 For an integration depth of 1/8th of the box. Top left: Dissipation field. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. Top right: CVI ε . Bottom left: δQ. Bottom right: δψ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 For an integration depth over the whole box. Top left: Dissipation field. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. Top right: CVI ε . Bottom left: δQ. Bottom right: δψ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Color maps of the dissipation field (Red: Ohmic, Green: viscous, Blue: ambipolar) overlaid with contours of the diagnostics CVI ε (white), δψ (blue), δU (red), δQ (green), for a lag l = 3 pxiels.Top left: Projection over 1/64th of the box. Top right: Projection over 1/8th of the box. Bottom left : Projection over the whole box. . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Efficiency of CVI ε as a function of the contour level. Red line: projection over the whole box. Green dotted line: projection over 1/8th of the box. Blue dashed line: projection oven 1/64th of the box. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Efficiency of CVI ε as a function of the lag l, for a contour level of 0.5. Red line: projection over the whole box. Green dotted line: projection over 1/8th of the box. Blue dashed line: projection over 1/64th of the box. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Efficiency of CVI ε for different types of dissipation and a lag l = 1. Left: Projection over 1/64th of the box. Right: Projection over the whole box. Blue line: total dissipation. Green dotted line: viscous dissipation. Red dashed line: Ohmic dissipation. Green dash-dotted line: ambipolar diffusion dissipation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Compared efficiency of different tracers of Ohmic dissipation for two integration depths. Left: Integration over 1/64th of the box. Right: Integration over the whole box. Blue line: CVI ε . Green dotted line: δU . Red dashed line: δQ. Greed dash-dotted line: δψ. . . . . . . . . . . . 6.11 Pdfs of the increments (normalised by their standard deviation) of diagnostic quantities integrated over the whole box, for various lags l. For clarity, they are artificially shifted vertically (their maximum value is set to 1/l). Top left: P l (CVI ε ). Top right: P l (δU ). Bottom left: P l (δQ). Bottom right: P l (δψ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 Structure functions S CVε p (l) as a function of lag (left), and computed using ESS, with S CVε 3 (l) taken as a reference structure function (right). Vertical lines in the left panel denote the extent of the intertial range in the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Intermittency exponents computed using ESS Left: OT incompressible runs. Right: ABC incompressible runs. Bottom: compressible runs (at 512 3 resolution). The black thin line is p/3, commonly referred to as K41. The diamonds correspond to the generalised She-Lévêque model with g = 3 and C = 1. The only exponents observed to date are those of Hily-Blant et al. (2008), displayed here as blue triangles. They happen to be very close to g = 3 and C = 2. 6.14 Color-map of the energy dissipation field in the case of weak phase mixing. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. A.D. dissipation dominates the rest (see text). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.15 Color-map of the viscous dissipation field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.16 Color-map of the ohmic dissipation field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.17 Color-map of the weighted centroid velocity increment CVI ε (x, l) for a lag l = 3, projected over the whole box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.18 Color-map of the weighted increment δψ(x, l) for a lag l = 3, projected over the whole box . . 6.19 Pdf of the weighted centroid velocity increment CVI ε (x, l) projected over the whole box for various lags l. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.20 Intermittency exponents for all variables in the case of weak phase mixing (dashed curves) and without phase mixing (full curves). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1

Figure 1 . 1 :

 11 Figure 1.1: Comparison of heating rates and temperature in two different types of shock: J-shock and C-shock. The x axis is the distance from the first integration point. As an example of how shocks can lead to molecule formation, consider a 10 km s -1 J-type steady state shock, figure 1.1. In such a shock, viscous friction raises the gas temperature to values as high as 4300 K, which is well above the 2980 K activation barrier of the reaction H 2 + O → OH + H. (1.3) The presence of OH allows the formation of other molecules through the reaction OH + H 2 → H 2 O + H (1.4) which has an activation barrier of 1490 K and C + + H 2 O → HCO + + H (1.5)

Figure 1 . 2 :

 12 Figure 1.2: Abundances of OH, H 2 O, HCO + and CO inside a J-shock.

Figure 2 . 1 :

 21 Figure 2.1: Examples of turbulent flows. Top left: Classical laboratory experiment of homogeneous isotropic turbulence generated behind a grid[START_REF] Corke | A new view on origin, role and manipulation of large scales in turbulent boundary layers[END_REF]. Top right: Phytoplankton bloom reveals turbulent ocean currents, an example of geophysical turbulence[START_REF] Signorini | Environmental factors controlling the barents sea spring-summer phytoplankton blooms[END_REF]. Bottom left: Turbulence in the interstellar medium as revealed by the gradient of linear polarization[START_REF] Gaensler | Low-mach-number turbulence in interstellar gas revealed by radio polarization gradients[END_REF]. Bottom right: Current density magnitude in a high resolution numerical experiment of homogeneous isotropic magnetohydrodynamic turbulence [G.Momferratos, unpublished].

Figure 2

 2 Figure 2.2: Streamwise velocity component recorded by a hot-wire anemometer in a classical wind tunnel experiment of homogeneous isotropic turbulence behind a grid[START_REF] Sreenivasan | Fractals and multifractals in fluid turbulence[END_REF].

Figure 2

 2 Figure 2.3: Schematic representation of the energy cascade process (from Ecke [2005]).

Figure 2

 2 Figure 2.4: Viscous dissipation rate along a fixed line from a 1024 3 incompressible magnetohydrodynamics simulation with initial condition based on the ABC flow. The snapshot shown is at the temporal peak of total (ohmic + viscous) dissipation [G. Momferratos, unpublished].

Figure 2

 2 Figure 2.5: A realization of one-dimensional Brownian motion. "Wiener process zoom". Licensed under Creative Commons Attribution-Share Alike 3.0 via Wikimedia Commons.

Figure 2

 2 Figure 2.6: Pdf of the longitudinal velocity increment from a high resolution incompressible hydrodynamical simulation for various separations r n = 2 n ∆x, ∆x being the distance between grid points in the simulations. (from Ishihara et al. [2009]).

Figure 2 0 e

 20 Figure 2.7: Comparison of experimental data with the predictions of the structure function scaling exponents of various intermittency models (from Anselmet et al. [2001]).

Figure 3 . 1 :

 31 Figure 3.1: Absolute value of the error in the conservation of kinetic energy and kinetic helicity in incompressible ideal hydrodynamics.

Figure 3

 3 Figure 3.3: Energy spectra from two identical hydrodynamic simulations, one well-resolved with k max l d = 2.0 (black line) and another under-resolved with k max l d = 0.5 (red dashed line).

Figure 3 . 4 :

 34 Figure 3.4: The wave fan.

Figure 3

 3 Figure 3.5: One-dimensional shock test results: DUMSES with Lax-Friedrichs solver. Blue: ρ. Green: B y .Red : v x . Cyan: v y . To be compared with[START_REF] Kim | A multidimensional code for isothermal magnetohydrodynamic flows in astrophysics[END_REF], figure1a(reproduced as figure3.6 here).

Figure 3

 3 Figure 3.6: Tilted shock test results: figure 1a from Kim et al. [1999]. Upper left: ρ. Upper middle: B y . Upper right: B z . Lower left: v x . Lower middle: v y . Lower right: v z .

Figure 3

 3 Figure 3.7: Tilted shock test results: DUMSES with Lax-Friedrichs solver. Upper left: ρ. Upper middle: B y . Upper right: B z . Lower left: v x . Lower middle: v y . Lower right: v z .

(

  u x , u y , u z ) = (-2 sin y, 2 sin x, 0) (b x , b y , b z ) = (-2 sin(2y) + sin z, 2 sin(x) + sin z, sin x + sin y)

Figure 4 . 1 :

 41 Figure 4.1: Compensated kinetic (top) and magnetic (bottom) energy spectra for OT runs 10, 12 and 14 at the temporal peak of dissipation. The assumed limits of the inertial and dissipation ranges are shown in dashed vertical lines. The vertical blue and red dashed-dotted lines with square symbols correspond to the expected AD critical wavenumber k a (see equation (4.5)) in runs 12 and 14, respectively.

Figure 4 . 2 :

 42 Figure 4.2: Time evolution of volume integrated dissipation rates for the OT runs 10, 12 and 14. Solid lines show the Ohmic plus viscous dissipation which we use to define peak dissipation.

Figure 4

 4 Figure 4.3: Log-normal core and power-law tail fit for the pdf of the total dissipation (Run 12 AD -OT).The core follows a log-normal distribution with mean µ l -4.27 and standard deviation σ l 1.03 while the tail follows a power-law with exponent -τ -2.61. Vertical lines show the mean value (red) and thresholds located at 1 (black), 2 (green) and 3 (blue) standard deviations above the mean value.

Figure 4

 4 Figure 4.4: Cumulative probability density function of the total dissipation for run 12 (AD -OT).

Figure 4 . 7 :

 47 Figure 4.7: Color maps of a slice of the dissipation fields, run 10 (MHD -OT). Red: Ohmic dissipation, green: viscous dissipation. All snapshots are taken at the peak of dissipation.

Figure 4

 4 Figure 4.8: Same as Figure 4.7, Run 12, (AD -OT, Re a = 100) with blue: ambipolar diffusion heating.

Figure 4

 4 Figure 4.9: Same as Figure 4.7, Run 14, (AD -OT, Re a = 10)

Figure 4 .

 4 Figure 4.10: Compensated plot S ± z (r)/r for the OT runs 10,12 and 14.

Figure 4 .

 4 Figure 4.11: ESS velocity field structure function exponents for ABC runs 9-11-13.

Figure 4 .

 4 Figure 4.12: ESS velocity field structure function exponents for OT runs 10-12-14.

Figure 4 . 13 :

 413 Figure 4.13: ESS magnetic field structure functions for ABC runs 9-11-13.

Figure 4 .

 4 Figure 4.14: ESS magnetic field structure functions exponents for OT runs 10-12-14.

Figure 4 .

 4 Figure 4.15: Inertial range structures extracted from the dataset corresponding to the peak of dissipation of Run 12 (AD -OT). They are defined as connected sets of points having values of the total dissipation three standard deviations above the mean value. Each little sphere has a 2-pixels diameter, ie: about the size of the viscous (or equivalently Ohmic) dissipation length.

Figure 4 .

 4 Figure 4.16: One of the largest structures extracted from run 12

Figure 4 .

 4 Figure 4.18: Scaling exponents for structures extracted based on the Ohmic dissipation (red circles), comparison with results of UR10 (blue squares). Upper panel: Inertial range exponents for our Run 10 (MHD-OT) compared with the corresponding Run III of UR10. Lower panel: Dissipative range exponents for the same runs. Thick error bars are three-sigma error brackets. Thin red error bars estimate the systematics related to the choice of the boundaries of the inertial and dissipative ranges as well as the bin size for the definition of the pdfs: see text for details.

Figure 4 . 13 Figure 4

 4134 Figure 4.20: Comparison of scaling exponents with three-sigma error bars between pure MHD (red circles) and AD MHD (green squares Re a = 100 and blue triangles Re a = 10) -ABC Runs 9,11,13

Figure 4 .

 4 Figure 4.22: Scatter plots of the ratio ε o /ε t versus the ratio ε a /ε t for AD-OT runs 12 and 14.

Figure 5 . 1 :

 51 Figure 5.1: Time evolution of the ratio of kinetic to magnetic energy E k /E m . Time is in units of l 0 /u r.m.s. . Blue: OT initial condition M s = 0 (incompressible). Green: OT initial condition M s = 1. Red: OT initial condition M s = 4. Dashed lines: same colors for ABC initial condition.

Figure 5 .

 5 Figure 5.1 presents the time evolution of the ratio of kinetic to magnetic energy E k /E m for ABC and OT initial conditions and three different Mach numbers: 0 (incompressible runs from the previous chapter), 1 and 4. Magnetic energy dominates over kinetic energy all along the time evolution, except in the case of the ABC M s = 1 and 4 runs, where kinetic energy dominates at early times, up to dissipation peak (see

Figure 5

 5 Figure 5.2: Time evolution of total physical and numerical heating rates. Blue: M s = 0 (incompressible). Green: M s = 1. Red: M s = 4. Solid lines: total heating rate -Ėtot (see text). Dashed lines: physical heating rate ε phys . Diamonds: total heating rate ε corr as recovered by method 1 of the appendix (see text). Left: ABC initial condition. Right: OT initial condition.

Figure 5 .

 5 Figure 5.2 presents the time evolution of total physical and numerical heating. The total energy density E tot is computed from

Figure 5 . 3 :

 53 Figure 5.3: Time evolution of the ratio of viscous to total physical dissipation rates ε v /( ε v + ε o ). Blue: M s = 0 (incompressible). Green: M s = 1. Red: M s = 4. Solid lines: OT initial conditions. Dashed lines: ABC initial conditions.

Figure 5 . 4 :

 54 Figure 5.4: Time evolution of the ratio of dilatational to total viscous dissipation rates (vortical plus dilatational) ε vc /( ε vc + ε vω ). Green: M s = 1. Red: M s = 4. Solid lines: OT initial conditions. Dashed lines: ABC initial conditions.

Figure 5 . 5 :

 55 Figure 5.5: Density cuts for runs 1-4. Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT).

Figure 5 . 6 :

 56 Figure 5.6: Column density (i.e. density field integrated over the whole cube) maps for runs 1-4. Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT).

Figure 5

 5 Figure 5.7: Pdfs of the density and column density. Green: OT initial condition, M s = 4. Red: ABC initial condition, M s = 1. Cyan: ABC initial condition, M s = 4. Left: Pdf of the density. Right: Pdf of the column density.

Figure 5 . 8 :

 58 Figure 5.8: Dissipation color-maps for runs 1-4, integrated over 1/64th of the box. Red: Ohmic dissipation. Green: viscous dissipation. Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT).

Figure 5 . 9 :

 59 Figure 5.9: Color map of Ohmic, viscous dilatational and viscous rotational dissipation integrated over 1/64th of the box. Red: Ohmic dissipation. Green: viscous rotational dissipation ε vω . Blue: viscous dilatational dissipation ε vc . Top left: Run 1 (M s = 1, ABC). Top right: Run 2 (M s = 1, OT). Bottom left: Run 3 (M s = 4, ABC). Bottom right: Run 4 (M s = 4, OT).

Figure 5 .

 5 Figure 5.10: Color map of Ohmic, viscous and numerical dissipation. Red: Ohmic physical dissipation. Green: viscous physical dissipation ε vc . Blue: viscous numerical dissipation ε corr.v

Figure 5 .

 5 Figure 5.11: Pdf P(log 10 ε corr. ) of the total dissipation. Blue: M s = 0 (incompressible). Green: M s = 1. Red: M s = 4. Left: ABC initial condition. Right: OT initial condition.

Figure 6 . 1 :

 61 Figure 6.1: CVIs based on the 12 CO(J = 2 -1) line observed in the Polaris Flare. Left: Normalized pdfs measured over variable lags expressed in units of 15 (upper right corner of each panel). Dashed red lines are Gaussian pdfs of same dispersion (given in km.s -1 at the bottom of each panel). Right: Spatial distribution of the azimuthally averaged CVI for a lag l = 4 × 15 . The color scale gives the values of the velocity increments in km.s -1 . The largest values (blue lanes) of CVI are associated to the tails of its pdf. From Hily-Blant et al. [2008a].

Figure 6

 6 Figure 6.2: Map of the polarization angle dispersion function S at 353 GHz with 1 deg resolution and for lag δ = 30 . The map is shown in log 10 scale over the range 0.1 deg < S < 70 deg. Only sky regions where the S/N on S is larger than 3 are shown (see Planck Collaboration Int. XIX, 2015, in press, A & A)

Figure 6 . 3 :

 63 Figure 6.3: For an integration depth of 1/64th of the box. Top left: Dissipation field. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. Top right: CVI ε . Bottom left: δQ. Bottom right: δψ.

Figure 6 . 4 :

 64 Figure 6.4: For an integration depth of 1/8th of the box. Top left: Dissipation field. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. Top right: CVI ε . Bottom left: δQ. Bottom right: δψ.

Figure 6 . 5 :

 65 Figure 6.5: For an integration depth over the whole box. Top left: Dissipation field. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. Top right: CVI ε . Bottom left: δQ. Bottom right: δψ.

Figure 6 . 6 :

 66 Figure 6.6: Color maps of the dissipation field (Red: Ohmic, Green: viscous, Blue: ambipolar) overlaid with contours of the diagnostics CVI ε (white), δψ (blue), δU (red), δQ (green), for a lag l = 3 pxiels.Top left: Projection over 1/64th of the box. Top right: Projection over 1/8th of the box. Bottom left : Projection over the whole box.

Figure 6 . 7 :

 67 Figure 6.7: Efficiency of CVI ε as a function of the contour level. Red line: projection over the whole box. Green dotted line: projection over 1/8th of the box. Blue dashed line: projection oven 1/64th of the box.

Figure 6

 6 Figure 6.8: Efficiency of CVI ε as a function of the lag l, for a contour level of 0.5. Red line: projection over the whole box. Green dotted line: projection over 1/8th of the box. Blue dashed line: projection over 1/64th of the box.

Figure 6 . 9 :

 69 Figure 6.9: Efficiency of CVI ε for different types of dissipation and a lag l = 1. Left: Projection over 1/64th of the box. Right: Projection over the whole box. Blue line: total dissipation. Green dotted line: viscous dissipation. Red dashed line: Ohmic dissipation. Green dash-dotted line: ambipolar diffusion dissipation.

Figure 6 .

 6 Figure 6.10: Compared efficiency of different tracers of Ohmic dissipation for two integration depths. Left: Integration over 1/64th of the box. Right: Integration over the whole box. Blue line: CVI ε . Green dotted line: δU . Red dashed line: δQ. Greed dash-dotted line: δψ.

Figure 6 .

 6 Figure 6.11: Pdfs of the increments (normalised by their standard deviation) of diagnostic quantities integrated over the whole box, for various lags l. For clarity, they are artificially shifted vertically (their maximum value is set to 1/l). Top left: P l (CVI ε ). Top right: P l (δU ). Bottom left: P l (δQ). Bottom right: P l (δψ).

Figure 6 .

 6 Figure 6.12 shows the structure functions S CVε p (l) computed using ESS, with S CVε 3 (l) taken as a reference structure function. The region of power-law dependence is significantly extended in comparison with the computation without ESS (shown in the left panel).

Figure 6 .

 6 Figure 6.12: Structure functions S CVε p (l) as a function of lag (left), and computed using ESS, with S CVε 3 (l) taken as a reference structure function (right). Vertical lines in the left panel denote the extent of the intertial range in the simulation.

Figure 6 . 13 :

 613 Figure 6.13: Intermittency exponents computed using ESS Left: OT incompressible runs. Right: ABC incompressible runs. Bottom: compressible runs (at 512 3 resolution). The black thin line is p/3, commonly referred to as K41. The diamonds correspond to the generalised She-Lévêque model with g = 3 and C = 1. The only exponents observed to date are those of Hily-Blant et al. (2008), displayed here as blue triangles. They happen to be very close to g = 3 and C = 2.

Figure 6 .

 6 Figure 6.14: Color-map of the energy dissipation field in the case of weak phase mixing. Red: Ohmic dissipation. Green: viscous dissipation. Blue: ambipolar diffusion dissipation. A.D. dissipation dominates the rest (see text).

Figure 6 .

 6 Figure 6.15: Color-map of the viscous dissipation field.

Figure 6 .

 6 Figure 6.16: Color-map of the ohmic dissipation field.

Figure 6 .

 6 Figure6.17: Color-map of the weighted centroid velocity increment CVI ε (x, l) for a lag l = 3, projected over the whole box

Figure 6 .

 6 Figure6.20: Intermittency exponents for all variables in the case of weak phase mixing (dashed curves) and without phase mixing (full curves).

Figure

  Figure A.1:

v=

  qε v /(ε v + ε o ) and conversely for ε corr.

  o

v=

  ε v and ε corr. o = ε o .

Table 1

 1

	-3 )

.2: Energy densities in the local interstellar medium (from

  11)which then produces CH + by dissociative recombination of CH + 3 . The comparison of the chemistry between the two types of shock at the same velocity shows (table1.3) that there is a clear composition difference depending on the nature of the dissipation, for a given dissipated energy: 10 13 3.36 × 10 13 3.48 × 10 13 7.77 × 10 10 4.36 × 10 13 1.69 × 10 11 4.80 × 10 12 N (C): 9.29 × 10 13 2.67 × 10 13 6.35 × 10 12 2.88 × 10 12 6.71 × 10 11 5.50 × 10 11 9.49 × 10 13 Table 1.3: Comparison of molecular column densities for different types of shocks (J-type and C-type) at the same level of energy dissipation (units: cm -2 ).

	species:	CH	CO	OH	CH +	H 2 O	HCO +	CH + 3
	N (J):	1.34 ×

Table 4 .

 4

	#	N	L	λ	l d	Re λ Re = Re m Re a Initial condition
	1	128 2.60 1.30 0.0607 189	219	-	ABC
	2	128 2.59 1.31 0.0605 210	219	-	OT
	3	128 2.70 1.44 0.0390 217	219	100	ABC
	4	128 2.72 1.44 0.0381 228	219	100	OT
	5	256 2.46 0.94 0.0375 353	551	-	ABC
	6	256 2.46 0.91 0.0389 377	551	-	OT
	7	256 2.62 1.08 0.0185 412	551	100	ABC
	8	256 2.66 1.09 0.0196 444	551	100	OT
	9	512 2.28 0.67 0.0237 591	1374	-	ABC
	10 512 2.12 0.58 0.0245 604	1374	-	OT
	11 512 2.49 0.83 0.0091 756	1374	100	ABC
	12 512 2.36 0.75 0.0095 750	1374	100	OT
	13 512 1.84 0.58 0.0074 640	1374	10	ABC
	14 512 2.80 1.00 0.0084 927	1374	10	OT

1: Parameters of the simulations. N : linear resolution, L: integral length scale at the peak of dissipation (pd) , λ: Taylor microscale (pd), l d : dissipative scale (pd), Re λ : Taylor microscale Reynolds number U λRe (pd), Re: kinetic Reynolds number, Re m : magnetic Reynolds number, Re a : AD Reynolds number.

  Figure 4.17: Scaling relations P i ∝ L D P

	All data
	Inertial data
	Inertial fit
	Uritsky slope
	Dissip. data
	Dissip. fit
	Uritsky slope
	i

Table 4 .

 4 2: Results of the structure extraction algorithm for all runs and structures defined as connected sets of points having a value of the total dissipation two standard deviations above mean value.

	#	N	Re λ Re a # of structures % of volume % of total dissipation
	1	128 189	-	87	3.90	23.85
	2	128 210	-	54	3.81	27.29
	3	128 217 100	72	3.78	23.65
	4	128 228 100	34	3.90	26.91
	5	256 353	-	215	3.26	27.32
	6	256 377	-	233	3.26	26.80
	7	256 412 100	144	2.86	26.40
	8	256 444 100	153	2.98	26.53
	9	512 591	-	790	2.70	31.30
	10 512 604	-	1166	2.77	30.33
	11 512 756 100	375	2.17	29.21
	12 512 750 100	418	2.57	29.48
	13 512 640	10	1167	2.94	28.01
	14 512 927	10	378	1.75	22.00

Table 4 .

 4 3: Characteristics of various components of the ISM. Dimensions are recovered from our simulations a where U a is the line-of-sight r.m.s. Alfvén velocity. The quantity γρ i is computed by assuming the ions are essentially C + ions with a number density n i = 10 -4 n H .

	by assuming L	2.5, * a	0.6 and u 0	√	3.U CNM	molecular clouds low-mass dense cores
	Density		n H (cm -3 )	30	200	10 4
	Length scale r.m.s. velocity Alfvén velocity	L.l 0 (pc) U.u 0 / √ 3 (km/s) U a = u 0 / √ 3 (km/s)	10 3.5 3.4	3 1 2	0.1 0.1 1
	AD Reynolds number AD heating length AD dynamical length	Re -1 a = 1/γρ i /(t 0 ) * a .l 0 (pc) a .l 0 = 2πl 0 Re -1 a .U 2 a /U 2 (pc)	1.2 10 -2 2.4 0.28	3.6 10 -3 0.72 0.10	1.1 10 -3 0.024 0.026

  Table5.1: Parameters of compressible numerical simulations. N : linear resolution. M s : rms sonic Mach number. M a : rms Alfvénic Mach number. ν: viscosity coefficient (in units l 0 u r.m.s. with L box = 2πl 0 and u r.m.s. is the initial r.m.s. velocity). η: physical resistivity. Re: kinetic Reynolds number. Re m : magnetic Reynolds number.

					M s M a Initial condition	ν = η	Re = Re m
			1 2 3 4 5 6 7 8	512 512 512 512 1020 1020 1020 1020	1 1 4 4 1 1 4 4	1 1 1 1 1 1 1 1	ABC OT ABC OT ABC OT ABC OT	1.0 × 10 -3 1.0 × 10 -3 1.0 × 10 -3 1.0 × 10 -3 7.0 × 10 -4 7.0 × 10 -4 7.0 × 10 -4 7.0 × 10 -4	6.2 × 10 3 6.2 × 10 3 6.2 × 10 3 6.2 × 10 3 8.9 × 10 3 8.9 × 10 3 8.9 × 10 3 8.9 × 10 3
	5.2 Physical dissipation in DUMSES
	In the numerical scheme implemented in DUMSES, physical viscosity and resistivity are included with
	two half-steps giving the intermediate state vectors U	n-1/2 i,j,k	and U	n+1/2 i,j,k	around the inviscid, non-resistive
	Godunov step giving the state vector U n i,j,k . The diffusive fluxes F d i,j,k which are used for the computation
	of U	n-1/2 i,j,k	n+1/2 i,j,k and U

k B 1 2 µ in (u iu n ) 2(1.9)

From the ancient greek word " τύρβη".

In the following, we will use the symbol ∝ to express proportionality whithin a constant of order unity.

The conservation of kinetic helicity and the law of conservation of circulation require, in the most general case, that ∇p/ρ can be written as a gradient of a function Q(p, ρ).

In the following we will use the velocity field as an example when we refer to general properties of turbulence.

In this chapter, we will follow the naming convention of moments in homogeneous turbulence where the indices of the components that correspond to the position x + r follow a coma.

We will use longitudinal structure functions for definiteness, although there is no distinction between longitudinal and transverse structure functions in the context of this chapter.

The code is named after Andrei Nikolaievich Kolmogorov, the Russian mathematician who was also one of the greatest contributors to the theory of turbulence.

http://www.fftw.org

This chapter was published on Monthly Notices of the Royal Astronomical Society[Momferratos et al.,

2014]

For example, with = la, we find that the AD Reynolds number of a C-shock is either R AD = 1 or R AD = M 2 a depending on whether we compute it with the post-shock or the pre-shock magnetic field strength.

http://www.lra.ens.fr/ ~giorgos/ank

where P (ε l ) is the probability density function of ε l and

(2.184) Thus, Q p (ε l ) is a weighted probability density function of which the intensity ε (p) l

is the expectation value. If P (ε l ) has a single maximum and an exponentially decaying tail, ε (p) l is a monotonically increasing sequence of p. It can be shown that if ε (∞) l < ∞, the variance of Q p (ε l ) tends to zero as p tends to infinity, thus for large p Q p (ε l ) is increasingly concentrated around ε (p) l and captures structures of this specific dissipation intensity. In addition to this, lim p→∞ ε

. It is known from numerical experiments [START_REF] Siggia | Numerical study of small-scale intermittency in three-dimensional turbulence[END_REF][START_REF] Kerr | Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF]] that dissipation of high intensity displays a certain degree of coherence in space, typically organizing itself in filamentary structures [Moisy and Jiménez, 2004]. The motivation behind the original formulation of the log-Poisson model is that these coherent structures can have scaling properties which are captured by the hierarchy ε

The first structure of the hierarchy ε 0 l is scale independent by definition, as it is equal to the mean energy dissipation rate ε used in the context of K41 to predict normal, non-anomalous scaling. Here, anomalous scaling is introduced by the divergent scaling dependence of ε (∞) l as l → 0 due to the presence of filamentary structures of high dissipation. The last structure in the hierarchy can be estimated dimensionally as

where

is the kinetic energy available to the cascade and t l is the cascade timescale. She and Leveque [1994b] choose to introduce anomalous scaling through E (∞) l , and thus t l scales just as in K41:

(2.186)

This assumes that the nonlinear processes that determine the timescale of the cascade are homogeneous in space, irrespective of the presence of structures of high dissipation. If we further assume that kinetic energy dissipated by the most intense structures E (∞) l is independent of scale and equal to the large-scale kinetic energy 1 2 U 2 , we obtain the scaling

where the K41 estimate ε = U 3 /L i has been used. This is equivalent to assuming an anomalous l 2/3 scaling law for the energy flux. Equations (2.183) and (2.187) imply that as p → ∞, τ p+1τ p → 2/3, or

(2.188)

By equation (2.161), C 0 can be interpreted as the co-dimension of the most intense dissipative structures.

Assuming that these structures are filaments, we have C 0 = 2 and τ p = 2/3p + 2. The dissipation intensities ε (p) l for p < ∞ characterize less intense dissipative structures, and the scaling of these intensities with l can provide a measure of how singular these structures are. The less intense structures are assumed to be the results of the same dynamics which form the most intense structures ε (∞) l , but fail to accomplish this high degree of coherence because of background disorder. In physical space, structures of higher intensity are likely to be surrounded by structures of lower intensity. Similarly in time, structures of lower intensity tend to come first. These facts suggest that there may be a relation between the intensities ε

. The relation proposed by She and Leveque [1994b] as a consequence of some hidden (statistical) symmetries of the Navier-Stokes equations reads

Chapter 6

Comparison with observations

In this chapter we present synthetic maps of diagnostic quantities that allow the comparison of numerical simulations with observational results, with emphasis on the dissipation field. In the context of incompressible ambipolar diffusion magnetohydrodynamics, we attempt to establish a link between structures of high dissipation and increments of projected diagnostic quantities which are accessible to observation. With the exception of subsection 6.4 (on intermittency exponents), all numerical results discussed in this section originate from a decaying simulation starting with an Orszag-Tang-based initial condition (simulation 8 in chapter 4). The Taylor microscale Reynolds number is Re λ = 750, the magnetic Prandtl number is equal to unity and the ambipolar diffusion Reynolds number Re AD = 100.

Introduction

In numerical simulations, one has access to all flow variables in any space-time point. In astrophysical observations, the projections severely hamper the quantitative information extracted from the observables. First, all the observables are obtained in one direction in the sky (x, y). Therefore any observed quantity is indeed a line-of-sight (LoS) integration of the corresponding variable. The gas velocity is measured using the Doppler effect of spectral lines of known rest frequency emitted by the gas, and the Doppler shift is a direct measure of the gas velocity projected on the LoS. The integration of the line emission across the medium that has some velocity dispersion σ provides a line profile T (v) of half-power width ∆v = 2 √ ln 2σ, neglecting radiative transfer. T (v) traces the emission of all gas cells along the LoS which have a given LoS-projected velocity u z = v. Any spectral line observed with high enough spectral resolution can be characterized by its moments in velocity space. In the following, we use the line centroid (or first moment) defined as

Given these projection constraints, one important observational diagnostic of the local structure of the interstellar velocity field is the centroid velocity increment (CVI), which is defined in terms of the centroid velocity CV. Below, we investigate the extent to which this quantity can be taken as a surrogate to viscous dissipation. In a numerical context the (synthetic) centroid velocity is defined as the integral of a velocity component along the corresponding line of sight:

where L d is an integration depth along the z direction. Then, the centroid velocity increment is defined, for each point on the plane of the sky x = (x, y), as the average difference between its value on x and the values on a circle C l of radius l centered at x:

Numerical dissipation in Godunov methods

We report here on the methods we used to recover the amount of numerical dissipation present in our compressible simulations.

In our compressible simulations, we have adopted twice the resolution of the incompressible ones (1020 vs 512), for the same dissipation coefficients (viscosity and resistivity). Indeed, there is a common belief that grid based methods need twice as many elements to obtain a resolving power equivalent to Fourier elements. This led us to run 1020 3 pixels simulations, which were about the most we could afford given our CPU time allowance.

In order to check whether our simulations are resolved with respect to dissipation, we run a 1D planar isothermal non-magnetized shock simulation with various resolutions. We adopt a viscosity ν = 10 -3 cL for a box length of unit L = 1, speed of sound c=1 and initial homogeneous left-directed inflow velocity u = 1 and homogeneous density ρ = 1. Boundary conditions are inflow on the right boundary and reflective wall on the left boundary. As the gas flows from right to left and hits the left boundary, a shock is produced at the wall and detaches from it as time proceeds. We look at the resulting shock at time t = 0.1L/c. In the shock frame, the resulting Mach number is about 1.61 for a γ = 5/3 polytropic gas.

Shocks have a viscous spread on the order of λ v = ν/ max(c, |u|). Our 3D simulation of decaying turbulence with 1020 3 pixels have box length L box = 2π and viscosity ν = 0.7 × 10 -3 : the corresponding resolution for the shock experiment is about n = 113 zones if we want the same number of zones in a viscous length λ v . We hence run 1D shocks simulations where we vary the resolution in powers of two from n = 128 up to n = 2048. In each of these runs, we compare the physical viscous heating term ε v = ν 4 3 ρ(∂ x u) 2 to the actual solution (integrated numerically from the 1D steady-state equations).

Figure A.1 shows the irreversible heating in non-dimensional units at a close-up of the shock front. It demonstrates that the resolution convergence for the heating rate is very slow and fully obtained only for n = 2048 (see red, blue and green solid curves vs. solid black curve), which corresponds to about 20 zones in a viscous length λ v . This would require a 18400 3 pixels box for the 3D turbulence simulation: way out of our capabilities. The situation closest to our 3D 1020 3 simulations is the red curve (n = 128): the viscous heating is largely underestimated and spread out.

Hence, we need to work a little bit if we want to use the results of our turbulence simulations to study dissipation. There are two sources of dissipation in our simulations: viscous (and Ohmic) dissipation due to the physical term we have introduced in Dumses, and numerical dissipation intrinsic to the scheme. The whole issue is to recover what amount of dissipation is produced by the scheme. We imagine here several methods to achieve this, in particular by considering variants of the energy conservation equation.

Method 0

We fit viscous isothermal shock models to the velocity and density profiles and we find a best fit value for the viscosity of ν = 3.3 × 10 -3 cL. Then, we compute the viscous heating in the best fit shock model (black

Method 1

Consider the evolution equation of kinetic and magnetic energy:

∂ t ( 1 2 ρu 2 + 1 2 B 2 ) + ∇.F 1 = -u.∇(p)q where q is the total irreversible heating and where the flux F 1 reads:

We compute the left hand side of equation (A) along a time step of the simulation, using the flux estimates of each dissipative half-step for the two dissipative contributions to F 1 and using a Lax-Friedrich estimate for its non-dissipative part (evaluated within the Godunov step). We estimate -u.∇(p) at the middle of the time step thanks to the same TVD (total variation diminishing) slopes used in the Godunov step. Finally, we recover q by subtraction of the left hand side. The heating measured thanks to this method 1 for n = 128 and u = 1 is the blue dashed line on figure A.1. It overestimates the actual total heating by about 60%, but provides roughly the correct shape for the heating: positive and symmetric, only too wide spread compared to the exact solution. Note that the accuracy for the total heating is much better in higher Mach number shocks (the relative error falls down to only 13% for an inflow speed u = 4).

Method 2 ∂ t ( 1 2 ρu 2 + 1 2 B 2 + p log ρ) + ∇.F 2 = -q where F 2 = F 1 + up(log ρ + 1). We compute the flux as in method 1 (the additional contribution is computed in the Godunov step using a Lax-Friedrich estimate). This method recovers the total heating through the shock with a much better accuracy (to within 4% for n = 128 and u = 1), but the shape of the heating is asymmetric, with strong negative values on the pre-shock side.