Giovanni Manfredi

Ozgür D Gürcan

Thierry Réveillé

Pierre Bertrand

Frédéric Brochard

Jean-Hugues Chatenet

Daniele Del Sarto

Erwan Deriaz

Stéphane Devaux

Thomas Drouot

Eric Faudot

Alain Ghizzo

Stéphane Heuraux

Anne-Sophie Latreille

Nicolas Lemoine

Jérome Moritz

Alex Aurelia

Chris

Diego

Greg

Huyen

Kathleen

Miguel Marion

Philippe, Sarah, Seb, Thibaud Vin- Cent Pierre

Merci À Abderzak

Jean-Loïs Dimitri

d'assister à ma soutenance. Ce fût aussi un plaisir de discuter avec eux et de répondre à leurs questions pertinentes.

Introduction "I would like nuclear fusion to become a practical power source. It would provide an inexhaustible supply of energy, without pollution or global warming." -Stephen Hawking

Global warming mitigation is the challenge of the XXI st century [START_REF] Cox | erratum: Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model[END_REF][START_REF] Thomas | Extinction risk from climate change[END_REF][START_REF] Pounds | Widespread amphibian extinctions from epidemic disease driven by global warming[END_REF][START_REF] Heimann | Terrestrial ecosystem carbon dynamics and climate feedbacks[END_REF][START_REF] Peters | The challenge to keep global warming below 2 c[END_REF][START_REF] Revesz | Global warming: Improve economic models of climate change[END_REF][START_REF] Berg | Landatmosphere feedbacks amplify aridity increase over land under global warming[END_REF][START_REF] King | Australian climate extremes at 1.5 c and 2 c of global warming[END_REF][START_REF] Russo | Half a degree and rapid socioeconomic development matter for heatwave risk[END_REF]. An international effort has been undertaken to limit the global warming below +2 • C in 2100, compared to pre-industrial era in the 1850s [START_REF] Rogelj | Paris agreement climate proposals need a boost to keep warming well below 2 c[END_REF]. However, it is likely that the +2 • C mark will be exceeded before 2100 [START_REF] Raftery | Less than 2 c warming by 2100 unlikely[END_REF][START_REF] Tollefson | Ipcc says limiting global warming to 1.5 c will require drastic action[END_REF]. Indeed, current predictions indicate a +3 • C to +5 • C warming in 2100 [START_REF] Knutti | A scientific critique of the two-degree climate change target[END_REF]. Humaninduced warming already reached +0.93 • C (±0.13 • C) in 2015 compared to the 1850s [START_REF] Millar | Emission budgets and pathways consistent with limiting warming to 1.5 c[END_REF][START_REF] Allen | Warming caused by cumulative carbon emissions towards the trillionth tonne[END_REF][START_REF] Fischer | Anthropogenic contribution to global occurrence of heavy-precipitation and high-temperature extremes[END_REF]. Recent climate studies have demonstrated the relation between average global temperature rises and cumulative emissions of greenhouse gases [START_REF] Mcglade | The geographical distribution of fossil fuels unused when limiting global warming to 2 c[END_REF][START_REF] Matthews | The proportionality of global warming to cumulative carbon emissions[END_REF]. For the +2 • C target, CO 2 emissions need to be reduced by 50 to 80% in 2050, compared to the 2010 emissions [START_REF] Knutti | A scientific critique of the two-degree climate change target[END_REF] while a +1.5 • C target requires reaching net Instead of breaking a heavy atomic nucleus into lighter parts like in nuclear fission, nuclear fusion combines two light nucleus into one to produce energy. The problematic of nuclear fusion is entirely different than nuclear fission, and is much harder to achieve on Earth. Nuclear fusion has been occurring in Nature for billions of years at the core of stars. At the center of the Sun, hydrogen nuclei, or protons, combine to create helium nuclei which are composed of 2 protons and 2 neutrons. More than 80 different types of fusion reactions between light nuclei have been observed and studied [START_REF] Rax | Physique des tokamaks[END_REF]. The most accessible one on Earth is the fusion between deuterium and tritium (D-T), two hydrogen isotopes, because it has the highest cross-section/energy ratio (see Fig. 1). Deuterium is composed of one proton and one neutron while tritium has one proton and two neutrons. This reaction produces a helium nucleus (an α particle) and a neutron,

D + T → α (3.52MeV) + n (14.06MeV) (1)
This reaction does not produce long-lasting radioactive product like in fission reactions. Also, the reaction presents zero risk of "runaway" or "nuclear meltdown". Deuterium is abundant on Earth as it is present in the ocean water: in every cubic meter of seawater, there are 33 grams of deuterium. At current rate of energy consumption, there is enough deuterium on Earth for 10 11 years.

Tritium is radioactive and has a half-life of ∼ 12 years, thus it is not present in nature and has to be produced. Luckily, it is possible to produce tritium in an elegant way using the following nuclear reactions [25, p. 26]:

n + 6 Li → α + T + (4.80MeV) (2)
and n

+ 7 Li → α + T -(2.47MeV) (3)
which uses the neutron from the nuclear reaction (1), combined with Lithium, which is available on Earth [25, p. 26].

Nuclear fusion is the promise of near unlimited energy, in vast quantity, with no CO 2 emission, with no long lasting nuclear waste and no possibility of "runaway" reaction.

Then why is nuclear fusion still not a commercial power source on Earth?

In order to fuse, the deuterium and tritium need to overcome the potential energy barrier originating from the fact that they are both positively charged (Coulomb repulsion). In order to achieve this, the positively charged particles need to reach extremely high velocity, or, equivalently, high energy. Fig. 1 shows the cross section of different fusion reactions. We see that for energies < 200 keV, the D-T reaction has, by one or two order of magnitude, the machine exploiting this concept is the tokamak, an acronym in Russian meaning "toroidal chamber with magnetic coils". A tokamak is a complex machine which mainly features a vacuum chamber in toroidal shape where the plasma is present, and magnetic coils generating a confining magnetic field of 1 -10 Tesla. To reach energies > 10 keV, the plasma needs to be heated. This is routinely done by a set of combined methods and the main ones are: ohmic heating, neutral beam injection and resonant absorption of radio frequency electromagnetic waves [26, p. 239].

For this approach to be worthwhile, the power produced from the nuclear fusion needs to be greater than the power used to heat the plasma. This leads to the Lawson criterion [27] [26, p. 11] which states:

nT τ E > 10 21 m -3 .keV.s [START_REF] Heimann | Terrestrial ecosystem carbon dynamics and climate feedbacks[END_REF] where n is the particle density, T is the plasma temperature in energy units (10-20 keV) and τ E is the energy confinement time, which is the time the energy dissipates from the plasma if all external heat sources are stopped.

In this triple product, the density is limited in tokamaks to a maximum of about 10 20 particles per cubic meter because of the Greenwald limit [21, p. 500]. The temperature is limited by the β limit where β is the ratio between kinetic pressure P k = nk B T and the magnetic pressure B 2 /2µ 0 . This ratio in tokamaks is usually about ∼ 1% and an increase in β can lead to reduced confinement. Thus, in current tokamaks, only the energy confinement time can be increased.

Throughout the end of the 20th century, numerous tokamaks were built all around the world, each one gathering precious data and gradually improving the triple product.

The Russian tokamak T-1 from the Kurchatov Institute in Moscow, which started operating at the end of 1958, is usually considered to be the first tokamak. It had a magnetic field of 1.5 Tesla, a plasma volume of 0.4m 3 inside a copper vacuum vessel and was far from achieving the Lawson criterion. The first controlled deuterium-tritium fusion, with megawatt-range power, was achieved on November 9 th , 1991, in the Joint European Torus (JET) in Culham (UK) [START_REF] Team | Fusion energy production from a deuterium-tritium plasma in the jet tokamak[END_REF], with a fusion amplification factor Q DT ≈ 0.15 (the ratio between the power injected in the tokamak and the power produced), and produced 1 MW of fusion energy over 2 s. In 1997, JET achieved a fusion power of 16.1 MW, with a Q DT ≈ 0.67, which is still the actual record of Q [START_REF] Keilhacker | High fusion performance from deuterium-tritium plasmas in jet[END_REF]. In its current configuration, JET has a plasma volume of 90m 3 . From all these data, an empirical law has been found, linking the energy confinement time τ E to the size of the tokamak. Fig. 2 shows the measured τ E in tokamaks, against the predicted τ E from the scaling law. This scaling law proved to be particularly accurate for the moment.

Figure 2: Measured τ E against predicted τ E from the scaling law, in log-log scale, for different tokamaks. Picture taken from [START_REF]Energy confinement time scaling[END_REF].

Thus, extrapolating this law, it is possible to estimate the required tokamak size to achieve a fusion amplification factor Q = 10, which would need τ E = 5 s. This is the goal of ITER, which results from an international cooperation between the European Union, USA, China, Russia, India, Japan and South Korea. It represents a leap in terms of size, and τ E compared to previous tokamaks. One goal is to achieve 500 MW of fusion power from D-T fusion with only 50 MW of power injected (Q = 10) for 400 to 600 s. ITER will have a plasma volume of 830m 3 and a magnetic field of 5.3 T and a core temperature of 13 keV.

It is currently under construction in Cadarache, France. It is scheduled to start its first plasma in 2025, and start the deuterium-tritium fusion experiments after 2035. Its total cost is currently estimated to be about 20 billions euros, which needs to be put in perspective knowing that ITER members represent half the population and 60% of the total wealth in the world.

If ITER is a success, its successor, DEMO, would be operational after 2050 and would serve as a demonstration that fusion power can produce energy at industrial level. DEMO would produce at least 1 GW of fusion power.

Nevertheless it is not certain whether or not the scaling law will saturate at some point because of physical processes amplified at higher tokamak sizes. Indeed there are multiple mechanisms that tend to reduce τ E . The main one is the transport of particles and heat from the center of the plasma (the core) toward the wall of the tokamak (the edge), across the magnetic field lines. Understanding transport in order to control it and increase τ E , is a subject of uttermost importance in the domain of magnetic fusion confinement. An early theory of transport, called neo-classical theory is an improvement of classical theory (which was based on the kinetic theory of gases, and modeled the transport as a random walk of particles caused by collisions [START_REF] Balescu | Transport processes in plasmas[END_REF]) as it takes into account the curvature of the magnetic field. Nevertheless neo-classical transport predictions are usually one order of magnitude lower than actual transport in tokamaks.

The difference between the predicted transport and the actual transport is called "anomalous transport".

Anomalous transport is now understood as being dominated by turbulent transport. Turbulence is a highly nonlinear process involving multi-scale phenomena in space and time [START_REF] Kadomtsev | Plasma turbulence[END_REF][START_REF] Davidson | Turbulence: an introduction for scientists and engineers[END_REF][START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF].

Comprehension of tokamak turbulence made tremendous progress in the last few decades thanks to a wealth of analytical, experimental and numerical results. Anomalous transport originates from microinstabilities (microscopic compared to the size of the device) in the plasma, generating turbulent transport [START_REF] Freidberg | Plasma physics and fusion energy[END_REF][START_REF] Wesson | [END_REF][START_REF] Garbet | Physics of transport in tokamaks[END_REF][START_REF] Doyle | Plasma confinement and transport[END_REF][START_REF] Yoshizawa | Plasma and Fluid Turbulence: Theory and Modelling[END_REF][START_REF] Balescu | Aspects of anomalous transport in plasmas[END_REF][START_REF] Diamond | Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas[END_REF]. They lead to electric potential fluctuations φ in the plasma and electric fields E. In the presence of an electric field and a magnetic field, the charged particles undergo the so called E × B drift, leading to the apparition of drift waves [START_REF] Horton | Drift waves and transport[END_REF][START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]. These electromagnetic (considered electrostatic when β 1) waves can be destabilized by density and temperature gradients. In this case, the electric potential φ grows in amplitude until it reaches a saturation level. The most common microinstabilities are the Ion Temperature Gradient (ITG) driven modes, the Electron Temperature Gradient (ETG) driven modes, the Trapped Electron Mode (TEM) and the Trapped Ion Mode (TIM) which belong to the family of drift waves. Plasma microinstabilities can also give rise to large scale structures such as zonal flows [START_REF] Lin | Turbulent transport reduction by zonal flows: Massively parallel simulations[END_REF][START_REF] Rogers | Generation and stability of zonal flows in ion-temperature-gradient mode turbulence[END_REF][START_REF] Kim | Dynamics of zonal flow saturation in strong collisionless drift wave turbulence[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF] or streamers [START_REF] Dimits | Formation of streamers in plasma with an ion temperature gradient[END_REF][START_REF] Champeaux | Streamer and zonal flow generation from envelope modulations in drift wave turbulence[END_REF][START_REF] Yamada | Anatomy of plasma turbulence[END_REF].

In order to investigate turbulent transport in tokamaks, we need to model the plasma, which can follow different approaches, for example:

• Describing the particles motion directly with Newton equation (6N dimensions), coupled to Maxwell equations with microscopic fields,

• Using a statistical, kinetic description of the plasma via a distribution function f (6 dimensions), coupled to Maxwell equations with macroscopic fields,

• Using a statistical fluid description (3 dimensions), considering the evolution of averaged quantities such as the density or the pressure.

Although the first approach is complete, it is unachievable considering that in ITER there will be of the order of 10 23 particles. The fluid approach requires, to be valid, that particles have a Gaussian distribution in velocities which is not necessarily the case in a tokamak plasma. Also, the resonant interactions between waves and particles are not fully taken into account.

We adopt the kinetic approach in this work, because we want to describe phenomena which tend to impact the velocity distribution. In this approach, the time evolution of the distribution function is governed by the Vlasov equation coupled to Maxwell equations with macroscopic fields. These coupled nonlinear equations present no analytical solutions for non trivial cases.

With the recent development of High Performance Computing (HPC), it is possible to obtain approximate solutions to these equations on supercomputers, using a reduced statistical approach. Indeed, the 6D distribution function evolution coupled to Maxwell equations proves too heavy to solve in general, even with state-of-theart supercomputers. The problem can be reduced from 6D to 5D (4 variables + one adiabatic invariant), by averaging the fast cyclotron motion of the particles. This is the concept of gyrokinetic theory [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF]. This approach has lead to a class of codes called gyrokinetic codes. Standard state-of-the-art gyrokinetic codes nonexhaustively include GYSELA [START_REF] Grandgirard | A drift-kinetic semi-lagrangian 4d code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Global full-f gyrokinetic simulations of plasma turbulence[END_REF][START_REF] Darmet | Intermittency in flux driven kinetic simulations of trapped ion turbulence[END_REF][START_REF] Dif-Pradalier | Interplay between gyrokinetic turbulence, flows, and collisions: Perspectives on transport and poloidal rotation[END_REF][START_REF] Sarazin | Large scale dynamics in flux driven gyrokinetic turbulence[END_REF], GENE [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF][START_REF] Dorland | Electron temperature gradient turbulence[END_REF][START_REF] Goerler | The global version of the gyrokinetic turbulence code gene[END_REF], GKW [START_REF] Peeters | The nonlinear gyro-kinetic flux tube code gkw[END_REF], ELMFIRE [START_REF] Heikkinen | Gyrokinetic simulation of particle and heat transport in the presence of wide orbits and strong profile variations in the edge plasma[END_REF], ORB5 [START_REF] Jolliet | A global collisionless pic code in magnetic coordinates[END_REF], GT5D [START_REF] Idomura | Conservative global gyrokinetic toroidal full-f five-dimensional vlasov simulation[END_REF][START_REF] Idomura | Study of ion turbulent transport and profile formations using global gyrokinetic full-f vlasov simulation[END_REF] and GYRO [START_REF] Candy | An eulerian gyrokinetic-maxwell solver[END_REF][START_REF] Candy | Anomalous transport scaling in the diii-d tokamak matched by supercomputer simulation[END_REF]. Each of these gyrokinetic codes focuses on particular ingredient of plasma turbulent transport. Indeed, there is no "super code" integrating every known ingredients and physical processes of the turbulence as it would require orders of magnitude more computational power. For example, some codes rely on an approximation to neglect the perturbation of the magnetic field induced by the particles. It is the electrostatic assumption, which allows to only describe the evolution of the electric field, using the Poisson equation. Even with this focused approach, these gyrokinetic codes are still very computationally intensive.

In a tokamak, the magnetic field intensity is higher toward the center of the torus than toward the edge. Thus there is a potential energy barrier as particles approach toward the high magnetic field side. This leads to a separation of the particles into two classes: the passing, and the trapped particles. The passing particles fully follow the winded, looped magnetic field lines, as they overcome the potential energy barrier toward the center of the torus. The trapped particles lack energy to pass the high field side potential barrier, thus at some point their parallel (to the magnetic field) velocity becomes 0 and they are repelled in the opposite direction. They are thus trapped on the outer region of the tokamak, the low field side, and their motion resembles the outline of a banana. This motion of trapped particles is called bounce motion. It has a lower frequency than the cyclotron frequency, and a larger spatial scale than the cyclotron radius, by at least one order of magnitude. The trapped particles ("bananas") slowly (compared to the bounce frequency) precess in the toroidal direction.

In TERESA [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF][START_REF] Ghizzo | Streamer-induced transport in the presence of trapped ion modes in tokamak plasmas[END_REF][START_REF] Cartier-Michaud | Optimizing the parallel scheme of the poisson solver for the reduced kinetic code TERESA[END_REF][START_REF] Cartier-Michaud | Staircase temperature profiles and plasma transport self-organisation in a minimum kinetic model of turbulence based on the trapped ion mode instability[END_REF][START_REF] Drouot | Étude de la turbulence liée aux particules piégées dans les plasmas de fusion[END_REF][START_REF] Médina | Radial density and heat fluxes description in the velocity space: Nonlinear simulations and quasi-linear calculations[END_REF], the trapped particle bounce motion is averaged out. Moreover, the passing particles respond adiabatically to the electric potential, thus, the passing particles density simply depends linearly on the electric potential. Thus, only turbulence on time-scale of the slow precession motion of the trapped particle is retained. The trapped particles, which are de facto "banana centers", are described using the Vlasov equation coupled to the Poisson equation. Focusing on turbulence driven by trapped particles is relevant because these modes play important role in transport [START_REF] Horton | Drift waves and transport[END_REF]. In TERESA, only Trapped Ion Mode and Trapped Electron Mode turbulence can appear.

The TERESA code underwent multiple improvements. It was first a 3D code (2 variables + one adiabatic invariant) with adiabatic trapped electrons developed in Nancy, in 2000 [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF], based on the 1999 Semi-Lagragian approach [START_REF] Sonnendrücker | The semilagrangian method for the numerical resolution of the vlasov equation[END_REF]. An adiabatic invariant was added in 2010 [START_REF] Ghizzo | Streamer-induced transport in the presence of trapped ion modes in tokamak plasmas[END_REF] therefore the code became 4D. In 2014, the code was rewritten at the CEA Cadarache to improve parallel computing efficiency. In 2015, kinetic trapped electrons were added, at IJL Nancy, with the possibility of N species kinetic description [START_REF] Drouot | Étude de la turbulence liée aux particules piégées dans les plasmas de fusion[END_REF][START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF][START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF][START_REF] Drouot | Self-generated zonal flows in the plasma turbulence driven by trapped-ion and trapped-electron instabilities[END_REF].

Part of the work presented in this thesis was to develop a test particle module to the TERESA code. Indeed, solving the distribution function from Vlasov-Poisson does not directly yield the particle trajectories. Knowing the particles trajectories leads to have better insights on various transport processes, such as diffusion [START_REF] Rechester | Electron heat transport in a tokamak with destroyed magnetic surfaces[END_REF][START_REF] Isichenko | Invariant measure and turbulent pinch in tokamaks[END_REF][START_REF] Harvey | Radial transport and electron-cyclotron-current drive in the tcv and diii-d tokamaks[END_REF][START_REF] Estève | Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport[END_REF], hyper-or sub-diffusion [START_REF] Misguich | Noble internal transport barriers and radial subdiffusion of toroidal magnetic lines[END_REF], advection [START_REF] Isichenko | Invariant measure and turbulent pinch in tokamaks[END_REF][START_REF] Boedo | Flow reversal, convection, and modeling in the diii-d divertor[END_REF][START_REF] Boedo | Transport by intermittent convection in the boundary of the diii-d tokamak[END_REF][START_REF] Dif-Pradalier | On the validity of the local diffusive paradigm in turbulent plasma transport[END_REF], the trapping of particles in potential wells, and ballistic events such as avalanches [START_REF] Sarazin | Large scale dynamics in flux driven gyrokinetic turbulence[END_REF][START_REF] Idomura | Study of ion turbulent transport and profile formations using global gyrokinetic full-f vlasov simulation[END_REF][START_REF] Dif-Pradalier | On the validity of the local diffusive paradigm in turbulent plasma transport[END_REF][START_REF] Sarazin | Predictions on heat transport and plasma rotation from global gyrokinetic simulations[END_REF].

The other part of the work was focused on analyzing the transport using the test particles, as well as to study the particles and heat fluxes with high accuracy in energy space, and to compare these fluxes to quasi-linear predictions.

In chapter I, we introduce the magnetic configuration of a tokamak, along with single charged particle trajectories in the magnetic field. We then describe the model of the TERESA code in section II.1 and we detail the implementation of test particles in the TERESA code in section II.2. The test particle trajectories are computed with the same order of accuracy as the Semi-Lagrangian solver for the Vlasov-Poisson system. In chapter III we use the newly developed test particle module to separate diffusive from advective contribution to the particle flux. For this we use millions of test particles and let them evolve in a turbulent plasma simulation. Studying their evolution in time, we calculate a random walk diffusion coefficient in energy space. We then estimate a diffusive flux and compare it to the total particle flux which takes into account diffusive and non-diffusive processes. We can thus estimate the importance of diffusion in the total transport processes.

In chapter IV, we study the particle and heat fluxes with a high accuracy in the energy dimension. The fluxes are dominated by a narrow resonance peak. We then proceed to compare the results with quasi-linear theory. We find that quasilinear predictions are in qualitative agreement with the nonlinear simulations from TERESA, but with almost a ∼ 50% discrepancy in the peak magnitude. A term by term study is made and we found that in our case, some neglected terms by quasi-linear theory are in fact of the same order of magnitude as other kept terms.

Chapter I

Magnetic confinement

In this chapter we start by describing the magnetic field B inside a tokamak, in section I.1. Its analytical expression in 3D space is not trivial as the magnetic field lines are winded around a torus. Nevertheless we make the approximation of axisymmetric magnetic field, so that magnetic surfaces become nested tori. It is worth to note that B can be described as a combination between a magnetic flux linked to the radial component, and two angles. The tokamak magnetic configuration is chosen to limit the particle drift, which degrades the confinement. Indeed the charged particle trajectories are non trivial as their velocity can be decomposed in multiple components: a cyclotron velocity, a parallel velocity (parallel to the magnetic field lines) and a drift velocity which we explicit in section I.2. Depending on the ratio between parallel velocity and perpendicular velocity, a charged particle can be passing or trapped. Trapped particles are confined to the low field side of the tokamak (toward the outside) and slowly (compared to other motions) drift in the toroidal direction. In section I.3 we introduce the statistical description of the plasma, using a distribution function which time evolution is governed by the Vlasov equation, coupled to Maxwell equations for the electromagnetic field, or only the Poisson equation if we make the electrostatic assumption. Then in section I.4 we introduce the concept of transport and why it is a challenge in tokamaks, then we recall that transport is mainly governed by turbulence driven by microinstabilities.

I.1 Magnetic configuration

The motion of a single particle of charge e thrown with velocity v in a uniform and time-independant magnetic field B 0 will endure the Laplace force F = e v × B 0 and 11 its trajectory can be described exactly by solving Newton's second law of motion

m d v dt = e v × B 0 . (I.1)
The particle will have a helical trajectory twisting around the magnetic field lines (see Fig. I.1). Its velocity can be decomposed into 2 components: one along the magnetic field lines v and one in the plane perpendicular to the magnetic field lines v ⊥ . Thus, the velocity can be written v = v + v ⊥ with the parallel velocity v constant. The perpendicular velocity is called the cyclotron velocity and can be written

v ⊥ = v ⊥ = v c = v c = ρ c ω c
with ω c the angular frequency called the cyclotron frequency and ρ c the cyclotron radius or Larmor radius that we will explicit later. In order to confine a population of charged particles, one could thus think to fold back the magnetic lines on themselves to make a loop, in a toroidal machine, a "donut". This simple approach is undermined because of the existence of a curvature drift v CD which occurs when a charged particle follows a curved magnetic line and slowly pushes the charged particle toward the roof or toward the floor, depending on the particle's charge, which in turn leads to charge separation and the creation of an electric field in the perpendicular plane. Furthermore, in this configuration, the magnetic field is more intense toward the center of the torus than toward the edge. This magnetic field gradient in the plane perpendicular to the magnetic field line also creates a drift of the particle called the ∇B drift.

Any deviation from the assumptions of an uniform magnetic field with straight field lines and no electric field leads to the apparition of drifts in the perpendicular plane to B or parallel accelerations. Later we will see that these parallel accelerations to B, or forces, will lead to the apparition of a class of particles which are "trapped" [26, p. 44]. Other main drift sources are the presence of an electric field perpendicular to B which leads to a drift called E × B drift, or the time dependency of this electric field which leads to a polarization drift. To compensate the curvature and ∇B drifts, the idea is to induce a poloidal field in addition to the toroidal field, so that the charged particle will now tend to follow a twisted magnetic field line, see Fig.

I.1.1 General idea of a tokamak magnetic field

The tokamak's magnetic field B can be seen as a sum of two components :

B = B θ e θ + B ϕ e ϕ (I.2)
where B θ and B ϕ are respectively the poloidal and the toroidal magnetic field and e θ and e ϕ are the unit vectors in the poloidal and toroidal directions respectively, see Fig. The toroidal field B ϕ is not homogeneous inside the torus. A very simple way to realize that is to use Ampere's law [26, p. 106]. It results that B ϕ is inversely proportional to R,

B ϕ = µ 0 N I 2πR (I.3)
where I is the current in the coil from Fig. I.4, N is the number of coils, µ 0 the vacuum permeability and R = R 0 + r cos θ. Thus,

B ϕ = B 0 1 + r R 0 cos θ (I.4)
On the other hand, the poloidal field B θ is mainly generated by an electrical current circulating in the plasma, in the toroidal direction. Its magnitude is usually an order of magnitude weaker than the toroidal magnetic field. The total magnetic field being helical, it is possible to define a safety factor q which measures the local ratio between the number of toroidal turns a magnetic field line has to do before making one poloidal turn,

q = B. ∇ϕ B. ∇θ (I.5)
which can be approximated [26, p. 112] in the case of an axisymmetric, large aspect ratio ξ = R 0 a system by :

q(r) ≈ r R B ϕ B θ (I.6)
In ITER, the safety factor at the edge (r = a) of the chamber will be 3 [25, p. 35] [87], meaning a magnetic field line has to do 3 toroidal turns before making one poloidal turn, and the aspect ratio will be ξ = 6.2 2 ≈ 3. Now let us find a way to describe B.

I.1.2 Axisymmetric magnetic field description

We define a set of general toroidal coordinates [88, p. 6] (ρ(x, y, z), η(x, y, z), ζ(x, y, z)) which are functions of the euclidean set of coordinates (x, y, z), where surfaces of constant ρ designate the surfaces of nested tori around the axis ρ = 0, η designates the angular position around this axis and ζ the angular position around the toroidal axis. Any tokamak magnetic field B can then be written [25, p. 163] as

B = ∇φ T × ∇η + ∇ζ × ∇φ P (I.7)
where φ T = φ T (ρ, η, ζ) and φ P = φ P (ρ, η, ζ) are respectively the toroidal and poloidal flux functions. If the field is axisymmetric and thus presents perfectly nested magnetic surfaces, the flux functions φ P and φ T are linked to the (internal) poloidal and toroidal magnetic fluxes [25, p. 165] as φ P = Φ Pi 2π and φ T = Φ T 2π where Φ Pi and Φ T depends only on ρ and are defined as :

Φ Pi (ρ) = S Pi (ρ,η) B. dS (I.8) Φ T (ρ) = S T (ρ,ζ) B. dS (I.9)
where S Pi designates the surface delimited by the magnetic axis ρ = 0 and the magnetic surface of constant ρ, and S T is the surface of a torus section at constant ζ.

The magnetic field in this particular axisymmetric case can thus be written [25, p. 169] [89, p. 121]:

B(Φ Pi , η, ζ) = 1 2π ∇Φ Pi (ρ) × ∇(q(ρ)η -ζ) (I.10)
where in this case, ρ = r, ξ = θ and ζ = ϕ with r, θ and ϕ respectively the standard radius, poloidal angle and toroidal angle in torus coordinates (see Fig . I.3). We also pose ψ = Φ Pi 2π so that in definitive :

B(ψ, θ, ϕ) = ∇ψ(r) × ∇(q(r)θ -ϕ) (I.11)
The poloidal flux normalized to 2π, ψ can be linked to the radial coordinate r [65, 66] [70, p. 17]. Indeed the internal poloidal surface vector can be written dS = -R(r)dϕdr e θ , and, B. dS = -B θ (r)R(r)dϕdr. Thus :

Φ Pi = 2πψ = - 2π φ=0 r r=0 B θ (r )R(r )dϕ dr = -2π r r=0 B θ (r )R(r )dr (I.12)
and therefore

ψ(r, θ) = -B 0 r r=0 r dr q(r ) 1 + r R 0 cos(θ) (I.13)
In the case of a high aspect ratio r R 0 1 we thus have

ψ(r) ≈ -B 0 r r=0 r dr q(r ) (I.14)
And when we neglect the radial dependence of the safety factor q, we obtain

ψ ∝ -r 2 (I.15)
In the TERESA code detailed in II.1, the poloidal flux ψ will be used as the radial variable. Now that we understand the magnetic configuration in a tokamak, let us see how a single charged particle moves inside this field.

I.2 Single charged particle trajectories

We saw that a single charged particle in a uniform time-independant magnetic field will have a helical trajectory around the field line with a velocity v = v + v ⊥ . But in a real tokamak the exact equations of motion of a particle of species s (s can be electron or ion s = e -, i) and charge e s = Z s e (Z s is the number of elementary charges of the species s, Z e -= -1 for electrons, and e is the elementary charge), are determined from Newton's second law of motion :

m s d v s dt = e s [ E( r s , t) + v s × B( r s , t)] (I.16) with d r s dt = v s (I.17)
with non uniform, time-dependent, electromagnetic fields, the single particle's velocity is rather decomposed as v = v ⊥ + v + v D with v D the vectorial sum of all particle drifts. As we will see, the single charged particle motion in a tokamak presents a clear separation of temporal and spatial scales.

I.2.1 Cyclotron motion

The cyclotron motion described in Chap. I.1 is still present even in more complex electromagnetic fields. This helical motion is characterized by its Larmor or cyclotron radius ρ c and frequency ω c defined respectively in Eqs. I.19 and I. [START_REF] Matthews | The proportionality of global warming to cumulative carbon emissions[END_REF]. We can evaluate these two quantities considering ITER's parameters [START_REF] Aymar | The iter design[END_REF] : an average plasma temperature T s of 150 millions Kelvin, a toroidal magnetic field In the following, we will always be under the adiabatic assumption.

ω cs = |e s |B m s (I.

I.2.1.1 The « guiding center » motion

Instead of describing the full cyclotron motion consisting of a fast rotation perpendicular to B superposed to a slower motion parallel to B motion, it is possible to average the motion over the cyclotron period, so that the position of the particle is only determined by its guiding center [25, p. 42] [90, p. 75] [21, p. 139] position r g . This is the gyro-cyclotron averaging. The position r s (t) of a single charged particle can thus be decomposed as

r s (t) = r gs (t) + ρ cs (t) (I.20)
where ρ cs (t) describes the fast rotating motion around the magnetic field line.

I.2.1.2 The magnetic moment

The fast rotating motion ρ cs (t) can be considered, on a time average, to be a current loop of intensity I = des dt = e s ωcs 2π . The current loop will create a magnetic field in the opposite direction of B. This phenomenon highlights that plasmas in magnetic fields are naturally diamagnetic. We thus introduce the magnetic moment µ so that

µ = µ = IS (I.21)
with S the surface enclosed by the current loop [91, p. 186]. Thus we have the magnetic moment of the species s :

µ s = 1 2 m s v 2 ⊥s B (I.22)
Under the adiabatic assumption, this magnetic moment is an ("the first") adiabatic invariant, which means it is an approximate constant of motion [25, p. 122

I.2.2 Motion parallel to B

In the introduction of this chapter it was stated that in a uniform, purely magnetostatic field, the particle's velocity along the magnetic field line, or parallel velocity v was constant, because Laplace's force only act in the plane perpendicular to B. In a tokamak, the charged particle can have a non-zero acceleration along B.

I.2.2.1 Acceleration due to an electric field parallel to B

If there is an electric field E with a parallel component (E ) to B, Newton's second law of motion along B simply becomes :

dv dt = e s m s E (I.23)
This parallel acceleration play a role in the generation of runaway electrons [START_REF] Rosenbluth | Theory for avalanche of runaway electrons in tokamaks[END_REF][START_REF] Helander | Runaway acceleration during magnetic reconnection in tokamaks[END_REF][START_REF] Eriksson | Simulation of runaway electrons during tokamak disruptions[END_REF] [21, p. 207] which are not in the scope of this manuscript.

I.2.2.2 Acceleration due to the non-uniformity of B

From I.3 we saw that the magnetic field in a tokamak was more intense toward the center of the torus than toward its edge. This fact is particularly important, because it divides the particle population into two different classes : the trapped particles and the passing particles.

A magnetic moment µ in a magnetic field B undergoes a force F = ∇( µ. B) [91, p. 189]. Thus there is a force applying to the magnetic moment introduced in I.22,

F = F b = -µ ∇ B (I.24)
where b = B B . It is to be noted that this force applies on the particle's guiding center, and not on its instantaneous position [21, p. 172]. Since µ s = The fact that µ s is an adiabatic invariant means that if B increases, v ⊥s must increase too. However, in a static magnetic field, the particle's kinetic energy [21, p. 173], so if v ⊥s increases, it means that v decreases. Eventually, if the increase in B is large enough, or equivalently, if the charged particle's parallel velocity over perpendicular velocity

1 2 m s (v s + v ⊥s ) 2 is an exact constant of motion
v v ⊥ ( v ⊥
v is called the pitch angle) is not high enough, the charged particle will eventually reach a point where v = 0 and will be reflected toward the zone with lower B . This is the concept of magnetic mirror in tokamak. This thus leads to two different populations of particles :

• The particles which have a pitch angle

|v s | |v ⊥s | < Bmax
B min -1, thus reaching a point called a mirror point at a critical poloidal angle θ crit where their parallel velocity v becomes zero, and are reflected toward the zone with lower B , are called the trapped particles [START_REF] Kadomtsev | Trapped particles in toroidal magnetic systems[END_REF], because they are trapped on the outer side of the tokamak and cannot reach the inner side.

• The particles which have a pitch angle

|v s | |v ⊥s | > Bmax
B min -1, so that their parallel velocity v is always finite ( = 0). These particles are called the passing particles because they can fully follow the helical magnetic field lines, making a whole poloidal turn, without being reflected toward the low B side.

I.2.2.3 Focus on the trapped particles population

The trapped particle population can represent an important part of the total particle population (trapped + passing particle population). The fraction of trapped particles can be written [25, p. 137] in the case of large aspect ratio ξ :

f T = 2r R 0 ∼ ξ -1 2 (I.25)
It is worth mentioning that the concept of trapped particles between magnetic mirrors formed the basis of mirror machines which are an earliest concept of fusion machine, which in experiment turned out to have a lot of particles lost in the walls of the machine [21, p. 175].

Banana trajectory of a trapped particle A priori, the trapped particle will thus approximately follow its magnetic field line, reach a mirror point at a critical poloidal angle θ crit , bounce back toward the low B side following its magnetic field line in the opposite direction, until it reaches the mirror point on the other end at a critical poloidal angle -θ crit . The guiding center of the trapped particle will thus describe a back and forth motion between 2 mirror points. This motion is called the bounce motion and is characterized by the bounce frequency ω bs which can be written [25, p. 135] under the assumption that v c v and that v D v c :

ω bs = µ s B 0 r m s q 2 R 3 0 ≈ ρ cs ω cs qR 0 r 2R 0 (I.26)
where B 0 is the intensity of the magnetic field at r = 0 and q is the safety factor. However, when reversing its trajectory at a mirror point (v = 0), the trapped particle will not exactly follow its previous path in the opposite direction. Thus, when looking at the projection of the trapped particle trajectory in a poloidal cross-section, one will not see a circle's arc on the low B zone, but rather a banana, with a certain width δ b at θ = 0 called banana width, see Fig. I.5, which can be written [25, p. 137]

δ bs = 2qρ cs R 0 r v s0 ω cs0 ≈ 2qρ cs 2R 0 r (I.27)
where q is the safety factor, v s0 and ω cs0 are the parallel velocity and cyclotron frequency at the poloidal angle θ = 0. For ITER's parameters [START_REF] Aymar | The iter design[END_REF], at r = a we obtain the numerical values : It is worth mentioning that the banana width is proportional to the Larmor radius.

For α particles which are highly energetic particles and typically have a Larmor radius more than one order of magnitude higher than the protons, the trapped orbit resembles more a potato than a banana [26, p. 131] [25, p. 142].

Precession in the toroidal direction of the banana trajectory Lastly, a banana orbit does not perfectly close on itself (see Fig. It is worth mentioning that this drift frequency does not depend on the particle's mass. The typical time and spatial scales of the cyclotron, bounce and toroidal drift motions, for ITER's parameters [START_REF] Aymar | The iter design[END_REF]

B ≈ B ϕ , v ⊥ = v ths = 2k B Ts ms
the thermal velocity of the species s, the safety factor q = 3, the minor radius a = 2 m and the major radius R 0 = 6.2 m.

I.2.3 Motions perpendicular to B called "drifts"

We saw that the velocity of a single charged particle in an electromagnetic field can be decomposed as follows:

v = v ⊥ + v D + v b (I.29)
This drift of the particle's guiding center at the velocity v D is slow compared to the parallel velocity v and the cyclotron velocity v ⊥ , and can be decomposed in multiple contributions. We will simply recall the main ones [21, p. 178] :

v D = v E + v ∇B + v curv + v p + ... (I.30) with v E the E × B drift v E = E ⊥ × B B 2 , (I.31) v ∇B the ∇B drift v ∇B = m s v 2 c e s B b × ∇B 2B , (I.32) v curv the curvature drift v curv = m s v 2 e s B b × ∇B B , (I.33)
and v p the polarization drift

v p = m s e s B 2 d E ⊥ dt . (I.34)
All of these drifts except the E × B drift create currents and charges separations because of their dependence on e s . As the diamagnetic drift is not a guiding center drift, it is not included in this list.

The ∇B drift and the curvature drift are called magnetic drifts and are directed upward or downward according to the sign of e s . The tokamak configuration cancels these drifts in the equilibrium.

We will now introduce a set of coordinates particularly adapted to the particle's motion in a tokamak.

I.2.4 Action-angle variables

When the motion is integrable, it is possible to build a set of canonical coordinates consisting of actions J and angles α in which the Hamiltonian of the system only depends of J. This is the case for the equilibrium motion in an axisymmetric tokamak: H eq = H eq (J 1 , J 2 , J 3 ). With the action-angle variables, Hamilton's equations thus become

d J dt = - ∂H eq ∂ α = 0 (I.35) d α dt = ∂H eq ∂ J ≡ ω(J 1 , J 2 , J 3 ) (I.36)
where J(t) = J 0 and α(t) = ωt + α 0 with J 0 = (J 1 , J 2 , J 3 ), α = (α 1 , α 2 , α 3 ) and ω = (ω 1 , ω 2 , ω 3 ). The actions J are conserved along the movement, while the angles α evolve linearly as a function of time.

For each periodic motion of the particles, an action is determined as J = 1 2π C p. dl, where C is the contour of the periodic motion, and is a constant of motion [25, p. 93]. The three quasi-periodic motions of a trapped particle have been explicited in I.2. They are the cyclotron motion, the bounce motion and the precession motion. We will now search the actions associated with these motions, knowing that p = m v + e A where A is the potential vector such that B = ∇ × A. Let us now find the actions of the trapped particle motion.

I.2.4.1 Action J 1s associated with the magnetic moment µ

The adiabatic invariant associated with the cyclotron motion can be linked to the magnetic moment µ introduced in I. This invariant is called the toroidal kinetic momentum and is usually written P ϕ .

In the case of trapped particles, we saw in I.2.2.3 that at a reflection point, v s = 0. Thus, at a reflection point, J 3s = e s ψ.

The three actions are thus

J =    J 1s = -ms es µ s J 2s = 1 2π C m s v. dl J 3s = e s ψ    (I.47)
and the associated frequency are simply the frequency associated to the quasiperiodic motions : [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF] where α 0 = 0. Now that we have described how a single particle moves in a tokamak magnetic field, let us see how we can describe the evolution of the whole plasma.

ω =    ω 1s = ω cs ω 2s = ω bs ω 3s = ω ds    (I.48) and α(t) =    α 1s = ω cs t α 2s = ω bs t α 3s = ω ds t    (I.

I.3 Kinetic description of the plasma I.3.1 The Vlasov-Maxwell system

In sections I.2 and I.1 we considered the Newton's second law of motion in order to describe a single charged particle trajectory. Coupled with Maxwell's equations it is possible to find the system's unknowns ( r j and v j ) for each particle j. In order to completely describe the microscopic state of the plasma at a specific time, one requires the 3 position and 3 velocity coordinates for each particle. However, one cubic meter of fusion plasma usually contains approximately 10 20 particles. The total number of particles in ITER's plasma will thus be around 10 23 . Describing the trajectory of each particles seems unnecessary and far out of reach even with recent advances in high performance computing (HPC). A more adapted approach to the problem is to use a statistical description of the system. Indeed, macroscopic measures like pressure, temperature, etc, rely very little on the microscopic behavior of single particles, but depend almost exclusively on the collective behavior of the particles. The plasma state at a time t can be described with a distribution function f s = f s ( r, v, t), which measures the probability to find a particle of species s in the infinitesimal phase space volume δ rδ v around ( r, v). In practice, δ r should be much smaller than the Debye length λ D but much greater than the mean distance between individual particles and δ v much smaller than the thermal velocity defined in I.2.1, or alternatively, it requires the graininess parameter g = 1 nλ 3 D → 0 which means that the plasma is considered to be a continuous medium governed by collective effects. The evolution of f s ( r, v, t) in the case of a collisionless, non-relativistic fusion plasma is described by the collisionless Boltzmann equation, also called Vlasov equation [START_REF] Vlasov | On the vibrational properties of an electronic gas[END_REF][START_REF] Vlasov | The vibrational properties of an electron gas[END_REF][START_REF] Hénon | Vlasov equation[END_REF]:

df s dt ( r, v, t) = ∂f s ∂t ( r, v, t) + v. ∂f s ∂ r ( r, v, t) + F m s . ∂f s ∂ v ( r, v, t) = 0 (I.50)
which can be written

df s dt ( r, v, t) = ∂f s ∂t + v. ∂f s ∂ r + e s m s [ E( r, t) + v × B( r, t)]. ∂f s ∂ v = 0 (I.51)
where the macroscopic fields E( r, t) and B( r, t) do not only contain the electromagnetic fields from external sources but also the (self-consistent) fields created from the moments of the distribution function itself. The charge density ρ s ,

ρ s ( r, t) = Z s e s f s ( r, v, t)d v (I.52)
and the current density J s ,

J s ( r, t) = Z s e s vf s ( r, v, t)d v (I.53)
both appears in the Maxwell's equations needed to close the system,

∇. E = s ρ s 0 (I.54) ∇. B = 0 (I.55) ∇ × E = - ∂ B ∂t (I.56) ∇ × B = µ 0   s J s + 0 ∂ E ∂t   (I.57)
with the electric field E = -∇φ -∂ A ∂t and the magnetic field B = ∇ × A, and where the sum s is on all the plasma species. The Vlasov equation is thus heavily coupled with the Maxwell's equations. This system is called the Vlasov-Maxwell system. The Vlasov equation states that f s is conserved along the trajectories dfs dt = 0 , or characteristics, defined by d r dt = v and d v dt = es ms [ E( r, t) + v × B( r, t)]. The TERESA code detailed in the next part (II.1) will take advantage of this property, which can be seen as a consequence of Liouville's theorem which states that a phase space volume advected along the characteristics is conserved.

I.3.2 The Vlasov-Poisson system

An approximation of the Vlasov-Maxwell system can be made in order to obtain a Vlasov-Poisson system. In the case of low β plasma,

β = kinetic pressure magnetic pressure = nk B T B 2 2µ 0 1, (I.58)
which is usually the case for tokamak's plasmas, we can make the assumption that the magnetic field varies slowly in time compared to the electric field, thus, ∂ B ∂t ≈ 0 and E ≈ -∇φ. It is the electrostatic assumption. Thus, only the Maxwell-Gauss equation is needed to close the system. Each density can be decomposed as follows :

∇. E =

n s = n eq,s + ñs (I.61)
where ñs is the density fluctuation of the species s and at equilibrium the sum of the ions densities equals the electron density thus we have n eq,e = n eq = N j=1 Z j n eq,j . Therefore we can rewrite the Poisson equation (I.60) as

λ 2 D ∆ eφ T eq,e = ñe n eq - 1 n eq N j=1 Z j ñj (I.62)
where we introduced the Debye length λ D = 0 Teq,e e 2 neq . In a space of canonical coordinates which we can call ( α, J), it is convenient to use a compact, Hamiltonian form of the Vlasov equation :

∂f s ∂t ( α, J) -[H s ( α, J), f s ( α, J)] ( α, J) = 0 (I.63) where [F, G] ( α, J) = 3 i=1 ∂F ∂α i ∂G ∂J i -∂G ∂α i
∂F ∂J i are the Poisson brackets and H s is the Hamiltonian of the system. This notation will be useful to average out the fast motions associated with the cyclotron and bounce frequencies, in the TERESA code described in II.1.

I.4 Transport

As we mentioned in the introduction, in order for the fusion plasma to reach ignition, the triple product nT τ E must be greater than a certain threshold. To reach this threshold, it is planned in ITER that the energy confinement time τ E will be of a few seconds, the core plasma density will be n core ≈ 10 20 m -3 and the temperature around T core ≈ 1.5 × 10 8 K. However, this density and temperature decrease rapidly toward the edge, by a few order of magnitude, so that temperature at the edge is around T edge ≈ 1×10 3 K, therefore creating extremely steep gradients of temperature, density and pressure. In this configuration, fluxes of particles and heat spontaneously appear in the radial direction. In order to have a self-sustaining fusion reactor, the heating by α-particles must balance the losses in energy toward the edge. Understanding the transport of particles and heat from the core plasma toward the edge in order to reduce it is a vast and active domain of research [START_REF] Garbet | Physics of transport in tokamaks[END_REF][START_REF] Doyle | Plasma confinement and transport[END_REF][START_REF] Balescu | Aspects of anomalous transport in plasmas[END_REF].

I.4.1 Classical and neo-classical theories

Considering these enormous gradients of temperature, density or pressure, it is thus tempting, in order to have a rough estimation of the density and heat fluxes, to use Fick's [START_REF] Fick | Ueber diffusion[END_REF] first law

Γ( r, t) = -D ∇n( r, t) (I.64)
where Γ is the density flux and D a diffusion coefficient, and Fourier's law

Q( r, t) = -nχ ∇T ( r, t) (I.65)
where Q is the heat flux and χ the heat diffusivity coefficient [START_REF] Garbet | Instabilités, turbulence et transport dans un plasma magnétisé[END_REF]. Classical theory and neo-classical theory give an estimation of these coefficients [START_REF] Balescu | Transport processes in plasmas[END_REF][START_REF] Balescu | Aspects of anomalous transport in plasmas[END_REF][START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity[END_REF]. These theories are based on the assumption that transport is mainly governed by binary Coulomb collisions between particles in a quiescent plasma. Classical theory gives a diffusion coefficient which corresponds to a random walk in a cylinder, D classical ∼ ρ 2 ce ν coll [START_REF] Balescu | Transport processes in plasmas[END_REF] where ρ ce is the electron Larmor radius and ν coll is the electron-ion collision frequency. When compared to the experimental coefficient, D classical appears to be several orders of magnitude lower than the one measured. A caveat of classical theory of transport is that is does not take into account the toroidal geometry of the magnetic field. Neo-classical theory of transport attempts to add these geometrical effects to classical theory, taking into account the magnetic field inhomogeneity and curvature. The neo-classical coefficient scales as D neo-classical ∼ q 2 ξ 3/2 D classical [START_REF] Balescu | Transport processes in plasmas[END_REF] with q the safety factor defined in Eq. (I.5) and ξ the aspect ratio. Although the neo-classical diffusion is larger than the classical one by two orders of magnitude, it is still far below the experimental measurements, see These theories are in fact crippled by the fact that collisions are not the dominant transport mechanisms in a tokamak and that a fusion plasma is not in a quiescent state.

I.4.2 Anomalous transport

Historically, the measured transport which was not explained by classical nor neoclassical theories has been named anomalous transport. Anomalous transport is now generally assumed to be dominated by a type of transport driven by highly non-linear turbulence processes occurring on multiple scales in space and time and originating from micro-instabilities [START_REF] Garbet | Physics of transport in tokamaks[END_REF][START_REF] Doyle | Plasma confinement and transport[END_REF][START_REF] Balescu | Aspects of anomalous transport in plasmas[END_REF][START_REF] Diamond | Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas[END_REF][START_REF] Horton | Drift waves and transport[END_REF][START_REF] Garbet | Turbulence in fusion plasmas: key issues and impact on transport modelling[END_REF]. Indeed, considering the collective nature of a fusion plasma, particles tend to organize themselves in more or less coherent structures, which can sometimes vastly enhance radial transport, or on the contrary, reduce it [START_REF] Lin | Turbulent transport reduction by zonal flows: Massively parallel simulations[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF][START_REF] Champeaux | Streamer and zonal flow generation from envelope modulations in drift wave turbulence[END_REF][START_REF] Yamada | Anatomy of plasma turbulence[END_REF]. Understanding anomalous, or turbulent transport is thus one of the top priority to achieve controlled thermonuclear fusion.

I.4.2.1 Turbulence

Turbulence can be defined as a physical nonlinear mechanism involving a fluid (or plasma) chaotic flow featuring a wide range of length and time scales, and possibly vortices [START_REF] Davidson | Turbulence: an introduction for scientists and engineers[END_REF][START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF].

3D hydrodynamic turbulence

In a 3D turbulent system, Kolmogorov 1941 theory (K41) [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF] assumes that vor- tices of size L = 2πk -1 (where k is a wavenumber) will transfer part of their kinetic energy to newly formed, smaller vortices. These smaller vortices will again transfer part of their kinetic energy to even small vortices and so on. Although the total kinetic energy of the system is conserved, the kinetic energy is transfered from the scale L = 2πk -1 to the smaller scales. This phenomenon is called direct cascade of energy and with E(k) the energy spectrum of the system we can write

E(k) ∝ k -5/3

2D hydrodynamic turbulence

Although, strictly speaking, 2D turbulence does not exist in nature, some physical systems can relevantly be described by 2D turbulence. It is the case for example of large scale turbulence of the atmosphere or the oceans because of their large ratio of lateral (∼ 10 4 km) to vertical length scales (∼ 10 km) [START_REF] Boffetta | Two-dimensional turbulence[END_REF][START_REF] Pandit | An overview of the statistical properties of two-dimensional turbulence in fluids with particles, conducting fluids, fluids with polymer additives, binary-fluid mixtures, and superfluids[END_REF]. As we will explain, turbulence in tokamaks can also be described as 2D [START_REF] Mckee | The beam emission spectroscopy diagnostic on the diii-d tokamak[END_REF]. In 2D turbulence, in addition to the kinetic energy, the enstrophy Ω is also conserved. The enstrophy is defined as the squared vorticity Ω = | ∇ × v| 2 [39, p. 10]. This results in a turbulent spectrum with the "double cascade" [39, p. 62] scenario of Kraichnan [START_REF] Kraichnan | Inertial-range spectrum of hydromagnetic turbulence[END_REF][START_REF] Kraichnan | Inertial-range transfer in two-and three-dimensional turbulence[END_REF] and Batchelor [START_REF] Batchelor | Computation of the energy spectrum in homogeneous twodimensional turbulence[END_REF]. When injecting kinetic energy at a scale L I = 2πk -1 I , the kinetic energy is now transfered to larger scales (smaller k < k I ) in an inverse cascade such that E(k) ∝ k -5/3 , while the enstrophy follows a direct cascade to the smaller scales (k

I < k) as E(k) ∝ k -3 , see Fig. I.9.

Turbulence in core tokamak plasma

In tokamaks, the turbulence is significantly different than in "3D turbulence". Indeed, the turbulence can mainly develop in a plane perpendicular to the magnetic field due to the strong magnetic field and the turbulence is thus similar to a 2D hydrodynamic turbulence with nevertheless key differences. In our code detailed in chapter II.1, the small dissipation scale is that of numerical dissipation, as the number of grid point is finite. Also, the spectral energy can be written

E k = kv 2 (k) where v(k) is the E × B (or B × ∇φ) drift velocity thus in Fourier space we have v(k) = kφ k and therefore E k = k 3 φ 2 k [65].
Core turbulence is considered to be triggered by microinstabilities, or microscopic instabilities. They are microscopic in the sense that they are of much smaller length scale than larger hydrodynamic, or MHD instabilities which typically happen at the equilibrium length scale. Microinstabilities are driven by density, temperature and sometimes velocity gradients [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF].

In the Vlasov-Poisson description of the plasma, only electrostatic waves are considered. An electrostatic plane wave is a purely longitudinal wave with no magnetic component. They can be written using complex form so that they are ∝ e i( k. r-ωt) with wave vector k and complex frequency ω = ω r + iγ, thus the wave can be written ∝ e γt e i( k. r-ωrt) . Therefore a mode is unstable when its growth rate is γ > 0.

The mode grows in amplitude and eventually reaches a saturation through diverse physical mechanisms.

A few of the most common unstable (considered electrostatic) modes in collisionless tokamak plasmas are:

• the Electron Temperature Gradient (ETG) -driven mode [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF][START_REF] Dorland | Electron temperature gradient turbulence[END_REF] with a typical scale ρ c,e λ ⊥ ρ c,i . At this scale, the ion response can be considered adiabatic, therefore, ETG driven modes do not induce significant ion particle, heat or momentum transport [START_REF] Waltz | Coupled ion temperature gradient and trapped electron mode to electron temperature gradient mode gyrokinetic simulations[END_REF].

• the Trapped Electron Mode (TEM) [START_REF] Dannert | Gyrokinetic simulation of collisionless trappedelectron mode turbulence[END_REF][START_REF] Terry | Nonlinear inward particle flux component in trapped electron mode turbulence[END_REF][START_REF] Lang | Nonlinear saturation of collisionless trapped electron mode turbulence: Zonal flows and zonal density[END_REF][START_REF] Vernay | Global gyrokinetic simulations of tem microturbulence[END_REF] which can overlap in scales with ITG and ETG-driven modes [START_REF] Garbet | Physics of transport in tokamaks[END_REF][START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF][START_REF] Arnichand | Identification of trapped electron modes in frequency fluctuation spectra[END_REF]. TEM microinstabilities are driven by resonance with trapped electrons.

• the Ion Temperature Gradient (ITG) -driven mode with a typical scale 0.1 < k ⊥ ρ c,i < 1 [START_REF] Garbet | Physics of transport in tokamaks[END_REF][START_REF] Arnichand | Identification of trapped electron modes in frequency fluctuation spectra[END_REF].

• the Trapped Ion Mode (TIM) which is analogous to the TEM but with large wavelength δ b,i < λ ⊥ /2π.

Table I.2 recalls the different scales in wavelength/wavenumber for the main microinstabilities.

In tokamaks, it is generally admitted that the low frequency modes (ω ω c ) like ITG driven modes and TEM dominate the transport [START_REF] Horton | Drift waves and transport[END_REF][START_REF] Tang | Microinstability theory in tokamaks[END_REF]. The TERESA code detailed in Chapter II.1 focuses on large spatial scale (λ ⊥ /2π > δ b,i ) and low frequency phenomena (ω < ω b ) and thus only allow Trapped Ion Mode and Trapped Electron Mode to develop.

Perpendicular wavelength Perpendicular wavenumber

ETG ρ c,e λ ⊥ ρ c,i 1 k ⊥ ρ c,i TEM overlaps ITG/ETG modes overlaps ITG/ETG modes ITG ρ c,i < λ ⊥ /2π < 10ρ c,i 0.1 < k ⊥ ρ c,i < 1 TIM δ b,i < λ ⊥ /2π k ⊥ δ b,i < 1
Table I.2: Summary of the different scales for the main microinstabilities. k ⊥ = 2π/λ ⊥ is the perpendicular microinstability wavenumber.

Trapped Ion Mode (TIM) and Trapped Electron Mode (TEM)

Trapped Ion Mode (TIM) [START_REF] Kadomtsev | Plasma instability due to particle trapping in a toroidal geometry[END_REF][START_REF] Kadomtsev | Turbulence in toroidal systems[END_REF][START_REF] Tagger | Trapped ion mode driven by ion magnetic drift resonance in a fat torus[END_REF][START_REF] Tang | Residual trapped-ion instabilities in tokamaks[END_REF][START_REF] Biglari | Toroidal ion-pressure-gradientdriven drift instabilities and transport revisited[END_REF][START_REF] Garbet | Ionic instability thresholds in tokamaks[END_REF][START_REF] Artun | Trapped ion mode in toroidally rotating plasmas[END_REF] are long wavelength (compared to the banana width δ b,i ), low frequency (compared to the ion bounce frequency ω b,i ), electrostatic, toroidal microinstabilities triggered by a radial ion temperature gradient ∇ ⊥ T i in the "unfavorable curvature region", similarly to the toroidal Ion Temperature Gradient (ITG) -driven instabilities [START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF]. TIMs are in some sense a type of ITG, considering they require ion pressure gradient, nevertheless we separate them from ITG, because they are driven through the resonant interaction between a wave and trapped ions at their precession frequency [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF].

The Trapped Ion Mode instability occurs because of a temperature gradient ∇ ⊥ T i and because of the ∇B drift, Eq. (I.32). This vertical drift depends on the particle's temperature, thus, when the temperature gradient is aligned with the magnetic field gradient, which is the case in a tokamak, particles down the temperature gradient drift slower than particles up the temperature gradient. If the temperature gradient were unperturbed (i.e. a perfect linear temperature gradient orthogonal to the ∇B drift velocity) the TIM instability would not occur. When the density gradient is perturbed, as in Fig. I.10, the difference in ∇B drift velocity creates charge separation at the interface of the perturbation. This charge separation creates an electric field E. On the low magnetic field side of the tokamak (the "unfavorable curvature region"), the ∇B is in the same direction as ∇ ⊥ T i .

The particles then undergo the E × B drift which amplifies the initial temperature gradient perturbation. It is worth mentioning that on the strong field side of the tokamak (the "favorable curvature region"), ∇B and ∇ ⊥ T i are in opposite direction. In this case, on Fig. I.10 the hotter region and the colder region would be swapped, and the electric field created from charge separation would be in the opposite direction, thus, in this case the E × B would suppress the perturbation instead of amplifying it. The TIM instability appears only above a certain threshold [START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF]. The linear threshold is determined (by linear analysis) by a parameter η i = d(log T i )/d(log n i ) [START_REF] Garbet | Ionic instability thresholds in tokamaks[END_REF] as well as the temperature ratio τ = T i /T e . Fig. I.11 shows an example of stability diagram of ITG and TEM. TEM and ITG (or TIM) appear only above a certain density and temperature threshold. In chapter III, we run numerical simulations of TIM/TEM by imposing an initial temperature and density gradient, so that the modes are out of the stable region.

I.4.2.2 The general "transport equation"

The general transport equation can be written as :

Q = F ( X ) (I.66)
where Q is a "supervector" grouping all the relevant fluxes [38, p. 10], X is grouping all the thermodynamic forces [128, p. 99], and F could be a nonlinear function.

The relevant fluxes in Q are usually the particles fluxes Γ s , the heat or energy fluxes q s , the electric current density J s or the dissipative pressure tensor π s (which is the traceless part of the total pressure tensor [31, p. 140]). These fluxes can be nonlinearly coupled to the thermodynamic forces such as ∇n( r, t) or ∇T ( r, t).

In order to have a tractable problem, it is generally assumed (and confirmed experimentally [START_REF] Moret | Tokamak transport phenomenology and plasma dynamic response[END_REF][START_REF] Gentle | Equilibrium and perturbed fluxes and turbulence levels in a tokamak: Implications for models[END_REF]) that each flux can be expressed as a linear combination of the thermodynamic forces [38, p. 12] [131, p. 340],

Qi = j M ij X j (I.67)
where Qi can be the particle, heat or current density fluxes and where M ij are the transport coefficients, which belong to the transport matrix M [132] [133, p. 12].

The first diagonal terms are the usual diffusion coefficient D linking the particle flux to the density gradient, the heat coefficient χ linking the heat flux to the temperature gradient and the electrical conductivity σ linking the electric current density to the electric field, while some of the off-diagonal terms can be linked to the ware pinch [START_REF] Ware | Pinch effect for trapped particles in a tokamak[END_REF], the thermodiffusive pinch, the curvature pinch or the bootstrap-current.

Perturbative experiments [START_REF] Cardozo | Perturbative transport studies in fusion plasmas[END_REF][START_REF] Gentle | Strong nonlocal effects in a tokamak perturbative transport experiment[END_REF][START_REF] Mantica | Perturbative studies of turbulent transport in fusion plasmas[END_REF][START_REF] Ryter | Perturbative studies of transport phenomena in fusion devices[END_REF] consist in applying a small perturbation to a steady state plasma in order to obtain transport coefficients. They allow to investigate the relation between the fluxes and the thermodynamic forces. For example it has been demonstrated that electron temperature gradient ∇T e has an effect on the particle flux Γ s [START_REF] Terry | Anomalous particle pinch for collisionless plasma[END_REF] through an off-diagonal term of the transport matrix, while the density gradient ∇n s has little effect on the heat flux Q s [132, p. 848]. However, most of the off-diagonal terms have not been measured experimentally.

Let us focus on the particle flux (or density flux). It can be written in terms of transport coefficients :

Γ s = -D ∇n s + M 12 ∇T e + M 13 ∇T i + ... (I.68)
or equivalently :

Γ s = -D ∇n s + n s V s (I.69)
where -D ∇n s is the diffusion term which comes from the diagonal term and n s V s is a convection term accounting for all the off-diagonal terms [132, p. 816].

The transport of particles is governed by the continuity equation :

∂n s ∂t + ∇. Γ s = S (I.70)
where S is a source or sink of particles, which leads in the simple case of a pure diagonal transport matrix, constant D and S = 0 to :

∂n s ∂t = D∆n s (I.71)
which is the usual diffusion equation.

Transport research in the second half of the XX th century relied heavily on experiments and on theoretical calculation. But in the last few decades, the massive progress in computational capabilities has opened new horizons for transport studies. It is now possible via numerical code to simulate some behaviors of a tokamak plasma, although each code still has to focus on particular physical aspects and cannot simulate the whole physics in a tokamak.

The TERESA code introduced in the next section II.1 is a numerical code which focuses on the physics of the trapped particles, averaging out the cyclotron and the bounce motion and keeping the slower toroidal precession motion. In chapter III we use the TERESA code with a newly implemented test particle module in order to separate the diffusive and non-diffusive contributions to the particle flux (see Eq. I.69).

Chapter II

The bounce-averaged gyrokinetic code TERESA

In section I.3, we mentioned that a plasma can be described using a statistical approach with a distribution function f s ( r, v, t) and that its time-evolution could be determined by the Vlasov-Maxwell or by the Vlasov-Poisson system. This description is called a kinetic description of the plasma. Although the Vlasov-Maxwell system is a 6D problem, it is still intractable either analytically or numerically. Fortunately, the field of High Performance Computing (HPC) combined with advances in gyrokinetic theory [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF] recently opened new horizons in the domain of numerical simulations of plasmas [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]. In essence, the gyrokinetic model allows the reduction of the problem from 6D to 5D, by averaging out the fast cyclotron motion of the charged particles. This reduction is relevant in the case of turbulent plasma studies, considering that turbulence typical frequencies are much lower than the cyclotron frequency. This approach has lead to a class of codes called gyrokinetic codes. Standard state-of-the-art gyrokinetic codes non-exhaustively include GYSELA [START_REF] Grandgirard | A drift-kinetic semi-lagrangian 4d code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Global full-f gyrokinetic simulations of plasma turbulence[END_REF][START_REF] Darmet | Intermittency in flux driven kinetic simulations of trapped ion turbulence[END_REF][START_REF] Dif-Pradalier | Interplay between gyrokinetic turbulence, flows, and collisions: Perspectives on transport and poloidal rotation[END_REF][START_REF] Sarazin | Large scale dynamics in flux driven gyrokinetic turbulence[END_REF], GENE [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF][START_REF] Dorland | Electron temperature gradient turbulence[END_REF][START_REF] Goerler | The global version of the gyrokinetic turbulence code gene[END_REF], GKW [START_REF] Peeters | The nonlinear gyro-kinetic flux tube code gkw[END_REF], ELMFIRE [START_REF] Heikkinen | Gyrokinetic simulation of particle and heat transport in the presence of wide orbits and strong profile variations in the edge plasma[END_REF], ORB5 [START_REF] Jolliet | A global collisionless pic code in magnetic coordinates[END_REF], GT5D [START_REF] Idomura | Conservative global gyrokinetic toroidal full-f five-dimensional vlasov simulation[END_REF][START_REF] Idomura | Study of ion turbulent transport and profile formations using global gyrokinetic full-f vlasov simulation[END_REF] and GYRO [START_REF] Candy | An eulerian gyrokinetic-maxwell solver[END_REF][START_REF] Candy | Anomalous transport scaling in the diii-d tokamak matched by supercomputer simulation[END_REF]. These codes require heavy computational capabilities and one well-resolved gyrokinetic simulation with many coupled physical processes in a specific case can take months even on the best supercomputers to this day. These codes aim to have predictive capabilities and give quantitative estimations of macroscopic quantities in tokamaks. Nevertheless this promising approach is limited by the cost (time and funds) of these simulations.

A complementary approach would be to have a reduced code focusing on qualitative results and would be an order of magnitude faster in order to clear up the path for heavier predictive codes. This is the approach taken by the TERESA [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF][START_REF] Ghizzo | Streamer-induced transport in the presence of trapped ion modes in tokamak plasmas[END_REF][START_REF] Cartier-Michaud | Optimizing the parallel scheme of the poisson solver for the reduced kinetic code TERESA[END_REF][START_REF] Cartier-Michaud | Staircase temperature profiles and plasma transport self-organisation in a minimum kinetic model of turbulence based on the trapped ion mode instability[END_REF][START_REF] Drouot | Étude de la turbulence liée aux particules piégées dans les plasmas de fusion[END_REF][START_REF] Médina | Radial density and heat fluxes description in the velocity space: Nonlinear simulations and quasi-linear calculations[END_REF] code which is a global, full f , collisionless code which solves the Vlasov-Poisson system in a 2D phase space + 2 parameters using a Semi-Lagrangian method [START_REF] Sonnendrücker | The semilagrangian method for the numerical resolution of the vlasov equation[END_REF] and making efficient use of parallel computing. The code treats kinetically the motion of trapped particles and focuses on slow (compared to the bounce motion of the trapped particles) phenomena.

In section II.1 we detail the pre-existing model of the TERESA code. Then in section II.2, we present the newly developed test particle module of TERESA.

II.1 Model

II.1.1 The Vlasov -Quasi-neutrality system

In I.3.1 we introduced the Vlasov-Poisson system which is an approximation of the Vlasov-Maxwell system in the case of low-β plasma. We recall that the Vlasov-Poisson system consists of the Vlasov equation :

∂f s ∂t ( α, J) -[H s ( α, J), f s ( α, J)] ( α, J) = 0 (II.1) where [F, G] ( α, J) = 3 i=1 ∂F ∂α i ∂G ∂J i -∂G ∂α i
∂F ∂J i are the Poisson brackets, H s is the Hamiltonian of the system, which we can decompose into an equilibrium and a perturbation parts : As we will see later, in our model we will assume ∆φ ≈ 0 and we will thus use the quasi-neutrality equation.

H s ( α, J) = H eq,s ( J) + Hs ( α, J) (II.
We also recalled in Table I.1 that the motion of a trapped particle can be decomposed in 3 motions of different time and space scales : the cyclotron, the bounce and the toroidal drift motions. This system of equations is 6D and has no analytical solutions and is impossible to solve numerically for non trivial cases even on state of the art super-computers. We will thus have to reduce the dimensionality of the problem.

In the context of low frequency phenomena (compared to the cyclotron and the bounce frequencies), it is possible to average out the cyclotron and bounce motions, which is equivalent to average the system of equation over the angle α 1 and α 2 , using a gyrobounce-gyroaverage operator J 0 (E) on all the quantities which depend on α 1 and α 2 .

II.1.1.1 The gyrobounce-gyroaverage operator J 0,s

Here we introduce the gyrobounce-gyroaverage operator J 0,s which is the development of a product of two Bessel functions in the case where nρ c 1 and kδ b 1 (where n and k are respectively the toroidal and the radial mode numbers), and has the following expression [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF]

J 0,s (E) ≈ 1 - E T eq,s δ 2 b0,s 4 ∂ 2 ∂ψ 2 -1 1 - E T eq,s q 0 ρ 2 c0,s 4a 2 ψ ∂ 2 ∂α 2 3 -1 (II.4)
where

E = 1 2 m s v 2 g + µ s B( r g ) (II.5)
is the kinetic energy of a particle and is an exact constant of motion in the case of a static magnetic field, T eq,s is the equilibrium temperature of species s and q 0 is the safety factor at a radius r 0 which is at the center of the banana. a ψ , ρ c0,s and δ b0,s are respectively the radial size of the system, the thermal cyclotron radius and thermal banana width, all three in units of ψ, so that J 0,s (E) is a dimensionless quantity. It is a differential operator which has to be applied to any function which depends on α 1 and α 2 .

II.1.1.2 The gyrobounce-gyroaverage Vlasov equation

After applying the J 0,s (E) operator on Eq. II.1, the Vlasov equation thus becomes:

∂ fs ∂t (α 3 , J) -[ Hs (α 3 , J), fs (α 3 , J)] ( α, J) = 0 (II.6)
where the double bar symbolizes the double average over the cyclotron and the bounce motions. fs thus represents the "banana-centers" distribution function of the species s. Now let us find an expression for the double averaged equilibrium Hamiltonian Heq (J 3 ).

Considering Eq. I.36, we have

dα 3 dt = ∂H eq ∂J 3 ( J) = 1 Z s e ∂H eq ∂ψ ( J) = ω ds , (II.7) thus, H eq,s ( J) = Z s eω ds ψ + cste (II.8)
when ω ds does not depend on ψ. It is possible to link [65] [70, p. 32] [73] the toroidal drift frequency drift ω ds to the kinetic energy E of a particle, thus:

ω ds = E Z s Ω d (II.9)
where

Ω d = q 0 r 0 B min eR 0 ωd = 1 eR 0 dψ dr ωd (II.10)
can be seen as a toroidal drift frequency which does not depend on the species s,

with ωd = 2E(κ 2 ) K(κ 2 ) -1 + 4s 0 E(κ 2 ) K(κ 2 ) + κ 2 -1 (II.11)
κ is defined as

κ 2 = sin 2 θ crit 2
where θ crit is the critical poloidal angle discussed in I.2.2.2, thus for trapped particle 0 ≤ κ ≤ 1. κ is linked to the first adiabatic invariant µ as κ = ξ 1-λ 2λ where λ = µB min (ψ) E and B min (ψ) is the minimum magnetic field intensity encountered along a particle's trajectory. s 0 = r q dq dr r=r 0 is the magnetic shear at the center of a banana. K(κ) and E(κ) are respectively the complete elliptic functions of the first and second kinds defined as

K(κ) = 1 0 dt √ 1-t 2 √ 1-κ 2 t 2 and E(κ) = 1 0 √ 1-k 2 t 2 √ 1-t 2 dt. The quantity Ω d is hence linked to κ.
The equilibrium Hamiltonian can thus be rewritten

H eq,s ( J) = E(1 + Ω d eψ) (II.12)
where the constant has been chosen equal to E. At equilibrium, the electric field is taken to be zero, but the electric field fluctuations can still occur and the electrostatic potential energy must be added to the Hamiltonian expression:

H s ( α, J) = E(1 + Ω d eψ) + Z s eφ( α, J) (II.13) thus, J 0,s (E)H s ( α, J) = Hs (α 3 , J) = E(1 + Ω d eψ) + Z s e φ(α 3 , J) (II.14)
where φ(α 3 , J) can be seen as the electric potential "felt" by a whole "banana" of charged particle. Vlasov equation thus becomes:

∂ fs ∂t (α 3 , J) -[ Hs (α 3 , J), fs (α 3 , J)] ( α, J) = 0 (II.15)
and can be rewritten:

∂ fs ∂t (α 3 , J) -[ φ(α 3 , J), fs (α 3 , J)] (α 3 ,ψ) + Ω d E Z s ∂ fs ∂α 3 (α 3 , J) = 0 (II.16)
The model can be even more reduced by taking into account that J 1 and J 2 are adiabatic invariants which does not appear in differential operators. J 1 and J 2 can be linked respectively to the trapping parameter κ and E. Thus, it can be written that fs = fs (α 3 , J 3 , κ, E) = fs (α 3 , ψ, κ, E) where κ and E are simple parameters.

II.1.1.3 The quasi-neutrality equation

As we saw in the section I.

λ 2 D ∆ eφ T eq,e = ñe n eq - 1 n eq N j=1 Z j ñj (II.17)
The fluctuations can be decomposed as :

ñs = ñP s + ñT s (II.18)
where ñP s and ñT s are respectively the passing and trapped particle density fluctuations. In the TERESA model, the passing particle density fluctuations are considered adiabatic which means that ñP s has a linear dependence with φ such as: Nevertheless, in Fourier space, the LHS can be written

ñP s = - e s n eq
-   λ 2 D k 2 + 1 + T eq,e n eq N j=1 Z 2 j n eq,j T eq,j   eφ k T eq,e (II.22)
Then if we consider an electric potential fluctuation φ of wavenumber k, we can assume that in a tokamak plasma λ D k 2 1, because generally, electric potential fluctuation length scales (2π/k) are greater than the ion Larmor radius ρ c,i , and typically λ D is smaller of at least an order of magnitude than ρ c0,i , or λ D ρ c,i < 2π/k. Thus we can neglect the term λ 2 D ∆ eφ Teq,e in Eq. (II.17) and we obtain the following quasi-neutrality equation :

ñe n eq = 1 n eq N j=1 Z j ñj (II.23)
Nevertheless, in the TERESA model, only the "banana-centers" distribution function fs appears in the system of equation. Thus, a link needs to be found between fs and ñs , in action angle variables.

Usually the densities are obtained by integrating over the whole velocity space v.

In our gyro-bounce average model we need to integrate over the ( α, J) space.

Assuming that H eq ≈ E [70, p. 177], an infinitesimal volume in velocity space can be written using the magnetic moment µ and E such as

d 3 v = 4π √ 2m -3 2 s √ EdE dµB min (ψ) 4ω b E . (II.24)
Furthermore we assume that the equilibrium distribution function is of the form [START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF]:

f eq,s = n eq,s T 3/2 eq,s e -E Teq,s (II.25)
The calculation of the relation between ñs and fs has been made in [70, p. 42-47] and is written :

ñs n eq,s = e s f T T eq,s q 2 0 a 2 ψ ρ c0,s ∂ 2 α 3 + δ b0,s ∂ 2 ψ φ - 1 -f T f T (φ -φ,s φ α 3 ) -f T + 4π √ 2m -3 2 s f T n eq,s ∞ 0 J 0,s (E) fs √ EdE (II.26)
where

-Z s e(1 -f T ) T eq (φ -φ,s φ α 3 ) = f P s f eq (II.27)
and where f P s is the perturbation in the passing particle distribution function of species s such that f P s = f eq + f P s and where f T is the fraction of trapped particles and has no dimension. Here only the integral over √ EdE is present. The integral over dµB min (ψ) 4ω b E is equal to 1 for the trapped particles.

φ,s is a simulation parameter [70, p. 46] which depends on the species and type of perturbation in the simulation. When φ,s = 0, the passing particles respond adiabatically to every toroidal modes. However when φ,s = 1, the passing particles do not respond at all to the 0 th toroidal mode φ α 3 , but respond adiabatically to every other toroidal modes. Finally, when the plasma contains only one ion species of charge +e and electrons, the quasi-neutrality of our model is written: which can also be written 4π √ 2m

e i T eq,i q 2 0 a 2 ψ ρ c0,i ∂ 2 α 3 + δ b0,i ∂ 2 ψ φ - 1 -f T f T (φ -φ,i φ α 3 ) + 4π √ 2m -3 2 i n eq,i ∞ 0 J 0,i (E) fi √ EdE = - e T eq,e q 2 0 a 2 ψ ρ c0,e ∂ 2 α 3 + δ b0,e ∂ 2 ψ φ - 1 -f T f T (φ -φ,e φ α 3 ) + 4π √ 2m
-3 2 i n eq,i ∞ 0 J 0,i (E) fi √ EdE - m e m i -3 2 ∞ 0 J 0,e (E) fe √ EdE = 1 T eq,i e ∆i φ + τ ∆e φ + C ad (φ -φ φ α 3 ) (II.29) where ∆s = q 0 a ψ ρ c0,s 2 ∂ 2 α 3 + δ 2 b0,s ∂ 2 ψ (II.30)
is a non-isotropic Laplacian operator and where the term 1 T eq,i e ∆i φ + τ ∆e φ represents a polarization density accounting for the difference between the "real" trapped particle density and the gyro-bounce averaged density [START_REF] Ghizzo | Streamer-induced transport in the presence of trapped ion modes in tokamak plasmas[END_REF],

τ = T i T e , (II.31) C ad = e 1 -f T f T (1 + τ ) (II.32) and φ = φ,i + τ φ,e 1 + τ (II.33)
is the simulation parameter discussed previously which, practically, can reduce or force the apparition of zonal flows in the simulation.

In order to numerically solve the system of equation constituted by the Vlasov and the quasi-neutrality equations, one must first nondimensionalize the equations.

II.1.2 Normalization

In this section we present the nondimensionalized equations of the system obtained by normalizing each variable of the system to a characteristic quantity [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF] [70, p. 50] [START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF].

Each nondimensionalized variable will be written with a hatˆ. The poloidal flux ψ used as a radial coordinate can be expressed in base S.I units as

[ψ] = kg.m 2 .A -1 .s -2 . ψ is normalized to L ψ = a dψ dr = aR 0 B θ
which is the radial length of the simulation box in units of ψ so that :

ψ = ψ L ψ (II.34)
α and κ are already dimensionless quantities so ˆ α = α and κ = κ. The energy E is normalized to a characteristic "temperature" T 0 in eV units, so that Ê = E T 0 (II.35)

The time t is normalized using a characteristic angular frequency ω 0 = T 0 eR 2 0 B θ which represents the toroidal precession angular frequency of the trapped particles at a reference temperature T 0 such that Ωd = Ω d T 0 ω 0 . Therefore:

t = ω 0 t (II.36)
Thus, for the simulations discussed in this thesis, t = 1 is equivalent to approximately t = 1 ms in real time for a tokamak. The electric potential's base S.I units are

[φ] = kg.m 2 .A -1 .s -3 and is normalized as φ = φ ω 0 L ψ = R 0 a eφ T 0 (II.37)
The gyro-bounce averaged distribution function is normalized using a typical density n 0,s of a species s in a 3D real space and a typical velocity space volume T 0,s ms 3 2 , so that:

fs = 1 n 0,s 2π T 0,s m s 3 2 fs (II.38)
The 2π factor appears so that the integral over the whole phase space of the equilibrium distribution function feq,s is equal to 1.

The following equations are dimensionless but now the hat ˆis omitted for readability purposes. The nondimensionalized Vlasov equation for a species s is thus written:

∂ fs ∂t (α 3 , J) -[J 0,s (E)φ( α, J), fs (α 3 , J)] (α 3 ,ψ) + Ω d E Z s ∂ fs ∂α 3 (α 3 , J) = 0 (II.39)
In the case of a plasma containing only one ion species of charge +e and electrons of charge -e, we thus write two nondimensionalized Vlasov equations:

∂ fi ∂t (α 3 , J) -[J 0,i (E)φ( α, J), fi (α 3 , J)] (α 3 ,ψ) + Ω d E ∂ fi ∂α 3 (α 3 , J) = 0 (II.40) ∂ fe ∂t (α 3 , J) -[J 0,e (E)φ( α, J), fe (α 3 , J)] (α 3 ,ψ) -Ω d E ∂ fe ∂α 3 (α 3 , J) = 0 (II.41)
which are closed by the quasi-neutrality equation treating kinetically the trapped ions and electrons while treating adiabatically the passing particles, written in nondimensionalized form as:

2 √ πn eq ∞ 0 J 0,i (E) fi (α 3 , J) √ EdE - ∞ 0 J 0,e (E) fe (α 3 , J) √ EdE = 1 T eq,i C ad (φ( α, J) -φ φ( α, J) α 3 ) -C pol ∆i+e φ( α, J) (II.42) with C pol = eω 0 L ψ T 0 = a R 0 (II.43) C ad = C pol 1 -f T f T (1 + τ ) (II.44) φ = φ,i + τ φ,e 1 + τ (II.45) and ∆i+e = ∆i + τ ∆e =   q 0 a ψ ρ c0,i 2 ∂ 2 α 3 + δ 2 b0,i ∂ 2 ψ   + T i T e   q 0 a ψ ρ c0,e 2 ∂ 2 α 3 + δ 2 b0,e ∂ 2 ψ   (II.46)
It is also possible to use this model to take into account multiple ions with different charges, which can be used in order to study impurities transport [START_REF] Idouakass | Impurity density gradient influence on trapped particle modes[END_REF][START_REF] Gravier | Diffusive impurity transport driven by trapped particle turbulence in tokamak plasmas[END_REF] [70, p. 52], but will not be considered in this work.

Electron dissipation modeling When considering trapped ion resonance-driven turbulence, it is possible to model the effects of electron-ion collisions (collisional electron dissipation) without solving the electron distribution function, nor adding a collision operator [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF]. It is done introducing an iδ m term in the quasi-neutrality, which accounts for a phase-shift between electron density and perturbed electric potential, which writes in Fourier space: ne = (1 + iδ m ) φ, where δ m = mδ with m the m th component of the Fourier decomposition of φ in the α direction and δ is a constant input parameter. Thus, the trapped electrons respond non-adiabatically to φ, although they are not considered "kinetic". When studying TIM turbulence we can thus decrease numerical computation cost and gain analytical tractability, in the philosophy of the reduced model of TERESA. The quasi-neutrality constraint, Eq. (II.42), can thus be re-written in the case of TIM focused turbulence as:

C ad φ + F -1 iδ m φm -C pol ∆i φ = 2 √ π ∞ 0 J 0,i (E) fi (α 3 , J) √ EdE -1 (II.47)
with i the imaginary number and F -1 the inverse Fourier transform. In this case, passing ions and electrons are still adiabatic, trapped ions are kinetic, and trapped electrons are a non-adiabatic neutralizing background.

This iδ m model is one of the three modeling options of the TERESA code:

• Kinetic trapped ions with no electron dissipation (no collision modeling)

• Kinetic trapped ions with electron dissipation, iδ m (collision modeling)

• Kinetic trapped ions and Kinetic trapped electrons (no collision modeling)

Now that we obtained a dimensionless closed set of equations, let us see how they are solved numerically.

II.1.3 Solving numerically the model equations

In this section we will briefly see how the system of equations is solved numerically. Solving the Vlasov equations and the quasi-neutrality constraint means finding fs (α 3 , J) and φ( α, J) at any time t. The numerical scheme accuracy is of second order in time (Appendix A).

II.1.3.1 Solving the Vlasov equations using the Semi-Lagrangian method

Each Vlasov equation will be solved using the Semi-Lagrangian method [START_REF] Ghizzo | Streamer-induced transport in the presence of trapped ion modes in tokamak plasmas[END_REF][START_REF] Sonnendrücker | The semilagrangian method for the numerical resolution of the vlasov equation[END_REF] which we will very briefly explain here.

As stated in section II.1.1.2, J 1 and J 2 which are respectively linked to E and κ, as well as α 1 and α 2 does not appear inside differential operators in our model. It is the case only for ψ (linked to J 3 ) and α 3 which we will now simply call α. α and ψ are thus the two variables of the model while E and κ are two parameters.

For the simulations discussed in this thesis, we impose Ω d = 1, which forces the trapped particles to be deeply trapped and is equivalent to using only a single value of κ = 0. Therefore it is possible to take

N E values of E such as E = {0 ≤ E ≤ E max } in order to separate fs (α, ψ, κ = 0, E) into N E distribution functions such as fs,E (α, ψ)
, where E is a fixed value.

We thus have to solve N E Vlasov equations for each species s, plus one quasineutrality equation.

As stated in section I.3.1, each distribution function fs,E (α, ψ) is conserved along the particle trajectories, or characteristics. These characteristics in (α, ψ) space are written as follows:

α = dα dt = ∂ Hs ∂J 3 = E Ω d Z s + ∂ φ ∂ψ (α, ψ, E , t) (II.48) and ψ = dψ dt = - 1 Z s e ∂ Hs ∂α = - ∂ φ ∂α (α, ψ, E , t) (II.49)
Consequently, the Vlasov equation can be written in its advective form [START_REF] Sonnendrücker | The semilagrangian method for the numerical resolution of the vlasov equation[END_REF]

∂ fs,E ∂t ( X(t), t) + U ( X(t), t). ∇ X(t) fs,E ( X(t), t) = 0 (II.50)
where X(t) = (α, ψ) represents the phase space coordinates. U is a vector field also called the advection field,

U = α ψ . (II.51) Therefore we can write d X dt = U ( X(t), t) (II.52)
The 2D phase space is discretized using respectively N α and N ψ points in the α and ψ directions. At a time t = j∆t, where j is the iteration number and ∆t the time step, f j s,E ( X, j∆t) is known for every points on the phase space grid. Thus we want to know the value of the distribution function at the next iteration f j+1 s,E ( X, (j + 1)∆t) for every phase space grid points. We call x j+1 the position of one point on the phase space grid at time t = t + ∆t where we want to know the value of f j+1 s,E ( x j+1 ). The distribution function is conserved along the characteristics, thus it means that somewhere in the phase space at time t and at a position x j which we want to know, we can find the same value for the distribution function such that f j+1 s,E ( x j+1 ) = f j s,E ( x j ). The difference of position between these two points can be written:

x j+1 = x j + δx (II.53)
where δx is obtained by integrating the advection field from t to t :

δx = t+∆t t U ( X x j+1 ,t x j ,t (s), s)ds (II.54)
and where X

x j+1 ,t x j ,t (s) (II. [START_REF] Jenko | Electron temperature gradient driven turbulence[END_REF] is the one characteristic which connects the points x j+1 at time t = t + ∆t and x j at time t. Note that while the point x j+1 is taken on a phase space grid point, the point x j can be anywhere in the phase space and not necessarily on the grid, see Fig.

II.1.
Therefore, f j s,E ( x j ) must be found by interpolation of the nearest 4 grid points using the bicubic spline method [START_REF] Boor | A practical guide to splines[END_REF] and finally the value of the distribution function at time t is equal to this interpolated value. This procedure now has to be repeated for every grid points in order to know f j+1 s,E on the whole grid at time t .

II.1.3.2 Solving the quasi-neutrality equation

To find the electric potential φ, we now need to solve the quasi-neutrality equation where

δn(α, ψ) = 2 √ πn eq ∞ 0 J 0,i (E) fi √ EdE - ∞ 0 J 0,e (E) fe √ EdE (II.58)
Thus the quasi-neutrality equation can be written for any m th mode of the α direction:

1 T eq,i C ad (1 -φ δ m,0 ) + C pol q 2 0 a 2 ψ (ρ 2 c0,i + τ ρ 2 c0,e )m 2 -C pol (δ 2 b0,i + τ δ 2 b0,e )∂ 2 ψ φ m (ψ) = δn m (ψ) (II.59) with C pol = eω 0 L ψ T 0 , C ad = C pol 1-f T f T (1 + τ ), φ = φ,i +τ φ,e 1+τ
and where δ m,0 is the Kronecker delta which is equal to 1 for the 0 th mode and equal to 0 for every other m modes. Hence, with the finite difference method in ψ, the number of equations to solve to obtain φ(α, ψ) is

N ψ × N α -1 2 + 1 (II.60)

II.1.4 Boundary and initial conditions

In the α direction we use 2π periodic boundary conditions for every quantity that have an α dependence. In the ψ direction we use Dirichlet boundary conditions so that φ = 0 at ψ = 0 and ψ = L ψ . Fixing the value of the electric potential at the edges of the box usually causes numerical instabilities which can distort the simulation results. In order to avoid this problem, it is possible to add buffer zones to both edges of the box in ψ which can soften the decrease in φ toward the edges.

In practice, we add a diffusion term to the right hand side of the Vlasov equation so that:

∂ fs ∂t -[J 0,s (E)φ( α, J), fs ] (α 3 ,ψ) + Ω d E Z s ∂ fs ∂α 3 = ∂ ψ (D ψ ∂ ψ fs ). (II.61)
This artificial diffusion is governed by a diffusion coefficient D ψ defined as:

D ψ (ψ) = D 0 2 -tanh ψ -ψ 0 L D + tanh L ψ -ψ 0 -ψ L D (II.62)
where D 0 , ψ 0 and L D are input parameters.

Motivations

We described the TERESA code in section II.1. Although TERESA solves the distribution function f and the electric potential φ, it does not yield individual particle trajectories. Particles, momentum and energy fluxes can be obtained from f and φ but discriminating diffusive and convective processes typically requires convoluted methods such as dedicated dynamical synthetic experiments. Investigating the particle trajectories in the turbulent plasma would lead to have better insights on diverse phenomena occurring in tokamaks such as diffusion [START_REF] Rechester | Electron heat transport in a tokamak with destroyed magnetic surfaces[END_REF][START_REF] Isichenko | Invariant measure and turbulent pinch in tokamaks[END_REF][START_REF] Harvey | Radial transport and electron-cyclotron-current drive in the tcv and diii-d tokamaks[END_REF][START_REF] Estève | Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport[END_REF], hyper-or sub-diffusion [START_REF] Misguich | Noble internal transport barriers and radial subdiffusion of toroidal magnetic lines[END_REF], advection [START_REF] Isichenko | Invariant measure and turbulent pinch in tokamaks[END_REF][START_REF] Boedo | Flow reversal, convection, and modeling in the diii-d divertor[END_REF][START_REF] Boedo | Transport by intermittent convection in the boundary of the diii-d tokamak[END_REF][START_REF] Dif-Pradalier | On the validity of the local diffusive paradigm in turbulent plasma transport[END_REF], the trapping of particles in potential wells, and ballistic events such as avalanches [START_REF] Sarazin | Large scale dynamics in flux driven gyrokinetic turbulence[END_REF][START_REF] Idomura | Study of ion turbulent transport and profile formations using global gyrokinetic full-f vlasov simulation[END_REF][START_REF] Dif-Pradalier | On the validity of the local diffusive paradigm in turbulent plasma transport[END_REF][START_REF] Sarazin | Predictions on heat transport and plasma rotation from global gyrokinetic simulations[END_REF]. Indeed, the analysis can be done locally in space, within a short timespan, and without ambiguity.

In order to have access to the particle trajectories, we added test particles to TERESA during this thesis. Test particles are particles advected by the electrostatic field, but which do not affect it. They can thus be used as markers in the turbulent plasma, representing exactly the motion of a single particle belonging to f . The test particle trajectories are computed directly in the TERESA code thus allowing the same order of accuracy as the solving of f and φ.

There are multiple approaches to use test particles and the main ones are: 1determine an electric potential map from analytical methods and let the test particles evolve in it [START_REF] Kwon | Global drift wave map test particle simulations[END_REF], or, 2-obtain the electric potential map either from a kinetic simulation (or experimental measurements [START_REF] Wünderlich | On the proton flux toward the plasma grid in a rf-driven negative hydrogen ion source for iter nbi[END_REF][START_REF] Minea | Kinetics of plasma species and their ionization in short-hipims by particle modeling[END_REF]) and determine the test particle trajectories in post processing, or, 3-solve the test particle trajectories directly in the kinetic simulation. Method 1 offers the main advantage that it does not require large computational power but it relies on a predetermined analytical description of the electric potential configuration. Method 2 is usually more computationally intensive as it generally requires a nonlinear kinetic simulation. One downside is that the test particle trajectories which are solved in post processing are not solved at the same order of precision as the simulation: ∆t dt where ∆t is the kinetic simulation time-step and dt the post processed test particle trajectories time-step. Although it is in principle possible to use method 2 with ∆t = dt and obtain method 3, we distinguish method 2 from method 3 because it would require the saving of the potential map at each dt and would thus be prohibitively expensive in term of numerical storage. In this work we will use method 3 which has the advantage of giving the test particle trajectories with the same order of precision as the numerical scheme of the TERESA simulation (∆t = dt), but with the downside of being more expensive in terms of computational time and numerical storage. With one million test particles, the TERESA simulation usually takes twice as much time than without test particles.

II.2.2 Implementation

The test particles, which are de facto "test banana-centers", follow the characteristics of fs which in (α, ψ) space are :

α = E Ω d Z s + ∂ φE ∂ψ (α, ψ, E , t) (II.63) and ψ = - ∂ φE ∂α (α, ψ, E , t) (II.64)
These positions are solved at each time-step using an explicit Runge-Kutta 4 algorithm method and they depend on the energy parameter

E = {0 ≤ E ≤ E max }.
The Runge-Kutta 4 algorithm requires the advection field

U = α ψ (II.65)
at the beginning of the time-step, at half of the time-step and at the end of the time-step.

The solving of the test particle trajectories is parallelized in E in the same way that the Vlasov-Quasi-neutrality solving is parallelized, as explained in section II.1.3.1.

II.2.2.1 Verification with an analytic electric potential

Let us first verify that the RK4 algorithm is properly implemented in the TERESA code. We choose the verification procedure as follows: we start a reference simulation with a constant time-step widths ∆t considered to be "small enough", we choose ∆t ref = 1 × 10 -5 . Then we start the same simulation but with different constant time-step width ∆t > ∆t ref , we take ∆t ∈ [5×10 -5 ; 1×10 -4 ; 5×10 -4 ; 1× 10 -3 ].

For each simulation, at a time 1 × 10 -3 , we calculate the relative error on α and ψ of one test particle, compared to the α ref and ψ ref of the reference simulation. We conduct the test with an analytic potential φ analytic instead of the self consistent field obtained from the Vlasov-Quasi-neutrality system. With an analytic potential, the accuracy of the particle trajectories is not limited by the Semi-Lagragian solver accuracy but by the Runge Kutta 4 algorithm with ∆t 4 accuracy. We choose φ analytic as:

φ analytic (α, ψ) = 2 sin(πψ) 100 m g m (α, ψ) (II.66) with g m (α, ψ) = cos(mα -ω m t)   a 1 (m)e - ψ-ψ 1 (m) ∆ 1 (m) 2 + a 2 (m)e - ψ-ψ 2 (m) ∆ 2 (m) 2   + sin(mα -ω m t)   a 3 (m)e - ψ-ψ 3 (m) ∆ 3 (m) 2 + a 4 (m)e - ψ-ψ 4 (m) ∆ 4 (m) 2   (II.67)
with m the mode number, a 1,2,3,4 (m), ψ 1,2,3,4 (m), ∆ 1,2,3,4 (m) constant parameters and φm (ψ) is Gaussian.

We conduct the test with this analytic electric potential, at t = 1 × 10

II.2.2.2 Verification in a self consistent electric potential

As detailed in Appendix A, the Vlasov-Quasi-neutrality solver in TERESA is of second order in time (∆t 2 ). It means a Runge Kutta 2 method would be acceptable since it has a local error of third order and a global error of second order, but we choose the RK4 method because the increase in computational time from RK2 to RK4 is negligible. Let us verify that test particle trajectories are computed with second order accuracy in time in the self consistent electric potential calculated by TERESA. We conduct the same test as in II.2.2.1 but with a self consistent field obtained from the Vlasov-Quasi-neutrality system instead of an analytic potential φ analytic . We choose the reference time-step width as ∆t ref = 5 × 10 -6 and other simulations are made with ∆t ∈ [1 × 10 -5 ; 5 × 10 -5 ; 1 × 10 -4 ; 5 × 10 -4 ; 1 × 10 -3 ; 5 × 10 -3 ; 1 × 10 -2 ]. II.1. This test confirms that test particle position accuracy is of second order in time as δα and δψ grows like ∆t 2 , which is the same order of accuracy as the Vlasov-Poisson solver for f and φ.

In terms of computational cost of the test particle module, the test conducted on a single test particle with the self consistent field showed a negligible increase in computation time compared to the same simulation without test particle. In the following, we usually use one million test particles, the increase in computation time is around 100% in this case.

II.3 Chapter conclusion

We developed, parallelized, integrated and tested a new test particle module in the bounce-averaged electrostatic gyrokinetic code TERESA. TERESA uses a reduced model, averaging out the cyclotron motion, as well as the trapped particle bounce motion, in order to focus on the trapped particle precession in the toroidal direction (the toroidal drift). Before this work, TERESA only solved the trapped particle distribution functions and the electric potential from the Vlasov-Quasineutrality system of equations, using a Semi-Lagrangian solver.

The new test particle module can now compute the test particle trajectories of millions of test particles. The test particle follow the Vlasov characteristics in phase space and the trajectories are solved using RK4 algorithm. We first verified the accuracy of the trajectories solving, using an analytic electric potential. It showed that the trajectories were computed with fourth order accuracy in time (∆t 4 ) accordingly to the RK4 algorithm. However, when using the self consistent electric potential obtained from the Semi-Lagragian solver, the accuracy was limited to second order in time (∆t 2 ), which is the order of accuracy of the Semi-Lagragian solver. We could use a RK2 algorithm to solve the trajectories but the difference in computation time between RK2 and RK4 is negligible therefore we use RK4.

Knowing the trajectories of the test particles, which are de facto "test banana centers", gives more insights on the transport phenomena such as diffusion, advection, trapping or ballistic motions.

In chapter III we exploit this new module using one million test particles on each energy grid point, in a turbulent simulation.

Chapter III

Transport characterization via test particles

In section II.2, we detailed the implementation of a new test particle module in TERESA and the motivations behind it.

In this chapter we study the evolution of millions of test particles in a turbulent plasma simulation, as a method to get insights on the type of transport governing the plasma. The goal is to compare the diffusive contribution to the particle flux with the non-diffusive contributions. We initialize the test particles, which are de facto "test banana-centers", at a time of the simulation when the plasma is turbulent. We impose an initial temperature and density gradients. We run a first simulation with fully kinetic trapped ions and "quasi-adiabatic" trapped electrons, and a second one with fully kinetic trapped particles.

The first simulation achieves a state of "global turbulence", while the second one is dominated by streamer-like structures. We calculate the Mean Squared Displacement (MSD) of the test particles as a function of time in order to obtain a random walk diffusion coefficient depending on the toroidal precession kinetic energy (E) of the particles. A radial particle diffusion flux is then calculated and compared to the total radial particle flux accounting for all the transport processes such as diffusion and advection which is obtained directly from the TERESA code. We find that for the first simulation, the total radial particle flux is essentially diffusive which is consistent with our simulation set-up aiming for "global turbulence". Both fluxes present a peak around a resonance energy E R ≈ 1.74 between the TIM and the particles. Both thermal and high-energy particles do not contribute significantly to radial transport. For the second simulation, which is streamer-like dominated, the total flux is the most intense for thermal particles (low E), and the diffusive processes at E = 0 accounts for only ∼ 10% of the total flux. We find that the transport of particles is dominated by non-diffusive processes.

59

Chapter III. Transport characterization via test particles

This chapter has been the subject of an article [START_REF] Médina | Test particle dynamics in low-frequency tokamak turbulence[END_REF].

III.1 Test particle transport in trapped ion turbulence III.1.1 Simulation configuration

We treat kinetically the trapped ions and the trapped electrons are assumed as a neutralizing background. Passing ions and electrons respond adiabatically to the electric potential φ. We introduce the term δ m which models the effects of electron-ion collisions as a phase shift between φ and the electron density [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF] The quasi-neutrality constraint, Eq. (II.42), can thus be re-written as:

C ad φ + F -1 iδ m φm -C pol ∆i φ = 2 √ π ∞ 0 J 0,i (E) fi (α 3 , J) √ EdE -1 (III.1)
where F -1 is the inverse Fourier transform, i is the imaginary number and we have taken φ = 0. We use a uniform grid in phase space: N ψ points in ψ ∈ [0; 1] where ψ = 1 is the center of a poloidal section and ψ = 0 is toward the edge (but still fulfilling core plasma conditions) and N α points in the toroidal precession angle α ∈ [0; 2π[. The number of points is N α × N ψ = 2045 × 1025. For the energy E, we choose a non uniform grid spacing with the introduction of a new parameter V = √ E, with N V = 96 points. The range E ∈ [0; 20] is chosen to allow good convergence of simulations results. We recall the grid configuration in table IV.1 and the normalization in table III.2.

Grid Number of grid points Value

α N α = 2049 α ∈ [0; 2π[ ψ N ψ = 1025 ψ ∈ [0; 1] E, V N E or N V = 96 E ∈ [0; 20]
Table III.1: Grid used for our simulation. α and ψ are the phase-space variables while E (or V ) is a parameter. For boundary conditions we use thermal baths on ψ = 0 and ψ = 1 thus we can impose an initial temperature gradient length κ T = 0.15 and an initial density gradient length κ n = 0.05.

The ion Larmor radius is ρ i = 0.001 and the ion banana width is δ bi = 0.01. The initial electrostatic potential is a sum of sines both in α and ψ with random phases and we choose the equilibrium ion distribution function f eq as locally Maxwellian (exponential in E):

f eq (ψ, E) = e -E [1 + ψ(κ T (E -3/2) + κ n )] (III.
Modes magnitude We recall the input parameters in table IV With the goal of studying test particle diffusion in a typical core plasma, we do not want zonal flows or streamers to be dominant because they would drastically enhance or reduce the radial transport of particles. We want a simulation with a "global turbulence" at the time of the study, where global means that there is no large electric potential structure either in ψ or in α. Then the modes reach the saturation level at t ≈ 2 and nonlinear interactions are dominant. This is the turbulent phase where the modes are Trapped Ion Modes (TIM). The 0 th mode is not dominant throughout the simulation.

A "global turbulence" should not be dominated by large scale modes such as k αm L α ∼ 1 and k ψn L ψ ∼ 1, where k αm and k ψn are the m th and n th wavenumbers of a φ wave in (α, ψ) directions and L α = 2π and L ψ = 1 are the sizes of the box in α and ψ. Moreover we should have a bulk of most intense α-modes (not one dominant mode over the others) so that ∆m ∼ m where ∆m is the mode-range of the bulk of most intense α-modes and m the mean mode of this bulk. The autocorrelation time of the α-modes is τ α ≈ 1. The spectrum of α-modes averaged between t = 6 and t = 7 is shown in Fig. III.4a. It is the absolute value of the Fourier transform in the α direction of the electric potential φ. We can see that the mode-range of the bulk of most intense α-modes is approximately ∆m ≈ 10 and the mean mode of this bulk of modes is m ≈ 10 so we do not have one very dominant mode but rather a collection of dominant modes of about the same level of intensity. Therefore we choose to study the test particle diffusion at time t = 6. At this time we have a ratio eφ/T ∼ 0.01 -0.03, Fig. III.5, which is typical in core fusion plasma. In addition, we observe that the electric potential map shows no large structure such as zonal flow or streamer.

We choose to initialize the test particles at time t = 6, with a Gaussian distribution in ψ centered in ψ = 0.5 and with a standard deviation ∆ψ = 0.022, in order to minimize the sensitivity to the radial variations of turbulence. In the α direction the test particles are distributed randomly. They have a fixed E. [START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF], from initialization at t = 6. For E = 0, the TPDF goes from a centered distribution around ψ = 0.5 at t = 6 and t = 6.1 to a skewed TPDF for later times t ≥ 6.4. This skewness in the TPDF indicates an advection of the test particles. This advection pushes the particles toward ψ < 0.5 at t = 6.4 and t = 8, then toward ψ > 0.5 at t = 10. This oscillatory motion continues so that on longer time scales, the TPDF stays in average around ψ = 0.5, with an increasing standard deviation (a flattening of the TPDF), thus, on long time scale (∆t 5) the TPDF undergoes a diffusion process and not an advection process. For E = 1.74, we can make a similar observation, with the difference that the TPDF flattens at a faster rate than at E = 0 and the skewness (advection) is less present. Fig. III.7 shows the trajectories in (α, ψ) space of 3 test particles with different energies, E = 0, E = 0.8 and E = 1.74, as a function of time. These particles are selected arbitrarily but each particle appears to be qualitatively representative of other particles with the same energy. The test particle with an energy E = 1.74 is transported to the edge of the simulation box faster than the other test particles. Section III.1.2 shows that the highest rate of test particle radial transport is observed for E = 1.74. For each E we use 10 6 test particles: 1000 × 1000 in the ψ and α directions. This simulation took about 9 days of computation on 96 cores.

III.1.2 Test particle dynamics in a turbulent plasma simulation

In this section, we aim at developing a robust method for dissociating radial diffusion and radial convection of the test particles. We first study the time-evolution of the test particles Mean Squared Displacement (MSD) in the radial direction ψ for each E as they evolve in a turbulent plasma simulation. It allows to estimate the turbulent radial diffusion coefficient in velocity space. We then find the radial diffusive flux of the test particles and we compare it to the total radial particle flux accounting for all the transport processes.

III.1.2.1 Estimation of a random walk radial diffusion coefficient in velocity space for the test particles

We let the test particles evolve in the turbulent simulation starting from time t 0 = 6. At each time-step of the simulation, TERESA calculate the MSD of the test particles in the radial direction, where the average is over all the test particles : (ψ(t)-ψ(t 0 )) 2 , for each energy E. Later we will see that around E ≈ 1.74 there is a resonance between trapped ions and TIM therefore we plot the time-evolution of this MSD for E = 0, E = 0.8, E = 1.74, E = 2.71 and E = 3.2 respectively on Figs. III.8a, III.8b, III.9a, III.9b and III.10. For each E we can distinguish different phenomena in time for the MSD. When the MSD grows linearly in time, we superpose the slope in red to the plot, and we will be able to calculate a diffusion coefficient.

In the general case we expect the MSD to have 4 different phases at different time and space scales:

1. Phase 1: a first, rapid (∼ 0.1Ω -1 d ) phase of local convection where the test particles reorganize themselves inside the local potential structure where they have been initialized in, typically on a space scale of 10 -2 in units of ψ.

2. Phase 2: a phase of fast-diffusion (∼ 0.2Ω -1 d ), on a space scale of 10 -2 -10 -1 in units of ψ 3. Phase 3: In this phase, the evolution of the MSD is somewhat complex, with strong fluctuations on a ∼ 0.1Ω -1 d timescale, indicating that the particle motion is not a simple combination of diffusion and convection on this timescale. However, on a timescale ∼ Ω -1 d , of the order of the turbulent autocorrelation time, the MSD grows roughly linearly in time, on the space scale of the simulation box. Therefore, transport may be modeled by a simple diffusive process on this longer timescale. We thus make a linear fit on phase 3, which we superpose in red to the plot, to find the diffusive coefficient The diffusion coefficient has a peak (D RW ≈ 5.2 × 10 -2 ) around the resonance energy E R ≈ 1.74 because at this E, the test particles tend to move simultaneously with the electric potential and thus diffuse in the radial direction much faster than at other E. In Chapter IV, we will introduce quasi-linear theory which allow, under certain conditions, to estimate the value of this resonance energy.

For high E, the test particles have a high velocity α compared to the evolution of the electric potential and tend to perceive only an average of φ along their trajectories, thus their radial diffusion coefficient is much smaller.

To confirm that the diffusion coefficient calculated from the MSD is not spuriously influenced by convection, we analyze the standard deviation of the ψ-distribution of test particles, which cannot be influenced by convection. We find that there is no significant difference between the time-evolution of the variance and that of the MSD, except for a constant shift due to a finite initial standard deviation, which has no impact on the slope. Therefore, this second method of analysis, which unambiguously discriminate convection and diffusion, confirms the results of the first method. This agreement indicates that, on a timescale of the order of the turbulence auto-correlation time, transport is predominantly diffusive in this simulation, with this set of parameters and this type of turbulence.

III.1.3 Comparison between the radial diffusion flux of the test particles and the total radial particle flux, in velocity space

From the random walk diffusion coefficient estimated in the previous section, we can estimate a radial diffusive flux for the test particles as which includes all radial transport processes such as diffusion or advection and which we average over ψ ∈ [0.4; 0.6] and over an auto-correlation time t ∈ [6; 7]. We find that radial particle transport is dominated by the resonant particle around the energy E R ≈ 1.74 where the trapped ions resonate with the TIM. As the peaks of the total flux and the diffusive flux are of the same magnitude, we can say that the transport of resonant particle is exclusively diffusive, following a random walk process, and moreover that the whole radial particle transport is dominated by diffusive processes. This is coherent with our choice of "global turbulence" as we chose to favour a turbulence driven by a bulk of dominant TIM, with non-dominant large potential structure.

Γ D RW = -D RW ∂ f α ∂ψ ψ∈[
Both flux peaks are negative (directed toward the edges) which is coherent with the gradient ∂ f α ∂ψ ψ∈[0.4;0.6] being positive at E = E R , recalling that ψ = 1 is toward the core and ψ = 0 is toward the edge.

For high E and thus high α, particle radial transport tends to be negligible and thus the two fluxes tend to 0, as explained in Sec. III.1.2.1.

For E < 1, the gradient ∂ f α ∂ψ ψ∈[0.4;0.6] is negative, thus the diffusive flux is slightly positive (directed toward the core). The total flux is of the same sign than the diffusive flux thus they are in the same direction. Between E = 0 and E ≈ 0.5, the total flux is a little less intense than the diffusive flux, indicating that the total flux may have a significant non-diffusive component directed toward the edge, although the discrepancy may be due to the uncertainty in measuring the slope of the MSD. Between E ≈ 0.5 and E ≈ 1 the total flux is more intense than the diffusive flux indicating that the total flux may have a non-diffusive component in the same direction than the diffusive flux, i.e. directed toward the core. The use of test particles thus allows us to estimate the diffusive part of the flux in the total particles flux.

III.2 Test particle transport in trapped ion and electron turbulence

In this section we use the complete Vlasov-Quasi-neutrality system described by eqs. II.40, II.41 and II.42 with fully kinetic trapped ions and electrons. Aiming to separate the diffusive contribution to the radial particle flux from the non-diffusive contributions, we proceed in the same way as in the previous section. We initialize one million test particles at each energy E grid point in a turbulent simulation. In terms of computational time, the simulation is heavier compared to the one with only kinetic ions and the simulation took 6 days on 256 cores. We give the grid configuration in Table III.4 and the input parameters in Table III.5. [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF] δ m = 0 Table III.5: Input parameters.

Grid Number of grid points Value

α N α = 2049 α ∈ [0; 2π[ ψ N ψ = 1025 ψ ∈ [0; 1] E, V N E or N V = 256 E ∈ [0; 20]
Ω d = 1 φ φ = 0 C ad C ad = 0.1 C pol C pol = 0.1 Electron dissipation

Fig. III

. [START_REF] Knutti | A scientific critique of the two-degree climate change target[END_REF] shows the time evolution of 5 modes in the α-direction along with the 0 th mode. The mode magnitudes grow exponentially from t = 0 to t ≈ 3 and then reach a saturation level. Fig. III.14 and III.15a show respectively the α-mode spectrum averaged over t ∈ [6; 7] and at ψ = 0.5, and the electric potential map at t = 6. We observe that the situation is different from the previous case. The turbulence is dominated by streamer-like, ψ-elongated α-modes and is not "global". Fig. III.15b shows the time evolution of φ as a function of α and t, for ψ = 0.5. When an α-mode propagates toward decreasing α, its real frequency is negative which indicates a trapped electron resonance-driven mode (TEM). We observe such mode for t 7. On the contrary, for 3 t 7, the modes propagates toward increasing α, which indicates that the turbulence is dominated by TIM. Fig. III. [START_REF] Fischer | Anthropogenic contribution to global occurrence of heavy-precipitation and high-temperature extremes[END_REF] shows the trajectories in (α, ψ) space of 2 test particles with different energies, E = 0 and E = 1.38, as a function of time. These particles are selected arbitrarily but each particle appears to be qualitatively representative of other particles with the same energy. The trajectories appear to present a non diffusive behavior. Using the same procedure as in the previous section, we study the MSD time evolution of one million test particles for each E. Function (TPDF) as a function of ψ and for different times. We observe that from t = 6 to t = 6.06 the test particles explored the potential structure they were initialized in, which in streamer-like turbulence, leads to a rapid flattening of the TPDF. From t = 6.2 to the end of the simulation, the test particles continue to be transported and the TPDF continues to be flattened, until saturation where the TPDF becomes totally flat and the test particles effectively explored the entirety of the simulation box in ψ. We choose to estimate the slope of the MSD time evolution between t = 6.5 and t = 10 in order to calculate a diffusion coefficient.

It is worth to note that the MSD points are plotted every ∆t = 0.02 and that with a lower ∆t, the MSD from t = 6 to t ≈ 6.04 should not grow linearly in time, but as t 2 , characteristic of an advection process.

For which is an order of magnitude lower than the first linear slope, plotted in red. Fig. III.18b shows the TPDF as a function of ψ and for different times. We can see that the test particles explore the potential structure they were initialized in from t = 6 to t ≈ 6.1. The test particles then continue to be transported and the TPDF continues to flatten. Between t = 7.5 and t = 10, the TPDF stays similar. We choose to estimate the slope of the MSD time evolution between t = 6.1 and t = 7 in order to calculate a diffusion coefficient.

For E = 2.02, Fig. III.19a, the test particle MSD again presents a first rapid growth (from t = 6 to t ≈ 6.1) to a MSD∼ 0.03. The MSD then grows linearly in time until the end of the simulation. Fig. III.19b shows the TPDF as a function of ψ and for different times. From t = 6 to t ≈ 6.1, the test particles explore the potential structure they were initialized in and the TPDF flattens rapidly. We choose to estimate the slope of the MSD time evolution between t = 6.1 and t = 10 in order to calculate a diffusion coefficient. Figure III.17: (III.17a) Time-evolution of the test particle MSD for E = 0, with the estimated slope for the interval t ∈ [6.5 -10] in red, and (III.17b) TPDF as a function of ψ, for E = 0, at times t = 6, t = 6.02, t = 6.04, t = 6.06, t = 6.2 and t = 10. nature of the turbulence. We find that both fluxes are fundamentally different from the previous simulation (III.12) and that the fluxes present different structures. The peak of the diffusion coefficient appears in the diffusive flux as a negative flux around E ≈ 1.8 and the diffusive flux is the most intense around E ≈ 0.1, Γ D RW ≈ 3 × 10 -4 . Nevertheless, the total flux for E = 0 is around Γ Total ≈ 2 × 10 -3 , one order of magnitude higher than the diffusive flux. For low E the total flux is positive, indicating a flux from the edge to the core. The total radial flux becomes negligible for E > 2.

Comparing the magnitude of Γ Total and Γ D RW , it is unambiguously clear that the radial particle flux for low energies E < 1 is, in this streamer-like turbulence, dominated by non-diffusive processes, which accounts for ∼ 90% of the flux at E = 0. In addition to this statement, we find that the maximum magnitude of the total radial flux regardless the energy is one order of magnitude higher in this second simulation (∼ 4.8 × 10 -4 vs ∼ 2 × 10 -3 ).

III.3 Chapter conclusion

We added a test particle module to the reduced bounce-averaged gyrokinetic code TERESA which is focused on investigating low frequency phenomena of the order of the trapped particles toroidal precession frequency. The code can henceforth solve, in addition to the distribution function f and the electric potential φ, the individual trajectories of millions of test particles. The test particle positions in phase space are computed directly in the code thus allowing the same order of accuracy as the TERESA numerical scheme. Test particles are particles which respond to the electrostatic field without contributing to it. The addition of test particles in our code gives us access to information which where not available before with only f and φ. It allows to have better insights on transport phenomena such as diffusion, advection or ballistic motions.

In this first work using test particles in TERESA we aimed at separating the contribution of the diffusive process of the particles in the radial direction, from the total radial transport. To proceed, we initialized one million test particles at t = 6, in a turbulent core plasma, in the center of our box (ψ = 0.5) and let them evolve in the electrostatic potential, in two simulations. In the first simulation, only the trapped ions were treated fully kinetically, while the trapped electrons were "quasi-adiabatic", using the parameter δ m [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF]. In the second simulation, both trapped ions and electrons were fully kinetic.

In the first simulation, the turbulence is TIM-driven and there is no dominant zonal flows, streamers or large potential structure which would drastically impact the transport. Instead, the α-mode spectrum presents a bulk of most intense modes ranging from m ≈ 1 to m ≈ 10 and the ratio eφ/T ∼ 0.01 -0.03 is typical of core turbulence. This simulation is rather representative of developed turbulence.

We then calculated the time evolution of the test particles Mean Squared Displacement in the radial (ψ) direction, for each E ∈ [0; 20] and observed that the MSD tends to first have a rapid growth indicating that the test particles reorganize themselves inside the potential structure where they were initialized in. Then the test particle MSD grows linearly indicating a radial diffusion process toward the edges of the box until the test particles start to undergo boundaries effects. At E R = 1.74 which is approximately when the particles resonate with the TIM, the MSD becomes constant in time at MSD≈ 0.08 indicating that the test particles fully explored the box in the radial direction. At the resonant energy, the particles tend to "see" constant potential structure and thus explore the simulation box faster than at non-resonant energies.

With the MSD obtained for each E, we calculated a radial random walk diffusion coefficient, which presents a peak around the resonant energy E R ≈ 1.74.

Then we estimated the radial diffusive flux of the particles which is the flux if there was only a diffusion process. We compared it to the total flux obtained directly from the TERESA simulation accounting all the radial transport processes. We found that the radial particle transport is clearly dominated by the resonant particles, as both fluxes present a peak around E R ≈ 1.74. Both peaks are negative and of the same intensity, indicating that radial transport of resonant particles is exclusively a diffusive process toward the edge. It is coherent with our choice of "global turbulence". The high energy particle radial transport tend to be negligible. Below the resonance, for E < 1, the gradient of f is negative and the diffusive flux is oriented toward the core. Between E = 0 and E ≈ 0.5, there might be a non-diffusive process, directed toward the edge so that the total flux is lower than the diffusive flux. Between E ≈ 0.5 and E ≈ 1, a non-diffusive process appears to induce a flux directed toward the edge, so that the total flux is more intense than the diffusive flux.

In the second simulation, we use kinetic trapped ions and electrons. In this case, the nature of the turbulence is streamer-like, with large, radially elongated structures of the size of the simulation box, instead of the "global turbulence" from the previous case. This simulation is rather representative of a streamer-dominated turbulence, which could occur in experiments of linear magnetized plasmas. From t ≈ 3 to t ≈ 7, the turbulence is dominated by TIM, then for t > 7, the TEM are dominant.

We performed the same procedure as mentioned above, in order to separate the diffusive flux from the total radial flux. In this case, the MSD and the TPDF have similar behavior than in the previous simulation. Indeed the MSD first experiences a rapid growth in a short time (∆t ∼ 0.1) which is due to the fact the the test particles first explore the potential structure they were initialized in, which in the case of streamer-like turbulence, leads to a faster rate of flattening of the TPDF, compared to the previous simulation with a "global turbulence" case. Then the MSD grows linearly in time, until the TPDF becomes totally flat and the test particles have explored all the simulation box in ψ.

From the slope estimated from the MSD, we calculated a random walk diffusion coefficient. This diffusion coefficient presents similar features than in the previous case. We observe a peak around E = 1.38 which indicates that the particles resonate with the modes at a lower energy than in the previous case (E = 1.74). For high energy particles, the diffusion coefficient tends to 0, but is finite for thermal particles.

From this diffusion coefficient we estimate a diffusion flux, which we compare to the total radial particle flux. In this second simulation, we find large discrepancies between the two fluxes. The total flux is the most intense for E = 0, where the diffusive contribution to the total flux is only ∼ 10% and both fluxes are positive (from the edge to the core). At E = 0, the non-diffusive processes account for ∼ 90% of the total radial particle flux. At E = 1.38, the total flux is positive while the diffusive flux is negative (from the core to the edge). For high energy particles, similarly to the previous simulation, the radial particle fluxes vanish. We thus find that in this second simulation, the transport is governed by highly non-diffusive

Chapter IV

Radial density and heat fluxes description in the velocity space: Nonlinear simulations and quasi-linear calculations

In chapter III, we studied the radial transport of one million test particles and compared the radial diffusive flux to the non-diffusive flux, in the energy direction.

In this chapter, we will investigate the structure of the flux in energy with ∼ 10 times more grid points in the energy dimension. We noticed for instance an error of about 6% on the density and heat fluxes when we choose 128 points in energy instead of 1024 points. With 64 points in energy the density flux has an error of 25% and the heat flux an error of 10%. We then compare the results to quasi-linear theory. The results in this chapter lead to the publication of an article in Physics of Plasmas [START_REF] Médina | Radial density and heat fluxes description in the velocity space: Nonlinear simulations and quasi-linear calculations[END_REF].

Exploring the details of the velocity space with brute force gyrokinetic simulation, can consume up to 10 8 core.hours on a modern supercomputer [START_REF] Maeyama | Cross-scale interactions between electron and ion scale turbulence in a tokamak plasma[END_REF], even with a local code. The global gyro-bounce averaged code TERESA allows one to study the details of the trapped particle dynamics in the velocity space, without degrading the accuracy in the real space, and at reasonable numerical cost.

For collisionless trapped-particles-driven modes, such as the trapped-electron mode (TEM) or the trapped-ion mode (TIM), we expect the mode-particle resonance to play an essential role for the particles and heat fluxes, as shown in the previous chapter. It was found recently that this poses an issue for the accuracy required to describe the mode structure, which has a strong, narrow peak around the resonant energy [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF]. Trapped-ion modes being weakly dispersive, all resonances in the turbulent case could be very localized in the velocity (or energy) space. This raises the issue of the precision required to accurately describe the fluxes in the velocity space in gyrokinetic simulation when weakly dispersive modes are present. Usually, conventional gyrokinetic simulation focus on fine description of the real space, but are limited to only a handful of points in the velocity dimensions.

In this chapter we aim to describe in details the velocity dependence of the heat and density radial fluxes in nonlinear simulations from the gyro-bounce averaged code TERESA, in order to give an estimate of the accuracy required in the velocity space to describe the resonant fluxes. Here, we focus on trapped-ion modes, although TERESA can treat TEM as well.

To provide some reference point, the results from TERESA's nonlinear simulations will be compared to the quasi-linear theory. More than five decades after the two pioneering papers [START_REF] Vedenov | Quasi-linear theory of plasma oscillations[END_REF][START_REF] Drummond | Non-linear stability of plasma oscillations[END_REF], quasi-linear theory [START_REF] Adam | Effets des interactions résonnantes ondesparticules en turbulence faible des plasmas[END_REF][START_REF] Laval | Breakdown of quasilinear theory for incoherent 1-d langmuir waves[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF] is still relevant nowadays because in spite of a priori crude simplifications, its estimations of turbulent fluxes remains in good agreement with experimental results [START_REF] Doveil | Two regimes of self-consistent heating of charged particles[END_REF][START_REF] Bourdelle | Core turbulent transport in tokamak plasmas: bridging theory and experiment with qualikiz[END_REF] as well as with nonlinear gyrokinetic simulations [START_REF] Casati | Validating a quasi-linear transport model versus nonlinear simulations[END_REF][START_REF] Jenko | Heat and particle transport in a tokamak: advances in nonlinear gyrokinetics[END_REF].

The principal assumptions of the quasi-linear theory are a low level of fluctuations of the distribution function (i.e. f = f eq + δf and |δf | f eq where f eq describes a slowly evolving background distribution that is changing due to the effects of the unstable waves themselves) and that there is no trapped particle inside "potential wells" [START_REF] Vedenov | Quasi-linear theory of plasma oscillations[END_REF][START_REF] Bourdelle | Turbulent transport in tokamak plasmas: bridging theory and experiment[END_REF]. These assumptions can be questioned in the case of trapped particles driven modes, where wave/particles resonances play a central role. Quasilinear theory also relies upon the assumption that the correlation time τ C of the electric field seen by a resonant particle is small as compared to the evolution time of averaged quantities [START_REF] Adam | Reconsideration of quasilinear theory[END_REF]. In section IV.1, we give the simulation configuration and input parameters. In section IV.2 we show the radial density and heat flux from the TERESA simulation in the real space (ψ, corresponding to the minor radius), then we separate the contribution to the heat flux from thermal (low energy) particles and resonant particles. We also introduce quasi-linear theory with the purpose to give an element of comparison for the fluxes. Then in section IV.3 we explore in details the energy dependence of the radial particle flux from TERESA, and we compare it to the quasi-linear estimation.

IV.1 Simulation configuration

We focus on the ion dynamics. The grid in phase space is uniform, with N ψ points in ψ ∈ [0; 1], with ψ = 1 being in the center of the fusion plasma and ψ = 0 being toward the edge (but still in the core) and N α points in the toroidal precession angle α ∈ [0; 2π[. We choose N α × N ψ = 256 × 257. For the 2 invariants, the grid has N E points in E ∈ [0; 20]. It should be noted that such a range of energy is required for the convergence of the simulation results. A single value in κ = 0 [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF] δ m = 0.02 Table IV.2: Input parameters.

Grid Number of grid points Value

α N α = 256 α ∈ [0; 2π[ ψ N ψ = 257 ψ ∈ [0; 1] E, V N E or N V = 1024 E ∈ [0; 20]
Ω d = 1 C1 C 1 = 0.1 C2 C 2 = 0.1 φ φ = 1 Electron dissipation
which corresponds to deeply trapped particles [START_REF] Cartier-Michaud | Optimizing the parallel scheme of the poisson solver for the reduced kinetic code TERESA[END_REF] is chosen. The energy grid is finer close to E = 0 and looser for higher E allowing greater precision for low E. This is done using a new parameter V = √ E instead of the parameter E itself, so that the new grid in V still have N E points but is finer in terms of E for low V , corresponding to low E. Here we choose N E = 1024.

Here, the TERESA simulations are performed with thermal baths at both outside boundary ψ = 0 and inside boundary ψ = 1. Artificial dissipation is imposed in buffer regions ψ < 0.15 and ψ > 0.85 to smooth out the transition between turbulent fluctuations of φ, and the constraint φ = 0 at ψ = 0 and ψ = 1.

Quantity

Normalization Poloidal magnetic flux ψ = ψ L ψ Energy Ê = E T 0 Time t = ω 0 t Electric potential φ = φ ω 0 L ψ = R 0 a eφ T 0 Distribution function fs = 1 n 0,s 2π T 0,s ms 3 2 f s Table IV.3:
Recall of the normalization from section II.1.2

In the figures represented with a ψ axis, the buffer zones are for ψ ∈ [0; 0.15] and ψ ∈ [0.85; 1], and will be shaded in grey. We choose the ion Larmor radius ρ i = 0.01 and the ion banana width δ bi = 0.1. The initial radial temperature gradient is κ T = 0.15. The equilibrium distribution function f eq is chosen as locally Maxwellian in V (so exponential in E) as

f eq (ψ, E) = e -E [1 + (κ T (E -3/2) + κ n )ψ] (IV.1)
where κ n ≡ ∂ log(neq) ∂ψ | ψ=0 is a constant input parameter, which measures the equilibrium density gradient and similarly κ T measures the temperature gradient. The grid configuration, the input parameters and the normalization are recalled in tables IV.1, IV.2 and IV.3 respectively. This simulation took about one day of computation on 96 cores.

IV.2 Trapped-ion-mode turbulence and radial transport

In this section we analyze the evolution of the turbulent field, and the resulting density and heat fluxes obtained from the TERESA simulation. We observe a phase of linear growth of the plasma potential from t = 0 until t ≈ 4 -5. Linearly, the most unstable mode is the mode number m = 9, where φ = m φ m exp(imα). Its linear frequency and growth rate are ω m=9 = 15.6 and γ m=9 = 0.987. After t = 5 the modes achieve their saturation level. The saturation amplitude in terms of the root mean square φ rms is of the order of φ rms ∼ 0.05. In physical units, since we have chosen an inverse aspect ratio = 0.1, we obtain eφ rms /T ∼ 5 × 10 -3 . At t ≈ 5, the most intense modes are the modes around m = 10.

IV.2.1 Time evolution of dominant modes

For the rest of this chapter we investigate the time t = 5, at the beginning of the turbulent phase, where the dominant modes have just achieved saturation.

IV.2.2 Particle, density, and heat fluxes

Averaging the Vlasov equation over the angle α yields where δf = f -f α . The quantity Λ NL ψ is the nonlinear radial particle flux. It is calculated directly within the TERESA simulation. Taking the moments of Eq. (IV.2) yields the following conservation equations:

∂ f α ∂t + ∂ ∂ψ Λ NL ψ (ψ, E, t) = 0, (IV.
∂ n α ∂t + ∂Γ NL ∂ψ = 0, (IV.4) ∂ nT α ∂t + ∂q NL ∂ψ = 0, (IV.5)
with the density flux

Γ N L (ψ, t) = E Λ N L ψ (ψ, E , t) √ E dE , (IV.6)
or the heat flux

q N L (ψ, t) = E Λ N L ψ (ψ, E , t)E √ E dE . (IV.7)
The velocity space is linked to the energy space as [70, p. 44], therefore we add a √ E factor inside the integral over energy space to physically integrate over the velocities. Note that we distinguish the particle (in the sense of "test particle") flux Λ, which is also the flux of phase-space density, and the density (in the sense of fluid density) flux Γ. Λ has an E dependence while Γ and q are integrated over E. The artificial density and heat fluxes due to the buffers diffusion in TERESA are included. The negative sign of both Γ and q is consistent with a flattening of the initial gradients. The fluxes are bell shaped because of the boundary conditions restraining the instability to a 0 value at ψ = 0 and ψ = 1. Simulation with larger box size yielded flatter flux profile, but with higher numerical cost, so we focus on a "small box". In order to obtain more information about the origin of turbulent transport, we are now going to separate the contribution from thermal particles and from resonant particles, by integrating over different ranges of energy E. Thermal particles and resonant particles will be discriminated as follows: thermal particle range of integration will be E ∈ [0; 1.4] These integration intervals are chosen to take into account all the particles below the resonance energy for the thermal particles, and all the particles in the resonance peak for the resonant particles. These choices are coherent considering The resonant particle contribution accounts for more than 90% of the total heat flux. Furthermore, the contribution from thermal particles has the opposite sign. This can be simply explained by the sign of ∂ f α ∂ψ which is positive for E > 3/2 but negative for E < 3/2. Here we focus on heat, but similar conclusions are reached for the density flux.

d 3 v = C √ EdE with C being a constant

IV.2.3 Quasi-linear theory

Quasi-linear theory describes the slow evolution ( ∂ log f ∂t γ k with γ k the linear growth rate of the most unstable mode) of the α-averaged particle distribution function f α , as the solution of a diffusion equation [START_REF] Diamond | Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas[END_REF],

∂ f α ∂t = ∂ ∂ψ D QL ∂ f α ∂ψ (IV.8) (a)

IV.2.3.1 Quasi-linear diffusion coefficient

In this section, we detail the quasi-linear diffusion coefficient derivation for the TERESA model, adopted from the procedure developed in the case of the 1D Vlasov-Poisson model [START_REF] Sagdeev | Nonlinear plasma theory[END_REF]. The starting point is the Vlasov equation in normalized units, The Fourier transform of Eq. (IV.9) is detailed in Appendix B and yields

∂ fm ∂t + im EΩ d Z fm -im φm ∂ f0 ∂ψ = - l il fl ∂ φm-l ∂ψ + l =0 i(m-l) φm-l ∂ fl ∂ψ (IV.12)
The RHS corresponds to nonlinear wave-wave and wave-particle interactions. In the framework of quasilinear theory, the nonlinear wave-wave and wave-particle interactions are assumed to be weak according to the hypotheses of weak turbulence and of no particle trapped in electrostatic potential. Thus we neglect these nonlinear interactions except for m = 0 or l = m. Therefore, for m = 0, Eq. (IV.12) is approximated by

∂ fm ∂t + im EΩ d Z fm -im φm ∂ f α ∂ψ + im fm ∂ φ0 ∂ψ = 0 (IV.13)
The latter equation is of the form

L fm = g(ψ, E, t), (IV.14)
where L is the linear operator

L = ∂ t + iω R,m (ψ, E, t), with ω R,m (ψ, E, t) = m   Ω d E Z + ∂ φ0 ∂ψ   (IV.15)
and g = im φm ∂ ψ f α . For ω R,m = 0, the solution is trivial. For ω R,m = 0, Eq. (IV.14) can be solved using Green's function

G E (t, s) = exp [iω R,m (s -t)],
which is such that

LG E = δ(s -t). The solution is fm (ψ, E, t) = t 0 e iω R,m (s-t) im φm (ψ, E, s) ∂ f α ∂ψ (ψ, E, s) ds (IV.16)
It turns out that the latter expression is also valid for ω R,l . Each mode is assumed to have a fixed frequency ω m (obtained from linear theory), and a time-dependent growth rate γ m , φm (ψ, E, t) = φm (ψ, E, 0) exp t 0 (-iω m + γ m (t )) dt (IV.17)

For m = 0, Eq. (IV.12) simplifies as

∂ f α ∂t + m im ∂( fm φ * m ) ∂ψ = 0 (IV.18)
Substituting Eqs. (IV. [START_REF] Fischer | Anthropogenic contribution to global occurrence of heavy-precipitation and high-temperature extremes[END_REF]) and (IV.17) yields

∂ f α ∂t = ∂ ∂ψ m m 2 φm (ψ, E, t) 2 t 0 e i(ω
∂ f α ∂t = ∂ ∂ψ D QL ∂ f α ∂ψ (IV.20)
In other words, Eq. (IV.20) is obtained by substituting the linear response for φ and δf in Eq. (IV.3), which implies neglecting some of the nonlinear coupling terms. The main hypotheses of quasi-linear theory are: weak turbulence (small fluctuations of the profiles compared to the equilibrium), no trapped particles in electrostatic potential 'wells' (with a very large number of electrostatic waves, the resonance region of each waves can overlap so that the particle motion becomes stochastic, and the particle can wander in the velocity phase space [START_REF] Diamond | Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas[END_REF][START_REF] Sagdeev | Nonlinear plasma theory[END_REF]) and a small auto-correlation time of the electric field compared to the evolution time of the profiles. The quasi-linear radial particle flux Λ QL ψ (ψ, E, t) is the product of the quasi-linear diffusion coefficient D QL and the mean gradient

Λ QL ψ (ψ, E, t) = -D QL (ψ, E, t) ∂ f α ∂ψ (ψ, E, t) (IV.21) with D QL (ψ, E, t) = m m 2 φm (ψ, E, t) 2 1 -e -i(ω R,m -ωm)t-|γm|t i(ω R,m -ω m ) + |γ m | (IV.22) and ω R,m (ψ, E, t) = m   Ω d E Z + ∂ φ0 ∂ψ   (IV.23)
where the sum m is over the m components of the Fourier decomposition in α of the gyro-bounce averaged electrostatic potential, noted φm (ψ, E, t). γ m and ω m are respectively the growth rate and linear frequency of each trapped-ion mode m, ω R,m (ψ, E, t) is homogeneous to a frequency and ∂ φ0 ∂ψ accounts for the Doppler effect from the mean flow (including zonal flow), and the trapped particles precession frequency Ω d and Z are equal to 1 in our normalized units.

The use of this version of quasi-linear theory is justified for a time ∼ 1/γ after the start of the saturation phase. For later times, renormalized quasi-linear theory is needed [38, p. 181]. However, in this simulation, the fastest growing mode has a growth rate γ ≈ 1 and the saturation phase begins at t ≈ 5, thus the use of quasi-linear theory is justified for our study around t ≈ 5. The moments of Eq. (IV.21) yields the QL density flux

Γ QL (ψ, t) = E Λ QL ψ (ψ, E , t) √ E dE (IV.24)
or the QL heat flux

q QL (ψ, t) = E Λ QL ψ (ψ, E , t)E √ E dE (IV.25)
To compute the quasi-linear fluxes, the electrostatic potential φ, the mean distribution function f α , the modes growth rate γ m and the modes frequency ω m are needed.

In order to test fundamental aspects of the QL formalism, we want to limit the sources of discrepancies between QL theory and NL simulations. We thus take φ and f α directly from the TERESA simulation, instead of using approximations such as the mixing length. γ m and ω m are obtained from the resolution of the linear dispersion relation [START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF][START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF].

D QL includes the contribution of the zonal flow. Because the zonal radial electric field ∂ φ0 ∂ψ fluctuates rapidly in ψ and in time, which can create artifact in the resonance considering its presence in the denominator of Eq. (IV.22), we average this quantity over time and ψ.

Considering we are only interested in the time t = 5, we thus perform a time averaging over t ∈ [4.8; 5.2]. This time window (0.4) corresponds to a typical oscillation time of ∂ φ0 ∂ψ . We then perform a moving average in ψ with a 100 points window. This window in ψ corresponds to approximately a third of the total box size in ψ.

After these averages, the quantity ∂ φ0 ∂ψ has only a ψ dependence and can take positive or negative values. For example, for the highest oscillation amplitude of The NL fluxes and QL predictions are in qualitative agreement. However, there are significant quantitative discrepancies, even as we have chosen the best case scenario for QL theory (a small auto-correlation time, φ and ∂f ∂ψ directly from simulation and no mixing length estimate). The quasi-linear estimation of the density flux, the heat flux, the thermal particles contribution to the heat flux and the resonant particle contribution to the heat flux respectively present a 16%, 16%, 50% and a 20% discrepancy compared to the nonlinear flux obtained from the TERESA simulation.

In appendix C, we present a comparison of the neglected or kept terms in the quasi-linear equation. We find that two key terms contribute to the radial transport : the linear term ∂φ ∂α ∂ f α ∂ψ which is kept in the quasi-linear equation and the nonlinear term ∂φ ∂α ∂δf ∂ψ which is neglected. A comparison term by term shows that the neglected nonlinear term which contributes to the radial transport is in fact greater than the kept linear term. This may explain the discrepancy between the nonlinear and quasi-linear fluxes.

The next section is focused on the investigation of the details in the velocity (or energy) space of theses fluxes.

IV.3 Anatomy of the radial fluxes in the energy space

TERESA allows us to investigate in detail the velocity (or energy) space at reduced computational cost. In this section we will investigate the energy dependence of the particle flux Λ NL (ψ, E, t), at ψ = 0.5 and at time t = 5, taking advantage of our large number of grid points in energy, N E = 1024. In order to have some reference point, we will then compare Λ NL to the quasi-linear particle flux Λ QL .

The quasi-linear diffusion, Eq. (IV.22), presents a resonance for each mode m,

when ω R,m (ψ, E R , t) = ω m , with E R = Z Ω d ωm m -∂ φ0
∂ψ . This peak is the result of the resonance between the low frequency trapped ions modes and the precession of trapped particles. The position of the resonance for each mode is mainly determined by the ratio ωm m . The mean flow contribution ∂ φ0 ∂ψ slightly shifts the position of the resonance as discussed in section IV.2. The contribution of the most intense modes to the fluxes will be greater, so that the peak position will be mainly determined by these most intense mode at t = 5 (see Fig. As expected, the nonlinear (from TERESA simulation) fluxes in the energy (or velocity space), at the middle of the box and at t = 5, present a narrow resonance negative peak. The resonance peak width is about ∆E ≈ 0.15 (obtained by a Gaussian fit of the resonance peak, Λ NL ∼ exp[(E -E 0 )/∆E]). Strikingly, the resonance peak accounts for most of the flux. Indeed the flux at the resonance is one order of magnitude stronger than the flux for thermal particles, while the high energy particle flux tends to 0. We must keep in mind that the position of this resonance peak is determined by the most intense electrostatic modes at this ψ and time, and can be slightly shifted by the intensity of the zonal radial electric field ∂ φ0 ∂ψ .

Fig. (IV.5) includes the QL particle flux Λ QL as well. The NL and QL fluxes are in qualitative agreement, in terms of sign, energy of sign reversal, position of peak, width of peak, and behavior at large energies. However, quantitatively, the discrepancy between the nonlinear fluxes and the quasi-linear estimation is of 41% at E = 0, 55% at E = 1, and 57% at the resonance peak. As discussed in subsection IV.2.3, appendix C shows that out of two terms contributing to the radial transport, only the linear term is kept in the quasi-linear theory while the nonlinear term is neglected. The term by term comparison shows that the neglected nonlinear term is greater than the kept linear term, which may explain this discrepancy.

IV.4 Chapter conclusion

We performed a numerical simulation of trapped ion mode turbulence with the bounce-averaged gyrokinetic TERESA code. It is based on a reduced model, which is meant to investigate fundamental mechanisms and trends, rather than provide realistic quantitative predictions for tokamaks. The trapped ions are treated kinetically, while the passing ions respond adiabatically. The electrons respond quasi-adiabatically. The simulation parameters correspond to a small aspect ratio (0.1), a radial box size of 10 banana widths, a flat density profile, and a moderate temperature gradient (gradient length of 67 banana widths). The electric potential in terms of root mean square φ rms saturates to a value eφ rms /T 0 ≈ 5 × 10 -3 .

To investigate how the radial turbulent transport depends on the energy dimension, the simulation was performed with a large number of grid points N E = 1024. We focused on a radial location in the middle of the simulation box, and on a time in the simulation where the dominant trapped-ion modes have just achieved saturation, corresponding to the beginning of the turbulent phase.

As expected, the radial flux is negative for E < 3/2 and positive for E > 3/2, consistent with a flattening of the radial gradient of the distribution function. We emphasize that the flux features a narrow peak in the resonant region, with a width ∆E ≈ 0.15, around the energy E ≈ 1.8. Strikingly, this resonant peak accounts for 90% of the density and heat fluxes. In contrast, the contribution from thermal particles is negligible. Based on these results, a fine mesh in the energy space (δE 0.1), in the resonant region, is required to accurately describe the radial transport of density and heat.

The quasilinear predictions for radial fluxes (including the effect of zonal flow) are in qualitative agreement with the simulation results, in terms of global structure in the radial direction, sign throughout the energy dimension, behavior at small and large energies, and for the resonant peak in terms of its shape, location and width in the energy dimension. However, quantitatively, there is a 57% overprediction at the peak, and a 55% overprediction at the thermal energy E = T 0 .

Since the flux is negative for resonant energies, and positive for thermal energies, these discrepancies can compensate each others. Indeed, the total density and heat fluxes predicted by quasi-linear theory are in good quantitative agreement with the nonlinear simulation result, with only 16% overprediction. However, this agreement is only a result of compensating errors. We showed that the non-zonal nonlinear part of the radial advection in the Vlasov equation is actually slightly larger than the linear part. This can explain the discrepancy, since quasi-linear theory discards the non-zonal nonlinear part to obtain the response of the perturbed distribution function f -f α . Here, we should emphasize that we have chosen to maximize the chances of success of quasilinear theory. Indeed, both electric potential and mean distribution were taken directly from the nonlinear simulation, and substituted into the quasilinear formula for the flux. We did not use theory-based assumptions for the potential spectrum and amplitude, nor did we use an analytic equilibrium for the mean distribution. In this sense, we did not really test a prediction, but rather we have tested the assumptions underlying quasilinear theory.

In this chapter we focused on the trapped ion mode (TIM), but we expect similar results for the trapped electron mode (TEM). However, our analysis does not include the effects of neither ion-temperature-gradient (ITG) nor electrontemperature-gradient (ETG) modes, since the kinetic description in TERESA is limited to trapped particles. We believe this is the main caveat to our analysis, since, for example in ITG turbulence, resonances may not play such a crucial role in the transport, and therefore the co-existence and/or the coupling between TEM and ITG may mitigate the importance of the resonant peak in the radial fluxes.

Another important caveat is that boundary conditions consist of thermal baths, which strongly restrict the evolution of the profiles. As a result, the heat flux remains small in amplitude, and the relaxation of the initial temperature gradient is only marginal. A future analysis based on flux-driven simulations may provide new information about the role of resonant particles in turbulent transport.

Outlooks

In this work we adopted a model focused on turbulence driven by trapped particles. We averaged out the fast cyclotron motion of the charged particles, as well as the bounce motion of the trapped particles, while the passing particles respond adiabatically to the electric potential φ. Using an adapted action-angle formalism, it allows us to reduce the dimensionality of the problem from 6D to 4D (2 variables and 2 adiabatic invariants). The bounce-averaged electrostatic gyrokinetic TERESA code is based on this model and thus focuses on low frequency phenomena, of the order of the trapped particle toroidal precession time-scale and is not aimed to quantitatively predict transport in existing or futures tokamaks. Instead its purpose is to investigate general trends and fundamental ingredients of turbulent transport in a qualitative way. TERESA solves the system of coupled Vlasov-Quasi-neutrality equations, for the banana-center distribution function of species s: f s and for the electric potential φ.

We presented in this thesis the implementation in the TERESA code of a newly developed test particle module. The test particle trajectories are computed at order 2 in time (∆t 2 ) which is the same level of accuracy as the Semi-Lagrangian solver for the Vlasov-Quasi-neutrality system. TERESA can now follow the trajectories of millions of test particles which can give useful insights on transport processes.

As a first goal we aimed to estimate the diffusive contribution to the total radial particle flux. To achieve this goal we initialized one million test particles on every energy grid point and let them evolve in two turbulent simulations.

The first simulation was a case of "global turbulence", without dominant large potential structure such as zonal flows or streamers, thus favoring a diffusive behavior of the transport. Here only the trapped ions were fully kinetic, while the trapped electrons were "quasi-adiabatic", using the parameter δ m [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF]. We calculated a random walk diffusion coefficient in energy space from the slopes of the test particles MSD and then estimated a diffusive flux, which we compared Outlooks to the total radial particle flux. We found that the diffusive flux accounts for almost the entirety of the total flux, except for thermal particles (low E).

For the second simulation, we treated kinetically the trapped ions as well as the trapped electrons. The resultant turbulence was dominated by TIM for 3 t 7 and then by TEM for t 7. Also, we observed the presence of streamer-like, large elongated radial structures. From the slopes of the MSD we estimated a diffusion coefficient, which presented similar features with the previous simulation: a large peak at a resonant energy, E = 1.38 and a vanishing for high energy particles. Although a difference with the previous simulation was that thermal particles presented a finite diffusion coefficient. The random walk diffusive radial particle flux estimated from this diffusion coefficient was then compared to the total radial flux. We found that for E = 0, the total flux was the most intense and that diffusive processes accounted only for ∼ 10% of the radial transport. This is coherent with our streamer-like dominated turbulence.

In the second part of this thesis work, we investigated the radial flux with a finer mesh in energy space (N E = 1024) and compared the results with quasi-linear theory predictions. We found that the accurate description of the radial particle flux diminished with lower number of grid points in the energy dimension, which can be a challenge for heavier, more complete, gyrokinetic codes. The quasi-linear prediction hold a qualitative (∼ 50% discrepancy) agreement with the nonlinear fluxes around t = 5.

On a follow up of this work, we recently found that the discrepancies between quasi-linear predictions and the nonlinear fluxes rise at longer times, up to the point where there is no more qualitative agreement. This is thought to be due to small-scale structures in energy space. Since, at later times, the flux can be dominated by the contribution of energies significantly larger than the resonant energies, we do not expect agreement even with renormalized quasi-linear theory. In fact, strong particle trapping implies a departure from quasi-linear theory hypotheses. This underlines the importance of an accurate phase space description in energy space.

To take into account this fine energy space structures we could increase the energy grid points in the test particle simulation. This increase could change the transport substantially. An increase from N E = 256 to N E = 1024 would mean 4 times the computational time, or 24 days on 256 cores for the fully kinetic trapped particle simulation with millions of test particles. Additionally, the turbulence from section III.2 with fully kinetic trapped particles was streamer-like dominated. An interesting investigation would be to find a "global turbulence" case similar to the first simulation in section III.1.2, but with both kinetic ions and kinetic electrons, which may include both diffusive and nondiffusive parts.

Further work on chapter IV would be to include fully kinetic trapped particles. Indeed, here we used fully kinetic ions and "quasi-adiabatic" electrons, and as we seen in section III.2, adding kinetic electrons can change the nature of the transport.

= max( n+1 f,i,j ) (A.5) with n+1 f,i,j = f n+1 i,j -f (α i , ψ j , t n+1 ) (A.6) Putting (A.3) into (A.6) we find n+1 f,i,j = f (α i -u(α i , t n )∆t, ψ j -v(ψ j , t n )∆t) -f (α, ψ, t n+1 ) (A.7)
which we can develop using Taylor expansion in

n+1 f,i,j = -∆t ∂f ∂t + u(α i , t n ) ∂f ∂α + v(ψ j , t n ) ∂f ∂ψ + ∆t 2 u 2 (α i , t n ) 2 ∂ 2 f ∂α 2 + v 2 (ψ j , t n ) 2 ∂ 2 f ∂ψ 2 + u(α i , t n )v(ψ j , t n ) ∂ ∂α ∂f ∂ψ - 1 2
∂ 2 f ∂t 2 + O(∆t 3 ) (A.8) and considering that ∂f ∂t + u(α, t n ) ∂f ∂α + v(ψ, t n ) ∂f ∂ψ = 0 (A.9)
we then obtain

n+1 f,i,j = ∆t 2 u 2 (α i , t n ) 2 ∂ 2 f ∂α 2 + v 2 (ψ j , t n ) 2 ∂ 2 f ∂ψ 2 + u(α i , t n )v(ψ j , t n ) ∂ ∂α ∂f ∂ψ - 1 2
∂ 2 f ∂t 2 + O(∆t 3 ) (A.10)
Analytically, we find that the local error in time of the numerical scheme is of second order ∆t 2 . Now we aim to find the same result with simulation data, using the adiabatic invariant E. In order to find the accuracy order of the numerical scheme we proceed as follows: we start a reference simulation with a constant time-step ∆t considered to be "small enough", typically we choose ∆t ref = 5 × 10 -6 . Then we start the same simulation but with different constant time-step ∆t > ∆t ref , typically we take ∆t ∈ [1 × 10 -5 ; 5 × 10 -5 ; 1 × 10 -4 ; 5 × 10 -4 ; 1 × 10 -3 ; 5 × 10 -3 ; 1 × 10 -2 ]. For each simulation, at a time 1 × 10 -2 , we calculate the relative error of a E, compared to the E ref of the reference simulation.

Appendix B Vlasov equation in Fourier space

In this appendix we derive the Fourier transform of the Vlasov equation, used in section IV.2.3.1. We start with the Vlasov equation:

∂ f ∂t +   Ω d E Z s + ∂ φ ∂ψ   ∂ f ∂α - ∂ φ ∂α ∂ f ∂ψ = 0 (B.1)
and we apply a Fourier transform in the α direction:

1 2π 2π 0   ∂ f ∂t +   Ω d E Z s + ∂ φ ∂ψ   ∂ f ∂α - ∂ φ ∂α ∂ f ∂ψ   e -imα dα = 0 (B.2)
The first term is simply

Θ 1 = 1 2π 2π 0 ∂ f ∂t e -imα dα = ∂ fm ∂t (B.3)
The second term is

Θ 2 = 1 2π 2π 0 Ω d E Z s ∂ f ∂α e -imα dα (B.4)
which becomes, after an integration by parts:

Θ 2 = 1 2π Ω d E Z s      f e -imα 2π 0 =0 + 2π 0 im f e -imα dα      = im Ω d E Z s fm (B.5)
The third term is

Θ 3 = 1 2π 2π 0 ∂ φ ∂ψ ∂ f ∂α e -imα dα (B.6) 109
Nevertheless,

Θ 3 = F   ∂ φ ∂ψ ∂ f ∂α   = F   ∂ φ ∂ψ   * F   ∂ f ∂α   (B.7)
where * is the convolution product operator and F is the Fourier transform. Thus,

F   ∂ φ ∂ψ   = ∂ φ ∂ψ (B.8)
and

F   ∂ f ∂α   = im fm (B.9) therefore: Θ 3 = i f ∂ φm- ∂ψ (B.10)
Finally the fourth term is

Θ 4 = F   - ∂ φ ∂α ∂ f ∂ψ   = -F   ∂ φ ∂α   * F   ∂ f ∂ψ   = -im φm *   ∂ f ∂ψ   = - i φ ∂ fm- ∂ψ (B.11) or Θ 4 = - i(m -) φm- ∂ f ∂ψ = - =0 i(m -) φm- ∂ f ∂ψ -im φm ∂ ¯ f α ∂ψ (B.12)
Thus, the Fourier transform of Eq. (IV.9) yields

∂ fm ∂t + im EΩ d Z fm -im φm ∂ f0 ∂ψ = - i f ∂ φm- ∂ψ + =0 i(m-) φm- ∂ f ∂ψ (B.13)
where we removed the double bar notation for clarity.

of the simulation and as a function of the particles energy.

L 1 is not represented because choosing an adequate timestep in order to evaluate ∂ fm ∂t can be a tricky task considering the timestep taken for the simulation. However it is possible to estimate the magnitude of this term.

From IV.12 we have L 1 + L 2 + L 3 = NL 1 + NL 2 + NL 3 . We observe on Figure C.1 that L 2 L 3 , NL 1 , NL 2 , NL 3 (except for E 1), thus, we deduce that L 1 is of the same order of magnitude as L 2 and that the two terms cancel each other out. L 2 corresponds to the banana precession around the toroidal direction and does not impact transport in the radial direction nor the QL radial fluxes. L 3 contributes directly to the radial QL fluxes since it involves the gradient of f α in the ψ direction. NL 1 involves the zonal radial electric field ∂ φ0 ∂ψ and is the one responsible for the resonance peak shift in the energy/velocity dimension. NL 2 and NL 3 correspond to all the other nonlinear coupling that do not involve the zonal flow. NL 3 contributes directly to the radial transport since it involves the gradient of f α in the ψ direction but is neglected in the framework of the quasi-linear theory.

The key terms influencing the radial fluxes are thus L 3 and NL 3 but only L 3 is taken into account. We observe that the key term influencing the radial transport, kept in the quasi-linear theory (L 3 ) is small compared to the neglected nonlinear coupling term also influencing the radial transport (NL 3 ), especially around the resonance energy. The other neglected term (NL 2 ) is in contrast small compared to the other non neglected term. From this, one can expect large discrepancy between nonlinear simulations and quasi-linear theory.

Résumé en français

La lutte contre le réchauffement climatique est le challenge du XXI e siècle [START_REF] Cox | erratum: Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model[END_REF][START_REF] Thomas | Extinction risk from climate change[END_REF][START_REF] Pounds | Widespread amphibian extinctions from epidemic disease driven by global warming[END_REF][START_REF] Heimann | Terrestrial ecosystem carbon dynamics and climate feedbacks[END_REF][START_REF] Peters | The challenge to keep global warming below 2 c[END_REF][START_REF] Revesz | Global warming: Improve economic models of climate change[END_REF][START_REF] Berg | Landatmosphere feedbacks amplify aridity increase over land under global warming[END_REF][START_REF] King | Australian climate extremes at 1.5 c and 2 c of global warming[END_REF][START_REF] Russo | Half a degree and rapid socioeconomic development matter for heatwave risk[END_REF]. Un effort international a été entreprit pour limiter le réchauffement climatique sous les +2 • C en 2100, comparé aux niveaux préindustriels des années 1850 [START_REF] Rogelj | Paris agreement climate proposals need a boost to keep warming well below 2 c[END_REF]. Cependant, il est probable que les +2 • C soient atteints bien avant 2100 [START_REF] Raftery | Less than 2 c warming by 2100 unlikely[END_REF][START_REF] Tollefson | Ipcc says limiting global warming to 1.5 c will require drastic action[END_REF]. En effet, les estimations actuelles tablent sur une augmentation de la température de +3 • C à +5 • C en 2100 [START_REF] Knutti | A scientific critique of the two-degree climate change target[END_REF]. Le réchauffement climatique causé par l'Homme a déjà atteint +0.93 • C (±0.13 • C) en 2015 comparés aux années 1850 [START_REF] Millar | Emission budgets and pathways consistent with limiting warming to 1.5 c[END_REF][START_REF] Allen | Warming caused by cumulative carbon emissions towards the trillionth tonne[END_REF][START_REF] Fischer | Anthropogenic contribution to global occurrence of heavy-precipitation and high-temperature extremes[END_REF]. Les études les plus récentes ont démontré la relation entre la hausse globale des températures et l'émission de gaz à effet de serre dans l'atmosphère [START_REF] Mcglade | The geographical distribution of fossil fuels unused when limiting global warming to 2 c[END_REF][START_REF] Matthews | The proportionality of global warming to cumulative carbon emissions[END_REF]. Pour limiter le réchauffement à +2 • C, les émissions de CO 2 doivent être réduites de 50 à 80% en 2050, comparé aux émissions de 2010 [START_REF] Knutti | A scientific critique of the two-degree climate change target[END_REF] tandis que la cible des +1.5 [START_REF] Horton | Drift waves and transport[END_REF][START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]. Ces ondes électromagnétiques (considérées électrostatiques lorsque β 1) peuvent être déstabilisées par des gradients de densité et de température. Dans ces cas, le potentiel électrique φ augmente en amplitude jusqu'à atteindre un niveau de saturation. Les microinstabilités les plus communes sont les instabilités de gradient ionique (ITG), de gradient électronique (ETG), les modes d'électrons piégés (TEM) et d'ions piégés (TIM) qui appartiennent à la famille des ondes de dérive. Les micro instabilités du plasma peuvent donner lieu à des structures à grandes échelles telles que les écoulements zonaux [START_REF] Lin | Turbulent transport reduction by zonal flows: Massively parallel simulations[END_REF][START_REF] Rogers | Generation and stability of zonal flows in ion-temperature-gradient mode turbulence[END_REF][START_REF] Kim | Dynamics of zonal flow saturation in strong collisionless drift wave turbulence[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF] ou les streamers [START_REF] Dimits | Formation of streamers in plasma with an ion temperature gradient[END_REF][START_REF] Champeaux | Streamer and zonal flow generation from envelope modulations in drift wave turbulence[END_REF][START_REF] Yamada | Anatomy of plasma turbulence[END_REF].

Afin d'étudier le transport turbulent dans les tokamaks, on peut utiliser un modèle cinétique qui décrit le plasma via une fonction de distribution f à 6 dimensions, couplée aux équations de Maxwell. Dans ce travail nous utilisons l'approche cinétique parce que nous voulons décrire des phénomènes qui impactent la distribution en vitesse. Dans cette approche, l'évolution temporelle de f est gouvernée par l'équation de Vlasov couplée aux équations de Maxwell pour les champs macroscopiques. Ces équations non linéaires couplées ne présentent pas de solution analytique pour les cas non triviaux.

Avec le développement du calcul haute performance (HPC), il est possible d'approximer les solutions de ces équations grâce à des supercalculateurs, en utilisant des modèles statistiques réduits. Une de ces approches est de garder uniquement les ingrédients fondamentaux de la turbulence. Le rôle de ces codes basés sur des modèles réduits est d'ouvrir la voie pour d'autres codes plus complets mais plus lourds. Dans ce travail nous utilisons le code TERESA qui est à mi-chemin entre les codes lourds gyrocinétiques et les modèles réduits. C'est un code 4D (2 variables + 2 invariants adiabatiques) moyenné sur le mouvement de rebond. Le code TERESA se focalise sur la dynamique des particules piégées.

Dans un tokamak, le champ magnétique est plus intense vers le centre que vers les bords. Ainsi, les particules ressentent une barrière de potentiel lorsqu'elles approchent du côté du tore où le champ est le plus fort. Cela sépare les particules en deux classes : les particules passantes, et les piégées. Les trajectoires des particules passantes s'enroulent complètement autour du tore en suivant les lignes de champ magnétiques. Les particules piégées n'ont pas assez d'énergie pour franchir le côté au fort champ magnétique et donc, finissent par être repoussées dans le sens opposé. Elles sont donc piégées du côté où le champ magnétique est faible et leur mouvement de rebond ressemble à celui d'une banane lorsqu'il est projeté dans le plan poloïdal. Ce mouvement de rebond a une fréquence plus faible que le mouvement cyclotron d'au moins un ordre de grandeur. Les particules piégées (les "bananes") précessent lentement (comparé à la fréquence de rebond) dans la direction toroïdale.

Dans TERESA [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF][START_REF] Ghizzo | Streamer-induced transport in the presence of trapped ion modes in tokamak plasmas[END_REF][START_REF] Cartier-Michaud | Optimizing the parallel scheme of the poisson solver for the reduced kinetic code TERESA[END_REF][START_REF] Cartier-Michaud | Staircase temperature profiles and plasma transport self-organisation in a minimum kinetic model of turbulence based on the trapped ion mode instability[END_REF][START_REF] Drouot | Étude de la turbulence liée aux particules piégées dans les plasmas de fusion[END_REF][START_REF] Médina | Radial density and heat fluxes description in the velocity space: Nonlinear simulations and quasi-linear calculations[END_REF], le mouvement de rebond des particules piégées et le mouvement cyclotron sont moyennés. De plus, on considère que les particules passantes réagissent de façon adiabatique au potentiel électrique, de sorte que le mouvement des particules piégées dépende linéairement du potentiel électrique. Ainsi, seule la turbulence aux échelles de temps du mouvement de précession des particules piégées est retenue. Le mouvement des particules piégées, qui sont donc de facto des "bananes piégées", est décrit par l'équation de Vlasov couplée à l'équation de Poisson. Se focaliser sur la turbulence générée par les particules piégées est pertinent puisque ces modes jouent un rôle important dans le transport [START_REF] Horton | Drift waves and transport[END_REF]. Dans TERESA, seuls les TIM et les TEM peuvent apparaître.

Une partie du travail présenté dans cette thèse était de développer un module de particules test et de l'implémenter dans le code TERESA. En effet, résoudre la fonction de distribution ne donne pas directement les trajectoires des particules. Connaitre ces trajectoires permet d'avoir une meilleure compréhension de diverses phénomènes de transport comme la diffusion [START_REF] Rechester | Electron heat transport in a tokamak with destroyed magnetic surfaces[END_REF][START_REF] Isichenko | Invariant measure and turbulent pinch in tokamaks[END_REF][START_REF] Harvey | Radial transport and electron-cyclotron-current drive in the tcv and diii-d tokamaks[END_REF][START_REF] Estève | Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport[END_REF], les sous ou hyperdiffusions [START_REF] Misguich | Noble internal transport barriers and radial subdiffusion of toroidal magnetic lines[END_REF], l'advection [START_REF] Isichenko | Invariant measure and turbulent pinch in tokamaks[END_REF][START_REF] Boedo | Flow reversal, convection, and modeling in the diii-d divertor[END_REF][START_REF] Boedo | Transport by intermittent convection in the boundary of the diii-d tokamak[END_REF][START_REF] Dif-Pradalier | On the validity of the local diffusive paradigm in turbulent plasma transport[END_REF], le piégeage de particules dans des puits de

Figure 1 :

 1 Figure 1: log-log diagram of the fusion cross section as a function of energy, for different reactions. Picture taken from [26, p. 5].

Figure I. 1 :

 1 Figure I.1: Helical trajectory of a charged particle with a non zero initial velocity parallel component, in a uniform magnetostatic field. Picture taken from [85].

  I.2.

Figure I. 2 :

 2 Figure I.2: A tokamak magnetic field line, in red, is twisted around the torus. Picture taken from [70, p. 15].

I. 3 .

 3 The magnetic field lines are thus helical, twisting around the torus.

Figure I. 3 :

 3 Figure I.3: Torus coordinates where ( e r , e θ , e ϕ ) form a direct basis. r is the radial coordinate, θ the poloidal angle, ϕ the toroidal angle, R 0 is the major radius of the torus (from the toroidal axis Z to the center of a poloidal section) and R is the major radius coordinate. Picture taken from[START_REF]Toroidal coordinates[END_REF].

Figure I. 4 :

 4 Figure I.4: The toroidal magnetic field is created from electric current flowing in external coils. Picture taken from [26, p. 16].

B ϕ = 5 . 3

 53 Tesla, a poloidal magnetic field B θ one tenth of the toroidal one so that B ≈ B ϕ , and v ⊥ = v ths = 2k B Ts ms the thermal velocity of the species s, k B the Boltzmann constant. Thus we obtain : proton electron ρ cs (m) 3.1 × 10 -3 7.3 × 10 -5 ω cs (rad/s) 5.1 × 10 8 9.3 × 10 11In a tokamak, we usually consider that the particles respond adiabatically to magnetic fluctuations (ω < ω c ), meaning that B is nearly constant over a cyclotron period ( 1 B ∂B ∂t |ω c |), and that its characteristic gradient length is very large compared to the Larmor radius (

  the force I.24 is directed toward the opposite direction of the parallel gradient of B.

Figure I. 5 :

 5 Figure I.5: The trapped particle trajectory forming a "banana" when projected in a poloidal section of the tokamak. Picture modified from [96, p. 8].

  proton electron δ bs (m) 4.6 × 10 -2 1.1 × 10 -3 ω bs (rad/s) 3.4 × 10 4 1.5 × 10 6

Figure I. 6 :

 6 Figure I.6: Toroidal precession of a trapped particle. The magnetic field is the lowest toward θ = 0 and the highest toward θ = ±π. The particle is trapped between the two magnetic mirrors m1 and m2 at the poloidal angles θ = ±θ crit . This bounce trajectory does not close on itself, as it slowly (compared to the bounce motion) precesses here toward the right which is the toroidal direction ϕ. Picture taken from [25, p. 140].

  Fig I.7.

Figure I. 7 :

 7 Figure I.7: Neo-classical estimate (dashed line) and experimental measurements (circles and squares) of the ion and electron heat diffusivities in the tokamak Tore Supra [103, p. 28].

Figure I. 8 :

 8 Figure I.8: Energy spectrum in 3D, homogeneous, isotropic turbulence. In the inertial range, the kinetic energy is transferred at constant rate to smaller and smaller vortices, until the size of the vortices reaches the dissipation scale where energy is lost due to viscosity effects. Picture taken from [106, p. 6].

  [34, p. 92]. As the energy transfer goes to smaller and smaller scales, the vortices eventually reach a critical scale where the viscosity of the fluid (or plasma) becomes dominant compared to inertial effects. This scale L d = 2πk -1 d is the dissipation scale. Fig.I.8 shows the energy spectrum E(k) in 3D, homogeneous, isotropic turbulence: the kinetic energy is injected at a scale k < k L (large scale) then the kinetic energy is transferred to smaller scales at a constant rate in the inertial range, and finally the kinetic energy is lost in the dissipation scale (or Kolmogorov scale).

Figure I. 9 : 1 I

 91 Figure I.9: Energy spectrum in 2D, homogeneous, isotropic turbulence. The kinetic energy is injected at a scale L I = 2πk -1 I . From k -1 I we observe a double cascade: one direct cascade of enstrophy toward the smaller scales (k I < k) as E(k) ∝ k -3 and one inverse cascade toward the larger scales (smaller k < k I ) as E(k) ∝ k -5/3 . Picture taken from [113].

Figure I. 10 :

 10 Figure I.10: Instability caused by a temperature gradient aligned with the magnetic field gradient, on the "unfavorable curvature region" of the tokamak. Here the particle following the ∇B drift are positively charged. Picture taken from [70, p. 75].

Figure I. 11 :

 11 Figure I.11: Example of stability diagram of ITG and TEM instabilities at τ = T i /T e = 1. R is the major radius of the tokamak. Picture taken from [126].

2 )

 2 and of the Poisson equation which closes the system :

  [70, p. 40] [73] :

  3.1, we need another equation in order to close the system. This equation is the quasi-neutrality equation[START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF] [70, p. 40] [73] which results from the fact that locally, the plasma stays almost neutral electrically. The derivation of the quasi-neutrality equation in the TERESA model is described extensively in [70, p. 40] and only a brief summary will be given in this section. The quasi-neutrality equation results from the Maxwell-Gauss equation I.54. In the electrostatic case this equation leads to the Poisson equation (I.62) which we recall:

  e (E) fe √ EdE (II.28)

  Figure II.1: Representation of the Semi-Lagrangian method used to solve the Vlasov equation. Picture taken from [70, p. 80].

Fig II. 2 Figure II. 2 : 1 II

 221 Figure II.2: Artificial diffusion in the buffer zone toward the edges of the simulation box. The parameters are D 0 = 0.001, ψ 0 = 0.1, L D = 0.02 and L ψ = 1

  -3 , with ∆t ref = 1 × 10 -5 and the different ∆t ∈ [5 × 10 -5 ; 1 × 10 -4 ; 5 × 10 -4 ; 1 × 10 -3 ].

Figure II. 3 :

 3 Figure II.3: Relative error in ψ of a single test particle in an analytic electric potential at time t = 1×10 -3 for simulations with different time-step ∆t, compared to the reference simulation with a time-step ∆t ref = 1 × 10 -5 . ψ 0 is the position of the test particle from the reference simulation.

Fig. II. 3

 3 Fig. II.3 shows that the relative error in ψ grows like ∆t 4 . Thus the test with the analytic electric potential indicates a proper implementation of the RK4 algorithm in the TERESA code.

  Figure II.4: In log-log scale, relative error in α (II.4a) and in ψ (II.4b) at time t = 1 × 10 -2 for simulations with different time-step ∆t, compared to a reference simulation with a time-step ∆t ref = 5 × 10 -6 , where α 0 and ψ 0 are the positions of the test particle from the reference simulation.

Figure III. 1 :Table III. 2 :

 12 Figure III.1: Ion density (III.1a) and temperature (III.1b) profiles at t = 0, t = 6 and t = 12.

2 )Figure III. 2 :

 22 Figure III.2: Equilibrium ion distribution function in E (III.2a) and in V (III.2b) at ψ = 0.5.

Figure III. 3 :Figure III. 4 :

 34 Figure III.3: Time-evolution of 5 α-modes, along with the 0 th mode, in semi-log at ψ = 0.5.

Figure III. 6 :

 6 Figure III.5: Electric potential φ map at t = 6.

  Fig. III.4b shows the energy spectrum E m of the α-modes, with E m = m 3 φ 2 m , superposed with the slope in m -5/3 and m -3 , specific to cascades detailed in chapter I.4.

  Fig. III.6 shows the Test Particle Distribution Function (TPDF) at different times in ψ, for E = 0 and E = 1.

Figure III. 7 :

 7 Figure III.7: Trajectories as a function of time of 3 test particles with different energies, E = 0, E = 0.8 and E = 1.74 in (α, ψ) space. The black cross is the initial position of the test particle at t = 6. Trajectories for E = 0.8 and E = 1.74 are cut at some time to preserve readability of the figure.

Figure III. 8 :

 8 Figure III.8: Time-evolution of the test particles MSD for E = 0 (III.8a) and E = 0.8 (III.8b). In red is the linear fit of the MSD in phase 3 (diffusive phase).

Figure III. 9 :

 9 Figure III.9: Time-evolution of the test particles MSD for E = 1.74 (III.9a) and E = 2.71 (III.9b). In red is the linear fit of the MSD in phase 3 (diffusive phase).

3 Figure III. 10 :Figure III. 11 :

 31011 Figure III.10: Time-evolution of the test particles MSD for E = 3.2. In red is the linear fit of the MSD in phase 3 (diffusive phase).

Figure III. 12 :

 12 Figure III.12: Radial diffusion flux of the test particles (in blue) in velocity space, compared to the total flux (in red) given by the gyrokinetic simulation in red, averaged over ψ ∈ [0.4; 0.6] and over an auto-correlation time t ∈ [6; 7].

For E = 0 ,Figure III. 13 :

 013 Figure III.13: Time-evolution of 5 α-modes, along with the 0 th mode, in semi-log at ψ = 0.5.

Figure III. 14 :

 14 Figure III.14: Log-log spectrum of the α-modes taken at ψ = 0.5 and averaged over t ∈ [6, 7].

FigureFigure III. 16 :

 16 Figure III.15: (III.15a) Electric potential φ map at t = 6, and (III.15b) Time evolution of the electric potential φ as a function of α, for ψ = 0.5.

Fig

  Fig. III.20 shows the random walk diffusion coefficient in energy space calculated from the slopes of the MSD. Similarly to the previous case, the coefficient present a peak around E = 1.38 (E = 1.74 before), but is 5 times smaller than the coefficient showed on Fig. III.11. The difference in peak positions indicates that the particles resonate with the modes at a different energy than in the previous section. From this diffusion coefficient we estimate a radial diffusion flux which we compare to the total radial flux obtained from TERESA, Fig. III.21. Here we averaged Γ Total and ∂ f α ∂ψ from t = 6 to t = 10 and from ψ = 0.4 to ψ = 0.6 to smooth out the important variations of these quantities due to the streamer-like

Figure III. 18 :

 18 Figure III.18: (III.18a) Time-evolution of the test particle MSD for E = 1.38, with the estimated slope for the interval t ∈ [6.1 -7] in red, and (III.18b) TPDF as a function of ψ, for E = 1.38, at times t = 6, t = 6.02, t = 6.04, t = 6.06, t = 6.2, t = 7.5 and t = 10.

FigureFigure III. 20 :Figure III. 21 :

 2021 Figure III.19: (III.19a) Time-evolution of the test particle MSD for E = 2.02, with the estimated slope for the interval t ∈ [6.1 -10] in red, and (III.19b) TPDF as a function of ψ, for E = 2.02, at times t = 6, t = 6.02, t = 6.1 and t = 10.

Fig. (IV. 1 )

 1 Fig. (IV.1) shows the obtained time evolution of a selection of dominant modes.We observe a phase of linear growth of the plasma potential from t = 0 until t ≈ 4 -5. Linearly, the most unstable mode is the mode number m = 9, where φ = m φ m exp(imα). Its linear frequency and growth rate are ω m=9 = 15.6 and γ m=9 = 0.987. After t = 5 the modes achieve their saturation level. The saturation amplitude in terms of the root mean square φ rms is of the order of φ rms ∼ 0.05. In physical units, since we have chosen an inverse aspect ratio = 0.1, we obtain eφ rms /T ∼ 5 × 10 -3 . At t ≈ 5, the most intense modes are the modes around m = 10.For the rest of this chapter we investigate the time t = 5, at the beginning of the turbulent phase, where the dominant modes have just achieved saturation.

Figure IV. 1 :

 1 Figure IV.1: Time-evolution of the amplitude of 10 Fourier modes of the electric potential φ. The modes m = 5 to m = 13 correspond to the 9 most intense modes and the m = 0 mode corresponding to the zonal flow. The intensity of each m = 0 mode grows in the linear phase from t = 0 to t ≈ 5, when they reach saturation level. The dominant modes at time t ≈ 5 are the modes 7, 8, 9, 10 and 11. This plot is for ψ = 0.5. For ψ = 0.3 or ψ = 0.7, the only significant difference is that the m = 0 mode corresponding to the zonal flow has an intensity of the same order of magnitude as the other modes.

  Fig. IV.2a shows the density flux Γ(ψ, t), and Fig. IV.2b shows the heat flux q(ψ, t), both at the given time t = 5. The quasi-linear fluxes will be discussed later.

  and for the resonant particles E ∈ [1.4; 2.4].

Figure IV. 2 :

 2 Figure IV.2: Nonlinear and quasi-linear radial density (IV.2a) and heat (IV.2b) fluxes as a function of ψ and at t = 5. These fluxes are integrated over the whole velocity space. The fluxes are negative: oriented from ψ = 1 (the center) toward ψ = 0 (toward the edge).

  Fig. IV.4 and IV.5, discussed later. Figs. IV.3a and IV.3b show the contribution to the heat flux of the thermal and the resonant particles, respectively.

Figure IV. 3 : 4 E 1 =0Λ 4 E 2

 34142 Figure IV.3: (IV.3a) q N L,QL th

9 )

 9 We transform f and φ in Fourier space in the α direction, fm (ψ, E, t)

∂

  Figures IV.2 and IV.3 include the QL prediction for the density and heat fluxes. The NL fluxes and QL predictions are in qualitative agreement. However, there are significant quantitative discrepancies, even as we have chosen the best case scenario for QL theory (a small auto-correlation time, φ and ∂f ∂ψ directly from simulation and no mixing length estimate). The quasi-linear estimation of the density flux, the heat flux, the thermal particles contribution to the heat flux and the resonant particle contribution to the heat flux respectively present a 16%, 16%, 50% and a 20% discrepancy compared to the nonlinear flux obtained from the TERESA simulation.

Figure IV. 4 :

 4 Figure IV.4: Resonance condition for each modes m with Ω d = 1 and Z = 1. The time averaging around t = 5 of ∂ φ0 ∂ψ is necessary and is explained in section IV.2. This figure, coupled with Fig. IV.1, shows that we can expect the resonance to be the most intense for an energy E ≈ 1.78 to E ≈ 1.83 at t = 5, because the most intense modes in ψ = 0.5 are the modes between m = 7 and m = 11.

  IV.1).

Fig. IV. 4 Figure IV. 5 :

 45 Fig. IV.4 shows the resonant energy E R as a function of the modes m. The dominant modes around t = 5 and at ψ = 0.5 (m = 7 -11) have their resonant

Fig. IV. 5

 5 Fig. IV.[START_REF] Peters | The challenge to keep global warming below 2 c[END_REF] shows the particle flux (or flux of phase-space density f ), Λ NL (ψ, E, t), at a fixed radius ψ = 0.5. The sign is everywhere consistent with a flattening in ψ of f α (since ∂ f α ∂ψ is negative for E < 3/2 and positive for E > 3/2, as explained in subsection IV.2.2).As expected, the nonlinear (from TERESA simulation) fluxes in the energy (or velocity space), at the middle of the box and at t = 5, present a narrow resonance negative peak. The resonance peak width is about ∆E ≈ 0.15 (obtained by a Gaussian fit of the resonance peak, Λ NL ∼ exp[(E -E 0 )/∆E]). Strikingly, the resonance peak accounts for most of the flux. Indeed the flux at the resonance is one order of magnitude stronger than the flux for thermal particles, while the high energy particle flux tends to 0. We must keep in mind that the position of this resonance peak is determined by the most intense electrostatic modes at this ψ and time, and can be slightly shifted by the intensity of the zonal radial electric field ∂ φ0 ∂ψ .

  Fig. A.1 shows the error on E relatively to E 0 = E(t = 0), as a function of time. Fig A.1a is from a simulation with buffers at the edges. We observe a relative error on the E conservation of a few percents. Fig A.1b is the same simulation but without the buffers at the edges. It shows a relative error on the E conservation of less than 0.1%.

Figure A. 1 :

 1 Figure A.1: Relative error on the E conservation as a function of time, compared to E 0 = E(t = 0), with buffers a the edges (A.1a) and without buffers (A.1b). The simulation time-step is ∆t = 1.10 -3 .

  Fig. A.2a and A.2b show the relative error in E, at time t = 1 × 10 -2 , of different simulations with different time-step ∆t, with and without buffers, compared to the reference simulation of time-step ∆t ref = 5 × 10 -6 . We observe on Fig. A.2b that the relative error on E grows linearly with the time-step, and when the time-step is reduced of one order of magnitude, the relative error on E is reduced by two order of magnitude. The same test but with buffers on the edges, Fig. A.2a, shows that the relative error on E stays of the same order of magnitude regardless of the time-step ∆t.

Figure C. 1 :

 1 Figure C.1: (C.1a) Comparison between linear and nonlinear terms from the system equations. In the quasi-linear theory the nonlinear coupling is neglected. (C.1b) Zoom around the resonance on the key terms contributing to radial transport : L 3 and NL 3 .

  are recalled in table I.1.

		s =proton s =electron
	Cyclotron motion :
	ρ cs (m)	3.1 × 10 -3 7.3 × 10 -5
	ω cs (rad/s)	5.1 × 10 8	9.3 × 10 11
	T cs (s)	1.2 × 10 -8 6.8 × 10 -12
	f cs (Hz)	8.1 × 10 7	1.5 × 10 11
	Bounce motion :
	δ bs (m)	4.6 × 10 -2 1.1 × 10 -3
	ω bs (rad/s)	3.4 × 10 4	1.5 × 10 6
	T bs (s)	1.9 × 10 -4 4.2 × 10 -6
	f bs (Hz)	5.4 × 10 3	2.4 × 10 5
	Toroidal drift motion :
	ω ds (rad/s)	96	96
	T ds (s)	6.7 × 10 -2 6.7 × 10 -2
	f ds (Hz)	15	15
	Table I.1: Typical frequencies and lengths of the trapped particles for ITER pa-
	rameters : a particle's temperature T s = 150 × 10 6 K, a toroidal magnetic field
	B ϕ = 5.3 T, a poloidal magnetic field B θ one tenth of the toroidal one so that

.4.2 Action J 2s associated with the bounce (or banana) motion

  is the contour of the cyclotron motion, dl s = ρ cs dθ e θ and v cs = -ρ cs ω cs e θ , thus Cs e s A. dl s to J 2s is negligible. Indeed we can write, using Stokes theorem :

	The contribution from C	A. dl =		∇ × A. dS =	B. dS	(I.41)
	where dS = dS n is the surface enclosed by the banana contour C. n which is
	the vector normal to the surface S is almost perpendicular to B, thus Cs e s A. dl s ≈
	0. Thus we have	J 2s =	1 2π Cs	m s v s . dl s .		(I.42)
	I.2
								2.1.2. Considering the axisymmetric assump-
	tion, we use the cylindrical coordinate system where B = B z e z and thus A = Bzr 2 e θ
			J 1s =	1 2π Cs	(m s v cs + e s A). dl s	(I.37)
	where C s J 1s = -	1 2π	0	2π	m s ρ 2 cs ω cs dθ +	e s 2π	0	2π	ρ 2 cs B z 2	dθ	(I.38)
	Using the definition of µ s we obtain for a species s :
							J 1s = -	m s e s	µ s	(I.39)
	I.2The second action is associated with the bounce motion of trapped particles,
			J 2s =	1 2π Cs	(m s v s + e s A). dl s .	(I.40)

.4.3 Action J 3s associated with the toroidal kinetic momentum

  Rdϕ e ϕ and v ϕs is the velocity in the toroidal direction ϕ, which can be approximated by v ϕs ≈ v s ,

	The third action J 3s is associated with the toroidal precession of the particles,
							J 3s =	1 2π Cs	(m s v ϕs + e s A). dl s	(I.43)
	where dl = J 3s =	1 2π	0	2π	(m s v s + e s A ϕ )Rdϕ = m s v s R + e s RA ϕ	(I.44)
	Using Stoke's theorem,
	RA ϕ =	1 2π	2π ϕ=0	A ϕ Rdϕ =	1 2π	A. dl =	1 2π	∇ × A. dS =	1 2π	B. dS = ψ
											(I.45)
	thus,
									J 3s = m s v s R + e s ψ	(I.46)

  .2.

	Quantity	Value
	Ion Larmor radius	ρ i = 0.001
	Ion banana width	δ bi = 0.01
	Initial temperature gradient	κ T = 0.15
	Initial density gradient	κ n = 0.05
	Trapped particles precession frequency	Ω d = 1
	C ad	C ad = 0.1
	C pol	C pol = 0.1
	Electron dissipation [141]	δ m = 0.02
	Table III.3: Input parameters.

  4. Phase 4: a phase of saturation due to nonlocal effects and boundary conditions, where the test particles have explored the whole simulation box in ψ and the MSD reaches a plateau at MSD≈ 0.08. For E = 0, see Fig. III.8a, the trapped particles have no kinetic energy in the α direction. Phase 1 is from t = 6 to t ≈ 6.1, where the MSD grows quadratically with time until MSD≈ 10 -3 . Phase 2 appears from t ≈ 6.1 to t ≈ 6.3 where the MSD grows to MSD≈ 2.3 × 10 -3 . Phase 3 is from t ≈ 6.3 to the end of the simulation. Phase 4 does not appear on this figure, but the MSD would reach the saturation phase if the total simulation time were approximately one order of magnitude longer. E = 0.8, Fig. III.8b, is an intermediary case. Phase 1 is present from t = 6 to t ≈ 6.1, where the MSD grows quadratically with time until MSD≈ 10 -4 . Phase 2 appears from t ≈ 6.1 to t ≈ 6.7 where the MSD grows to MSD≈ 1.5 × 10 -2 . Phase 3 is from t ≈ 6.3 to the end of the simulation, and is where we fit linearly the MSD. Phase 4 is again not present, for the same reason as before. Then the MSD enters phase 4 as the test particles are subject to boundary effects, finally reaching a plateau at MSD≈ 0.08 at the end of the simulation. For E = 3.2, Fig. III.10, the test particles first explore the potential structure they where initialized in, in phase 1, from t = 6 to t ≈ 6.1. Then the MSD directly enters phase 3 as the test particles follow Brownian motion and the MSD grows linearly, until the end of the simulation. Phase 4 does not appear during the simulation time although the saturation would appear with a longer simulation (around t ≈ 20 -30). From the MSD at each E, we calculate the slope of the MSD in phase 3, and thus we can estimate a radial random walk diffusion coefficient of the test particles in E space (or velocity space), Fig. III.11.

For E = 1.74, see Fig III

.

9a, the test particles resonate with the TIM. Phase 1 is from t ≈ 6 to t ≈ 6.4. Phase 2 does not appear as the MSD transitions directly to phase 3. Phase 3 is from t ≈ 6.4 to t ≈ 6.8 where the MSD grows, linearly in time, from MSD≈ 0.02 to MSD≈ 0.06, as the test particles diffuse rapidly in ψ. Phase 4 appears from t ≈ 6.8 to the end of the simulation, where the test particles have explored the whole simulation box in ψ and the MSD reaches a plateau at MSD≈ 0.08. At E = 2.7, Fig. III.9b, the test particles are above the resonance energy and have a higher α than the precedent cases. Phase 1 is from t ≈ 6 to t ≈ 6.1. Phase 2 does not appear as the MSD directly enters phase 3 from t ≈ 6.1 to t ≈ 8, with the MSD growing linearly from MSD≈ 10 -3 to MSD≈ 0.04.

Table III .

 III 4: Grid used for our simulation. α and ψ are the phase-space variables while E (or V ) is a parameter.

	Quantity	Value
	Ion Larmor radius	ρ i = 0.01
	Ion banana width	δ bi = 0.1
	Electron Larmor radius	ρ e = 0.001
	Electron banana width	δ be = 0.01
	Initial temperature gradient	κ T i = κ T e = 0.15
	Initial density gradient	κ ni = κ ne = 0.05
	Trapped particles precession frequency

  E = 1.38, Fig. III.18a, the test particle MSD again presents a first rapid growth (from t = 6 to t ≈ 6.1) to a MSD ∼ 0.04. Then the MSD grows linearly in time a first time from t ≈ 6.1 to t ≈ 7.5 and a second time, slower, from t ≈ 7.5 to the end of the simulation. The slope of the second linear growth is 2.51 × 10 -3

Table IV .

 IV 1: Grid used for our simulations. α and ψ are the phase-space variables while E (or V ) is a parameter.

	Quantity	Value
	Ion Larmor radius	ρ i = 0.01
	Ion banana width	δ bi = 0.1
	Initial temperature gradient	κ T = 0.15
	Initial density gradient	κ n = 0
	Trapped particles precession frequency

  Since the RHS phase-mixes for t -s larger than a typical growth time γ -1 and a typical timescale of relaxation of f α , we can approximate s t γ(t )dt by γ(t)(s -t) and ∂ ψ f α | s by ∂ ψ f α | t . With this approximation, the time evolution of the α-averaged distribution function simplifies to

	R,m -ωm)(s-t) ∂ f α ∂ψ s	exp	t	s	γ(t )dt ds
					(IV.19)

  • C requiert une neutralité carbone mondiale en 2050[START_REF] Tollefson | Ipcc says limiting global warming to 1.5 c will require drastic action[END_REF][START_REF] Rogelj | Energy system transformations for limiting end-of-century warming to below 1.5 c[END_REF]. En 2015, 42% des émissions de CO 2 dans le monde provenait de la combustion d'énergies fossiles à des fins de génération d'électricité et de chaleur, tandis que le transport comptait pour 24%[START_REF]Co2 emissions from fuel combustion highlights[END_REF] des émissions. Bien que les énergies renouvelables soient en développement, ce sont des sources d'énergie intermittentes qui doivent encore être associées à d'autres sources d'énergies[21, p. 13] [START_REF] Steinke | Grid vs. storage in a 100% renewable europe[END_REF][START_REF] Gowrisankaran | Intermittency and the value of renewable energy[END_REF]. La fission nucléaire produit de l'énergie en cassant un noyau atomique lourd comme celui de l'uranium et donnant des éléments plus légers. La fission nucléaire ne produit pas directement de CO 2 mais génère des déchets radioactifs avec une durée de vie de plusieurs centaines de milliers d'années. Ces déchets doivent être stockés dans des bunkers souterrains. En France, où 75% de l'électricité est générée par la fission, 1.54 × 10 6 m 3 de déchets radioactifs sont stockés en 2016, en augmentation de +8.5 × 10 4 m 3 depuis 2013[24, p. 14].Dans ce contexte, il existe une autre réaction nucléaire semblant très prometteuse pour l'avenir qui est la fusion nucléaire. Au lieu de casser un noyau atomique lourd en éléments légers, la fusion nucléaire combine deux noyaux atomiques légers en un plus lourd afin de produire de l'énergie. La problématique de la fusion nucléaire est totalement différente de celle de la fission, et est bien plus difficile à réaliser sur Terre. La fusion nucléaire a lieu 115 dans la Nature depuis des milliards d'années, au coeur des étoiles. Au centre du soleil, des noyaux d'hydrogène, ou protons, se combinent afin de créer des noyaux d'hélium, composés de deux protons et de deux neutrons. Plus de 80 types de réactions de fusion entre éléments légers ont été observés et étudiés[START_REF] Rax | Physique des tokamaks[END_REF]. La plus accessible sur Terre est la fusion du deutérium et du tritium (D-T), deux isotopes de l'hydrogène, parce qu'ils présentent le plus grand ratio section efficace/énergie (voir Fig.1). Le deutérium est composé d'un proton et d'un neutron tandis que la tritium est composé d'un proton et des deux neutrons. Cette réaction de fusion produit un noyau d'hélium (une particule α) et un neutron (Eq. 1). Cette réaction ne produit pas de déchet radioactif à longue durée de vie comme les réactions de fission. De plus, cette réaction ne présente pas de risque "d'emballement" ou de "fusion du coeur". Le deutérium est abondant sur Terre et est présent dans l'eau des océans : dans chaque mètre cube d'eau de mer, on trouve 33 grammes de deutérium. Au rythme actuel de consommation d'énergie, il y aurait assez de deutérium sur Terre pour 10 11 ans. Le tritium est radioactif et a une demi-vie de ∼ 12 ans, il n'est donc pas présent dans la nature et doit être produit. Par chance, il est possible de produire du tritium de façon élégante via les réactions [25, p. 26] : Eq. 2 et Eq. 3 qui utilisent le neutron de la réaction (Eq. 1), combiné avec du Lithium, qui est disponible sur Terre[25, p. 26]. La fusion nucléaire est la promesse d'une énergie quasi-illimitée, en grande quantité, sans rejet de gaz à effet de serre, sans déchet radioactif à longue durée de vie et sans possibilité d'emballement.Dans ce cas pourquoi la fusion nucléaire n'est toujours pas une source d'énergie commerciale sur Terre ?Afin de fusioner, le deutérium et le tritium doivent vaincre la barrière de potentiel provenant du fait que les deux noyaux sont chargés positivement (répulsion Coulombienne). Afin de réaliser cela, les particules chargées positivement doivent atteindre une vitesse très élevée ou de façon équivalente, une très grande énergie. La Fig.1montre la section efficace de différentes réactions de fusion. On peut voir que pour les énergies < 200 keV, la réaction D-T a la plus grande section efficace, ainsi, la réaction D-T est la plus simple à effectuer. La réaction de fusion du D-T a un maximum de section efficace à 100 keV [26, p. 2] puis retombe pour les énergies plus élevées. Malgré cela, la température optimale se situe plutôt autour de 10 -20 keV [21, p. 51] [26, p. 2]. Cependant, 10 keV sont équivalent à plus de 100 millions de Kelvin, ce qui est près de 10 fois la température au coeur du Soleil.A ces températures, la matière se trouve sous forme de plasma. L'état plasma est considéré comme le quatrième état de la matière. Il peut être compris de façon simple avec le raisonnement suivant : chauffer un solide le transforme en liquide, chauffer un liquide le transforme en gaz et chauffer un gaz le transforme en plasma. Plus précisement, un plasma est un gaz qui est partiellement ou totalement ionisé. On considère que 99% de la matière baryonique de l'univers se trouve sous forme de plasma. A des énergies > 10 keV, le deutérium et le tritium sont totalement ionisés. Etant donné qu'un plasma est constitué de particules chargées (ions positifs ou électrons), il peut être confiné dans une chambre à vide via de puissants champs magnétiques. Cela est très pratique et mène au concept de Fusion par Confinement Magnétique (FCM). La première machine à confinement magnétique que l'on pourrait imaginer serait basé sur un champ magnétique bouclé sur lui même, avec le plasma en forme de tore, de façon à ce que les particules chargées suivent le champ magnétique bouclé sur lui même et soient confinées. Cela ne marche pas, puisqu'il existe des dérives à travers les lignes de champs, ce qui fait que le plasma n'est pas confiné. Afin de pallier cela, les lignes de champ magnétique doivent être enroulées aussi sur elles-mêmes à la façon d'un ressort. Une classe de machine exploitant cette idée est le tokamak, un acronyme russe de "chambre toroïdale avec bobines magnétiques". Un tokamak est une machine complexe qui présente une chambre à vide de forme torique où le plasma est présent, et des champs magnétiques de confinement de 1 -10 Tesla. Afin d'atteindre des énergies > 10 keV, le plasma doit être chauffé, ce qui se fait de façon routinière via des méthodes comme le chauffage ohmique, l'injection de neutres ou l'absorption résonante d'ondes électromagnétiques radiofréquences[26, p. 239].Pour que cette approche soit digne d'intérêt, la puissance produite par le réacteur à fusion doit être supérieure à la puissance utilisée pour chauffer le plasma. Cela mène au critère de Lawson[START_REF] Lawson | Some criteria for a power producing thermonuclear reactor[END_REF] [26, p. 11], nT τ E > 10 21 m -3 .keV.s (1) où n est la densité de particule, T la température du plasma en unité d'énergie (10-20 keV) et τ E est le temps de confinement de l'énergie, qui est le temps que mettrait l'énergie du plasma à se dissiper si l'on coupait tous les systèmes de chauffage extérieurs. Dans le produit triple, la densité est limitée dans les tokamaks à un maximum d'environ 10 20 particules par mètre cube à cause de la limite de Greenwald [21, p. 500]. La température est limitée par la limite de β où β est le ratio entre la pression cinétique P k = nk B T et la pression magnétique B 2 /2µ 0 . Ce ratio dans les tokamaks est autour de ∼ 1% et une augmentation de β peut mener à une réduction du confinement. Ainsi, dans les tokamaks actuels, seul le temps de confinement de l'énergie τ E peut être augmenté. Cependant, de nombreux mécanismes tendent à réduire τ E . Un des principaux est le transport de particules et de chaleur depuis le centre du plasma (le coeur) vers les parois du tokamak (les bords), à travers les lignes de champ magnétique. Comprendre le transport afin de le contrôler et d'augmenter τ E est un sujet de la plus haute importance dans le domaine de la fusion par confinement magnétique. Le transport dit "anormal" est maintenant reconnu comme étant la principale composante du transport. Le transport anormal est dominé par du transport turbulent. La turbulence est un processus hautement non linéaire, multi-échelles en temps et en espace [32-34]. La compréhension de la turbulence dans les plasmas de tokamak a fait d'énormes progrès dans les dernières décennies grâce à une abondance de résultats analytiques, expérimentaux et numériques. Le transport anormal prend son origine dans des micro instabilités (microscopiques comparées à la taille de la machine) dans le plasma, ce qui génère du transport turbulent [21, 26, 35-39]. Elles créent des perturbations du potentiel électrique φ dans le plasma et mènent à des champs électriques E. Lorsqu'un champ électrique et un champ magnétique sont présents, les particules chargées subissent un mouvement de dérive appelé E × B, et ce qui mène à l'apparition d'ondes de dérive

processes, which is coherent with the observed streamer-like turbulence.

Acknowledgments

Acknowledgments potentiels et les phénomènes balistiques tels que les avalanches [54,62,83,84]. Une autre partie de ce travail était focalisé sur l'étude du transport grâce à ces particules test, ainsi qu'étudier les flux de particules et de chaleur avec une haute résolution dans l'espace des énergies, et de comparer ces flux aux estimations de la théorie quasi-linéaire.

Dans le chapitre I, nous introduisons la configuration magnétique d'un tokamak, et les trajectoires de particules chargées dans ce champ magnétique. Nous décrivons ensuite le modèle puis le code TERESA dans la section II.1 et nous détaillons l'implémentation des particules test dans le code TERESA dans la section II.2. Dans le chapitre III nous utilisons ce nouveau module de particules test pour séparer les contributions advectives des contributions diffusives aux flux de particules. Pour cela nous utilisons des millions de particules test et les laissons évoluer dans une simulation de plasma turbulent. En étudiant leur évolution temporelle, on calcule un coefficient de marche aléatoire dans l'espace des énergies. On peut ainsi estimer l'importance de la diffusion dans le flux total qui prend en compte tous les processus de transport. Dans le chapitre IV, nous étudions les flux de particules et de chaleur avec une haute résolution dans la dimension en énergie. Ces flux sont dominés par un fin pic de résonance. On compare ensuite ces résultats à la théorie quasi-linéaire. Nous trouvons que les prédictions quasi-linéaires sont en accord qualitatif avec les simulations non linéaires de TERESA, mais avec près de ∼ 50% d'écart dans l'amplitude du pic. Nous faisons ensuite une étude terme par terme et nous trouvons que dans notre cas, des termes négligés par la théorie quasi-linéaire sont en fait du même ordre de grandeur que les autres termes non négligés.

Appendix A

Order of accuracy of the numerical scheme

In the TERESA code, the Vlasov-Quasi-neutrality system of equations is solved with a second order accuracy in time.

We aims to verify this, first analytically and then via simulation results. First let us find the analytical expression of the local error n+1 φ on the electric potential generated after each time-step of the simulation. In the TERESA code, the electric potential φ is solved using the distribution function f via a matrix relationship, thus we can link n+1 φ to n+1 f , the local error on the distribution function, as:

where ϑ is a constant. The analytical solution of the distribution function after one time-step ∆t is written :

where u and v represent the advection in, respectively, the α and ψ directions. The distribution function is calculated at each time-step as:

where F n is interpolating f n i,j in the (α, ψ) grid, using a 2D spline. Because we focus on the error in time, we can make the assumption, which is later verified, that the mesh (α, ψ) does not affect the error in time, which is equivalent to take the limit (∆α → 0, ∆ψ → 0), where ∆α and ∆ψ are the grid space. Thus we can write

Appendix A. Order of accuracy of the numerical scheme

The local error on the distribution function f generated from the numerical scheme is written:

Appendix C

QL approximation terms comparison

Quasi-linear theory neglects some of the nonlinear terms in the Vlasov-Poisson system.

In this appendix we give details about the neglected or kept linear or nonlinear terms in the quasi-linear approximation.

The linear terms from Equation IV.12 are:

and the nonlinear terms are:

In order to compare their magnitude, we take their absolute value and sum them over all non zonal modes, which gives us the terms L 1 , L 2 , L 3 , NL 1 , NL 2 and NL 3 , where

The terms kept in the quasi-linear equations are the terms L 1,2,3 and NL 1 . We plot on Figure C.1 the linear and nonlinear terms of the system at a time t = 5

111

Abstract

One of the most promising approach to controlled nuclear fusion is the tokamak. It is a toroidal machine confining a fusion plasma using magnetic fields. Transport of particles and heat, from the core toward the edges happens spontaneously, degrades the efficiency of the tokamak, and is driven by turbulence. We use a bounce-averaged 4D gyrokinetic code which solves the Vlasov-Quasi-neutrality system. The code is based on a reduced model which averages out the cyclotron and the bounce motion of the trapped particles to reduce the dimensionality. In this work we developed and tested a new module for the code, allowing to track test particle trajectories in phase space. As a first result obtained with test particles, we achieved to separate the diffusive contribution to the radial particle flux in energy space, from the non-diffusive contributions. Both fluxes present an intense peak indicating resonant particles dominate transport. On short period of time the test particles undergo a small scale advection, but on longer times, they follow a random walk process. We then explored with greater accuracy the fluxes in energy space. Furthermore we compared the obtained fluxes with quasi-linear predictions and found a qualitative agreement, although there was a ∼ 50% discrepancy in the peak magnitude.

Keywords: Transport, Turbulence, Trapped particles, Simulations

Résumé

Une des approches les plus prometteuses pour réaliser la fusion nucléaire est le tokamak. Un tokamak est une machine toroïdale qui confine le plasma de fusion via des champs magnétiques. Le transport de particules et de chaleur du centre vers les bords apparaît spontanément, dégrade le confinement et est produit par la turbulence. Nous utilisons un code 4D gyrocinétique qui résout le système Vlasov-Quasi-neutralité. Le code est basé sur un modèle réduit qui moyenne le mouvement cyclotron et de rebond des particules piégées pour réduire la dimensionnalité. Dans ce travail, nous avons développé et testé un nouveau module pour le code permettant de suivre des particules test dans l'espace des phases. En premier résultat, grâce aux particules test nous avons pu séparer la contribution diffusive au flux de particules, des contributions nondiffusives. Sur des temps longs, elles suivent un processus de marche aléatoire. Les deux flux présentent un important pic indiquant que les particules résonantes dominent le transport. Nous avons ensuite exploré avec une plus grande précision les flux dans l'espace des énergies. Enfin nous avons comparé les flux obtenus aux prédictions quasi-linéaires et avons trouvé un accord qualitatif, avec cependant un écart de ∼ 50% dans l'intensité du pic.

Mots-clés:

Transport, Turbulence, Particules piégées, Simulations