D ' Abord

Roland

Merci Pedro

Melissa

Özgür Shaokang

Alessandro Pierre

Nicole Olivier

Ana Abhyu

Victor Montagud-Camps
email: victor.montagud-camps@lpp.polytechnique.fr

Roland Grappin
email: roland.grappin@lpp.polytechnique.fr

Andrea Verdini
email: verdini@arcetri.astro.it

Turbulent Heating between 0.2 and 1 au: A Numerical Study

Keywords: magnetohydrodynamics (MHD)methods, numericalsolar wind -Sun, heliosphereturbulence

qui est présenté dans les pages suivants a été le résultat de trois ans de thèse. Mes remerciements sont dirigés à ceux qui ont fait de mon séjour en France et de mon travail, une expérience tellement bonne et inoubliable.

C Résumé en français

List of Figures

1.1 Radial structure of atmosphere and wind: (a) density (dashed line) and temperature (solide line) profiles close to solar surface [START_REF] Withbroe | Mechanical heating in the transition region[END_REF]; (b) radio-scintillation measurements (points in the image) of the solar wind bulk speed in the first 50 solar radii (from [START_REF] Scott | Solar wind observations near the sun using interplanetary scintillation[END_REF]). Solid lines correspond to theoretical models for the bulk speed evolution. Each kind of symbol for the measured points correspond to a different author. The references for these authors can be found in the cited paper. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 (Left) Eclipse white-light image of the corona [START_REF] Sr Habbal | Localized Enhancements of Fe +10Density in the Corona as Observed in Fe xi789.2 nm during the 2006 March 29 Total Solar Eclipse[END_REF]. (Right) Coronal magnetic field-line configuration on 6 February 2004, as derived by applying a PFSS (Potential Field Source Surface) extrapolation with source surface at 2.5R s to a photospheric magnetic map from Mount Wilson observatory [START_REF] Wang | Observations of Low-Latitude Coronal Plumes[END_REF]. . . . 1.3 Three Ulysses spacecraft orbits of the Sun. Figure shows radial solar wind velocity and images of the Sun at varying degrees of solar activity [START_REF] Mccomas | Weaker solar wind from the polar coronal holes and the whole Sun[END_REF]. . . . 1.4 Scheme of an electron velocity distribution function at 1AU for high speed solar wind electrons. Bottom panel: isocontours in the plane of velocities parallel and perpendicular to the mean magnetic field (Green circle represents the limit of the Maxwellian core). Top panel: parallel (solid black line) and perpendicular (dashed blue line) cross section of the VDF in the bottom panel. Dashed red line is a Maxwellian distribution that fits the core of the VDF.

From [START_REF] Pierrard | Solar Wind Electron Transport: Interplanetary Electric Field and Heat Conduction[END_REF] without distinction between fast and slow winds that form the so called Maltese cross (from: [START_REF] William H Matthaeus | Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind[END_REF]). Center: same kind of measurements done at higher frequencies and setting appart slow (up) and fast wind (down) data (from: [START_REF] Dasso | Anisotropy in Fast and Slow Solar Wind Fluctuations[END_REF]). Right: reduced two dimensional spectra computed from the integration of the 3D magnetic spectra in the azimuthal direction to B 0 . They have been computed from slow (up) and fast winds (down) at 1AU using multi-spacecraft techniques (from: [START_REF] Owen W Roberts | Three-dimensional density and compressible magnetic structure in solar wind turbulence[END_REF]) . . averaged in the octave around the two waveperiods of one hour (two upper curves) and 3 minutes (two lower curves). Thick line: proton thermal velocity in arbitrary units. From [START_REF] Grappin | Alfvenic' versus 'standard' turbulence in the solar wind[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Measurement of proton temperature as a function of radial distance with Helios 1 and 2 data for winds with mean speed over 600km s -1 . Solid line corresponds to the fitted function T p ≈ 2.5 • 10 5 (R/R 0 ) -0.74 K. From [START_REF] Hellinger | Heating and cooling of protons in the fast solar wind between 0.3 and 1 AU: Helios revisited[END_REF]. . . 4.7 Kolmogorov-like estimate of turbulent heating Q K41 (eq. 3.13) at frequency f = 10 -2 Hz (red symbols) and critical heating rate Q c (blue symbols), vs proton temperature at 1 AU ( [START_REF] Bernard | Evaluation of the turbulent energy cascade rates from the upper inertial range in the solar wind at 1 AU[END_REF]). Solid lines for each data population correspond to fitted power laws with proton temperature. . . . . . . . . . . . 4.8 Helios mission, solar minimum. Kolmogorov-like estimate of turbulent heating Q K41 (eq. 3.13) evaluated at frequency f = 4 10 -3 Hz (crosses) and critical heating rate Q c (lines), vs proton temperature. From left to right panels: increasing distance range from 0.3 to 1 AU. This is to be compared with the previous figure 4. [START_REF] Belcher | Large-amplitude Alfvén waves in the interplanetary medium, 2[END_REF] which gives a similar analysis using a much larger dataset at 1 AU. Vertical units: [10 at 1AU for slow and fast winds. Turbulent heating is computed via thirdorder moments laws that have assumed an isotropic hydrodynamic turbulence (green triangles), isotropic MHD turbulence (blue squares) or anisotropic MHD turbulence (red circles). From [START_REF] Stawarz | THE TURBULENT CASCADE AND PROTON HEATING IN THE SOLAR WIND AT 1 AU[END_REF]. . . . . . . . . . . . . . . . . . . . . 6.1 Initial total energy spectra, 2D slice of the angular spectrum in the plane k z = 0. Left: Isotropic (ISO) initial spectrum; center: Gyrotropic (GYRO) spectrum; right: Gyro-Alfvén (G+A) spectrum. . . . . . . . . . . . . . . . . 7.1 2D slice of the angular spectrum in the plane k z = 0 at the end of an ISO simulation, i.e. E 3D (k x , k y , k z = 0). The solid straight line corresponds to the largest vector in modulus from the center to the isocontour of level 10 -8 , k maj , while the dotted straight line is the smaller vector in modulus from the center to the same isocontour, k min . The mean magnetic field axis (not showed) forms an angle of 45 • with respect to the radial at the end of all simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Temporal evolution of the 1D total energy spectrum in the x direction for an ISO simulation. Each thin solid line correspond to the reduced spectra compensated by k 5/3 and averaged over one fifth of duration of the simulation. The arrow indicates the temporal evolution for the spectra: first time interval on top and the last one at the bottom. The solid thick line below the spectra helps to signal the length of the inertial range at the last time interval. Values in the vertical axis are arbitrarily chosen to help visualize the spectral evolution and are not representative of the real decrease of the amplitude of the fluctuations. . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Varying : ISO initial spectrum, σ 0 c = 0, spectral evolution of runs R12 ( = 0, first row), R2 ( = 0.2, second row), R5 ( = 0.4, third row). Rows 1 to 3: left panel, 1D energy spectra averaged in five time intervals, compensated by k -5/3 , vs k x ; mid panel, same vs k z ; right panel, 2D slices of the final 3D spectra E 3D (k x , k y , k z = 0). Last row: spectral aspect ratio A(k) and angle α(k) between main axis and radial (R12: solid line, R2: dotted, R5: dashed). 8.2 Varying : ISO initial spectrum, σ 0 c = 0.8, R x = 1, spectral evolution of runs R13 ( = 0, first row), R4 ( = 0.2, second row), R6 ( = 0.4, third row). Same caption as in fig. 8 (c,d) RB3, = 100; (e,f) RB5, = 0.5. The top row (a,c,e) shows the evolution of the energy spectra compensated by k -2 , line styles vary according to the increase in time, from a thin solid line at t = 0, followed be the dotted and dashed line and so on. The bottom row (b,d,f) shows the evolution of the dissipation terms: total dissipation dE/dt (solid line), viscous dissipation ν(∂ x u) 2 (dotted line), damping due to the expansion 2n a u 2 (dashed line), prediction for turbulent dissipation k 0 u 3 /a α (dot-dashed line). The turbulent dissipation always matches the viscous dissipation term after a nonlinear time

List of Tables

1. [START_REF] Alexandrova | Comment on "Evidence of a Cascade and Dissipation of Solar-Wind Turbulence at the Electron Gyroscale[END_REF] Average properties of fast and slow solar winds measured at 1AU made by [START_REF] Hundhausen | Direct Observations of Solar-Wind Particles[END_REF] and [START_REF] Milan Maksimovic | Etude d'un modele cinetique du vent solaire et spectroscopie du bruit thermique sur ulysses[END_REF] Adapted from a table from [START_REF] Lemaire | Kinetic Models of Solar and Polar Winds[END_REF]. . . . . . . . . . . . . . . . . . .

6.1 List of parameters for the initial conditions. Run: name of run; Type: class of initial conditions (see figure 6.1); R x : initial aspect ratio of the numerical domain; k max : the maximum wavenumber in directions perpendicular to radial; k min x : the minimum wavenumber in the x direction (k min y =k min z =1 always); m 0 : initial index of the reduced spectrum; σ c : initial cross-helicity; : initial value of the expansion parameter; M ach = u rms /c s , with c s the sound speed; B 0 : initial mean magnetic field; µ 0 fixes the initial values of the three diffusive parameters µ, η and κ (see eq. 2.30). . . . . . . . . . . . . . . . . . 10.1 Run column contains the names used for the simulations; Spectrum designates the initial energy spectrum; = (U 0 /R 0 )/(k 0 u rms ) denotes the expansion parameter; k max stands for the maximum wave-number containing energy at the beginning of the simulation; x/y defines the direction of the wavevector, radial or transverse; the pair (α, n) varies, depending on the wavevector direction: (0, 1/3) for the radial direction, (1, 5/6) for the transverse direction and (0, 0) without expansion; µ 0 is the common initial value of the diffusive parameters (viscosity and conduction). . . . . . . . . . . . . . . . . . . . . .

12.1 List of parameters for the initial conditions. Run: name of run; Type: class of initial conditions (see figure 6.1); R x : initial aspect ratio of the numerical domain; k max : the maximum wavenumber in directions perpendicular to radial; m 0 : initial slope of the reduced spectrum; σ c : initial cross-helicity; : initial value of the expansion parameter; M = u rms /c s , with c s the sound; ν 0 : initial value of the kinematic viscosity (equal to η and κ). . . . . . . . . .

Part I

Introduction 13

The solar system is embedded in an extension of the Sun's high atmosphere, known as the solar wind. An early hint that the medium is indeed a wind of ionized particles (a plasma) blowing from the Sun, comes from observations of particle acceleration and ionization in comet tails.

Plasmas are composed by a quasi-neutral set of charged particles, that is, with approximately the same amount of positive and negative charges. The movement of these particles is determined by the possible presence of external electromagnetic fields combined with the fields generated by the movement of the particles. This co-dependency between particles and fields gives rise to a collective behavior for particles and fields that characterizes plasmas and differentiates them from other states of matter.

The composition of the solar plasma is made mainly of protons and electrons, and a small concentration of Helium ions (alpha particles). In the following, we will describe the solar wind plasma via the one-fluid equations of magnetohydrodynamic (MHD). Right after that, we will present the modified version of these equations that we use in our numerical simulations, the so-called Expanding Box Model.

This work is centered on the study of turbulence in the solar wind plasma. Once the equations describing the solar wind are introduced we will describe fluid turbulence and some of its features that will be studied in this work.

After this general introduction to plasmas and fluid turbulence, we then focus on the medium that the previous theoretical framework is going to describe, the solar wind. We will start with an overview of its origins close to the Sun's surface. Then we will explain how in-situ measurements of the solar wind are done and what those measurements tell us about the two properties of the solar wind that we have studied: turbulent anisotropy in the inertial range and turbulent heating.

Chapter 1

Solar Wind

The first hydrodynamic theory of the origin of the solar wind was published by Parker [START_REF] Parker | Dynamics of the Interplanetary Gas and Magnetic Fields[END_REF]. He remarked that a strong difference in thermal pressure between the high atmosphere (corona) and interstellar plasma could lead to a global supersonic wind at a few solar radii. Parker's theory required a source of internal (thermal) energy which had to be transported from surface to the high corona and at larger distances. Direct in situ spacecraft observations then confirmed the existence of the supersonic solar wind which cools with distance more slowly than predicted by the fluid equations.

One of the sources of the solar wind heating required by Parker's theory is the thermal conduction: heat is thus transported from a 10 6 K corona through the solar wind predominantly by electrons. Another source of heat, which we shall study in detail in this thesis, is provided by turbulent dissipation of fluctuations transported by the solar wind.

Sources of the wind and solar cycle

The solar atmosphere, mainly composed of protons, electrons and alpha particles, is strongly stratified in the vertical (radial) direction: density drops from n 10 17 cm -3 at the surface down to about 5 -15 cm -3 at the Earth's orbit. The strongest stratification occurs close to the Sun. One identifies at least two layers, the thin, dense chromosphere and the hot, rarefied corona which extends up to several solar radii. The chromosphere is separated from the corona by an abrupt transition (chromospheric transition), at about 2000 km above the photosphere. While the chromosphere is relatively cold (between 6000K and 10000K) the corona has a proton temperature around 2MK. High in the corona, at several solar radii, the plasma begins to be accelerated to outer distances from the Sun and to form the solar wind.

Fig. 1.1 shows the average density and temperature profiles around the chromospheric transition (panel a) and bulk velocity measurements using radio source scintillation measurements ( [START_REF] Scott | Solar wind observations near the sun using interplanetary scintillation[END_REF], panel b).

Not all the solar surface is a source of wind. In fact, sources of the wind are largely controlled by the topology of the magnetic field. The atmosphere shows in white light density contrasts which actually trace the magnetic structures, more specifically, closed and open structures. An example of density structures during an eclipse is shown in fig. 1.2 (left), while a reconstruction of magnetic field lines is shown in fig. 1.2 (right).

In the whole corona (up to about 2.5 R s ), the magnetic field energy dominates, and the magnetic field is thought to be organized in flux tubes that force the flow to follow the Figure 1.1: Radial structure of atmosphere and wind: (a) density (dashed line) and temperature (solide line) profiles close to solar surface [START_REF] Withbroe | Mechanical heating in the transition region[END_REF]; (b) radio-scintillation measurements (points in the image) of the solar wind bulk speed in the first 50 solar radii (from [START_REF] Scott | Solar wind observations near the sun using interplanetary scintillation[END_REF]). Solid lines correspond to theoretical models for the bulk speed evolution. Each kind of symbol for the measured points correspond to a different author. The references for these authors can be found in the cited paper. This global magnetic pattern changes in a quasi-regular way during the solar cycle (see fig 1.3). The Sun's magnetic field adopts a quasi-dipolar pattern, with two magnetic poles with opposite polarities, during the so-called Solar minimum activity, and a multipolar pattern during Solar maximum activity. The quasi-dipole associated with minimum activity reverses each eleven years in average. At the time of maximum activity on the contrary, no clear polarity is defined, Solar spots, cold regions of the photosphere with a high concentration of magnetic field lines are frequently seen. Strong solar events such as Coronal Mass Ejections (CME) are abundant during this period and perturb properties of the pristine solar wind.

Independently of the period of solar activity, the intensity of the magnetic field decreases with distance. From about 10 solar radii onwards, the magnetic field energy becomes subdominant compared to the kinetic energy of the flow (which defines the Alfvén point or Alfvén surface), and the magnetic field lines cease to control the flow.

From the Alfvén point onwards, the mean magnetic field lines no longer co-rotate with the Sun, but spiral around it. This forms what is known as Parker's spiral.

Fast and slow winds

One can see important variations in the properties of the plasma depending on its source region, that is, the place where the solar wind originates. The winds originating from sources near the closed magnetic regions are characterized by a mean speed around 350 km s -1 when measured at 1AU . At the same distance, the particle density per unit volume is in average 15 cm -3 , the proton temperature is around 5 • 10 4 K and the mean magnetic field modulus is 6 nT ( [START_REF] Bruno | The Solar Wind as a Turbulence Laboratory[END_REF]).

On the other hand, winds from sources well within open magnetic regions reach speeds close to 600 kms -1 , lower particle densities, 4 cm -3 , larger proton temperatures, 2 • 10 5 K, but the mean magnetic field is the same as when the source is close to closed regions. Due to these differences we shall group solar winds into slow and fast winds, or into cold and hot winds respectively. Note that this difference between cold and hot winds is only based on proton temperature: electrons have similar temperature in fast and slow winds T e ≈ 1 • 10 5 K, (see table 1.1). The radial variation of electron temperature between 0.3 and 1 AU, T e ∝ R -0.59±0.32 [START_REF] Maksimovic | Radial evolution of the electron distribution functions in the fast solar wind between 0.3 and 1.5 AU[END_REF] is also slower than that of protons, T p ∝ R -0.9±0.1 [START_REF] Lynn Totten | An empirical determination of the polytropic index for the free-streaming solar wind using Helios 1 data[END_REF] (see section 4.4.1 for a further discussion on the measurements of proton temperature evolution).

Other ion species, such as He +2 (alpha particles) also show some differences between slow and fast solar wind streams. At 1AU, their relative abundance with respect to protons is around the 4% and the temperature ratio, T α /T p , is around 1.2 for the more collisional slow winds and around 4.5 for fast winds [START_REF] Maruca | Collisional Thermalization of Hydrogen and Helium in Solar-Wind Plasma[END_REF]. Bulk speed of alpha particles also shows differences between slow and fast streams, as it is smaller than that of protons for slow streams and higher for fast winds, the difference reaching 170km s -1 [START_REF] Marsch | Solar wind helium ions -Observations of the HELIOS solar probes between 0.3 and 1 AU[END_REF].

Until now we have just specified the amplitudes of macroscopic (i.e. fluid) properties in fast and slow winds. There are also differences in microscopic quantities, such as velocity distribution functions. While proton velocity distribution functions (VDF) for slow winds are similar to a Maxwellian distribution, fast winds VDFs are gyrotropic with respect to the mean magnetic field, as can be seen in figure 1.5. In contrast, electron VDFs do not present strong temperature anisotropies such as proton VDFs (see table 1.1). Protons, electrons and alpha particles present departures from a gaussian distribution in the form of non-thermal tails, as the one shown in figure 1.4. The generation of non-thermal tails has been linked to the wave-particle interactions of high frequency waves in plasmas such as Kinetic Alfvén Waves (KAWs) and whistler waves: KAWs are invoked in the generation of the non-thermal tail of proton distributions ( [START_REF] Pierrard | Modification of Proton Velocity Distributions by Alfvénic Turbulence in the Solar Wind[END_REF] and references within), and whistler waves in the generation of the non-thermal tail of electron distributions [105, 104][78].

In this thesis, we restrict ourselves to the study of solar wind plasma at low frequencies (correspondingly large scales), that can be well described by MHD equations (see section 2.1). Thus, high frequency waves and their effects on the particles VDFs are excluded from our description of the plasma. Also, the deviations of the VDFs from Maxwellian distributions will be minimized: in particular, perpendicular and parallel temperatures will be considered equal, as well as electron and proton temperatures.

From corona to Earth

Starting from the corona, the solar wind speed becomes larger than the sound speed beyond one solar radius, and larger than the Alfvén speed beyond ten solar radii (see fig. 1.6). After this distance, the wind flow is quasi-radial. The wind speed then reaches a value close to its cruise speed, e.g., it doesn't change much between 0.2 and 1 astronomical units (AU), Figure 1.4: Scheme of an electron velocity distribution function at 1AU for high speed solar wind electrons. Bottom panel: isocontours in the plane of velocities parallel and perpendicular to the mean magnetic field (Green circle represents the limit of the Maxwellian core). Top panel: parallel (solid black line) and perpendicular (dashed blue line) cross section of the VDF in the bottom panel. Dashed red line is a Maxwellian distribution that fits the core of the VDF. From [START_REF] Pierrard | Solar Wind Electron Transport: Interplanetary Electric Field and Heat Conduction[END_REF]. [START_REF] Hundhausen | Direct Observations of Solar-Wind Particles[END_REF] and [START_REF] Milan Maksimovic | Etude d'un modele cinetique du vent solaire et spectroscopie du bruit thermique sur ulysses[END_REF] Adapted from a table from [START_REF] Lemaire | Kinetic Models of Solar and Polar Winds[END_REF]. i.e., between 42 and 210 solar radii: this will allow us to use the expanding box model to be described below (see section 2.2). The quasi-radial flow and the quasi-constant speed forces a plasma volume to expand in the two directions perpendicular to the radial and not in the radial direction (see fig. 2.1), so that the volume expands as R 2 instead of R 3 . This increase of a plasma volume embedded in the wind has two consequences (i) density should decrease as 1/R 2 , which is indeed observed; (ii) applying the one-fluid adiabatic law, i.e., assuming pressure forces are the only ones against which the expanding plasma volume is working, with no extra internal energy source, one obtains for the temperature

T ∝ ρ γ-1 ∝ ρ 2/3 ∝ R -4/3 (1.1)
However, as we shall see in the heating section 4.4.1 below, the proton temperature gradient is flatter than that, which requires a heat source to slow down the plasma cooling during expansion.

Chapter 2

Plasma description

The MHD equations

We adopt a one-fluid description (Magnetohydrodynamics or MHD) of the Solar Wind plasma. A fluid description describes the plasma in terms of the evolution of its macroscopic variables: density, velocity, pressure, heat flux,... and the evolution of the magnetic and electric fields, B and E, given by Maxwell's equations. A fluid description can be derived in several ways, either from a kinetic description or starting directly with a fluid description, imposing a series of conservation laws (see for instance [START_REF] Chen | Introduction to Plasma Physics[END_REF]). In all cases, fluid equations are a hierarchy of equations that need to be closed by imposing the evolution of certain macroscopic variables. For an MHD description of the plasma, it is assumed quasineutrality (approximately the same density of negative and positive charges), a Gaussian velocity distribution function for all species and slow variations of the plasma fluctuations in space and time 1 , that allows to consider ∇ × B = µ 0 j. The evolution of the current j is given by assuming Ohm's law: E + u × B = ηj.

For the MHD description presented below, we also assume that the flux of thermal energy is q = -κ∇T , that the plasma follows the equation of state for ideal gases P = ρT and that the polytropic index is γ = 5/3 (thus, in the absence of local sources of heating, the system is adiabatic).

The previous assumptions made to obtain the MHD description allow to describe the plasma as a conductive fluid whose evolution is coupled to that of the magnetic field. Several observations however argue against this choice to describe the Solar Wind. The main ones are as follows: (1) electron and proton temperatures differ; (2) mostly in fast winds: the velocity distribution function is not gaussian, resulting in different values for the P ⊥ (pressure perpendicular to the mean magnetic field ) and P (pressure parallel to the mean field); (3) other species (Alpha particles in particular) show a velocity drift along the mean field.

These features all derive from the low rate of Coulomb collisions (as measured by the Knudsen number 2 ) starting from the corona. The justification for neglecting the effects of low collisionality in a first approach is because no direct simulations of the 3D turbulent cascade and the associated energy flux has been done yet properly for the solar wind (see next section 2.2). So it is valuable to investigate if, in spite of the simplifications of the MHD description, turbulent heating and spectral anisotropy properties that will be obtained from simulations are already close or not to the solar wind observations. Finally, we measure the energy cascade rate (to be defined later in Section 3.1.3) and the associated visco-resistive dissipation rate, defined in eq. 2.2 below. The actual dissipation in the solar wind occurs at kinetic scales. Our MHD description replaces these scales by visco-resistive scales: this will be justified later in Section 3.1.2).

The MHD equations read as follows,

∂ t ρ + ∇ • (nu) = 0 ρ(∂ t u + u • ∇u) + ∇P -(∇ × B) × B = µ(∆u + 1 3 ∇(∇ • u)) ∂ t B + u • ∇B = B • ∇u -B(∇ • u) + η∆B ∂ t P + u.∇P + γP (∇.u) = κ∆T + (γ -1)Q ∇ • B = 0 (2.1) Q = µ((∇ × u) 2 + 4 3 (∇ • u) 2 ) + η(∇ × B) 2 (2.2)
The last term at the left-hand side of the momentum equation in (2.1) can be decomposed into

(∇ × B) × B = -∇(B 2 /2) + (B • ∇)B 2 (2.3)
where the first term contains the magnetic pressure and the second one is the magnetic tension. In front of a perturbation of the mean magnetic field, these terms will try to restore straight magnetic field lines. Due to this, kinetic and magnetic oscillations propagate in the plasma. From the linearized MHD equations 2.1 one can obtain the dispersion equation of the MHD waves.

Apart from the entropy mode that does not propagate (ω = 0), propagating modes are Alfvén waves and the fast and slow magnetosonic waves. Alfvén waves are of particular interest, as they are directly identified in the solar wind, via their eigenmode relation, in the so-called Alfvénic streams. They are obtained from the linearized MHD equations by imposing that the density, pressure, magnetic and velocity fluctuation parallel to the mean magnetic field are zero. The remaining equations for the fluctuating amplitudes u and b

= (B -B 0 )/ √ ρ 0 = δB/ √ ρ 0 are ∂ t u = (V a • ∇)b (2.4) ∂ t b = (V a • ∇)u (2.5)
with phase speed being the Alfvén speed

V a = B 0 √ ρ 0 (2.6) and dispersion relation ω = ±k • V a (2.7)
The previous propagation equation may be rewritten as

∂ t z ± ∓ (V a • ∇)z ± = 0 (2.8) with z ± = u ± b (2.9)
For a mean magnetic field line pointing towards the Sun, z + propagates outwards the Sun while z -propagates inwards. The Alfvénic streams mentioned above thus have either with z + z -or z - z + , depending on the polarity of the mean field B 0 . The dispersion relation of magnetosonic modes is:

(ω/k) 2 = 1 2 [(v 2 A + c 2 s ) ± (v 2 A + c 2 s ) -4v 2 A c 2 s cos 2 θ] (2.10)
where c s = γP ρ is the sound speed and θ is the angle between k and B 0 .

MHD equations with large scale radial flow (EBM)

We describe here the expanding box model (EBM) defined in ( [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF]), which is the method that will be used in this thesis to simulate turbulent evolution in the solar wind. More precisely, the model will allow us to follow a turbulent plasma volume transported by the mean solar wind between say 0.2 and 1 AU. As we will show, it consists in subtracting the mean radial flow from the velocity field in the MHD equations presented in section 2.1. We begin by explaining the motivation. As explained in section 3, turbulence redistributes kinetic and magnetic energies among scales, usually from large to small scales. The problem is that to follow the turbulent cascade process between 0.2 and 1 AU in an Eulerian frame requires locally (at a fixed distance, with no wind) a large number of mesh points (typically 512 3 or 1024 3 ) but to follow the cascade from 0.2 AU up to the Earth's orbit (1 AU) requires in principle much more, as (i) the largest energy containing eddies have a frequency of f=4 10 -4 Hz i.e., a period τ = 2500s; (ii) the time during which we want to follow the plasma evolution is about 4 days. The ratio between the two times is about 144, which means either that we leave any hope to follow turbulence correctly, or that we ignore the effect of the transport by the wind.

A solution to this problem consists in employing a technique often used to simulate the flow evolution in a wind tunnel: it amounts to subtract the average flow, changing to the Galilean frame transported by the mean flow. This transforms the distance evolution into time evolution. Of course, since we ignore the boundary conditions of our advected domain, we must choose the only possible method in this case, that is, adopt periodic boundary conditions.

The EBM approach relies on the idea that the change of Galilean frame described above is but the first step to eliminate the mean flow. The reason is that the mean flow, contrary to that of a wind tunnel, is not uniform, it is radial:

U 0 = U 0 êr (2.11)
As a consequence, after our Galilean frame change, the plasma (as described with local cartesian coordinates) actually still expands in directions perpendicular to the radial (see fig. 2.1). To get rid of this residual flow, we subtract this transverse expansion by adopting coordinates comobile with it. We obtain in this way a set of modified MHD equations, with both linear and nonlinear modifications, which predict the evolution of the plasma volume in a cartesian system attached to the central radial line. We now are free to use periodic boundary conditions, which closes our definition of the expanding box model, using a local cartesian coordinate system.

Note finally that, since the model relies on a mean radial flow with constant speed, it can be used only far from the acceleration range of the solar wind and, a fortiori, at distances larger than the Alfvén point. In practice, we will use it between 0.2 and 1 AU in this thesis. This choice is also supported by in-situ measurements between 0.3 and 1AU, using data from Helios missions: they have shown that the evolution of the particle density with heliospheric distance is close to n ∝ R -2 [START_REF] Lynn Totten | An empirical determination of the polytropic index for the free-streaming solar wind using Helios 1 data[END_REF] [START_REF] Hellinger | Heating and cooling of protons in the fast solar wind between 0.3 and 1 AU: Helios revisited[END_REF][START_REF] Hellinger | Proton thermal energetics in the solar wind: Helios reloaded[END_REF], thus close to the prediction of radial constant expansion obtained from the conservation of mass.

The radius R at which the plasma box is located varies with time τ as

R(τ ) = R 0 + U 0 τ (2.12)
where R 0 is the initial position of the box. Space, time, velocity, temperature and density are measured in the following units:

L 0 /(2π) (2.13) t 0 N L = L 0 /(2πu 0 rms ) (2.14) u 0 rms (2.15) m p (u 0 rms ) 2 /(2k B ) (2.16) ρ 0 (2.17)
where ρ 0 is the initial average density of the plasma, u 0 rms is the initial rms velocity of the fluctuations, t 0 N L is the initial nonlinear time based on the initial rms velocity, and L 0 is the initial size of the box perpendicular to the radial direction.

Let us choose a cartesian coordinate system (X,Y,Z), with the X direction aligned with the Sun-Earth radius passing by the center of the initial plasma volume. Y and Z are taken to complete the Radial Tangential Normal coordinate system. The set of coordinates comobile with the transverse expansion is then

t = τ (2.18) x = (X -U 0 τ )/R x (2.19) y = Y /a(t) (2.20) z = Z/a(t) (2.21)
We choose the unit length to be the initial transverse size (i.e., L 0 Y = L 0 Z ) of the domain divided by 2π) (by transverse we will always mean "perpendicular to the radial" in the following, while "perpendicular" and "parallel" will mean "with respect to the mean magnetic field"). Note that we have chosen a regular box as our physical domain. Curvature terms from the spherical expansion of the wind can be neglected if we assume L X,Y,Z /R 1 ( [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF]). The parameter R x is the initial aspect ratio of the domain:

R x = L X /L 0 Y = L X /L 0 Z (2.22)
Remark that L X is constant as the wind speed is assumed constant, and therefore L X = L 0 X . The parameter a(t) is the normalized heliospheric distance:

a = R(t)/R 0 = 1 + t (2.23)
where = da/dt is the initial expansion rate defined as the initial ratio between the characteristic expansion and turnover times in the transverse directions (perpendicular to the From [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF].

radial):

= τ N L τ exp = U 0 /R 0 k 0 u 0 rms (2.24)
with k 0 the minimum wavenumber in the transverse direction. At a given distance R(t), the domain aspect ratio is thus decreasing as:

L X /L Y = R x /a(t) (2.25)
The choice of the initial aspect ratio R x of the domain is somewhat arbitrary: in the following, we will consider R x = 1 and R x = 5. If we choose R x = 1, the aspect ratio varies from unity to 1/5 (this is the choice made in Dong et al 2014). If we choose R x = 5, then the aspect ratio (when traveling from 0.2 to 1 au) varies from 5 to unity (this is the choice made in [START_REF] Verdini | IMPRINTS OF EXPANSION ON THE LOCAL ANISOTROPY OF SOLAR WIND TURBULENCE[END_REF], [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF] and in [START_REF] Montagud-Camps | Turbulent Heating between 0.2 and 1 au: A Numerical Study[END_REF]).

The EBM equations read

∂ t ρ + ∇(ρu) = -2ρ a ∂ t u + u.∇u + 1 ρ ∇(P + B 2 2 ) - B.∇B ρ = µ ρ ( ∇2 u + 1 3 ∇( ∇.u)) -U a ∂ t B + u.∇B -B.∇u + B(∇.u) = η ∇2 B -B a ∇ • B = 0 ∂ t P + u.∇P + γP (∇.u) = -2γP a + ρκ ∆T + (γ -1)ρQ ν Q ν = µ(ω 2 + 4/3 ( ∇ • u) 2 ) + η J2 P = ρT (2.26)
where U = (0, u y , u z ), B = (2B x , B y , B z ) are obtained from terms involving the gradient and the divergence of the mean perpendicular flow, ∇ • U 0⊥ ≈ 2 /a (see for instance [START_REF] Dong | EVOLUTION OF TURBULENCE IN THE EXPANDING SOLAR WIND, A NUMERICAL STUDY[END_REF]).

As in section 2.1, we continue to assume that in the absence of dissipation and heat flux, the system is adiabatic, γ = 5/3. In the heating term, the vorticity and current are defined in terms of the modified nabla operator ∇: ω = ∇ × u and J = ∇ × B (see eq. 2.29 for the definitions of the modified nabla operator).

In these equations, ρ is the density normalized by its average:

ρ = ρ/ρ = a 2 ρ (2.27)
and ρ = 1/a 2 is the average density of the plasma. The standard nabla operator is related to gradients expressed in terms of comobile coordinates x,y,z as:

∇ = ((1/R x )∂ x , (1/a)∂ y , (1/a)∂ z ) (2.28)
The nabla operator appearing in the dissipative terms in the right hand side of eqs. 2.26 is defined so as to adapt to the comobile coordinates x,y,z, in order to prevent the excessive damping of fluctuations perpendicular to the radial:

∇ = (∂ x , ∂ y , ∂ z ) (2.29)
In order to prevent a too fast decrease with time of the Reynolds number (defined in section 3.1.2), we adopt the following temporal dependance of the diffusion coefficients:

µ = η = κ = µ 0 /a(t) (2.30)
The definitions in equations (2.29) and (2.30) do not respond to physical but to numerical reasons. On the one hand, the loss of small perpendicular scales due to expansion and the limitations imposed by the use of an MHD description impedes us to properly describe small scale physics, in particular the dissipative process. Since our focus is not put on the description of the dissipation process itself, but on the amount of energy dissipated by turbulence, we can prescribe the aforementioned definitions. On the other hand, in any turbulence simulation it is necessary to avoid the artificial accumulation of energy at the smallest mesh-size, in order to prevent an unphysical equipartition of energy between all degrees of freedom. Equations (2.29) and (2.30) have been defined to do so, and at the same time maintain a high Reynolds number.

Chapter 3

Turbulence

Homogeneous turbulence

A qualitative picture of turbulence in every day life is given by random patterns (usually eddies) forming in the wake of a solid obstacle facing a flow, or air turbulence forming for instance above a tea cup, as in fig. 3.1(Left). Turbulence can be found in a variety of contexts, such as geophysical flows ( [START_REF] Lilly | Two Dimensional Turbulence Generated by Energy Sources at Two Scales[END_REF]), gravitational waves ( [START_REF] Galtier | Turbulence of Weak Gravitational Waves in the Early Universe[END_REF]), interstellar medium ( [START_REF] Hennebelle | Turbulent molecular clouds[END_REF]), or the Solar Wind ( [START_REF] Chen | INTERPRET-ING POWER ANISOTROPY MEASUREMENTS IN PLASMA TURBULENCE[END_REF]). Here, we shall concentrate on fluid turbulence, with fluids that can be described either by Navier-Stokes or MHD equations.

The simplest case is that of homogeneous and isotropic fluid turbulence, that is, with nothing destroying a priori homogeneity and isotropy of space. It shows two main properties (i) a quasi-constant flux of energy (or related quantity) flowing in general from large to small scales (in that case one speaks of direct cascade); (ii) a systematic randomness which prevents one to exhibit any analytic solution of a turbulent system and asks for a statistical description of the phenomenon, a situation analogous to that of thermodynamical equilibrium for instance.

Shock formation as a simplified model of turbulence

The flux of energy from large to small scales is well illustrated by a related phenomenon: shock formation. It is illustrated in fig. 3.1(right) by a numerical solution of the 1D compressible Navier-stokes equations,

∂ t ρ + ∂ x (ρu) = 0 (3.1) ∂ t (ρu) + ∂ x (ρu 2 ) + ∂ x P = µ∂ xx u (3.
2) Kinetic energy is thus finally distributed on a whole range of wave numbers, while the initial energy was concentrated on a single wavenumber: this evolution is due to the nonlinear term, namely the advection term. Moreover, as shown in panel (d), while during the first part of the evolution the total kinetic energy is more or less constant, it begins to systematically decrease when t > 2, after the shock has formed (panel (a)). The shock formation process may be considered a deterministic version of what is called "turbulent dissipation" which begins when energy has been distributed among various scales, including the so-called dissipation scales where the viscous term is effective.

∂ t P + u x ∂ x P + γP ∂ x u x = κ∂ xx T (3.

3D hydrodynamic turbulence

In the previous 1D shock formation, the eddies visible in fig. 3.1(Left) are lacking completely, since only 1D movements are possible. In this case, the random character of turbulent evolution is absent, the evolution is deterministic. As a consequence, another property of turbulence which we explain now is lacking: the locality of the kinetic energy flux which allows to predict how energy is distributed among scales.

Assume indeed that we deal now with incompressible 3D hydrodynamic equations. The equations are similar to the incompressible 3D MHD equations, but with no magnetic field and with the pressure equation replaced by a divergenceless condition on the velocity field:

∂ t u + u • ∇u + ∇P = ν∆u (3.5) ∇ • u = 0 (3.6)
The divergenceless (incompressibility) condition is obtained from the compressible equations (see MHD eqs.2.1 with B = 0) by taking the limit of infinite sound speed (or very subsonic flow). The situation is not very different from the previous 1D compressible Navier-Stokes equations (3.4), except for two things: (i) the incompressibility conditions implies that the total kinetic energy of the system is conserved by nonlinear couplings, while in the 1D compressible case this conservation was only approximately satisfied; (ii) the 3D configuration in practice transforms completely the evolution by making it chaotic: the precise spatial and temporal distribution of large gradients, contrary to the shock case seen previously, is unpredictable. A consequence is that no analytic (or even semi-analytic) solution of the problem exists presently. This implies in particular that in 3D turbulence, no simple typical profile can be produced such as the one shown for a 1D shock at time t = 2 in fig. 3.1a.

Fourier space description: K41 phenomenology

Happily enough, the random character of turbulence does not prevent one to find average quantities which allow a description of the statistical equilibrium. This approach started with [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF]. Consider first the hydrodynamic system with no viscosity or, which is equivalent, with only large scale velocity fluctuations, for which the viscous term is negligible. This means that the Reynolds number, which is the ratio of the nonlinear term evaluated at the energy containing scales, L, and of the typical viscous term, is much higher than unity:

Re = (u • ∇)u/(ν∆u) 1 (3.7) or Re = U L/ν 1 (3.8)
In this situation, the kinetic energy is an invariant, which nonlinear couplings (associated to the advection term) redistribute among the different Fourier modes of the system. Since at small scales (high wavenumbers) modes are more and more numerous, these degrees of freedom rapidly gain energy, taking it from large scale ones. But, again, kinetic energy as a whole is conserved, at least during a first phase in which viscous dissipation can be neglected. So, if we would have strictly ν = 0, a Gibbs-like equilibrium with equal distribution of energy between all modes would rapidly be reached, thus with the small scales sharing the largest amount of energy.

In the real physical world with non zero viscosity, the transfer of energy from the large to the small scales is moderated by the fact that finally, as small scales become populated, the viscous terms become active in a range of scales called dissipative range.

In between the largest scales and the dissipative scales, a qualitative description of this turbulent equilibrium has been proposed by Richardson [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF]. According to his description, energy is injected at large scales in the form of whirls of characteristic size L and turnover velocity U . The large size whirls would split into smaller whirls due to the non-linear interactions imposed by the term (u∇)u in the Navier-Stokes equations. Then, the kinetic energy stored in the large scale whirls would pass to the smaller ones of size l and turnover velocities u.

It has been then proposed by Kolmogorov [48] that the process of energy redistribution between scales is local in scale space, that is, with the energy being transferred between neighboring scales l and l/2, thus making possible to estimate it quantitatively. This process of transferring energy from scale to scale was named turbulent cascade. Assuming that the dissipation of energy only takes place at small scales, this process would continue until the size of the whirls reaches the dissipation scales, where the kinetic energy would finally be transmitted to the internal energy of the system. In this picture we can differentiate three scale ranges: the injection scale (with largest energy content), the inertial range (that of the turbulent cascade), and the dissipation range, where the energy flux is ultimately dissipated.

In the inertial range, viscosity can be neglected, so the energy is conserved during the cascade, so that the energy flux t between neighboring scales is scale-independent:

t = dE dt ≈ u 2 t N L ≈ u 3 l -1 = ku 3 (3.9)
where

t N L = l/u = 1/(ku) (3.10)
is the characteristic time for nonlinear interactions, or also turnover time, a term appropriate for eddies. By writing that the energy flux is scale invariant, it can then be deduced that the isotropic energy spectrum at the inertial range scales is a power law of the wavenumber

k, E k = C K 2/3 t k -5/3 (3.11)
where E k is the energy spectrum, related to the velocity as

u 2 = k √ 2 k/ √ 2 E k dk kE k . The k -5/3
scaling is thus the 3D counterpart of the k -2 energy spectrum which holds in the shock formation case considered earlier. It is known as the Kolmogorov spectrum and C K as the Kolmogorov constant, an integration constant that depends on the system in which turbulence takes place. For turbulence in 3D incompressible neutral fluids, one measures C K ≈ 1.6.

The end of the inertial range begins with scales where viscous terms begin to be comparable to nonlinear terms, that is for which u • ∇u ν∆u or ku νk 2 , also u νk. Using the Kolmogorov spectrum (eq. 3.11) and u 2 kE k , we finally obtain the dissipative scale

k D 1/4 t ν -3/4 = k 0 Re 3/4 (3.12)
which is seen to go to infinity when viscosity goes to zero (or, equivalently, when the Reynolds number goes to infinity), which means that the dissipation scale goes to zero in this limit, and the inertial range where the k -5/3 law holds extends without limit. Eq. 3.12 will be used to normalize the spectra when analyzing our numerical results. Note that ν = 0 is not equivalent to a vanishing viscosity: the spectral slope k -5/3 holds only in the latter (physical) case.

An important point is the clear separation between the inertial range and the dissipation range, which ensures that one can indeed measure the energy flux , using eq. 3.9. This is possible only if the Reynolds number is large enough, i.e., when viscosity is small. Such a situation, that is, a large inertial range allowing to measure the energy flux via eq. 3.9), will be assumed to hold in more general situations, whatever the precise form of the dissipation term which may differ from the standard viscous term. The only assumption is that dissipation is completely non negligible only at the smallest scales, as are the viscoresistive terms, which scale with wavenumber k as k 2 u. This assumption is basic in our following work on turbulent heating, where the dissipation process is still unknown, but yet believed to be active at small scales (corresponding to frequencies larger than 1 Hz, see fig. 4.1). This assumption will allow us to integrate numerically the MHD equations with visco-resistive terms and then to use the visco-resistive dissipation rate in the simulations to estimate the turbulent dissipation in the solar wind.

3D MHD turbulence

We again consider here the incompressible limit of the equations, as it is only in this limit that a turbulent phenomenology is easily built. In this limit, one can prove that there are three inviscid invariants (by this we mean a quantity conserved by nonlinear couplings, in the limit of zero viscosity): (i) "total" energy, that is, the sum of kinetic and magnetic energy (in Alfvén speed units)

E v = u 2 /2 and E b = B 2 /(2ρ); (ii) cross helicity v.B/ √ ρ; (iii)
magnetic helicity a.B where a is the magnetic potential (B = ∇ × a).

Let us define the fluctuating magnetic field as δB = B -B 0 , where B 0 is the mean magnetic field. The standard view is that the turbulent spectrum, in the MHD case, is formed as in the hydrodynamic case by the cascade of the total energy. Thus, applying the Kolmogorov phenomenology to the total energy, i.e., with E(k) replacing the kinetic energy spectrum, we obtain again E(k) ∝ k -5/3 . More precisely, the total energy flux between scales reads

Q K41 = dE dt ≈ u 2 + δB 2 /ρ t N L ≈ k(u 2 + δB 2 /ρ) 3/2 (3.13)
where the non-linear time is

t N L = 1/(k(u 2 + δB 2 /ρ) 1/2 ) (3.14)

Incompressible MHD using Elsasser variables, cross-helicity

It is useful to recast MHD equations in a new form, still in the incompressible limit ∇•u = 0. We assume here (which is standard in the incompressible limit) ρ = constant and define u and b = (B -B 0 )/ √ ρ to be the kinetic and magnetic fluctuations in Alfvén speed unit, with the Alfvén speed being v A = B 0 / √ ρ. We now define the Elsasser variables as:

z ± = u ± b (3.15)
One can then write the MHD equations in the incompressible limit as:

∂ t z ± ∓ (v A • ∇)z ± + (z ∓ • ∇)z ± + 1 ρ ∇P = 0 (3.16) ∇ • z ± = 0 (3.17)
where the pressure P is the total pressure (actually, pressure is no longer an independent degree of freedom, since the incompressible condition ∇ • z ± = 0 implies that the pressure is a functional of the z ± amplitudes). For the sake of clarity, we have left the diffusive terms aside.

The linearized equations read:

∂ t z ± ∓ (v A • ∇)z ± = 0 (3.18)
These equations describe the propagation of Alfvén waves along the mean field B 0 , the mode z + propagating opposite to the mean field and z -in the direction of the mean field. Both are Alfvén waves (that is, with dispersion relation ω = ±k • v A ), corresponding to two of the three modes of compressible MHD, namely true Alfvén waves and slow magnetosonic waves (the fast modes have disappeared since the incompressible limit corresponds to infinite sound speed). Energies E ± =< z 2 ± /2 > in each mode are conserved separately by nonlinear terms (this is equivalent to the separate conservation mentioned above of total energy and v.b crosshelicity). This should lead to two separate cascades, one for each of the two energies E ± , with each its own energy flux. Each energy flux may be written formally in the Kolmogorov line:

Q ± K41 = dE ± dt ≈ (z ± ) 2 t ± (3.19)
with the definitions of the nonlinear times, obtained by visual inspection of eq. 3.16:

t ± = 1/(kz ∓ ) (3.20)
The nonlinear time (that is, the cascade time) of each Alfvén species z ± thus goes to infinity when the other species is absent. This leads to the well known paradox of fast winds (or Alfvénic winds) which are almost devoid of one of the two species, but nevertheless show well-developed energy spectra (see section 5.1). Last, we define the standard measure of the imbalance between the two (z ± ) 2 energies, called normalized cross-helicity:

σ c = |z + | 2 -|z -| 2 |z + | 2 + |z -| 2 (3.21)
Anisotropy of the cascade with mean magnetic field B 0

In general, that is, in the absence of specific forcing or boundary conditions, the 3D MHD cascade is expected to be isotropic in Fourier space. However, when a non zero mean field B 0 is present with intensity comparable or larger than the large eddies' amplitudes, then one expects Alfvén waves to introduce a global anisotropy with respect to the mean magnetic field axis. The origin of this important property lies in the following points: (i) nonlinear terms in the MHD equations couple Alfvén waves propagating in opposite directions along the mean field; (ii) the Alfvén waves frequency varies from zero to kB 0 / √ ρ when wave vectors k pass from perpendicular to parallel to the mean field B 0 ; (iii) the efficiency of the nonlinear coupling averages to zero when the nonlinear time becomes longer than the oscillation period 1/ω. As a consequence, only wave vectors with directions close to perpendicular to the mean field contribute to an efficient turbulent cascade, while the cascade is negligible in other directions (e.g., [START_REF] John V Shebalin | Anisotropy in MHD turbulence due to a mean magnetic field[END_REF], [START_REF] Grappin | Onset and decay of two-dimensional magnetohydrodynamic turbulence with velocity-magnetic field interaction[END_REF], [START_REF] Belmont | Collisionless Plasmas in Astrophysics[END_REF]). In other words, as soon as the mean field is large enough, the cascade should be mainly perpendicular, and the resulting energy spectrum should adopt an anisotropic pattern with aspect ratio given by

k /k ⊥ = b rms /B 0 < 1 (3.22)

Local and nonlocal interactions

It will be useful to comment shortly the equations in Fourier space. For the sake of simplicity, we still consider the incompressible limit:

∂ t ẑ± i (k) = ±i(k • B 0 )ẑ ± i (k) + M ijm (k) dpdq ẑ∓ j (p)ẑ ± m (q)δ(k -p -q) -νk 2 ẑ± i (k) (3.23) k • ẑ± (k) = 0 (3.24)
where ẑ± are the Fourier coefficients of the Elsasser variables. The kernel M ijm (k) stands for

M ijm = (-i/2)P ijm (3.25)
P ijm = k m P ij + k j P im (3.26)
P ij = δ ij - k i k j k 2 (3.27)
The pressure has disappeared from the equations, as, using the zero divergence condition, the pressure gradient can be expressed in terms of quadratic terms involving the Elsasser fields. Note that the non-linear terms in the right-hand side of equation (3.23) only contribute to the evolution of ẑ± i as the triad (k, p, q) verifies k = p + q (3.28)

Triadic interactions of modes are therefore responsible for the development of turbulence, at least in the incompressible limit. We can distinguish two kinds of interactions. When |k| ≈ |p| ≈ |q| we talk about local interactions since it concerns modes with wave-numbers close to each other. Conversely, when |p| << |k| |q|, the large scale wavenumber p directly drives excitation from q to k. Such nonlocal interactions dominate in general only during the early phase of turbulence onset, before the spectrum is well developed. Only later on, when the local couplings dominate, the simple K41 expression for the energy flux (eq. 3.9) can be used.

Spectra, autocorrelations and structure functions

Let us start with the definition of the velocity vector u and its Fourier transform u k :

u = u k exp ik•x dk. (3.29)
Conversely:

u k = 1 (2π) 2 u exp -ik•x dx. (3.30)
The 3D spectral density is the modulus of the Fourier coefficient,

E 3D k = |u k | 2 (3.31)
The 3D spectral density tells us how the energy is distributed among the different wavenumbers k (or scales 1/k) in the three directions of space.

Note that, most of the time, the 3D spectral density cannot be measured directly in situ. Only 1D measurements along the Sun-Earth direction can be done, which are equivalent to measuring the 1D reduced energy spectrum, E(k x ), defined as

E(k x ) = E 3D k dk y dk z (3.32)
x being the coordinate along the radial direction. Consider the case of homogeneous turbulence in neutral fluids, turbulence develops an isotropic energy spectrum that scales as a power law of the wavenumber, E(k) ∝ k -5/3 .

A similar information can be obtained from autocorrelations

R uu (l) = R 3 u(x + l)u(x)dx =< u(x + l)u(x) > (3.33)
where l is the distance between two points in real space. It can be related to spectra by

E 3D k = R uu (l) exp -ik•l dl
The structure functions (SF) of order n are defined as

S n (l) = R 3 (u(x + l) -u(x)) n dx =< (u(x + l) -u(x)) n > (3.34)
with SF of order 2 satisfying to

S 2 (l) 2 =< u 2 > -R uu (l). (3.35)
Note that in the above lines we have presented several ways to measure at a given scale the energy of the system and other quantities (autocorrelations and second order structure functions) bearing a similar physical meaning. It is possible to derive the conservation equation of these quantities from the primitive equations with just a few assumptions.

Exact laws can be derived for the evolution of second moment (autocorrelations, structure functions and/or their Fourier transforms, namely energy spectra) in terms of corresponding third order moments. However, these equations are not directly solvable, since the evolution of third order moments require to solve their evolution equations in terms of fourth order moments, etc... This is clearly not feasible, and forces us to solve direct simulations of the primitive equations, and to compute the above defined averages (structure functions, spectra) as a last step, once the primitive equations have been numerically solved.

Turbulence with large scale radial flow

Imposed large scale flows or magnetic fields can lead to new or modified turbulent regimes. An example in hydrodynamics is the case of rigid rotation which modifies both the scaling (the spectrum steepens) and anisotropy (the cascade is mainly active in the Fourier plane perpendicular to the rotation axis), and in MHD with mean field as seen in subsection 3.1.3.

In general, large scale gradients are supposed to play the role of an energy reservoir and/or help to trigger the onset of the energy cascade. In the solar wind context, for instance, the shear between fast and slow streams has been considered in this way, and, also, to trigger the decrease of cross-helicity with heliocentric distance ( [START_REF] Roberts | Magnetohydrodynamic simulation of the radial evolution and stream structure of solarwind turbulence[END_REF]).

Inhibition of the turbulent cascade by expansion

Contrary to plain (divergenceless) shear flows which trigger the cascade, the mean radial flow, due to the associated anisotropic expansion, should instead prevent or at least delay the turbulent cascade.

An analogy with gravitational collapse and cosmological expansion helps to understand the case of turbulence onset within the solar wind. Gravitational collapse is known to occur at scales larger than Jeans' scale, while smaller scales are stable, but cosmological expansion shifts the Jean's scale at larger scales, thus inhibiting gravitational collapse, a rather intuitive result (eg, [START_REF] Holian | Fragmentation by molecular dynamics: The microscopic "big bang[END_REF]). It happens that the plasma expansion during solar wind transport has a similar effect on the turbulence onset. The analogy can be understood easily if we think of quasi-stationary turbulence as a permanent formation of quasi-singularities and/or gradients. As expansion is continuously damping the existing gradients, it is clearly fighting against the cascade. The clearest example of turbulence inhibition by expansion is given by acoustic waves, for which one can demonstrate analytically the inhibition of shock formation when the wave is transverse to the radial direction ( [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF], see also part III of this thesis). In the two-dimensional case, inhibition by expansion of the turbulent cascade starting with a large scale spectrum was found numerically by [START_REF] Grappin | Onset of anisotropy and Alfvén waves turbulence in the expanding solar wind[END_REF].

The inhibition of turbulent cascade by expansion is complicated by its anisotropy: as the radial size of a plasma volume doesn't expand far enough from the acceleration region, the inhibition should be absent for a radial cascade. However, such a conclusion holds only in the special case of 1D purely compressive waves, for which indeed wavevector triads (see, eq. 3.24) are colinear: in quasi-solenoidal flows dominated by currents and vortices, triads are not colinear, so that the detailed expansion effect is difficult to predict.

The general problem is well summarized by considering the two extreme ideal cases shown in fig. 3.3, depending on the strength of the expansion compared to the nonlinear terms. We show in panel (a) the initial (isotropic) spectrum. Final spectra are shown in panels b and c: panel b shows the case (i) of negligible expansion and panel c that (ii) of a very strong expansion. At the beginning of (i), the mean field is at 45 0 with the x direction; at the beginning of (ii), it is close to the radial, so that at the end, the strong damping of the B x component due to expansion (see section 3.2.2) will have turned the mean field to 45 0 . Therefore, the mean field is at 45 0 at the end of both cases. However, the spectra E(k x , k y , k z = 0) show very different geometries.

In case (i) (panel (b)), expansion of the wind can be neglected, so that the cascade proceeds in the standard way, namely perpendicular to the mean field, as described earlier in Section 3.1.3. In second (ii) (panel (c)), nonlinear couplings are negligible, and one sees only the linear effect of the expansion, which consists in expanding all structures perpendicular to the radial, corresponding, in Fourier space, to a plain collapse of energy isocontours towards the radial axis.

In conclusion, in these two extreme ideal cases, well identifiable (and different) spectral anisotropies will be found, with respectively the mean field or the radial directions as a symmetry axis. As a corrolary, one can expect in general that the spectra will exhibit intermediate patterns mixing the two symmetry axis ( [START_REF] Saur | Geometry of low-frequency solar wind magnetic turbulence: Evidence for radially aligned Alfvénic fluctuations[END_REF], [START_REF] Grappin | Waves and streams in the expanding solar wind[END_REF], [START_REF] Dong | EVOLUTION OF TURBULENCE IN THE EXPANDING SOLAR WIND, A NUMERICAL STUDY[END_REF], [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF]).

Fluctuations decay: strong versus weak expansion

When expansion alone is active, that is, nonlinear couplings can be neglected, then two important simple cases arise in the framework of a radial flow with constant speed.

The first case occurs when the typical frequency ω of the fluctuation (e.g., for Alfvén waves the frequency given by eq. 2.7) is small compared to the expansion rate U 0 /R. Then the invariants of the problem are: mass, radial momentum, angular momentum and magnetic flux. By taking the EBM equations and assuming that the temporal evolution of the previous invariants only depends on the damping terms due to expansion on the right-hand side of eq. (2.26), we obtain the decay laws corresponding to the previous invariants:

ρ ∝ 1/R 2 (3.36) u x = constant (3.37) u y,z ∝ 1/R (3.38) B x ∝ 1/R 2 (3.39) B y,z ∝ 1/R (3.40)
Another set of decay laws is found when the Alfvén frequency is larger than the expansion rate. This implies u and B/ √ ρ adopt a common (WKB) decay rate, or in other words, that wave-action is conserved for Alfvén waves (see [START_REF] Grappin | Waves and streams in the expanding solar wind[END_REF] for more details on the derivation of these power laws)

u x,y,z B x,y,z / √ ρ ∝ 1/ √ R (3.41)

Cascade rate and turbulent heating

How much heat is necessary in the one-fluid MHD approach to achieve a temperature profile T p ∝ 1/R α ? We reproduce here the analysis given e.g., in [START_REF] Bernard | Evaluation of the turbulent energy cascade rates from the upper inertial range in the solar wind at 1 AU[END_REF] (see also chapter 11 in part III). We take the spatial average of the temperature evolution equation, replace the temporal derivative by the radial derivative, ∂ t → U 0 ∂ R , and as well the damping term /a by its dimensional expression U 0 /R(t). This gives:

U 0 ∂ R T + 2(γ -1)T (U 0 /R) (γ -1)Q ν (3.42)
We define now Q α as the heating necessary to obtain a temperature profile of the form T = 1/R α . Thus we replace in eq. 3.42 U 0 d T /dR = -αU 0 T /R; this gives:

Q α = 4-3α 2 T U 0 /R (3.43)
In the following we will insist on the special value α = 1, which leads to what we will call the critical heating Q c :

Q c = Q 1 = (1/2)T U 0 /R (3.44)
Finally, it is important to note that, to obtain a temperature profile such as T = 1/R α in a given distance range, heating must be continuously equal to eq. 3.43 in this range.

Chapter 4

Solar Wind turbulence

In situ observations

Thanks to the deployment of satellite probes in the heliosphere since the 90s (Ulysses, Helios, ACE, Wind...), solar wind is one of the few examples of astrophysical plasmas in which in-situ measurements can be made. Since the launch of these probes, it has been possible to recover data from solar wind plasma using several types of instruments among which search-coils and mass spectrometers are specially important in the measurement of turbulence in the solar wind. The first ones consist on three coils, each one aligned along three orthogonal axis. Magnetic fluctuations of the plasma induce currents in the coils that are then amplified and pretreated by an on-board analyser. This finally allows to measure the three components of the magnetic field at the position of the satellite.

The mass spectrometer on the other hand consists in a particle detector that measures the kinetic energy of particles in several directions, allowing to recover the distribution function f (v) for different particles species. Once the distribution function is obtained, it is possible to compute its moments and thus to obtain the fluid variables: density, velocity and temperature of the plasma.

There are however some limitations with respect to the measurements that can be done with satellite probes. The range of values to be measured by the on-board instruments is limited. Most importantly, the number of spacecraft is limited, since probes can only provide information on in situ properties of plasma at the point where they are placed. Finally, spacecraft can only sample a limited part of the heliosphere. At the points where the probes are present, the Solar Wind flows radially at super-Alfvénic speed. Since the flow is super-Alfvénic, one can assume that Solar Wind fluctuations are frozen-in and thus transform time into space (or frequency into wavevectors) for fluctuations propagating as fast as the Alfvén speed (i.e. Alfvénic fluctuations). This is known as the Taylor hypothesis and is valid only as long as the velocity of the wind is much larger than the phase velocity of fluctuations. It allows to relate directly the Doppler frequency f (the one seen by the spacecraft) to the (radial) wavenumber and the wind bulk speed, neglecting the frequency f 0 seen in the plasma frame:

f = f 0 + kU sw 2π ≈ kU sw 2π (4.1)
At MHD scales, fluctuations propagate at Alfvén or magnetosonic speeds, which are one order of magnitude below the wind bulk speed, so the Taylor hypothesis is valid. At subion scales, Taylor's hypothesis might be no longer valid for certain modes, such as whistler modes ( [START_REF] Howes | VALIDITY OF THE TAYLOR HYPOTHE-SIS FOR LINEAR KINETIC WAVES IN THE WEAKLY COLLISIONAL SOLAR WIND[END_REF]).

A single spacecraft can only give information about fluctuations propagating along the direction of propagation of the wind, namely the radial direction. In particular, the 1D energy spectrum of magnetic field fluctuations can be recovered from the the autocorrelations function R BB (l) =< B(x + l)B(x) >. Let us define a reference frame where x is the unitary vector in the direction of propagation of the wind, ŷ is perpendicular to x and B 0 and ẑ is the orthogonal to the other two. Then, the autocorrelation function measured by a satellite is the autocorrelation in the x direction,

R BB (l) = R BB ((τ U sw , 0, 0)) (4.2)
where τ is the sampling time. Using the conversion relations stated in section 3 it becomes the reduced spectrum for k x wave-numbers, E(k x ), as defined in eq. 3.32. The reduced energy spectrum in terms of wave-numbers is related to the energy spectrum in terms of frequency by

E(k x )dk x = P (f )df (4.3)
By Taylor's hypothesis, f = k x U sw /2π, and therefore the reduced energy spectrum in terms of the satellite frequency is

P (f sat ) = 2π
Usw E(k x ). If one assumes symmetry of the 3D energy spectrum around the mean magnetic field axis, a single spacecraft can provide the reduced energy spectrums in the parallel and perpendicular directions to the mean magnetic field ( [START_REF] Saur | Geometry of low-frequency solar wind magnetic turbulence: Evidence for radially aligned Alfvénic fluctuations[END_REF]), E and E ⊥ (= E ky in our reference frame ) respectively. To do so, it is enough to project the distance l = τ U sw into the parallel and perpendicular directions to the mean magnetic field l = τ U sw cos(θ) and l ⊥ = τ U sw sin(θ), where θ is the angle between B 0 and U sw . Thus, the measured autocorrelation becomes

R BB (l) = R BB ((l , l ⊥ , 0)) (4.4)
There are other techniques that do not rely on the assumption of gyrotropy of the energy spectrum around the B 0 axis. The so-called k-filtering technique ( [START_REF] Lefeuvre | Local characterization of homogeneous turbulence in a space plasma from simultaneous Measurements of field components at several points in space[END_REF]) uses multi-spacecraft measurements to give the spectral energy density in k space for every satellite frequency ω sat . Data is obtained from ensembles of N satellites forming a polyhedral geometry in space, as is done by the 4 spacecraft of CLUSTER mission that form a tetrahedron. The k-filtering technique is capable of measuring the energy spectra with low error for frequencies between 10 -3 Hz and 10 -1 Hz, at the lower half of the inertial range (see figure 4.1). Nevertheless, even within this range the validity of the results depend on maintaining the same distance between the satellites, the regularity of the tetrahedron and the stationarity of the measured signal ( [START_REF] Sahraoui | Three Dimensional Anisotropic kSpectra of Turbulence at Subproton Scales in the Solar Wind[END_REF]). These constraints have to be imposed to the selection of data, that makes more difficult the generalized application of this technique. Thus, we will rely more on the results obtained by single spacecraft measurements since a higher number of studies with this technique are at our disposal. We will further address this issue in section (4.3.1).

MHD inertial range

In previous sections, we have already mentioned some of the peculiarities that distinguish solar wind turbulence from neutral fluid turbulence. In particular, what stands out is the presence of multiple spectral slopes as can be seen in fig. 4.1. The spectrum shown is the magnetic energy spectrum, which can be measured up to high frequencies, while the kinetic Figure 4.1: Measurements of magnetic energy spectrum at 1AU. From: [START_REF] Kiyani | Dissipation and heating in solar wind turbulence: from the macro to the micro and back again[END_REF]. energy spectrum cannot be measured at frequencies higher than 1 Hz. In this spectrum, one identifies several frequency ranges. At the lower frequencies, the spectrum scales as f -1 . The origin of this spectrum is still an open question (e.g. [START_REF] Verdini | ON THE ORIGIN OF THE 1/f SPECTRUM IN THE SOLAR WIND MAGNETIC FIELD[END_REF]). It is mostly observed in fast solar winds and its frequency range decreases with heliospheric distance ([6], [START_REF] Bruno | The Solar Wind as a Turbulence Laboratory[END_REF]).

Between 10 -4 and 10 -1 Hz the spectral index of the total energy (kinetic + magnetic energy) is close to -5/3 in average (see the distribution of the index in fig. 4.2, left panel). This is generally interpreted to be due to a turbulent quasi-incompressible cascade with constant flux of the total energy (sum of kinetic and magnetic energy) in the framework of MHD. Note however that the spectral index shown is only an average. In reality, it is varying systematically with the wind velocity, passing for the total energy from -1. Coming back to fig. 4.1, at frequencies larger than 1 Hz, around the gyration frequency of protons, the magnetic spectrum becomes steeper, with a scaling close to f -2.8 . It is thought that at these frequencies, linear perturbations in the plasma such as kinetic Alfvén waves (KAWs), cyclotron waves and/or fast magnetosonic waves are present ( [START_REF] Howes | A prescription for the turbulent heating of astrophysical plasmas[END_REF]). Some authors propose that at still higher frequencies, turbulence cascades into the dissipation scales around the electron gyration frequency, where the spectra have been found to decrease exponentially ([2], [START_REF] Alexandrova | Comment on "Evidence of a Cascade and Dissipation of Solar-Wind Turbulence at the Electron Gyroscale[END_REF]). Other works report a power law f -4.1 instead of an exponential decrease ( [START_REF] Sahraoui | Evidence of a Cascade and Dissipation of Solar-Wind Turbulence at the Electron Gyroscale[END_REF], [89]).

Spectral anisotropy

As it has already been explained in sections (3.1.3) and (3.2.1), turbulence properties can vary with wave-vector direction. In the context of turbulence, anisotropy may actually refer to either spectral index anisotropy, power anisotropy or wave-vector anisotropy ( [START_REF] Chen | INTERPRET-ING POWER ANISOTROPY MEASUREMENTS IN PLASMA TURBULENCE[END_REF][START_REF] Chen | Anisotropy of Alfvénic turbulence in the solar wind and numerical simulations[END_REF]). The variation of the spectral indices of the reduced spectrum with the directions of measurement is known as spectral index anisotropy. Power anisotropy is defined as the ratio of the reduced power spectra in two different directions, typically the parallel and one perpendicular radiation to the mean magnetic field, E(k )/E(k ⊥ ). In particular, it provides information about the end of the inertial range and the relative weight of the energy flux in each direction.

Here we are interested by the wave-vector anisotropy (that will be called simply spectral anisotropy in the following). This anisotropy can be measured most easily in direct simulations by plotting energy isocontours of 3D spectra, or related 2D spectra. Similar analysis can be made of autocorrelations of the magnetic and velocity fluctuations. However, such 3D or 2D measurements are not directly possible with single spacecraft data, which produce only temporal records, from which 1D radial records can be extracted via the Taylor hypothesis.

Nevertheless, 2D turbulent structures of the Solar Wind can be recovered from the 1D records, if we assume that turbulence in the solar wind is gyrotropic, namely axisymmetric with respect to the mean magnetic field B 0 . In this case, thanks to the fact that the mean magnetic field is changing direction constantly, one may attach a frame to the mean field direction: in this case one can produce autocorrelation figures which depend on distance, but in different directions with respect to the mean field, i.e., 2D autocorrelation figures.

Two methods can be used for that, which use different definitions of the mean magnetic field, thus leading to different definitions of the fluctuation amplitudes which are to be correlated:

B LOC = 1 ∆t ∆t Bdt (4.5) B GLOB = 1 T T Bdt (4.6)
The local mean field (B LOC ) is computed on the same time interval ∆t as the one on which the correlation is computed. The global mean field (B GLOB ) is computed on a fixed interval of length T , the time intervals over which the correlation is computed being all smaller than T.

Remember from section 4.1 that with the frozen-in hypothesis this is equivalent to com-pute the mean magnetic field at a scale l ≈ U SW 2πf = U SW ∆t 2π . When considering the local magnetic field, magnetic fluctuations are computed at each scale with respect to B LOC . As a consequence, the resulting anisotropy changes with respect to the one we would find when using the global mean magnetic field.

We will consider here only the global method which has two advantages on the local method: (i) it has been used extensively in the solar wind, showing in particular different geometries in the slow and fast winds, as we will see; (ii) it can be directly related to the spectra and hence to the turbulent cascade rate and the resulting heating, which is not possible with the local method [START_REF] William H Matthaeus | Local Anisotropy, Higher Order Statistics, and Turbulence Spectra[END_REF].

The Maltese Cross

Measurements of turbulent anisotropy are made at 1 AU in order to ensure the stationarity of the fluctuations ( [START_REF] William H Matthaeus | Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind[END_REF]). At this distance from the Sun, in situ measurements of turbulent anisotropies have been done using the techniques already introduced in section 4.1: singlespacecraft autocorrelations and multiple-spacecraft k-filtering/wave telescope techniques.

In slow winds turbulent anisotropy has been measured using k-filtering techniques at 1AU between 10 -3 and 1Hz ( [START_REF] Narita | Wave-Vector Dependence of Magnetic-Turbulence Spectra in the Solar Wind[END_REF], [START_REF] Narita | Three-dimensional spatial structures of solar wind turbulence from 10 000-km to 100km scales[END_REF]). In fast winds, measurements have been done at higher frequencies, between 10 -2 and 1Hz ( [START_REF] O W Roberts | Direct Measurement of Anisotropic and Asymmetric Wave Vector Spectrum in Ion-scale Solar Wind Turbulence[END_REF], [START_REF] Owen W Roberts | Three-dimensional density and compressible magnetic structure in solar wind turbulence[END_REF]). For all cases observed using this technique, turbulence was found to be developed in the plane perpendicular to the mean magnetic field.

Single spacecraft techniques have also been used to measure turbulence anisotropy in fast and slow winds but, in contrast with k-filtering techniques, they need to assume axial symmetry around B 0 . Their advantage with respect to multi-spacecraft techniques is that they are able to measure anisotropy at lower frequencies. It is at these low frequencies, between 10 -4 and 10 -3 Hz, that the autocorrelations of the magnetic fluctuations show an structure known as the Maltese cross ( [START_REF] William H Matthaeus | Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind[END_REF]).

The interpretation of Matthaeus et al. was that the observed autocorrelation is composed of two components: one aligned with the mean magnetic field and the other in the directions perpendicular to B 0 . Since autocorrelations are computed in real space, the Fourier space counterparts of these components would correspond to spectral isocontours mainly perpendicular to B 0 , known as the "2D" component of the Maltese cross; and the other component would correspond to spectral isocontours extended along the mean magnetic field axis, the "slab" component. Later analysis of the power anisotropy showed that the 2D component contains about 80% of the energy ( [START_REF] Bieber | Dominant two-dimensional solar wind turbulence with implications for cosmic ray transport[END_REF]).

In 2005, [START_REF] Dasso | Anisotropy in Fast and Slow Solar Wind Fluctuations[END_REF] found distinct structures for slow winds U SW < 400km s -1 and fast winds U SW > 500km s -1 : the 2D geometry was observed in slow winds and the SLAB geometry in fast winds. This observation was made in the middle of the inertial range, from 10 -3 Hz (entering Matthaeus et al's interval of measurements) and 10 -2 Hz. In the same range of frequencies as that of Dasso's measurements, it was observed that for intermediate wind speeds, the wave-vector anisotropy was also between the "2D" and the Slab geometries, suggesting a continuous transition from one to the other as wind speed increases ( [START_REF] James M Weygand | Correlation and Taylor scale variability in the interplanetary magnetic field fluctuations as a function of solar wind speed[END_REF]).

This change of the spectral anisotropy in fast winds is beyond the standard description of MHD turbulent cascade, which is generally believed to be mainly developed in the plane perpendicular to the mean magnetic field. Nevertheless at still higher frequencies, between 10 -2 and 10 -1 Hz, the "2D" pattern has been found to be dominant ( [START_REF] Owen W Roberts | Three-dimensional density and compressible magnetic structure in solar wind turbulence[END_REF]).

Thus, in fast winds, anisotropy varies along the inertial range from an slab geometry at low frequencies to a 2D one in the high frequency range of the inertial range. On the other Figure 4.3: Left: autocorrelations of the magnetic field fluctuations computed at 1AU without distinction between fast and slow winds that form the so called Maltese cross (from: [START_REF] William H Matthaeus | Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind[END_REF]). Center: same kind of measurements done at higher frequencies and setting appart slow (up) and fast wind (down) data (from: [START_REF] Dasso | Anisotropy in Fast and Slow Solar Wind Fluctuations[END_REF]). Right: reduced two dimensional spectra computed from the integration of the 3D magnetic spectra in the azimuthal direction to B 0 . They have been computed from slow (up) and fast winds (down) at 1AU using multi-spacecraft techniques (from: [START_REF] Owen W Roberts | Three-dimensional density and compressible magnetic structure in solar wind turbulence[END_REF]) hand, slow winds anisotropy is more robust, showing the same 2D geometry perpendicular to B 0 all along the inertial range. A clear correlation is visible between velocity and proton temperature variations. It is not absurd to think that this correlation begins not too far from the accelerating region, and thus is somewhat a signature of the acceleration process. Then, as the plasma flows away, this initial correlation works as an initial condition, that is, the velocity and temperature live their lives independently. In other words, during the advection of plasma volumes by the solar wind, the initial correlation should in principle be progressively destroyed due to the plasma cooling during transport, and, to a lesser degree, to a residual additional acceleration of the wind -which act in opposite ways. Can we separate the two effects?

Let us try, ignoring for simplicity any systematic speed variations with distance. We multiply the temperature by R α , and measure the correlation between this compensated temperature and bulk velocity. We find in this way that the correlation rises from 71% with no temperature compensation (α = 0) up to the much higher value 91% when compensating by α within [0.8,1.3]. If we use α outside this interval, the correlation coefficient decreases.

However, although heating, (that is, a cooling slower than adiabatic), is revealed in this simple way, the interval of optimal α values is large, too large to conclude with certainty that substantial heating is at work (no heating corresponds to α = 4/3 = 1.33). The problem is that the systematic radial dependence is hidden by the much larger (correlated) variations of bulk velocity and proton temperature. We give in fig. 4.4b-c the scatter plots of velocity vs compensated temperature R α T p , with α = 0 and with α = 1.

Note that the turbulent amplitude as well is highly correlated with proton temperature, as shown in fig. 4.5. Hence, the radial variation of the turbulent amplitude is also largely masked by the variation due to the wind source.

Temperature gradient: a long series of studies

The methodology used to measure the index α of T p ∝ R α has changed and evolved with years. The initial work [START_REF] Schwenn | Two States of the Solar Wind at the Time of Solar Activity Minimum -Part Two -Radial Gradients of Plasma Parameters in Fast and Slow Streams[END_REF] took the data from radial line-up configurations of the two Helios satellites. They conclude to a small departure from the adiabatic prediction in the case of slow solar winds (V sw < 400km/s), while fast winds (V sw > 600km/s) show a flatter temperature profile.

In the following years, [START_REF] Marsch | Solar wind protons: Three-dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU[END_REF] improved the statistics of the previous method by measuring the bulk velocity distributions at different heliospheric distances between 0.3 and 1AU for fast and slow winds. Temperature was then computed for bulk velocity bins of 100 km/s. As a consequence, the measured radial dependance of proton temperature for slow winds becomes R -1 , and flatter temperature profiles obtain for fast winds. Note that this work neglected the possible mixing of slow and fast wind populations due to the acceleration of the winds.

Measuring of the velocity gradients between 0.3 and 1AU done by [START_REF] Arya | Estimates of solar wind velocity gradients between 0.3 and 1 AU based on velocity probability distributions from HELIOS 1 at perihelion and aphelion[END_REF] was implemented into Marsch's methodology by Totten et al. [START_REF] Lynn Totten | An empirical determination of the polytropic index for the free-streaming solar wind using Helios 1 data[END_REF]. This resulted in fast and slow wind proton temperatures decreasing at the same rate, namely as T p ∝ R -0.9±0.1 .

Other authors have continued to measure the temperature gradient in the inner heliosphere assuming constant radial speed. Such is the case of [START_REF] Hellinger | Proton thermal energetics in the solar wind: Helios reloaded[END_REF], who finds in wind populations with speed below 300 km s -1 that T p ∝ R -0.58 . For winds with speed over 600 km s -1 , it was found that T p ∝ R -0.74 [START_REF] Hellinger | Heating and cooling of protons in the fast solar wind between 0.3 and 1 AU: Helios revisited[END_REF]. The difference with respect to Totten et al's measurements is larger for the case of slow winds, since they suffer a stronger acceleration in the inner heliosphere.

More recently, [START_REF] Leslie J Lamarche | Proton temperature change with heliocentric distance from 0.3 to 1 AU according to relative temperatures[END_REF] has found a non-monotonic variation of α with the temperature associated to the closest distance to the Sun. The hotter and cooler winds shown presented an evolution similar to Totten's, R -0.8 , while intermediate temperatures reached a faster decrease, R -1 . The difference with respect to Totten et al. consisted in dividing the temperature distributions into different classes instead of working with the velocity distributions. The possible mixing of temperature classes was also accounted for in this work.

Thus, different methodologies used in the past have led to different results in the mea-surement of the evolution of proton temperature with heliospheric distance. In this thesis, we choose for the sake of simplicity to take as a reference the simple profile T p 1/R, compatible with Totten et al's result T p 1/R 0.9±0.1 [START_REF] Lynn Totten | An empirical determination of the polytropic index for the free-streaming solar wind using Helios 1 data[END_REF].

Possible explanations of the slow proton cooling in the inner heliosphere

As discussed in [START_REF] Scudder | Radial Variation of the Solar Wind Proton Temperature: Heat Flow or Addition[END_REF] there are four possible ways to slow down the proton cooling below 1AU: (i) the compression of the plasma volume, (ii) the deposition of coherent energy at dissipative scales, (iii) the divergence of the heat flux and (iv) the energy exchange with other plasma particles. The contribution of each of these elements to the variation of proton temperature can vary along the heliosphere. For instance, from 10 AU onwards, the temperature gradient observed by [START_REF] Richardson | Radial evolution of the solar wind from IMP 8 to Voyager 2[END_REF] is strongly affected by the presence of pick-up ions [START_REF] Isenberg | Turbulent Heating of the Distant Solar Wind by Interstellar Pickup Protons[END_REF][74], ionized particles by the passage of the solar system through the Galaxy's interstellar medium. Other works also suggest a certain variability of the contribution of the turbulent cascade to solar winds energy between 3 and 4 AU [START_REF] Marino | Heating the Solar Wind by a Magnetohydrodynamic Turbulent Energy Cascade[END_REF].

A significant part of the research done by past and future missions such as Parker Probe and Solar Orbiter is dedicated to the proton temperature evolution within the inner heliosphere, below 1AU. Within this region, the observed proton temperature results from the combination of the four mechanisms mentioned in the previous paragraph. Their contributions however appear to be unequal: (i) Solar wind plasma is compressed mainly by the interaction of different wind streams.

It has been argued that the compression of streams contributes no more than a 15% to proton heating [START_REF] Burlaga | Solar wind temperature and speed[END_REF]. [START_REF] Miyake | Stream interaction as a heat source in the solar wind[END_REF] showed that this heating source is localized at the compression regions between interacting jets (around 0.6 AU in their study). Since these compression regions are not abundant and very localized between 0.3 and 1 AU [START_REF] Richardson | Solar wind stream interaction regions throughout the heliosphere[END_REF], and their estimated contribution to the total heating is low, it can be neglected as a possible source of the low temperature decrease.

(ii) The energy spectrum has been measured by several authors within the inner heliosphere [6] [59]. As reported in the following sections (4.4.2 and 4.4.3), several methods for measuring, directly or indirectly the energy cascade rate associated with the turbulent fluctuations in the inertial range of the solar wind indicate that it provides a source term for heating large enough to explain the observed slow proton cooling, at least at 1 AU.

(iii) The possibility that the divergence of heat flux is the main contributor to the proton temperature gradient is defended by Scudder [START_REF] Scudder | Radial Variation of the Solar Wind Proton Temperature: Heat Flow or Addition[END_REF]. The measurement of the proton heat flux from Helios data via the third order cumulant (skewness) of the proton VDF is not unique. The reason is that Helios particle detector selected particles in terms of their energy per electric charge. Since alpha particles are positive charges with a mean drift with respect to the proton VDF, the resulting VDF detected by Helios shows a double peak, accounting for the protons and alphas contributions. The method used to separate the VDFs was developed by Marsch et al. [START_REF] Marsch | Solar wind protons: Three-dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU[END_REF] and has been the most widely used to evaluate proton heat flux. With Marsch's method [START_REF] Marsch | Solar wind protons: Three-dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU[END_REF], the contribution of the divergence of the proton heat flux can be neglected. On the contrary, Scudder shows in his paper [START_REF] Scudder | Radial Variation of the Solar Wind Proton Temperature: Heat Flow or Addition[END_REF] that with his method to compute the skewness of the proton VDF measured by Helios, the divergence of the proton heat flux becomes the main contributor to proton heating of the solar wind. As it is acknowledged by Scudder [START_REF] Scudder | Radial Variation of the Solar Wind Proton Temperature: Heat Flow or Addition[END_REF], the discrepancies between the results based on his method and that of Marsch will be solved only by the data gathered by the 3D mass-resolved spectrum of Parker Solar Probe mission, in which the VDF of protons will be measured separately from any other positive ion.

(iv) The exchange of energy between different particles species can also alter the evolution of proton temperature in the heliosphere. Part of the contribution of (i), (ii) and (iii) to proton temperature can be directed to the heating of electrons and the heavy ions in the solar wind plasma [START_REF] Neugebauer | Observations of solar-wind helium[END_REF], [START_REF] Anders | Abundances of the elements: Meteoritic and solar[END_REF]. The opposite, electrons and heavy ions heating protons, is also possible.

Considering the temperature and density of protons and alpha particles (see section 1.1.1), [START_REF] Bernard | Evaluation of the turbulent energy cascade rates from the upper inertial range in the solar wind at 1 AU[END_REF] estimates that alpha particles represent, at most, 12% to 16% of proton turbulent heating and refers to previous bibliography that sustains this idea. The dependency of the temperature and thermal velocity ratio of heavy ions to protons with the collisional age for each species1 also suggests that the thermalization of heavy ions (and therefore energy exchange by Coulomb collision with protons) occurs mostly at the lower Corona [START_REF] Maruca | Collisional Thermalization of Hydrogen and Helium in Solar-Wind Plasma[END_REF], [START_REF] Tracy | Thermalization of Heavy Ions in the Solar Wind[END_REF]. Until a direct measurement of Helium to proton temperature ratio distribution at distances below 0.3 AU is done (by Solar Orbiter or Parker Solar Probe missions), we consider that the role of heavy ions in the process of turbulent heating between 0.2 and 1AU is negligible.

For what concerns the electrons, [START_REF] Stawarz | THE TURBULENT CASCADE AND PROTON HEATING IN THE SOLAR WIND AT 1 AU[END_REF] argues that the measured excess of turbulent heating goes to other ions and electrons (at most 1/3 of the total turbulent heating [START_REF] Robert J Leamon | Observational constraints on the dynamics of the interplanetary magnetic field dissipation range[END_REF][21] [START_REF] Breech | Electron and proton heating by solar wind turbulence[END_REF]). In contrast, recent measurements by [START_REF] Štěpán Štverák | Electron energetics in the expanding solar wind via Helios observations[END_REF] show that the high electron heat flux for both fast and slow winds is enough to explain the heating of this species. Energy exchange by Coulomb collisions between electrons and protons can be neglected in the inner heliosphere, according to this work. The observation of a different set of solar wind data and the use of temperature measurements done by [START_REF] Hellinger | Heating and cooling of protons in the fast solar wind between 0.3 and 1 AU: Helios revisited[END_REF][START_REF] Hellinger | Proton thermal energetics in the solar wind: Helios reloaded[END_REF] might explain part of the differences with respect to the authors defending that part of the turbulent heating is needed to heat the electrons and other species.

Hence, the role of electrons cannot be easily neglected in the energy balance equation of the solar wind: either electrons might act as a source of heating for protons or as a sink for part of the turbulent heating. Due to the lack of this ingredient in our numerical model, we will not be able to establish that turbulent heating is the only one responsible for the proton temperature gradient. In our work, we make the assumption that all turbulent heating contributes only to the heating of protons. In future works, the contribution of electrons to proton heating might be implemented into our work thanks to a method based on the papers of [START_REF] Cranmer | Empirical Constraints on Proton and Electron Heating in the Fast Solar Wind[END_REF] and [START_REF] Breech | Electron and proton heating by solar wind turbulence[END_REF].

Measuring turbulent heating via second order moments

A first idea to measure directly the turbulent heating possibly at the origin of the flat proton temperature profile is to use the Kolmogorov phenomenological estimate Q K41 (eq. 3.13) to evaluate the turbulent energy flux and thus the turbulent dissipation. This has been done by [START_REF] Bernard | Evaluation of the turbulent energy cascade rates from the upper inertial range in the solar wind at 1 AU[END_REF], who find that the Kolmogorov energy flux is a factor 10 above the critical heating necessary to obtain Totten et al's profile (T p ∝ 1/R 0.9 ), in cold winds (see fig4.7). A factor close to 10 between the simple expression of Q K41 and the true turbulent dissipation rate is also found in direct simulations of MHD turbulence, as noted by [START_REF] Bernard | Evaluation of the turbulent energy cascade rates from the upper inertial range in the solar wind at 1 AU[END_REF]. Thus one may consider that this result is a good indication that turbulent heating is close to critical. When in hot winds, the ratio becomes close to a factor 40, however, the latter result depends on the frequency, due to the flat spectral slopes in fast winds which makes the Kolmogorov's expression to be non constant in the inertial range (the authors made this estimate in the frequency band around f = 0.01 Hz).

We have done the same exercice by using the first 118 of Helios I data, as specified in the second section of [START_REF] Grappin | On the origin of solar wind MHD turbulence -HELIOS data revisited[END_REF]. Thus, we computed the Kolmogorov energy flux as

Q K41 = 2πf U 0 (kE V (k) + kE M (k)) 3/2 (4.7)
where f = 4 • 10 -3 Hz and the energy of magnetic and kinetic fluctuations is

δu 2 ≈ kE V (k) and δB 2 ≈ kE M (k).
We compared our computation of the Kolmogorov energy flux to the critical dissipation for 0.3 ≤ R ≤ 1 AU (see fig. 4.8), while the previous work considered only data at 1 AU. If we find that the result doesn't change much during the transport between 0.3 and 1 AU, this will be a good indication that the turbulent heating hypothesis is at the origin of the proton temperature profile T p = 1/R 0.9 .

It is seen that in cold winds, the Kolmogorov energy flux is at all distances a factor ten above the critical heating, in the frequency interval [5 10 -4 ,5 10 -3 ] Hz. With the same data used to obtain the curves in fig. 4.8, we have computed averages during four distance intervals. Figure 4.8 shows that at all distances the period with minimum solar activity and cold winds is comparable to the period studied by Vasquez et al 2007 at 1 AU (note however that the latter period is between 1998 and 2002, which includes a peak of solar activity).

Figure 4.9 shows values which ratio Q K41 /Q c grow quasi-linearly in the log-log plots with frequency for hot winds. These show spectra systematically flatter than for cold winds (see right panel of figure 4.2), that is, flatter than the k -5/3 of Kolmogorov's phenomenology. As a result, in fig. 4.9, the high frequency values show normalized values of the Kolmogorov flux estimate that reach 50 at the high frequency end. The level is not constant, which is not 13) evaluated at frequency f = 4 10 -3 Hz (crosses) and critical heating rate Q c (lines), vs proton temperature. From left to right panels: increasing distance range from 0.3 to 1 AU. This is to be compared with the previous figure 4. [START_REF] Belcher | Large-amplitude Alfvén waves in the interplanetary medium, 2[END_REF] which gives a similar analysis using a much larger dataset at 1 AU. Vertical units: [10 6 J/(kgs)].

contradictory with the possibility that the true normalized heating is distance-independent, as the Kolmogorov estimate is clearly not correct in hot flows.

Measuring turbulent heating via third order moments

Another, more direct method to check the ratio between turbulent heating and critical heating rate (here Q , corresponding to an index α = 0.9 for the temperature gradient), consists in measuring third order moments using the Politano-Pouquet expression for the energy flux in the incompressible limit of MHD ( [START_REF] Politano | Dynamical length scales for turbulent magnetized flows[END_REF]), using several different hypothesis to reconstruct the 3D moments from the 1D measurements. Fig. 4.10 shows the ratio between the heating estimate via 3-third order moments Q 3 and the critical heating Q , vs the product U T p . It is seen that for most cases except perhaps for the coldest-slowest winds, heating is sufficient to generate the temperature profile T p ∝ R -0.9 observed by Totten et al [START_REF] Lynn Totten | An empirical determination of the polytropic index for the free-streaming solar wind using Helios 1 data[END_REF]. This has been shown to be true also for expressions of the energy flux valid for compressible isothermal MHD( [START_REF] Carbone | Scaling Laws of Turbulence and Heating of Fast Solar Wind: The Role of Density Fluctuations[END_REF][START_REF] Banerjee | Exact relation with two-point correlation functions and phenomenological approach for compressible magnetohydrodynamic turbulence[END_REF]): the measurement of the energy cascade rate using these laws is equal or higher than for the Politano at 1AU for slow and fast winds. Turbulent heating is computed via third-order moments laws that have assumed an isotropic hydrodynamic turbulence (green triangles), isotropic MHD turbulence (blue squares) or anisotropic MHD turbulence (red circles). From [START_REF] Stawarz | THE TURBULENT CASCADE AND PROTON HEATING IN THE SOLAR WIND AT 1 AU[END_REF].

Chapter 5

Plan of this thesis

We consider successively two issues related to solar wind turbulence: (i) finding the origin of spectral anisotropy (in Part II) and (ii) understanding turbulent heating (in Part III). This resembles a challenge in that observational measurements (eg, [START_REF] Lynn Totten | An empirical determination of the polytropic index for the free-streaming solar wind using Helios 1 data[END_REF]) support a single proton temperature gradient, pointing to a single turbulent heating regime in both slow and fast winds, while the spectral anisotropy is known to differ in fast and slow streams, which may on the contrary be interpreted as pointing to different turbulent regimes.

Spectral anisotropy: understanding the Maltese Cross

In part II we study the generation of the spectral anisotropy in both the slow and fast wind. As we have seen in section 3.2.1, spectral anisotropy should provide a test for the idea that turbulence and expansion are in competition in the Solar Wind. In this competition, turbulence forces energy to cascade in (wave vector) directions perpendicular to the mean field, and expansion forces wave vectors perpendicular to the radial directions to mimic an inverse cascade that actually is a simple linear consequence of expansion. Such a direction of research has been advocated a long time ago in the first 1D and 2D works by [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF] and [START_REF] Grappin | Waves and streams in the expanding solar wind[END_REF], and more recently by [START_REF] Dong | EVOLUTION OF TURBULENCE IN THE EXPANDING SOLAR WIND, A NUMERICAL STUDY[END_REF].

The work by [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF] (VG16) provides our starting point. In VG16, one shows using the EBM equations that initial conditions naturally evolve towards distinct spectral shapes, reproducing convincingly the two branches of the maltese cross ( [START_REF] William H Matthaeus | Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind[END_REF]). First, when the initial spectrum (at 0.2 AU) has its wave vectors lying mainly in the plan perpendicular to the mean field, the perpendicular cascade proceeds, with the spectrum remaining perpendicular to the mean field while it rotates. This provides the first branch of the maltese cross. Second, when the initial spectrum is isotropically distributed in all directions of Fourier space, the evolution proceeds differently, with the wave-vectors being progressively aligned in the radial direction. The latter structure may be interpreted as showing a dominant effect of the expansion (leading to a collapse of wave-vectors parallel to the radial direction), while the previous structure may be interpreted as showing a standard cascade with dominant turbulent couplings.

We will denote here the two final relaxed regimes respectively as "2D" or gyrotropic, corresponding to the structure observed in slow winds by Dasso et al 2005, and the radial parallel regime as "R-Slab", corresponding to the "Slab" structures observed in fast winds. The reason of the change in terminology from "SLAB" to "R-SLAB" is that the Slab structure, being found by assuming axisymmetry around the mean magnetic field, was believed to be indeed axisymmetric around the mean field, while in the EBM simulations they are found to be axisymmetric around the radial direction. There is however no contradiction, since the sampling method automatically transforms the latter into the former. In the VG16 simulations, to obtain a 2D structure at 1 AU, one must start at 0.2 AU with a 2D structure, while, to obtain an R-SLAB structure at 1 AU, one has to start with an isotropic structure at 0.2 AU.

Our goal here is to see how the previous picture changes when we include cross helicity. The motivation to include cross-helicity is that it characterizes fast winds ( [START_REF] Belcher | Large-amplitude Alfvén waves in the interplanetary medium, 2[END_REF]). This leads possibly to purely linear waves, not to well-developed spectra, while in fact well-defined power-laws are indeed observed, a paradox already remarked 40 years ago ( [START_REF] Dobrowolny | Fully Developed Anisotropic Hydromagnetic Turbulence in Interplanetary Space[END_REF]).

We first present the initial conditions for the simulations and the method. We then study how the combination of expansion and cross-helicity affects the spectral anisotropy, thus extending the work of VG16 to high cross-helicity simulations.

Turbulent heating in slow and fast winds

Previous attempts to cope numerically with turbulent heating in the solar wind ( [START_REF] Breech | Turbulence transport throughout the heliosphere[END_REF]) have all considered a set of closure hypothesis, which allowed to write closed equations for second order moments as total (kinetic and magnetic) turbulent energy, residual energy (magnetic minus kinetic energy) etc... These unknowns were determined not only by the choice of initial conditions, but also by the choice of a set of parameters that appeared in the closure equations. Our aim here is to compute the heating rate in the plasma volume transported by the wind, solving the primitive MHD equations, thus relying only on the initial conditions of the plasma.

We will begin by considering the simplified case of 1D hydrodynamic simulations with expansion. Then we give a copy of the paper we published in 2018 on turbulent heating of slow winds. We then present the extension of these results to fast winds. To attack turbulent heating, it will be necessary to consider rms Mach numbers close to unity, that is, larger than the one considered in part II, so we will give an account of how spectral properties (in particular anisotropy) are modified due to the increase of the Mach number.

Part II

The Maltese Cross revisited Chapter 6

Parameters and initial conditions

We use here EBM simulations to study the generation of the spectral anisotropy in the distance range [0.2,1] AU, which as seen previously (section 4.3.1) is different in the fast and slow wind at 1 AU. All runs have a resolution of 512 3 grid-points. The different parameters are given in Table 6.1.3

Physical parameters and initial spectra

Expansion parameter, cross helicity, Mach number

We model the low frequency range anisotropy (around f = 6.6 10 -3 Hz) by adopting an expansion parameter = 0.2. The simpler high frequency anisotropy around 10 -1 Hz should in principle be modeled by using a much lower expansion rate. However, since this would be highly demanding in computing time, we consider sufficient to consider zero expansion simulations in this case.

Following the observations of the dependance of normalized cross helicity with wind speed at the solar minimum of the 23rd solar cycle [START_REF] Raffaella | ON THE ORIGIN OF HIGHLY ALFVÉNIC SLOW SOLAR WIND[END_REF], we consider the slow wind to be characterized by a zero initial cross-helicity (as defined in eq. (3.21)), and the fast wind to be characterized by σ 0 c ≈ 0.8. Note that in Dasso et al [START_REF] Dasso | Anisotropy in Fast and Slow Solar Wind Fluctuations[END_REF] the slow wind cross-helicity at 1 AU is mentioned to be 0.5 -0.6. This is not contradictory as cross-helicity is known to increase with heliocentric distance.

Simulations have a low Mach number to minimize the possible effects of compressible fluctuations on the development of the anisotropies and enable to reach reasonably high Reynolds numbers. Higher Mach numbers and their influence on spectral anisotropy will be treated in part III. Since they do not show strong differences with respect to low Mach simulations, they will not be considered here.

Initial spectra

We have chosen three initial spectral forms at 0.2 AU: (i) isotropic spectra (ISO simulations), (ii) spectra perpendicular to the radial directions, that is in practice quasi-perpendicular to the mean magnetic field B 0 (gyrotropic or GYRO simulations) and finally (iii) a composite spectrum we have called the Gyro-Alfvén spectrum (G+A simulations).

The G+A form has been inspired from [START_REF] Oughton | Transport of solar wind fluctuations: A two-component model[END_REF]. It consist in the sum of the GYRO spectrum and a 1D component along the radial direction made of Alfvén modes with maximum crosshelicity, that is, propagating outward from the Sun. The relative amplitude of the 1D component fixes the cross-helicity of the whole spectrum: e.g., the 1D component energy is equal to the GYRO component energy to build a global spectrum with negligible crosshelicity σ 0 c ≈ 0, while with an amplitude of the 1D component equal to 5 times the amplitude of the GYRO component one obtains an initial σ 0 c ≈ 0.8. A 2D slice of the three dimensional spectra, i.e. E 3D (k x , k y , k z = 0) is shown in fig. 6.1. We want kinetic and magnetic energy equipartition with amplitude 1, divergenceless velocity fluctuations, an initial 1D spectrum k -m 0 and to set the initial cross helicity. This is obtained by choosing randomly and independently z + and z -in 3D Fourier space with an energy scaling as k -m 0 +2 4 and then reconstructing the u and b fields from the z + and z -fields. It is also worth noticing that, since the random number generator algorithm that we use is not perfectly random, z + and z -cannot be completely uncorrelated. This implies initial energy equipartition is not perfect and also that the initial cross helicity is slightly different from the value listed in table 6.1.3.

Domain aspect ratio

As explained in the introduction, the transverse sizes of the transported plasma volume L Y and L Z increase linearly with distance while the radial size remains constant. As a consequence, the aspect ratio of the domain L X /L Y decreases with distance as 1/R.

In a first numerical study of spectral anisotropy evolution between 0.2 and 1 AU, made by [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF] (VG16 from now on), the choice had been made of an initial aspect ratio at 0.2 AU of R x = L X /L 0 Y = 5, implying that the final domain aspect ratio L X /L Y is unity. Such a choice was made in part to satisfy the critical balance condition at large scales.

In homogeneous MHD simulations with mean magnetic field, it is generally admitted that the standard strong cascade requires a (constant) aspect ratio with the size L of the domain parallel to the mean field B 0 being larger than the perpendicular sizes L ⊥ by the "critical balance" ratio:

L ⊥ /L = δB/B 0 (6.1)
Satisfying to this condition at the beginning of the computation (0.2 AU) would require L ⊥ /L 1/2, i.e., an initial aspect ratio R x = 2, as the initial mean field is close to the radial direction. Note however that the critical balance condition cannot be satisfied during the whole computation time. Indeed, at the end of the computation, the mean field has rotated so as to be at 45 0 from the radial. There is thus no simple way to satisfy the critical balance condition during the whole transport. =1 always); m 0 : initial index of the reduced spectrum; σ c : initial cross-helicity; : initial value of the expansion parameter; M ach = u rms /c s , with c s the sound speed; B 0 : initial mean magnetic field; µ 0 fixes the initial values of the three diffusive parameters µ, η and κ (see eq. 2.30).

Run Type R x k max k min x m 0 σ 0 c M ach B 0 µ 0 R1 ISO 5
Another viewpoint consists in remarking that the aspect ratio of the domain controls the proportion of available "local" triads (that is, quasi-equilateral triads) at a given time of the simulation: an aspect ratio different from unity decreases the proportion of local triads, thus somewhat inhibiting the standard strong cascade. So, choosing R x = 1 will favor the strong cascade at the beginning of the simulation, but not at the end (due to the expansion), while the reverse is true if R x = 5. So, in general, there is no very strong reason to prefer

R x = 1 to R x = 5.
However, in practice, this is not completely true. The choice of an initial aspect ratio R x can alter the distribution of the initial energy spectrum for R x > 1 if no precaution is taken beforehand. This will be discussed in detail later in the case of an initial isotropic spectrum (see section 9.2).

Chapter 7

Defining spectral properties in EBM simulations

We consider first the procedures used to measure properties of the 3D spectra and then of 1D spectra.

Anisotropy of 3D spectra

We consider in this work the spectral aspect ratio A and spectral angle α with respect to the radial axis as defined in [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF] (VG16 from now on).

We take the same frame of reference than in section 2.2 with the x axis in the direction of propagation of the radial solar wind, pointing anti-sunward, ŷ axis orthogonal to x, coplanar to it and to the mean magnetic field, and ẑ defined to complete the right-hand coordinate system. Turbulent anisotropy will be shown qualitatively with a 2D slice of the 3D spectra, E 3D (k x , k y , k z = 0) at the end of the simulation, which corresponds to 1AU. Remark that for all simulations, the mean magnetic field has an angle of 45 • with respect to the radial at the end of the simulation.

Let k maj be the vector from the center of one isocontour to the most distanced point from the center and k min the vector from the center of the same isocontour to the closest point from the center. We then define the aspect ratio of the isocontour as

A = |k min | |k maj | (7.1)
and the characteristic angle of the main axis with respect the radial direction as

α = | arctan( k maj,y k maj,x )| (7.2)
The anisotropy index A gives information about the shape of the isocontours, ranging from 0 to 1, that is, from a straight line when A=0 to isotropic isocontours for A=1 and taking an ellipsoid shape for intermediate values. On the other hand, α ranges from 0 • , when the isocontours are completely aligned to the radial, to 90 • , perpendicular to the radial. We expect intermediate values for A and inclinations close to 0 • or 45 • , the orientation of the mean magnetic field axis at 1AU. We remark that when the isocontours do not have a clear ellipsoidal shape, the method used to define α is not able to accurately recover the major symmetry axis of the isocontour, ). The solid straight line corresponds to the largest vector in modulus from the center to the isocontour of level 10 -8 , k maj , while the dotted straight line is the smaller vector in modulus from the center to the same isocontour, k min . The mean magnetic field axis (not showed) forms an angle of 45 • with respect to the radial at the end of all simulations.

as it can be seen in figure 7.1. Therefore, in order to complement the information provided by α, we will also show the 2D slice of the 3D spectra along with the other information about turbulence anisotropy. We also signal two particularities that appear throughout several figures showing A and α. On the one hand, for k maj < 2 the isocontours just have a few points due to resolution limitations. As a consequence, measurements of A and α below this threshold are not trustworthy. On the other hand, for simulations with aspect ratio of 5, the isocontour tend to reach the limit of the numerical box at about k maj ≈ 10. Beyond this distance the periodic boundary conditions at the limit of the numerical box open the isocontours, and the computed A and α parameters loose all physical meaning for those isocontours.

1D spectral slopes

We now consider defining the reduced energy spectra in the radial direction and in the ẑ direction (always perpendicular to the radial and to B 0 by construction), where the cascade is expected to develop. An example is shown in figure 7.2.

The time of the simulation is divided uniformly into into five intervals. Then, the reduced spectra have been averaged on each of these intervals. Since we are interested on turbulent anisotropy at 1AU, i.e., the end of the simulations, we focus on analyzing the reduced of the last time interval, and use the others to check the stationarity of the inertial range. Note that during the first non-linear times, the system has not yet entered the turbulent regime, so the first interval is not representative of a developed turbulent cascade. The difference of amplitude between spectra shown in figure 7.2 does not represent the real decrease of the amplitude of the fluctuations. For each averaged spectrum, the amplitude has been multiplied by a constant to obtain a more or less constant separation between the spectra.

Due to the dissipation of energy, the amplitude of the fluctuations diminish with time, so the lower spectra correspond to the last intervals of the simulation. To better perceive the spectral slope, all reduced spectra have been compensated by k 5/3 . In the same spirit of clarity for the spectral slope, we have normalized the wave-numbers in the abscissa axis by an estimation of the dissipative wave-number already defined in equation (3.12),

k d = ( ν k k 2 ( û2 2 + b2 2 ) ν 3 ) 1/4 (7.3)
This normalization gives as a result the dissipative part of the spectra for the normalized wave-numbers larger than 1 and the inertial range for values between zero and 1. A visual inspection in figure 7.2 show that this normalization is not always optimal: part of the dissipation ranges may lay in a range of normalized wave-numbers below of 1. Nevertheless, this normalization gives a qualitative indication of the beginning of the inertial range that is useful to compare the variation in time of the inertial range length. Chapter 8

Results

In the following, except explicitly mentioned, we fix the initial domain aspect ratio to be most of the time R x = 1, and the initial spectral index to be m 0 = -1. We have checked that after an initial phase, the spectra don't depend on the choice of these initial parameters. In the zero cross-helicity case, spectral anisotropy bas been analyzed by VG16. Our results will be compared with those of VG16 in the discussion. We consider successively starting from an isotropic spectrum (ISO), a spectrum with main axis perpendicular to the mean magnetic field (GYRO), and a spectrum combining a GYRO symmetry with a radial component, namely the Gyro-Alfvén spectrum (G+A).

Initially ISO spectrum: expansion and cross helicity effects

We study in this section how an isotropic spectrum transforms with time/distance, depending on the initial cross helicity σ 0 c and expansion rate . We compare the evolution with increasing values of the initial expansion rate and with low and large initial cross helicity σ 0 c . Then we consider the evolution of initially gyrotropic and Gyro-Alfvén spectra, with fixed non zero expansion rate = 0.2.

Varying expansion at zero cross-helicity

Increasing the initial expansion rate corresponds physically to focusing on larger scales in the solar wind turbulent spectrum. When doing so, one expects to observe growing departure from the known properties of the standard turbulent cascade perpendicular to the mean field ( [START_REF] Dong | EVOLUTION OF TURBULENCE IN THE EXPANDING SOLAR WIND, A NUMERICAL STUDY[END_REF], VG16).

The evolution starting from the ISO symmetry with zero cross-helicity and initial spectral index m 0 = -1 is shown in fig. 8.1. The figure shows three runs with growing initial expansion rates = 0, 0.2, 0.4 respectively in rows one, two and three (respectively runs R12, R2, R5). The first three rows show the 1D spectra (first two columns, vs k x and k z , with normalization by the dissipation wavenumber k d ), and the last column shows energy isocontours in 2D slices at k z = 0 of the final 3D energy spectrum. The properties of the energy isocontours are summarized in the bottom row, which shows their aspect ratio A(k) and the angle α(k) between the radial and the main axis: the solid line shows the case = 0, the dotted line = 0.2, the dashed line = 0.4. With zero expansion, in the interval 3 ≤ k ≤ 20 (which encompasses the inertial range), one sees a complete loss of the initial isotropy: one finds a perfect alignement of the energy spectral isocontours with the perpendicular direction to the mean magnetic field. With = 0.2 or 0.4, one sees a partial loss of isotropy: the isocontour's main axis adopt an angle α in between the radial (0 • ) and the mean field perpendicular (45 • ), in average around 30 • . This angle of 30 • indicates a regime not far from the gyrotropic, which is in complete contrast with the result of VG16 (Run C, with nearly the same initial conditions than R6), that passed from an isotropic initial spectrum to energy isocontours aligned with the radial at the end of the simulation. The cause of this discrepancy will be discussed in section 9.2.

Varying expansion at large cross-helicity

When initial cross-helicity is large, σ 0 c = 0.8, the next fig. 8.2 shows a significantly different result: while with zero expansion the angle α oscillates again around 45 • , the non zero expansion rates lead to α 15 • (for = 0.2) and α 10 • (for = 0.4). In summary, increasing σ 0 c increases the basic effect of expansion, turning more the main cascade axis from the direction perpendicular to the mean field to the radial axis.

Systematic comparison of zero vs large cross helicity

ISO

As it can be seen in figure 8.3, the value of the anisotropic index A is independent of σ 0 c , giving an almost constant value of A = 0.4 for simulations R2 and R4. The angle α of energy isocontours however, shows a dependence on σ 0 c . Isocontours are more aligned with the radial (α ≈ 15 • ) with high initial cross-helicity, while when σ 0 c = 0, one gets α ≈ 30 • . Simulation R2 (σ 0 c = 0) maintains low cross helicity values, |σ c | < 0.1 during the whole simulation. On the other hand, for simulation R4 (σ 0 c = 0.8) cross helicity decreases slowly, until it reaches σ c = 0.7 at the end of the simulation. The evolution of cross helicity is qualitatively the same for GYRO and G+A simulations: for the cases with high initial cross helicity, σ c decreases slowly and for low initial cross helicity, σ c remains low during the whole simulation.

All reduced spectra in figure 8.3 show a developed inertial range with a 5/3 slope, although its length depends on the initial cross-helicity, being larger for zero cross-helicity. While for R2 the reduced spectra in the x and z directions have an inertial range from 0.015 to 0.06 k x /k d and from 0.05 to 0.2 k z /k d respectively, the corresponding inertial ranges when the initial cross-helicity is σ 0 c ≈ 0.8 spans from 0.025 to 0.05 k x /k d and from 0.05 to 0.15 k z /k d for the reduced spectra in the x and z directions.

GYRO

For GYRO runs, the effect of cross-helicity in turbulent anisotropy is not as remarkable as for ISO simulations. In figure 8.4, runs R7 (σ 0 c = 0) and R9 (σ 0 c ≈ 0.8) show a small difference in the isocontour aspect ratio A: for σ 0 c = 0, it takes values closer to 0.4 for low k maj and slowly increases up to 0.55, whereas it stays at 0.55 for σ 0 c ≈ 0.8. This is similar to the ISO case (fig. 8.3). The angle α is also changed by cross helicity. For run R7, α takes values between 30 • and 35 • with a peak of α = 75 • for k maj = 4.

The reduced spectra of figure 8.4 show a developed spectrum with a 5/3 for x and z direction when the initial cross-helicity is zero. The inertial ranges spans from 0.015 to 0.04 k x /k d in the x direction and 0.04 to 0.1 k z /k d for the z direction. On the other hand, for the high σ 0 c , the reduced spectra in the x direction are steeper, with m = -2.

Gyro-Alfvén model

For Gyro-Alfvén runs (see figure 8.5 ), high cross-helicity simulations become strongly aligned with the radial. For the σ 0 c = 0 run, R9, the spectral aspect ratio A ≈ 0.6, while the runs with σ 0 c ≈ 0.8, the spectra become very flat, with A ≈ 0.4. For σ 0 c = 0, α ≈ 30 • , while for large σ 0 c , it lies between 0 • and 10 • . Reduced spectra in the x direction show slopes steeper than 5/3, k -2 for R9 and k -3 for R10, low and high cross-helicity respectively. In the z direction, reduced spectra are closer to k -5/3 , with inertial ranges from 0.08 to 0.15 k z /k d for R10 to 0.02 to 0.07 k z /k d for R10.

Chapter 9

Discussion

Summary

In conclusion, we find again, as in VG16, two final stable attractors: the GYRO and the R-SLAB.

The main factors that lead to the R-SLAB geometry, with the main axis closest to the radial are: first expansion, then cross helicity, and last the initial geometry. This is at variance with the VG16 results in which expansion and an ISO initial geometry with zero cross helicity were sufficient.

The case of the initial ISO geometry is a good example (see fig. 9.1): (i) with zero expansion, the ISO spectrum transforms into a gyrotropic spectrum, namely with α = 45 • , whatever the value of cross helicity; (ii) with finite expansion, the main spectral axis starts turning towards the radial, the larger the cross-helicity, the better the alignement with the radial.

Explaining the discrepancy with VG16

As explained in section 8.1.1, the spectral anisotropy obtained for run C in VG16 differs from the one obtained in the mentioned section for R5, despite the apparent similarity in their initial conditions. Recall that the 30 • inclination of the isocontours for run R6 is in Figure 9.1: ISO initial conditions, summarizing how the final spectral main axis varies with and σ 0 c : passing from the GYRO structure perpendicular to the mean field (when = 0) to the R-SLAB (quasi-radial) structure (when both and σ 0 c are large).

complete contrast with the result of VG16 (Run C, with nearly the same initial conditions than R6), that passed from an isotropic initial spectrum to energy isocontours aligned with the radial at the end of the simulation. The cause of this discrepancy will be discussed in section 9.2. .We found that the origin of the difference between our work and VG16 comes from the fact that the ISO initial geometry in VG16 is not truly isotropic. Indeed, in VG16, the initial domain aspect ratio is larger than unity (R x = 5). In that case, there is a large scale range 0.2 ≤ k x ≤ 1 which does not exist in the two other directions k y,z for which the excited range is k y,z ≥ 1. In VG16, this range was initially excited, while in our simulations, it is not, precisely to obtain true isotropy at all scales.

In VG16, the large scales are not at all isotropic in the range 0.2 ≤ k ≤ 1: only radial wave vectors are excited in this range. Since this is a large scale range, it contains much energy, so this may trigger a large parallel component in the whole wave number range after some time. This is indeed what happens, leading to an exaggerate dominance of the radial scales at the end. In the following, we denote by k min x the smallest excited wavenumber in the radial direction (see table 6.1.3).

Our argument is illustrated in fig. 9.2 by three runs: runs R11, R3 and R1. They all have an initial spectral index m 0 = -5/3. Runs R11 and R3 are Isotropic runs already considered here, with respectively R x = 5 and R x = 1, with k min x = 1 (see table 6.1.3). Run R1 has R x = 5, but now with k min x = 0.2, i.e., non isotropic in the range 0.2 ≤ k ≤ 1. One sees that, while runs R11 and R3 exhibit the same main symmetry axis rather far from the radial axis, run R1 main axis is seen to be close to the radial. than the others, actually close to what we have obtained in our previous results by adding a large cross helicity from start to our ISO initial pattern.

Part III

Can the Maltese Cross heat?

In this third part of the thesis, we simulate the turbulent evolution of the plasma and measure the visco-resistive dissipation of turbulence via the turbulent cascade, following the plasma between 0.2 and 1 AU. Measuring the visco-resistive dissipation allows us to bypass our ignorance of the exact mechanism that dissipates turbulent energy in the solar wind. Indeed, inasmuch as the Reynolds number is large enough, the energy flux flowing from large to small dissipative scales leads in principle to a universal measure of the dissipation, independent of the precise value of viscosity and resistivity. Moreover, since equations transmit the amount of energy lost by visco-resistive dissipation to the internal energy, we will then be able to measure the true cooling rate of the plasma.

We first report in the following a study of the heating produced by 1D shocks transported by the wind. The numerical and analytical results obtained with this first model will give us hints to understand the results of the more complex 3D MHD case.

We then consider the 3D MHD in the case of the slow wind, published in an ApJ paper on 2018, and finally the case of the fast wind.

Chapter 10

1D turbulent heating

We consider here the evolution of a 1D spectrum in the expanding (radial) wind, that is, with wavevectors all aligned in a given direction. We consider two cases: wavevectors in the radial direction, and wavevectors in a transverse direction (perpendicular to the radial).

1D HD equations with expansion

The code used in this section, BOX1D from now on, solves a hydrodynamic version of the EBM equations (eqs. 2.26):

∂ t ρ + (ρu) a α = -2ρ a ∂ t u + uu a α + 1 ρ P a α = µ ρ u -αu a ∂ t P + u P a α + γP ( u a α ) = -2γP a + ρκT P = ρT (10.1)
where u = ∂ x u, u = ∂ xx u, α = 0 for wave vectors in the radial direction and α = 1 for wave vectors in a transverse direction. We have used the same notations as for the EBM equations (see section 2.2), including the units used to measure space, time, velocity, temperature and density.

There are three different time-scales in this system: the expansion time, t exp = R 0 /U 0 , the non-linear time, t N L = 1/(ku) and the acoustic time t L = 1/(kc s ). Note that if we assume a Kolmogorov's scaling for the velocity fluctuations, we can write down how the characteristic times scale with wave numbers:

t L ∝ k -1 , t exp ∝ 1, t N L ∝ k -2/3 . (see fig. 10.1).
Consider first the case where expansion dominates all other terms (Section 3.2.2), i.e. t exp < (t L , t N L ), then, for the particular case of a wave in the radial direction: This we will call the non-WKB regime. It can be easily seen that the former relations are obtained by neglecting all terms of equations 10.1 with the exception of the temporal variation and the linear damping terms proportional to /a. If instead the acoustic term dominates in a particular range of scales, t L < (t exp , t N L ), the scaling laws change along; this is the so-called WKB regime which reads:

δρ/ρ ∝ 1 (10.2) δu/c s ∝ R 2/3 (10.3) c s = γT ∝ R -2/3 (10.
δρ/ρ = δu/c s ∝ R 1/3 (10.6) c s ∝ R -2/3 (10.7) δu ∝ R -1/3 (10.8)
These equations can be obtained analytically from the conservation of wave-action1 , an invariant in the WKB regime. Similar equations can be derived for the case of a wave with wavevector perpendicular to the radial. Finally, when non-linear interactions dominate, t N L < (t L , t exp ), the WKB-theory cannot predict the evolution of variables.

Modified Burgers equation and semi-analytical solutions

We now derive a simplified model (a modification of Burgers' equation [START_REF] Burgers | A Mathematical Model Illustrating the Theory of Turbulence[END_REF]) equation with non-linear terms taken from the HD 1D equations at zero pressure, but depending only on the variable u, the amplitude of the velocity perturbation [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF]. In order to do so, we assume the WKB limit. We thus obtain an equation similar to the well known Burgers equation, but with a new linear term that accounts for damping due to expansion:

∂ t u + 1 a α uu = νu -nu a (10.10)
The pair (α, n) is determined by the direction of propagation of the modes: α = 0 and n = 1/3 for radial propagation; α = 1 and n = 5/6 for the case of perpendicular propagation. These parameters are recovered from the WKB solutions of equations 10.1. Neglecting the viscous and nonlinear terms of the equations and imposing the WKB predictions (see eq. 10.8 for the radial case), the values of α and n are obtained.

It is worth noticing that in the radial direction, although the plasma doesn't expand in this direction, linear modes are still affected by the damping caused by expansion. The cause of this counter-intuitive phenomenon is that expansion always acts to damp average quantities as density and sound speed, which has an effect on fluctuations as well.

The WKB term in the modified Burgers equation makes difficult the analytical treatment of equation 10.10. However, a change of variables allows to recast the equations in the form of a standard Burgers' equation. This is done by adopting a new time variable, τ , and amplitude, v:

u = v a n (10.11) t = [1 + (1 -β) τ ] 1/(1-β) -1 (10.12)
with n = α + β. Using the new variables v and τ , equ. 10.10 takes the form of a standard Burgers equation:

∂ t v + v∂ x v = ν∂ 2 x v (10.13)
We can now predict the evolution of the amplitude. To do so, we assume that an inertial range has developed, starting from the largest scales. We then derive from eq. 10.15 an equation for the evolution of the average energy content in the domain, integrating the kinetic energy in the periodic domain of length L = 1/k 0 . Assuming periodic boundary conditions, this reads:

∂ t (1/L) L 0 v 2 2 dx = -ν(1/L) L 0 v 2 dx (10.14)
The nonlinear term doesn't appear, as it is the average of a total derivative, which is zero in the periodic case. We obtain (eg, [START_REF] Belmont | Collisionless Plasmas in Astrophysics[END_REF]), changing notation as

√ < v 2 > → v: ∂ t v 2 2 = -k 0 v 3 (10.15)
The contribution of the dissipative term in eq. 10.14 has been transformed into the term -k 0 v 3 by evaluating the shock width in terms of the viscosity and wave amplitude as l ν/v. The solution for the amplitude in terms of the new variables is:

v = v 0 1 + k 0 v 0 τ (10.16)
We express the solution in terms of the original variables using v 0 = u 0 (10.17)

τ = a 1-β -1 (1 + β) (10.18)
where we recall that a(t) = R/R 0 = 1 + t (eq. 2.23). Thus, the general evolution for the amplitude of the fluctuations in terms of the original variables is

v = u 0 1 + k 0 u 0 ( a 1-β -1 (1+β) ) (10.19)
It is worth to express separately the parallel and the perpendicular cases. For the parallel case, α = 0,

β = n = 1/3 u = v a n = u 0 a 1/3 [1 + k 0 u 0 ( a 2/3 -1 4/3 )] (10.20)
For the perpendicular case, α = 1, β = -1/6, n = 5/6:

u = v a n = u 0 a 5/6 [1 + k 0 u 0 ( a 7/6 -1 5/6 )]
(10.21)

Note that in the limit of a = R/R 0 1 we have u ∼ a ∝ 1/R for the parallel case and u ∼ a 2 ∝ 1/R 2 for the perpendicular one. As expected, the amplitude of the fluctuations decreases faster in the direction transverse to the radial, but it also decreases in the radial direction, thanks to the couplings of magnitudes in the WKB regime.

From this point onwards, we have to present several remarks. First, the amplitude of the fluctuations diminishes in both cases as a power law of the heliospheric distance. Second, from 10.12 it can be seen that the formation of a shock wave, or, what is the same, the beginning of the turbulent cascade, t * , depends on the direction and the value of . In the standard Burger's equation, this critical time can be estimated by τ * ≈ 1/(k 0 v 0 ) = 1. Thus, for the original variables,

t * = [1 + (1 -β) ] 1/(1-β) -1 (10.22)
In the parallel direction, even in the limit of high expansion parameter, 1, the critical time can be estimated by t * ≈ β 1-β . On the contrary, in the perpendicular direction, there is an upper limit for the expansion parameter, beyond which the critical time takes negative values. In those cases, no shock is developed and the system is completely dominated by expansion.

Simulations of shock turbulence with transverse waves

Initial conditions

Numerical simulation will be performed for both radial and perpendicular directions using a monochromatic wave as initial condition, u = 10 -2 cos(x) in real space, or an initial spectrum with scaling in k -1 . In this way, we pretend to observe the effects that expansion applies on a particular mode and what happens in a more realistic situation when several Run Spectrum

k max x/y (α, n) µ 0 RB1 δ(k -1) 0 1 (0, 0) 5 • 10 -4 RB2 k -1 0 16 (0, 0) 5 • 10 -4 RB3 δ(k -1) 100 1 Radial (0, 1/3) 5 • 10 -4 RB4 k -1 100 16 Radial (0, 1/3) 5 • 10 -4 RB5 δ(k -1) 0.5 1 Transverse (1, 5/6) 10 -4 RB6 k -1 0.5 16
Transverse (1, 5/6) 10 -4

Table 10.1: Run column contains the names used for the simulations; Spectrum designates the initial energy spectrum; = (U 0 /R 0 )/(k 0 u rms ) denotes the expansion parameter; k max stands for the maximum wave-number containing energy at the beginning of the simulation;

x/y defines the direction of the wavevector, radial or transverse; the pair (α, n) varies, depending on the wavevector direction: (0, 1/3) for the radial direction, (1, 5/6) for the transverse direction and (0, 0) without expansion; µ 0 is the common initial value of the diffusive parameters (viscosity and conduction).

modes interact with each other. The choice of a k -1 scaling is made by reference to the large-scale 1/f spectrum in the solar wind (see fig. 4.1). For each initial conditions, we will vary the expansion parameter, in order to show three different cases: first the case with zero expansion ( = 0), one where the expansion effects dominate non-linear couplings ( = 100), and a last run with ( = 0.5).

Simulations initial conditions are presented in table 10.1. In all cases, the space resolution is N x = 512 and initial temperature is T=15. Numerical simulations will be performed for both radial (α = 0) and perpendicular directions (α = 1). For each direction, we take two types of initial energy distribution: either a monochromatic wave of the form u = 10 -2 cos(x) or a spectrum with an initial scaling in k -1 . The expansion parameter is taken according to the direction: = 100 for the radial direction, = 0.5 for the perpendicular one and = 0, that is common to both of them. These values of the expansion parameter have been taken so turbulence can develop according to the predictions made in the previous section.

Results

In the figures of this section, we show the evolution in time of the energy spectra compensated by k -2 (the scaling expected for a shock) and the evolution of the energy dissipation of the complete system described by equations (10.1) compared to the right-hand side terms (damping and viscous terms) of equation (10.10). The evolution in time of the energy spectra show the energy spectra at six instants, equally spaced in the simulation, beginning by the thin solid line and ending by the dashed triple-dotted line. For what concerns the dissipation curves, we also compare them to the energy cascade rate predicted by the expression k 0 u 3 /a α , directly obtained from Burgers' equation. We expect that the energy cascade rate equals the dissipation term ν(∂ x u) 2 once turbulence is fully developed.

In figure 10.2 we see the results corresponding to single mode initial conditions. For the non expanding case, turbulence develops an inertial range characterized by a spectral index of 2, as expected. The dissipation terms show that the energy dissipation from the BOX1D equations and ν(∂ x u) 2 match perfectly, since in the absence of expansion this is the only term responsible for energy dissipation. The energy cascade rate converge towards both of them at the peak of dissipation, t * ≈ 100 that corresponds to the shock formation time, which is nicely given also by the nonlinear time which reads t N L ≈ 1/(k 0 u 0 ) = 100.

When the initial expansion parameter is = 100 and modes propagate in the radial direction, we again observe a well developed inertial range in k -2 . However, energy dissipation becomes remarkably different with respect the non-expanding case: during the first half of the simulation the main contribution to the total dissipation is given by the damping term due to expansion (2/3) a u 2 , while during the second half the damping term contributes as much as ν(∂ x u) 2 . It is during this second half that shocks are formed, approximately at t * ≈ 600 (as indicated by the local peak of energy dissipation) and the energy cascade rate k 0 u 3 converges with the viscous dissipative term from this moment onwards. We remark that the shock formation time does not match anymore with the estimation of the non-linear time obtained from Burgers' equation, t N L ≈ 100. Instead, it is closer to t ≈ 555.7, the value given by equation (10.12), which is derived from the burgers-like model (eq. 10.10).

It is also worth mentioning the fluctuations in the energy dissipation appearing in the second half of the simulation, which can also be observed in the energy dissipation of run RB4 (see figure 10.3). In both cases, the energy damping due to expansion predicted by the WKB theory (dotted line in the bottom panels of figures 10.2 and 10.3) match well the dissipation of energy during the first half of the simulation (before we enter the turbulent regime). As it is indicated by equation 10.7, the WKB theory predicts an increase with time of the Mach number as M = δu/c s ∝ R 1/6 , close to the Mach evolution in the simulations (not shown). Hence, compressibility in the simulation increases with time 2 . As a consequence, kinetic energy is no longer a proper approximation for the second order invariant (instead of the sum of kinetic and magnetic energy) and energy exchanges between the kinetic and the internal energy are produced. The fluctuations of the kinetic energy are then amplified when its derivative is computed. Since we are interested in the trend of energy dissipation we have smoothed the energy dissipation curve taking local averages of the function on a four points width.

Seemingly to the two previous cases, when = 0.5 and the mode propagates in the same direction of expansion, the spectral index tends to 2, but in this case the inertial range is about one wave-number decade larger than in the previous case. This can be explained by the choice of a smaller viscosity and the subsequent shift of the dissipative region to higher wave-numbers.

On the other hand, the main contribution to total energy dissipation comes from the damping due to expansion (5/3) a u 2 during the whole simulation. Contrary to the radial case, for which the shock was formed after the time interval during which expansion dominates the dissipation, now the shock is formed during this time interval at t * ≈ 200. From that moment onwards, viscous dissipation and the energy cascade rate prediction converge, although the latter has been modified to account for the expansion of the domain, becoming k 0 u 3 /a. Again, we remark that this does not match the Burgers' equation prediction of t N L ≈ 100, but it is closer to the value given by equation (10.12), that is, t ≈ 123. [START_REF] Anders | Abundances of the elements: Meteoritic and solar[END_REF].

In figure 10.3 we observe that when several modes are initially involved in the development of turbulence, we obtain similar results to the ones obtained with a single mode. There are however three differences. First, spectra is more irregular; however, despite that, when averaged over long times, they conserve a scaling in k -2 . Second, there is no clear peak of viscous dissipation that can point to the shock formation time, t * . Third, the predicted energy cascade rate k 0 u 3 /a α does not converge as well toward the viscous dissipation but their values stay close. Dissipation due to the expansion terms 2n a u 2 behave in the same , line styles vary according to the increase in time, from a thin solid line at t = 0, followed be the dotted and dashed line and so on. The bottom row (b,d,f) shows the evolution of the dissipation terms: total dissipation dE/dt (solid line), viscous dissipation ν(∂ x u) 2 (dotted line), damping due to the expansion 2n a u 2 (dashed line), prediction for turbulent dissipation k 0 u 3 /a α (dot-dashed line). The turbulent dissipation always matches the viscous dissipation term after a nonlinear time (see text). way as their single mode counterparts: when modes propagate in the radial direction, the expansion term contributes the most to the total dissipation only during the first half of the simulation, whereas when the modes propagate in the transverse direction, the expansion term dominates all the time.

Discussion

From the simulations aforementioned we can state that the approximate expression for the cascade rate k 0 u 3 /a α is a good approximation of the viscous dissipation for the monochromatic mode simulations. For runs with a fully excited spectrum at t = 0, the previous expression for the energy cascade rate can still be used as an approximation of the viscous dissipation, but it is less reliable. Let us try to extrapolate the results of this section to have an idea of what we are going to see in EBM. We expect that expansion-driven damping will be the maximum contributor to total energy dissipation, even for low values of the expansion parameter, as it is shown by runs RB5 and RB6. The runs in the present chapter confirm that the effect of expansion on the modes propagating in the perpendicular direction is more important than in the parallel direction. This is because the damping due to expansion dominates the simulation a longer time with transverse waves, although the value of is three orders of magnitude lower than the one used for the run with parallel wave. From the runs with a continuous spectrum as initial condition we can also infer that the time at which turbulence develops, t * will not be observed as a clear maximum of the local dissipation. This could be caused by the development of the cascades for the modes with higher wave-number before the ones with lower k, leading to multiple instants at which viscous dissipation could increase.

For what concerns the predictions done with the Burgers-like model of equation (10.10), we have seen that they provide a good approximation of the time needed to trigger turbulence for the initially monochromatic simulations. We have not been capable to predict this time for the simulations with a continuous spectrum as initial condition, as there is no clear peak for the viscous dissipation.

Introduction

The radial evolution of the proton temperature in the solar wind has not yet been fully explained. From 0.2 to 1 au, the wind expands radially and with a low-varying mean velocity. Observations have shown that the rate at which the proton temperature decreases in this interval does not correspond to an adiabatic cooling of a spherically expanding flow; instead of a temperature evolution given by T R 4 3 µ -, a slower cooling rate is observed, T R µ x -for 0. We know that in quasi-stationary incompressible turbulence, the turbulent energy (i.e., the sum of the kinetic and magnetic energies) cascades along the inertial range without being dissipated. When energy reaches the dissipation scales, it is transformed into heat. If the process is quasi-stationary, the heating rate can be obtained either from the dissipation rate (if the small-scale dissipation process is known) or via the energy cascade rate. The dissipation mechanisms take place at sub-ion scales that are not yet reachable by current measurements; therefore, little can be said about the associated heating. On the contrary, the cascade rate can be computed at much larger (inertial range) scales, and this is the approach commonly used in analyzing solar wind data to obtain the turbulent heating.

To understand if turbulent heating is responsible for the nonadiabatic decrease of temperature in the solar wind, one needs to compute the heating rate required to produce the observed temperature profile. Such heating is obtained by exploiting an argument based on the internal energy equation (e.g., Vasquez et al. 2007; see Section 2.3 below). In particular, it shows that the existence of a radial power law implies a direct relation between the energy cascade rate Q, wind velocity U, proton temperature T p , and heliocentric distance R,

Q k mU T R 1 2 , 1 B p = ( )( ) ( )
the coefficient 1/2 being associated with the scaling T R 1 p µ , which we adopt here as a representative scaling. Equation (1) expresses a balance between two decaying quantities: the energy cascade rate, which is a turbulent quantity, and the proton temperature, which is not.

The most accurate measurements of the cascade rate are obtained by computing third-order moments of the distribution of the magnetic field (Politano & Pouquet 1998). Such measurements have shown that Equation (1) holds at 1 au (Stawarz et al. 2009;Coburn et al. 2012) and at larger distances (Marino et al. 2008), although the precise value of the cascade rate somewhat depends on the hypothesis made on the 3D geometry of the angular spectra (Verdini et al. 2015).

The existence of the 1/R profile between 0.3 and 1 au suggests that the balance in Equation (1) is also realized during this whole distance range, but there is presently no measure of the cascade rate and no proof (either theoretical or numerical) that such an equilibrium is indeed achieved in the inner heliosphere. In fact, previous attempts to verify Equation (1) used models of solar wind turbulence with simplified nonlinear couplings (Tu 1988;Tu & Marsch 1997;Smith et al. 2001;Breech et al. 2009). These models actually managed to reproduce the wind temperature decrease; however, they relied on the choice of free parameters to fit their results to the observations. We aim here to examine at what conditions such a relation between turbulent heating and temperature, as well as the associated temperature profile, can be found using direct numerical simulations. It is indeed most probable that not all turbulent conditions close to the Sun are able to generate such a close adjustment between the turbulent cascade rate and the temperature in the distance range 0.3-1 au, so that the numerical solution to this problem will provide constraints on turbulent properties close to the Sun.

We directly compute the temperature and heating evolution versus distance, adopting successively different initial states of the plasma at the minimum heliocentric distance (here 0.2 au). We use the MHD equations modified by expansion, as given by the expanding box model (EBM; Grappin et al. 1993). We study in this work the case of the slow solar wind, where turbulence is mainly in a 2D geometry (Dasso et al. 2005; Verdini & Grappin 2016). The initial conditions found to lead to Equation (1) between 0.2 and 1 au will be characterized by (i) spectral properties and (ii) global plasma properties such as expansion parameter, Mach number, plasma β, and mean magnetic field angle with radial.

Equations, Control Parameters, Diagnostic Tools, and Initial Conditions

EBM Equations (Ideal)

We give here a short description of the EBM equations The wind is assumed to be radial and have a uniform speed (U 0 =const). The radius R at which the numerical domain is located varies with time τ as

R R U , 2 0 0 t t = + ( ) ( )
where R 0 is the initial position of the box. Space, time, velocity, temperature, and density are measured in the following units:

L 2 , 3 0 p ( ) ( ) t L u 2 , 4 NL 0 0 rms 0 p = ( ) ( ) u , 5 rms 0 ( ) m u k 2 , 6 p B rms 0 2 ( ) ( ) ( ) and
, 7 0 r ( )
where 0 r is the initial average density of the plasma, u rms 0 is the initial rms velocity of the fluctuations, t NL 0 is the initial nonlinear time based on the initial rms velocity, and L 0 is the initial size of the numerical domain perpendicular to the radial direction.

The EBM approach relies on the idea that a simple change of Galilean frame is not sufficient to eliminate the expansion. After such a frame change, the plasma still expands: it is stretched in directions perpendicular to the radial. In other words, a systematic velocity field perpendicular to the mean radial direction remains. To recover the usual theoretical setup where the fluctuating quantities are homogeneous (i.e., have zero average) in the plasma volume, we need to subtract this transverse expansion. This is done by using coordinates comobile with this transverse expansion:

t , 8 t = ( ) x X U a , 9 x 0 t = - ( ) ( ) y Y a t , 1 0 = ( ) ( ) and z Z a t . 1 1 = ( ) ( )
The parameter a L L L L 5

x x y x z = = = is the initial aspect ratio of the domain. The parameter a is defined as the normalized heliospheric distance,

a R t R t 1 , 12 0  = = + ( ) ( )
where da dt  = is the expansion parameter defined as the initial ratio between the characteristic expansion and turnover times in the transverse directions (perpendicular to the radial),

U R k u , 1 3 NL exp 0 0 0 rms  t t = = ( )
with k 0 the minimum wavenumber in the transverse direction. At a given distance R(t), the domain thus has an aspect ratio L L a t a x radial = ^( ) . The EBM equations with dissipation terms omitted (but see Equations ( 20) and (21) below) read

u a 2 , 14 t  r r r ¶ +  = - ( ) ( ) ( )
u u P P P P a . . 2 , 1 5 t  g g ¶ +  +  = - ( ) ( ) ( ) u u u B B P B a . 2 . , 16
t 2   r r ¶ +  +  + - = -  ( ) ( ) ( ) B u B B u B u a . . . , 1 7 t   ¶ +  - +  = -  ( ) ( )
and P T. 1 8 r = ( )
In these equations, ρ is the density, P is the total pressure, B is the magnetic field, and u is the velocity fluctuation u U U e r 0 = -ˆ, U being the total velocity. The pressure equation with 5 3 g = is the perfect gas equation, and T T T i e = = is the proton (and electron) temperature, m p being the proton mass.

The above equations are standard MHD equations with, however, two modifications. First, additional linear terms involving the constant average speed U 0 appear on the righthand side: u u 0, ,

y z  = ( ) and B B B 2 , , x y z  = (
) . Hence, depending on the component, the right-hand side damping term differs, as is well known. Note that in the following, we will sometimes use the dimensional form of these damping terms, namely, with

a t U R t 0   ( ) ( ). Second, a new
expression for the gradients is used, accounting for the increasing lateral stretching of the plasma volume with time/ distance:

a at at 1 , 1 , 1 , 19 x x y z  = ¶ ¶ ¶ (( ) ( ( )) ( ( )) ) ( )
With Ox being along the radial, expansion acts only in the other two directions. All fields x y z , , r ( ), u x y z , , ( ), etc. are then considered to be periodic in all three directions of the domain comobile with the mean expansion. This allows to use a pseudo-spectral method for the spatial scheme, as in the standard Lagrangian approach. The temporal scheme is a third-order Runge-Kutta method.

Defining Visco-resistive Terms and Dissipation Rate

The diffusive terms read

u u u 1 3 , 20 t dis 2 m r ¶ = ~ ~+  ~ | ( )( ( ) ( • )) ( ) B B , 2 1 t dis 2 h ¶ =  | ˜( ) and P T , 2 2
t dis rk ¶ = D | ¯˜( )
where r ˜is the density normalized by its average,

a , 2 3 2 r r r r = = ˜¯( )
and a 1 2 r = ¯is the average density of the plasma. Note that the normalized density r ˜is unity, on average. Note that density enters in two different ways in the dissipative terms: as the normalized density r ˜in the momentum equation and as the average r ¯in the pressure equation. In the homogeneous (nonexpanding) case, the dissipative terms given by Equations ( 20) and ( 21) are identical to the standard ones in MHD. In the expanding case, these dissipation terms have two specific differences. First, the ∇ operator components are comobile derivatives, without the anisotropic prefactors a a a

 = ¶ ¶ ¶ ˜( ) ( ) 1 , 1 , 1 x ( ): , , . 24 x y z
This allows the dissipation to be isotropic in comobile coordinates, that is, to use the available Fourier domain most efficiently. Second, we impose that viscosity, resistivity, and conductivity κ decrease as time/distance increases (see Equation ( 12)):

a. 2 5 0 k m h m = = = ( )
This choice allows us to somewhat moderate the decrease of the Reynolds number associated with the fast damping of the turbulent amplitude. From Equations (20) and (21), one derives the contribution of dissipative terms to the turbulent energy evolution,

u B Q F 2 2 , 2 6
d dt 2 2 diss r r + =- + n ( ) | ¯˜• ( )
where F is a flux that does not change the average turbulent energy. The turbulent heating is given by the visco-resistive damping term Q ν , which is always positive,

Q u J 4 3 , 27 2 2 2 m w h = +  + n ( ˜( ˜• ) ) ˜( ) where u w =  ´and J B = 
´. Finally, what is lost by turbulent energy is transmitted to internal energy; this reads

u u P P P P a T Q . . 2 1 , 2 8 t  g g r k g r ¶ +  +  + = D + - n ( ) ( ) ¯˜( )¯( )
where κ is the thermal conductivity. Since we are in the following more directly interested in the temperature T, we also write down its equation: We here rederive Equation (1), which expresses the critical heating leading to a temperature decrease as T R 1 p µ

u u T T T T a T Q . 1 . 2 1 1 . 29 t  g g k r g r ¶ +  + -  + - = D + - n ( ) ( ) ( ) ( ) ( ˜) ˜( ) ˜( ) 2
. For this, we need to take the spatial average of Equation (29), so eliminating the thermal conductive term. The other terms are checked to be negligible, e.g., using the dimensional factor U R

0 instead of a t  ( ), u u T T T U R . . 2 , 3 0 0 d g  =   ¯( )
which gives for the average temperature (the spatial average being denoted either by a bar or by angular brackets)

T T a Q Q 2 1 1 1 , 31 t  g g r g ¶ + - = -á ñ - n n  ( ) ( ) ( ) ˜( ) ( )
where the last equality is obtained assuming 1 dr r  . Now we replace the temporal derivative by the radial derivative, U l R 0 ¶  ¶ , and the damping term

a  by its dimensional expression U R t 0 ( ), U T T U R Q 2 1 1 . 32 R 0 0 g g ¶ + - - n  ( ) ( ) ( ) ( )
We now define Q α as the heating necessary to obtain a temperature profile of the form T R 1 = a . Thus, we replace in Equation (32) U dT dR U T R;

0 0 a = - ¯¯this gives Q T U R . 3 3 4 3 2 0 = a a - ( )
In the following, we will insist on the special value 1 a = , which leads to what we will call the critical heating Q c :

Q Q TU R 1 2
. 34

c 1 0 = = ( ) ( )
Using simple phenomenology, we now derive a condensed formula for such a critical heating involving basic parameters of the turbulent wind: the Mach number and the expansion parameter ò. We first define the Kolmogorov cascade rate, generalized to MHD, in two ways. A first definition (Vasquez et al. 2007) is given by

Q k u B ku 3 , 35
K41 2 2 32 3 d r +   ( ) ( )
with B B B d = -á ñ being the magnetic field fluctuation and both u and B d being evaluated at the Taylor's scale, thus in the inertial range. In Equation (35), u 2 and B 2 d r are, respectively, the kinetic and magnetic energy content in the wavenumber range

k k 2 , 2 [ ] ,
where k lies in the inertial range. Note that to derive the last approximate equality, we have assumed equipartition between kinetic and magnetic energy. When analyzing our simulation results, we shall use the equivalent but more precise definition,

Q E k k , 3 6
K41 3 2 5 2 = l l ( ) ( )
where E k l ( ) is twice the 1D total (kinetic+magnetic) energy spectrum, integrated on directions perpendicular to the radial and depending on the radial wavenumber, evaluated at the Taylor's wavenumber k λ . The Taylor wavenumber is defined as

k k E k dk E k dk . 3 7 x x x x x 0 2 1 2 0 ò ò = l ¥ ¥ ( ) ( ) ( ) ( )
More precisely, k l marks the middle of the inertial zone of the reduced spectrum in the radial direction. Now, we know from Vasquez et al. (2007) that in cold winds, the Kolmogorov cascade rate Q K 41 overestimates by a factor of 10 the true average dissipation rate Q n :

Q Q 10 .
38

K41 n  ( )
From the latter equation, Equations [START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF] and (35), we rewrite the critical condition

Q Q = n a to produce a R 1 a temperature profile as ku T U R 3 10 . 39 3 4 3 2 0 a -  ( )
We want to express this condition in terms of the rms turbulent Mach number,

M u c , 4 0 s rms = ( )
where c P s 1 2 g r = (
) is the sound speed and the expansion parameter ò (Equation ( 13)), which also includes the rms velocity in its definition. Assuming that the inertial range begins at the largest scale, we evaluate the kinetic energy content at the largest scale u 2 and, for the particular value

1 a = , M 4.4. 42 2   ( )
Equation (42) will guide us in selecting "critical" couples of initial Mach number and expansion parameter in our simulations.

Turbulent Energy, Expansion, and Viscous Damping

We define the average turbulent energy per unit mass:

e u B 1 2 . 4 3 2 2 r d r = á + -ñ ˜¯( )
From Equations ( 16) and (17), one finds the expression of the energy damping rate due to expansion,

Q u u B a 1 2 2 , 44 y z x exp 2 2 2  r d r = á + + ñ ˜( ) ¯( ) ( )
or, in dimensional terms,

Q u u B U R 1 2 2 . 45 y z x exp 2 2 2 0 r d r = á + + ñ ˜( ) ¯( ) ( )
The two main terms leading to the damping of turbulent energy are thus (i) the small-scale diffusive (viscous and resistive) damping term Q ν , which is fed by the nonlinear turbulent energy cascade, and (ii) the expansion damping term Q exp . However, here turbulent energy conservation is subject to further effects: (i) compressibility breaks the turbulent energy invariance, the remaining invariant being the sum of turbulent energy and internal energy; and (ii) expansion modifies the nonlinear terms as well, thus breaking down all inviscid invariants of homogeneous MHD. These two additional effects are gathered in a single term denoted by Q NL . It is found by subtracting turbulent dissipation and linear expansion decay from the time derivative of turbulent energy:

de dt Q Q Q . 4 6 exp NL = -- - n ( )
In the previous equation and from now on, Q ν (without average) denotes the spatial average Q á ñ n . Note that the very existence of a Kolmogorov-like spectral scaling k 5 3 implies that the term Q NL should be subdominant compared to Q ν .

Numerics, Initial Conditions, and Parameters

All simulations are computed in a numerical box with resolution N N N 512

x y z = = = . We start at 0.2 au with a numerical box elongated with an aspect ratio a x =5 along the radial of dimensions L L L 5 5 5 2

x y z p = ´= ´= ´. The domain is then stretched by expansion in directions perpendicular to the radial so that at 1 au the domain becomes a cube, and the mean field rotates accordingly (Figure 1).

We set up energy equipartition: u b 1 rms rms = = and about zero correlation between magnetic and velocity fluctuations (zero cross-helicity). Energy isocontours are spatially anisotropic, having the same aspect ratio as the numerical box. As the mean magnetic field in all runs makes a small angle θ with the radial (between 11°and 20°), the initial spectrum has isocontours not far from perpendicular to the radial direction with an aspect ratio equal to 5 (see Figure 9(a)). This corresponds to the so-called "2D" configuration characteristic of slow winds (Dasso et al. 2005

; Verdini & Grappin 2016).
The parameters of the simulations are listed in Table 1, which lists all runs described here. The main parameters are the Mach number and expansion parameter ò. In Table 1, ò goes from 0.12 to 0.4, M from 0.3 to 1, and M 2  from 0.45 to 5. Varying ò by a factor of 2 as done here corresponds to varying the initial turnover time of the largest scale by such a factor, so the wavenumber by a factor of 2 2 3 , adopting a k 5 3 scaling for the initial spectrum. Also, the simulation duration is longer when ò is smaller, as the transport distance (from 0.2 to 1 au) is fixed.

Each simulation provides a test window of wavenumbers (reduced to a little more than two decades) on the much larger solar wind inertial range. When choosing the couple M,  ( ) for a given run, we choose a wind regime, and at the same time we place our simulation range on the solar wind wavenumber range. Typical values, as taken in the runs A, B, C, and E, are M=1 and 0.2  = . Such values are not far from the values found in cold winds in Helios data at the scale of several hours (Grappin et al. 1991).

As seen in Section 2.3, the expression M 2  was found to be approximately 4.4 in the inertial range at 1 au by Vasquez et al. (2007) for cold/slow winds. This is in the middle of the range of values listed in Table 1.

Other parameters appearing in Table 1 are the plasma β,

P P c v 2 B s a th 2 b g = =( )(
) , which is varied from 0.29 to 1.48, and the mean field B 0 , which varies between 0.86 and 2.04, meaning that, since b 1 rms = , b B rms 0 varies from 1.2 to 0.5. The initial viscosity 0 m is about the same in all calculations except for run Z: the origin of this exception will be discussed in Section 4.2. The parameters k max and m are, respectively, the extent of the initial spectrum and its 1D spectral slope.

Results

Run A: Extended Initial Spectrum

We start with the case of run A, with M=1, 0.2  = , and thus M 5 2  = . The initial spectral extent in the direction perpendicular to the radial is k 64 max = , and the spectral slope is m 5 3 = . Figure 2 shows the evolution of several quantities versus heliocentric distance. Panel (a) shows the evolution of the rms turbulent quantities u rms , B rms 1 2 d r ¯, and u c rms , where u c is the compressible part of the velocity field. The quantities shown are multiplied by R R 0 0.6 ( ) . One sees that the compensated profiles of velocity and magnetic field (in velocity units) are close to a plateau, indicating a decay close to R 1 0.6 , somewhat faster than the Wentzel-Kramers-Brillouin (WKB) prediction

u b R 1 rms rms 1 2 1 2 r   ¯.
The compressive rms velocity amplitude, initially zero, rapidly reaches about half that of the total rms velocity and then becomes closer to 1/3 of it. Panel (b) shows the average turbulent dissipation of total (solid line), solenoidal kinetic ( 2 mw ˜; dotted line), magnetic ( J 2 h ˜;

dashed line), and compressible ( u 4 3 2 m  ( • ) ; dot-dashed line) energy. While the kinetic and magnetic dissipation are comparable as expected, the compressible dissipation decreases rapidly and becomes 1/10 of the total dissipation at the end.

Panel (c) of Figure 2 shows the total energy (Equation ( 43)) decay rate de dt | | (solid line) and its components: the expansion decay rate Q exp (dotted line), the turbulent decay rate Q ν (dashed line), and the residual term Q NL (dot-dashed line).

One can distinguish two phases: (a) a short initial transient, during which the turbulent dissipation dominates the expansion decay, in agreement with the small value of the expansion parameter 0.2  = and a very large residual term comparable to the turbulent dissipation,

Q Q Q ; exp NL
< n  and (b) the rest of the evolution, during which the turbulent decay is smaller than the expansion decay, and the residual decay is the smallest,

Q Q Q NL exp < < n .
These two points will be clarified in the discussion. We also remark that the sign of Q NL varies: it is an energy loss (thus increasing de dt | |), denoted by a thick line, but during the beginning phase ( R R 1 2 0   ; thin line), it is an energy gain, thus decreasing de dt | |.

Finally, panel (d) gives the resulting temperature curve, compensated by a 1/R law. One sees that a power-law regime appears for R 0.5 au  , with an index between 4/3 and 1. The turbulent energy reservoir is clearly used in two phases: (i) an early phase with rapid and strong dissipation that almost stops the plasma cooling and (ii) a long-lasting phase with reduced dissipation that only mildly delays the cooling of the plasma.

Varying Spectral Extent

Is the set of initial conditions made for run A the most efficient in terms of resulting temperature curve, or can we achieve a resulting curve closer to the observed 1/R decrease? Note. In the table, R is the name of the run; M u c s rms = , with c s the sound speed; ò is the initial expansion parameter; M 2  , see Section 2.3; B 0 is the initial magnetic field amplitude (and Alfvén speed); β is the ratio of thermal over magnetic pressure; k max is the maximum wavenumber in directions perpendicular to radial (nota bene largest perpendicular scale corresponds to unit wavenumber); m is the 1D spectral slope; and 0 m is the initial value of the diffusive parameters (viscosity, resistivity, and conduction). As we will see later in the discussion, an important characteristic of run A is the exaggerated importance of small scales, compared to the one found in the quasi-stationary turbulent state. This is at the origin of the excessive heating occurring in the early phase of the run. To reduce the energy of small scales, we now change one or two of the following parameters: (i) the initial spectral slope m and (ii) the initial power-law extent, as measured by k max .

In Figure 3 Run A shows an initial large overheating phase with

Q Q 10 c n  , followed by insufficient heating Q Q 1 c < n for R 0.3 au > (panel (a)
, solid line). This explains the different phases of temperature evolution considered earlier: a large part of the turbulent energy is lost during the first phase; so, in the second phase, the remaining turbulent energy is too small to substantially heat the plasma, leading to a temperature decrease between adiabatic and 1/R (panel (b), solid line).

Due to the reduced importance of small scales, runs B, C, and E show a different behavior. Initially, for R 1.2  , all of them show comparable critical heating ratios (panel (a)), too small to lead to observable heating. This phase corresponds to a quasi-adiabatic decrease of temperature (panel (b)). This is followed by a quasi-stationary regime (R 1.2  ) in which the heating is close to critical and leads to a common temperature decrease, with all three temperature curves showing a quasiplateau, thus close to a 1/R decrease (panel (b)).

Decreasing the importance of small scales in the initial spectrum thus succeeds in suppressing the too-large energy loss of the first phase. Note that run E shows the profile closest to 1/R during the whole nonadiabatic phase. In the following runs, we thus fix the spectral parameters as in run E: spectral slope 2.2 and a short spectral extent with k 4 max = .

Mach Number, ò, and M 2 

Decreasing the initial Mach number intuitively decreases the turbulent energy reservoir compared to the internal energy, so it should also decrease the heating ratio Q Q c n . In order to check this conjecture, we compare two runs, C and D, both with

k 4 max = , 0.2  =
, and, respectively, M=1 and M=0.77. Figure 4(a) shows that our conjecture for the heating ratio is correct. As a consequence, the average radial slope of the temperature profile changes substantially (Figure 4(b)). Note, however, that the temperature curve shows a break and decreases at a faster rate in the end.

The expansion parameter ò measures the expansion rate normalized by the nonlinear shearing rate. Intuitively again, a low expansion parameter should favor heating to the detriment of cooling. In order to check this second conjecture, we compare two values of the expansion parameter, 0.12  = (run F) and 0.2  = (run K), with M=0.77 in both cases. The result is as expected (Figure 5); i.e., the heating is larger for the run with lower ò during the first part of the transport for R 0.6 au < (R R 3 0 < ). However, the reverse is true for the second half of the travel. This happens because, for smaller values of ò (run F, solid line), the same travel distance corresponds to a larger number of nonlinear times (i.e., larger "age"; see Grappin et al. 1991), which may easily result in a too-fast decrease of the energy reservoir and thus of the heating rate.

Last, we test the parameter M 2  as a possible control parameter for the heating and temperature profile (see Section 2.3). We choose two pairs of runs: E and F have the largest parameter value, M 5 2  = , and G and H have the lowest one, M 2.5 2   (see Table 1). Figure 6(a) shows that the ordering of critical heating by M 2  is approximately verified for R 0.5 au < but not at larger distances. Nevertheless, the temperature curves (Figure 6(b)) ; dotted line). Same caption as Figure 3.

appear to be gathered in two groups according to the parameter value.

Varying the Plasma b and Angle VB q

Due to variations of the initial temperature, with the mean initial magnetic field remaining constant, the β of the different runs considered up to now (runs A to K) has been varied in the interval 0.3 0.8   b (Table 1). Since, in the slow solar wind, the β of the plasma can be larger than unity, we now consider runs with larger values,

(runs M and N), and compare with run K with 0.49 b =

. The Mach number is 0.77 in the three runs. The overall effect of the β variation appears to be small (Figure 7), with, however, a slight advantage (stronger heating) to the two runs with larger β (runs M and N; solid and dotted lines with quasi-superposed curves).

A last parameter is of interest: it is well known that, although the angle VB q between the magnetic field and the radial direction is, on average, 45°at 1 au (which corresponds to 11 q =  at 0.2 au), its distribution actually varies widely around the average. To test the effect of small-angle variation, we considered doubling the initial angle VB 0 q . We found that when passing from 11 VB 0 q =  to 20 VB 0 q = , the critical heating ratio, as well as the temperature profile, shows no variation at all (not shown).

Discussion

Summary

Our numerical results support the possibility that MHD turbulence can drive a proton temperature profile that is decreasing significantly more slowly than the adiabatic prediction in the distance range R 0.2 1 au < <

. We started with a spectrum initially having a 2D configuration, corresponding in principle to the slow wind regime as observed by Dasso et al. (2005). This led, with an rms Mach number close to unity and expansion parameter 0.2  = , to a temperature profile significantly steeper than observed; however, when considering a strong reduction of the initial spectral inertial range, we obtained a temperature profile close in average to a 1/R law, thus not far from the average R 1 0.9 profile measured by Totten et al. (1995).

We found that the parameters regulating the heating rate are the rms Mach number and the expansion parameter ò, combined as M 2  , while other parameters, such as plasma β and angle VB q between the mean field and the radial (for small initial values) have a minor effect.

Spectral Properties versus Mach Number

To understand the necessity of reducing the initial smallscale energy content, we examine here the spectral evolution.

We first examine the spectral anisotropy. Figure 8 shows the case of run E; it is representative of the other runs that all show comparable evolutions. The 2D spectra shown are cuts though the plane k z =0 of the 3D spectrum for total energy u B

2 2 d r +
. As explained in Section 2, the initial energy isocontours (panel (a); R=0.2 au) are quasi-perpendicular to the mean field direction (straight line), which is close to the radial direction. At 1 au, however (panel (b)), the mean magnetic field has an angle of 4 p with the radial, and the main symmetry axis of the isocontours is now quasiperpendicular to this mean field. In other words, the cascade is not only initially perpendicular to the mean field, it remains so during transport, following the rotation of the mean field. This corresponds nicely to the so-called 2D spectrum dominant in the slow wind (Dasso et al. 2005), which was first found numerically in Verdini & Grappin (2016) to be one of the two robust attractors in the wind.

We now consider the 1D reduced total energy spectra at 0.2 and 1 au for runs A, B, C, and E, shown in Figure 9. The initial spectra have dotted lines, and the final ones have solid lines (1 au). The final spectra depend on either radial (thick solid line) or perpendicular (thin solid line) wavevectors. The wavenumber is normalized by the Kolmogorov dissipation wavenumber k Q

d 3 1 4 m = n ( )
, where the dissipation rate Q n is defined in Equation [START_REF] Galtier | Turbulence of Weak Gravitational Waves in the Early Universe[END_REF]. Each final spectrum is then obtained by averaging the spectra so normalized in the distance interval R 0.6 1   . The four final spectra are all comparable (either along the radial or transverse direction), showing a very reduced spectral extent of about half a decade with a slope m 5 3  . The origin of the small extent of the inertial range in the final spectra actually lies in the high Mach number (M=1) adopted in these runs. A high Mach number leads to large intermittent variations of density associated with shocks, thus requiring large viscosities to prevent the occurrence of unresolved gradients in high-density regions. As a matter of comparison, when dealing with run Z (M=0.33), we could use a viscosity 10 times smaller than that used for runs with M=1 (see Table 1). Due to its much lower viscosity, run Z follows a k 5 3 -scaling on more than one decade, thus significantly larger than that for run A (Figure 10).

The short spectral extent seen previously for runs with M=1 thus results from the necessity to increase the viscosity with such large Mach numbers to prevent a catastrophic (unphysical) evolution of the run at a given numerical resolution. This also explains why it is necessary to start with a small spectral extent: otherwise, during a transient phase, one obtains excessive heating produced by the artificial initial excess of energy at visco-resistive scales (see the evolution of Q Q c n for run A in Figure 3(a)).

Dissipation Rate and Kolmogorov Rate

The parameter M 2  (Equation ( 42)), derived in Section 2.3, has been used to specify conditions allowing us to approach critical heating. In particular, in runs B, C, and E, the value M 5 2  = (not too far from the nominal value of 4.4 derived in Equation ( 42)) allowed us to obtain temperature profiles with a power law close to 1/R. The argument used in Section 2.3 leading to this prescription relies on the assumption that the ratio R Q Q

V K 41
= n is 0.1  . This very low value of R V (corresponding to an effective very high Kolmogorov constant) has been found to hold in cold winds by Vasquez et al. (2007). It would be satisfying to check whether or not the runs studied in this paper show the same low value of R V .

Figure 11(a) gives R V versus distance in runs A, B, C, and E. The curve varies wildly for run A, with R V passing from larger than 1 to lower than 0.1, while runs B, C, and E all show a ratio clustering around the value 0.1 on the whole distance range. Panel (b) shows the evolution of the Taylor wavenumber (Equation ( 37)) for the four runs, summarizing the spectral width evolution for the different runs. By comparing panels (a) and (b), one sees again that the initial oversized spectral width of run A leads to an anomalously large value of R V . In contrast, runs B, C, and E, which have a spectral width adapted to their viscosity, show R V values clustering around 0.1.

The fact that in our simulations B, C, and E leading to a critical heating, we find a value for Q Q K 41 n close to that observed indicates that our numerical setup, in spite of the previous remarks on the limited spectral range, leads to turbulent properties close to those of the actual solar wind turbulence.

This small value of R V means that the characteristic turbulent decay time is about 10 times longer than the plain nonlinear time t k u 1 NL =

. This allows us to interpret the relative importance in our runs of the expansion decay rate and turbulent decay rate. Indeed, with an expansion parameter 0.2  = , as in runs A, B, C, and E, one expects a priori an expansion decay rate Q exp smaller than the turbulent dissipation rate, since ò is the expansion rate normalized by the inverse of the turnover time at large scales (Equation ( 13)). However, due to R V being 0.1, the effective turbulent decay timescale is 10 times longer than the nonlinear turnover time. The effective expansion decay rate is thus finally not smaller, but rather larger than the turbulent dissipation rate, in spite of 0.2  = . This is true in the whole distance range [0.2, 1] au for runs B, C, and E, and also for run A, except during the early transient where R V is close to unity.

A corollary to the overall dominance of expansion damping in runs A, B, C, and E should be that the turbulent fluctuation amplitude decays close to the WKB prediction, which is, for Alfvén waves,

Z u B R 1 . 47 2 2 1 2 1 2 d r = +  ( ) ( )
This is indeed the case: Figure 11(d) shows that the turbulent amplitude in runs A, B, C, and E decays as Z R 1 0.6  , thus close to the WKB prediction.

Loss of Energy Conservation during Cascade

The deviation from turbulent energy conservation during cascade has been measured by the residual term Q NL , defined as the difference between the total turbulent energy decay and the sum of turbulent dissipation Q ν and linear expansion decay Q exp (Equation [START_REF] Burgers | A Mathematical Model Illustrating the Theory of Turbulence[END_REF], run A, Figure 2(c)). Other runs show that with a fixed Mach number, the residual term Q NL is proportional to the expansion parameter ò. In order to eliminate the contribution of compressibility and so to determine without ambiguity the contribution of expansion alone, we consider the deviation of total energy conservation instead of just turbulent energy. We denote the new residual term by Q NL¢ :

d dt e T Q T a Q 1 2 . 48 exp NL  r g g r + á ñ - = - -á ñ - ¢ ( ˜ (
)) ˜( )

Figure 12 shows the total (turbulent + internal) energy time derivative and the associated residual term Q NL¢ for run A. The residual dissipation is limited to 1% or 2% of the total variation. This is substantially smaller than Q NL (see Figure 11(c)). It shows that compressible exchanges between turbulent and internal energy are the dominant contribution to the deviation of turbulent energy conservation during the cascade, especially during the beginning of the evolution. The same remarks can be made for the other runs in Table 1. However, in run H, with a larger expansion parameter ( 0.4  =

), we find that compression and expansion contribute more equally (not shown).

Conclusion

In conclusion, using complete nonlinear couplings of MHD equations, we have shown that radial temperature profiles as 1/R simply result from the combination of adiabatic decrease and turbulent dissipation. This has been done starting at 0.2 au with an rms Mach number 1 and an expansion parameter 0.2  =

. With these parameters, the decrease of rms turbulent amplitude is not much faster than the Alfvén WKB prediction, actually as R 1 0.6 . This demonstration has been done by starting with a spectral anisotropy characteristic of slow winds, namely mainly perpendicular to the mean magnetic field. This included showing that the Q Q K 41 n ratio is close to 0.1 in our simulations, as in the solar wind. Finally, we measured the deviation from the conservation of turbulent energy during the cascade (residual energy loss Q NL ). We found that expansion was a minor cause of deviation, the main cause being compressible exchanges between turbulent and internal energy.

Future work includes (i) providing a clearer signature (i.e., with larger Reynolds) of the 5/3 power-law index characteristic of slow winds by lowering the Mach number and the expansion parameter and (ii) considering the case of fast winds, which, despite having a different spectral anisotropy, produce a similar radial dependence of the temperature. This work was performed using HPC resources from GENCI-IDRIS (grant 2017-040219). It has been supported by Programme National Soleil-Terre (PNST/INSU/CNRS).

ORCID iDs

Victor Montagud-Camps https:/ /orcid.org/0000-0002-7848-9200 Roland Grappin https:/ /orcid.org/0000-0001-7847-3586 Andrea Verdini https:/ /orcid.org/0000-0003-4380-4837

Heating fast winds

In the previous chapter we saw that under reasonable initial conditions at 0.2 AU it is possible to generate a proton temperature profile close to R -1 in the [0.2, 1] AU interval, solely from the contribution of turbulent heating, starting with a gyrotropic initial spectrum with M = 1, = 0.2. This spectral geometry is robust: it remains gyrotropic at 1 AU. It is thus interesting to generalize this first study by examining whether one can reproduce the 1/R temperature profile also when starting with different spectral geometries, in particular the one specific of fast winds.

Since wind heating probably needs as well are characterized by a large cross helicity, we will introduce the initial cross helicity as a parameter.

Initial conditions

The parameters are listed in table 12.1. Runs all have a resolution of 512 3 grid-points, an initial expansion parameter of = 0.2, initial spectral slope m 0 = -5/3, an initial mean magnetic field close to the radial direction, B 0 = (2, 2/5, 0). The initial spectra and the initial cross helicity have been constructed as described in section 6.1.2.

The Mach number does not differ substantially in average in fast and slow winds (cf fig. 13.1), thus we adopt M = 1 initially in our fast solar wind simulations. We wish to obtain the R-slab geometry characteristic of the fast winds, so starting with an isotropic energy spectrum is appropriate (see part II); nevertheless, we will also, for completeness, consider starting with Gyro of Gyro-Alfvén initial conditions (see section 6.1.2). We consider either σ 0 c = 0 or a large initial cross helicity. The initial domain aspect ratio is R x = 1 or R x = 5. (See table 12.1) As already remarked, a Mach number M = 1 demands a high viscosity, which shifts the dissipation wave-numbers to large values. An initial spectrum with a large part of the initial turbulent energy attributed to quasi-dissipative wave-numbers would be a waste of energy, probably leading to a high heating in a first phase, thus spoiling the temperature profile (see Chapter 11). For this reason, a strong truncation k max = 8 is adopted. The only exception has been run RC7, with a stronger truncation of k max = 4. 12.2 Results

Run

Type R x k max m 0 σ c M µ 0 RC1 ISO 1

Spectral anisotropy with Mach=1

In this section we present the spectral evolution of runs with large starting turbulent Mach number (M = 1), thus complementing the results of part II which was dealing with lower initial Mach number (M = 0.3). The three figures in this section start with successively Gyrotropic, Isotropic and Gyro-Alfvén geometries. For each initial model, we vary the initial cross-helicity, comparing σ 0 c = 0 and 0.8 and also as well 0.6 for the G+A model. Figure 12.1 shows the case of the ISO model. The spectral anisotropy is similar to that observed in M = 0.3 runs in part II: gyrotropic or 2D for the zero cross helicity run, and R-slab for high cross helicity. The isocontours' aspect ratio does not differ much between the former series and the present one, with values ranging between 0.4 and 0.6. The tilt angle α is evaluated in the range k maj = 2 and k maj = 10: |α| ≈ 30 • for run RC1, whilst for RC2, |α| ≈ 15 • . This is close to the results shown in figure 8. [START_REF] Anders | Abundances of the elements: Meteoritic and solar[END_REF].

The 1D spectra in figure 12.1 are best developed in the direction z, namely perpendicularly to the mean magnetic field: this is true whatever the initial cross-helicity. The inertial range in the radial direction is very short, with an approximate index m = -2, while it is m = -5/3 in the z direction.

The initial gyrotropic anisotropy of RC3 and RC4 is robust, whatever the cross helicity, (no figure is shown).

For the G+A simulations RC5, RC6 and RC7, the differences in anisotropy between the low and high cross-helicity are largest. Figure 12.2 shows that for the zero cross helicity run, RC5, the anisotropic index A takes values around 0.7 below k maj = 5 and then decrease continuously until A = 0.4. For the high cross-helicity run, RC6, A shows a smaller variations, ranging from 0.2 to 0.3. On the other hand, RC7 stays around A = 0.6, with a sudden increase beyond 10k maj . As it is explained in chapter 7, this increase is just caused by the isocontours reaching the boundaries of the domain and opening. Thus, the increase is just an artifact of the boundary conditions.

The tilt α stays below 10 • for the high cross-helicity run. On the other hand, for zero cross helicity σ 0 c , α also stays below 10 • only for k maj > 5. Below this mark, the inclination varies between 30 • and 60 • , as the rectangular shape of the isocontours forces to take the The reduced energy spectra scale as k -3 in the radial direction indenpently of crosshelicity. In the z direction however, the reduced spectra for zero cross helicity shows a spectral index m -2 from 0.05 to 0.2 k z /k d . For high cross-helicity σ 0 c , a short inertial range with the same index spans from 0.1 to 0.2 k z /k d .

Turbulent heating

In figure 12.3, one sees the evolution with distance of some important properties of the simulations listed in table 12.1. First the rms Mach number, then the cross helicity, the heating normalized to critical heating Q/Q c , then the compensated temperature. Heating is generally subcritical in the beginning, then reaches a maximum at about R/R 0 = 1.5. After that, both ISO behave alike in terms of heating and independently of the initial cross-helicity, the ratio Q/Q c decreasing slowly. For both runs, Q/Q c 1 at R/R 0 = 3.

Mach numbers decrease for both isotropic runs (RC1 and RC2) and for the gyro both gyro-Alfvén runs with non-zero initial cross helicity (RC6 and RC7). The decrease is faster for both isotropic runs and for the gyro-Alfvén run with the highest cross helicity (RC6), taking values around 0.6 at the end of the simulation. The other runs, maintain the Mach number around 1 during the hall simulation. The more weight has the gyrotropic component of initial anisotropy, the higher the Mach stays.

Low initial cross helicity stays below 0.1 and high initial cross helicity slowly decreases, taking values between 0.3 and 0.7 at the end of the simulation. We have not found any way to explain this evolution of cross helicity in terms of the other considered parameters.

If we look now at how the variation in Q/Q c corresponds to the evolution of proton temperature we see that between R/R 0 = 1 and R/R 0 = 2 there is a fast increase of temperature, which was expected considering that the ratio Q/Q c obtains its maximum values in this interval. From R/R 0 = 2 to R/R 0 = 4, the proton temperature profile is closer to R -1 , although not as constantly as in the paper (Chapter 11) for simulations with a similar initial set-up (e.g. run E, fig. 3), with the exception of runs RC3 and RC7.

For both runs with isotropic initial spectrum, in the interval (2 ≤ R/R 0 ≤ 4), Q/Q c diminishes from 1.5 to 0.8 approximately. After that, Q/Q c continues to diminish until 0.7, and the temperature cools faster than 1/R.

The Gyro run with low cross helicity reproduces the results of run C in chapter 11, obtaining a ratio Q/Q c around 1 in the interval 2 ≤ R/R 0 ≤ 5 and an evolution of proton temperature T p ∝ R -1 in the same interval. We remark that in this case, the initial truncation has been relaxed from k max = 4 to k max = 8, and the initial aspect ratio has also changed from 5 to 1. Both changes imply giving more energy to modes higher wave-numbers than for run C in chapter 11. This implies that the variations in k max do not affect the turbulent heating, as long as the dissipative scales are not reached.

The Gyro simulation with initial cross helicity 0.8 shows remarkable differences with respect its zero cross helicity version. Although from R/R 0 = 1.2 to R/R 0 = 2 turbulent heating attains the critical heating level, it rapidly diminishes and stays around Q/Q c = 0.8 for R/R 0 > 2. Correspondingly, proton temperature follows a power law intermediate between the adiabatic prediction and the observed R -1 for that interval.

For the run with gyro-Alfvén initial anisotropy and σ 0 c = 0, the peak of heating at R/R 0 = 1.5 is about 1.3. After the peak, the ratio Q/Q c stays close to 1 for the rest of the simulation. Conversely, after a small increase in temperature at R/R 0 = 1.5, T p evolves almost perfectly as a power law, with an index slightly higher than -1.

On the other hand, for RC6, the G+A run with σ 0 c ≈ 0.8, the ratio Q/Q c stays between 2 and 1.5 for the interval R/R 0 ∈ (1.5, 3). After that, the ratio stays between 1 and 0.8 for R/R 0 > 3. Proton temperature makes visible these two intervals, showing a power law evolution T p ∝ R α with α > -1 for the first interval and close to the inverse of heliospheric distance fo the second interval.

The last Gyro-Alfvén run, with σ 0 c is characterized by stable value of the ratio Q/Q c close to 1 and a proton temperature evolution following the power law R -1 in the interval R/R 0 ∈ (1.5, 4) and a small departure from that behavior for R/R 0 > 4. We associate the more stable behavior and the absence of overheating with respect the other gyro-Alfvén simulations, to the stronger truncation of energy in run RC7, as it is the main variation with respect the other Gyro-Alfvén simulations.

Chapter 13

Discussion

In this third part of the thesis, we have focused on the study of turbulent heating between 0.2 and 1AU and its relation to turbulent anisotropy.

We have started with a chapter dedicated to the 1D hydrodynamic turbulence with expansion. The results obtained in this chapter cannot be directly linked to solar wind observation, since solar wind plasma turbulence cannot be described with such a model (e.g. the spectral index -2 of the energy spectra is nowhere to be found in the solar wind). However, it has anticipated some of the results obtained in the following chapters: 1) the need of low values for the expansion parameter to allow the development of turbulence in the transverse directions; 2) the major contribution of expansion in the damping of coherent energy; and 3) the Kolmogorov estimation for the cascade rate Q k41 = k 0 u 3 /a(t) is a good approximation for the viscous dissipation (although for the EBM simulations they are just proportional).

Chapter 11 was devoted to the study of turbulent heating in slow winds with EBM simulations. In this chapter, it is shown for the first time that turbulent heating is capable to attain the critical heating level and make proton temperature evolve as the inverse of heliospheric distance. Other important contributions from this chapter have been the parameters needed to obtain a critical heating and a proton temperature evolution like R -1 : i) the ratio M 2 / ≈ 4.4 for the initial Mach and expansion parameters allows to approach the critical heating; ii) a strong initial truncation of the energy spectra is necessary to avoid an early dissipation of all the energy; iii) despite the use of relatively high Mach numbers (M = 1), turbulent heating is mainly incompressible.

The aforementioned chapter can be seen as a particular case since all runs had gyrotropic initial spectra and zero cross helicity. In chapter 12, the last of this part, we have studied the variation of turbulent heating in runs with different cross helicity and initial geometry. Our aim has been to reproduce a proton temperature evolution T p ∝ R -1 for slow and fast wind conditions (low and high cross helicity respectively) and to determine the possible influence of turbulent anisotropy. For the last, we have also verified whether the results obtained in part II stayed true for high Mach number simulations.

The results regarding turbulent anisotropy at high Mach number show similar properties to those observed for low Mach number simulation in part II, at least for the inertial range scales. Beyond the inertial range, the isocontours follow the stretching of the box in the directions transverse to the radial. Thus, for low k maj , we have verified that increasing the initial cross helicity turns the inclination of the isocontours towards the radial and if σ 0 c = 0, the cascade is mostly perpendicular to B 0 .

The images of the reduced spectra display shorter radial inertial lengths and steeper slopes compared to simulations in part II, specially in the simulations with gyrotropic and gyro-Alfvén initial conditions. Note that for the GYRO simulations and the GYRO component of the G+A simulations, k max = 8/5 = 1.6 by construction of these initial geometries. As a consequence, the former two initial geometries are quasi-2D at t = 0. This make difficult the development of turbulence in the x direction as it limits the nonlinear interactions with modes in the other two directions. We assume that changing the initial aspect ratio of the isocontours in the GYRO and G+A initial geometries (imposing for example k max = 8/2 in the x direction) could have led to similar reduced spectra to those seen in part II.

For what concerns the turbulent heating results of chapter 12, we have managed to confirm that turbulent heating, obtained in slow and fast wind conditions, is able to reach and maintain the critical heating level so proton temperature can evolve approximately as R -1 . The information provided by the simulations goes beyond this achievement.

Despite the parallel inertial ranges for gyrotropic and gyro-Alfvén runs with zero initial cross helicity were very short and steep, their turbulent heating remained at the critical level for most of the time. All simulations that have shown the same feature also had a well developed inertial range in the z direction, perpendicular to the mean magnetic field. Hence, a radial cascade is not necessary to provide enough turbulent heating to the solar wind.

Simulations with isotropic initial conditions have shown that proton temperature can evolve as 1/R even though Q/Q c = 1 is not exactly maintained during the whole simulation. Variations of the ratio between 0.8 and 1.2 also allow to obtain a fair approach to the 1/R temperature gradient, even if the value of the ratio is not constant in time.

The influence of σ 0 c in turbulent heating has proven to be linked to the turbulent anisotropy. Isotropic runs have shown ratios Q/Q c and proton temperature evolutions almost superposed, despite the difference in cross helicity. In contrast, with gyrotropic initial spectra and the other same initial parameters of the ISO runs, low cross helicity is necessary to maintain turbulent heating at the critical heating level. Note that the temporal evolution of cross-helicity for each geometry cannot justify these differences, since it maintains larger values for cross-helicity in the ISO simulation than in the gyro one. For now, we have not found any explanation for the strong influence of cross helicity in the heating of gyrotropic initial geometries.

Finally, gyro-Alfvén initial geometry has showed to be a plausible initial anisotropy for fast and slow solar wind in terms of heating. The runs with this initial geometry showed that they can achieve and maintain similar or even higher levels of turbulent heating than with the other initial geometries.

Numerical parameters vs Helios data

In fig. 13.1, we have plotted data from the Helios mission: rms Mach number (top panel) and expansion rate (bottom panel) vs day of year, during four months, at solar minimum. Rms quantities are computed in two ways, summing energies starting (i) with frequency f = 3.310 -5 Hz (thick solid lines); (ii) with frequency f = 5.510 -5 Hz (dotted lines). One has also plotted the wind bulk speed (plain solid line) and the heliocentric distance in arbitrary units (varying from 0.3 to 1 AU, dashed line). It is seen that the Mach number oscillates around unity and around 0.2, as considered in this work.

Here it is important to remark an important difference in the definition of the Mach number in this figure and in our numerical simulations. In this figure, the rms Mach number is defined each day by computing the ratio of the rms velocity fluctuation at the required frequency, normalized by the sound speed associated to the proton temperature only, with no attempt to take electrons into account:

M helios = u rms /c s = u rms /( 5 3 T p ) 1/2 (13.1)
while in the simulations it is

M = u rms /c = u rms /( 5 3 P/ρ) 1/2 (13.2)
where the pressure is P = 2nk B T , T being identified to the proton temperature and, as we know, the factor 2 coming from the fact that T e = T p in the MHD system of equations that we integrate. This remark deserved to be made, although it is probably not basic, in the absence of a realistic model of the electron dynamics in this work.

Chapter 14

Conclusions

Obtention of Maltese Cross components at 1AU

The results concerning turbulent anisotropy have shown that the increase of Mach number has a moderate effect on turbulent anisotropy. Runs with high Mach number (M=1) have shown a slight tilt of the isocontours towards the radial in comparison with the simulations with low Mach number (M=0.3). The influence of the other parameters on turbulent anisotropy is independent of the choice of the initial Mach number.

For wave-numbers within the inertial range, we have observed that the numerical simulations with zero initial cross-helicity develop a 2D cascade perpendicular to the mean magnetic field, regardless of the variation of other initial parameters. In contrast, when the inhibition of non-linear couplings by expansion and cross helicity are present, the energy isocontours become more aligned with the radial direction. In this case, the tilt of the isocontours depends on the initial geometry, as we detail now.

An initial gyrotropic distribution of energy in Fourier space (GYRO simulations) favors the perpendicular cascade from the beginning of the simulation. On the other hand, simulations with isotropic and gyro-Alfvén initial spectra (ISO and G+A ), develop turbulent cascades more aligned with the radial direction, that is, a radial-slab anisotropy ([102], VG16). Among these, the G+A simulation has the more radially aligned anisotropy, mimicking the situation of a system completely dominated by expansion (as the one shown in section 3.2.1).

The above results for ISO simulations are also observed for high values of the expansion parameter. In contrast, in the limit of no expansion, ISO simulations develop a 2D cascade perpendicular to B 0 , independently of the initial cross helicity. With the ensemble of ISO simulations we managed to reproduce for the first time the observations of wave-vector anisotropy observed at 1AU in the solar wind at MHD scales, including the components that form the Maltese Cross. On more fundamental level, we showed that expansion is physically needed to obtain the observed anisotropies.

Can the Maltese Cross heat the solar wind?

The use of the EBM equations have allowed us to verify for the first time that turbulent heating is able to produce a power law evolution of proton temperature, R -1 , between 0.3 and 1AU for fast and slow winds. Thus, the work presented here renforces the idea that MHD energy cascade rate could be the main contributor to the Solar Wind dissipation rate.

Turbulent anisotropies at the end of the simulations, that is, the "2D" and the R-slab anisotropy do not determine the proton temperature evolution. It has been possible to obtain both anisotropies and achieve the critical level of turbulent heating necessary to obtain the proton temperature evolution T p ∝ R -1 .

In contrast, the initial anisotropy that we have used (isotropic, gyrotropic and gyro-Alfén) have more influence on turbulent heating. This influence however, depends upon the initial value of cross helicity. Both isotropic and gyro-Alfvén initial conditions can arrive and maintain the critical heating level with similar initial conditions. In the gyrotropic case, with the same initial parameters, this is only achieved for zero initial cross-helicity.

In addition, we have found that the condition M 2 / ≈ 4.4 was derived in chapter 11 to obtain the critical level heating. This condition has proven to be useful for the choice of our initial Mach number and expansion parameter. Both values (M = 1 and = 0.2) are close to in-situ measurements in fast and slow winds. Despite the use of this Mach number, incompressible turbulence had the major contribution to turbulent heating.

Finally, our results have shown that the development of an inertial range in the directions perpendicular to the mean magnetic field is key to obtain an acceptable turbulent heating level. Regardless of the initial conditions, all simulations that obtain a proton temperature evolution close to observations has shown a well developed inertial range in the perpendicular direction to B 0 , even when the reduced spectra in the parallel direction were shorter or steeper than k -5/3 .

Expectations from Solar Orbiter and Parker Solar Probe

The validation and the refutation of physical models is achieved via experimental data, by improving the sensibility of previous measurements and by extending the reach of these measurements to scales and places never explored before. In the particular case of Heliospheric and Solar physics, multiple questions and theories strongly depend on the properties of the plasma at the lower solar corona (including our work). It is in this context that the missions Parker Solar Probe (launched on August of 2018) and Solar Orbiter (launch expected at 2020) were designed to obtain in-situ measurements at the closest distance from the solar surface ( [START_REF] Fox | The Solar Probe Plus Mission: Humanity's First Visit to Our Star[END_REF], [START_REF] Müller | Solar Orbiter . Exploring the Sun-Heliosphere Connection[END_REF])

Parker Solar Probe (PSP) will be the first spacecraft to approach to less than 10 solar radii at the perihelium of its orbits, which are constrained to the ecliptic plane. Solar Orbiter (SO) on the other hand, will reach distances at its perihelium already explored by previous mission, (0.28AU at its maximal approach to the Sun) but will also provide measurements at these close distances outside the ecliptic plane (latitudes up to 25 • at the end of the nominal mission), allowing to directly measure the properties of fast streams. Remote observations will be also combined to the in-situ measurements carried by both probes.

The programmed orbits of PSP and its multiple passages at 10 solar radii allow the insitu measurements of plasma properties at the heating region of the solar corona, responsible for the solar wind acceleration. Hence, we expect that the data collected by PSP can clarify the role of turbulence in the question of coronal heating via the measurement of plasma properties at this distances. We are particularly interested in the measurement of the turbulent heating rate in the accelerating region of the solar wind (not possible to describe with the EBM equations), since it could improve the choice of our initial conditions: for instance, more information about the formation of the large scale spectral range in f -1 could make us treat it as a turbulence forcing instead of a "fossile" energy reservoir. In addition, it could improve the estimation of the radial distance at which collisions between solar wind particles become negligible, a crucial parameter for models of solar wind alternatives to an MHD description, the exospheric models [START_REF] Lemaire | Kinetic Models of Solar and Polar Winds[END_REF], [START_REF] Lamy | A kinetic exospheric model of the solar wind with a nonmonotonic potential energy for the protons[END_REF], [START_REF] Pierrard | Solar Wind Electron Transport: Interplanetary Electric Field and Heat Conduction[END_REF].

As stated in the previous paragraphs, the strong points of SO mission with respect to former missions will be the improvement of measurements at already explored distances (manly explored between 0.3 and 0.5 AU) and the direct recollection of data from coronal holes, situated out of the ecliptic plane. SO mission will provide new measurements of the velocity distribution functions of protons, electrons and heavy ions, and of the turbulent levels and heating rates possibly before the cascade is fully developed. Hence, by comparison with the measurements at other distances, it will be possible to give a better estimation of the distribution of MHD turbulent heating among the different species of the solar wind, and by extension, re-evaluate the importance of turbulence as a local source of proton heating. In addition, we also expect that the data collected by SO magnetometers at the multiple passage of the spacecraft by 0.3 AU, will supply enough data to analyze the spectral anisotropy at this distance (specially for fast winds), as done previously at 1AU ( [65][23]).

Open questions and (partial) answers

We know from [START_REF] Dasso | Anisotropy in Fast and Slow Solar Wind Fluctuations[END_REF] that the fast and the slow wind have different spectral anisotropies at 1 AU. What is then the most probable pattern at, say, 0.3 AU? The data that we will be able to obtain from Solar Orbiter mission will tell us more about the turbulent anisotropy of the solar winds at 0.28 AU. But until the mission has gathered enough data, we can extract some predictions from our results.

Before our work we would have said as found in [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF]: at 0.3 AU, the most probable geometry for the slow wind is the gyrotropic pattern (also 2D), and the ISO pattern for the fast wind.

We know there are several ways to obtain the R-Slab anisotropy for fast solar winds: ISO initial spectrum and high initial cross helicity or initial excitation of large scale modes in one direction and zero initial cross helicity (VG16 [START_REF] Verdini | BEYOND THE MALTESE CROSS: GEOMETRY OF TURBULENCE BETWEEN 0.2 AND 1 au[END_REF]). Both situations cause the inhibition of nonlinear couplings in the directions transverse to the radial necessary to obtain the R-Slab anisotropy. However, we consider that the initial anisotropy proposed by VG16 is difficult to justify, as it would need an unknown physical mechanism to act before 0.2AU. Thus, our preferred initial conditions are: (i) Gyro with σ 0 c = 0 for the slow wind, and (ii) the ISO with large cross helicity for the fast wind. We could also propose as an alternative the Gyro-Alfvén model with large σ 0 c for the fast wind but, remembering the very poor scaling of these runs in the radial direction, we prefer the ISO model.

Then, what can this tell us on the turbulent regime in the fast wind? Before our work, a long-lasting time idea was that perhaps large σ c could lead to a special turbulent regime (e.g. [START_REF] Grappin | Dependence of MHD turbulence spectra on the velocity field-magnetic field correlation[END_REF])

Now we know that without expansion a large cross helicity changes nothing to the gyrotropic cascade: even a value as high as σ 0 c = 0.8 leads to gyrotropy whatever the initial spectrum. But we know also that this changes completely with some expansion: then the cross helicity has a big effect, leading to the strange R-SLAB, with a spectrum elongated in the radial direction.

Last, what can we say reasonably about turbulent heating ? Before our work, one had only tried models devoid of detailed description of the nonlinear couplings (of course, one could integrate more sophisticated models as MHD HALL or hybrid models, but then expansion was absent, or space dimension was 2 or 1). Now, we have true 3D simulations including the anisotropic wind expansion, and we succeeded in heating reasonably the plasma. However, we are forced to severely truncate the initial spectrum, otherwise all energy is lost in the first moments. What the initial conditions really are e.g. at 0.2 AU, we have an idea, because after all the turbulent cascade relaxes relatively rapidly, but we don't really understand what the relaxation is. We only know that the most important role is probably played in all winds (fast and slow) by the cascade in direction z, that is perpendicular to both the radial and the mean field direction, otherwise obtaining the same heating rate and/ or temperature profile in both winds would be strange.

Chapter 15

Future work: Anisotropy temperature description

The EBM equations assume the MHD closure, that is, that proton temperature is isotropic. This assumption is not always valid. As it can be seen from 2D cuts of the 3D velocity distribution functions of protons ( [START_REF] Marsch | Solar wind protons: Three-dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU[END_REF]), while the core of the velocity distribution function for slow winds is isotropic, the distribution for fast winds is perpendicular to the mean magnetic field axis. Consequently, temperature anisotropy is more important for fast winds than for slow winds. For instance, in situ measurements have confirmed that temperature perpendicular to B 0 , T p ⊥ , decreases slower with heliospheric distance than the parallel temperature, T p ( [START_REF] Matteini | Evolution of the solar wind proton temperature anisotropy from 0.3 to 2.5 AU[END_REF], [START_REF] Hellinger | Heating and cooling of protons in the fast solar wind between 0.3 and 1 AU: Helios revisited[END_REF][START_REF] Hellinger | Proton thermal energetics in the solar wind: Helios reloaded[END_REF] ).

By neglecting temperature anisotropy in the EBM description, we are not only unable to reproduce the different evolution of T p and T p ⊥ . In addition, we are also leaving aside all the instabilities triggered by temperature anisotrpy, such as the beam, ion cyclotron or parametric instabilities ( [START_REF] Manfred | Stability analysis of double-peaked proton distribution functions in the solar wind[END_REF], [START_REF] Tenerani | The Parametric Instability of Alfvén Waves: Effects of Temperature Anisotropy[END_REF]).

Following a previous work by Fabrizio Cametti, Roland Grappin and André Mangeney (not published), we will implement temperature anisotropy to the EBM equations. This new description of solar wind plasma will use a modified version of the classical Chew-Goldberger-Low closure, also known as CGL closure or double-adiabatic closure ( [START_REF] Chew | The Boltzmann equation an d the onefluid hydromagnetic equations in the absence of particle collisions[END_REF]).

It has already been shown that the invariants from the classical CGL equations are broken in the solar wind ([61]). In a new CGL-EBM description for the solar wind plasma, the CGL invariants will be altered and possibly closer to reality. Moreover, proton temperature in the parallel and perpendicular directions to the mean magnetic field will be described separately in two different equations, allowing to take into account the temperature anisotropy observed in the solar wind and its possible effects on turbulent heating and turbulent anisotropy. The second term on the right-hand side accounts for all the (k + N m)th modes whose value is the same as the k-mode on the points of the grid. These modes are said to "alias" the kth-mode on the grid, giving rise to the so called "aliasing errors". Usually, there are two possibilities to "de-alias", that is, to avoid these errors in the representation of variables in the discretized Fourier space. The first way is to truncate the third of the modes with largest wave-number at each time-step of the simulation. The second one consists in computing the transformation on the original grid and in another mesh that is half a grid-cell of distance shifted from the original one (see "Spectral methods in Fluid dynamics" by [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF] for an extensive description of each de-aliasing method).

We use the following strategy in our code to minimize the numerical errors. Instead of the de-aliasing methods previously described, we adopt values for the viscosity, resistivity and thermal conductivity high enough to dissipate the energy at mesh scale and thus, most of the aliasing errors. In this way, we enforce the conservation of energy which is necessary to study the turbulent heating.

At the frequencies at which we study solar wind plasma, Mach number has values close to 1. At these values, density fluctuations can create regions where the kinematic viscosity is not high enough to dissipate the energy accumulated at the mesh-size scales.

In all simulations, we set equal viscosity, resistivity, and thermal conduction at the initial time. For simulations with plasma expansion, all these coefficients diminish as the inverse of heliospheric distance, in order to maintain a constant Reynolds number (see an explanation of this in the attached paper in part III). Following a suggestion by Prof. Thierry Passot we have split the kinematic viscous term into a compressible and a solenoidal component. The kinematic dissipation term in equation (2.26) becomes,

∂ t u| diss = - µ s ρ (∆u) - µ c ρ ( 1 3 ∇ • u) (B.3)
where µ s and µ c are the solenoidal and the compressible viscosity, respectively. The intention was to increase the compressible dissipation to moderate the steepness of shock-waves and, at the same time, diminish the solenoidal viscosity so as to achieve higher Reynolds numbers for the Alfvénic terms. Unfortunately, the improvement of the Reynolds number was minor with respect to the use of a single viscosity because the required solenoidal viscosity to dissipate energy at mesh scales was at most 20% smaller than the original one. The several runs that had to be launched to find the optimal choice for the compressible and solenoidal dissipation for each new simulation were not worth such a small decrease in the Reynolds number. As a consequence, we have chosen to continue using only one viscosity.

B.2 Time integration method

Along with the pseudo-spectral methods that evaluate the spatial gradients in the EBM equations, it is also necessary to use a numerical method to integrate in time the same equations. The time integration method we use is the third order Runge-Kutta method. The time step is chosen accordingly to the Courant-Friedrich-Lax (CFL) condition, which assures the stability of the time integration.

Let ∆x be the distance between two points of the mesh, dt the time-step in the simulation and dt ph the characteristic time at which perturbations propagate in the system. According to the CFL condition, for the third order Runge-Kutta method to be stable, it is necessary that dt < dt ph (B.4)

For instance, in the simulation of the Navier-Stokes equations, perturbations propagate at sound speed c s = ∆x/dt s , where dt s is the the characteristic time for sound waves.

In the EBM simulations, several characteristic times coexist in the simulation: the sound time, the Alfvén time, the dissipation time and the expansion time, defined as, Appendix C

Résumé en français

Le but de cette thèse est l'étude du développement de la turbulence dans le vent solaire entre 0.2 et 1 unité astronomique du soleil (l'orbite terrestre). L'étude est faite en résolvant numériquement les équations de la MHD après soustraction de l'écoulement moyen radial (Expanding Box Model ou EBM). Les deux aspects de la turbulence qui nous intéressent sont la structure 3D des spectres d'énergie et le chauffage du plasma qui résultent de la dissipation turbulente des tourbillons et couches de courant emportés par le vent. Il s'agit de problèmes inverses: on cherche à déterminer quelles sont les conditions initiales du plasma à 0.2 UA qui permettent d'aboutir à ce qu'on observe à 1 unité astronomique. Un but important de mon travail est aussi de déterminer si la physique qui est présente dans les équations que j'intègre (la Magnétohydrodynamique) suffit pour arriver à reproduire ce qu'on observe dans l'intervalle où se font les mesures in situ, entre 0.3 et 1 UA.

Notre thèse contient quatre parties. Dans la première, nous donnons les équations de base, une introduction à la turbulence et au vent solaire, ainsi qu'un exposé de nos motivations dans cette thèse. La deuxième partie est consacrée à l'étude de l'anisotropie de la cascade turbulente, et plus précisément à l'anisotropie du spectre 3d, considérée comme un diagnostic de la nature de la turbulence, et en même temps un moyen de remonter aux conditions plus près du soleil. La troisième partie est consacrée au chauffage turbulent, qui apparaît comme une conséquence universelle de la turbulence, quel que soit le régime.

Partie I: Introduction

Le vent solaire est un plasma accéléré dans la direction radiale à partir d'environ dix rayons solaires de la surface du Soleil. Dans cette première partie, nous présentons les propriétés suivantes: température des différentes espèces, densité, fonction de distribution des vitesses . . . Nous continuons avec la description du plasma du vent solaire dans une approche fluide, en particulier celle de la boite en expansion ou EBM.

La turbulence et ses propriétés sont définies de façon générale successivement dans le cadre des fluides neutres et celui des plasmas. Dés que ces concepts de base sont établis, les deux problèmes principaux traités dans cette thèse sont présentés: l'anisotropie de la turbulence à grandes échelles et le chauffage turbulent du vent solaire. Dans les deux cas, nous rappelons les travaux observationels qui ont mené à ces questions et les différentes tentatives d'y répondre avant nous. Cette partie finit par la présentation du plan que nous avons suivi pour faire face aux deux problèmes proposés.

Partie 2: Anisotropie de la turbulence dans le vent solaire Dans la zone inertielle, les mesures in situ vers 1 UA montrent des figures turbulentes complexes représentées ici dans l'espace de Fourier par des spectres d'énergie. Dans le vents lents, on observe une géométrie spectrale qui correspond à une cascade d'énergie dans le plan perpendiculaire au champ magnétique moyen. Pour les vents rapides, à fréquence pas trop élevée, les isocontours d'énergie se présentent alignés le long de l'axe magnétique. Cette deuxième géométrie est connue dans la littérature sous le nom de "Slab turbulence" et la première géométrie, caractéristique des vents lents, sous le nom de cascade perpendiculaire ou simplement "2D".

Ces observations peuvent être interprétées de plusieurs façons, parce que la méthode utilisée pour obtenir les autocorrélations est basée sur une hypothèse de gyrotropie (symétrie par rapport au champ magnétique moyen), ne permet pas de reconstituer une géométrie qui ne serait pas à symétrie axiale par rapport au champ moyen. Mais les simulations numériques nous ont permis de lever l'ambiguïté et de reconstituer complètement la structure 3D.

La question est en effet de savoir quand intervient l'axe soleil-terre, et quand intervient l'axe du champ magnétique moyen dans la définition de la géométrie 3D, c'est-à-dire la forme des structures turbulentes.

Dans ce cadre, nous avons montré les résultats d'une série de simulations où nous avons étudié l'effet de trois paramètre reliés au développement de la turbulence: la pente spectrale initiale du spectre 1D (réduit à sa dépendance par rapport à l'échelle dans la direction soleil-terre) , le paramètre d'expansion et l'hélicité croisée.

Pente initiale

La pente initiale du spectre 1D (réduit à sa dépendance par rapport au nombre d'onde radial) caractérise la distribution initiale de l'énergie des fluctuations magnétiques et cinétiques entre les grandes et le petites échelles le long de la direction soleil-terre. Avec quelques modes seulement excités au début de la simulations, une turbulence ordinaire (dans l'atmosphère terrestre par exemple) devrait être aveugle à un changement de la pente spectrale initiale. Pourtant, nous avons trouvé que dans certains cas où les interactions non-linéaires sont affaiblies, la distribution initiale de l'énergie entre les différentes échelles compte et affecte fortement le développement ultérieur de la turbulence.

Paramètre d'expansion

L'expansion transverse du vent solaire affaiblit les couplages non-linéaires dans le système. Cette expansion est plus intense à plus basses fréquences. Le réglage du paramètre d'expansion dans nos simulations nous a permis d'étudier l'inhibition de la cascade turbulent par l'expansion à plusieurs fréquences. En prenant aussi des valeurs initiales pour l'hélicité croisé proches des valeurs mesurées dans les vents lents et les vents rapides, nous avons réussi à reproduire pour la première fois l'anisotropie turbulente observée le long des fréquences de la zone inertielle.

Hélicité croisée

L'hélicité croisée mesure l'efficacité des couplages non-linéaires présents dans les équations de la MHD dans sa limite incompressible, simplification utile pour décrire les régimes turbulents. Dans cette limite, les couplages se font entre ondes d'Alfvén se propageant en directions opposées. Les couplages non-linéaires (et donc la turbulence) meurent si un seul des deux modes d'ondes est présent. Dans le vent solaire, on trouve dans les vents rapides une situation de ce genre, avec une amplitude très réduite d'un des deux modes. Pourtant les spectres montrent une loi de puissance caractéristique des régimes turbulents, ce qui est contradictoire: cette contradiction a été relevée très tôt par [START_REF] Dobrowolny | Fully Developed Anisotropic Hydromagnetic Turbulence in Interplanetary Space[END_REF].

Dans notre thèse, nous retrouvons la relation observée entre la géométrie des spectres et l'hélicité croisée qui apparaît dans les vents rapides. Par ailleurs, on sait qu'à petite échelle, la géométrie redevient 2D même dans les vents rapides: nous expliquons cela par la diminution du taux d'expansion qui rend l'hélicité croisée moins efficace pour tuer la cascade.

Partie 3: Chauffage turbulent

La troisième partie est centrée sur le chauffage turbulent dans les vents rapides et lents. Entre 0.3 et 1 UA, la température des protons diminue anormalement lentement, ce qui indique une source de chauffage, qu'on suppose ici être la dissipation des tourbillons et couches de courant emportés par le vent depuis la couronne solaire. Nous commençons par exposer nos résultats utilisant un modèle 1D de turbulence (équation de Burgers) qui reproduit le raidissement des ondes sonores et la formation de chocs. Nous montrons comment le taux d'expansion du volume de plasma emporté par le vent permet de tout prédire analytiquement: la décroissance de la turbulence et le taux de chauffage.

Ensuite, nous étudions l'évolution d'un volume tri-dimensionnel de plasma emporté par le vent. Nous examinons en premier le cas de conditions initiales correspondant aux vents lents. Le cas des vents rapides est abordé ensuite. Dans les deux cas, on adopte des anisotropies spectrales différentes, tenant compte de nos résultats de la partie 2 sur l'anisotropie.

Dans les deux cas, et avec des paramètres (nombre de Mach rms proche de 1, taux d'expansion égal à 0.2) représentatifs des valeurs observées, nous observons numériquement que la turbulence transfère l'énergie des fluctuations à l'énergie interne du système à un taux satisfaisant, c'est-à-dire suffisant pour créer une évolution de la température des protons en loi de puissance de la distance R au soleil du type T p 1/R.

Un tel gradient de température est proche de celui observé par plusieurs auteurs. Ce résultat nous permet d'affirmer que la turbulence est une source possible de chauffage local du vent solaire qui contribue majoritairement au ralentissement de la diminution de la température des protons avec la distance héliosphérique entre 0.3 et 1 UA. Ceci nous permet de répondre à la question initiale: oui, des régimes turbulents différents (et des géométries turbulentes différentes) peuvent aboutir à un même gradient de température des protons dans le vent solaire.

Partie 4: Conclusions et travail futur

Dans cette dernière partie, nous exposons les conclusions de l'ensemble de nos résultats concernant l'anisotropie turbulente du vent solaire et le chauffage turbulent. Ces résultats ont tous été obtenus avec les équations de la MHD, qui supposent une température isotrope, qu'on ne trouve pas dans les vents rapides, dans lesquels la température parallèle diminue plus rapidement avec la distance héliosphérique que la température perpendiculaire. Nous proposons une modification à implémenter dans les équations EBM pour tenir compte de cet effet dans des travaux futurs. We introduce the context of our work in the first part. We give a summary of the physics concerning the solar wind and the solar corona, and the basic equations used to describe the solar wind plasma and an introduction to turbulence. Part 2 is dedicated to the study of anisotropy in the turbulent cascade, which characterizes 3D spectra. In the inertial range, in-situ measurements at 1 AU show complex figures for these spectra that we can interpret in several ways : numerical simulations allow to clear ambiguities. An important question is to know whether the Earth-Sun symmetry axis or the mean magnetic field axis is dominant. The third part focuses on turbulent heating in fast and slow winds. Between 0.3 and 1 AU, proton temperature decreases more slowly than expected, which requires a heating source. This source is supposed to be the continuous dissipation of eddies and current layers transported by the wind. To start with, we consider the simple case of Burgers equation, which is a one-dimensional model for shock formation. Thereupon, we switch to the 3-dimensional case, where we consider initial conditions appropriate for slow and fast winds.

In the last part we expose our conclusions and propose the implementation of temperature anisotropy as future work.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Velocity distribution functions of protons in the Solar Wind at 1AU (top row) and 0.3 AU (Bottom row). Left panels correspond to slow winds and right panels to fast winds. The dashed line corresponds to the mean magnetic field axis. From: [60]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Model profiles of the solar wind speed (U ) and the Alfvén wave speed (V a ) with distance from the Sun. The vertical bar separates the source, or sub-Alfvénic region, from the super-Alfvénic flow. "Previous Missions" marks the region explored by the Helios mission. From [26]. . . . . . . . . . . . . . . . 2.1 Sketch of the evolution of a plasma volume advected by a radial wind with constant speed. (a) exact evolution; (b) approximate evolution adopted in the EBM model. From [33]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Introducing turbulence. Left: Oblique light reveals turbulent eddies forming spontaneously above a warm cup of tea (courtesy of R. Grappin). Right: Shock steepening in 1D compressible Navier-Stokes equations ([8]). (a) Shock formation (formation of a steep gradient starting with a smooth cosine profile; (b) growth of successive harmonics with time; (c) successive energy spectra; (d) energy decay with time of the first Fourier mode (dotted) and of the total kinetic energy (solid line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Schematics of a local triad interaction (top) and a non-local one (bottom) and the corresponding position for each mode in the inertial range of a well developed turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Turbulent cascade vs linear expansion. Two possible evolutions of an initial isotropic spectrum (a), in two limiting cases: (b) standard cascade perpendicular to the mean magnetic field with negligible expansion; (c) no cascade, with only expansion. (From [25]). . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Measurements of magnetic energy spectrum at 1AU. From: [47]. . . . . . . . 4.2 Spectral indices at 1AU for magnetic energy fluctuations (red dots), kinetic energy fluctuations (blue dots), their sum (black dots), and their difference (green dots). From: [16]. Left panel: distribution of spectral indices; Right panel: spectral indices vs wind bulk speed. . . . . . . . . . . . . . . . . . . . 4.3 Left: autocorrelations of the magnetic field fluctuations computed at 1AU

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 First 118 days of Helios 1 mission starting from 1975, during solar minimum. (a) Synoptic view of daily averages. Solid line: bulk speed; dotted line: thermal proton velocity (arbitrary units); thick solid line: heliocentric distance (between 0.3 and 1 au); (b) Scatter plot of bulk velocity vs proton temperature T p ; (c) Scatter plot of bulk velocity vs compensated temperature RT p . . . 4.5 Same Helios 1 subset as in fig 4.4. Synoptic view of the daily fluctuations of energy densities eV and eB of u (dotted lines), b = δB/ √ ρ fields (solid lines)

. 1 .= 1 ,= 1 ,

 111 Last row: R13: solid line, R4: dotted, R6: dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Varying σ 0 c : ISO initial spectrum, R x = 1, m 0 = -1, = 0.2, spectral evolution of runs R2 (σ 0 c = 0, first row) and R4 (σ 0 c = 0.8, second row); same caption as in fig. 8.1. Left panel of the bottom row also shows the temporal evolution of the normalized cross helicity. Last row: R2: solid line, R4: dotted line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4 Varying σ 0 c : GYRO symmetry, R x = 1, m 0 = -1, = 0.2, spectral evolution of runs R7 (σ 0 c = 0, first row) and R8 (σ 0 c = 0.8, second row); same caption as in fig. 8.3; last row: R7: solid line, R8: dotted line. . . . . . . . . . . . . . 8.5 Varying σ 0 c , Gyro+Alfvén model, R x = 1, m 0 = -1, = 0.2, spectral evolution of runs R9 (σ 0 c = 0, first row) and R10 (σ 0 c = 0.8, second row): same caption as in fig. 8.3; last row: R9: solid line, R10: dotted line. . . . . . . . . 9.1 ISO initial conditions, summarizing how the final spectral main axis varies with and σ 0 c : passing from the GYRO structure perpendicular to the mean field (when = 0) to the R-SLAB (quasi-radial) structure (when both and σ 0 c are large). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Comparing runs with complete and incomplete initial isotropy, σ 0 c = 0, spectral evolution of runs R11 (R x = 5, k min x top row), R3 (R x = 1 and k min x second row), and R1 (R x = 5 and k min x = 0.2, third row). Same caption as in fig. 8.1. Last row: R11: solid line, R3: dotted, R1: dashed line. 10.1 Scheme of the different characteristic times for the Hydrodynamic expanding Box as fonction of the wave-number. The characteristic time of propagation of linear sound waves is denoted by t L , the characteristic time for non-linear interactions is defined by t N L and the characteristic expansion time by t exp . 10.2 Runs with monochromatic sinusoidal initial conditions: (a,b): RB1, = 0;

  (see text). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.3 Runs with non monochromatic initial conditions: (a,b): RB2, = 0; (c,d) RB4, = 100; (e,f) RB6, = 0.5. Same caption as previous figure. . . . . . . 12.1 Spectral evolution with initial ISO symmetry and Mach=1, varying σ 0 c : RC1 (σ 0 c = 0, first row) and RC2 (σ 0 c = 0.83, second row); same caption as in fig. 8.1; last row: RC1: solid line, RC2: dotted line. . . . . . . . . . . . . . . 12.2 Spectral evolution with initial Gyro-Alfvén symmetry and Mach=1, varying σ 0 c : RC5 (σ 0 c = 0, first row), RC7 (σ 0 c = 0.6, second row) and RC6 (σ 0 c = 0.83, third row); same caption as in fig. 8.1; last row: RC5: solid line, RC7: dotted line, RC6: dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.3 Turbulent heating. First column: Mach number vs distance; Second column: σ c vs distance, thicker lines correspond to positive values and thiner lines to the absolute value of negative ones; third column: normalized turbulent heating Q/Q c ; last column: normalized temperature profile RT (R). Oblique straight line corresponds to the adiabatic cooling prediction T ∝ R -4/3 . Top: initial ISO geometry (RC1: solid line, RC2: dotted); Middle row: GYRO geometry (RC3: solid, RC4: dotted); Bottom row: Gyro-Alfvén geometry (RC5: solid, RC7: dotted, RC6: dashed ). . . . . . . . . . . . . . . . . . . . 13.1 Helios mission, first four months of solar minimum; (a) rms turbulent Mach number M = u rms /c (c = (5/3 k B T p /m) 1/2 ); (b) expansion parameter (eq. 2.24). M and are computed at two frequencies: f= 3.3 10 -5 Hz (solid bold line) and 5.5 10 -5 Hz (dotted line). Both (a) and (b) show the bulk velocity U (plain line, arbitrary units); Panel (a) shows the heliocentric distance R × 10 [AU] (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 2 :

 12 Figure 1.2: (Left) Eclipse white-light image of the corona [35]. (Right) Coronal magnetic field-line configuration on 6 February 2004, as derived by applying a PFSS (Potential Field Source Surface) extrapolation with source surface at 2.5R s to a photospheric magnetic map from Mount Wilson observatory [106].

Figure 1 . 3 :

 13 Figure 1.3: Three Ulysses spacecraft orbits of the Sun. Figure shows radial solar wind velocity and images of the Sun at varying degrees of solar activity [67].

Figure 1 . 5 :

 15 Figure 1.5: Velocity distribution functions of protons in the Solar Wind at 1AU (top row) and 0.3 AU (Bottom row). Left panels correspond to slow winds and right panels to fast winds. The dashed line corresponds to the mean magnetic field axis. From: [60].

Figure 1 . 6 :

 16 Figure 1.6: Model profiles of the solar wind speed (U ) and the Alfvén wave speed (V a ) with distance from the Sun. The vertical bar separates the source, or sub-Alfvénic region, from the super-Alfvénic flow. "Previous Missions" marks the region explored by the Helios mission. From [26].

Figure 2 . 1 :

 21 Figure 2.1: Sketch of the evolution of a plasma volume advected by a radial wind with constant speed. (a) exact evolution; (b) approximate evolution adopted in the EBM model.From[START_REF] Grappin | Nonlinear wave evolution in the expanding solar wind[END_REF].

3 )

 3 shows how the initial smooth (cosine) profile evolves finally into a shock with a large gradient in the center of the figure. Panel (b) shows the growth of harmonics of the successive Fourier modes with time. Panel (c) shows the corresponding evolution of the energy spectrum with time. Panel (d) shows finally the kinetic energy evolution: the energy in the first Fourier mode first decreases (solid line), followed by the decrease of the overall energy content (dotted line).

Figure 3 . 1 :

 31 Figure 3.1: Introducing turbulence. Left: Oblique light reveals turbulent eddies forming spontaneously above a warm cup of tea (courtesy of R. Grappin). Right: Shock steepening in 1D compressible Navier-Stokes equations ([8]). (a) Shock formation (formation of a steep gradient starting with a smooth cosine profile; (b) growth of successive harmonics with time; (c) successive energy spectra; (d) energy decay with time of the first Fourier mode (dotted) and of the total kinetic energy (solid line).

Figure 3 . 2 :

 32 Figure 3.2: Schematics of a local triad interaction (top) and a non-local one (bottom) and the corresponding position for each mode in the inertial range of a well developed turbulence.

Figure 3 . 3 :

 33 Figure 3.3: Turbulent cascade vs linear expansion. Two possible evolutions of an initial isotropic spectrum (a), in two limiting cases: (b) standard cascade perpendicular to the mean magnetic field with negligible expansion; (c) no cascade, with only expansion. (From [25]).

7

 7 in the slowest wind to -1.55 in the fastest wind (cf fig. 4.2, right panel))

Figure 4 . 2 :

 42 Figure 4.2: Spectral indices at 1AU for magnetic energy fluctuations (red dots), kinetic energy fluctuations (blue dots), their sum (black dots), and their difference (green dots). From: [16]. Left panel: distribution of spectral indices; Right panel: spectral indices vs wind bulk speed.

Figure 4 .

 4 3 summarizes this variation of anisotropy depending on wind speed and frequency.4.4 Proton temperature gradient and turbulent heating4.4.1 Turbulent amplitude and proton temperature variations in the inner heliosphereAs we shall see in a moment, measurements of the proton temperature gradient in the inner heliosphere, say (Helios mission) between 0.3 and 1 AU are somewhat controversial. To understand why, let us examine the left panel in fig.4.4, which shows the daily variations of bulk velocity and proton thermal velocity, during the first four months of Helios 1 mission.

Figure 4 . 4 :

 44 Figure 4.4: First 118 days of Helios 1 mission starting from 1975, during solar minimum. (a) Synoptic view of daily averages. Solid line: bulk speed; dotted line: thermal proton velocity (arbitrary units); thick solid line: heliocentric distance (between 0.3 and 1 au); (b) Scatter plot of bulk velocity vs proton temperature T p ; (c) Scatter plot of bulk velocity vs compensated temperature RT p .

Figure 4 . 5 :

 45 Figure 4.5: Same Helios 1 subset as in fig 4.4. Synoptic view of the daily fluctuations of energy densities eV and eB of u (dotted lines), b = δB/ √ ρ fields (solid lines) averaged in the octave around the two waveperiods of one hour (two upper curves) and 3 minutes (two lower curves). Thick line: proton thermal velocity in arbitrary units. From [32].

Figure 4 . 6 :

 46 Figure 4.6: Measurement of proton temperature as a function of radial distance with Helios 1 and 2 data for winds with mean speed over 600km s -1 . Solid line corresponds to the fitted function T p ≈ 2.5 • 10 5 (R/R 0 ) -0.74 K. From [37].

Figure 4 . 7 :

 47 Figure 4.7: Kolmogorov-like estimate of turbulent heating Q K41 (eq. 3.13) at frequency f = 10 -2 Hz (red symbols) and critical heating rate Q c (blue symbols), vs proton temperature at 1 AU ([100]). Solid lines for each data population correspond to fitted power laws with proton temperature.

Figure 4 . 8 :

 48 Figure 4.8: Helios mission, solar minimum. Kolmogorov-like estimate of turbulent heating Q K41 (eq. 3.13) evaluated at frequency f = 4 10 -3 Hz (crosses) and critical heating rate Q c (lines), vs proton temperature. From left to right panels: increasing distance range from 0.3 to 1 AU. This is to be compared with the previous figure4.[START_REF] Belcher | Large-amplitude Alfvén waves in the interplanetary medium, 2[END_REF] which gives a similar analysis using a much larger dataset at 1 AU. Vertical units: [10 6 J/(kgs)].

Figure 4 . 9 :

 49 Figure 4.9: Helios mission, solar minimum. Kolmogorov estimate Q K41 (eq. 3.13) of the turbulent heating, normalized by the critical heating Q c vs frequency. Curves show averages in four distance ranges from 0.3 to 1 AU. The distance interval [0.3, 1] is divided into four equal subsets starting with 0.3 AU. Solid line: closest to 0.3 AU; then dotted, dashed, with dotted-dashed for the interval closest to R = 1AU . Left: Low proton temperatures; right: high proton temperatures.

Figure 4 . 10 :

 410 Figure 4.10: Estimations of turbulent heating normalized by the critical heating (eq. 3.44)at 1AU for slow and fast winds. Turbulent heating is computed via third-order moments laws that have assumed an isotropic hydrodynamic turbulence (green triangles), isotropic MHD turbulence (blue squares) or anisotropic MHD turbulence (red circles). From[START_REF] Stawarz | THE TURBULENT CASCADE AND PROTON HEATING IN THE SOLAR WIND AT 1 AU[END_REF].

Figure 6 . 1 :

 61 Figure 6.1: Initial total energy spectra, 2D slice of the angular spectrum in the plane k z = 0. Left: Isotropic (ISO) initial spectrum; center: Gyrotropic (GYRO) spectrum; right: Gyro-Alfvén (G+A) spectrum.

Figure 7 . 1 :

 71 Figure 7.1: 2D slice of the angular spectrum in the plane k z = 0 at the end of an ISO simulation, i.e. E 3D (k x , k y , k z = 0). The solid straight line corresponds to the largest vector in modulus from the center to the isocontour of level 10 -8 , k maj , while the dotted straight line is the smaller vector in modulus from the center to the same isocontour, k min . The mean magnetic field axis (not showed) forms an angle of 45 • with respect to the radial at the end of all simulations.

Figure 7 . 2 :

 72 Figure 7.2: Temporal evolution of the 1D total energy spectrum in the x direction for an ISO simulation. Each thin solid line correspond to the reduced spectra compensated by k 5/3 and averaged over one fifth of duration of the simulation. The arrow indicates the temporal evolution for the spectra: first time interval on top and the last one at the bottom. The solid thick line below the spectra helps to signal the length of the inertial range at the last time interval. Values in the vertical axis are arbitrarily chosen to help visualize the spectral evolution and are not representative of the real decrease of the amplitude of the fluctuations.

Figure 8 . 1 :

 81 Figure 8.1: Varying : ISO initial spectrum, σ 0 c = 0, spectral evolution of runs R12 ( = 0, first row), R2 ( = 0.2, second row), R5 ( = 0.4, third row). Rows 1 to 3: left panel, 1D energy spectra averaged in five time intervals, compensated by k -5/3 , vs k x ; mid panel, same vs k z ; right panel, 2D slices of the final 3D spectra E 3D (k x , k y , k z = 0). Last row: spectral aspect ratio A(k) and angle α(k) between main axis and radial (R12: solid line, R2: dotted, R5: dashed).

Figure 8 . 2 :

 82 Figure 8.2: Varying : ISO initial spectrum, σ 0 c = 0.8, R x = 1, spectral evolution of runs R13 ( = 0, first row), R4 ( = 0.2, second row), R6 ( = 0.4, third row). Same caption as in fig. 8.1. Last row: R13: solid line, R4: dotted, R6: dashed line.

Figure 8 . 3 :

 83 Figure 8.3: Varying σ 0 c : ISO initial spectrum, R x = 1, m 0 = -1, = 0.2, spectral evolution of runs R2 (σ 0 c = 0, first row) and R4 (σ 0 c = 0.8, second row); same caption as in fig. 8.1. Left panel of the bottom row also shows the temporal evolution of the normalized cross helicity. Last row: R2: solid line, R4: dotted line.

Figure 8 . 4 :

 84 Figure 8.4: Varying σ 0 c : GYRO symmetry, R x = 1, m 0 = -1, = 0.2, spectral evolution of runs R7 (σ 0 c = 0, first row) and R8 (σ 0 c = 0.8, second row); same caption as in fig. 8.3; last row: R7: solid line, R8: dotted line.

Figure 8 . 5 :

 85 Figure 8.5: Varying σ 0 c , Gyro+Alfvén model, R x = 1, m 0 = -1, = 0.2, spectral evolution of runs R9 (σ 0 c = 0, first row) and R10 (σ 0 c = 0.8, second row): same caption as in fig. 8.3; last row: R9: solid line, R10: dotted line.

Figure 9 . 2 := 1 ,= 1 ,

 9211 Figure 9.2: Comparing runs with complete and incomplete initial isotropy, σ 0 c = 0, spectral evolution of runs R11 (R x = 5, k min x = 1, top row), R3 (R x = 1 and k min x = 1, second row), and R1 (R x = 5 and k min x = 0.2, third row). Same caption as in fig. 8.1. Last row: R11: solid line, R3: dotted, R1: dashed line.

4 )( 10 . 5 )Figure 10 . 1 :

 4105101 Figure 10.1: Scheme of the different characteristic times for the Hydrodynamic expanding Box as fonction of the wave-number. The characteristic time of propagation of linear sound waves is denoted by t L , the characteristic time for non-linear interactions is defined by t N L and the characteristic expansion time by t exp .

Figure 10 . 2 :

 102 Figure 10.2: Runs with monochromatic sinusoidal initial conditions: (a,b): RB1, = 0; (c,d) RB3, = 100; (e,f) RB5, = 0.5. The top row (a,c,e) shows the evolution of the energy spectra compensated by k -2 , line styles vary according to the increase in time, from a thin solid line at t = 0, followed be the dotted and dashed line and so on. The bottom row (b,d,f) shows the evolution of the dissipation terms: total dissipation dE/dt (solid line), viscous dissipation ν(∂ x u) 2 (dotted line), damping due to the expansion 2n a u 2 (dashed line), prediction for turbulent dissipation k 0 u 3 /a α (dot-dashed line). The turbulent dissipation always matches the viscous dissipation term after a nonlinear time (see text).

Figure 10 . 3 :

 103 Figure 10.3: Runs with non monochromatic initial conditions: (a,b): RB2, = 0; (c,d) RB4, = 100; (e,f) RB6, = 0.5. Same caption as previous figure.

(

  Grappin et al. 1993; Grappin & Velli 1996; Dong et al. 2014) that allow us to follow the turbulent evolution transported by the radial wind. Let us denote by U e

Figure 1 .

 1 Figure 1. Initial and final domains of simulation (and plasma volume) in the ecliptic plane. Thin lines: direction of mean magnetic field. For all runs, the aspect ratio of the domain varies from 1/5 to unity. In the figure, the mean magnetic field angle with the radial varies varies from tan 1 5 11 . 3 1

. 3 .

 3 Critical Heating, Cascade Rate, and Parameter M 2 

4

 4 us to express Equation[START_REF] Hennebelle | Turbulent molecular clouds[END_REF] in terms of quantities based on rms velocity amplitude, namely the Mach number M and expansion parameter ò. One finds M

Figure 2 .

 2 Run A: evolution of basic quantities vs. heliocentric distance R R 0 . (a) Velocity amplitude u rms (solid line), compressible velocity u c rms (dotted line), and magnetic field fluctuation b rms 1 2 r ¯(dashed line). (b) Visco-resistive dissipation Q ν (solid line) decomposed as the sum of u 2 m  | ˜| (dotted line), dot-dashed line). (c) Dissipation rates per unit mass: dE/dt (solid line), expansion-driven damping Q exp (dotted line), viscoresistive dissipation Q ν (dashed line), and nonlinear loss during cascade Q NL (dot-dashed line, thick when increasing the decay rate and thin when decreasing the decay rate). (d) Temperature compensated by 1/R decrease. Distance is normalized by the initial distance R 0 =0.2 au.

  , we compare runs A (solid line), B (dotted line), C (dashed line), and E (dot-dashed line) with the following values of m and k max : run C, and m=2.2 and k 4 max = for run E. Panel (a) shows the critical heating ratio Q Q c n that in principle reveals how close we are from critical heating; panel (b) shows the temperature curve compensated for a 1/R decrease.

Figure 3 ..

 3 Figure 3. Heating ratio Q Q c n and temperature profiles, all with M 1, = 0.2  = but varying small-scale initial excitation. Runs A (solid line), B (dotted line), C (dashed line), and E (dot-dashed line). (a) Heating ratio Q Q c n vs. heliospheric distance R. (b) Average temperature (normalized by its initial value) compensated by R R 0 . Distance is normalized by the initial distance R 0 = 0.2 au. The thin solid line in panel (b) corresponds to T T R R 0 0 4 3 =

Figure 4 .

 4 Figure 4. Heating ratio Q Q c n and temperature profiles, with same 0.2  = but varying initial Mach number, for runs C (M=1; solid line) and D (M=0.77; dotted lines). Same caption as Figure 3.

Figure 5 .

 5 Figure 5. Heating ratio Q Q c n and temperature profiles, with same M=0.77 but varying expansion parameter, for runs F ( 0.12  = ; solid line) and K ( 0.2  = ; dotted line). Same caption as Figure 3.

Figure 6 .

 6 Figure 6. Heating ratio Q Q c n and temperature profiles, with varying M 2  , for runs E (solid line), F (dotted line), G (dashed line), and H (dot-dashed line). Same caption as Figure 3.

Figure 7 .

 7 Figure 7. Heating ratio Q Q c n and temperature profiles, with M=0.77, 0.2  = , and varying β, for runs M (solid line; β=2.75), N (dotted line; β=1.48), and K (dashed line; β=0.49). Same caption as Figure 3.

Figure 8 .

 8 Figure 8. Run E: rotation of the energy spectrum (kinetic + magnetic) with the mean magnetic field. The angular energy spectrum E(k x , k y ) is in the plane k z =0. (a) At R=0.2 au. (b) At R=1 au. The mean magnetic field direction is represented by a straight line in each panel.

Figure 9 . 3 -

 93 Figure 9. Runs A, B, C, and E: evolution of the 1D energy spectrum (kinetic + magnetic), compensated by the Kolmogorov scaling k 5 3 -. Abscissa: wavenumber normalized by the Kolmogorov wavenumber k Q d 3 1 4 m = n ( ) .

Figure 10 . 3 -

 103 Figure 10. Runs Z (solid line) with M=0.3 and A (dotted line) with M=1. The total 1D energy spectrum is averaged during the last phase of transport, compensated by the k 5 3 scaling. Abscissa: radial wavenumber normalized by the Kolmogorov wavenumber.

Figure 11 .

 11 Figure 11. Checking the Vasquez law and consequences. (a) Ratio Q Q K 41 n between visco-resistive dissipation and Kolmogorov's energy cascade rate (Equation (36)). (b) Radial Taylor wavenumber. (c) Ratio Q Q exp n between visco-resistive dissipation and expansion decay rate. (d) Turbulent rms amplitude Z u B 2 2 12 d r = + ( ) evolution compensated by R 1 0.6 . Shown are runs A (solid line), B (dotted line), C (dashed line), and E (dot-dashed line).

Figure 12 .

 12 Figure 12. Run A: deviation from total (turbulent + internal) energy conservation during cascade. Solid line: time derivative of total energy (turbulent + internal, Equation (48)). Dot-dashed line (thick when increasing the decay rate, thin when decreasing the decay rate): residual term Q NL¢ measuring the deviation from total energy conservation during the cascade apart from linear effects.

Figure 12 . 1 :

 121 Figure 12.1: Spectral evolution with initial ISO symmetry and Mach=1, varying σ 0 c : RC1 (σ 0 c = 0, first row) and RC2 (σ 0 c = 0.83, second row); same caption as in fig. 8.1; last row: RC1: solid line, RC2: dotted line.

Figure 12 . 2 :

 122 Figure 12.2: Spectral evolution with initial Gyro-Alfvén symmetry and Mach=1, varying σ 0 c : RC5 (σ 0 c = 0, first row), RC7 (σ 0 c = 0.6, second row) and RC6 (σ 0 c = 0.83, third row); same caption as in fig. 8.1; last row: RC5: solid line, RC7: dotted line, RC6: dashed line.

Figure 13 . 1 :

 131 Figure 13.1: Helios mission, first four months of solar minimum; (a) rms turbulent Mach number M = u rms /c (c = (5/3 k B T p /m) 1/2 ); (b) expansion parameter (eq. 2.24). M and are computed at two frequencies: f= 3.3 10 -5 Hz (solid bold line) and 5.5 10 -5 Hz (dotted line). Both (a) and (b) show the bulk velocity U (plain line, arbitrary units); Panel (a) shows the heliocentric distance R × 10 [AU] (dashed line).

U 0 ,

 0 U sw radial component of mean solar wind speed T, T p proton temperature T e electron temperature T α alpha particles (He +2 ) temperature T , T p, proton temperature parallel with respect to the mean magnetic field T ⊥ , T p,⊥ proton temperature perpendicular with respect to the mean magnetic field T e, , T e,⊥ electron temperature parallel and perpendicular to the mean magnetic field P = ρT isotropic pressure P component of the pressure tensor parallel to the mean magnetic field P ⊥ component of the pressure tensor perpendicular to the mean magnetic field γ polytropic index E electric field B magnetic field B 0 mean magnetic field δB magnetic field perturbation V a = B 0 / √ ρ 0 Alfvén speed b = δB/ √ ρ 0 z ± = u ± b Alfvén modes written in Elsasser variablesσ c = ((z + ) 2 -(z -) 2 )/((z + ) 2 + (z -) 2 ) normalized cross helicityσ 0 c initial value for the normalized cross helicity µ, η, κ dynamic viscosity, magnetic and thermal conductivity ν = µ/ρ kinematic viscosity L integral length of a system t ν characteristic time of viscous dissipation. In HD t ν = 1/(νL -2 ) t N L characteristic time of nonlinear interaction. In HD t N L = 1/(L -1 u) Re = t ν /t N L Reynolds number ũk

  dt s = ∆x/c s (B.5) dt A = ∆x/v A (B.6)dt diss = ∆x/max(µ/min(ρ), η, κ) (B.7)As a consequence, in the CFL condition becomes dt < min(dt s , dt A , dt diss ) (B.8)

Titre:

  Dynamique Turbulente du Vent Solaire Mots cl és : Turbulence, Magn étohydrodynamique (MHD), Vent Solaire R ésum é : Le but de cette th èse est l' étude du d éveloppement de la turbulence dans le vent solaire entre 0.2 et 1 unit é astronomique (UA) du soleil (i.e. l'orbite terrestre). L' étude est faite en r ésolvant num ériquement les équations de la MHD apr ès soustraction de l' écoulement moyen radial. Les deux aspects de la turbulence qui nous int éressent sont la structure 3D des spectres d' énergie et le chauffage du plasma qui r ésulte de la dissipation turbulente des tourbillons et couches de courant emport és par le vent. On cherche à d éterminer quelles sont les conditions du plasma pr ès du soleil qui permettent d'aboutir à ce qu'on observe à 1 UA. Un but important de mon travail est aussi de d éterminer si la physique qui est pr ésente dans les équations que j'int ègre (la MHD) suffit pour arriver à reproduire ce qu'on a d éj à observ é dans cet intervalle de distance. Nous introduisons le contexte de notre travail dans la premier partie. On y trouve les équations de base, une introduction à la turbulence, un r ésum é sur la physique du vent solaire et de la couronne solaire. La partie 2 sera consacr ée à l' étude de l'anisotropie de la cascade turbulente, et plus pr écis ément du spectre 3D. Dans la zone inertielle, les mesures in-situ vers 1 UA montrent des figures complexes pour ces spectres qu'on peut interpr éter de plusieurs fac ¸ons : nos simulations num ériques permettent de lever toute ambiguït é. Plus pr écis ément, la question est de savoir quand intervient l'axe soleil-terre, et quand intervient l'axe du champ magn étique moyen. La partie trois est centr ée sur le chauffage turbulent dans les vents rapides et lents. Entre 0.3 et 1 UA, la temp érature des protons diminue anormalement lentement, ce qui indique une source de chauffage, qu'on suppose ici être la dissipation des tourbillons et couches de courant emport és par le vent. Pour d émontrer que cette hypoth èse est raisonnable, nous consid érons d'abord le mod èle de Burgers qui est un mod èle pour l' évolution d'ondes sonores. Ensuite, nous passons à l' étude du cas plus complexe d'un volume de plasma 3D. Nous examinerons les conditions initiales correspondant aux vents lents et rapides. Dans les deux cas, on adoptera des anisotropies spectrales diff érentes. Dans la derni ère partie, nous exposerons les conclusions de notre travail et proposerons d'introduire l'anisotropie de la temp érature dans un travail futur. Title : Turbulent Dynamics of the Solar Wind Keywords : Turbulence, Magnetohydrodynamics, Solar Wind Abstract : The aim of this thesis is the study of the development of turbulence in the solar wind between 0.2 and 1 astronomical unit (AU) from the Sun (i.e. Earth's orbit). The study is done by solving the magnetohydrodynamics equations (MHD) after subtracting the mean radial flow. The two aspects of turbulence that interest us are the 3D structure of the energy spectra and the heating of plasma that results from the turbulent dissipation of eddies and current layers transported by the wind. We want to determine which conditions of the plasma close to the Sun can result into what we observe at 1 AU. We have relatively detailed measurements of what happens between 0.3 and 1 AU. One important goal of this work is to determine if the physics present in the equations that are integrated (MHD) is sufficient to reproduce what is observed in this interval of distances.

  13 Discussion Numerical parameters vs Helios data . . . . . . . . . . . . . . . . . . . Obtention of Maltese Cross components at 1AU . . . . . . . . . . . . . . . . . 14.2 Can the Maltese Cross heat the solar wind? . . . . . . . . . . . . . . . . . . . 14.3 Expectations from Solar Orbiter and Parker Solar Probe . . . . . . . . . . . . 14.4 Open questions and (partial) answers . . . . . . . . . . . . . . . . . . . . . . .

	IV Conclusions and future work
	14 Conclusions
	14.1 15 Future work: Anisotropy temperature description
	A Liste of symbols
	B Numerical method

B.1 Computation of spatial gradients . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Time integration method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  6 J/(kgs)]. . . . . . . . . . . . . . . . . . 4.9 Helios mission, solar minimum.Kolmogorov estimate Q K41 (eq. 3.13) of the turbulent heating, normalized by the critical heating Q c vs frequency. Curves show averages in four distance ranges from 0.3 to 1 AU. The distance interval [0.3, 1] is divided into four equal subsets starting with 0.3 AU. Solid line: closest to 0.3 AU; then dotted, dashed, with dotted-dashed for the interval closest to R = 1AU . Left: Low proton temperatures; right: high proton temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Estimations of turbulent heating normalized by the critical heating (eq. 3.44)

Table 1 .

 1 1: Average properties of fast and slow solar winds measured at 1AU made by

		Slow winds Observations Fast winds Observations
	Bulk velocity (km s -1 )	320	667
	Number density (cm -3 )	5.4	3
	Proton Temperature (K)	4.8 • 10 4	2.8 • 10 4
	Electron Temperature (K)	1.1 • 10 5	1.3 • 10 5
	Proton anisotropy T p /T p⊥	3.4	1.2
	Electron anisotropy T e /T e⊥	1.2	1.2

Table 6 .

 6 1: List of parameters for the initial conditions. Run: name of run; Type: class of initial conditions (see figure6.1); R x : initial aspect ratio of the numerical domain; k max : the maximum wavenumber in directions perpendicular to radial; k min x : the minimum wavenumber in the x direction (k min

	y	=k min z

Table 1

 1 List of Parameters for the Initial Conditions

	R	M	ò	M 2 	B 0	β	k max	m	m	0
	A	1	0.2	5	2.04	0.29	64	5/3	2.4 10 -3
	B	1	0.2	5	2.04	0.29	64	3	1.7 10 -3
	C	1	0.2	5	2.04	0.29	4	5/3	2.1 10 -3
	D	0.77	0.2	3	2.04	0.49	4	5/3	1.8 10 -3
	E	1	0.2	5	2.04	0.29	4	2.2	1.7 10 -3
	F	0.77	0.12	5	2.04	0.49	4	2.2	1.8 10 -3
	G	0.6	0.14	2.6	2.04	0.8	4	2.2	1.5 10 -3
	H	1	0.4	2.5	2.04	0.29	4	2.2	1.3 10 -3
	K	0.77	0.2	3	2.04	0.49	4	2.2	1.8 10 -3
	M	0.77	0.2	3	0.86	2.75	4	2.2	1.8 10 -3
	N	0.77	0.2	3	1.17	1.48	4	2.2	1.8 10 -3
	Z	0.3	0.2	0.45	2.04	0.29	64	5/3	1.5 10 -4

Table 12 .

 12 1: List of parameters for the initial conditions. Run: name of run; Type: class of initial conditions (see figure6.1); R x : initial aspect ratio of the numerical domain; k max : the maximum wavenumber in directions perpendicular to radial; m 0 : initial slope of the reduced spectrum; σ c : initial cross-helicity; : initial value of the expansion parameter; M = u rms /c s , with c s the sound; ν 0 : initial value of the kinematic viscosity (equal to η and κ).

The frequency of the fluctuations should be below the giration frequency of the protons, fp = (2π) -1 |e||B|/mp, and space fluctuations larger than the Larmor radius of the protons, ρL = mp|u ⊥ |/|e||B|, where u ⊥ is the modulus of the velocity components perpendicular to the mean magnetic field line B

The Knudsen number is a dimensionless parameter that indicates the "collisionality" for a population of particles. λ mf p being the mean free path of particles and ω col the frequency of collision between particles, the Knudsen number can be defined as Kn = kλ mf p or as Kn = ω/ω col , to characterize the collisionality respectively at the spatial scale k -1 or at the frequency ω.

The collisional age is defined as the number of collisional timescale for a given species in the characteristic expansion time for the wind. See for instance[START_REF] Maruca | Collisional Thermalization of Hydrogen and Helium in Solar-Wind Plasma[END_REF] for more details on the definition.

Wave-action is defined as the energy of a wave over the frequency of the wave, E/ω. In the case of sound waves, E = ρδu

and ω = kcs = ω0cs/(U0 + cs). Note that in the absolute frame of reference, the sound waves transported in a fluid with mean velocity U0 have a frequency ω0 = k(U0 + cs).The conservation of wave-action states ∂t(E/ω) + ∇ • ((U0 + cs)(E/ω)) = 0 (10.9)

The opposite happens for the transverse propagating waves, as the WKB theory predicts a decrease of Mach number as M = δu/cs ∝ R -1/12 . In the non-expanding simulations compressibility also decreases in time, as velocity fluctuations are dissipated by the viscous terms

The Astrophysical Journal, 853:153 (10pp), 2018 February 1 https://doi.org/10.3847/1538-4357/aaa1ea © 2018. The American Astronomical Society. All rights reserved.

The Astrophysical Journal, 853:153 (10pp), 2018 February 1 Montagud-Camps, Grappin, & Verdini

Remerciements

Chapter 11

Paper ApJ 2018: "Turbulent Heating between 0.2 and 1 au: A Numerical Study"

Part IV

Conclusions and future work

Rs Solar radii (1Rs = 6.96

Numerical method

B.1 Computation of spatial gradients

The numerical resolution of the EBM equations uses a pseudo-spectral method to solve the spatial gradients. For physical problems where it is possible to assume periodic conditions at the boundaries of the numerical domain, pseudo-spectral methods can be easily implemented. This is one of the reasons, along with the simplicity of the equations, to consider a cubic box instead of a shell section as the numerical domain for our simulations. This approximation done by the EBM description is acceptable for a box small enough, that is,

). Pseudo spectral methods are called in this way because one computes gradients in Fourier space . Non-linear terms are computed in real space.

Our code uses the algorithm FFTW (Fast Fourier Transformation of the West) to compute the Fourier and the inverse Fourier transformation of variables. This algorithm is specially efficient with arrays whose dimension can be expressed as a power of 2. Because of this, the resolution of the numerical domain is chosen as a power law of 2, typically taking a uniform mesh of 512 3 points to discretize the numerical domain.

In the absence of sharp discontinuities (e.g. Heavyside-like functions), the evaluation of the spatial gradients by pseudo-spectral methods is highly accurate, reducing the numerical errors to the order of the number representation in the machine. Because of that, energy conservation can be assured in the simulation or, more precisely, the numerical dissipation can be neglected in comparison with the contribution of the dissipative terms in the equations. In the particular case of the EBM equations, in the absence of expansion, the sum of the kinetic and magnetic field energy of the fluctuations and the internal energy is conserved.

When the expansion is non zero, then total energy is not conserved in the plasma volume, but this is physically correct and not due to the numerical method.

Despite of the advantages of the pseudo-spectral methods, their use has some limitations. One of the drawbacks of these methods comes from the need to represent the variables in a discretize Fourier space. Let u(x) be a function on R. It can be represented in Fourier space over a discretized grid of N points, j = 0, . . . , N -1, as a sum of trigonometric polynomials: where the k-mode coefficients, ũk , are related as follows to the k-modes of the continuous Fourier transformation, ûk ,