Résumé

Le contrôle du transport turbulent est l'une des clés pour l'amélioration du temps de confinement nécessaire à la réalisation de la fusion thermonucléaire contrôlée. La description de la turbulence cinétique du plasma est un problème à 3 coordonnées spatiales et 3 coordonnées en vitesse.

En se basant sur le fort champ magnétique de confinement, l'approche gyrocinétique, largement utilisée, consiste à moyenner le mouvement cyclotron, qui est beaucoup plus rapide que le phénomène de turbulence. Une telle réduction permet de simplifier le problème à trois coordonnées spatiales, relatives aux centres-guides des particules, une vitesse parallèle ou énergie et une vitesse perpendiculaire apparaissant comme l'invariant adiabatique. La description gyrocinétique non linéaire requiert des simulations numériques de haute performance massivement parallèles. Toute la difficulté est dûe aux termes non linéaires (crochets de Poisson) qui décrivent les interactions multi-échelles, ce qui constitue un défi tant pour la théorie que pour la simulation. Toute approche réduite, basée sur des hypothèses bien contrôlées, est donc intéressante à développer.

Sur la base de cette ambition, cette thèse concerne la turbulence des particules piégées dans le plasma magnétisé : ces particules disposent en effet d'une énergie cinétique parallèle insuffisante pour décrire toute une surface magnétique et décrivent un mouvement de rebond côté faible champ magnétique. En considérant des échelles de temps de la turbulence grandes à la fois devant la fréquence cyclotron et la fréquence de rebond, une double moyenne de la fonction de distribution des particules sur les mouvements cyclotron et de rebond est utilisée, et un système 4D est obtenu qui peut être considéré comme une forme réduite de la théorie gyrocinétique standard. Nous l'avons appelé "bounce averaged gyrokinetics" pendant ce travail. Même si cette description est grandement réduite par rapport à la théorie gyrocinétique, sa simulation numérique directe reste un challenge.

Une description des termes non linéaires en coordonnées polaires est choisie, avec une grille logarithmique en norme du vecteur d'onde, tandis que les angles sont discrétisés sur une grille régulière. L'utilisation d'une grille logarithmique permet de prendre en compte une large gamme de vecteurs d'ondes, donc la physique à très petite échelle. De manière analogue aux modèles en couches en turbulence fluide et afin de simplifier le système, seules les interactions entre couches voisines sont considérées.

iii Dans un premier temps, l'étude du système linéaire est présentée, en particulier les dépendances paramétriques des seuil linéaire d'instabilité et des solutions de l'équation de dispersion permettent de retrouver la forte anisotropie des taux de croissance des modes d'ions piégés (ou TIM) et des modes d'électrons piégés (ou TEM). Ces études permettent également de valider les codes numériques non-linéaires vis-à-vis d'un solveur aux valeurs propres développé indépendamment.

Dans un second temps, l'hypothèse d'isotropie dans l'espace de Fourier est utilisée. Ainsi il n'y a pas d'information de phase exacte pour de tels modèles en couches 1D, ce qui laisse un paramètre libre dans les coefficients d'interaction et un choix de phase arbitraire pour les inconnues permettant d'utiliser deux modèles : Gledzer-Ohkitani-Yamada (GOY) et Sabra. Les spectres en nombre d'onde révèlent des lois de puissance originales, approximativement ∝ k -4 pour l'énergie potentielle électrostatique (E φ ), et ∝ k -1 pour l'énergie libre (E f ). Ces exposants ne s'avèrent pas affectés par la valeur du paramètre libre, mesurant l'intensité des effets non-linéaires relativement aux termes linéaires. La comparaison entre les modèles GOY et Sabra montre de grandes différences quant au niveau de saturation et au comportement temporel de la phase turbulente, en particulier des oscillations à basse fréquence sont obersvées à l'aide du modèle GOY, mais ne sont pas retrouvées dans les simulations basées sur un modèle Sabra.

L'information de phase apparait donc comme très importante. De plus, puisque le taux de croissance d'instabilité est anisotrope, la simulation du modèle anisotrope doit être réalisée dans un troisième temps. Le système résolu numériquement est réduit à une espèce cinétique, en supposant que les autres espèces sont adiabatiques. Deux systèmes différents peuvent ainsi être étudiés : ions cinétiques et électrons adiabatiques (IC) d'une part, ou électrons cinétiques et ions adiabatiques (EC) d'autre part. Les spectres en nombre d'onde de l'énergie potentielle électrostatique et de l'entropie montrent des exposants similaires dans les deux cas : ∝ k -5 pour l'énergie potentielle électrostatique et ∝ k -7/3 pour l'énergie libre. Par ailleurs, dans le cas d'électrons cinétiques une dynamique de type prédateur-proie entre le flux zonal et la turbulence est observée, où le flux zonal présente une oscillation forte et régulière. Il faut noter que le niveau de flux zonaux dans le système avec ions cinétiques est significativement plus elevé que dans le cas d'électrons cinétiques.

Enfin, un système complètement cinétique (CC) est étudié. La comparaison avec les deux systèmes (IC) et (EC) montre des différences significatives iv dans les pentes des spectres, en particulier le spectre en énergie potentielle électrostatique du système complètement cinétique (CC) montre une pente en ∝ k -4 , différente de celles en ∝ k -5 observées quand une des deux espèces est considérée adiabatique (IC) ou (EC). Le système (CC) s'avère suivre le cas d'électrons cinétiques avec ions adiabatiques (EC) en ce qui concerne le mécanisme de saturation : le maximum du spectre est observé glisser progressivement du maximum de taux de croissance linéaire (avec k ≈ 20) vers des nombres d'ondes plus faibles (typiquement k ≈ 3). En même temps, l'anisotropie du spectre observée dans le cas (IC) est retrouvée dans le cas (CC), avec un rôle important joué par les streamers (modes correspondant à des tourbillons allongés radialement). Ceci interroge sur la pertinence de l'approximation d'une réponse adiabatique que ce soit pour les ions ou, plus communément effectuée, pour les électrons.

Introduction: turbulence and transport in Tokamaks

Energy is the most important issue in modern society. As a consequence of the development of civilization, the worldwide demand for energy has been increasing rapidly. In most of the countries around the world, energy is produced mainly by fossil fuels, even though many new sources and techniques exist. As a resource that can not be reproduced, the fossil fuels stored in the earth is fixed, which will certainly reach the end one day in the future, especially when the growth of world population is considered. Another big problem that is more urgent, is the pollution and the greenhouse effect, caused by the combustion of fossil fuels, which endanger the health of everyone, including the planet earth. So the research of new energy sources, especially those that are clean and abundant, is becoming more and more important.

Nuclear fusion 1.Nuclear fusion reaction

One of the possibilities, inspired by the reaction in the stars, is the transformation of fusion energy to electric energy, which has been envisaged at least from 50s. The method is to confine the light nuclides such as deuterium (D), tritium (T ), helium-3 (He 3 ), and lithium (Li ) as shown in figure (1.1), by a strong magnetic field in a reactor like Tokamak such that the nuclides can produce energy in the process of fusion reactions: Li 6 + n → T + α + 4.8Mev

D + T → α (3.
Li 7 + n → T + α + n -2.5Mev , (1.1) Figure 1.1: Binding energy per nucleon [START_REF] Azevedo | Selection of fuel cladding material for nuclear fission reactors[END_REF].

Since the energy is produced by nuclear reactions, the greenhouse effect that is linked specially to fossil fuels, is not an issue. Deuterium exists abundantly in sea water with an amount that could supply the energy consumption for more than millions of years, while tritium can be "breeded" using lithium, another abundant element. Furthermore fusion reactions do not leave long lasting radioactive products, and the problem of nuclear waste is much less serious than that for fission reactors.

The energy produced in this reaction is so high that the power of the designed fusion reactor, for example 500 Mw for ITER (International Thermonuclear Experimental Reactor), is usually very high comparing to that of the other power station in the world. If one day the nuclear fusion becomes a reality, energy will not be a limit for human beings for a long time.

Lawson criterion

The binding energy per nucleon is relatively small in very light or very heavy nuclides, as shown in figure (1.1) [START_REF] Azevedo | Selection of fuel cladding material for nuclear fission reactors[END_REF] [START_REF] Miyamoto | Plasma physics and controlled nuclear fusion[END_REF], while heavy nuclides provide the fuel for fission reactors, the lighter ones (which have especially low binding energy, such as H , D and T ) provide the potential fuel for fusion. The difficulty in extracting this energy is that, in order to achieve fusion reaction, one has to supply enough energy to heat the plasma to temperatures, which will allow positively charged particles to overcome their repulsive Coulomb potentials with a sufficient rate so that fusion reactions can occur. For example the reaction rate of one deuterium due to a large number of tritium atoms is given by:

d N d t
= n T 〈σv〉 , Figure 1.2: Fusion cross section as a function of the kinetic energy for different reactions [START_REF] Azevedo | Selection of fuel cladding material for nuclear fission reactors[END_REF].

where N is the number of reaction and d N d t is the reaction rate, n T is the number of tritium per volume, σ is the cross section and v is the relative velocity between the particles, the average 〈•〉 is performed on the velocity distribution of the particles. There are many possible fusion reactions, however D -T reaction has the largest cross section as shown in figure (1.2) [START_REF] Arnold | Cross Sections for the Reactions D(d , p)T, D(d , n)He 3 , T(d , n)He 4 , and He 3 ([END_REF] and the released energy is also the largest among the reactions given in equation (1.1), which is the main reason to choose the deuterium and tritium as the fuel in modern fusion reactor like the ITER machine.

In the case of many particles, one can obtain the reaction rate as follows:

d N d t = n D n T 〈σv〉 ,
thus the total fusion power P fusion should be :

P fusion = d E t ot d t = d (N E fusion ) d t = n D n T 〈σv〉E fusion ,
where E t ot is the total fusion energy, which is the sum of the energy E fusion produced in each fusion reaction. In order to understand the power budget, the fusion power should be compared to the power loss from the plasma, for which one can introduce the confinement time τ E , such that:

dW t h d t = - W t h τ E + P l ,
where P l is the loss of power, which can be seen as the sum of the radiation power of α particles defined as P a = P f us E α E α +E n and the external power P ext supplied to heat the plasma. W t h is the thermal energy of plasma, which can be written as:

W t h ≈ 3 2 n D T D + 3 2 n T T T + 3
2 n e T e V ≈ 3n e T e V , with the assumption of a D -T mixture such that n D = n T = n e 2 . In stationary regime, the confinement time can be defined as the ratio of the plasma thermal energy to the loss of power [START_REF] Miyamoto | Plasma physics and controlled nuclear fusion[END_REF]:

τ E = W t h P l .
The confinement time can be understood as the period τ E for which the plasma is confined.

Another key parameter of a fusion reactor is the amplification factor of the plasma Q , which is defined as:

Q = P f us P ext .
If Q = ∞, the reaction is self sustained by the fusion power, which is called ignition. If Q = 1, it is called "break-even", since one gets as much power from the device as one puts in.

An important and general criterion to measure the system to be able to produce energy, is that the generated fusion power should be able to maintain the temperature of the plasma against the loss of power without external power input such that the temperature is high enough to heat the plasma and to continue the fusion reaction [START_REF] Drouot | Study of turbulence associated with trapped particles in fusion plasmas[END_REF]:

P f us > P l ,
which, after a detailed calculation, gives a minimum required value for the product of the plasma (electron) density n e , the temperature of the plasma T (10 ∼ 20 kev) and the energy confinement time τ E , called the Lawson criterion [START_REF] Lawson | Some Criteria for a Power Producing Thermonuclear Reactor[END_REF]. For the system of deuterium and tritium, it should satisfy: n e T τ E > 3 * 10 21 kev.m -3 .s . This relation gives two different ways of fusion: the inertial fusion, which works with a very short confinement time (τ E ∼ 10 -11 s) and a very high density (n ∼ 10 31 m -3 ) and the magnetic confinement fusion, where the density is reduced to n ∼ 10 20 m -3 and the confinement time is τ E ∼ 1s.

In the magnetic confinement fusion, the temperature of the system should be high enough (~15 kev ) to continue the reaction. The density of the plasma n e is limited by the parameter β = p ki n p mag of the plasma with p ki n , the kinetic pressure of the plasma defined as p ki n = n e kT and p mag , the magnetic pressure of the plasma defined as p mag = B 2 2µ 0 , where B is the magnetic field and µ is the magnetic moment. The plasma β parameter must be smaller than 1 such that the magnetic pressure is higher than the kinetic pressure to avoid the disruption of the plasma.

Due to these conditions, the choice of the plasma density n e and the temperature T is limited, to satisfy Lawson criterion, it is necessary to make the energy confinement time τ E as long as possible. However the confinement of the particle and the heat in the Tokamak is not stable, due to the diffusive transport of the density, as well as the heat from the center to the boundary of machine, which can be formulated as:

∂ t n = D△n , ∂ t T = χ△T ,
with n and T , the density and the temperature, D and χ the diffusion coefficients of density and heat, respectively. This diffusive transport is what limits the confinement and strongly depends on the level of turbulence energy. So the control of turbulent transport is one of the key requirements for the improvement of the energy confinement time to make fusion energy a reality.

Kinetic description of plasma turbulence

In order to control the turbulent transport in a Tokamak, it is necessary to understand the mechanism of the turbulence and transport of the magnetized plasma. The understanding of the turbulence and transport in a Tokamak is mainly via the description of the temporal behavior of the plasma. A plasma consist of many ions and electrons, but the individual behavior of each particle can hardly be observed. What can be observed instead are statistical average, thus it is necessary to define a distribution function in the phase space to describe the evolution of a plasma. Note that if the density of a plasma is very low, it is seen as a collection of individual particles. If the density as well as the collision rate are high enough, it is seen as a fluid, between these two cases, it is seen as a kinetic plasma. The distribution function of a kinetic plasma can be introduced as f (r, v, t ), where r and v present the particle positions and the velocities, respectively.

Vlasov equation

The evolution of the particle distribution function f (r, v, t ) is described by the Vlasov equation, that is:

∂ f ∂t -H , f r,v = 0 , (1.2)
with H , the Hamiltonian of the particles, which can be generally defined as H ≡ 1 2m s Pq s A 2 + q s Φ, with P , the canonical momentum; A, the magnetic vector potential; Φ, the electric potential; m s and q s , the mass and the charge of the particle, respectively. H , f is the Poisson bracket defined as

H , f = d H d r d f d v - d f d r d H d v .
The right hand side of this formula is 0 for a collisionless plasma. Note that the particle distribution function f of a kinetic plasma depends on 3 spatial coordinates and 3 kinetic coordinates (or velocity). So it is a 6D system.

Quasi-neutrality equation

Vlasov equation is coupled to the electromagnetic fields, described by Maxwell's equations. When magnetic fluctuations are small as that observed in a tokamak and for scales that are larger than Debye length λ D , the system can be seen as quasi-neutral, which means the fluctuation of ion density δn i should be equal to the fluctuation of electron density δn e :

δn i = δn e , (1.3)
where the density δn s , can be calculated by the integration of the particle distribution function f s over the kinetic coordinates, where the index s = i , e represents the species of interest.

Difficulty of nonlinear simulation

A 6D problem is a challenge for theoretical analysis, so the understanding of this high dimensional system is based mainly on numerical simulations. Considering the capability of today's computers, simplification is useful for this system. One possible technique is to average out the cyclotron motion, since the cyclotron frequency ω c of the charged particle in a strong magnetic field is much higher than the typical frequency of plasma turbulence. This widely used technique, the so called gyro average [START_REF] Catto | Linearized gyro-kinetics[END_REF] in the field of gyrokinetic simulations, can reduce the 3D velocity coordinate to 2D, since the phase information of the cyclotron motion is averaged out:

v || , v ⊥ , θ =⇒ v || , v ⊥ ,
where v || is the velocity that is parallel to the magnetic field, v ⊥ is the perpendicular velocity and θ is the phase of the cyclotron motion. After the gyro average (figure (1.3)), the 6D kinetic system is reduced to 5D (typically three spatial coordinates of the guiding-centers of particles, one Figure 1.3: The 6D phase space is reduced to 5D by gyro average, where the phase of cyclotron motion is lost [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF].

parallel velocity or energy and one perpendicular velocity appearing as the adiabatic invariant, the cyclotron phase is lost in gyro average).

Nonlinear gyrokinetic [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF][START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF] simulations are usually massively parallel high performance numerical simulations. Note that being a 5D problem, the computation of nonlinear dynamics in gyrokinetic description that includes self-consistent multi-scale interactions usually requires millions of CPU hours and makes anything beyond medium nonlinear runs dominated by a single type of instability over a small range of scales, impractical. This, coupled with the complexity of numerical implementation can also make it harder to understand and isolate important physical mechanisms. In this regard, reduced models appear as intermediate tools, that are useful for isolating important physical mechanisms, to provide guidance to large scale gyrokinetic simulation efforts as well as comparison with experiments.

Reduced model

Caused by the complexity and the difficulty of the nonlinear gyrokinetic simulation, reduced models are potentially very promising.

Trapped particle model

Gyro average allows to reduce the 6D kinetic system to a 5D gyrokinetic system. Considering the motion of a charged particle, further simplification is possible. In a magnetic configuration of a Tokamak, where the magnetic field is stronger on the inner side than on the outer side, a charged particle can be trapped or passing depending on its kinetic energy due to the mirror effect. So in addition to the cyclotron motion, it displays also the bounce motion and the toroidal precession motion. Since the cyclotron frequency ω c and the bounce frequency ω b are much higher than the toroidal precession motion ω p , it allows to average out the cyclotron motion as well as the bounce motion as shown in figure (1.4):

This double average (gyro average plus bounce average) is called gyro bounce average [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF][START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF], which is usually used to simplify the discussion of the trapped particle (i.e, Trapped Electron Mode (TEM) and Trapped Ion Mode (TIM)) [START_REF] Kadomtsev | Trapped particles in toroidal magnetic systems[END_REF][START_REF] Cohen | Non-linear saturation of the dissipative trapped-ion mode by mode coupling[END_REF][START_REF] Tang | Residual trapped-ion instabilities in tokamaks[END_REF] turbulence and results in a 4D system [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] (i.e. two spatial coordinates: the toroidal angle α and the poloidal flux ψ playing the role of radial coordinate r , and two kinetic parameters: particle energy E and an adiabatic invariant µ appearing as the trappeed particle ratio f t ):

f ⇒ f α, ψ (r ) , E , µ .
Note that the gyro average reduces the 3D kinetic coordinates to 2D and the bounce average reduces, in some sense, the 3D spatial coordinates to 2D. So this is a 2D kinetic system, which is used to describe the TEM and TIM turbulence in the toroidal plane of a Tokamak. This 4D model is initially developed by IJL (Institut Jean Lamour) of Lorraine University [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] Gravier | Stimulated zonal flow generation in the case of tem and tim microturbulence[END_REF][START_REF] Gravier | Transport hysteresis and zonal flow stimulation in magnetized plasmas[END_REF][START_REF] Drouot | Study of turbulence associated with trapped particles in fusion plasmas[END_REF] in collaboration with CEA IRFM Cadarache as a reduced form of the 5D gyrokinetic model [START_REF] Garbet | Entropy production rate in tokamaks with nonaxisymmetric magnetic fields[END_REF]. A quick and detailed introduction of this model is given in chapter 2 section 2.2.

Simplifications in fluid turbulence

Even though the trapped particle model is much simpler than 5D gyrokinetics, the 4D simulation is still expensive, especially if a wide range of scales are to be considered, due to the nonlinear nature of this model. Thus, further simplifications, especially on the nonlinear interactions, without losing any essential features of the initial model is useful.

The study of kinetic plasma turbulence is relatively new, however turbulence is everywhere in nature and its study dates back to ancient Greeks, yet important developments have been made in this field over the years. In plasmas, there exist a number of descriptions, from the full Klimontovich description to simple reduced fluid models such as Hasegawa-Mima [START_REF] Hasegawa | Pseudo-three-dimensional turbulence in magnetized nonuniform plasma[END_REF], or Hasegawa-Wakatani equations [START_REF] Hasegawa | Plasma edge turbulence[END_REF]. Furthermore, different variants of reduced fluid models such as MHD [START_REF] Plunian | A non-local shell model of hydrodynamic and magnetohydrodynamic turbulence[END_REF][START_REF] Gürcan | Nested polyhedra model of turbulence[END_REF], reduced MHD, Hall MHD etc. have been studied into details in the context of solar wind turbulence, and of the solar dynamo problem.

Inspired by previous works, especially the description of nonlinear effects in fluid turbulence, a kinetic reduced form of trapped particle model has been worked out during this thesis. Different from the generic gyrokinetic simulation in real space, the method here is to resolve the kinetic Vlasov-Poisson system in Fourier space, where the Poisson bracket becomes a convolution after the Fourier transform. The convolution is associated to a constraint on the wave vectors: k+p+q=0, which means that the evolution of mode k is determined by the mode p and the mode q such that p+q=-k. Note that this standard constraint gives rise to the so called nonlinear multiscale interactions, which is also the most challenging part of the turbulence problem: all the difficulty and complexity comes from this term.

A promising technique used in fluid turbulence is the logarithmic discretization for the wave number k:

k = k n = k 0 g n ,
where k 0 is a free parameter, g > 1 is the logarithmic spacing factor and n is the shell number. This method allows to consider a very large range of scales in numerical simulations, which is practically impossible for the gyrokinetic simulation in real space.

In order to run the simulation, an analytic expression of the nonlinear terms should be given. Note that in general the relation of the wave vector k+p+q = 0 corresponds exactly to a triangle. Based on the existence condition of a triangle, which is principally determined by the moduli of the wave vector (i.e. wave number), a selection of all the possible triads in a given range of wave numbers has been worked out, which allows a full description of the nonlinear terms based on the triad information. In practice only a subset of all the possible interactions is taken into account in simulation, to simplify the nonlinear expression and to save time. For example the popular subset of the local interactions that is used in GOY [START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF] or Sabra models [START_REF] Vov | Differential model for 2D turbulence[END_REF]:

{n -2 , n -1 , n , n + 1 , n + 2} ,
which means that the nonlinear interactions are between the nearest neighboring modes. Simulation presents that the logarithmical discretization of the wave number space and the assumption of local interactions are rather efficient and some important results have been found based on this method, which is comparable to the direct nonlinear simulation in real space.

The aims of this thesis is to give a detailed and analytical formulation of the nonlinear multiscale interactions by using previous methods that exist in fluid or MHD turbulence problems and then apply these new nonlinear descriptions to the 4D trapped particle turbulence model.

Even though kinetic plasma turbulence is different from the fluid turbulence, the nonlinear description here doesn't destroy the kinetic effects, for example toroidal precession resonance is well observed in these simulations.

Moreover, inspired by the logarithmically discretized model (LDM) [START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF], which has been developed to resolve the 2D passive scalar equation, a further description of the nonlinear terms respecting the conservation of entropy and electrostatic potential energy has been worked out in a 2D polar coordinates based on the triad condition. This description is more general and complete, where different approaches connected to the previous models like GOY, LDM, etc can be found: it can be seen as a development of the previous ones.

This report is organized as follows: in chapter II the trapped particle model is introduced quickly, as well as the Fourier transformation of the model equations. An analytic expression of the nonlinear terms as well as their connections to the previous models has been presented too. Chapter III is contributed to discuss the linear phase of the system, such as the linear dispersion relation, the threshold of parameters and the linear instability. A numerical method based on the argument principle has been developed to resolve the linear dispersion relation, where singularity and multi roots co-exist. In chapter IV and chapter V the nonlinear simulation is presented for the isotropic model and the anisotropic model, respectively. A conclusion and a perspective is given in chapter VI as the last chapter.

Chapter 2

Bounce averaged gyrokinetics δ f equations in Fourier space.

Introduction

Nonlinearity is the most challenging aspect of kinetic plasma turbulence, yet in order to make fusion a reality, it is necessary to understand the nonlinear mechanisms of the turbulent phenomena in fusion reactor. Our ambition in this thesis is to try some new methods to represent the nonlinearity of the kinetic plasma turbulence, since direct nonlinear simulations of the 5D gyrokinetic system is expensive. This kind of large simulation is largely impractical for the study of multi-scale physics for most of the usual laboratories or institutes, which, as a consequence, limits the understanding of this problem.

Inspired by previous efforts to treat the nonlinearity in fluid turbulence, for example the Hasegawa-Wakatani model [START_REF] Hasegawa | Plasma edge turbulence[END_REF], the simple shell models such as GOY [START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF], Sabra [START_REF] Vov | Improved shell model of turbulence[END_REF], as well as the more advanced and recently developed 2D LDM [START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF] and the 3D Nested polyhedra model of turbulence [START_REF] Gürcan | Nested polyhedra model of turbulence[END_REF] etc, a new description of the nonlinear kinetic plasma turbulence will be built to interpret the multiscale interaction between the particle distribution function and the electrostatic potential fluctuations. As a first step, this description will be applied to the trapped particle turbulence [START_REF] Kadomtsev | Trapped particles in toroidal magnetic systems[END_REF][START_REF] Cohen | Non-linear saturation of the dissipative trapped-ion mode by mode coupling[END_REF][START_REF] Tang | Residual trapped-ion instabilities in tokamaks[END_REF] model, which, as a reduction of the standard 5D gyrokinetics [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF], provides a simple testbed.

Plasma turbulence is usually more complicated than the fluid turbulence, because in addition to the spatial phenomenon, the kinetic effects, such as the kinetic resonance, is also contained in kinetic plasma turbulence. Note that this work does not do any simplification on the kinetic aspects of plasma turbulence, so all the models beyond this work are totally kinetic even though we use ideas of fluid turbulence.

This chapter is organized as follows: in section (2.2) the trapped particle model [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF], or as it is more generally called bounce averaged gyrokinetics, will be introduced, where we will detail the Vlasov equation, the quasineutrality equation as well as the normalization of this model. In section (2.3) we will develop a formulation to describe the nonlinear terms of the bounce averaged gyrokinetic model in a special log polar coordinate system. This new developed formulation is able to describe all the nonlinear couplings in a very large range of scales as well as handle anisotropy and it can be applied to other turbulence problems even though we begin with the trapped particle model. Since the formulation developed here contains all the nonlinear couplings, different reduced models can be found based on different assumptions, which shows a clear connection to the reduced models of fluid turbulence. The conservation of the quadratic quantities, such as the electrostatic potential energy and the entropy, is demonstrated in section (2.4). Note that the conservation laws, as an essential nonlinear property of the Poisson bracket must be valid by our models.

Bounce averaged gyrokinetics

The bounce averaged gyrokinetics model [START_REF] Depret | Trapped-ion driven turbulence in tokamak plasmas[END_REF][START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF] is initially developed by the Institute of Jean Lamour (IJL) in University of Lorraine, in collaboration with IRFM (Institut de Recherche sur la Fusion par confinement Magnetique) in CEA Cadarache, as a reduced form of the 5D gyrokinetics system, to study the TEM and TIM turbulence in Tokamak [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] Gravier | Stimulated zonal flow generation in the case of tem and tim microturbulence[END_REF][START_REF] Gravier | Transport hysteresis and zonal flow stimulation in magnetized plasmas[END_REF]. In addition to the gyro average operator that is widely used in gyrokinetics, the bounce average operator is also used in this model to further simplify the system, which finally results in a 4D (i.e. 2 spatial coordinates + 2 kinetic coordinates ) kinetic plasma turbulence model.

Model equations

A full description of a kinetic ion-electron plasma requires solving the 6D Vlasov equation for electrons and ions, which can be given in the action-angle coordinates (J i , α i ) as follows:

∂ f s ∂t -H , f s α i , J i = 0 , (2.1)
where f s is the full particle distribution function corresponding to species s (i.e. s = i , for ions, s = e for electrons), and H is the Hamiltonian, which can be generally defined as H ≡ 1 2m s Pq s A 2 +q s Φ, with P , the canonical momentum; A, the magnetic vector potential; Φ, the electric potential; m s and q s , the mass and the charge of the particle, respectively. H , f α i ,J i is the Poisson bracket defined as:

H , f α i ,J i ≡ i ∂H ∂α i ∂ f ∂J i - ∂H ∂J i ∂ f ∂α i , ( 2.2)
where α i and J i correspond respectively to the angles (or phases) and actions associated to the cyclotron motion (i = 1), the bounce motion (i = 2) and the toroidal precession motion (i = 3), which are defined as [START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF]:

J =    J 1 = - m s µ q s J 2 = ¸ms v g ∥ 2π d s J 3 = m s R v g ∥ + q s ψ    α =   α 1 = ω c t + cnt α 2 = ω b t + cnt α 3 = ω d t + cnt   (2.3)
Table 2.1: Standard spatial and temporal parameters of ions and electrons for the ITER machine [START_REF] Drouot | Study of turbulence associated with trapped particles in fusion plasmas[END_REF].

where µ is the usual adiabatic invariant defined as µ ≡

m s v 2 ⊥
2B with v ⊥ , the perpendicular velocity, v g ∥ is the parallel velocity of the guiding center. ω c , ω b and ω d correspond respectively to the frequencies of cyclotron motion, bounce motion and toroidal precession motion. ψ is the poloidal flux and stands for the radial coordinate since d ψ = -R 0 B θ d r . Note that the action variables J i can be considered as adiabatic invariants because the variation of the actions J i are much slower than the time scales of the periodic movements. In this case the action variables J i can be seen as conserved during the movement:

d J i d t = - ∂H ∂α i = 0 , d α i d t = ∂H ∂J i = ω i . (2.4)
This spatial, temporal and kinetic system is 6 + 1D and is equivalent to a description in real space and velocity space Vlasov equation:

f s = f s (t , α 1 , α 2 , α 3 , J 1 , J 2 , J 3 ) .
(2.5)

Considering the motion of particles, allows further simplification. Since in a standard tokamak configuration, for example the standard parameters of ions and electrons for the ITER machine as shown in table (2.1), the cyclotron frequency ω c and the bounce frequency ω b are much higher than the toroidal precession frequency ω d (i.e. (ω c , ω b ) ≫ ω d ), it is possible to average the Vlasov equation over both the cyclotron motion and the bounce motion, if we are interested mainly by time scale of the toroidal precession motion, which is also the characteristic time scale of the TEM/TIM turbulence [START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF]. This double average is called gyro bounce average in the following. Comparing to the general gyro average, the average over the bounce motion is also taken into account, which results in a 4D system, since the phase information of both the cyclotron motion and the bounce motion are averaged out. This means that when we perform the double average, we can write:

f (t , α 1 , α 2 , α 3 , J 1 , J 2 , J 3 ) ⇒ f (t , α 3 , J 1 , J 2 , J 3 ) . (2.6)
Inserting the definitions of the action variables in equation (2.3), the distribution function can be represented with the usual variables:

f (t , α 3 , J 1 , J 2 , J 3 ) ⇒ f t , α, µ, E , ψ , (2.7)
here α = α 3 is the toroidal precession angle, µ is the adiabatic invariant associated with the action of cyclotron motion J 1 , E is the particle energy associated with the action of bounce motion J 2 and ψ is the poloidal flux associated with the action of precession motion J 3 . The resulting Vlasov equation after this double average, depends only on 2 variables (precession angle α and poloidal flux ψ, which plays the role of the radial coordinate r in this configuration.), parametrized by the particle energy E and the adiabatic invariant µ, which gives 2 spatial coordinates and 2 coordinates in velocity space, and thus a complete problem in 4D.

Considering the dynamics of each motion in the case of electrostatic fluctuations, the Hamiltonian H of the particle can be given as follows:

H = J 3 Ω d + q s φ + ... , (2.8)
where J 3 Ω d is the kinetic energy connected to the toroidal precession motion with Ω d , the precession frequency of the particle. q s φ is the electrostatic potential energy with φ, the electrostatic potential. There should be other potentials associated to the bounce motion and the cyclotron motion, however after the gyro bounce average, these potentials will not pass to the further calculation, so they are not detailed here. Substituting the Hamiltonian into Vlasov equation and after some manipulations, it results in [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]:

∂ f s ∂t -J 0s φ, f s α,ψ + E Ω d Z s ∂ f s ∂α = 0 , (2.9)
where f s are the full (equilibrium and fluctuation) particle distribution functions associated with the species s that are averaged over bounce and cyclotron motions. Ω d E /Z s is the toroidal precession frequency of a particle with the energy E and the atomic number Z s . J 0s is the gyro bounce average operator (which will be defined in Fourier space in the next section). Note that in the calculation of equation (2.9) we use d H

d α i = 0 and d H d J 1,2 = 0 .

Scale separation

It is common to make an assumption of scale separation between background profiles and small scale fluctuations that make up the turbulence in Tokamaks. Considering the scale separation between the equilibrium and the fluctuations, one can write the bounce center distribution function as:

f s E , µ, α, ψ, t = F s E , µ, ψ + δ f s E , µ, α, ψ, t . (2.10)
where F s is the equilibrium distribution function of the species s, which, in this model, is assumed to be Maxwellian and independent of time:

F s E , µ, ψ = n eq,s ψ T 1.5 eq,s ψ e - µB + 1 2 m s v 2 ∥ T eq,s .
(2.11)

Inserting this definition into the Vlasov equation, one can obtain the equation of the fluctuating part δ f s :

∂δ f s ∂t - ∂J 0s φ ∂α ∂F s ∂ψ + E Ω d Z s ∂δ f s ∂α -J 0s φ, δ f s α,ψ = 0 . (2.12)
Note that this is the Vlasov δ f s equation of the bounce averaged gyrokinetics model. In order to close the system, this equation is coupled to Maxwell equations with self-generated and imposed fields, which in the case of electrostatic fluctuations with no externally applied electric fields becomes the quasi-neutrality equation:

δn i = δn e , (2.13)
for scales that are larger than the Debye length. In tokamak geometry, a charged particle can be either "passing" or "trapped" in the low field side depending on its kinetics energy, as shown in figure (1.4), so the density of the particle can be seen as the sum of two parts:

δn s = δn t + δn p , (2.14)
where δn t is the density of the trapped particles, with the index t standing for "trapped", and δn p is the density of the passing particles, with the index p standing for "passing". For the trapped particles, dynamics is kept kinetic, which can be obtained from the integration of the particle distribution function over the velocity coordinates:

δn t = ˆqs J 0s δ f s E d E d µ .
(2.15)

Note that the integration with respect to the adiabatic invariant µ gives the trapped particle ratio f t , which should be between 0 and 1. Only one value of f t is passed to simulation, such that this 4D system can be simplified further to 3D. The correction of the gyro bounce average should also be considered (the difference of the real density and the density of the bounce center), which means:

δn t = f t ˆqs J 0s δ f s E d E - q s T s ˆ 1 -J 2 0s φF s E d E . (2.16)
The details of the gyro correction can be found in appendix C.

If the ratio of the trapped particle is f t , then the ratio of the passing particle is 1f t . In order to simplify the system, the passing particles are assumed to be adiabatic in this model, which means:

δn p = -1 -f t q s φ -ǫ φ φ α T s , (2.17)
where φǫ φ < φ > α is due to the response from zonal flow and ǫ φ defines this response, which takes values between 0 and 1 for the zonal modes and is zero for the other modes, because zonal flows have no direct influence on them. Inserting all the ingredients into the density equation, the final quasineutrality equation can be written as follows:

2 π 1 T e a R 0 ˆ F e 1 -J 2 0e φ + 1 τ F i 1 -J 2 0i φ E d E + (1 -f t )n 0 f t T e a R 0 1 + 1 τ φ -ǫ e,φ + ǫ i ,φ τ φ α = 2 π ˆ J 0i f i -J 0e f e E d E . (2.18)
where τ is the ion to electron temperature ratio defined as τ = T i T e . ǫ i ,φ and ǫ e,φ define respectively the response from zonal flow for ions and electrons. In order to simplify, we use q i = 1 and q e = -1 respectively for the charge of ions and electrons.

Note that the primary focus of this thesis is the development and the discussion of the reduced nonlinear descriptions, so we do not discuss the derivation of this model in great detail, we recommend the interested readers to find the details in [START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF] .

Normalization

It is important to provide proper normalization in order to identify spatial and temporal scales associated with the phenomena that can be studied in bounce averaged gyrokinetics. Here the energy is normalized by a temperature T 0 , which in this thesis refers to the ion temperature T i :

Ê = E T 0 . (2.19)
The time is normalized by a frequency ω 0 typically in the time scale of the toroidal precession motion:

ω 0 = T 0 eR 2 0 B θ , (2.20)
with R 0 , the major radius of the Tokamak and B θ , the strength of the poloidal magnetic field. This frequency represents the ionic frequency of the toroidal precession motion with the temperature T 0 , such that:

t = ω 0 t . (2.21)
Note that this frequency is also the characteristic frequency of the TIM turbulence. The distance is normalized by a scale defined as:

L ψ = a | d ψ d r |= aR 0 B θ , (2.22)
with a, the minor radius of the Tokamak. This length can be seen as the radial scale of the problem:

ψ = ψ L ψ . ( 2

.23)

Since the poloidal flux plays the role of radial coordinate: ψ = ψ (r ). The electrostatic potential is normalized by ω 0 L ψ :

Φ = Φ ω 0 L ψ = R 0 a eφ T 0 .
(2.24)

The particle distribution function f s is normalized as follows:

f = f s n 0s 2πT 0s m s - 3 2 .
(2.25)

The bounce averaged gyrokinetics model after these normalizations can be written as follows:

∂ f i ∂t - ∂J 0i φ ∂α ∂F i ψ ∂ψ + E Ω d Z i ∂ f i ∂α -J 0i φ, f i α,ψ = 0 , ∂ f e ∂t - ∂J 0e φ ∂α ∂F e ψ ∂ψ - E Ω d Z e
∂ f e ∂α -J 0e φ, f e α,ψ = 0 .

(2.26)

c n φ = ˆ J 0i f i -J 0e f e E d E . (2.27)
In this equation and from here on, we use f s instead of δ f s to simplify the notation. Here, the prefactor c n in the quasi-neutrality relation is defined as:

c n = a R 0 1 T i ˆ 1 -J 2 0i e -E E d E + a R 0 1 T e ˆ 1 -J 2 0e e -E E d E + π 2 (τ + 1) a R 0 (1 -f t ) f t 1 -ǫ φ , (2.28)
which is the sum of the polarization (i.e. gyro bounce correction) and the adiabatic response (i.e. passing particle contribution).

δ f equations in Fourier space

The nonlinear gyrokinetics is the domain of massively parallel high performance numerical simulations through the computation of nonlinear dynamics within a description, which includes self-consistent multi-scale interactions. It usually requires millions of CPU hours and makes anything beyond medium nonlinear runs dominated by a single type of instability over a small range of scales, impractical. As a result of the complexity of numerical implementation, sometimes it is hard to understand and isolate important physical mechanisms. In this regard, reduced models appear as intermediate tools, that are necessary for isolating important physical mechanisms, to provide guidance to large scale gyrokinetic simulation efforts as well as comparison with experiments.

Based on these ambitions, the following reduced model, which stands on the simplification of the nonlinear dynamics, without losing any essential features, has been worked out for the trapped particle turbulence model. In this thesis the simulation is implemented in Fourier space, which means:

f s α, ψ, µ, E = k f s k µ, E e i k α α+i k ψ ψ , ( 2.29)
Inserting this definition into the Vlasov equation, one can obtain the Vlasov equation in the wave numbers plane k= k ψ , k α as follows:

∂ f s k ∂t -i k α J 0s φ k ∂F s ∂ψ + E Ω d Z s i k α f s k - k+p+q=0 p ψ q α -p α q ψ J 0s φ * p f * s q = 0 . (2.30)
The last term of equation (2.30) is the convolution, which corresponds to the Fourier transform of the Poisson bracket. It means that the evolution of the mode k is determined by the mode p and the mode q that should satisfy the wave vector constraint p + q + k = 0. This constraint defines the nonlinear interactions across scales in turbulence, which presents infinite possibilities and thus a challenge for both theory and simulation.

Nonlinear interactions among scales has been discussed since the beginning in the context of fluid turbulence [START_REF] Hasegawa | Pseudo-three-dimensional turbulence in magnetized nonuniform plasma[END_REF], as well as in MHD turbulence [START_REF] Plunian | A non-local shell model of hydrodynamic and magnetohydrodynamic turbulence[END_REF]. Various reduced models that describe these scale by scale interactions have been proposed such as the shell models or LDM, or Nested polyhedra model [START_REF] Gürcan | Nested polyhedra model of turbulence[END_REF]. Inspired by these previous ideas, a reduced representation of the convolution will be given based on the description of the wave vector constraint.

Another important point in above equation is the double average operator J 0s . In real space, the bounce average, as well as the gyroaverage, correspond to operators acting on the unknowns, and the double averaging procedure has been symbolized in equation (2.9) by J 0s φ. In Fourier space, the bounce-average and gyro-average can be obtained as simple zeroth order Bessel function of the first kind, so that the total average, symbolized by J 0s for the sake of simplicity in equations (2.18, 2.30), is given by the following product:

J 0s ≡ J 0 (k ψ δ s0 E ) J 0 (k α ρ s0 E ) .
(2.31)

In order to simplify the simulation, different approximations of the average operator has been used based on different scales in previous works [START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF]. In appendix D we will discuss the difference between the exact average operator and its different approximations, especially the validation of its approximations in very small scales.

Description of nonlinear interactions

In theory, if the medium is infinite, the wave vector k is a continuum quantity, so in the Fourier plane there are infinite number of waves and in the absence of quantization due to boundary conditions, the number of triads defined by k+p+q = 0 is also infinite, which results in an infinite number of nonlinear interactions. The nonlinearity is the most complicated part in turbulence and analytic solution for such a problem is probably impossible.

However in numerical simulation, the wave vector should be discretized in a cartesian or polar coordinate, linearly or logarithmically, and the number of waves is not infinite. For a given range of discretized waves, the number of interacting triads is fixed, which can be found out easily by a program based on the conditions fixed by the wave vector constraint: k+p+q = 0. The triad information then can be written in a table as a coupling card. Based on this coupling card, it is possible to give an analytic expression of the nonlinear terms. Of course even with an analytic expression, the system is still very complicated and it is difficult to find a solution by theoretical analysis, due to the average operator (Bessel function) and the integration over the kinetic coordinates, etc, but it can be resolved numerically.

In order to give an analytic expression of the nonlinear terms, it is necessary to start from the coupling information defined by the wave vector constraint: k+p+q = 0. Figure 2.1: k + p + q = 0 in a plane, which in fact forms a triangle with sides of length k, p and q and angles of α k , α p and α q , respectively, as shown by the left figure. The right figure shows the existence conditions of a triangle: the maximum of (k, p, q) should be smaller or equal to the sum of the two smaller ones; the minimum of (k, p, q) should be larger or equal to the difference of the two larger ones.

Description of k+p+q = 0

The wave vector constraint k+p+q = 0 is very common in fluid turbulence [START_REF] Hasegawa | Plasma edge turbulence[END_REF][START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF], kinetic plasma turbulence [START_REF] Xu | A turbulent cascade model of bounce averaged gyrokinetics[END_REF], as well as in MHD [START_REF] Plunian | A non-local shell model of hydrodynamic and magnetohydrodynamic turbulence[END_REF]. A detaileddescription of the Poisson bracket in Fourier space ( i.e. equation (2.30)) will be given in this section. The vectors are presented in polar coordinates, which means the wave vector k is defined as k = (k, θ k ), with k, the modulus of the wave and θ k , its orientation. In this coordinate system k α and k ψ can be expressed as:

k α = ksin (θ k ) , k ψ = kcos (θ k ) .
(2.32)

So k = k ψ , k α forms a cartesian coordinate and k = (k, θ k ) forms a polar coordinate one. The nonlinear terms can be represented in vector form, independent of components as follows:

p+q+k=0 p ψ q α -p α q ψ J 0s φ * p f * s,q = p+q+k=0 ẑ × p • qJ 0s φ * p f * s,q , (2.33)
where p = p, θ p and q = q, θ q . ẑ is the direction that is perpendicular to the wave number plane k, where k ψ , k α , ẑ forms a 3D cartesian coordinate. The wave vector constraint p + q + k = 0 in the Fourier plane is shown in figure(2.1), which in fact describes triangles, with the sides of length k, p and q and angles α k , α p and α q , respectively. The interaction coefficient is defined as M kpq = ẑ × p • q , which can be calculated as:

M kpq = ẑ • p × q = pq sin θ p -θ q .
(2.34)

The absolute value of the interaction coefficient M kpq equals the surface of the triangle formed by k, p and q. In reduced models, like GOY, Sabra and LDM, only local interactions (i.e. p, q ∼ k) are considered. The obvious motivation for this is to simplify the nonlinear description, however a more important and concrete reason for this is that the surface of the triad formed by neighboring waves is maximum, which, as a reduced model, can qualitatively represent the maximum nonlinearity.

Another important situation is the disparate scale interaction [24][47], where p ∼ k and q ∼ 0. This nonlinear interaction represents the direct transfer between the large scale modes, like the zonal flow (q = 0) and the small scale modes, which is another important mechanism in turbulence.

Note that M kpq and M kqp is different due to their opposite sign, as will be shown later that the opposite sign is very important to satisfy the conservation laws.

In general, the constraint p + q + k = 0 defines a triangle. In order to compose a triangle, the moduli (or length) of the vectors should satisfy some constraints as shown in figure(2.1), which are: i) the sum of any two edges must be larger than (or equal to, in the limit case) the third one.

ii) the difference of any two edges must be smaller than (or equal to, in the limit case) the third one.

These two criteria are the necessary and sufficient conditions to determine the existence (or not) of a triangle if the lengths of the three sides are given. Based on these criteria, it is possible to select all the couples p, q that can interact with the wave number k for a given range of wave numbers. Note that the selection process acts only on the wave moduli k, p and q. In a triangle, if the three sides are given, the angles can be fixed by the law of the cosines, then this triangle is fixed, which in fact means the phase information is fixed.

The phase information θ p and θ q of the wave p and q can be obtained from the wave k by the following relations ( see figure(2.1)):

θ p = π + θ k + α q , θ q = π + θ k -α p , (2.35)
here π + θ k means the inverse direction of the wave k, α p (α q ) is the angle between the wave q (p) and -k (see figure(2.1)), which can be calculated by the cosines relation of a triangle:

α q = arccos q 2 -p 2 -k 2 2pk , α p = arccos p 2 -q 2 -k 2 2qk . (2.36)
From equation (2.35), we can obtain θ p -θ q = α q +α p , since α p +α q +α k = π, so θ p -θ q = πα k , which will be used in the calculation of interaction coefficients M kpq . Finally the information of the waves k, p, q as well as the interaction coefficient M kpq in a polar coordinate are given as follows:

k = (k, θ k ) , p = p , π + θ k + α q , q = q , π + θ k -α p , M kpq = pq sin (π -α k ) .
(2.37)

In order to know the coupling information of a triad, we should have k, p, α p , q, α q , M kpq , or at least the wave numbers k, p, q , since the other elements, such as the phase and the interaction coefficients can be calculated simply from the wave numbers.

Description of the nonlinear terms

When the coupling information is known, a description of the nonlinear terms can be written. A formulation of Poisson bracket in Fourier space will be laid out in this section based on the definition of the triad.

Considering the permutation of p and q (i.e. φ p f q ↔φ q f p , because the evolution of f k can be determined by φ p and f q , it can be also determined by φ q and f p ), one can write (i.e., see figure(2.2)) the nonlinear terms as:

ẑ × p • qJ 0s φ * p f * s,q = M kpq J 0s φ * p f * s,q + M kqp J 0s φ * q f * s,p . (2.38)
Considering the vectors that are mirror symmetric with respect to the vector k, for example the reflection vector p ′ (q ′ ) of p (q) with respect to k , as shown in figure(2.2), is:

p = p , π + θ k + α q k ← → p ′ = p , π + θ k -α q , q = q , π + θ k -α p k ← → q ′ = q , π + θ k + α p . (2.39)
This allows us to write the nonlinear term as: .40) Note that this expression gives a full description of the nonlinear terms for one triad k, p, q in the wave plane, which includes the permutation and the reflection of the wave vectors with respect to k, and there are no other possibilities for the same triangle.

ẑ × p • qJ 0s φ * p f * s,q =M kpq J 0s φ * p f * s,q + M kqp J 0s φ * q f * s,p + M kp ′ q ′ J 0s φ * p ′ f * s,q ′ + M kq ′ p ′ J 0s φ * q ′ f * s,p ′ . ( 2
From the definition of M kpq and the relation between the wave vectors p, q and p ′ and q ′ , we have the following relations:

k p q α k α p α q Φ p f q f p p q k α k α p α q Φ q α q α k k α p f q Φ p Φ q f p α q p q α k k α p permutation mirror reflection q p
M kpq = -M kqp , M kpq = -M kp ′ q ′ .
(2.41)

Inserting the definitions of the waves p, q, p ′ , q ′ and the interaction coefficients into equation (2.40), the nonlinear terms of one triad in polar coordinate can be represented as:

ẑ × p • qφ * p h * q = pq sin (π -α k ) (φ * π+θ k +α q p h * π+θ k -α p q -φ * π+θ k -α p q h * π+θ k +α q p -φ * π+θ k -α q p h * π+θ k +α p q + φ * π+θ k +α p q h * π+θ k -α q p
) , (2.42) With these, the Vlasov equation can be expressed in polar coordinates as follows:

∂ f s,k ∂t =i k α J 0s φ k ∂F s ∂ψ - E Ω d Z s i k α f s,k + k(p, q) pq sin (π -α k ) φ θ k +α q p h θ k -α p q -φ θ k -α p q h θ k +α q p -φ θ k -α q p h θ k +α p q + φ θ k +α p q h θ k -α q p , (2.43)
here k p, q means all the couples p, q that can interact with the mode k . So the nonlinear terms should be the sum of all the possible couplings, which is represented by in equation (2.43). pq sin (πα k ) is the interaction coefficient of the corresponding triad. Note that this expression gives a full description of the Fourier transform of 2D Poisson bracket in polar coordinates, which includes all (not only a subset of ) the possible couplings (i.e. local couplings, disparate scale couplings and other types) in a given scale. Before doing the simulation, one needs to find out all the couples k p, q that gives the coupling card of system. After having the coupling card, the Vlasov equation can be simulated directly.

Different approaches and connections to previous models

In the previous section, a description of the nonlinear terms (i.e. Fourier transform of the Poisson bracket) is worked out in the wave number plane. In simulation, the wave numbers k, p, q can be linear, logarithmic or even random, which is not a limit in this description. However in practice, one must consider the capacity (i.e. memory, threads,etc) of the machine and the intention of the research. For example if the nonlinear physics in small scale is of interest, it may be necessary to use logarithmic discretization of the wave number magnitudes, to describe a large range of scales easily. The number of the couplings in logarithmic discretization is much less than that in linear discretization, which will take less time and consume less computer resources to run the simulation. This was the main motivation for the development of "shell models". Logarithmic discretization means that the wave number k is discretized as:

k = k n = k 0 g n , (2.44)
where k 0 is a free parameter that defines the smallest wave numbers and g > 1 is the logarithmic scaling factor and n is the shell number. Equation (2.43) gives a full expression of the Poisson bracket in polar coordinates. Based on this expression, different reduced models can be derived, which are connected to the previous models developed in fluid turbulence.

If only local interactions are considered, for example the widely used reduction:

{(n -2, n -1, n) , (n -1, n, n + 1) , (n, n + 1, n + 2)} , (2.45)
where there are only three coupling triads in nonlinear interactions, the Vlasov equation can be reduced to the form:

∂ f j s,n ∂t = i k α J 0s φ j n ∂F s ∂ψ - E Ω d Z s i k α f j s,n + 1 2 k 2 n g -4 µ 0 J 0s φ * j +r 0 n-2 f * j -s 0 s,n-1 -J 0s φ * j -s 0 n-1 f * j +r 0 s,n-2 + J 0s φ * j +s 0 n-1 f * j -r 0 s,n-2 -J 0s φ * j -r 0 n-2 f * j +s 0 s,n-1 + 1 2 k 2 n g -2 µ 0 J 0s φ * j +l 0 n-1 f * j -s 0 s,n+1 -J 0s φ * j -s 0 n+1 f * j +l 0 s,n-1 + J 0s φ * j +s 0 n+1 f * j -l 0 s,n-1 -J 0s φ * j -l 0 n-1 f * j +s 0 s,n+1 + 1 2 k 2 n µ 0 J 0s φ * j +l 0 n+1 f * j -r 0 s,n+2 -J 0s φ * j -r 0 n+2 f * j +l 0 s,n+1 + J 0s φ * j +r 0 n+2 f * j -l 0 s,n+1 -J 0s φ * j -l 0 n+1 f * j +r 0 s,n+2
.

(2.46)

This has the same nonlinear form as the LDM [START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF], which has been recently developed to resolve the 2D passive scalar equations in fluid turbulence. In this model the wave vector is defined as k = k n , θ j with k n , logarithmically discretized and θ j , linearly discretized as θ j = j * 2π M , which means that the unit circle in wave number space is divided into M regular sections here. n is the index of the shell number and j is the index of the angle, which should be an integer between 1 and M .

In equation (2.46) r 0 , l 0 and s 0 are the coupling angles (or rather integer index offsets that correspond to these angles) originated from α p and α q , which should be normalized by 2π M in simulation. Note that the coupling angles gives the exact phase information of this 2D cascade model.

Since the particle distribution function f as well as the electrostatic potential φ are real physical quantities, their Fourier components should satisfy the following relations:

f -k = f * k , φ -k = φ * k , (2.47)
which allows to simulate the system only in the upper half wave plane [0; π] in Fourier space, the quantities in the other half plane [π; 2π] can be obtained automatically from equation (2.47).

The nonlinear simulation of this model will be presented in chapter (V).

Considering the isotropic fluctuations such that k α = k ψ , the 2D anisotropic turbulence problem can be simplified to a 1D problem. In this case the phase information θ j is lost and one has the freedom to choose different phase approximations. Note that the fact that there is no exact phase information, is a general problem in 1D shell models. Shell models also have a free parameter α in the nonlinear interaction coefficients. Since the model becomes 1D , the coupling angles r 0 , l 0 ,s 0 , etc, of this model will be 0 or π, depending on different phase approximation, which finally results in different models. For example if we use the well known GOY version phase approximation, the isotropic cascade model [START_REF] Gürcan | Wave-number spectrum of drift-wave turbulence[END_REF][START_REF] Xu | A turbulent cascade model of bounce averaged gyrokinetics[END_REF] of the bounce averaged gyrokinetics can be written as follows:

∂ f s,n ∂t =i k n J 0s φ n ∂F s ψ ∂ψ -i E Ω d Z s k n f s,n + αk 2 n g -3 J 0s φ * n-2 f * s,n-1 -J 0s φ * n-1 f * s,n-2 -αk 2 n g -1 J 0s φ * n-1 f * s,n+1 -J 0s φ * n+1 f * s,n-1 + αk 2 n g J 0s φ * n+1 f * s,n+2 -J 0s φ * n+2 f * s,n+1 , (2.48)
If we use the Sabra version [START_REF] Xu | A turbulent cascade model of bounce averaged gyrokinetics[END_REF] phase approximation, the Vlasov equation can be presented as:

∂ f s,n ∂t =i k n J n 0s φ n ∂F s ψ ∂ψ -i E Ω d Z s k n f s,n + αk 2 n g -3 J n-2 0s φ n-2 f s,n-1 -J n-1 0s φ n-1 f s,n-2 -αk 2 n g -1 J n-1 0s φ * n-1 f s,n+1 -J n+1 0s φ n+1 f * s,n-1 + αk 2 n g J n+1 0s φ * n+1 f s,n+2 -J n+2 0s φ n+2 f * s,n+1 , (2.49)
where α is a free parameter, which represents the magnitude of the linear to nonlinear terms. Since α is not a physical parameter, as will be shown later the physical properties of the system do not depend on this parameter. The nonlinear coefficients in equation (2.49) and equation (2.48) is g -3 , g -1 , g , which is exactly same as g -4 , g -2 , 1 in equation (2.46), because the isotropic models such as GOY and Sabra and the anisotropic LDM takes the same triads (i.e. equation(2.45)) in the nonlinear interactions. Since the nonlinear coefficients can be seen as the surface area of the triad, the ratio should be the same in GOY, Sabra and LDM. The interaction coefficient of GOY model is calculated from the conservation of electrostatic potential energy and entropy, which is presented in appendix E. The nonlinear coefficients of Sabra model are exactly the same as that of GOY model. The nonlinear simulation of the isotropic models will be discussed in chapter (IV) Considering the simple limit f = k 2 φ and neglecting the Bessel function, we find the simplest form of these models:

∂φ n ∂t =F -D + a n φ * n-2 φ * n-1 + b n φ * n-1 φ * n+1 + c n φ * n+1 φ * n+2 , (2.50)
where F is the forcing and D is the dissipation. This is the so called GOY model [START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF], which can be seen as the original of all the above models.

Conserved quantities.

In turbulence, it is of key importance, to identify the conserved quantities. The integration of the multiplication of Poisson bracket and any function of its arguments over the full phase space should be zero, which means:

ˆ J 0s φ, f s f s d Λ = 0 , ˆ J 0s φ, f s φd Λ = 0 , (2.51)
where d Λ refers the whole phase space, which is defined as d Λ = d αd ψd µd E . These two relations correspond to the conservation of the entropy and the conservation of the electrostatic potential energy, respectively (appendix E). These essential nonlinear properties should also be satisfied in the reduced models. The demonstration will be given in this section.

In fact in the system of triad interactions, the energy and the entropy are conserved nonlinearly among the triad, which means in quasi-stationary state:

∂ | f k | 2 + | f p | 2 + | f q | 2 ∂t = 0 , (2.52) here | f k | 2 = f k * f * k with f * k , the conjugate of f k , using the nonlinear terms of ∂ f k
∂t defined in equation(2.33) and separating the terms with and without * , we obtain:

f k ∂ f * k ∂t + f p ∂ f * p ∂t + f q ∂ f * q ∂t = + ẑ • p × q φ p f q f k -φ q f p f k + ẑ • k × q φ k f q f p -φ q f k f p + ẑ • p × k φ p f k f q -φ k f p f q = + φ p f q f k p × q + p × k + φ q f p f k -p × q -k × q + φ k f q f p k × q -p × k =0 . (2.53)
Similar calculations can be done with the others terms, finally the entropy is demonstrated to be conserved by the nonlinear terms. Due to the symmetry of φ and f in Poisson bracket, the demonstration of electrostatic potential energy conservation follows the same logic as that of entropy conservation. Note that the convolution constraint p + q + k = 0 is the key to define conservation properties of quadratic quantities in Fourier space.

Conclusion

In the present chapter, a short introduction of the 4D trapped particle model has been presented in section II [START_REF] Arnold | Cross Sections for the Reactions D(d , p)T, D(d , n)He 3 , T(d , n)He 4 , and He 3 ([END_REF], consisting of the derivation of the Vlasov δ f equation, the quasineutrality relation, as well as the normalization of this system. Since this work is focused on the reduction of the description of nonlinear physics in kinetic plasma turbulence, for which this model provides a very good testbed, the introduction was kept concise. In section II(2), the Fourier transform on the spatial coordinates α and ψ of this system has been given, where an explanation of the multiscale interactions defined by the nonlinear terms, as well as the difficulty in simulation is detailed. In section II(3), a description of the nonlinear terms in the polar coordinate of the wave number plane has been worked out, based on the nonlinear wave interaction constraint k+p+q = 0. This expression gives a full description of the nonlinear interaction, which can be seen as a development of the previous models in fluid turbulence, because from this description, the previous models like GOY , Sabra and LDM can be obtained. Finally a demonstration concerning the nonlinear conservation property, imposed by the Poisson bracket, is presented, in order to show that this nonlinear description, as well as its reduced models, respect the same conservation laws.

The Fourier transform is applied to the spatial coordinates α and ψ, the dependence on the kinetic coordinates (i.e. particle energy E and the adiabatic invariant µ appearing as the trapped particle ratio f t ) are not modified. If the reduced model works well, it should be able to keep the kinetic effect, such as the kinetic resonance. This is a challenge for the shell models, because these models are firstly developed to resolve the fluid turbulence, which is independent of the kinetic effect. As will be shown later that these models conserve all the kinetic effects. So this work gives the first working example of a kinetic shell model.

Chapter 3 Linear Phase

Before studying the nonlinear dynamics of trapped ions and electrons, it is useful to consider the characterization of the linear response of the plasma to electrostatic perturbations. Such a response provides an understanding of the role of the principal external parameters. Moreover, it allows to distinguish in the external parameter space, the nature of the instability driving the turbulence such as Trapped Ion Modes (TIM) and Trapped Electron Modes (TEM) [START_REF] Kadomtsev | Trapped particles in toroidal magnetic systems[END_REF][START_REF] Cohen | Non-linear saturation of the dissipative trapped-ion mode by mode coupling[END_REF][START_REF] Tang | Residual trapped-ion instabilities in tokamaks[END_REF].

In this model, the free physical parameters are: temperature gradient κ T , density gradient κ n , the ion to electron temperature ratio τ, the trapped particle ratio f t , and the different scales such as the electron Larmor radius ρ e , the ion Larmor radius ρ i , the electron banana width δ e and the ion banana width δ i , etc. So the free physics parameters are rich. In linear study, to have a clear understanding of the effect of one parameter, we scan the values of this special parameter and keep the values of the other parameters unchanged.

Linear dispersion relation

Linear dispersion relation, which is obtained by setting the plasma dielectric function to zero in plasma physics, gives the relationship between the linear mode frequency and the free physical parameters, in particular, the wave number. As a general convention in plasma physics, the frequency ω is a complex quantity that is defined as ω = ω r + i γ, where the real part ω r gives the wave propagation frequency and the imaginary part γ presents the growth rate ( if γ > 0) or damping rate ( if γ < 0) of the system and γ = 0 can be used to distinguish the stable (γ < 0) and unstable modes (γ > 0). So the linear dispersion relation allows to learn the effect of the external physical parameters to wave propagation frequency, as well as their influences on the instability of the linear system. In other words, all the properties of the linear system can be characterized by the linear dispersion relation.

Furthermore, the linear instability or the linear injection is very important in magnetized plasmas. As a self-organized system, the stationary state of the plasma is a balance between the linear injection, the nonlinear self-consistent multiscale interaction and the dissipation applied to the system. A quantitative measurement of the linear effect allows a better understanding of the phenomena in nonlinear phase as well as choosing the proper parameters in nonlinear simulation, such as the injection strength that is mainly determined by the source of instability, the dissipation in large scale ν L , the dissipation in small scale ν s , etc.

Plasma dielectric function: ǫ (ω)

The linear dispersion relation can be obtained from the linearized Vlasov-Poisson system (i.e, equations (2.30) and (2.18)), by setting the nonlinear terms to zero and assuming plane wave solutions for the fluctuations, i.e.

f s E , k ψ , k α , t ∝ f s k (E ) e -i ωt +i k α α+i k ψ ψ and φ k ψ , k α , t ∝ φ k e -i ωt +i k α α+i k ψ ψ : i ω f s k = -i k α J 0s φ k ∂F s ψ ∂ψ + E Ω d Z i k α f s k , (3.1)
c k φ k = 2 π f t ˆ J 0i f i k -J 0e f e k E d E , (3.2)
where the prefactor:

c k = a R 0 ˆ 1 -J 2 0i e (-ξ) ξd ξ (3.3) + τ 5 2 a R 0 ˆ 1 -J 2 0e e (-τξ) ξd ξ + π 2 (τ + 1) a R 0 (1 -f t ) f t
is introduced and the normalized energy ξ ≡ E /T i is defined for each species. k is the wave vector defined as k= k ψ , k α . By eliminating the fluctuating part of the distribution function associated with each species f s k , as well as assuming a Maxwellian form for the background distribution function

F s ψ = n s (ψ) T s (ψ)
e [-E/Ts(ψ)] , the plasma dielectric function ǫ (ω) reads as follows:

ǫ (ω) = c k - ˆ κ n + κ T ξ -3 2 Ω d (ξ -ω) J 2 0i e (-ξ) ξd ξ -τ 3 2 ˆ κ n + κ T τξ -3 2 Ω d (ξ + ω) J 2 0e e (-τξ) ξd ξ, ( 3.4)
where κ n = d ψ ln n and κ T = d ψ ln T stand respectively for logarithmic gradients of background plasma density and temperature and ω is the complex mode frequency normalized by the toroidal wave number k α .

Singularity and residue

The integrals in equation (3.4) admit singularities when the frequency is purely real (ℑ(ω)=0).

Since the Fourier-Laplace transform as given in (3.4) is only valid for ℑ(ω) > 0, residue contributions for both integrals have to be taken into account for the analytic continuation of the function ǫ (ω) into the lower half complex plane ℑ (ω) ≤ 0. Since the energy coordinate ξ is always positive, the relevant poles correspond to ξ = ω for ℜ (ω) > 0 and to ξ = -ω for ℜ (ω) < 0.

The connection with TIM and TEM becomes apparent here with the origin of these two different modes coming from the signs of electron and ion charge, inducing opposite mode frequencies [START_REF] Kadomtsev | Trapped particles in toroidal magnetic systems[END_REF]. Thus, the plasma dielectric function, which is valid everywhere in the complex plane, can be written in the form:

ǫ (ω) → ǫ (ω) + λ j i πg (ω) , (3.5)
where the ǫ (ω) on the right corresponds to the definition given in (3.4), valid only for ℑ (ω) > 0, and

λ j = 0, 1, 2 if respectively ℑ (ω) > 0, ℑ (ω) = 0 or ℑ (ω) < 0. The function g (ω) is residue contribution, which is g i (ω)
and g e (ω) for ions and electrons respectively:

g i (ω) = - κ n + κ T ω -3 2 Ω d J 2 0i e (-ω) ω ∀ ℜ (ω) > 0 , (3.6
)

g e (ω) = τ κ n + κ T -τω -3 2 Ω d J 2 0e e (τω) -ω ∀ ℜ (ω) < 0 . (3.7)
The dispersion relation, which is obtained by setting ǫ (ω, k) = 0, can then be solved numerically using a method based on the argument principle [START_REF] Johnson | Enclosing all zeros of an analytic function -a rigorous approach[END_REF], which allows us to find all the roots ω = ω (k) in the complex plane of the frequency, for any given set of parameters κ n ,κ T ,τ, f t ,k ψ ,k α . Before solving the linear dispersion relation, the parameter threshold is discussed.

The detailed calculation of the linear dispersion relation ǫ (ω) as well as its residues g (ω) can be found in appendix B.

Threshold of the temperature gradient driven instability

The instability threshold is defined as the point where the imaginary part of the roots of the dispersion relation becomes exactly equal to zero: ℑ (ω) = 0, and therefore it allows to distinguish between stable (γ < 0) and unstable (γ > 0) modes. In the multispecies model of interest here, the two branches can coexist at the threshold, corresponding to TIM or TEM. Depending on the sign of the real frequency, the TIM or TEM thresholds can be defined by setting ℑ (ω) = 0 in (3.5), (3.6) and (3.7) in the dispersion relation. Since the frequency is exactly real at the threshold, the plasma dispersion function can be seen as two parts, which are the real part given by ǫ (ω) and the imaginary part given by g (ω). Note that this is only valid for the frequency at the threshold (i.e, ℑ (ω) = 0, ω is purely real). In order that the plasma dielectric function vanish, the real and imaginary parts should both vanish, that is:

δ i δ e ρ i ρ e f t τ a/R 0 Ω d 0.1 0.01 0.03 0.01 2/3 1.0 0.1 1.0
Table 3.1: Standard parameters used for TIM and TEM studies, similar to the Ref. [START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF].

ǫ (ω) = 0 , g i (ω) = 0 , g e (ω) = 0 . (3.8)
From the residues, two explicit solutions corresponding to TIM and TEM can be found as follows:

κ n = -κ T ω - 3 2 , (3.9
)

κ n =κ T τω + 3 2 . ( 3.10)
Inserting these two into ǫ (ω) = 0, two expressions for the critical temperature gradients for TIM and TEM can be obtained:

κ T I M T = c k ´J 2 0i exp(-ξ) Ω d ξd ξ + τ 3 2
´(τξ-ω)

Ω d (ξ+ω) J 2 0e e (-τξ) ξd ξ , (3.11) κ T E M T = c k
´(ξ+τω)

Ω d (ξ-ω) J 2 0i exp (-ξ) ξd ξ + τ 5 2 ´J 2 0e Ω d e (-τξ) ξd ξ . ( 3.12)
In nonlinear phase we only discuss the system driven by the temperature gradient, which means the density gradient is 0. In the absence of the density gradient, two solutions of the frequency for TIM and TEM can be obtained:

ω T I M = 3 2 , ω T E M = - 1 τ 3 2 . (3.13)
Note that these expressions give the real normalized frequencies of TIM and TEM at the threshold (i.e, ℑ (ω) = 0) and it is only valid in the absence of density gradient. When the density gradient exists, the real frequency at the threshold is determined by both density and temperature gradients (i.e, equations(3.9) and 3.10).

For κ n = 0, the critical temperature gradient can be written explicitly as:

κ T I M T = c k ´J 2 0i exp(-ξ) Ω d ξd ξ + τ 3 2
´ τξ- 3 2

Ω d ξ+ 3 2 J 2 0e e (-τξ) ξd ξ , (3.14) κ T E M T = c k ´ ξ+ 3 2 Ω d ξ-3 2τ J 2 0i exp (-ξ) ξd ξ + τ 5 2 ´J 2 0e Ω d e (-τξ) ξd ξ . (3.15)
For a set of given parameters κ n , κ T , τ, f t , k ψ , k α , et c , the TEM and TIM critical temperature gradient can be calculated from above relations.

Threshold of κ T and the wave number k

The dependence of the critical temperature gradients with respect to the wave vector as implied by (3.11) and (3.12) is shown in figure 3.1. Here a logarithmic representation of the wave vectors has been used, such that

k ψ = k α = k n = k 0 g n , ∀n = 0, • • • , N , with k 0 = 0.1, g = 1.2, N = 20
anticipating the later use of a similar representation in the bounce averaged gyrokinetic cascade model. The other parameters chosen here are given in table 3.1.

It is shown in figure 3.1 that the critical value of the temperature gradient increases with the wave vector both for TIM as well as TEM. It can be seen in figure 3.1 also that, for these parameters, for a given wave vector k n , the critical temperature gradient is larger for TIM than for TEM, meaning that the TEM instability appears easier to trigger than the TIM. It can also be noticed that the threshold increase with wave number, meaning that large scale modes are easier to destabilize. Accordingly the lowest value of the threshold is obtained with the mode k n = 0. However when the growth rate is considered as a function of the wave number, the k n = 0 mode is not the most unstable one (see for example figure 3.5). In any case, since in the sense that it is the first mode to become unstable, the following analysis will be considered with this value for the wave vector.

Threshold of κ T and ion to electron temperature ratio τ

Even though the TIM linear threshold is found higher than the TEM linear threshold for the set of parameters given in figure 3.2, this is in fact only valid for values such that: τ ≥ 1.0. The effect of the ion to electron temperature ratio is investigated in figure 3.2 in order to fully understand the relationship between TIM and TEM linear thresholds. Note that the interval τ ∈ [0.5 ; 2] is clearly sufficient to cover most of the temperature ratios observed experimentally [START_REF] Casati | Temperature ratio dependence of ion temperature gradient and trapped electron mode instability thresholds[END_REF].

Figure 3.2 shows the TIM and TEM critical temperature gradients as functions of the ion to electron temperature ratio. Four regions appear, where TIM and TEM become alternatively stable or unstable. The TEM are found easier to destabilize while increasing τ, on the other hand the TIM instability is harder to destabilize for increasing values of τ.

Threshold of κ T and the trapped particle ratio f t

Finally, we consider the effect of trapped particle ratio f t on TIM and TEM stability. The result is shown in figure (3.3), for three different values of the ion to electron temperature ratio τ, where one can clearly see that the threshold decreases with trapped particle ratio f t for both TEM and TIM, which means that the system will be more unstable if more particles are trapped. This of course is valid within the constraints of the model where f t should be in a reasonable range (for example f t = 1 being meaningless). Also, the effect of the temperature ratio can be observed again in this figure : increasing τ leads to a stabilization of TIM modes, and on the contrary to a destabilization of the TEM modes.

Linear instability

As is the convention in linear stability analysis in fusion plasma the imaginary part of the eigenfrequency ℑ (ω) gives the linear growth rate of the associated instability if positive (ℑ (ω) > 0), or linear damping if negative (ℑ (ω) < 0), while the sign of real part ℜ (ω)allows to discriminate between TIM (if positive) and TEM (if negative).

Take the wave oscillating solution for the electrostatic potential φ, which means

φ (t ) = φ k e -i ωt .
As the frequency is complex (i.e. ω = ω r + i γ), the electrostatic potential φ can be expressed as:

φ (t ) = φ k e -i (ωr +i γ)t = φ k e -i ω r t +γt .
So the differentiation of φ with respect to the time t gives:

∂φ (t ) ∂t = ∂ φ k e -i ω r t +γt ∂t = -i ω r + γ φ (t ) .
This solution means that in the evolution of electrostatic potential φ as a function of time, one should be able to see the growing (or damping) process due to γ as well as the oscillation due to ω r .

The linear dispersion relation is a multi-root problem, which means that there are many possible roots in the complex plane ω = ω r +i γ. So the final solution of the electrostatic potential should be the sum of all the possible modes:

φ (t ) = φ 1 e -i ω 1 t + φ 2 e -i ω 2 t + ... + φ n e -i ω n t ,
with ω 1 ,ω 2 ,...,ω n , the eigenfrequencies and φ 1 ,φ 2 ,...,φ n the initial state of the corresponding eigenmode. However this solution is always dominated by the mode with the maximum growth (or damping) rate due to the exponential e γt when time tends to infinity, which means:

φ 1 e -i ω 1 t + φ 2 e -i ω 2 t + ... + φ n e -i ω n t t → ∞ ----→ φe i ωt ,
with ω equals to the mode of max [ℑ (ω 1 ) , ℑ (ω 2 ) , ..., ℑ (ω n )], so it is only necessary to find the eigenfrequency with the maximum imaginary part ℑ (ω).

In the eigenvalue problem, since the system is dominated by the mode with maximum γ after a short time evolution, we are only interested by one root for both TIM and TEM, whose imaginary part is maximum, even though the numerical method allows to calculate all the possible roots.

Numerical method: argument principle

As there is no compact analytical solution, the linear dispersion relation(equation (3.5)) must be resolved numerically. Since two kinds of instability TIM and TEM coexist in linear system, it is necessary to find the two roots separately. A method from complex analysis [START_REF] Johnson | Enclosing all zeros of an analytic function -a rigorous approach[END_REF] called argument principle has been introduced and numerically developed here. This principle states that for a given function f (z) the number of roots in the fixed rectangular domain R can be obtained by the integral where, f ′ (z) is the derivative of f with respect to the argument z. N and P are the numbers of zeros and poles of f in the rectangle R, respectively. For an analytical function like the plasma dielectric function there is no pole, i.e. P = 0. So this integral gives the number of roots.

´R ∂ z f (z) f (z) d z, which means: 1 i 2π ˆR f ′ (z) f (z) d z = N -P , (3.16)
If the number of roots is zero in a rectangle R, the calculation is terminated in this region. If there is only one root in R , some modules developed for optimization problem, for example, lmfit, scipy.optimize.root, etc,(usually based on Newton-Raphson methods) will be used to find the root more efficiently. If more than one root exist in this region, it is necessary to divide this region into some subregions and check the number of roots in each subregion. This division process will be repeated until only one root remains in one region. Then one can go back to the second step to find the root.

This method asks the differentiation of the function f with respect to the argument z, which in this special problem, means the differentiation of the full linear dispersion relation ǫ (ω) + g (ω) (i.e, equation(B.6)) with respect to the frequency ω. This calculation is presented in appendix B, where the differentiation of the Bessel function of the first kind of integral order 0 is used: J ′ 0 (x) = -J 1 (x), with J 1 (x), the Bessel function of the first kind of integral order 1. The numerical implementation is realized in Python. Quadrature method is used for the energy integration in the dispersion relation. The division of the rectangular domain and the verification of the number of root are written in a separate module. After the resulting rectangular region where only one root exists is identified, a method named 'dogbox' [START_REF] Voglis | A rectangular trust region dogleg approach for unconstrained and bound constrained nonlinear optimization[END_REF] is used to find the root. Note that the 'dogbox' module can find the root in a given region, which is very suitable for the linear dispersion relation.

Since TIM and TEM instability are distinguished by the sign of the real part of the frequency, initially the complex plane of the frequency ω = ω r + i γ can be divided into two rectangular TEM TIM TEM TIM Figure 3.5: Linear growth rate as a function of the wave vector for two different temperature gradients κ T = 0.25 and κ T = 0.5 for TIM(red) and TEM(blue). One can see that the growth rate for the higher temperature gradient case is almost an order of magnitude higher than the lower gradient case. domains to find the eigenfrequency for TIM and TEM separately: one domain is defined as ω r > 0 and another one is defined as ω r < 0.

In order to verify this numerical method, it is tested for some randomly selected wave numbers and the result is shown in figure 3.4. In this example the left figure presents the test for the wave number k = k α = k ψ = 20, while the left one presents the test for the wave number k = k α = k ψ = 100. The values of the plasma dielectric function ǫ (ω) in the complex plane ω = ω r + i γ has been calculated and the points where the imaginary part ℑ [ǫ (ω)] and the real part ℜ [ǫ (ω)] equal to zero have been disclosed and marked by the green line and the black line, respectively, in both figures. So the regions where the two lines cross should be the regions where the roots exist. The stars (for TEM) and squares (for TIM) are the roots found by the algorithm. It has been shown clearly that the roots found by the algorithm fall exactly in the region where the lines crosses, which means this root finding algorithm resolves very well the linear dispersion relation.

These figures present numerically the linear dispersion relation as a multi-root problem. Comparing the left one to the right one, it can be seen that there are many damped modes in small scale, i.e. especially for TIM. In the region that is close to ℜ (ω) ∼ 0, a discontinuity in the linear dispersion relation has been found, which is in fact branch cut. In order to avoid the problem of the integration in the discontinuity region, it is necessary and useful to define a small TEM TEM TEM Figure 3.6: Growth rate for TIM(red) and TEM(blue) as a function of the wave vector for three different values of the trapped particle ratio: f t = 1/3, f t = 1/2 and f t = 2/3. It can be seen that, within the validity of the model, increasing f t increases the growth rate. distance (∼ 0.01) between the edge of the initial rectangle and the axis ℜ (ω) ∼ 0.

Linear instability for isotropy system

: γ (k) with k = k α = k ψ
In the previous sections, the critical temperature gradient κ T at the linear stability threshold has been studied with respect to ion to electron temperature ratio τ, trapped particle ratio f t and the wave number k. However, the threshold by itself does not allow to characterize the linear response of the plasma completely : the strength of the instability, which determines the energy injection for the nonlinear cascade model, indeed requires the knowledge of the linear growth rate γ. This section is contributed to resolve the linear dispersion relation and obtain the growth rate by using the developed numerical algorithm.

In the nonlinear parts of this thesis (i.e. chapter IV and chapter V), two different models will be discussed, corresponding to the isotropy assumption (i.e. k α = k ψ = k n , chapter IV ) and the anisotropy assumption (i.e. k = k ψ , k α = (k n , θ k ), chapter V). So linear instability of these two different cases should be discussed. The linear phase of the isotropy assumption will be discussed first, which can identify the effect of the free physical parameters in linear system. The TIM and TEM linear growth rate as a function of the wave number will be presented as follows for different temperature gradient κ T , different ion to electron temperature ratio τ and different trapped particle ratio f t .

In figure 3.5, the linear growth rate is shown as a function of the wave vector k n , for two values of the temperature gradient κ T . The other parameters used are the same as those given in table 3.1. For these parameters, the TEM linear growth rate is observed to be higher than the TIM one. The maximal growth rates can be observed to have higher values of the wave vector in the case of TEM than TIM, due to the difference between the trapped electron and trapped ion characteristic length scales. More importantly, increasing the temperature gradient increases the linear growth rates of the modes strongly, as discussed previously in the literature for TIM/TEM as well as ion temperature gradient driven mode (ITG) [START_REF] Dimits | Comparisons and physics basis of tokamak transport models and turbulence simulations[END_REF].

The effect of the trapped particles ratio has been investigated in figure 3.6, where the growth rate γ is given as a function of the wave vector k n for three different values of the trapped particles ratio f t , for both TIM (red lines) and TEM (blue lines). The rest of the parameters are given in table 3.1. Here the TEM present a stronger growth rate than TIM for the three values of the trapped particles ratio considered. It can also be noticed from figure 3.6 that increasing the ratio of trapped particles, the linear growth rate as well as the wavenumber of the most unstable mode is found to increase consequently. A last observation concerns the case f t = 1/3, which are found almost completely stable.

The effect of the ion to electron temperature ratio has been advertised in figure (3.7), which in fact is used to validate the numerical schemes of the nonlinear model with three values of τ. However it can present clearly the role of τ in linear phase. In figure 3.7, the growth rates obtained with the IVP code (black disks) are compared with the values given by the solution of the dispersion relation, for three values of τ, while the other parameters are kept constant as given in table 3.1, and κ T = 0.25. Note that numerical solution of the dispersion relation can make the distinction between TIM and TEM branches : the associated growth rates are respectively represented by red and blue crosses in figure 3.7, while the IVP can only show the maximum of them. Here three values for the parameter τ ∈ {1.0 , 0.6 , 0.2} are chosen in order to verify that the method is valid for TEM dominated, TIM dominated, as well as TEM/TIM mixed regimes. Figure (3.7) shows that the linearly dominating mode can be changed by changing the ion to electron temperature ratio.

Linear instability for anisotropic system: γ k ψ , k α

Since the linear terms in Vlasov equation depends on k α , which is defined as k α = k n sin (θ k ), the linear system should also have a dependence on the wave-vector angle θ k , which means that it is highly anisotropic. This section will be contributed to the discussion of the linear instability and the issue of anisotropy.

Since the eigenfrequency ω in linear dispersion relation (i.e., equation(3.4)) is normalized by k α ≡ k n sin(θ k ), it should be symmetric with respect to π 2 , which means that we only need to solve the eigenvalue problem in the plane θ ∈ 0; π 2 and the results in the rest plane θ ∈ π 2 ; 2π can be obtained automatically by the symmetry. The eigenfrequency of the linear dispersion relation is resolved by the algorithm based on the argument principle introduced in section(3.2).

The linear growth rate for three different cases will be shown: adiabatic electron, adiabatic ion and full kinetic system. They are also the cases that will be discussed in nonlinear simulation. In gyrokinetic study adiabatic electron assumption is widely used to decrease the work of Figure 3.7: Comparison between the growth rates are shown as functions of the wave vectors k n obtained by solving the linear dispersion relation (red and blue crosses for respectively TIM and TEM) and from the initial value problem (IVP, black squares). Three values for the temperature ratio, τ = 1.0, τ = 0.6, and τ = 0.2 (with f t = 2/3 and κ T = 0.25 in all cases) were considered. These three cases, which are presented from top to bottom for decreasing τ, correspond to TEM dominated, mixed TEM/TIM, and TIM dominated regimes respectively. simulation. In this case only ions are kinetic and only trapped ion mode(TIM) exist in the system. However in trapped particle problem TEM is usually much more important than TIM, so it would also be interesting and useful to solve the system for kinetic electron. The linear growth rate of kinetic ion and electron will be shown to compare the difference between the fully kinetic system and the system with one adiabatic species. In order to discuss the adiabatic electron(ion), one needs to set the parameter f t of electron(ion) to be zero, so only ions(electrons) are trapped and kinetic and finally there is only TIM(TEM) in the linear system. In figure (3.8) the TEM and TIM growth rates are presented and the parameters used in these two simulations are given in table (3.2), where the upper set of parameters is used in the simulation of adiabatic electron and the bottom one is used in the simulation of adiabatic ion. One can see that all the other parameters are the same except for the Larmor radius ρ, the banana width δ, and the trapped particle ratio f t .

ρ i δ i f t ,e f t ,i τ a/R 0 Ω d κ T 0.03 0.
From figure (3.8) one can see that in a given direction (i.e θ k fixed), the growth rate γ increases with the wave modulus k n , when it arrives to the maximum value, it begins to fall due to Landau damping. The maximum growth rate γ max in the full wave vector plane is located in the direction θ k = π 2 , which corresponds to a "streamer", and it is always zero in the direction θ k = 0, which is the zonal flow direction. So in linear system streamers are the most unstable modes Figure 3.9: the TIM (left) and the TEM (right) growth rate in the wave plane k = k ψ , k α for the full kinetic system (i.e, kinetic ions + kinetic electrons). and zonal flows are stable. These features are similar for both TIM and TEM.

The white region presents the modes with γ = 0, so the inside zone (i.e.γ > 0) is unstable and the outside zone (i.e.γ < 0) is stable. Comparing the left figure and the right figure, one can find that the unstable zone of TEM is much larger than the unstable zone of TIM and the unstable zone of TIM is like an ellipse, but the unstable zone of TEM is close to a circle. This result is caused by the difference between the Larmor radius ρ s and the banana width δ s : if ρ s = δ s , the unstable zone is a circle, otherwise it is always ellipse. We should mention that for real physical parameters the banana width δ s is always much larger than the Larmor radius ρ s , so the result should be more like a ellipse in real physics. The maximum growth rate of TEM is γ T E M max ∼ 3.5, but the maximum value of TIM is γ T I M max ∼ 1.5: it is the main reason that drive us to discuss the unusual adiabatic ion (TEM) case in nonlinear simulation.

The growth rate of the full kinetic system is shown in figure(3.9), since the ions and electrons are both kinetic, so the TEM and TIM co-exist in the system. The left figure shows the TIM linear growth rate and the right figure shows the TEM linear growth rate. Comparing with the adiabatic system, their form doesn't change much and the strongest mode is located around the same region, but the growth rate is much more stronger than the adiabatic system for both TIM and TEM, which means the full kinetic system is much more unstable than the system with one adiabatic species. For example the maximum growth rate of the full kinetic system are γ max T I M ∼ 3.5 and γ max T E M ∼ 7 , however the maximum growth rate of the adiabatic system are: γ max T I M ∼ 1.5 and γ max T E M ∼ 3.5 . The linear energy injection of the full kinetic system is much stronger than the adiabatic system.

Conclusion

In this chapter, the linear system has been studied. The linear dispersion relation has been calculated from the linear Vlasov and Poisson equations by assuming a plane wave solution of the particle distribution function. Due to the singularity along the real axis of the frequency, the residue in the lower half plane of the complex frequency has been added to the linear dispersion based on the residue definition in Fourier-Laplace transform. Since the TEM and TIM instabilities co-exist in the system and the linear dispersion relation is a multiple roots equation, it is necessary to find the eigenfrequency for both TEM and TIM. However a system with singularity and multiple roots is a challenge for most of the existing root finders, an algorithm based on the argument principle has been developed to calculate the eigenfrequency of the linear dispersion relation.

In this model, there are two sources of instability: the temperature gradient and the density gradient, the threshold for the external parameters and the sources of instability have been discussed before solving the eigenvalue problem, which allows a direct understanding of the effect of these parameters, such as the wave number k,the ion to electron temperature ratio τ and the trapped particle ratio f t .

After having analyzed into details the dependence of the linear stability threshold with respect to the external parameters τ, κ T and f t , a characterization of the linear response in terms of the growth rate has been carried out, with the help of an eigenvalue solver based on the argument principle.

Again, the influence of the three main parameters τ, κ T and f t has been studied with respect to the linear growth rates associated to TIM and TEM, showing that the linear system becomes more unstable with higher temperature gradient (increasing κ T ) and more trapped (kinetic) particles (increasing f t ). The linearly dominated mode can be changed by varying the ion to electron temperature ratio τ.

A further study shows that the linear energy injection in the wave number plane k = k ψ , k α is strongly anisotropic: the zonal flows (i.e. k α = 0) are always stable and the most unstable mode is always located at k ψ = 0, which is a streamer. Comparing the linear growth rate of three different cases: kinetic ions (+adiabatic electrons), kinetic electrons (+adiabatic ions) and fully kinetic system, it has been observed that under the same external parameters, the fully kinetic system is much more unstable than the other two. The nonlinear simulations about these cases will be discussed in chapter V.

Chapter 4 Isotropic model: Sabra and GOY

Before the direct nonlinear simulation of the Sabra (2.49) and GOY (2.48) models, the linear phase, such as the stability threshold of the external parameters, the linear instability, etc, have been discussed in chapter III, which gives a linear picture of the trapped particle model. However the complex, turbulent behavior originates from the nonlinearity of the equations, which is much more interesting and complicated compared to the linear phase. The nonlinear term is defined by the Poisson bracket in real space, which after the Fourier transform, becomes a convolution in wave number plane.

In chapter II, in order to simplify the nonlinear structure of turbulence problem, a description of nonlinear interactions has been worked out, inspired by the fact that the wave number is discretized in simulation such that the number of wave as well as the number of interaction couldn't be infinite. To study the physics at small scales, a large range of wave numbers is necessary. Thus a logarithmic grid is used to discretize the wave number space. This idea is mainly inspired by LDM [START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF] that has been recently developed to resolve the passive scalar equation in fluid turbulence and this formulation (2.42) gives a complete description of the nonlinear interactions, which is the main difference compared to previous shell models, so it can be seen as a development or a progress of the shell models in the field of turbulence.

However in practice, limited by the computer resource, usually only a subset of all the nonlinear couplings is included in nonlinear simulation. In kinetic plasma turbulence, in addition to the spatial coordinate, the kinetic coordinates also exist in the equation, so the dimensionality is higher than that in fluid turbulence, which means more computer resources are required to resolve the kinetic plasma system than that in fluid system. Therefore as a first application of shell model to kinetic plasma turbulence, we will follow the approach from fluid turbulence, which means only a subset of all the couplings will be considered in nonlinear simulation. Even in this simplified case, it still takes some days to finish the simulation of the anisotropic model.

The nonlinear phase concerns the simulation of the isotropic (k = k α = k ψ ) model and the anisotropic model (k = k ψ , k α ). This chapter starts with the simplest one, where the dimension of the wave number space is reduced by imposing the isotropy assumption. Together with the kinetic coordinate, it is a 2D ODE system, which takes around 20 hours to finish the simulation (not parallelized ). This chapter is organized as follows: in section (4.1), the nonlinear Vlasov equation of the isotropic model will be presented, where a reduced analytic expression of the nonlinear terms (i.e. Poisson Bracket) as well as the calculation of the nonlinear interaction coefficients based on the conservation of the electrostatic potential energy and the entropy will be detailed. The phase problem that is very general in the 1D shell model will be explained in the beginning of section (4.1) and based on different phase approximations, two different versions called GOY model and Sabra model will be given. In section (4.2), the results of nonlinear simulation of Sabra model, such as the k spectra of electrostatic potential energy and of entropy, the effect of free parameters as well as the kinetic effects that are especially linked to the kinetic plasma turbulence will be presented. In section(4.3), the nonlinear simulation of the GOY model will presented, where the role of the phase in shell models will be better understood by a comparison between the models (GOY vs Sabra) based on different phase approximations.

Model equations

In shell models [START_REF] Vov | Improved shell model of turbulence[END_REF][START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF][START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF] the wave number can be divided into shells with k = k n = k 0 g n , where g is the logarithmic spacing factor and n is the shell number. Shell model, which uses a logarithmic spacing between wave numbers can be used to describe a very large interval in the wave number space, something that is virtually impossible with linear discretization.

In the description of the nonlinear interaction in chapter II, a formulation of the complete interaction has been given, where the nonlinear couplings can be obtained easily by a program. For example in a range of wave number k n = k 0 g n with g = 1.56, k 0 = g -2 , and N = 40, for each wave number, there are around 42 ∼ 45 nonlinear triads, note that the number of triads are a little bit different for different wave numbers, because of the wave numbers in the boundaries (beginning and ending). Even though this description is very simple and can present the real nonlinear physics in some sense, the numerical implementation of such a nonlinear system is still very expensive and the physics with the small scales should be treated carefully. As a first application of the shell models to the kinetic plasma turbulence, we will start with the simplest one, where the isotropy assumption is imposed to the system and only the local interactions will be taken into account. One of the simplest choices that achieve such a reduction used in GOY [START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF] or Sabra [START_REF] Vov | Improved shell model of turbulence[END_REF] shell models is the restriction of the interactions to the closest neighbors of the n t h shell, which means the following triads:

{(n -2, n -1, n) , (n -1, n, n + 1) , (n, n + 1, n + 2)} .

Phase approximation

When the system pass to 1D, the phase information is lost. So the nonlinear Vlasov equation (2.30) in wave number space can be written in different forms, based on different phase approximation. For example the two widely used cases : Sabra and GOY, where the nonlinear expressions are given respectively as follows:

f , f sabr a n = α f n-2 f n-1 + β f * n-1 f n+1 + γ f * n+1 f n+2 , (4.1) f , f GOY n = α f * n-2 f * n-1 + β f * n-1 f * n+1 + γ f * n+1 f * n+2 , (4.2)
where f , f n means the nonlinear terms of the n t h shell. α, β and γ are the interaction coefficients corresponding to the triad (n -2, n -1, n) , (n -1, n, n + 1) and (n, n + 1, n + 2) , respectively. Since the shell models in general does not really give reliable information about the complex phases of the fields, because they only rely on the conservation of absolute quadratic quantities, which are real, so the phase information is not fixed in these models, one only need to satisfy the following relation of the phase:

θ n-2 + θ n-1 + θ n = θ n-1 + θ n+1 + θ n = θ n+1 + θ n+2 + θ n = 0 ,
where θ n is the phase for the n t h shell. Based on this relation, different phase approximations exist. In a 1D coordinate, considering the direction of the axis, only two choices of the phase exist: 0 or π, corresponding to the positive direction and the negative direction of the axis, respectively. So the wave number k in a 1D reference can exist in two states: k (phase θ k = 0) or -k(θ k = π), where the corresponding physics variables are f k and f * k , with f * k , the conjugate of f k , due to the fact that the variable φ and f are real quantities, so f -k = f * k . In GOY model, all the phase information roughly equals to 0, which means the direction of the wave is neglected. Sabra model treats the phase information a bit better. In order to satisfy rigorously the relation k +p +q = 0, the direction of the wave must be taken into account. In the triad (n -2, n -1, n), since k n-2 , k n-1 < k n , the modes k n-2 and k n-1 must be positive. However in the triad (n -1, n + 1, n), since k n+1 > k n , the mode k n-1 must be negative, such that k n+1k n-1 = k n , which means that f for the mode k n-1 should appear as conjugate. Using the same logic to analyze the triad (n + 1, n + 2, n), it can be found that f n+1 should appear as conjugate.

Model equations

Based on these two phase approximations, the expression of the nonlinear term can be written in two forms, which result in two different models, called GOY and Sabra, respectively. The discussion of the nonlinear simulation begins with the Sabra model. Applying the Sabra truncation to the trapped particle model, the nonlinear terms can be written as follows:

J 0s φ, f s n = αJ n-2 0s φ n-2 f s n-1 + α ′ J n-1 0s φ n-1 f s n-2 +βJ n-1 0s φ * n-1 f s n+1 + β ′ J n+1 0s φ n+1 f s * n-1 +γJ n+1 0s φ * n+1 f s n+2 + γ ′ J n+2 0s φ n+2 f s * n+1 . (4.
3)

The exchange of φ and f should be considered compared to the simple Sabra expression given in equation (4.1). Since in the triad p + q + k = 0, the evolution of f k can be determined by f p and φ q , it can be also determined by φ p and f q .

The coefficients α, β, γ α ′ , β ′ and γ ′ are the interaction coefficients, which can be calculated via the conservation laws. Note that the Poisson bracket has the general property that the phase space integration of its multiplication by any function of either of its arguments is zero. Since the shell model that we want to develop has to preserve the conservation properties of the original system, the following two relations have to be satisfied also by the shell model:

n ˆ∞ 0 J 0s φ n J 0s φ, f s n E d E = 0 , (4.4
)

n ˆ∞ 0 f s n J 0s φ, f s n E d E = 0 , (4.5)
which corresponds to the conservation of electrostatic potential energy and entropy respectively (Appendix E). In a very general manner, such conservation properties are used as a basis for spectral analysis of turbulence: energy transfer for trapped ion modes have, for instance been studied using fluid models [START_REF] Newman | Energy transfer dynamics of dissipative trapped ion convective cell turbulence[END_REF]. These constraints together with the quadratic k dependence of the Poisson bracket allow us to determine the coefficients α, β ,γ, etc in equation (4.6).

The detail of this calculation can be found in appendix (A), here we give directly the results as follows:

α = -α ′ =α 0 k 2 n g -3 , β = -β ′ = -α 0 k 2 n g -1 , γ = -γ ′ =α 0 k 2 n g ,
where α 0 is a free parameter (we will simply call α from this point on for simplicity), which remains unknown in general for shell models since the number of constraints is always lower than the number of free parameters. Note that the problem of free parameter is not present in the 2D model, where the phase in a triad can be exactly calculated [START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF]. Inserting the coefficients in equation (4.6), the final Vlasov equation is presented as follows:

∂ f s n ∂t =i k n J n 0s φ n ∂F s ψ ∂ψ -i E Ω d Z s k n f s n + αk 2 n g -3 J n-2 0s φ n-2 f s n-1 -J n-1 0s φ n-1 f s n-2 -αk 2 n g -1 J n-1 0s φ * n-1 f s n+1 -J n+1 0s φ n+1 f s * n-1 + αk 2 n g J n+1 0s φ * n+1 f s n+2 -J n+2 0s φ n+2 f s * n+1 , (4.6)
and the quasi-neutrality relation remains unchanged, since it is linear with respect to the unknowns:

c n φ n = ˆ J n 0i f i n -J n 0e f e n E d E . (4.7)
At this point the model equations are clear. Note that all the transformation, simplification and assumption, etc. are imposed on the spatial coordinates, the dependence on the kinetic coordinates, such as the integration over the energy in Poisson equation, is not modified, so the kinetic effect should be kept by our model.

Nonlinear simulation of the Sabra model

A nonlinear expression of the Vlasov equation of trapped particle model has been given explicitly in previous section, based on the phase approximation in Sabra truncation. The nonlinear simulation of this simple model will be discussed, where we will show some general results in turbulence study, such as the k spectra of the energy and the entropy, as well as some new results of the kinetic plasma due to the kinetic effect, for example the resonance.

This nonlinear Vlasov-Poisson system (4.6,4.7) is a complex IVP (initial value problem) ODE (ordinary differential equation) system. For such a general problem, there exists a set of ODE solvers, such as the Vode [START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF], Zvode, dopris5, dop853 [START_REF] Oliphant | Python for scientific computing[END_REF], etc. In this thesis we use dop853 to resolve the system since it is more accurate and well adapted to stiff problems after a comparison between the results obtained by different solvers. Before the nonlinear simulation, it is necessary to validate the numerical schemes.

Code verification

In order to validate the numerical schemes that are used in our model, we consider the linear limit by setting α = 0 (which sets nonlinear couplings to zero) and compare its results (which we call the initial value problem -IVP) against the growth rates that are obtained by solving the dispersion relation by the root finding algorithm previously discussed. To be certain that the code works well in any situation, three different values of the ion to electron temperature ratio τ, corresponding to different linear regimes, such as TEM dominated, TIM&TEM mixed and TIM dominated, are simulated in the code. In the initial value simulation, the Gauss Legendre quadrature is chosen to resolve the energy integration in Poisson equation, after a comparison against the trapezoidal method and the Gauss Laguerre quadrature (Appendix F).

In figure 4.1, the growth rates obtained with the IVP code (black disks) are compared with the values given by the solution of the linear dispersion relation, for three values of τ, while the other parameters are kept constant as given in table 3.1, and the temperature gradient is κ T = 0.25. Note that numerical solution of the dispersion relation can make the distinction between TIM and TEM branches : the associated growth rates are respectively represented by red and blue crosses in figure 4.1, while the IVP can only show the maximum. Here we choose three values for the parameter τ ∈ {1.0 , 0.6 , 0.2} in order to verify that the method is valid for TEM dominated, TIM dominated, as well as TEM/TIM mixed regimes. The comparisons given in figure (4.1) display very good agreement between the linear growth rates computed from the dispersion relation and the linear solutions of the model, since the difference between γ I V P and max γ T E M , γ T I M is less than 1%. This remains valid with positive growth rate as well as with negative growth rate. The small difference observed could be explained by the discrete energy grid chosen (E max = 20, N ℓ = 200), which may cause an error in the Landau damping mechanism. It is interesting to notice that by varying τ, the dominant instability changes from TEM with high values of τ, to TIM with low values of τ, which agrees well with the linear threshold study. Note that the maximum energy is determined by the convergence of the linear growth rate. Note that the linear growth rate of the whole system is always closed to max γ T E M , γ T I M , which can be explained by the simple qualitative analysis given in chapter 3 section 3.3.

k spectra of the electrostatic potential energy E φ

A reference nonlinear simulation is presented here, with the standard parameters given in table 4.1. Energy injection is ensured self-consistently by the linear terms in the Vlasov equation 4.6. This means that in order to reach a quasi-stationary turbulent state, dissipations have to be

E f i (k n ) = ´Cn | f i | 2 /F i /k n E d E and and electron entropy E f i (k n ) = ´Cn | f i | 2 /F i /k n E d E
as a function of k n , with a time average over the turbulent quasi- stationary phase for the parameters given in table 4.1. The blue line is the entropy spectra of electrons, while the red one is the entropy spectra of ions. The k -4/3 and k -1 present qualitatively the power law scalings of ion entropy and electron entropy. added to the model. Here we choose a form that localizes the dissipations in "buffer" regions in k-space using a hyper-viscosity of the form D = ν L k 6 + ν s k 4 , where large scale dissipation is controlled by the value of ν L , while the small scale dissipation is given by ν s . Note that in the nonlinear simulations of this chapter, τ = 1.0 and f t = 2/3 and the other parameters are presented in table 3.1.

As the principal result of the reference simulation, the wave number spectrum of the electrostatic potential energy is presented in figure 4.2. Note that it is the electrostatic potential energy in the shell E φn = C n |φ n | 2 which is conserved by (4.6), which means that the spectral electrostatic energy density (equivalent of E (k), the spectral energy density in Kolmogorov theory) can be written as

E φ (k n ) = C n |φ n | 2 /k n so that we have ´Eφ (k) d k = n E φn , with d k = k n g d n in shell N ℓ E max N g k 0 α κ T ν s ν L τ 200 20 30
E f i = | f i | 2 F i
averaged over the turbulent phase as a function of energy E which is linear, and the wave number k which is logarithmic. Note that the effect of dissipation is most evident at high k and large E . model.

In the spectrum of electrostatic potential energy, a k -4 (a linear regression gives 4.08) power law scaling has been observed in the range where nonlinear transfer dominate the system: the so called inertial range in Kolmogorov's notation. The reason that k -4 slope is different from the classical result k -3 is probably due to the existence of Bessel functions in the nonlinear term, since

Γ 0 k 2 ⊥ ∼ k -1
, where Γ 0 is the exponentially scaled modified Bessel function of first kind. Comparing the spectrum of different large scale dissipations, it can be found that a stronger large scale dissipation results in a lower electrostatic potential energy in large scale. It should be noticed that in the turbulence research, the power law is usually assumed to be determined by the form of the nonlinear term, since the analytic result of the power law scaling can be calculated directly from the nonlinear expression, if the nonlinear expression is not very complicated. While the difficulties , such as the Bessel function and the energy integration complicate the analytic work here. The dissipation and the external parameters have no qualitative effect on the power law scaling, which has been demonstrated many times by simulations.

k spectra of the entropy E f and kinetic effect

In this model, the entropy

E f (k n ) = ´| f n | 2 /F i /k n E d E
is another nonlinearly conserved quan- tity. The k spectra of the entropy for ions as well as for electrons are shown as follows:

In figure (4.3) the k spectra of the entropy of ions and the entropy of electrons are presented. It can be seen that the differences of the ion entropy and electron entropy are mainly located at very large and very small scales. In the range where the linear injection and the dissipations are negligible, the difference between E f i and E f e is very small. we give two different power laws scalings k -4/3 and k -1 respectively for ions and electrons. Note that the small difference between the power law scalings is due to the behavior in the beginning and ending of the inertial range.

While the model that we propose, reduces the complexity of the nonlinear problem by the use of truncation based on shells in k-space, the equations that are solved are still kinetic as the electron and ion distribution functions depend on the energy coordinate E ℓ , and the kinetic dependence of this system is not simplified, so this model should keep the kinetic effect, such as the precession resonance. To present the kinetic phenomenon of this model, the free energy for ions 4.4 as a function of energy E (linear scale) and the wave number k n (log scale), averaged in time over the turbulent phase.

E f i = | f i nl | 2 /F i is shown (log scale) in figure
It can be observed from figure 4.4, that the ion free energy is rather high for small energy E < 10, while the contribution for E > 20 is relatively negligible. This phenomenon can be used to reduce the size of the energy grids, which will be very useful in the anisotropic 3D (2D spectral +1D energy) model simulation. The effect of small scale dissipation is mostly visible at high k n and large values of E . A sudden drop in free energy can be observed around E ∼ 3/2 for all k. This is caused by the resonance ω r /k = E /T s ∼ 3/2, as given in equations (3.9, 3.10), which somehow persists in the nonlinear phase. Note that, it is unclear whether the loss of resonance at high-k is due to the dispersion relation deviating from linear or due to the effects of dissipation. Here we only observe the resonance around E ∼ 3/2, because only one advection frequency Ω d associated to the toroidal precession motion is chosen.

The result in figure (4.4) shows clearly that this is a kinetic shell model (also the first shell model of kinetic plasma) and the kinetic resonance is well observed, which means the shell models developed in fluid turbulence can be applied to kinetic plasma turbulence.

The effect of free parameter α

In shell models, a free parameter α remains in the equation, which can not be determined from physical parameters. This means that the results that we obtain, should be qualitatively insensitive to the value of α. In fact, since it is a factor in front of the nonlinear term, it can be seen as a ratio (or magnitude) of the linear effect compared to the nonlinear effect. One can guess that the saturation level will be determined by a balance between α and γ, the growth rate. In order to verify this idea, simulations for different values of α are presented in figure 4.5. Parameters that are used in these simulations are the same with the previous simulation except that we used ν L = 1.0. It is seen in figure 4.5 that the spectra are parallel to each other and follow exactly the same trends, so that the slopes of the spectra in the range of interest (i.e. equiva-Figure 4.5: Test for free parameter α. The spectra from top to bottom correspond to values of α equal to 1, g 2 , 10 and 20 respectively. It can be seen that the spectra are perfectly parallel to each other. lent to the inertial range in Kolmogorov's theory) do not change with α, which is a key feature of this model. These simulations also show that the saturation level increase with decreasing α as expected from a basic mixing length argument of the type γφ = αφ 2 ⇒ φ 2 ∝ γ α

2 . In the k-range where the production and dissipation rates are small compared to the nonlinear transfer (for k n ∈ [56 ; 3000], while k n < 56 corresponds to the injection range, and k n > 3000 is the dissipation range) the slope obtained by a linear regression is -4.08, suggesting a k -4 power law scaling, which is indicated in figure 4.5 by a black line.

Influence of phase information: GOY vs Sabra

In principle one has the freedom to choose different combinations of phases that nonetheless respect the conservation laws of the original system, as explained in the section on phase approximation. Therefore, we have developed two different models: Sabra version (equation(4.6)) and GOY version (equation(4.8)). The nonlinear results of Sabra version has been discussed in previous sections, this section is contributed to the discussion of the nonlinear simulation of GOY version as well as the comparison and the difference between the models based on different phase approximation.

comparison of Sabra and GOY GOY Sabra , which is a value of k n a little higher than the value associated to the maximum in the kspectrum. Interestingly, while the GOY model shows distinct, large amplitude oscillations, the Sabra model, seems to evolve more randomly.

Oscillatory dynamics of the GOY Model

The GOY version is obtained by using complex conjugates of the fields in all the nonlinear terms appearing in equation (4.6). Applied to the trapped particle models, the Vlasov equation in GOY version is written as follows:

∂ f s n ∂t =i k n J n 0s φ n ∂F s ψ ∂ψ -i E Ω d Z s k n f s n + αk 2 n g -3 J n-2 0s φ * n-2 f s * n-1 -J n-1 0s φ * n-1 f s * n-2 -αk 2 n g -1 J n-1 0s φ * n-1 f s * n+1 -J n+1 0s φ * n+1 f s * n-1 + αk 2 n g J n+1 0s φ * n+1 f s * n+2 -J n+2 0s φ * n+2 f s * n+1 . (4.8)
It can be found that the prefactors of the nonlinear terms in GOY model are the same as these in Sabra model. As explained in chapter 2 section (2.3) the derivation of the model, the interaction coefficients equal the surface of the triads, since the triads considered in these two models are same, the interactions coefficients should be the same. In order to understand the influence of the phase approximation, we define all the other parameters , such as the external physical parameters, the wave numbers as well as the free parameter α (i.e, α = 20 for both) and the initial conditions to be exactly the same as those used in the simulation of Sabra model.

In figure (4.6) the time evolution of the absolute value of electrostatic potential φ for mode n = 6 are presented for GOY version and Sabra version. In the linear phase (t = 0 ∼ 5), these two models follow exactly the same curves, since linear terms dominate this phase and they are the same in both models. However, when the system reaches saturation, the nonlinear terms become important, and the difference between the forms of the nonlinear terms in the two models start to play a role. For the mode n = 6, with the same parameters, the GOY version saturates at |φ| ∼ 10, while the Sabra saturates around 0.1. One may guess that the nonlinear effect in Sabra model is more important than that in GOY model, because the linear terms are really the same for both. The temporal evolution of the Sabra model is quite smooth compared to the GOY model (which is true also for other large scale modes). However except these differences, the power law scalings of the electrostatic potential energy are the same for both GOY and Sabra model.

The periodic oscillation that is observed in the turbulent phase for the GOY model is a characteristic feature specific to this model. Different parameters were tested in the simulations of the Sabra model, but no such regular oscillation was ever observed. When the oscillations of different shells are compared with n = 6, which is taken as the reference shell, one sees a little phase shift as a function of the shell number n (suggesting a propagation towards higher k). For some special shells (like n = 12 and n = 18 ) the phase shift is actually about π/2, it means the maximum electrostatic potential φ for mode n = 12 corresponds to the minimum value for mode n = 18, as shown in figure (4.7). This kind of dynamics may present a similarity with the predator prey phenomenon, well known between zonal flows and drift waves in fusion plasma [START_REF] Malkov | On the nature of bursting in transport and turbulence in drift wave-zonal flow systems[END_REF][START_REF] Kim | Zonal flows and transient dynamics of the L-H transition[END_REF], which is believed to play a key role in particular in the low to high confinement transition in Tokamaks [START_REF] Miki | Spatio-temporal evolution of the l -> i -> h transition[END_REF]. Note that here since we do not have a proper treatment of zonal flows, the largest scale may end up playing the role of the predator, while the small scales may play the role of the prey [START_REF] Berionni | Predator prey oscillations in a simple cascade model of drift wave turbulence[END_REF][START_REF] Morel | Characterization of predator-prey dynamics, using the evolution of free energy in plasma turbulence[END_REF].

Comparison of entropy

The entropy of these two models are compared in figure (4.8). Due to the different saturation mechanisms, the saturation level of entropy is very different between GOY and Sabra model. For the parameters used here, the maximum values of the entropy in GOY model is around 1.0, while the maximum value of the entropy in Sabra model is only 1.0e -3: the difference is quite big, which means that the phase approximation has a very important influence into the final saturation level. In the k range, the spectra of the same species (i.e, ions and electrons) are almost parallel to each others, which means that the power law scalings should be similar for the same species in both case. Here a k -4/3 scaling is given in the figure, which is valid for the entropy of ions. As shown in figure(4.3), the k -1 is a better scaling for the entropy of electrons. It can be also found that a little bit difference exists in the entropy spectrum of ions and electrons, which produce the different power law scaling (i.e, k -4/3 vs k -1 ). In detail this difference appears in the inertial range (∼ 1.0e2) in GOY model, where the entropy of ions is a little bit higher than the entropy of electrons, however in Sabra model this difference appears at very large scale (< 10), where the entropy of electrons is higher than the entropy of ions. In the inertial range of Sabra model, the entropy of ions and electrons is almost the same, as shown in figure(4.8).

Figure 4.9: The two phases of the oscillations of the electrostatic potential energy spectrum of GOY. Here the blue line represents the mean spectrum, while the red one represents the spectrum around the maximum of mode n = 6 , and the green line represents the spectrum around the minimum. It can be seen that the oscillations are roughly out of phase between the large and the small scale regions.

Oscillation of the k-spectra

Interestingly, some qualitative difference have been observed in the behavior of the simulations of the GOY version of our model with respect to the reference Sabra version. In order to compare these two models, the same parameters are used as in table 4.1 but with ν L = 2.0 as in the case above.

In fact, a shorter time averaging of the wave number spectrum, where the average can be computed around the maximum and the minimum of one oscillation can be used in the turbulent phase to see the effect of these oscillations on different wave numbers. Here the time windows that define the minima and the maxima of the oscillations are defined with respect to mode n = 6, which is used as reference. The results are presented in figure 4.9, which shows that the oscillations are roughly out of phase between the large and the small scale regions.

In the standard predator prey picture, assuming that the growth rate is constant, the frequency of the oscillation is determined by the large scale dissipation. In order to understand how the oscillation in the GOY version of our model varies with large scale dissipation we perform a scan of ν L ∈ [2.0, • • • , 10.0]. Simulation result shows that the oscillation amplitude of the mode does not change much for different ν L , while the oscillation frequency changes. Using fast Fourier transform in the turbulent phase, the nonlinear oscillation frequency ω can be de- termined. In fact, it was observed that in the turbulent phase, all the modes oscillate with the same frequency ω, which can be shown to vary by changing the large scale dissipation ν L as shown in figure 4.10.

Conclusion

The first nonlinear simulation of the kinetic shell model is presented in this chapter, in order to discussed the trapped particle turbulence in Tokamak. Here we use the isotropic assumption, which is very general in turbulence study, to reduce the system, as well as to simplify the simulation. Based on different phase approximations, two models, called GOY and Sabra, equivalent respectively to the phase of GOY and Sabra models in fluid turbulence, have been derived, where the Sabra model has been discussed first, since its phase approximation is more rigorous. The k spectrum of the conserved quantities is commonly discussed in turbulence study, which in our system refer to the electrostatic potential energy density and the free energy of the particles (i.e.entropy). For the spectrum of the electrostatic potential energy, the simulations shows a k -4 power law scaling (linear regression gives 4.08), which is similar for both Sabra and GOY models.

The entropy spectrum appears to be a slightly different for ions and electrons in both models: the ion entropy spectrum shows a k -4/3 power law scaling against the k -1 scaling of electron entropy. Another important result that is especially associated with the kinetic plasma system is the kinetic resonance, which is very well observed in this work. This means that the shell models that have been developed for the fluid turbulence, can be applied to kinetic plasma turbulence and can be used to describe the kinetic effects. The simulations concerning a scan of α, present that the k-spectra is not changed qualitatively by this free parameter. But a larger α results in a lower saturation level, which can be qualitatively explained by a balance between the linear and nonlinear terms. The temporal evolution of the electrostatic potential in GOY model shows oscillating dynamics, which can not be observed in Sabra model. A further study shows that the oscillation frequency in GOY model is influenced by the large scale dissipation. Moreover the entropy in GOY model is observed much higher than that in Sabra model under the same parameters, even though their entropy spectra are almost parallel (i.e. same power law scaling).

The linear physics is anisotropic for kinetic plasma. In the linear study (i.e., chapter III), the TEM and TIM linear instability shows a strong anisotropy in the wave plane: the zonal modes are always stable and the streamers are the most unstable. The comparison between the isotropic models based on different phase approximations (i.e., GOY vs Sabra) presents that the phase information may change the nonlinear dynamics. Thus, a model with exact phase information is necessary to be worked out for the bounce averaged gyrokinetics. All these differences, caused by the phase approximations, as well as the fact that the linear physics is anisotropic as shown in chapter III, motivate us to work on the nonlinear anisotropic model, that we present in the next chapter.

Chapter 5 Anisotropic model: LDM

Logarithmically discretized model (LDM) has been recently developed in the context of fluid turbulence [START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF], which allows keeping track of exact phase information of each triad and can handle anisotropy. The basic formulation of LDM can be applied to the bounce averaged gyrokinetics, since the nonlinear term of the passive scalar equation is very similar to the 2D Poisson bracket in Fourier space if gyro bounce average operator (Bessel function) and the kinetic dependences are neglected. Since the spatial coordinates and the kinetic coordinates are independent, using the idea of LDM to simplify the nonlinear description on the wave number plane k = k ψ , k α while keeping the kinetic part as it is, a logarithmic discretization model of bounce averaged gyrokinetics has been obtained. After the introduction of the modelling details of this anisotropic model, the results of nonlinear simulations will be discussed in this chapter.

LDM of bounce averaged gyrokinetics

In this section we introduce the LDM grid and the notations associated and we recall the expressions of Vlasov and Poisson equations under the LDM notation that has been developed in chapter II section (2.3).

LDM Grid

The LDM grid is based on a polar representation of the Fourier space k = (k, θ k ), where a logarithmic grid is used with respect to the wave-vector modulus k n = k 0 g n , n = 0, ..., N -1, where g is the logarithmic spacing factor and n is the shell number. The angles θ k are discretized onto a regular grid θ j ≡ j 2π M , with j = 0, ..., M -1. The logarithmic grid with respect to the wave-vector modulus allows a clear link with shell models [START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF], the denomination of "shell" will be used in the following. The log polar coordinates system of LDM is presented in figure(5.1).

Vlasov-Poisson equation

The Vlasov equation and Poisson equation of trapped particle model are introduced in section (2.1). The nonlinear Vlasov equation based on the LDM grid has been given in section (2.3). We will give more details about this nonlinear Vlasov-Poisson system before discussing the results of nonlinear simulations.

The nonlinear interactions are located between the nearest neighboring modes as represented in figure (5.1), which imply in terms of triads:

k, p, q ∈ (n, n -1, n -2) ; (n, n -1, n + 1) ; (n, n + 1, n + 2) .
(5.1)

Given the locality assumption, the derivation of the LDM from Fourier transformed trapped particles model is straightforward, the details and discussions about the derivation have been given in chapter II. The Vlasov LDM equation describing trapped particles dynamics can be given in an explicit analytic form as follows:

∂ f j s,n ∂t = i k α J 0s φ j n ∂F s ∂ψ - E Ω d Z s i k α f j s,n + 1 2 k 2 n g -4 µ 0 J 0s φ * j +r 0 n-2 f * j -s 0 s,n-1 -J 0s φ * j -s 0 n-1 f * j +r 0 s,n-2 + J 0s φ * j +s 0 n-1 f * j -r 0 s,n-2 -J 0s φ * j -r 0 n-2 f * j +s 0 s,n-1 + 1 2 k 2 n g -2 µ 0 J 0s φ * j +l 0 n-1 f * j -s 0 s,n+1 -J 0s φ * j -s 0 n+1 f * j +l 0 s,n-1 + J 0s φ * j +s 0 n+1 f * j -l 0 s,n-1 -J 0s φ * j -l 0 n-1 f * j +s 0 s,n+1 + 1 2 k 2 n µ 0 J 0s φ * j +l 0 n+1 f * j -r 0 s,n+2 -J 0s φ * j -r 0 n+2 f * j +l 0 s,n+1 + J 0s φ * j +r 0 n+2 f * j -l 0 s,n+1 -J 0s φ * j -l 0 n+1 f * j +r 0 s,n+2 ,(5.2)
This is the Vlasov equation that will be solved in simulations of anisotropic model hereafter. With the notation used f j s,n , the bottom index n is the shell number, while the upper index j stands for the index of the Fourier space angle normalized by 2π M and it should accordingly be an integer. r 0 , l 0 and s 0 are the angular coupling indices between the wave-vectors. Due to self-similarity, only three values of coupling angle exist in LDM, which are given as follows:

r 0 = ⌈ M 2π arccos g 2 -g 4 -1 2g 2 ⌋ , s 0 = ⌈ M 2π arccos 1 -g 2 -g 4 2g 2 ⌋ , l 0 = ⌈ M 2π arccos g 4 -g 2 -1 2g 2 ⌋ , (5.3)
where the bracket ⌈•⌋ stands for rounding to closest integer. Note that the coupling angles r 0 , l 0 and s 0 that are linked to α p and α q give the exact phase information, which is the crucial difference between this 2D LDM model and the 1D shell models such as GOY or Sabra [START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF][START_REF] Vov | Improved shell model of turbulence[END_REF]. Moreover since LDM is obtained from a systematic reduction, it has no free parameters. LDM can be seen as an extension of GOY or Sabra model to two dimensions, so one can also find GOY or Sabra from LDM using the isotropy assumption. The prefactor 1 2 k 2 n g -4 µ 0 ,etc. are the coupling coefficients corresponding to the surface of the triads.

The representation of the nonlinear terms in equation( 5.2) automatically satisfies the general property of the Poisson bracket:

n, j ˆ∞ 0 J 0s φ j n J 0s φ, f s n, j E d E = 0 , (5.4) n, j ˆ∞ 0 f j s,n J 0s φ, f s n, j E d E = 0 , (5.5)
which represent the conservation of electrostatic potential energy and entropy, respectively. It means that energy and entropy are transferred among the triad such that they are conserved nonlinearly.

The quasineutrality equation is left unchanged:

c j n φ j n = ˆ J 0i f j i ,n -J 0e f j e,n E d E , (5.6)
with the prefactor c j n defined as:

c j n = a R 0 ˆ 1 -J 2 0i e (-ξ) ξd ξ + τ 5 2 a R 0 ˆ 1 -J 2 0e e (-τξ) ξd ξ + π 2 (τ + 1) a R 0 (1 -f t ) f t 1 -ǫ φ , (5.7)
where ξ ≡ E T i and τ ≡ T i T e and ǫ φ takes values between 0 and 1 for zonal flows and is 0 for the other modes. f t is the trapped particle ratio, which should take the values between 0 (no particles are trapped) and 1 (all the particles are trapped). In this model the passing particles are assumed to be adiabatic, thus for example if f t = 0 for electrons, it means that the electrons are not trapped and adiabatic. So we can take the value f t to be 0 to control the system. In this case the Poisson equation as well as the coefficient c j n should be redefined. Here we will give the Poisson equations and c j n for kinetic ions ( + adiabatic electrons) and kinetic electrons (+adiabatic ions). In kinetic ions system, f t = 0 for electrons, thus electrons are assumed to be adiabatic and there is only kinetic Vlasov equation for ions, the Poisson equation is given as follows:

c j n φ j n = ˆ J 0i f j i ,n E d E ,
with c j n defined as:

c j n = a R 0 ˆ 1 -J 2 0i e (-ξ) ξd ξ + π 2 a R 0 (1 + τ -f t ) f t 1 -ǫ φ .
In the kinetic electrons system, f t = 0 for ions, ions are all passing and assumed to be adiabatic according to this model, in this case there is only kinetic equation for electrons, and the Poisson equation should be:

c j n φ j n = ˆ -J 0e f j e,n E d E , with c j n defined as: c j n =τ 5 2 a R 0 ˆ 1 -J 2 0e e (-τξ) ξd ξ + π 2 a R 0 1 + τ -τ f t f t 1 -ǫ φ .

Numerical scheme

In section (3.4), three different linear cases have been discussed, corresponding to the pure TEM (kinetic electrons and adiabatic ions), pure TIM (kinetic ions and adiabatic electrons), and mixed TEM/TIM (kinetic ions and kinetic electrons) cases respectively. The nonlinear simulation of these three cases will be discussed in the sections hereafter. The simulations are performed in Python with help of Scipy ODE library (dop853) [START_REF] Oliphant | Python for scientific computing[END_REF] coupled to Fortran to Python interface (f2py) [START_REF] Peterson | F2PY: a Tool for Connecting Fortran and Python Programs[END_REF] for improving the performance of the time integration. The calculation of d f d t that is the most complicated part of the simulation, is implemented in a Fortran subroutine using Openmp [START_REF] Dagum | Openmp: An industry-standard api for shared-memory programming[END_REF] for parallel processing. A Gauss-Legendre quadrature scheme is used for energy integrals on typical intervals E ℓ ∈ [0 ; 8] with about 80 energy grid points (see Appendix F), resulting in typical simulation times of about 0.8 CPU hours per normalized time units in the simulations of the adiabatic system (which is usually integrated up to 100 normalized time units). The simulation of the full kinetic system usually takes much more time in simulation, because one needs to resolve the kinetic equations for both ions and electrons. However due to stronger linear injection as shown in figure (3.9), the full kinetic system saturate significantly earlier than the adiabatic system as will be shown later, thus the integration up to 30 normalized time units is considered to be enough for the simulation of the fully kinetic system, which finally takes about the same total simulation time as the adiabatic system in simulation (usually run up to 100 time units ). Note that the energy interval is determined by the convergence of the linear growth rate, with respect to the value given by the eigenvalue solver (see Appendix F).

Due to the multiscale nature of the TIM/TEM interaction problem, we use N = 24 shells for the fully kinetic system, while for the system with one adiabatic species, we use only N = 20 shells, which is enough for the nonlinear dynamics. The angles θ k are linearly discretized by dividing the wave number plane θ ∈ (0, 2π) into M regular sections with M = 54. With such parameters, the coupling angles defined in equation ( 5.3) take the values: r 0 = 21, l 0 = 10, s 0 = 23.

Definition of zonal flow and dissipation

Energy injection is ensured self-consistently by the linear terms in the Vlasov equation 5.2, which means that dissipations have to be added to this model in order to reach a quasi-stationary turbulent state. We again choose a form that localizes the dissipations in "buffer" regions in kspace using a hyper-viscosity of the form D = ν s k 4 + ν L k 6 , while the small scale dissipation is given by ν s and the large scale dissipation is given by ν L .

Zonal flows [START_REF] Diamond | Zonal flows in plasma -a review[END_REF][START_REF] Gürcan | Zonal flows and pattern formation[END_REF] are axisymmetric large scale flow structures that correspond to k y = 0 with k y being the binormal wave vector in tokamak 3D geometry (i.e. k α = 0 in this notation). In our simulations, where the discretization is log polar, there is no single well defined k y = 0 mode. In polar coordinates, the mode with j = 0 corresponds roughly to the zonal modes, however since it is defined as k = k n , θ j = -π M , + π M , it actually corresponds to a region around k α = 0, which becomes larger at smaller scales as shown in figure (5.1) and most of this region can not be considered as zonal flow modes (since k y can be quite large for many of the modes in this region, especially for large |k|).

As a result we have also defined a threshold on the shell number: n t h , to limit the zonal flow modes. For the modes in the direction j = 0, if the shell number n < n t h , it is seen as zonal flow mode, otherwise it is treated as fluctuations. The primary difference between the zonal flow modes and the turbulence modes in this model is that the zonal modes get no response from the passing particles as the parameter ǫ φ becomes 1 (since the effect of passing particle is canceled by the response from zonal flow). In this work we use n t h = 5 and ǫ φ = 1 . A constant dissipation ν 0 is applied to the zonal flow, which represents a constant linear drag.

E f i = ´| f j i ,n | 2 /F i /k n E d E (right)
averaged over the angular index j , as a function of the time t and the wave number k n . This is the case of TIM with adiabatic electrons.

(k n ) = c j n |Φ j n | 2 /k n (blue) and ion entropy E f i = | f j i ,n | 2 /F i /k n (
blue) averaged over the turbulent phase and the angles for the case of TIM with adiabatic electrons. The red and blue lines represent the k -5 and k -7/3 power law scalings.

Kinetic ions

Nonlinear simulations concern three different cases that have been identified in the linear phase, that are kinetic ions and adiabatic electrons system, kinetic electrons and adiabatic ions system, and the kinetic ions and kinetic electrons system. In order to simplify the notation, from now on, we will use kinetic ions, kinetic electrons and fully kinetic system to denote these three systems respectively.

In gyrokinetics, adiabatic electrons is a common simplifying assumption and is widely used [START_REF] Parker | Gyrokinetic simulation of ion temperature gradient driven turbulence in 3d toroidal geometry[END_REF][START_REF] Hammett | Developments in the gyrofluid approach to tokamak turbulence simulations[END_REF][START_REF] Lin | Turbulent transport reduction by zonal flows: Massively parallel simulations[END_REF][START_REF] Grandgirard | Computing itg turbulence with a full-f semi-lagrangian code[END_REF] in many codes. We present first the simulation for this case. and the ion entropy (also called the free energy) ) averaged over the angular index j , as a function of time t and the wave number k n . It can be seen that before t < 20, the linear growth dominates the dynamics. The most unstable modes (i.e. k n ∼ 10 as shown in figure 3.8) starts to saturate first, driving up the other modes via nonlinear transfer, which in turn begin to saturate via nonlinear transfer to dissipation regions. Note that the strongest mode is unchanged during the nonlinear phase. It is observed that the spectrum in the nonlinear phase is anisotropic, especially at large scales, and the strongest mode is located in the direction k ψ = 0, which is, as expected, a streamer. Comparing this 2D nonlinear spectrum and the 2D linear growth rate map that is shown in figure (3.8), it can be seen that the strongest mode in the nonlinear phase (i.e. k n ∼ 10 in figure 5.3) of the TIM corresponds exactly to the most unstable mode in the linear phase (i.e. k n ∼ 10 in figure 3.8). It means that the linear injection is still very important in the nonlinear phase of the TIM with adiabatic electrons.

Temporal spectrum of E φ and E f i

E f i = ´| f j i ,n | 2 /F i /k n E d E (right figure

Spectrum of electrostatic potential E φ in wave number plane

k spectra of E φ and E f i

The k-spectra of the conserved quantities is of main interest in turbulence study. In this kinetic ion system, the conserved quantities are the electrostatic potential energy E φ and the ion free energy E f i , or the ion entropy. The electrostatic potential energy E φ is defined as

E φ = c j n |Φ j n | 2 /k n ,
which should be averaged over the angles j to plot the k-spectra. The ion entropy is defined as

E f i (k n ) = ´| f j i ,n | 2 /F i /k n E d E .
Since the electrons are assumed to be adiabatic, the electron spectrum is determined by the spectrum of the electrostatic potential φ.

Figure (5.4) shows the k-spectra of the conserved quantities of this system: the electrostatic potential energy E φ = c j n |Φ j n | 2 /k n (blue) and the ion free energy

E f i = | f j i ,n | 2 /F i /k n (red)
, which are averaged over the angles and time in the nonlinear phase. The large scale dissipation defined by ν L = 2.0 and the small scale dissipation defined by ν s = 5.0e -8 are applied respectively to large and small scales, which cause the fall of the k-spectra in these two regions. Note that strong dissipation was necessary in the kinetic ions system in order to obtain a quasi stationary nonlinear phase. It is well known that in a two dimensional flow, the energy is transferred to large scales while the enstrophy is transferred to small scales [START_REF] Kraichnan | Inertial-range transfer in two-and three-dimensional turbulence[END_REF]. However in the kinetic ions system, the self-consistent energy injection is also located near the large scale (i.e. k ∼ 10). Thus two different physical mechanisms: energy transfer and linear energy injection co-exist at large scale, which finally leads to an accumulation. In order to obtain a stationary state, the dissipation at large scale had to be chosen as indicated.

In a range of wave-numbers where the linear injection and the dissipation can be ignored (which is called inertial range in Komogorov's theory), a power law scaling can be expected. Even though the existence of such an inertial range is questionable for the problem at hand, in some mid range far from the obvious energy injection and dissipation regions, the electrostatic potential energy, seems to follow a power law scaling of the form k -5 , while entropy seems to follow a k -7/3 scaling. Note that these are rough observational estimates rather than theoretical predictions.

Kinetic electrons

In gyrokinetics the assumption of adiabatic ion is not very common, however in TEM/TIM turbulence, TEM is usually more important than TIM, as shown in figure (3.8). Thus, nonlinear simulations of kinetic electron system may be more relevant than those of kinetic ions. In order to simplify the complexity of the kinetic problem, we start with the simpler case of kinetic electrons with the (admittedly nonphysical) assumption that ions are adiabatic.

Temporal spectrum of E φ and E f e

Since the grid are similar to the previous kinetic ion system, the number of ODEs in this simulation equals that of the previous one, therefore the simulation takes almost the same time. However it should be mentioned that the system of kinetic electron is easier to be saturated in simulation than that of the kinetic ion system. The problem to achieve a stationary state in such a model is mainly due to the dynamics at large scales, where the linear injection is important, even in nonlinear phase. For a two dimensional flow, the energy is transferred to the large scales and the enstrophy is transferred to the small scales. Since the strongest modes of the kinetic electron system is farther (i.e. k ∼ 30) from the large scale (i.e. the scales of zonal flow) than that in kinetic ion system (i.e. k ∼ 10), the injection and the transfer of energy will not accumulate at the large scale, so the kinetic electron simulation is easier to be stationary with a dissipation at large scales: ν L = 0.1 while in kinetic ion system, the large scale dissipation is ν L = 2.0 to achieve the stationary state.

The spectrum of electrostatic potential energy E φ = c 5.2), it can be seen that the nonlinear phase of the kinetic electrons case starts at around t ∼ 5, which is much earlier than that of the kinetic ions case. As a result the kinetic electrons system should saturate earlier than the kinetic ions system.

The first mode that starts to saturate is again the most unstable mode of the linear phase (i.e. γ = γ max ), which is true in kinetic electrons as well as in kinetic ions cases. The striking difference is that a shift of the strongest mode (i.e. the mode with highest energy and entropy) to larger scales is observed in the kinetic electrons case. The temporal behavior of the electrostatic potential energy and the electron free energy in figure (3.8) shows that the strongest mode is located at about k n ∼ 30 in the beginning of nonlinear phase (i.e. t ∼ 5), but as the energy is transferred to large scales in the process of nonlinear interactions, the strongest mode moves to k n ∼ 3 when the system is nonlinearly saturated at t ∼ 30. In contrast in the kinetic ions case (section (5.2)), the strongest mode always equals the most unstable mode in the linear phase, which is never changed during all the nonlinear phase.

Note that between the linear phase and the nonlinear phase, there is a phase of saturation, where the system is passing to a quasi-stationary nonlinear state, due to nonlinear energy transfer between different scales. This saturation phase is roughly from t ∼ 5 to t ∼ 30 in the kinetic electrons case, while it is between t ∼ 13 and t ∼ 20 in the kinetic ions case: the saturation phase of the kinetic electrons system is much longer than that of the kinetic ions system. The reason for this difference is unclear, but may be linked to the existence or not of available large scales in the system.

Predator-prey dynamics between zonal flow and turbulence modes

Predator prey dynamics is very well known in fusion plasmas [START_REF] Berionni | Predator prey oscillations in a simple cascade model of drift wave turbulence[END_REF]. Analyzing the temporal spectrum of the zonal flow and the turbulence, the predator prey like dynamics is also observed in the LDM model.

Figure (5.6) presents the temporal spectra of the absolute value of the electrostatic potential Φ. In upper figure, the blue line shows the mean value of the zonal flow modes while the green line shows the mean value of the rest of the modes, which can be seen as turbulence. We can say that the system is quasi stationary, nonetheless with strong oscillations in the spectrum of zonal flows, which is somewhat quasi periodic. Comparing the saturation levels, we can see that the strength of the zonal flow (∼ 10 -1 ) is less than the strength of turbulence(∼ 1). This is also an important difference with the kinetic ions simulation, where the mean value of zonal flow and the turbulence were almost at the same level. One possible explanation may be through the properties of linear injection: in the kinetic ions system, the zonal flow modes (i.e n = 1, 2, 3, 4, 5, j = 0) are close to the most unstable region (i.e n ∼ 7). This means that the zonal flows can gain energy almost directly from the most unstable modes by the nonlinear transfer. In contrast, for TEM, zonal flow is far from the the most unstable region (i.e n ∼ 11), which means that the energy injection near zonal flow scales is very low, therefore there is no direct interaction between most unstable modes and zonal flows, resulting in a low level of zonal flow and therefore a saturation through shifting towards smaller k. Note that since the growth rate of zonal flow is zero, the energy can never be injected directly to zonal flow. It can only gain energy indirectly from turbulence by nonlinear interactions.

Figure (5.6) also shows that the oscillations of zonal flow have a direct influence on the turbulence spectrum even though zonal flow is one order smaller than turbulence. It should be Figure 5.6: Time evolution of the average amplitude of electrostatic potential Φ is shown in semi log plot as a function of time t for the TEM case with adiabatic ions. In the upper figure the green line is the mean value of turbulence, the blue line is the mean value of the zonal flow (i.e n, j = n = 1, 2, 3, 4, 5; j = 0 ). The subfigure shows a zoomed in version of the deviation from the mean divided by the root mean square average of this same deviation, which brings out the oscillations in the fluctuations. In bottom figure the green line is the spectrum of one zonal flow mode with n, j = (1, 0), the blue line is the spectrum of one turbulence mode with n, j = (2, 10). Here the subfigure is simply a zoom in of the original figure .   noticed that the maximum (minimum) of zonal flow corresponds to the minimum (maximum) of turbulence, even though the response of the turbulence is very weak. This phenomenon is very similar to the predator prey dynamics. A comparison of a single zonal flow mode together with a single turbulence mode may demonstrate this phenomenon in more detail. Bottom figure shows the temporal evolution of one zonal flow mode (blue) with n, j = (1, 0) and one turbulence mode (green) with n, j = (2, 10). Note that these two modes can interact directly in the model due to the value of the coupling angles given in equation (5.3). A clear correlation in time can be observed between them, showing that when the zonal flow mode display a strong increase, the turbulence has a response at the same time usually with a drop. For example at t ∼ 50, a strong oscillation is observed in the zonal flow mode. Detailed time evolution of these two modes are presented in the sub figure. A phase shift of roughly π 2 can be seen between these two modes, which means that the maximum (minimum) of the zonal flow mode corresponds to the minimum of the turbulence mode: this is in accordance with the predator prey evolution, which can be observed in our model. This kind of phenomenon exists also in other kinetic LDM simulations. . It shows that the system is isotropic in small scales, while due to the anisotropic linear injection in large scales and the interplay with zonal flows, the system remains anisotropic in large scales. It can nonetheless be seen as isotropic because statistically the anisotropy is not significant and in a particular direction (i.e. predator-prey evolution keeps changing the direction of anisotropy not preferring a particular direction over long time). This result is different from the TIM case where the system was significantly anisotropic in the nonlinear phase.

Spectrum of E φ in wave number plane

Note that the isotropy or not of the linear phase may also have something to do with that final result. For example in figure 3.8, it was shown that the unstable zone of the kinetic ions system has the form of an ellipse, while the unstable zone of the kinetic electrons system is circular. In our model the nonlinear interaction takes place between neighboring shells. Therefore, while in kinetic ions system, (i.e. ellipse) the unstable modes (γ > 0 ) can interact directly with the stable modes (γ < 0) in the same region, in kinetic electrons system, (i.e. circle) the unstable modes mainly interact with other unstable modes. This difference in the interaction mechanisms may cause the differences that are observed in the nonlinear results. In order to achieve a stationary state, large scale dissipation defined by ν L = 0.1 and the small scale dissipation defined by ν s = 5.0e -8 as well as a constant dissipation ν 0 = 2.0 on zonal flows were applied. Note that the dissipations used in kinetic electrons case are much smaller than those used in kinetic ions case, especially the large scale dissipation : ν L = 0.1 of kinetic electrons vs ν L = 2.0 of kinetic ions. Since the strongest linear injection is located at log 10 k n ∼ 1.5 , which is far from the large scale, thus the energy production and energy transfer do not lead to a large accumulation at large scales, which explains why different large scale dissipations are needed to stabilize the system in different cases. Despite this difference, the power law scalings of the kinetic electrons system are the same as those of the kinetic ions system, which means the k-spectra of the electrostatic potential energy is close to k -5 , while the entropy is close to k -7/3 . Note that the inertial ranges are different in these two system, due to the different scales for electrons and ions.

k spectra of electrostatic potential energy E φ and entropy E f e

Fully kinetic system

The simulation of the fully kinetic system evolving both ion and electron distribution functions, is naturally more complicated than the simulation with one adiabatic species. So we left the detailed study of this case to the last part of the chapter, armed with the results from the simpler cases (section (5.2) & (5.3)). However, since with the same temperature gradient, the linear instability of the full kinetic system, as shown in figure (3.9) is much stronger than the others, the nonlinear phase arrives earlier. Thus we think that running the simulation from t = 0 to t = 30 is enough to resolve the saturated state.

Temporal spectrum of E φ , E f i and E f e

Three conserved quantities exist in the fully kinetic system, they are the electrostatic potential energy

E φ (k n ) = c j n |Φ j n | 2 /k n , the entropy of the ions E f i = | f j i ,n | 2 /F i /k n ,
and the entropy of the electrons E f e = | f j e,n | 2 /F e /k n . These quantities, averaged over the angles, are presented in figure (5.9) as a function of time t and the wave number k. The saturation phase begins at t ∼ 2, which is much earlier than those of the kinetic electrons (t ∼ 5) and the kinetic ions (t ∼ 13) cases, and is very short: at t ∼ 3 or 4, the system passes to the nonlinear phase, which is qualitatively similar to the kinetic ions system. In contrast the saturation takes place, through the shift of the strongest mode to larger scales, a phenomenon that has been observed in kinetic electrons simulations in section IV-B. As expected, the TEM dominates the linear phase in this fully kinetic simulation.

Spectrum of E φ in wave number plane

E φ (k n ) = c j n |Φ j n | 2 /k n (red), ion en- tropy E f i = | f j i ,n | 2 /F i /k n (blue) and electron entropy E f e = | f j e,n | 2 /F i /k n
averaged over the turbulent phase and the angles in the fully kinetic case. The red, blue and black lines present the k -4 , k -1.2 and k -1.8 power law scalings, respectively. 20 for the fully kinetic case. This nonlinear system is isotropic in small scale, as has been observed in the previous simulations. The system at large scale is anisotropic, and the form is like an ellipse, which is very similar to what has been observed in the kinetic ions simulation that is shown in figure (5.3). This means that the effect of kinetic ions at large scales is important in the nonlinear phase, even though it is the TEM that dominate the system in the linear phase.

k spectra of E φ , E f i and E f e

As the key result of the turbulence cascade model, the power law scaling of the quadratic invariants, that represent the nonlinear dynamics of the system, is of great interest. As shown in figure (5.11) in quasi-stationary nonlinear state, we find k -4 , k -1.2 and k -1.8 power law scalings for the electrostatic potential energy

E φ (k n ) = c j n |Φ j n | 2 /k n (red), ion entropy E f i = | f j i ,n | 2 /F i /k n (blue) and electron entropy E f e = | f j e,n | 2 /F i /k n (black)
respectively, which are quite different from the kinetic ions and kinetic electrons simulations. This significant discrepancy means that there is an important difference in the nonlinear dynamics of the fully kinetic system and the system with one adiabatic species. Note that in the simulation of the isotropic system [START_REF] Xu | A turbulent cascade model of bounce averaged gyrokinetics[END_REF], where ions and electrons are both kinetic as the simulation here, the power law scaling of E φ was also Figure 5.12: Electrostatic potential energy density E φ of LDM and GOY, suggests that the average spectrum is well captured by the isotropic model, and that the k -4 scaling that is observed at the end of the LDM range probably continues further in to smaller scales. found to be k -4 . The electrostatic potential energy density E φ of the anisotropic (2D LDM) and the isotropic (1D shell ) models are compared in figure (5.12). The spectra of LDM and GOY are almost parallel, which means the power law scalings of these two models are the same. The saturation level of LDM is higher than that of GOY, because the free parameter α in GOY model can change the saturation level. Since LDM gives the exact nonlinear interaction coefficients, thus the saturation level of LDM is more persuasive.

We have also compared the electrostatic potential energy spectra of the kinetic electrons and kinetic ions in the LDM and GOY, as shown in figure (5.13) In kinetic electrons simulations, the spectra are parallel between LDM and GOY, while in kinetic ions simulations, the spectra are not exactly parallel, due to the strong anisotropy of kinetic ions system as shown in figure (5.3). Hence the isotropy assumption is not very suitable for kinetic ions system.

Comparison of the k-spectra: adiabaticity

At this point, we have studied the k-spectra of three different systems: system with kinetic ions and kinetic electrons, system with kinetic ions and adiabatic electrons, and system with kinetic electrons and adiabatic ions. These spectra have been presented separately in different sections. In order to compare their difference as well as to understand the role of adiabaticity in the trapped particle turbulence, we plot the previous spectra in figure (5.14).

Figure (5.14) shows all the k-spectra of the conserved quantities, such as the electrostatic potential energy E φ (left figure), the ion free energy E f i (right figure) and the electron free energy E f e (right figure) of the kinetic system (presented by solid lines) and the systems with one adiabatic species (presented by dashed lines) in simulation. For E φ , we find k -5 in the simulations with adiabatic species, which is true for both kinetic ions and kinetic electrons. while in the full kinetic system, we find k -4 scaling for E φ . The entropy are very different for different systems: in fully kinetic system the power law scalings of ion entropy and electron entropy are respectively k -1.2 and k -1.8 , while in the system with one adiabatic species, the power law scalings of the entropy are both k -7/3 . It should be seen that the electrostatic potential energy is higher than the entropy at large scales and is lower than the entropy at small scales, which is valid for all these systems.

To get more information about the adiabaticity, we have run the simulations for the isotropic system that has been detailed in chapter (IV) and compare the difference between the kinetic and adiabatic system, the result is given in figure (5.15).

It is observed that the saturation levels are different for different regime. Comparing the blue one and the green one, the saturation level of kinetic electron (blue) is always higher than that of kinetic ion (green). Because the growth rate γ for kinetic electron is higher than the growth rate of kinetic ion. If we consider qualitatively the balance between the linear and nonlinear terms in stable state, which is γφ + αφ 2 ≈ 0, the system with stronger growth rate γ ( or smaller α ) will result in a higher saturation level, which has already been demonstrated by figure (4.5). And these two cases we find the same power law scaling k -5 , which is also the power law scaling of the 2D anisotropic simulation. However the power law scaling is different in very small scale, it shows a k -4 scaling for kinetic electron and a k -3 scaling for kinetic ion.

The interesting point is the spectrum of the case where ion and electron are both kinetic, which is presented by the red line with markers in figure (5.15). It is also the result presented in [START_REF] Xu | A turbulent cascade model of bounce averaged gyrokinetics[END_REF]. We find that the saturation level of the fully kinetic system (red line with markers) is between the cases with one species being adiabatic(blue and green). The detail is that in very Figure 5.14: The k-spectra as well as their power laws scalings of three different systems: the left figure presents the k spectra of the electrostatic potential energy E φ and the right figure presents the entropy spectra: E f i and E f e . Dashed blue lines show the results of kinetic electron system. Dashed red lines show the results of kinetic ion system. The solid lines present the results of fully kinetic system. large scale the saturation level is close to that of kinetic ion(green), but in very small scale it is close to the result of kinetic electron(blue). In other regions it is between the two. From this feature one may imagine that in nonlinear phase the saturation level of the fully kinetic system is not determined by TEM or TIM, which is dominated linearly, the nonlinear phase may be decided by ion in very large scale and by electron in very small scale due to the different scale of ion and electron. For the intermediate region it is determined by both ion and electron, which means in this region the interaction of ion behavior and electron behavior may be important and the system is passed from ionic physics to electronic physics. One may change the ion to electron Larmor radius ratio and the ion to electron banana width ratio to see this phenomenon, which can lead to a very interesting further study.

Conclusion

A logarithmically discretized model (LDM) of bounced averaged gyrokinetics has been derived and nonlinear, multiscale physics of TIM and TEM have been studied using this model. The model is found to be suitable for the study of nonlinear physics over a very large range of wavenumbers using a simplified kinetic formulation, similar in spirit to shell models, but with the ability to resolve anisotropy and thus a special treatment of zonal flows etc. This demonstration for bounce averaged gyrokinetics allows us to assess that similar formulation would be possible for full 5D gyrokinetics for the study of multi-scale physics of ions and electrons in Tokamak geometry.

The linear study of the model in a log-polar representation consistent with the LDM shows Figure 5.15: Blue line with square markers is the electrostatic potential energy spectrum for the system with kinetic electron and adiabatic ion. For this line the power law are k -5 (red line in the top) for large scale and k -4 (black line) for small scale. Green line with square markers is the spectrum for the case that ion is kinetic and electron is adiabatic, where the power law are k -5 (red line in the bottom) for large scale and k -3 (pink line) for small scale. The red line with markers is the result of previous work where ions and electrons are both kinetic.

that energy injection in TIM and TEM are both anisotropic but in slightly different ways such that while the peak of TIM with adiabatic electrons is rather narrow and localized around θ k = π/2 showing an elliptic form, the peak of the TEM with adiabatic ions is more diffuse and have almost a circular form. Fully kinetic calculations (where both ions and electrons are treated kinetically), show elevated levels of linear instability beyond a simple superposition of these two pure modes (i.e. where the other species were treated adiabatically). We argue that this difference propagates into the nonlinear regime, with TIM saturating via coupling to nearby scale damped modes (which exist because the instability boundary has the form of an ellipse), while the TEM couples to other unstable modes and therefore transfer energy to larger and smaller scales and finally saturate via dissipation resulting in a clear shift of the maximum of the wavenumber spectrum towards smaller wave-numbers. The fully kinetic simulations show characteristics of both modes: saturation via energy transfer to dissipation scales similar to the pure TEM case, with persistence of overall anisotropy similar to the pure TIM case. It was also shown that while in the kinetic ions case (section (5.2)), the zonal flow amplitude was about the same level with that of turbulence, in the kinetic electrons case (section (5.3)), it was much smaller. Strong oscillations of zonal flows, as well as their influence on turbulence spectra can be observed clearly in the kinetic electrons case, which is shown in figure (5.6). The time correlation between a zonal flow mode and a turbulence mode that is coupled to it, is clearly presented in figure (5.6), which is similar to the predator prey dynamics.

In the two cases where one species was treated kinetically, it has been observed that the wave-number spectrum is anisotropic only at large scales, as the energy injection is clearly anisotropic, and is isotropized rapidly, leading generally towards an isotropic k -5 spectrum for spectral potential energy density, which is different from the k -4 spectrum observed earlier using shell models [START_REF] Xu | A turbulent cascade model of bounce averaged gyrokinetics[END_REF]. A detailed study shows this difference is mainly caused by the adiabaticity (section (5.4)). Indeed, the wave-number spectrum is k -4 in the case where both species are kinetic, while it shows a k -5 scaling for adiabatic ions or adiabatic electrons.

From chapter 4, we understood that the nonlinear saturation level can be seen as a balance between the linear and nonlinear effect. Since the wave number of LDM is defined as k n = k 0 g n , it leaves 2 free parameters: k 0 and g , which may change the prefactors of the nonlinear terms and result in different nonlinear saturation level.

Note that the linear aspects of the model are validated against a separate linear eigenvalue solver, while nonlinear dynamics was observed to be qualitatively similar to previous results from the TERESA code [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] Gravier | Stimulated zonal flow generation in the case of tem and tim microturbulence[END_REF][START_REF] Gravier | Transport hysteresis and zonal flow stimulation in magnetized plasmas[END_REF], even though a direct comparison was not deemed particularly meaningful due to different assumptions that go into the two models. Detailed nonlinear comparison with TERESA should be considered for further studies.

In addition to the derivation of the model that is studied in this chapter, a general expression of the nonlinear convolution sum based on triad interaction condition (i.e. k + p + q = 0) is also presented in log polar coordinates in chapter 2, where different limits of this expression lead to different models such as the LDM (2D, local coupling), GOY and Sabra (1D, local coupling). The simulations that are discussed in this chapter concern the LDM, where the interactions are assumed to be local {n -2, n -1, n, n + 1, n + 2} but no assumption of isotropy was made.

The model presented in this chapter can be ameliorated in various ways: an extension to more complete kinetic descriptions such as 5D gyrokinetics or a better representation of zonal flows are two key directions. However understanding the details of the reduction procedure that is inherent in the LDM formulation is another interesting direction, which may allow the development of superior reduced models. For instance, a similar model that is being developed based on spiral decimation (where only one wave-vector direction is sampled at each scale) which has a performance between those of shell models and LDMs, provides a promising new direction [START_REF] Gürcan | Spiral chains in wavenumber space of two dimensional turbulence[END_REF]. One may be able to study the effect of reduction in detail by starting from a complete pseudo-spectral formulation and zeroing out Fourier modes in a systematic way until one arrives at the LDM or the spiral model at least over the small wave-number range where the pseudo-spectral model can operate. Such studies of the effect of reduction is necessary in order to understand how these models correspond to the original system, and whether or not a renormalization of the interaction coefficients is necessary.

Chapter 6 Conclusion and perspective

Turbulent transport is one of the keys to improve the energy confinement time required for thermonuclear fusion reactors. The description of the the kinetic turbulence of the plasma is a problem with 3 spatial and 3 velocity coordinates. Both theory and simulation of a problem with such high dimensionality are very difficult, and therefore reduced models are helpful to understand turbulence in Tokamaks. The common technique consists of averaging over the cyclotron motion, which is much faster than the turbulence time scale, which is called gyrokinetics. Such a reduction makes it possible to simplify the problem to three spatial coordinates of the particle guiding centers, a parallel velocity or energy and a perpendicular velocity appearing as the adiabatic invariant. Nonlinear gyrokinetic description requires massively parallel high performance numerical simulations. The difficulty lies in the nonlinear terms that describe the multi-scale interactions.

In a Tokamak the magnetic field is strong on the inner side and weak on the outer side. Thus, charged particles can be trapped or passing in such a magnetic configuration depending on their kinetic energy, therefore in addition to the cyclotron motion, a trapped particle also has bounce motion and toroidal precession motion. Since the characteristic time scale of the trapped electron mode (TEM) and the trapped ion mode (TIM) turbulence is in the order of the toroidal precession frequency, it allows us to average out both the cyclotron motion and the bounce motion, which finally results in a 4D kinetic system to describe the TEM and TIM turbulence in Tokamak. Even though this 4D description is greatly reduced compared to the 5D gyrokinetics, nonlinear simulation remains a challenge.

We have developed a new description of the TEM/TIM turbulence, which is mainly inspired by the nonlinear descriptions previously developed in the context of fluid turbulence [START_REF] Ohkitani | Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence[END_REF][START_REF] Vov | Improved shell model of turbulence[END_REF][START_REF] Gürcan | Logarithmic discretization and systematic derivation of shell models in two-dimensional turbulence[END_REF]. This new formulation of the Poisson bracket in Fourier space is obtained with a log polar coordinates system k = (k n , θ k ), where the wave numbers k n are logarithmically discretized and the angles θ k are regularly discretized. The logarithmical discretization with respect to the wave numbers presents a clear link to the models known as "shell models". The conservation of the electrostatic potential energy and the entropy, which are the essential nonlinear features of the Vlasov-Poisson system, is also preserved by our nonlinear description.

We have studied the linear system initially, which concerns the instability thresholds as a function of the external physics parameters and the linear growth rate of TEM and TIM. The linear growth rate is obtained by solving the linear dispersion relation, for which an eigenvalue solver has been developed based on the argument principle for this special problem where multiple roots and singularities co-exist.

It was observed that the linear dominated mode can be changed by varying the ion to electron temperature ratio τ. As the source of instability, a stronger temperature gradient κ T results in a more unstable system. Scanning the trapped particle ratio f t , we observe that the system is more unstable with more trapped particles. Note that according to this model, the trapped particles are kinetic while the passing particles are assumed to be adiabatic. The system is stable at small scale due to Landau damping.

It was also observed that a strong anisotropy of the linear growth rate exist in the wave number plane, which is true for both TEM and TIM: zonal flows are always stable and the most unstable modes are always the streamers. Under the same temperature gradient, the fully kinetic system is much more unstable than the system with one adiabatic species. Note that the form of the unstable zone in the wave number plane can be changed by Larmor radius ρ s and the banana width δ s , which was probably due to the definition of the gyro bounce average operator. Since the linear instability is very important even in nonlinear phase, the Larmor radius ρ s and the banana width δ s also have an important influence on the nonlinear system. The effect of these two parameters on the linear and the nonlinear system need more detailed understanding in future. For example we can scan the values of these two parameters to learn the instability threshold and see the difference on the linear as well as the nonlinear results.

The nonlinear studies consist of the simulation of the isotropic system and the anisotropic system. The isotropic system that consume less computer resource has been studied first. Similar to the 1D shell models, a free parameter α is left in the interaction coefficients and there is no exact phase information, which allows different phase approximations, which finally results in two different versions of the model, called Sabra and GOY. The simulations with different α gives the same power scaling in the range where the linear injection and the dissipation is negligible (i.e. the inertial range in Kolmogorov's theory), which means the free parameter α does not change the nonlinear dynamics. It was observed that the power law scaling is k -4 for the isotropic kinetic system (note that ions and electrons are both kinetic), which is different from the result that has been found in fluid and MHD turbulence, we guess that this is due to the average operator (Bessel function) and the kinetic integration, the analytic work on this subject is still in progress.

The toroidal precession resonance that is linked to kinetic plasma is clearly observed by our models. Note that the simplification in our models concerns only the spatial coordinates, the kinetic dependence rest unchanged, thus the kinetic effect should be kept by our description. This also means that the method of shell model that is initially developed for fluid turbulence can be applied to the kinetic plasma turbulence.

We have observed that the nonlinear temporal behavior of GOY and Sabra are totally different: GOY model shows an oscillating dynamics while such a phenomenon is not observed in Sabra model. Yet, the spectrum of the quadratic conserved quantities of GOY and Sabra follow exactly the same power law scaling, except that the saturation level of GOY model is higher than that of Sabra model under the same external parameters. These phenomena suggest that the correct description of the phase information may be very important.

Based on the understanding of GOY and Sabra, the 2D (+1D energy coordinate) anisotropic model in the wave number plane has been developed as well. Three different cases have been simulated by this model, that are: kinetic electrons (+ adiabatic ions ); kinetic ions(+adiabatic electrons) and fully kinetic system (ions and electrons are both kinetic).

It was observed that among these simulations, the system of kinetic electrons is the easiest one to achieve the quasi-stationary state with the minimum large scale dissipation. This is due to the different scales of electrons and ions: since the strongest linear injection of TIM (k ∼ 10) is close to the large scale and the energy is transferred to the large scale as has been observed in a 2D flow, the energy will accumulate in the large scale, which makes the kinetic ions system difficult to become stationary.

It was also observed that in kinetic ion system, the strongest mode in nonlinear phase equals the one that is the most unstable in linear phase, which means that the self-consistent linear injection is still important even in nonlinear phase. However in the other two cases, a transition of the strongest modes to large scales has been observed during the saturation phase since the energy is transferred to large scales.

The crucial difference between these simulations is that the fully kinetic system follows a k -4 power law scaling for the electrostatic potential energy, which is the same as that has been obtained in the isotropic system. While in the systems with one adiabatic species, the results are k -5 . This difference means that the fully kinetic system and the adiabatic system have very different nonlinear dynamics. Thus adiabaticity could be another interesting subject for future studies.

It was observed that the mean values of the zonal flow modes in kinetic ions system is much higher than that in kinetic electrons system. This was probably due to the nonlinear transfer mechanism of our model: in kinetic ions system, the strongest linear injection (k ∼ 10) is close to the large scale, thus zonal flow can gain energy easily from these modes, however in kinetic electrons, the strongest linear injection ( k ∼ 30) is far from the large scales, thus zonal flows modes can not gain energy by coupling directly to the linear injection, which results in a lower zonal flow.

We have also observed the predator-prey like dynamic between the zonal flows and the turbulence in our models: this is general and should be possible to study in our model, since in quasi stationary state, the energy, which can be seen as the sum of zonal flow modes and turbulence modes, is conserved. If the energy of zonal flows increases (decrease), then the energy of turbulence should decrease (increase), this is the mechanism of predator-prey dynamics. Furthermore regular oscillations of zonal flow modes have been observed in kinetic electron simulations.

Comparing the k-spectra of the electrostatic potential energy in isotropic system and anisotropic system, we found that the k spectra are parallel in the kinetic electron system and in the fully kinetic system. However due to the strong anisotropy of the kinetic ions system, the spectrum are different between the 2D and the 1D kinetic ion system. This means that the isotropy assumption is not very suitable in the study of the kinetic ion system.

We have run also some simulations with the 1D isotropic models to understand the adiabaticity. It has been found that the saturation level of the full kinetic system is between the cases with one species being adiabatic. The detail is that in very large scale the saturation level is close to that of kinetic ion system, but in very small scale it is close to the result of kinetic electron system. In other (intermediate) regions it is between the two. From this feature one may imagine that in nonlinear phase the saturation level of the full kinetic system is not determined by TEM or TIM, which is dominant linearly. The nonlinear phase may be controlled by ions in very large scale and by electrons in very small scale due to the different scale of ions and electrons. For the intermediate region it is determined by both ions and electrons, which means in this region the interaction of ion behavior and electron behavior may be important and the system is passed from ionic physics to electronic physics. One may change the ion to electron Larmor radius ratio and the ion to electron banana width ratio to see this phenomenon, which can lead to a very interesting further study.

This work gives the first application of shell models and LDM to kinetic plasmas, where some important results have been found based on these models. The idea here may be potentially used to discuss the nonlinearity of the other kinetic plasma systems, as well as fluid ones. The comparison between this reduction and the TERESA Code (where the simulations are performed in real space) is left to future studies. In the isotropic models one may insert the disparate scale interactions, where the nonlinear dynamics are probably different in this case. The adiabaticity has been demonstrated to be important by this work, which we hope will guide other kinetic or gyrokinetic simulation efforts.

Appendix A

Calculating the coefficients of GOY model

The coefficients of nonlinear terms can be obtained directly from equation (2.43) by the isotropy assumption. However a simpler approach consists of a direct calculation based on the conservation of electrostatic energy and entropy, which are given as follows:

n ˆ∞ 0 J 0s φ * n J 0s φ, f s n E d E = 0 , (A.1) n ˆ∞ 0 f s * n J 0s φ, f s n E d E = 0 . (A.2)
Note that these relations are obtained directly from the general property of Poisson bracket, that is the integration of Poisson bracket multiplied by one of its elements in the phase space should be zero.

Vlasov equation of GOY models is thus assumed to have the following general form:

∂ f s nl ∂t -i 2k n 2 ∂F s ψ ∂ψ J 0s φ n + i E l Ω d Z 2k n 2 f s nl +a n J 0s φ * n-2 f s * n-1,l + a ′ n J 0s φ * n-1 f s * n-2,l +b n J 0s φ * n-1 f s * n+1,l + b ′ n J 0s φ * n+1 f s * n-1,l +c n J 0s φ * n+1 f s * n+2,l + c ′ n J 0s φ * n+2 f s * n+1,l . with some coefficients a n ,b n , c n , a ′ n ,b ′ n , c ′ n .
In order to simplify the notation, the electrostatic potential φ n corrected by the Bessel function J 0s will be treated as a new quantity defined as Φ n : Φ n = J 0s φ * n . Then the conservation of electrostatic energy and entropy by the nonlinear terms of this model can be written as:

n (a n Φ n Φ n-2 f s * n-1,l + b n Φ n Φ n-1 f s * n+1,l c n Φ n Φ n+1 f s * n+2,l + a ′ n Φ n Φ n-1 f s * n-2,l b ′ n Φ * n Φ * n+1 f s * n-1,l + c ′ n Φ n Φ n+2 f s * n+1,l ) = 0 , n,l (a n Φ n-2 f s * nl f s * n-1,l + b n Φ n-1 f s * nl f s * n+1,l c n Φ n+1 f s * nl f s * n+2,l + a ′ n Φ n-1 f s * nl f s * n-2,l b ′ n Φ n+1 f s * nl f s * n-1,l + c ′ n Φ n+2 f s * nl f s * n+1,l ) = 0
, where the first constraint corresponds to the conservation of electrostatic potential energy and the second one corresponds to the conservation of entropy. Take the first constraint as an example, we can write the detail of the sum in this relation. One should care about the coefficients a n±2 , a n±1 etc for the shells in the boundary. For example when n = N , the nonlinear coefficients a N +1,N +2 , etc should be zero, because the modes f N +1 and f N +2 do not exist in the system. The detailed relation is given as follows:

c 0 Φ 0 Φ 1 f * 2 + b 1 Φ 1 Φ 0 f * 2 + c 1 Φ 1 Φ 2 f * 3 + a 2 Φ 2 Φ 0 f * 1 + b 2 Φ 2 Φ 1 f * 3 + c 2 Φ 2 Φ 3 f * 4 + ... ... ... a n-2 Φ n-2 Φ n-4 f * n-3 + b n-2 Φ n-2 Φ n-3 f * n-1 + c n-2 Φ n-2 Φ n-1 f * n + a n-1 Φ n-1 Φ n-3 f * n-2 + b n-1 Φ n-1 Φ n-2 f * n + a n Φ n Φ n-2 f * n-1 c ′ 0 Φ 0 Φ 2 f * 1 + b ′ 1 Φ 1 Φ 2 f * 0 + c ′ 1 Φ 1 Φ 3 f * 2 + a ′ 2 Φ 2 Φ 1 f * 0 + b ′ 2 Φ 2 Φ 3 f * 1 + c ′ 2 Φ 2 Φ 4 f * 3 + ... ... ... a ′ n-2 Φ n-2 Φ n-3 f * n-4 + b ′ n-2 Φ n-2 Φ n-1 f * n-3 + c ′ n-2 Φ n-2 Φ n f * n-1 + a ′ n-1 Φ n-1 Φ n-2 f * n-3 + b ′ n-1 Φ n-1 Φ n f * n-2 + c ′ n-1 Φ n-1 Φ n+1 f * n + a ′ n Φ n Φ n-1 f * n-2 = 0 .
here the first block contains the terms associated to a n ,b n and c n , while the second block is the terms associated to a

n ˆ∞ 0 J 0s φ * n J 0s φ, f s n = 0≤n≤N -2 a n+2 + c ′ n Φ n Φ n+2 f n+1 + 0≤n≤N -2 (c n + b n+1 ) Φ n Φ n+1 f n+2 + 2≤n≤N a ′ n + b ′ n-1 Φ n Φ n-1 f n-2 .
For the conservation of entropy, the calculation is very similar, so the detail is not shown here. After some manipulations, it can be represented as follows:

n ˆ∞ 0 f s * n J 0s φ, f s n = 0≤n≤N -2 a ′ n+2 + c n f * n f * n+2 Φ n+1 + 0≤n≤N -2 c ′ n + b ′ n+1 f * n f * n+1 Φ n+2 + 2≤n≤N (a n + b n-1 ) f * n f * n-1 Φ n-2 .
In order to preserve the conservations in equation (A.1) and equation (A.2), the following relations of the coefficients should be satisfied:

a n = -c ′ n-2 ,c ′ n = -b ′ n+1 , c n = -b n+1 a ′ n = -b ′ n-1 ,a ′ n = -c n-2 ,a n = -b n-1
which means:

a n = -b n-1 = -c ′ n-2 = b ′ n-1 = -a ′ n = c n-2 . (A.3)
Take into account the wave number k n = k 0 g n in the coefficients, it means:

a n = αk n-2 k n-1 = -a ′ n , b n = βk n-1 k n+1 = -b ′ n , c n = γk n+1 k n+2 = -c ′ n . (A.4)
From equation(A.3) and equation(A.4), it is found:

αk 2 0 g 2n-3 = -βk 2 0 g 2n-2 = -γk 2 0 g 2n-1 =⇒ α = -βg = -γg 2 .
It means that all the coefficients can be represented by the wave number k n and a free parameter α. a n ,b n ,c n , etc. are givens as follows:

a n = -a ′ n = αk 2 n g -3 , b n = -b ′ n = -αk 2 n g -1 , c n = -c ′ n = -αk 2 n g .
Then the Vlasov equation of this model is:

∂ f s n ∂t =i k n J n 0s φ n ∂F s ψ ∂ψ -i E Ω d Z s k n f s n + αk 2 n g -3 J n-2 0s φ * n-2 f s * n-1 -J n-1 0s φ * n-1 f s * n-2 -αk 2 n g -1 J n-1 0s φ * n-1 f s * n+1 -J n+1 0s φ * n+1 f s * n-1 + αk 2 n g J n+1 0s φ * n+1 f s * n+2 -J n+2 0s φ * n+2 f s * n+1 . (A.5)
Since the nonlinearity is truncated to the same couplings, the coefficients of Sabra model are exactly the same as that of GOY model.

Appendix B The linear dispersion relation solver B.1 Linear dispersion relation ǫ (ω)

In this appendix, we want to calculate the linear dispersion relation ǫ (ω) from the linear Vlasov-Poisson system. The linear Vlasov equation is obtained by setting the nonlinear terms to zero:

∂ f s k ∂t = i k α J 0s φ k ∂F s ψ ∂ψ -i E Ω d Z s k α f s k ,
Assuming plane wave solutions for the fluctuations, i.e. f s E , k α , k ψ , t ∝ f s (E ) e -i ωt +i k α α+i k ψ ψ with k = k ψ , k α , the linear Vlasov equation can be represented as follows:

i ω f s k = -i k α J 0s φ k ∂F s ψ ∂ψ + E Ω d Z s i k α f s k , (B.1) c k φ k = ˆ J 0i f i k -J 0e f e k E d E , (B.2)
where f s k can be expressed as a function of by the equilibrium state F s and of the electrostatic potential φ as follows:

f s k = k α J 0s φ k ∂F s (ψ) ∂ψ -ω + E Ω d Z s k α .
Inserting this relation of f s k in the Poisson equation, it can be found:

c k φ k = ˆ   J 0i   k α J 0i φ k ∂F i (ψ) ∂ψ -ω + E l Ω d Z i k α   -J 0e   k α J 0e φ k ∂F e (ψ) ∂ψ -ω + E l Ω d Z e k α      E d E .
Since the electrostatic potential φ k doesn't depend on the energy E , it can be removed from the energy integration, which gives:

c k = ˆ   J 0i   k α J 0i ∂F i (ψ) ∂ψ -ω + E l Ω d Z i k α   -J 0e   k α J 0e ∂F e (ψ) ∂ψ -ω + E l Ω d Z e k α      E d E ,
where the prefactor c k :

c k = a R 0 ˆ 1 -J 2 0i exp (-ξ) ξd ξ (B.3) + τ 5 2 a R 0 ˆ 1 -J 2 0e exp (-τξ) ξd ξ + π 2 (τ + 1) a R 0 (1 -f t )
f t is introduced and the normalized energy ξ ≡ E /T i is defined for each species. By eliminating the fluctuating part of the distribution function associated with each species f s k , as well as assuming a Maxwellian form for the background distribution function

F s ψ = n s (ψ) T s (ψ)
e [-E/Ts(ψ)] , the plasma dielectric function ǫ (ω) reads as follows:

ǫ (ω) = c k - ˆ κ n + κ T ξ -3 2 Ω d (ξ -ω) J 2 0i exp (-ξ) ξd ξ -τ 3 2 ˆ κ n + κ T τξ -3 2 Ω d (ξ + ω) J 2 0e exp (-τξ) ξd ξ, (B.4)
where κ n = d ψ ln n and κ T = d ψ ln T stand respectively for logarithmic gradients of background plasma density and temperature and ω is the mode frequency normalized by the toroidal wave number k α .

The integrals in equation (B.4) admit singularities (i.e.ξ = ±ω) when the frequency is purely real (ℑ(ω)=0). Since the Fourier-Laplace transformation as given in (3.4) is only valid for ℑ(ω) > 0, residue contributions for both integrals have to be taken into account for the analytic continuation of the ǫ(ω) into the lower half of the complex plane ℑ (ω) ≤ 0. Since the energy coordinate ξ is always positive, the relevant poles correspond to ξ = ω for ℜ (ω) > 0 and to ξ = -ω for ℜ (ω) < 0. The connection with TIM and TEM becomes apparent here with the origin of these two different modes coming from the signs of electron and ion charge, inducing opposite mode frequencies [START_REF] Kadomtsev | Trapped particles in toroidal magnetic systems[END_REF].

For a given integral

I (ω/k) = ´+∞ -∞ d u g (u)
u-w/k the integral in the full plane is defined (for k > 0) as:

I (ω/k) ≡        ´+∞ -∞ d u g (u) u-ω/k ℑ (ω/k) > 0 ´+∞ -∞ d u g (u) u-ω/k + j πg (ω/k) ℑ (ω/k) = 0 ´+∞ -∞ d u g (u) u-ω/k + j 2πg (ω/k) ℑ (ω/k) < 0 (B.5)
Thus, the plasma dielectric function, which is valid everywhere in the complex plane, can be written in the form:

ǫ (ω) → ǫ (ω) + λ j i πg (ω) , (B.6)
where the ǫ(ω) on the right corresponds to the definition given in (3.4), valid only for ℑ(ω) > 0, and λ j = 0, 1, 2 if respectively ℑ(ω) > 0, ℑ(ω) = 0 or ℑ(ω) < 0. The function g (ω) is residue contribution, which are g i (ω) and g e (ω) for ions and electrons respectively: At this point, the linear dispersion relation as well as its residues in the complex plane is clear. The frequency ω exist in a complicated form in the linear dispersion relation, in other words, it can not be found analytically. So we have to apply numerical methods in order to solve this eigenvalue problem. Since two kinds of instability, TIM and TEM coexist in the linear system, it is necessary to find all the possible roots. Here we use a method based on the argument principle, which is introduced in complex analysis [START_REF] Johnson | Enclosing all zeros of an analytic function -a rigorous approach[END_REF]. This principle presents for a given function f (z), the number of roots N and the number of poles P in a fixed rectangular domain R can be obtained by the integral

g i (ω) = - κ n + κ T ω -3 2 Ω d J 2 0i exp (-ω) ω , ∀ ℜ (ω) > 0 , (B.7) g e (ω) = τ κ n + κ T -τω -3 2 Ω d J 2 0e exp (τω) -ω , ∀ ℜ (ω) < 0 . (B.8)

B.2 Eigenvalue solver

´R ∂ z f (z) f (z) d z, which means: 1 i 2π ˆR f ′ (z) f (z) d z = N -P , (B.9)
For an analytic function like the linear dispersion relation introduced here, there is no pole, and so P = 0, which means that this integral gives exactly the number of roots in the rectangle. If there is no root in R, the calculation will stop for this region. If there is only one root in R some root finding solver based on Newton-Raphson method will be used to find the exact root more efficiently. If more than one root exist in R, the rectangle domain R should be divided into subregions and then the number of roots in each subregion should be checked. The division process should be repeated until only one root exists in the associated subregion. Then the root finding solver can be used to find the exact location of the root. In practice, it is always necessary to choose the size of the final subregion as small as possible such that the root finding solver can finish the last step faster. For example in this work the size of the final rectangle is 1, so the root exist in a small squared region.

To use this method, the derivative of plasma dispersion function with respect to the argument ω should be given explicitly ( f ′ (z) in equation (B.9)). The derivation of the linear dispersion relation as well as its residues are given as follows:

∂ǫ (ω) ∂ω = - κ n + κ T ξ -3 2 Ω d (ξ -ω) 2 J 2 0i exp (-ξ) ξ + τ 3 2 κ n + κ T τξ -3 2 Ω d (ξ + ω) 2 J 2 0e exp(-τξ) ξ , ∂g i (ω) ∂ω = - κ T Ω d J 2 0i (ω) exp (-ω) ω -κ n + κ T ω - 3 2 2 Ω d J 0i (ω) ∂J 0i (ω) ∂ω exp (-ω) ω -κ n + κ T ω - 3 2 1 Ω d J 2 0i (ω) (-1) exp (-ω) ω -κ n + κ T ω - 3 2 1 Ω d J 2 0i (ω) exp (-ω) 1 2 ω , ∂g e (ω) ∂ω = -τ 5 2 κ T Ω d J 2 0e (ω) exp (τω) -ω + τ 3 2 κ n + κ T -τω - 3 2 2 Ω d J 0e (-ω) ∂J 0e (-ω) ∂ω (-1) exp (τω) -ω + τ 3 2 κ n + κ T -τω - 3 2
1

Ω d J 2 0e (-ω) exp (τω) -ω + τ 3 2 κ n + κ T -τω - 3 2
1

Ω d J 2 0e (-ω) exp (τω) -1 2 -ω ,
where the derivation of the gyro bounce averaged operator J 0s (ω) with respect to the frequency is given as:

∂J 0s (ω) ∂ω = -J 1 τωδ s k ψ δ s k ψ τ 2 ω J 0 τωρ s k α -J 0 τωδ s k ψ J 1 τωρ s k α ρ s k α τ 2 ω , for J 0i (ω), τ = 1, while for electron τ = τ = T i
T e . J 0 stands for the Bessel function of first kind of integral order 0, J 1 stands for the Bessel function of first kind of integral order 1. The differentiation rule of Bessel function J ′ 0 (x) = -J 1 (x) is used here. The numerical method is implemented in Python. Quadrature schemes are used for energy integration. The process of integration along the rectangle and the division of a rectangle into subregions are realized in a module, which is called in main function. In figure B.1 this numerical method based on the argument principle is tested with a randomly selected wave number k (i.e take the case of isotropy k ψ = k α = 20 for simplicity).

The red and black lines correspond respectively to the zeros of the real part and imaginary part of the linear dispersion relation in the complex plane of the frequency ω = ω r + j ω i . So the points where the two lines cross should be the roots (or eigenvalues) of the linear dispersion relation. The points marked by red squares and stars are the roots found by the numerical method introduced here. It shows an agreement between the numerical method and the roots presented by the linear dispersion relation which means the numerical method finds roots correctly and it can be used to solve the eigenvalue problem.

In this figure it can also be seen that the linear dispersion relation has more than one root in the complex plane of the frequency. In physics one is interested mainly by the mode with the maximum imaginary frequency w i , it is the maximum value of ω i that dominate the linear system when the time t is long, as has been detailed in section (3.1). which is the Bessel function of first kind of integral order 0. The average over the bounce motion follows the same logic as the gyro average, which means the banana orbit can be approximately seen as the cyclotron motion where the Larmor radius ρ is replaced by the banana width δ [START_REF] Drouot | A gyro-kinetic model for trapped electron and ion modes[END_REF]. Then we will discuss the implications on the quasi neutrality equation. For the first order, the density does not depend on the gyro phase:

n(r) = ˆf (r, µ, v ∥ )d 2 v .
For the second order, one need to consider the adiabatic response, which is due to the fact that the adiabatic invariant µ is not an exact invariant. The adiabatic invariant holds for the case where the temporal variation of the field is slower than the cyclotron motion and the spatial variation is large compared to the cyclotron radius (∂ t l og B ≪ ω c and ∇B B ≪ ρ c ). So one need to consider the variation of µ in the scale of the cyclotron radius:

f (R, µ, v ∥ , φ) = f (R, µ, v ∥ ) + ∂F ∂µ d µ .
The second term represent the second order contribution or the adiabatic perturbation with F the equilibrium distribution. For a population with constant energy E = µB + qφ = c st , d µ can be obtained as follows:

d µ = - q B d φ = - q B φ (r) -φ (R) ,
where φ r and φ R represent respectively the electrostatic potential of the real trajectory and the guiding center. Take F as Maxwellian: F = n eq T 3/2 eq e (-E T eq

) we can obtain: ∂F ∂µ = -B T F thus:

f (R, V) = f (R, µ, v ∥ ) + q T F φ (r) -φ (R) .
The density should be:

n(r) = ˆJ0 f (R)d 2 v + ˆq T F φ (r) -φ (R) d 3 v .
Applying the gyro average operator for the potential, we can obtain

n(r) = ˆJ0 f (R)d 2 v + ˆq T F φ (r) -J 0 φ (r) d 2 v = ˆJ0 f (R)d 2 v + ˆq T F φ (r) -J 2 0 φ (R) d 2 v ,
here we use the fact that the potential of the real coordinate can be seen as an average over all the possible guiding centers. Since the electrostatic potential depend only on the spatial coordinate, it can be simplified as:

n(r) = ˆJ0 f (R)d 2 v + ˆq T F 1 -J 2 0 φ (r) d 2 v .
So the second term, which is called gyro correction, can be seen as the polarization of the plasma caused by the error of the potential between the guiding center and the real coordinate.

Appendix D Comparison of J 0s and its approximations

The gyro bounce average operator defined as:

J 0s = J 0 (k ψ δ s0 E T s ) J 0 (k α ρ s0 E T s ) , (D.1)
where J 0 stands for the Bessel function of first kind, which depends on the banana width δ s0 of the bounce motion, the Larmor radius ρ s0 , the particle energy E and the temperature T s . J 0 (k ψ δ s0 E T s ) and J 0 (k α ρ s0 E T s ) present the corrections of the bounce motion and the cyclotron motion, respectively. Considering the lowest order Padé approximation of Bessel function, J 0s can be simplified as:

J ap 0s = 1 1 + ξ k 2 α ρ 2 s0 4 1 1 + ξ k 2 ψ δ 2 s0 4 , (D.2)
which is widely used in gyrokinetics when the scale is not very large (i.e max(kρ)~10 ). Here ξ ≡ E T i

, to simplify we use T i = T e . However during this thesis the wave number is discretized logarithmically, so a large range of wave number was accessible. The difference of the exact Bessel function and its approximations, especially the correction of the coefficients c k in Poisson equation as well as its corrections on the linear instability caused by its approximation will be discussed here.

In figure D.1 the exact gyro bounce average operator J 0s and its approximation J ap 0s with respect to the wave number k is compared. It can be seen that the gyro bounce average operator is close to to 1 when the wave number k is small , while it is close to 0 when the wave number k is large. The exact operator has some oscillation around 0, so it can be positive and negative. The oscillations become more and more weak when the wave number increase. However in the approximation, it is a smooth line without any oscillation and it is always positive even though it is close to zero.

In Poisson equation the coefficient c k is the sum of the polarization (also called gyro correction) c pol and the adiabatic response c ad . The gyro correction c pol is a function of the gyro bounce operator, which is given as:

c pol = ˆ 1 -J 2 0i exp (-ξ) ξd ξ + ˆ 1 -J 2 0e exp (-ξ) ξd ξ , (D.3)
which is usually simplified as Laplace operator △ in gyrokinetics [START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF][START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF], that is:

c pol = ˆ 1 -J 2 0i exp (-ξ) ξd ξ + ˆ 1 -J 2 0e exp (-ξ) ξd ξ = △ i + △ e . (D .4)
This can be seen as a Taylor approximation, there exist another approximation, which is: 3). Red line presents the approximation where the gyro bounce averaged operator J 0s is replaced by its approximate form in equation (D.2), which can be seen as the Padé approximation. The black line is the result that is approximated directly to △ in gyrokinetics, this can be seen as Taylor approximation.

c pol = ˆ 1 -J 2 0i exp (-ξ) ξd ξ + ˆ 1 -J 2 0e exp (-ξ) ξd ξ = ˆ     1 - 1 1 + ξ k 2 α ρ 2 i 0 4 2 1 1 + ξ k 2 ψ δ 2 i 0 4 2      exp (-ξ) ξd ξ + ˆ     1 - 1 1 + ξ k 2 α ρ 2 e0 4 2 1 1 + ξ k 2 ψ δ 2
The blue line in figure (D.2) presents the exact gyro bounce correction c pol presented in equation (D.3). The red line presents the approximation in equation (D.5) where the Bessel function J 0s is replaced by its approximation given by equation(D.2), which can be seen as Padé approximation. The black line is the result that is approximated to a Laplace operator △ in gyrokinetics, where the integration over the kinetic coordinates disappears. The one is the simplest and the most widely used form in gyrokinetics. It is obtained by using Taylor approximation. Comparing these results, it can be seen clearly that these two approximations works well at large scale (i.e. k < 10). However at small scale only the Padé approximation rest somewhat valid. It can be seen that the difference between the Laplace operator △ and the exact gyro bounce correction c pol is huge at small scales, which means this widely used technique can not be used in the system where small scale is taken into account.

In shell models, the range of wave number k is very large. In order to obtain the correct result, it is better to use the exact Bessel function without any approximation. Or at least use the approximation given in equation (D.5). Since the growth rate γ obtained from the exact gyro bounce correction and the Laplace operator △ are very different.

From figure(D.3), it can be seen that the linear instability of TEM as well as TIM obtained with exact gyro bounce correction is higher than that obtained with the Laplace operator △ The result with exact c pol is marked by squares, while the result obtained with the Laplace operator approximation is marked by points, which is also the result given in [START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF].

. γ max is ∼ 7 compared to γ max ∼ 4 in case of approximation. The system becomes stable for k > 45 in the exact case, while with the Laplace operator, the system becomes stable while k > 30. The difference between the Laplace operator approximation △ i +△ e and the exact gyro bounce correction c pol is very huge, and affect the reliability of a code especially at large values of k. However, some numerical schemes, resolving the unknowns in direct space, requires Laplace operator or Padé approximations in order to compute the gyro corrections as derivatives space. In our study, the use of Fourier representation allows the use of the exact gyro bounce average operators, in terms of Bessel functions, which is useful since we require a valid description up to very small scales.

Appendix E Electrostatic potential energy E φ and entropy E f s

Definition of the conserved quantities

Let us remind the general definition of a Poisson bracket associated to coordinates q i ; p i , in a phase space of dimension 2N : f ; g p,q ≡ N i =1 ∂ p i f ∂ q i g -∂ p i g ∂ q i f , from which it comes after integrations by parts: ˆ f ; g p,q d pd q = 0 , ˆf f ; g p,q d pd q = ˆ f 2 2 ; g p,q

d pd q = 0 , ˆg f ; g p,q d pd q = ˆ f ; g 2 2 p,q d pd q = 0 .

In our problem, the nonlinear term in the gyro-bounce averaged Vlasov equation can be written as a Poisson bracket with respect to the two coordinates ψ, α :

∂ t f s ψ, α, E , t = - ∂ J 0s φ ∂α ∂F s ∂ψ - E Ω d Z s ∂ f s ∂α + J 0s φ ; f s ψ, α -D f s ,
where J 0s φ = J 0s φ ψ, α, E , t , and f s = f s ψ, α, E , t , and a dissipation operator -D f s can be required for numerical satbility reasons.

The last term, written in an explicit Poisson bracket form, should accordingly respect the conservation properties: In our formulation, the unknowns are expressed in Fourier space k ψ , k α , and we remind the Vlasov and quasi-neutrality equations :

∂ t f s, k = i k α ∂F s ∂ψ J 0s φ k + i k α E Ω d Z s f s, k -D f s, k + k+p+q=0 p ψ q α -p α q ψ J 0s φ ⋆ p f ⋆ s, q , c k φ k = s Z s ˆJ0s f s, k E d E .
The budget associated to the first nonlinearly conserved quantity is obtained by multiplying the Vlasov equation by

f ⋆ s k
F s and integrating along k ψ , k α and energy:

∂ t E f s = L f s -D f s , (E.3)
where the cancellation of the nonlinear term, inherited from the property of the Poisson bracket (E.1), has been used, so that only remains the contribution of the linear terms corresponding to energy injection:

L f s = i k k α φ k ˆ∂lnF s ∂ψ J 0s f ⋆ s ,k E d E ,
where the second linear term, associated to the advection by the precession drift Ω d cancels due to the symmetry with respect to ±k α :

k f ⋆ s, k i k α E Ω d Z s f s, k = i E Ω d Z s k k α f 2 s, k = 0 .
The dissipation operator reads as follows:

D f s = k d k 2 ˆ f s, k 2 F s E d E ,
where the dissipation operator is symbolized by the function d k 2 , usually taking the form: d k 2 = ν L k -6 + ν s k 4 , in order to dissipated large scales with the rate ν L and small scales with the rate ν s . The first conserved quantity then reads:

E f s ≡ k ˆ f s, k 2 2F s E d E ,
where it is important to notice that the energy integration is not mandatory with respect to the cancellation of the nonlinear term, so that the previous budget could also be expressed locally in E .

The second conservation property (E.2) suggests to multiply the Vlasov equation by Z s J 0s φ ⋆ k before integrating along Fourier modes and energy, in order to use the quasi-neutrality equation and to finally obtain:

∂ t E φ = L φ -D φ , (E.4)
with the definition of the injection due to the precession frequency advection:

L φ ≡ i Ω d s k k α φ k ˆJ0s f ⋆ s, k E E d E ,
while the dissipation operator reads:

D φ = k d k 2 φ k 2 ,
and the second conserved quantity reads, by use of the quasi-neutrality equation:

E φ ≡ k c k φ k 2 .

Link with the fluxes

In the two previous budgets (E.3, E.4), the right hand sides contains dissipation terms-D f s and -D φ that are by construction negatively definite, and appear due to the need of dissipation for numerical stability. On the other hand, the contributions from the linear terms can be written as follows: It is interesting to notice that the terms injecting energy, onto the natural quadratic conserved quantities of the system, are directly linked to the volume averaged radial particle and heat fluxes, that are of prime importance for radial transport in Tokamak plasmas.

L f s = k v E×B.

Link with the free energy

A general definition of the free energy states that this quantity represents the total available energy from which a thermodynamical system can produce some work, in other words, it is given by the difference between the total, internal energy U , and the non-available energy, given by the produce of entropy S and temperature T :

E ≡ U -T δS = ˆH F d v + T ˆF ln F d v ,
where we separate each unknown F = F 0 + δ f , and H = H 0 + δH , into an equilibrium and some fluctuations such that δ f /F 0 ∼ δH /H 0 << 1, so that we obtain for each species s:

ˆFs ln F s d v = ˆFs0 ln F s0 d v + ˆ(1 + ln F s0 ) δ f s d v + ˆδ f 2 s 2F s0 d v , ˆHs F s d v = ˆHs0 F s0 d v + ˆHs0 δ f s d v + ˆδH s F s0 d v + ˆδH s δ f s d v .
With the help of the definition of the equilibrium distribution function F s0 = e -H s0 /T s , the free energy can then be written for each species s:

E s = ´Fs0 δH s d v+T s ´δ f s d v+ ´δH s δ f s d v+T s ´δ f 2 s 2F s0 d v ,
where the zeroth order is automatically cancelled due to the equilibrium. In this work, we are interested into the characterization of the electrostatic gyro bounce averaged fluctuations δ f s and δH s = q s J 0s φ, so that the first order is of no interest and we will study mainly the spatial averages of the two quantities that are quadratic with respect to the fluctuations:

E f s ≡ k ˆTs δ f s, k 2 2F s0 d v , E φ ≡ k c k φ k 2 ,
where a summation over the species s has allowed to make use of the quasi-neutrality equation, in order to write the second quantity as proportional to the square root of the electrostatic potential.

It appears now clearly why we refer to the first quantity as an entropy (the temperature being usually neglected since it corresponds in our studies to a constant prefactor), and to the second quantity as the electrostatic potential energy. In contrast, the Gauss Legendre quadrature is defined for the integral in the interval [-1 , 1], for any given function. Thus, before using this method one needs to change the variables such that the interval for the integration is transformed to [-1 , 1]. For example the integral of f (x) in any given interval [a , b], we can change variables as follows:

ˆb a f (x) d x = b -a 2 ˆ1 -1 f b -a 2 x + a + b 2 d x .
Thus take a = 0, the integral is defined in the interval [0, b], we have: If we define the energy grid in the n t h interval as :

E n ℓ = a n -a n-1 2 x ℓ + a n + a n-1 2 ,
we can get:

ˆan a n-1 f (x) d x = a n -a n-1 2 ℓ ω n ℓ f E n ℓ ,
where x n ℓ and ω n ℓ are abscissa and the weights as above. The energy grid E n ℓ can be calculated from x l , a n and a n-1 . The integral in the full interval [a, b] can be obtained as the sum over all subintervals:

ˆb a f (x) d x = N n=1 (a n -a n-1 ) 2 ℓ ω n ℓ f E n ℓ .
In pratice we divide the interval [0, ∞] into some subregions regularly as [0, 2], [START_REF] Azevedo | Selection of fuel cladding material for nuclear fission reactors[END_REF][START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF], [START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF][START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF], [START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF][START_REF] Catto | Linearized gyro-kinetics[END_REF]..., etc. If we define ∆E = a n -a n-1 , the energy grid in the each subinterval can be obtained regularly as:

E n ℓ = a n-1 + ∆E 2 x ℓ + ∆E 2 ,
then the integral in the full interval can be written as:

ˆb a f (x) d x = ∆E 2 N n=1 ℓ ω ℓ f ∆E 2 x ℓ + a n-1 + ∆E 2 .
Take all the subinterval as one interval, then the number of the ℓ t h energy grid in the n t h subinterval

E n ℓ is m = ℓ + (n -1) * N r ,
where N r is the number of grids in each subinterval, thus the sum can be simply written as:

ˆb a f (x) d x = ∆E 2 m ω m f (E m ) ,
where ω m = ω nℓ = ω ℓ . E m = E n ℓ = ∆E 2 x ℓ + ∆E 2 + n∆E . ω l and E l of Gauss Legendre quadrature are presented in figure (F.1).

In order to verify the validity of these approaches, and choose the best one, we will compare the growth rates obtained with these three numerical integral method against the value from linear dispersion relation. the result is shown as follows:

The results are shown in figure (F.2), where the left figure shows the results of the isotropic case and the right figure shows the results of anisotropic cases with the plane 0, π 2 divided into M = 9 regular sections. From the left figure we can see that trapezoidal and Gauss Laguerre quadrature works better than the Gauss Legendre quadrature, however Legendre quadrature is much more efficient than the trapezoidal method. The small difference between the value of the The integral method is written in the figures. For trapezoidal method we use 500 regular points for E max = 10. For Laguerre quadrature, we use 96 points for energy grid. For Legendre quadrature, we use 80 points for E max = 8 in the left figure (i.e. 20 points in each subinterval [0, 2], [START_REF] Azevedo | Selection of fuel cladding material for nuclear fission reactors[END_REF][START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF], [START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF][START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF], [START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF][START_REF] Catto | Linearized gyro-kinetics[END_REF]), while in right figure we use 100 points for E max = 8, with 40 points in the interval [0 , 2]. IVP method and that from the eigenvalue solver is mainly due to the Landau damping mechanism, which is difficult for numerical implementation. In the right figure we use the Gauss Legendre quadrature, the maximum energy is 8, which is determined by the convergence. In the interval [0, 2] where the resonance exists, we use 40 points. For the other intervals, we use 20 points. It can be seen that the convergence is better than that in left figure .  Mots clés : Turbulence, gyrocinétique, Tokamak, TIM, TEM, Physique Nonlinéaire Résumé : La turbulence et ses mécanismes d'autorégulation régissent le transport dans les dispositifs de confinement magnétique. La turbulence cinétique du plasma en présence de champs magnétiques puissants, telle qu'elle est présente dans ces dispositifs, est décrite à l'aide de la formulation gyrocinétique. En raison de l'inhomogénéité du champ magnétique dans le tokamak, il est possible de simplifier encore la description en moyennant sur le mouvement de rebond, pour obtenir un système 4D décrivant la turbulence des seules particules piégées TIM (trapped ion mode) et TEM (trapped electron mode) dans le plan toroïdal du Tokamak. Au cours de cette thèse, une hiérarchie de modèles réduits a été établie, basée sur des reformulations et réductions successives du produit de convolution régissant la dynamique non-linéaire entre modes de Fourier. On distingue en particulier des modèles anisotropes, basés sur l'approche LDM, et des modèles isotropes de type modèles en couches. Abstract : Turbulence and its self-regulating mechanisms govern transport in magnetic confinement devices. The kinetic turbulence of the plasma in the presence of strong magnetic fields, as present in these devices, is described using the gyrokinetic formulation. Due to the inhomogeneity of the magnetic field in the tokamak, it is possible to further simplify the description by averaging the rebound motion, to obtain a 4D system describing the turbulence of trapped particles: trapped ion mode (TIM) and trapped electron mode (TEM) in the toroidal plane of the Tokamak. In this thesis, a hierarchy of reduced models has been established, based on successive reformulations and reductions of the convolution product governing nonlinear dynamics between Fourier modes. In particular, anisotropic models, based on the LDM approach, and isotropic models inspired by shell models, have been developed.

The resulting models are cheaper implementations of the nonlinear multi-scale physics involving trapped electron mode (TEM), trapped ion mode (TIM) and large scale structures (i.e. zonal flow). The non-linear codes are firstly validated by comparing the linear growth rates with the eigenvalue of the linear dispersion equation. The parametric dependencies as well as the linear instability of TIM and TEM are studied in great detail. The entropy and electrostatic potential energy spectra obtained show relatively similar exponents whether the description chosen is isotropic or anisotropic, except for kinetic ions, where strong streamers exist in Fourier Plane. However, the saturation level and the temporal dynamics are highly dependent on whether the wave number anisotropy is taken into account or not. Finally, the adiabaticity hypothesis for the response of a species turns out to strongly modify the turbulent spectra.

v vi 6 F

 6 Conclusion and perspective A Calculating the coefficients of GOY model B The linear dispersion relation solver B.1 Linear dispersion relation ǫ (ω) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 B.2 Eigenvalue solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 C Formulation of gyro correction D Comparison of J 0s and its approximations E Electrostatic potential energy E φ and entropy E f s Method of numerical integration Chapter 1

  5Mev) + n (14.1Mev)D + D → T (1.01Mev) + p (3.03Mev) D + D → He 3 (0.82Mev) + n (2.45Mev) D + He 3 → α (3.67Mev) + p(14.67Mev)

Figure 1 . 4 :

 14 Figure 1.4: The motion of a charged particle in a Tokamak, due to the gradient of the magnetic field, a charged particle can be passing (right upper figure) or trapped (right lower figure) depending on its kinetic energy. The left upper figure shows the parallel motion of the passing particle and the left lower figure shows the bounce motion and the toroidal precession motion of the trapped particle[18].

Figure 2 . 2 :

 22 Figure 2.2: Expansion of the nonlinear terms (Poisson bracket) in the wave number plane by the permutation (right vs left) and the reflection with respect to k(up vs bottom)

Figure 3 . 1 :

 31 Figure 3.1: Threshold for TIM (red) and TEM (blue) as a function of κ T and the wave number k.

Figure 3 . 2 :

 32 Figure 3.2: Threshold for TIM (red) and TEM (blue) as a function of κ T and the ion to electron temperature ratio τ ≡ T i /T e .

Figure 3 . 3 :

 33 Figure 3.3: Threshold for TIM (top) and TEM (bottom) as a function of κ T and the trapped particle ratio f t , for three different values of the temperature ratio τ.

Figure 3 . 4 :

 34 Figure 3.4: Validation of root finding algorithm: the blue line and the green line present the points where the imaginary part ℑ [ǫ (ω)] and the real part ℜ [ǫ (ω)] of the linear dispersion relation equal to zero, the start and the squares are the roots found by the algorithm. The left one is tested for k = 20 and the right one is tested for k = 100.

2 :

 2 t ,e f t ,i τ a/R 0 Ω d κ T 0.01 0.01 2Standard parameters used for TIM (top) and TEM(bottom) studies.

Figure 3 . 8 :

 38 Figure 3.8: TIM (left) and TEM (right) linear growth rates in the wave vector plane k = k ψ , k α = (k n cos(θ k ), k n sin (θ k )) for the kinetic ions (adiabatic electrons) and kinetic electrons (adiabatic ions), respectively.

Figure 4 . 1 :

 41 Figure 4.1: Comparison between the growth rates are shown as functions of the wave vectors k n obtained by solving the dispersion relation (red and blue crosses for respectively TIM and TEM) and from the initial value problem (IVP, black squares). Three values for the temperature ratio, τ = 1.0, τ = 0.6, and τ = 0.2 (with f t = 2/3 and κ T = 0.25 in all cases) were considered. These three cases, which are presented from top to bottom for decreasing τ, correspond to TEM dominated, mixed TEM/TIM, and TIM dominated regimes respectively.

Figure 4 . 2 :

 42 Figure 4.2: E φ (k n ) = c n |Φ n | 2 /k n as a function of k n , with a time average over the turbulent quasistationary phase for the parameters given in table 4.1 with two different values of large scale dissipation [i.e. ν L = 0.1 (blue) and ν L = 1.0 (green)]. Note that ν L determines how rapidly the electrostatic potential energy density falls at large scales, but does not change the rest of the spectrum. However due to this effect, the maximum of the fluctuation level is in fact determined by ν L .

Figure 4 . 3 :

 43 Figure 4.3: Ion entropyE f i (k n ) = ´Cn | f i | 2 /F i /k n E d E and and electron entropy E f i (k n ) = ´Cn | f i | 2 /F i /k n E d Eas a function of k n , with a time average over the turbulent quasi-

Figure 4 . 4 :

 44 Figure 4.4:E f i =

Figure 4 . 6 :

 46 Figure 4.6: Comparison of the electrostatic potential |Φ n | between GOY and Sabra models for n = 6, which is a value of k n a little higher than the value associated to the maximum in the kspectrum. Interestingly, while the GOY model shows distinct, large amplitude oscillations, the Sabra model, seems to evolve more randomly.

Figure 4 . 7 :

 47 Figure 4.7: Comparison of different modes (n = 6, 12, 18) in GOY model.

Figure 4 . 8 :

 48 Figure 4.8: comparison of entropy in GOY version and Sabra version, the blue lines are the entropy of electrons, while the red lines are the entropy of ions. The results of GOY model (top) are marked by points and the results of Sabra model (bottom) are marked by squares.

Figure 4 . 10 :

 410 Figure 4.10: Nonlinear oscillation frequency obtained by computing the fast Fourier transform in the turbulent phase as a function of the large scale dissipation ν L , showing that the nonlinear frequency increases as the large scale dissipation is increased, while the other parameters are unchanged.

Figure 5 . 1 :

 51 Figure 5.1: The log polar coordinate of LDM: the wave-vector modulus k is logarithmically discretized as k = k n = k 0 g n and the angles are discretized onto a regular grid θ j ≡ j 2π M with g = 1.424, k 0 = g -2 and M = 54 used in simulation.

Figure 5 . 2 :

 52 Figure 5.2: Spectrum of the log of electrostatic potential energy E φ = c j n |Φ j n | 2 /k n (left) and ion

Figure 5 . 3 :

 53 Figure 5.3: The log of electrostatic potential energy E φ = c j n |Φ j n | 2 /k n in log-log wave vector space [i.e. log k n sin θ j , log k n cos θ j ], averaged over the turbulent phase t > 50 for the case of TIM with adiabatic electrons.

Figure 5 . 4 :

 54 Figure 5.4: k-spectra of the electrostatic potential energyE φ (k n ) = c

Figure 5 .

 5 Figure 5.2 shows the spectra of the electrostatic potential energy E φ ≡ c j n |Φ j n | 2 /k n (left figure)

Figure 5 .

 5 Figure 5.3 shows the electrostatic potential energy E φ = c j n |Φ j n | 2 /k n in the full wave number plane, averaged over the turbulent phase.It is observed that the spectrum in the nonlinear phase is anisotropic, especially at large scales, and the strongest mode is located in the direction k ψ = 0, which is, as expected, a streamer. Comparing this 2D nonlinear spectrum and the 2D linear growth rate map that is shown in figure(3.8), it can be seen that the strongest mode in the nonlinear phase (i.e. k n ∼ 10 in figure5.3) of the TIM corresponds exactly to the most unstable mode in the linear phase (i.e. k n ∼ 10 in figure 3.8). It means that the linear injection is still very important in the nonlinear phase of the TIM with adiabatic electrons.

  k n (left figure) and the electron free energy or entropy E f e = ´| f j e,n | 2 /F e /k n E d E (right figure) averaged over the angular index j , are presented in figure (5.5), as a function of time t and the wave number k n . Compared to the kinetic ions case that has been discussed in figure (

Figure 5 . 5 :

 55 Figure 5.5: Spectrum of electrostatic potential energy E φ = c j n |Φ j n | 2 /k n (left) and electron free

Figure 5 . 7 :

 57 Figure 5.7: Electrostatic potential energy E φ (k) = c j n |Φ j n | 2 /k n in the wave number plane k = log 10 (k n ) cos (θ k ) , log 10 (k n ) sin (θ k ) , averaged over the turbulent phase t > 60.

Figure 5 .

 5 Figure 5.7 gives the spectrum of electrostatic potential energy E φ (k) = c j n |Φ j n | 2 /k n in the full wave number space averaged over the turbulent phase (i.e t > 50). It shows that the system is isotropic in small scales, while due to the anisotropic linear injection in large scales and the interplay with zonal flows, the system remains anisotropic in large scales. It can nonetheless be seen as isotropic because statistically the anisotropy is not significant and in a particular direction (i.e. predator-prey evolution keeps changing the direction of anisotropy not preferring a particular direction over long time). This result is different from the TIM case where the system was significantly anisotropic in the nonlinear phase.Note that the isotropy or not of the linear phase may also have something to do with that final result. For example in figure3.8, it was shown that the unstable zone of the kinetic ions system has the form of an ellipse, while the unstable zone of the kinetic electrons system is circular. In our model the nonlinear interaction takes place between neighboring shells. Therefore, while in kinetic ions system, (i.e. ellipse) the unstable modes (γ > 0 ) can interact directly with the stable

Figure 5 . 8 :

 58 Figure 5.8: k-spectra of electrostatic potential energy E φ (k n ) = c j n |Φ j n | 2 /k n (blue) and electron entropy E f e = | f j e,n | 2 /F i /k n (black) averaged over the turbulent phase and angles for the TEM case with adiabatic ions.. The blue and black lines represent the k -5 and k -7/3 power law scalings, respectively.

Figure 5 .

 5 Figure 5.8 shows the k-spectra of the electrostatic potential energy E φ = c j n |Φ j n | 2 /k n (blue) and the electron free energy E f e = | f j e,n | 2 /F e /k n(black). In order to achieve a stationary state, large scale dissipation defined by ν L = 0.1 and the small scale dissipation defined by ν s = 5.0e -8 as well as a constant dissipation ν 0 = 2.0 on zonal flows were applied. Note that the dissipations used in kinetic electrons case are much smaller than those used in kinetic ions case, especially the large scale dissipation : ν L = 0.1 of kinetic electrons vs ν L = 2.0 of kinetic ions. Since the strongest linear injection is located at log 10 k n ∼ 1.5 , which is far from the large scale, thus the energy production and energy transfer do not lead to a large accumulation at large scales, which explains why different large scale dissipations are needed to stabilize the system in different cases.

Figure 5 . 9 :

 59 Figure 5.9: Spectrum of the electrostatic potential energy E φ = c j n |Φ j n | 2 /k n (top), ion entropy E f i = | f j i ,n | 2 /F i /k n (middle) and electron entropy E f e = | f j e,n | 2 /F e /k n (bottom) are presented as a function of time t and the wave number k n for the fully kinetic case.

Figure 5 . 10 :

 510 Figure 5.10: Electrostatic potential energy E φ (k) = c j n |Φ j n | 2 /k n in the wave number plane k = log 10 (k n ) cos (θ k ) , log 10 (k n ) sin (θ k ) , averaged over the turbulent phase t > 20 for the fully kinetic case.

Figure ( 5 .Figure 5 .

 55 Figure (5.10) presents the electrostatic potential energy E φ (k) = c j n |Φ j n | 2 /k n in the log-log wave number plane k = log 10 (k n ) cos (θ k ) , log 10 (k n ) sin (θ k ) , averaged over the turbulent phase t >

Figure 5 . 13 :

 513 Figure 5.13: Electrostatic potential energy density E φ of LDM and GOY. The left figure present the result of kinetic ions simulation and the right figure present that of kinetic electrons simulation.

′

  n ,b ′ n and c ′ n . From the detailed relation we can find same terms with different coefficients, which have been marked by the same color in above relation. After arrangement, this equation can be rewritten in a simple form as follows:

Figure B. 1 :

 1 Figure B.1: Validation of root finding algorithm: the blue line and the green line present the points where the imaginary part ℑ [ǫ (ω)] and the real part ℜ [ǫ (ω)] of the linear dispersion relation equal to zero, the start and the squares are the roots found by the algorithm, tested for k = 20.

Figure D. 1 :

 1 Figure D.1: comparison of gyro bounce average operator and its approximation for ion and electron. The left figure is for J 0i and the right figure is for J 0e . The blue line is the Bessel function and the red line is the approximation in both figures.

  can be seen as the Padé approximation. The exact value of the coefficient c pol in equation (D.3) and its two approximations in equation (D.4) and equation (D.5), with respect to the wave number k have been compared in figure (D.2).

Figure D. 2 :

 2 Figure D.2:Comparison of the gyro bounce correction c pol and its approximations. Blue line presents exact value without any approximation as given by equation (D.3). Red line presents the approximation where the gyro bounce averaged operator J 0s is replaced by its approximate form in equation (D.2), which can be seen as the Padé approximation. The black line is the result that is approximated directly to △ in gyrokinetics, this can be seen as Taylor approximation.

Figure D. 3 :

 3 Figure D.3: Linear instability γ obtained with different c pol in Poisson equation: the blue and red lines presented respectively the linear growth rate of TEM and TIM.The result with exact c pol is marked by squares, while the result obtained with the Laplace operator approximation is marked by points, which is also the result given in[START_REF] Drouot | Global gyrokinetic simulations of trapped-electron mode and trapped-ion mode microturbulence[END_REF].

  ˆfs J 0s φ ; f s ψ, α d ψ d α =0 . (E.1) ˆJ0s φ J 0s φ ; f s ψ, α d ψ d α =0 , (E.2)

  T s L φ = Ω d s 〈Q ψ, T s 〉 , where we noticed that the Fourier transformed pseudo-radial component of the electric drift velocity reads v E×B.e ψ = i k α φ k , and we introduced the volume averaged particles and heat fluxes: 〈Γ ψ, n s 〉 = k v E×B.e ψ δn s , 〈Q ψ, T s 〉 = k v E×B.e ψ δn s δT s .

ˆ+∞ 0 ff

 0 (x) * e -x = ℓ ω ℓ f (x ℓ ) ,where x ℓ and ω ℓ are the abscissa, which are the root of Laguerre polynomials and the weights, respectively, ℓ is the number of grid . However for any given function f (x), in order to use Gauss Laguerre quadrature, we need to multiply the function by e x , since: (x) e x * e -x = ℓ ω ℓ f (x ℓ ) e x ℓ .

1 f 2 ˆ1 - 1 f

 121 ℓ f (E ℓ ) , withx ℓ and ω ℓ are the abscissa, which are the root of of Legendre polynomials and the weights, respectively and ℓ is the grid number. Note that the energy grid is not x l , it should be E ℓ = b 2 (x ℓ + 1). Now we can use Gauss Legendre quadrature with b = E max . For the integral in a large interval such as [0, +∞], the interval can be divided into N multiple regions, such as (take [a, b] for example): (x) d x , here we use a 0 = a and a N = b. a n-1 and a n are the first and last point of each subinterval, respectively. The integral in the n t h subinterval is given as follows:ˆan a n-1 f (x) d x =a na n-1 a na n-1 2x + a n + a n-1 2 d x .

Figure F. 1 :

 1 Figure F.1: The energy grid E l and the corresponding weight ω l of Gauss Legendre quadrature.

Figure F. 2 :

 2 Figure F.2: Verification of the initial value code (IVP) with different energy integral methods. Hereκ T = 0.25, f t = 2/3, τ = 1, TEM dominate the linear phase, thus the results should be compared to the TEM growth rate. Left figure presents the linear growth rate of the isotropic model (i.e.k α = k ψ ). The right figure present the linear growth rate of anisotropic model, the plane 0, π 2 is divided into M = 9 regular sections.The integral method is written in the figures. For trapezoidal method we use 500 regular points for E max = 10. For Laguerre quadrature, we use 96 points for energy grid. For Legendre quadrature, we use 80 points for E max = 8 in the left figure (i.e. 20 points in each subinterval [0, 2],[START_REF] Azevedo | Selection of fuel cladding material for nuclear fission reactors[END_REF][START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF],[START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF][START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF],[START_REF] Brown | Vode: A variable-coefficient ode solver[END_REF][START_REF] Catto | Linearized gyro-kinetics[END_REF]), while in right figure we use 100 points for E max = 8, with 40 points in the interval [0 , 2].

  Les modèles obtenus permettent une description numériquement moins coûteuse, et couvrant plusieurs décades en nombre d'onde, de la physique non-linéaire voyant interagir des modes TEM, TIM, et des structures à grande échelle. Les codes nonlinéaires sont dans un premier temps validés par comparaison des taux de croissance linéaire avec ceux donnés par la résolution aux valeurs propres de l'équation de dispersion. Les dépendances paramètriques ainsi que les instabilités linéaires des modes TIM et TEM sont étudiées en détail. Les spectres en énergie libre et en énergie potentielle électrostatique obtenus montrent des exposants relativement similaires, que la description retenue soit isotrope ou anisotrope. Cependant, le niveau de saturation et la dynamique temporelle s'avèrent dépendre fortement de la prise en compte ou non de l'anisotropie en nombre d'onde. Enfin, l'hypothèse d'adiabaticité pour la réponse d'une espèce s'avère modifier fortement les spectres turbulents. Title : Study of reduced models for kinetic plasma turbulence Keywords : Turbulence, gyrokinetic, Tokamak, TIM, TEM, Nonlinear Physics

Table 4 .

 4

	1.4 1/g 2 g 2 0.25 10 -12 0.1 1

1: Standard parameters for a reference nonlinear simulation of a shell model of trapped particle kinetic turbulence. Others parameters are the same as in table 3.1.

  e ψ κ n s δn s -κ T s 〈Γ ψ, n s 〉 + κ T s 〈Q ψ, T s 〉

			3 2	δn s -δn s	δT s T s	,
	= κ n s -	3 2	κ T s

Appendix C Formulation of gyro correction

The effect of gyro average is related to the simple fact that particles do not sustain exactly the same fields as their guiding centers. In the same way a field is affected not by the guiding center but by the distribution of particles that the guiding center represents. The gyro average operator is defined as:

where r and R represent respectively the real coordinate and guiding center coordinate, which satisfy the relation r = R + ρ with ρ the cyclotron radius and φ the gyro phase that is between 0 and 2π. The Fourier transform of the guiding center f (R, µ, v ∥ ) is :

) * e i ka , we can obtain:

where we can define the gyro average operator as follows:

Method of numerical integration

In Poisson equation we have the integral over the energy from 0 to ∞, which is:

In simulation we need to find a method to compute this integral numerically. We have tried three different methods: Gauss Legendre quadrature, Gauss Laguerre quadrature and trapezoidal method. The code based on different methods need to be verified by comparing the linear growth rate against that obtained from the linear dispersion relation. Even though the interval of the kinetic energy is [0 , +∞], the contribution from the high energy values is not significant, which means that we can use the convergence of the linear growth rate, with respect to the value of linear dispersion relation, to determine the interval of energy (i.e. find E max for which the value of the integral converges).

Trapezoidal method is the simplest one in numerical analysis. For an integral ´b a f (x) d x, the trapezoidal rule works by approximating the region under the graph of the function f (x) as a trapezoid and calculating its area s i . The integral can be seen as the sum of all the trapezoid areas s i , which means:

This is the basic expression of trapezoidal rule. This rule is written as a module in Python numpy library. One only needs to define a discrete range of energy and use this module to compute the integral.

Gauss Laguerre quadrature is defined in the interval [0, +∞] for the function with a special form like: ´+∞ 0 f (x) e -x , The integral result can be witten as: