Keywords: angle variables and motion invariants, 2, 2, 3 Passing and trapped particles: orbit width, 44, 58 Nonlinear interaction between energetic particle modes and turbulence 69

The multi-scale problem in complex systems 1.1 Complex systems: from the smallest entity to the global behaviour 1.2 Some examples of complex systems . . . . . . . . . . .

Abstract

On the route towards nuclear fusion as a clean, safe and potentially limitless worldwide source of energy, many challenges remain to be met within the framework of the ongoing ITER project. Nuclear reactions will produce a significant fraction of energetic particles, which must be sufficiently well confinement in order to transfer their energy to the thermal plasma, ensuring this way the steadystate production of energy. In the presence of energetic particles, plasma instabilities can be excited, degrading the confinement and the overall performance of the device. In addition, other instabilities are at play, eventually leading to turbulent transport. When energetic particles and plasma instabilities co-exist, their mutual interaction should not be neglected. In addition, the analysis of energetic particles must be done within a kinetic framework, describing the trajectories in a 5D phase-space. The self-consistent study of this problematic is the goal of the present dissertation, which focuses on the kinetic description of energetic particles in the presence of ion temperature gradient driven turbulence and macro-tearing modes. Special attention is paid to describing the analytical framework employed within a Hamiltonian approach based on action-angle variables, together with the state-of-the-art gyro-kinetic simulations.

iii

Foreword

The results presented in this dissertation are the outcome of my research enriched by strong international collaborations with my colleagues at Oak Ridge National Laboratory (USA), the University of Saskatchewan (Canada), the University of Seville (Spain), the National Institute for Fusion Science (Japan) and the Max-Planck Institute for Plasma Physics (Germany). Also, these collaborations have been reinforced by my active participation in EUROfusion projects on the verification of gyro-kinetic codes, on the nonlinear energetic particle dynamics and on the use of reduced models for the multi-scale transport of energetic particles. An essential part of the research presented here has been possible thanks to a project I have coordinated, funded by A*Midex, the Initiative of Excellence of the Aix-Marseille University. Within the framework of my project, I had the opportunity to supervise master students, PhD candidates and one post-doctoral researcher.

I would like to stress that in the past years I have interacted with other researchers outside the fusion community, mainly in the fields of industrial plasma processing and medical plasma applications (for the treatment of cancer). The discussions I had with these colleagues were extremely productive and led me to some digressions to think about my research in the context of a broader field, which is that of dynamical complex systems. Part of these reflections are given in chapter 1, for I consider it necessary to give a brief overview of other complex systems which can help us understand a bit more the one I analyse in detail: the fusion plasmas. At the end of that chapter, I present the main topic of the dissertation as a special type of dynamical complex system and I introduce the different instabilities I focus on in the remainder of the manuscript.

It is important to emphasize that much of the calculations that I present in this dissertion, especially those of chapter 2, are the outcome of many thoughts resulting from interactions and discussions motivated by questions coming from my master and PhD students as well as my post-doc. These questions have helped me acquire new perspectives on the way towards the understanding of the fundamental mechanisms in fusion plasmas. In that respect, I have tried to write down in a single document a theory which is sometimes missing in most of the manuals or that is scattered in many separate articles, books or proceedings. I hope this effort will be useful for the next generations of students, especially when looking for the mathematical background to understand the origin of plasma instabilities. Many questions are addressed, such as what a resonance is, what the drive of an instability means, how the general expression of the drive can be applied to different types of instabilities, how particle trajectories are described or what the v gyro-kinetic approach represents. Although one can always try to include more material in one dissertation, I must have found a compromise in order to include also new results from my project with my students and post-doc as well as from the interaction with my international collaborators.

Chapters 3 and 4 present the results of my research initiated to understand the real origin of the so-called EGAM, the existence of the EGAM channel between energetic and thermal particles and how the EGAM can directly impact the trajectories of particles, producing anomalous transport. Part of this work has been done in collaboration with Oak Ridge National Laboratory. Some considerations regarding the integration of particle trajectories were the result of the supervision of a master student and discussions with my collaborators at CEA (France).

Chapter 5 presents detailed analysis of the nonlinear coupling between non axi-symmetric modes and EGAMs, enabling a channel for the energy to flow from EGAMs to turbulence. Part of these results were obtained during the supervision of a post-doctoral researcher.

In chapter 6 I present results on the kinetic effects on the tearing instability, which has been historically analysed using fluid approaches. Part of this work is the result of the supervision of a PhD student and the collaboration with the University of Saskatchewan. At the end of that chapter I present results on the particle trajectories in the presence of tearing modes, which are obtained within the framework of an ongoing EUROfusion project and motivated by fruitful discussions with my collaborators at the University of Seville.

Last but not least, this dissertation does not contain all my research activities. Perspectives and new directions are explained in chapter 7. For instance, the analysis of the Alfvénic modes in turbulent regimes is the subject of a PhD I currently supervise. Also, the obtention of reduced models for the transport of energetic particles and the use of new numerical techniques such as those of the Artificial Intelligence constitute some of my current activities.

Evolution of the zeros of the dispersion relation with increasing energetic

particle concentration with q = 3. . . . . . . . . . . . . . . . . . . . . . . 3.3 Frequency (left) and growth rate (right) of EGAM as predicted in [1] (lines) and as measured from GYSELA simulations (symbols).

The scan is performed using the same species for the thermal population and modifying the mass and the charge for the energetic particle population. . 5.12 For the three simulations used in this chapter, plotted is the turbulent intensity at the end of the simulations as measured by 5.24 for a sufficiently long time window where the system exhibits a statistically steady state. For reference, the case without EP source is given (dashed thin line). The shadowed region represents the radial localisation of the EGAM. . The region inside the yellow line represents the area around the Opoint used to calculate the averaged quantity δI turb . For reference, the separatrix is given by a black line. . . . . . . . . . . . . . . . . . 6.7 Scan on island rotation frequency (ω isl ) and island width (W ). The vertical axis represents the modification of turbulence around the O-point of the island, normalized to the reference turbulence value of a simulation without island. A negative (resp. positive) value implies reduction (resp. increase) of turbulence. . . . . . . . . . . . 6.8 Normalized density gradient R/L n averaged around the O-point region as a function of the island rotation frequency and the island width. For the sake of clarity only one error bar is provided, measured as the standard deviation in the radial direction of the normalized density gradient. . . . . . . . . . . . . . . . . . . . . . . 6.9 Normalized ion (left) and electron (right) temperature gradients R/L T,i , R/L T,e averaged around the O-point region as a function of the island rotation frequency and the island width. For the sake of clarity only one error bar is provided, measured as the standard deviation in the radial direction of the normalized ion temperature gradient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xv

LIST OF FIGURES 6.10 Time evolution of the turbulent intensity for the CBC simulation (solid line) and the simulation without island but taking into account the profiles given by the presence of the island (dashed line). 110 6.11 Radial autocorrelation functions for the ITG parameters used in this paper in the absence of magnetic island. . . . . . . . . . . . . . 112 6.12 Radial profile through the O-point of the shear of the radial electric field for two island rotation directions: electron diamagnetic direction (dashed line) and ion diamagnetic direction (solid line). The case of an island width W = 6 is considered. . . . . . . . . . . . . . 113 6.13 Poincaré map of the magnetic field perturbed by a tearing-like instability, leading to the presence of m = 2 and m = 3 islands (left panel), which are located on the resonant surfaces q = 2 and q = 3 (right panel). .

Complex systems: from the smallest entity to the global behaviour

Reality is sometimes much more complex than we expect, but what does make real world so complex? Or in other words, what is the definition of complexity? An ant moving on the ground towards the nest exhibits apparently no complexity: in general, its motion will be directed following the shortest path. Nevertheless, the interior of an ant's nest exhibits a much more complex behaviour. It is apparently a messy space where the ants live and this is due to the fact that the motion of each single ant is dictated not only by its function in the ant's nest, but also by the interaction with the closest neighbours, and the whole ant's nest behaves in an apparently unpredictable way This is basically the definition of complex systems: systems made up of many individual differentiated or un-differentiated entities interacting with each other following local (and sometimes also nonlocal) rules, leading to a global behaviour which is difficult to predict or model, characterised by collective phenomena and pattern formation. Given the previous definition, it is obvious that complex systems can have many different natures. One can find them in physics, which is the main scope of the present manuscript, but also in many other fields such as engineering, biology, chemistry, economics or sociology. The fact that the different entities interact with each other not only locally, but also through long distances, introduces the nonlinear and nonlocal behaviours characteristic of these systems.

The fact that complex systems exhibit apparent unpredictable behaviours imply that modelling is required on the route towards the understanding of such systems. Modelling means that we must accept that we are not able to find neither to solve the exact equations dictating the evolution of complex systems. The model can be more or less complete, in the sense that increasing the completeness of the model makes it closer to the real system whose behaviour we want to predict or understand. In general, the more the model is complete, the more its mathematical treatment reveals challenging and sometimes even impossible. More specifically, let us assume the existence of a system represented mathematically by the N-tuple S = {u i } 1<i<N , where u i is the i-th entity of the system, made up of N entities. Let us assume that the system S is fully described by the time evolution of each entity u i , which is in its turn fully described by the time evolution of a vector x i belonging to a finite dimensional space W of dimension d i . Then if we are able to find a model for the system S we can write its evolution by means of the equations

ẋi = F i (x 1 , . . . , x N , t) , 1 ≤ i ≤ N (1.1)
with the initial conditions

x i (t 0 ) = x i,0 , 1 ≤ i ≤ N (1.2)
where F i is a rule to be applied and that depends on time and all the vectors x i , for 1 ≤ i ≤ N . We can give a mathematical sense to the rule F i . Probably the most

CHAPTER 1. MULTI-SCALE PROBLEM

unrealistic but also the most convenient case is that F i be a function

F i ∈ C ∞ .
This representation is unrealistic because by definition nature is discontinuous.

It is nevertheless convenient because this allows one to either perform analytic calculations or use numerical methods to solve the differential equation 1.1. But in general F i is not continuous and can even be a mix of functions, distributions, boolean operators and random processes. In this case, analytic calculations are in general intractable and only performing numerical simulations can the system 1.1 be solved. A schematic representation of what we have in mind when talking about complex systems is given in figure 1.1. The magenta circles represent the entities u i . Each entity can interact with a given number of neighbours

N int ≤ N - 1.
For the schematic view of figure 1.1 we decided that each entity would interact with any other entity with a probability p int = 0.05. The interaction is represented by red arrows and the strength of the interaction is represented by their width. This kind of system can represent for instance the flow of information between all the users in a social network, the bank transfers between all the customers in a financial market or the motion of players on a field during a soccer game. As explained earlier, a complex system is characterised by self-organisation and generation of patterns. A pattern is nothing else but a regularity which allows us to partially characterise and predict the behaviour of the complex sys-

SOME EXAMPLES OF COMPLEX SYSTEMS

tem by means of a function

Q {x i } 1≤i≤N , t : Π i W d i × R → S, where S ⊂ Π i W d i , dim (S) = n.
Mathematically, a pattern can be defined by the couple

P = G, Qj , k j , ω j j∈N
, where Qj , k j , ω j ∈ C n+d j +1 and G is such that

G ({Q j , k j , ω j }) = Q.
The function Q contains some physical information. For instance, the mean velocity of a fluid or the time evolution of the prices in a stoke market. The less is n with respect to i d i , the more reduced the pattern -and therefore the prediction -is. Basically, the function Q represents a dimensionality reduction, i.e. we go from a rather complete description of the complex system to a more restricted one. For example, in the kind of complex systems we are interested in the present dissertation, the distribution of the position of all the particles in space and time (what we will call kinetic description) can be reduced to the velocity of the fluid, density and temperature (what we will call fluid description).

The regularity is reflected by the wave vectors k j and the frequencies ω j , representing regularities in space Π i W d i and time, respectively. One can define the spatio-temporal scales as {δ j , τ j ,

γ j } = k j -1 , ℜ (ω j ) -1 , ℑ (ω j ) -1 .
The weight of each scale to the characterisation and eventual prediction of the system is determined by the amplitude of Q j . Note that the frequency is a complex number: the real part is the frequency characterising the periodicity in time and the imaginary part is the growth (ℑ (ω j ) > 0) or damping (ℑ (ω j ) < 0) rate characterising the growth or damping of a given pattern.

In the definition of the pattern, G is a kind of mathematical transformation. For the sake of clarity, the most common choice for G is the Fourier transform, i.e. G ≡ F, but other choices like the wavelet transform can be made, as we will do in chapter 5.

Some examples of complex systems

Before analysing the complex system that has been at the centre of my research in the last ten years, I would like to emphasize the fact that all complex systems share some common properties, which are therefore found in many fields. I will give hereafter two examples of complex systems very close to our daily lives that can be studied using approaches similar to those that I use in my research.

Modelling the traffic in a round-about or how to get a jam without any bottleneck

In our daily lives we encounter, sometimes without being aware of, many complex systems. One of the most well-known and sometimes hated complex systems is the flow of vehicles in our roads. It is indeed a complex system, for it is composed of many entities (the vehicles) interacting with each other locally or at least at short distances, characterised by a collective behaviour and a selforganisation, making it possible that the local interactions lead to the evolution of CHAPTER 1. MULTI-SCALE PROBLEM the global system with the formation of patterns. Although apparently very technical and complicated to understand, it basically means that because traffic flow is a complex system, we have the chance to experience traffic jams every morning on our way to work. In other words, a traffic jam is nothing else but the signature of a complex system. The analysis of traffic flow has become an extremely important topic in the last decades, especially due to the increasing use of cars and therefore to the increasing impact on the environment. For a detailed overview on the topic, the reader is encouraged to go through Ref. [START_REF] Helbing | Traffic and related self-driven many-particle systems[END_REF] and references therein, in particular the starting works of Biham [START_REF] Biham | Self-organization and a dynamical transition in traffic-flow models[END_REF], Nagel and Schrekenberg [START_REF] Nagel | A cellular automaton model for freeway traffic[END_REF], and Kerner and Konhäuser [START_REF] Boris | Cluster effect in initially homogeneous traffic flow[END_REF]. It has attracted also the work of many scientist who wanted to explain surprising effects like for instance how it is possible that in a road without any bottle-neck a traffic jam is formed. We experience sometimes this phenomenon called phantom traffic jams. These traffic jams are characterised by some vehicles stopping and going in a quasi-periodic way with no apparent reasons. This was for the first time evidenced by an experience where 22 initially uniformly-spaced drivers were asked to drive at a constant speed of 30 km/h inside a round-about of 230 m circumference [START_REF] Sugiyama | Traffic jams without bottlenecks. experimental evidence for the physical mechanism of the formation of a jam[END_REF]. After some time, small fluctuations in the velocity of the vehicles appeared, which grew in time, breaking the flow and leading to a stop-and-go wave propagating in the direction opposite to the motion of vehicles. We provide here a quick explanation. For this purpose we must assume a model for the traffic flow, or in other words, for the behaviour of the drivers, which is the same as giving the rule of Eq. (1.1). We can use different models like the optimal velocity one proposed in [START_REF] Frank | Nonlinear effects in the dynamics of car following[END_REF] using the optimal velocity suggested in [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF][START_REF] Bando | Phenomenological study of dynamical model of traffic flow[END_REF][START_REF] Bando | Analysis of optimal velocity model with explicit delay[END_REF]. The rule can be summarised as follows. When the car in front of us accelerates, we can accelerate, but always keeping in mind the speed limitation v ≤ v max and a safety distance between cars ∆s i = s i -s i+1 ≤ s min . Basically, we accelerate as long as we can. Mathematically, this can be modelled by the following monotonically increasing function [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF]

F i = v max 2 [tanh (∆s i (t) -s min ) + tanh (s min )] , (1.3)
which exhibits a turning point, essential to obtain the traffic jam. This model has been called in the literature optimal velocity model (OVM). Other models exist, but further insight into the whole traffic flow analysis is beyond the scope of the present manuscript. Let us for the moment simply solve Eq. (1.1) using the model given by the Eq. (1.3). We do this for a round-about exhibiting a circumference of 500 m, with 69 cars, each car with a length L c = 4.5m, v max = 2 m/s and s min = 2 m. As a first case, we will consider that all the cars are equivalent. Therefore, the standard deviation of their lengths vanishes. The result is shown in Fig. 1.2. The left panel gives the speed of the cars at the last instant of the simulation. It is observed that there are cars that are completely stopped, whereas other cars are moving at a speed close to the upper limit. The middle panel shows the time evolution of the speed of the first car (car number 2 is ahead of car number 1, car number 3 is ahead of car number 4, and so on). It is observed that the velocity of this car exhibits oscillations in such a way that sometimes the car is completely stopped and right afterwards the car exhibits the maximum velocity. This

SOME EXAMPLES OF COMPLEX SYSTEMS

is clearly seen in the magnified view of the red rectangle area. It can be observed also the that this phenomenon is not related to a low-resolution problem, since each point in that magnified view corresponds to a time step. The right panel shows the same pictures in 2D, i.e. the speed of all the cars as a function of time, exhibiting a self-organisation pattern characterised by perturbation fronts propagating at v ph = -2.7 m/s, i.e. in the direction opposite to the motion of vehicles, as was observed experimentally [START_REF] Sugiyama | Traffic jams without bottlenecks. experimental evidence for the physical mechanism of the formation of a jam[END_REF]. This result might be understood as the outcome of succesive and mutually coupled Stackelberg-like equilibria [START_REF] Von | Marktform und gleichgewicht[END_REF]. Indeed, a leader drives and the follower moves sequentially. This process is repeated, but this time the previous follower becomes now the leader with a new follower who will become the leader afterwards and so on. It can also be observed the presence of different spatio-temporal scales. In this example we have considered only cars with the same length. But in reality, there are different types of vehicles (lorries, cars, motorbikes, bycicles) together with pedestrians, having all different maximum velocities. Such models for the traffic-flow were for the first time introduced in [START_REF] Anthony | Car-following model of multispecies systems of road traffic[END_REF] to analyse the dynamics of cars and trucks interacting with each other.

Other models can be introduced to take into account in a more realistic way the long-distance interactions [START_REF] Delitala | Mathematical modeling of vehicular traffic: a discrete kinetic theory approach[END_REF]. As a result different spatio-temporal scales are at play. In our example, there are three spatial scales: the smallest one or micro-scale (related to the vehicles), the intermediate one or meso-scale (related to the size of the fronts propagating backwards) and the biggest one of macro-scale (the size fo the whole system, in this case the round-about). Similarly, there are also three time scales associated to the motion of the vehicles, the propagation of the fronts and the slow evolution of the whole system. Similar self-organisation patterns with the existence of disparate spatio-temporal scales are found in the complex systems that we analysed in the present manuscript.

Modelling a tumour growth

The previous subsection gives an example of a complex system that is modelled by a rule F i ∈ C ∞ . However, many complex systems cannot be described in such a

CHAPTER 1. MULTI-SCALE PROBLEM

way, especially those whose entities are characterised by discrete (and sometimes binary) states. In that case, the rule F i is not a classical function and one had to abandon the idea of regularity and continuity. This is the starting point of what has been called cellular automata. An example of a complex system modelled by cell automata is the human body, for it is composed of many individual entities called cells, which reproduce by means of the mechanism called mitosis. The cells can be differentiated if they are devoted to performing a very specific task. But a more general classification of cells can be done depending on whether the cells are healthy or not and in particular whether the cells are cancerous or not. The behaviour of the cancerous cells can be modelled as a first approximation using the so-called hallmarks of cancer [START_REF] Hanahan | The hallmarks of cancer[END_REF][START_REF] Hanahan | Hallmarks of cancer: the next generation[END_REF][START_REF] Natalya | The emerging hallmarks of cancer metabolism[END_REF]. A cancerous cell is that which reproduce itself without any control, killing the healthy cells around it and eventually leading the death of the biological tissue, transferring to the contiguous tissues and resulting in the death of the body. The state of a cell can then be binary: it is an healthy cell or it is a cancerous cell. A cell is characterised by the length of the telomere l tel . After a mitosis, the telomere length is reduced by one unit. When the length reaches l tel = 0, the cell dies. Also, a cell that has experienced a mitosis cannot experience a new one before a certain time ∆t mit . When this minimum time is elapsed, the healthy cell reproduces with a probability p. In addition, a healthy cell can reproduce only if there is available place in its neighbourhood, but also a minimum number of neighbours n nb . This limits naturally the growth of the healthy cells. However, there is a probability p mut that a healthy cell, during the mitosis, be mutated into a cancerous cell. A cancerous cell will evade the programmed cell death when its telomere length reaches l tel,max and will reproduce even when there is no place in the neighbourhood. For this purpose, a cancerous cell will kill one of its neighbours and will place a new cancerous cell at its place. Although very simple with respect to the real mechanism behind the evolution of biological tissues, the system that we have just described contains all the ingredients of a complex system and actually, an analytical treatment would be intractable. For this reason, following the work in [START_REF] Robert G Abbott | Simulating the hallmarks of cancer[END_REF][START_REF] Monteagudo | A cellular automaton model for tumor growth simulation[END_REF][START_REF] Santos | Study of cancer hallmarks relevance using a cellular automaton tumor growth model[END_REF][START_REF] Monteagudo | Studying the capability of different cancer hallmarks to initiate tumor growth using a cellular automaton simulation. application in a cancer stem cell context[END_REF][START_REF] Monteagudo | Treatment analysis in a cancer stem cell context using a tumor growth model based on cellular automata[END_REF], we perform a 3D numerical simulation of the cancer growth, starting the simulation with one single healthy cell that is born in the middle of the simulation box. The result is shown in figure 1.3, where the number of healthy (solid blue line) and cancerous (dashed red line) cells are plotted as a function of time for two cases depending on the mutation probability: p mut = 10 -4 (left panel) and p mut = 10 -2 (right panel). Three aspects related to instabilities analysed in the present dissertation can be observed in this system. First, the presence of a threshold for the cancer cells to develop. Second, the existence of an overshoot for the population of healthy cells. This overshoot is also observed in the Ion Temperature Gradient driven turbulence that we analyse in chapter 5. Third, the fact that the cancerous cells act as a predator of healthy cells, which leads to the decrease of healthy cells and eventually the saturation of the tumor.

NUCLEAR FUSION PLASMAS AS COMPLEX SYSTEMS WITH MULTI-SCALE INSTABILITIES

Nuclear fusion plasmas as complex systems with multi-scale instabilities

Many other complex systems can be found in everyday life. Humans thinking, interacting and communicating, the different stock markets in the world or the atmosphere on Earth are some more examples. In this manuscript I focus on a special case of complex system: a magnetised fusion plasma. The research I have carried out during my career centres around the study of plasmas confined inside toroidal devices, called tokamaks, in order to produce energy from nuclear fusion reactions. The nuclear fusion reaction that is envisioned for energy production occurs between two isotopes of hydrogen resulting in alpha particles with 3.5 MeV and neutrons with 14.1 MeV. This reaction is only possible if the hydrogen isotopes reach a sufficiently high energy. Therefore, the plasma needs to be heated up to 10 8 K. Ideally, the hydrogen isotopes should be heated by collisions with the alpha particles produced during the nuclear reactions in such a way that the reaction is self-sustained. This is the reason why the plasma and the fusion products need to be well confined and the confinement is possible by means of an externally applied toroidal magnetic field and a poloidal magnetic field, which constitutes the principle of tokamaks. Nuclear fusion as a source of energy is virtually safe, clean and limitless, and represents therefore the future of energy production at worldwide scale. The way to the steady-state production of this energy is not straightforward and a first step is the realization of the first experimental nuclear fusion reactor. This is the goal of the international ITER project, one of the most ambitious technological and scientific projects for energy production to date. The neutrons, which are not charged, can escape from the plasma and are subsequently used to react with a lithium blanket to produce energy which will be coupled to the electrical network. However, reality is much more complicated than this simple extraction of neutrons. Indeed, when understanding and modelling a plasma we have to face all the problems of a complex system, where nonlinearity, self-organisation, collective behaviour and multi-scale interactions CHAPTER 1. MULTI-SCALE PROBLEM are at play. Complex systems represent an incredible challenge for modellers, because strictly speaking all the scales and interactions must be solved without separation. In addition, a fusion plasma is an extreme environment with high temperatures in the core in order to fulfil the condition for nuclear fusion reactions to occur. However, the plasma is confined inside a vacuum vessel composed of materials which cannot support the same temperatures as those of the core. Therefore, the temperature must be reduced from the core to the edge, resulting in high gradients, not only of temperature, but also of density. The strong gradients lead to instabilities that saturate nonlinearly and are able to induce transport of energy and particles from the core to the edge, limiting therefore the confinement. Although there are many classes of instabilities [START_REF] Mo | Microinstability theory in tokamaks[END_REF] depending for instance on the type of particles, the location in the tokamak, the source of free energy and the electrostatic/electromagnetic nature, in my research career, I have focused on three types of instabilities, characteristic of current and future fusion devices.

Ion-temperature-gradient (ITG) driven instability

Initially evidenced to belong to the group of instabilities due to temperature gradients [START_REF] Coppi | Instabilities due to temperature gradients in complex magnetic field configurations[END_REF], the ITG instability is a Rayleigh-Bénard-like instability developed between layers exhibiting different temperatures and densities. It is similar to the atmospheric instability, which develops when the density gradient is opposed to the gravity. In a tokamak, the counterpart of the gravity is the curvature of the magnetic field and the free energy comes from the gradient of the ion temperature, hence its name. The ITG instability saturates leading to a turbulent behaviour [START_REF] Lee | Theory of ion-temperature-gradient-driven turbulence in tokamaks[END_REF], characterised by eddies which transport energy and particles from the core to the edge of the tokamak, even in the presence of a confining magnetic field. Independently of the nature of turbulent transport in tokamaks driven by ITG modes [START_REF] Ra Úl Sánchez | Nature of transport across sheared zonal flows in electrostatic ion-temperature-gradient gyrokinetic plasma turbulence[END_REF], turbulence must absolutely be controlled, for it limits the confinement of energy and particles, decreasing this way dramatically the performance of the machine. Much effort has been done in this direction, especially within the framework of numerical simulations including kinetic effects. This has resulted in an extremely productive field of research centred around the gyrokinetic simulations of turbulent transport [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]. Within the framework of turbulence driven my temperature gradient, an important role is played by the zoncal flows [START_REF] Diamond | Zonal flows in plasmas -a review[END_REF] and the oscillatory component of the zonal flows, which is referred to as geodesic acoustic modes, or GAMs [START_REF] Winsor | Geodesic acoustic waves in hydromagnetic systems[END_REF], because of their potential role in suppresing turbulence, which might have an impact on the transition towards the high-confinement regimes in tokamaks [START_REF] Gd Conway | Mean and oscillating plasma flows and turbulence interactions across the lh confinement transition[END_REF]. As we will see, GAMs play a crucial role in the present dissertation in the absence and in the presence of ITG-driven electrostatic turbulence, when they are driven by energetic particles.

Energetic particle driven modes (EPMs)

Energetic particles (EP) are ubiquitous in both laboratory and astrophysical plasmas. By definition, they exhibit velocities much larger than the thermal velocity

NUCLEAR FUSION PLASMAS AS COMPLEX SYSTEMS WITH MULTI-SCALE INSTABILITIES

of the bulk plasma, which is characterised by a Maxwellian distribution function. EP, such as the alpha particles, must be sufficiently well confined in order to transfer their energy to the bulk plasma through Coulomb collisions or to ensure the current drive efficiency [START_REF] Heidbrink | The behaviour of fast ions in tokamak experiments[END_REF][START_REF] Se Sharapov | Energetic particle physics in jet[END_REF][START_REF] Sd Pinches | The role of energetic particles in fusion plasmas[END_REF]. Nevertheless, because they exhibit orbits that can be as large as the size of the device, they tend to be de-confined when they approach the wall of the tokamak. In addition, the presence of a substantial population of particles at high energies leads to gradients in phase space, which may result in instabilities called energetic particle modes (EPM) (see for instance [START_REF] Chen | Theory of alfvén waves and energetic particle physics in burning plasmas[END_REF][START_REF] Ww Heidbrink | Basic physics of alfvén instabilities driven by energetic particles in toroidally confined plasmas[END_REF]35,36] and references therein). EPM tend to increase the transport of energetic particles, reducing inevitably the tokamak performance. Therefore, understanding and controlling the EPMs is also of prime importance for the future of ITER. During my career, I have focused on a special class of EPMs, called energetic geodesic acoustic modes (EGAMs), characteristic of toroidal devices, such as tokamaks, but it can be extended to other complex systems. EGAMs are born from the standard geodesic acoustic modes (GAMs) [START_REF] Winsor | Geodesic acoustic waves in hydromagnetic systems[END_REF] that are excited in the presence of a population of energetic particles [START_REF] Fu | Energetic-particle-induced geodesic acoustic mode[END_REF]. As the standard GAMs, EGAMs have the particularity to be axi-symmetric modes, meaning that there is no dependence on the toroidal angle, and they correspond to the oscillating component of the zonal flows. Because of the beneficial impact that zonal flows can have on turbulence, GAMs have received much attention in the last years, and EGAMs have been seen as a way to control turbulence by injection of energetic particles. GAMs are usually stable when only thermal particles exist, but in the presence of energetic particles, they can be driven unstable, leading to the socalled EGAMs. The EGAM instability is a bump-on-tail-like instability, i.e. it is driven by an inversion of the slope of the distribution function in velocity space. EGAMs have been observed in many experimental devices [38,39,[START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF]41], modeled in numerical simulations [42,43,44,45,46,47,48] and explained in analytic theory [START_REF] Fu | Energetic-particle-induced geodesic acoustic mode[END_REF][START_REF] Hl Berk | Fast excitation of egam by nbi[END_REF]50,51,42,1,44,52,[START_REF] Qu | Energetic geodesic acoustic modes associated with two-stream-like instabilities in tokamak plasmas[END_REF]54,55,56]. There has been an increasing interest in these modes mainly due to the rather optimistic possibility of using them as an external knob to suppress turbulence [38]. Before my work, these modes were considered as candidates for the suppression of turbulence and in addition they were assumed to have negligible impact on particle. I show in this thesis not only the evidence that these modes can couple nonlinearly to ITG modes [48] having a deleterious impact on turbulence [43,57,48], but also how they result in anomalous transport of energetic particles even in the absence of turbulence [START_REF] Zarzoso | Particle transport due to energetic-particle-driven geodesic acoustic modes[END_REF].

The tearing mode (TM) instability

The externally applied confining magnetic field can be modified by magnetohydrodynamic (MHD) instabilities, exhibiting different scales and leading to an important phenomenon called magnetic re-connection and subsequent formation of magnetic islands [59,60,61], which strongly limits the confinement. These instabilities exist in both astrophysical and laboratory plasmas [62]. Beacause the tearing mode with the formation of magnetic islands can have have a negative impact on the performance of tokamaks [63], this instability must be detected in

CHAPTER 1. MULTI-SCALE PROBLEM

ITER before it is triggered in order to mitigate it or suppress it. In the tearing instability, the available free energy finds its origin in the radial gradient of the parallel current, as we will see later on in this dissertation. During my career, I have analysed the stability of the tearing modes with gyro-kinetic theory and simulations [START_REF] Zarzoso | Gyro-kinetic theory and global simulations of the collisionless tearing instability: the impact of trapped particles through the magnetic field curvature[END_REF]. I have also developed detailed work on the numerical verification of the implementation of magnetic islands in gyro-kinetic codes [START_REF] Zarzoso | Verification of a magnetic island in gyro-kinetics by comparison with analytic theory[END_REF], together with their impact on ITG-driven turbulence and the radial profiles [START_REF] Zarzoso | Impact of rotating magnetic islands on density profile flattening and turbulent transport[END_REF]. Tearing modes can be of course coupled to energetic particles. As we will see in this dissertation, similarly to EGAMs, tearing modes can have a direct impact on particle trajectories, and especially on energetic particles due to their significant orbit widths compared to that of thermal particles. The transport of energetic particles induced by low-frequency MHD instabilities such as tearing modes has been analysed theoretically in the past [START_REF] Harry E Mynick | Transport of energetic ions by low-n magnetic perturbations[END_REF][START_REF] He Mynick | Stochastic transport of mev ions by low-n magnetic perturbations[END_REF] and also experimentally [START_REF] Em Carolipio | Simulations of beam ion transport during tearing modes in the diii-d tokamak[END_REF][START_REF] García-Mu Ñoz | Ntm induced fast ion losses in asdex upgrade[END_REF]. In addition to the fact that tearing modes can modify the trajectories, energetic particles can also have an impact on the linear and nonlinear stabilities of tearing modes. This was for the first time analysed using a kinetic theory and leading to the puzzling result that energetic particles can stabilize nonlinearly the tearing mode if the energetic ion density peaks just outside the rational surface [START_REF] Chris | Suppression of the tearing mode by energetic ions[END_REF][START_REF] Chris | Suppression of magnetic islands by energetic ions in toroidal plasmas[END_REF]. Later works followed with kinetic simulations coupled to resistive MHD codes showing how energetic ions can have a stabilizing effect at experimentally relevant parameters [START_REF] Takahashi | Kinetic effects of energetic particles on resistive mhd stability[END_REF]. Theoretical analysis proved that the role of energetic ions might be related to the circulating direction [START_REF] Cai | Influence of energetic ions on tearing modes[END_REF], which was later confirmed by hybrid simulations [START_REF] Cai | Hybrid simulation of energetic particle effects on tearing modes in tokamak plasmas[END_REF]. Also energetic ions might be responsible for the generation of an uncompensated current resulting in a nonlinear stabilization [START_REF] Cai | Influence of energetic ions on neoclassical tearing modes[END_REF]. Recent experimental measurements and analyses do not show any definitive phase-space dependent effect on the tearing mode stability, but only a possible diffusive transport of energetic ions [START_REF] William W Heidbrink | The phase-space dependence of fast-ion interaction with tearing modes[END_REF]. Therefore, so far no clear agreement has been reached regarding the mutual interplay between energetic particles and the tearing mode. This is the reason why part of my current research activities are focused on the analysing of the confinement of energetic particles in the presence of tearing modes and how the energetic particles can impact the stability of tearing modes. I develop this area of studies by means of reduced models within the framework of my EUROfusion participation.

The multi-scale problem in my research

All these instabilities are characterised by spatio-temporal scales exhibiting a wide separation. The ITG instabilities exhibit spatial scale lengths of the order of the ion Larmor radius (∼ 1 -2 cm in current tokamaks for thermal ions) and temporal scales of the order of the inverse diamagnetic frequency (∼ 10 -6 -10 -5 s). The spatial scale lengths of EPM are generally related to the orbit width (see chapter 2 for a detailed explanation) of energetic particles exciting the instability (ranging from ∼ 10 cm to the size of the device), whereas the temporal scales are of the order of the transit frequency of energetic particles inside the tokamak (∼ 10 -5 s). Finally, the tearing instability is characterised by two spatial

THE MULTI-SCALE PROBLEM IN MY RESEARCH

scales: close to the resonant surface we find the electron Larmor radius scale and far from the resonant surface the spatial scale is of the order of the size of the machine. This instability leads to magnetic islands of the order of ∼ 10 -20 cm. The tearing modes are characterised by low frequencies, evolving within resistive times of the order of ∼ 1 s and therefore slowly compared to the other instabilities.

Much work has been carried out and published in the literature concerning each of these issues separately, but several difficulties have been encountered. For instance, using high performance computers (HPC) allowing the simulation of such complex systems has been possible only very recently. In addition, energetic particle physics must be analysed in the so-called kinetic framework, where particles are described in velocity and position space, which increases the computational cost of the simulations. Finally, the mutual interplay between these instabilities has been little explored. In particular, the impact of energetic particles on turbulence or the interaction between energetic particles and tearing modes. Therefore, the goal of my past research has been the analysis of the selfconsistent evolution of all these instabilities when they co-exist without any scale separation. This has been an extremely ambitious project, but it has been possible using the state-of-the-art numerical simulations taking advantage of the most recent progress in HPC together with analytic theory based on Hamiltonian formulation.

The remainder of this dissertation is structured as follows. In chapter 2, the reader will be guided through a introduction to the fundamental mathematics and physics necessary for the understanding of the dissertation. The main results concerning the EPM analysed in my research are reported in chapter 3. I quantify in chapter 4 the impact of the EPM on the particle trajectories, providing for the first time evidences of anomalous diffusion. The nonlinear interaction between these modes and ITG dominated turbulence is shown in chapter 5. The theory of tearing modes, their interaction with particles and ITG-dominated turbulence is presented in chapter 6, together with ongoing studies of the stochastic orbits in the presence of tearing modes. I give perspectives and future directions of my research in chapter 7.

Chapter 2

Gyro-kinetic formalism for the description of instabilities in a tokamak

2.2.1

From the equations of motion to the kinetic approach . . .

THE TOKAMAK AND THE TOROIDAL GEOMETRY

Understanding the complex behaviour of a magnetized plasma to achieve the steady-state production of energy by means of nuclear fusion is a hard task. Therefore, this introduction aims at giving a brief description of the fundamentals that will be used throughout this manuscript to model the behaviour of magnetized plasmas.

The tokamak and the toroidal geometry

A plasma is a ionised gas, therefore composed of charged particles which can move in the presence of an electro-magnetic field. It is well-known that a charged particle with a velocity v ⊥ perpendicular to a magnetic field B experiences a circular motion (called cyclotron motion) around the magnetic field line, and the radius of the circle described by the particle is proportional to v 2 ⊥ . Therefore, it is intuitive to think that introducing a magnetic field in a plasma will confine the particles, since they will be constrained to move around the magnetic field lines. Particles with low (resp. high) energy will be closer (resp. further) to (resp. from) the magnetic field. The particles with a velocity v parallel to the magnetic field will simply describe helical trajectories characterised by a parallel motion and a circular motion around the magnetic field lines, i.e. they will follow the magnetic field lines. We can therefore have a cylinder confining the plasma by means of a vacuum magnetic field. Of course, since the confining device must have a finite size, the idea is to connect the end with the beginning of the cylinder, providing therefore a periodic condition in the parallel direction, which becomes the toroidal angle ϕ. Such a geometrical structure is called torus and the magnetic field resulting from this transformation (cylinder→torus) is called toroidal magnetic field. The unit vector in the toroidal direction is noted e ϕ . It can be shown that a device with a purely toroidal magnetic field cannot confine particles. An additional component is required, which is called poloidal component, for it is supported by the unit vector in the poloidal direction e θ . The resulting device, with a magnetic field with toroidal and poloidal components, is called tokamak and is the device around which the present manuscrit is centred. In a tokamak, there are therefore two natural periodicity angles, namely θ and ϕ. This periodicity makes it natural to decompose any quantity (periodic on these variables) onto Fourier series, which results in the so-called poloidal and toroidal numbers, m and n, respectively. The motion of particles around the magnetic field lines introduces a third periodicity angle called gyro-phase, ϕ c , where the c subscript stands for cyclotron.

In the following, the magnetic field of the tokamak is written in the simplified form

B = B 0 R 0 R (ζe θ + e ϕ ) (2.1)
with B 0 the magnetic field at the magnetic axis, R 0 the major radius at the magnetic axis and R the major radius at a given position R = R 0 + r cos θ. The unit vectors in the poloidal and toroidal directions are written as e θ = r∇θ, e ϕ = R∇ϕ, respectively. The magnetic topology is characterised by circular and concentric

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

flux-surfaces, which makes it natural to introduce the variable r, called minor radius. The variable r indicates the distance from the magnetic axis of the tokamak to the magnetic flux surface and therefore identifies in a unique way the magnetic surface. In the expression of the magnetic field, we have introduced the parameter

ζ = r qR 0 (2.2)
where q is the so-called safety factor, which is a function of r and measures the helicity of the magnetic field lines belonging to the magnetic surface at the distance r. The safety factor is extremely important for the operation of a tokamak and introduces necessary conditions for the development of some of the instabilities analysed in this manuscript. To give an idea of the importance of the safety factor, it is important to realise that for a rational (including integer) numbers it quantifies the number of toroidal turns a magnetic field lines must complete before completing a poloidal turn. If the safety factor is rational (or integer), it can be written, without any loss of generality, as q = m/n. This means that the magnetic field will close on itself after having completed m toroidal turns and n poloidal turns. This means that any perturbation characterised by these mode numbers m and n introduced on a magnetic field line, will feedback itself and can therefore grow, leading to an instability, as we will see later. Another important point of the safety factor is the fact that if q on a given magnetic surface characterised by the minor radius r is an irrational number, the magnetic field lines on that magnetic surface will never close on themselves, which means that a magnetic field can fill ergodically the magnetic surface. We can introduce another coordinate to measure the distance to the magnetic axis. This coordinate is the poloidal flux ψ, defined as

ψ (r) = r 0 r ′ q (r ′ ) dr ′ (2.3)
where r and r ′ are normalized distances. Note that with this definition dψ/dr > 0.

For circular concentric flux surfaces, working with ψ and r is almost identical, since a magnetic surface is unambigously determined by either ψ or r. However, when working with more complicated geometries, the magnetic flux surface is the right coordinate to measure the distance to the magnetic axis, since a magnetic surface cannot be characterised unambigously by r, but only by ψ. In the same way the factor factor describes the helicity of the magnetic field lines, the ratio r/R 0 represents how far from a cylinder the tokamak is. Strictly speaking, the limit case of a cylinder is found only when r/R 0 → 0. In practice, a tokamak is delimited by the last closed flux surface (LCFS), inside which all the magnetic field lines are closed and outside which the magnetic field lines end at the either limiter or divertor (depending on the configuration). The minor radius of the LCFS is called a. The quantity ǫ = a/R 0 is called inverse aspect ratio, and it is considered to be a small number, i.e. ǫ ≪ 1. We will not discuss here all the basics on geometry that are required in this manuscript. For further detail, the reader can refer to the Appendix. The described toroidal geometry, which is used for the analysis of instabilities in a toka-

PARTICLE TRAJECTORIES

mak within the framework of the present manuscript, is represented in figure 2.1. For the sake of clarity only two magnetic surfaces are represented here by blue and black surfaces. It is clearly observed that the magnetic surfaces are concentric. The coordinate system is also represented, with the minor radius denoting in this case the black surface.

Particle trajectories

From the equations of motion to the kinetic approach

Without any loss of generality, we can say that the plasma is composed of N particles, belonging to n s species. Strictly speaking, this is a complex system as the one described schematically in figure 1.1 and the set of equations governing the evolution of the plasma is therefore the system 1.1 with the initial conditions 1.2, where x i = (r i , v i ) T , with r i the spatial position and v i the velocity of the i-th particle, and the rule F i is simply

F i = v i eZ i m i (E (r i , t) + v i × B (r i , t))
(2.4)

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

where the second term is nothing else but the Lorentz force. Here, E is the electric field, deriving from an electrostatic potential φ and a magnetic vector potential A, i.e. E = -∇φ -∂ t A and B is the magnetic field, derived as B = ∇ × A. In the previous expression, e is the elementary charge, Z i is the charge number of the i-th particle and m i is its mass. The electric and magnetic field are obtained from Maxwell's equations together with the constitutive relations of the medium. Eq. 1.1 using the Lorentz force is equivalent to the Hamilton's equations

q = ∂H ∂p (2.5a) ṗ = - ∂H ∂q (2.5b)
where q and p are the canonically conjugated variables for position and momentum, respectively Ideally, one could solve the Hamilton's equations numerically for each particle inside a tokamak using the self-consistent Hamiltonian, as we have done for the examples of complex systems in the previous chapter. Nevertheless, in order to get accurate predictions of the behaviour of a plasma, one should integrate the equations of motion for ∼ 10 23 particles, which remains so far numerically prohibited. This is the reason why we use in the present manuscript a statistical approach, alos called a kinetic approach. We are not interested in knowing the exact motion of a particle, we are interested in knowing the probability for a particles to be in the volume q+dq, with a canonical momentum p+dp within a time window t + dt. This probability is given by the distribution function F (q, p, t). In the absence of collisions, the distribution function is conserved along the trajectories and therefore satifies the Vlasov equation

dF dt = ∂F ∂t -[H, F ] = 0 (2.6)
where [A, B] is the Poisson bracket between A and B, defined as

[A, B] = ∂ q A • ∂ p B -∂ p A • ∂ q B.

Action-angle variables and motion invariants

One can take advantage of the fact that the motion of a charged particle in a tokamak is quasi-periodic around the magnetic field lines and in the poloidal and toroidal directions to define the set of action-angle variables (α, J). These variables are canonically conjugated and satisfy therefore Hamilton's equations

J = - ∂H ∂α (2.7a) α = ∂H ∂J (2.7b)

PARTICLE TRAJECTORIES

In the absence of any perturbation and assuming that particles do not interact with each other, the motion of charged particles inside a tokamak can be decomposed into three motions: the cyclotron, poloidal and toroidal motions, characterised by the cyclotron frequency Ω 1 , the poloidal frequency Ω 2 and the toroidal frequency Ω 3 , respectively. These are called the eigen-frequencies and their expressions are derived in detail in Appendix. These frequencies describe the motion of all the particles in a tokamak in the absence of any perturbation. In the following, we will define the equilibrium as H eq ≡ dαH, i.e. the equilibrium does not depend on the angular variables. This implies in particular that the equilibrium satisfies the equations

J = 0 (2.8a
)

Ω = ∂H eq ∂J (2.8b)
in other words, the actions are the motion invariants associated to the three eigenfrequencies. Using the invariance of J in the absence of perturbations, important properties of the trajectories in a tokamak can be derived, but for this purpose one needs to determine first of all the motion invariants. The first one, J 1 is based on the so-called adiabatic limit, which represents the fact that the length scale of the equilibrium is much longer than the length scale of the cyclotron motion, namely the Larmor radius ρ, and also on the fact the the characteristic time scale of the equilibrium is much longer than the characteristic time scale of the cyclotron motion, namely Ω -1 1 . Mathematically, this is represented by

Ω 1 ≫ ∂ log B eq ∂t (2.9a) ρ ≪ d log B eq dr (2.9b)
Within the adiabatic limit and using the Euler-Lagrange equation, it is straightforward to show that

J 1 = - mµ eZ (2.10)
is the first motion invariant, with µ = mv 2 ⊥ / (2B) the magnetic moment. Since the mass and the charge are usually conserved, the magnetic moment is also a motion invariant. It is to be noted that J 1 is not an exact invariant, but only an invariant within the adiabatic limit. Therefore, µ is also called adiabatic invariant.

The third motion invariant J 3 is related to the toroidal geometry. It is possible to show that in the absence of any perturbation and if the equilibrium magnetic field does not depend on the toroidal angle, which is the case for a tokamak, the toroidal canonical momentum

J 3 ≡ P ϕ = -eZψ + m B ϕ B v (2.11)
is an exact invariant. In this expression, B ϕ is the co-variant component of the magnetic field in the toroidal direction.

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

The second motion invariant does not have a very useful expression

J 2 = Φ (J 3 ) + 1 2π 2π 0 mv ds (2.12)
where Φ is a periodic function of J 3 and s is the parallel coordinate. For our purposes, it is more convenient to use the kinetic energy instead. It is easy to prove that the kinetic energy E = mv 2 /2 + µB, must be conserved in the absence of any perturbation.

Passing and trapped particles: orbit width

We will use the conservation of the kinetic energy and the toroidal canonical momentum to derive two characteristics of the trajectories in a tokamak. The first one is the fact that there are two types of particles. Depending on the topology of their trajectories they are called passing (or circulating) and trapped particles.

The second one is the fact that the orbits do not lie on the flux surfaces that we have already defined. This introduces a length scale that is usually called orbit width.

As a first approximation, the separation between passing and trapped particles can be made based on the conservation of the kinetic energy. For a particle starting its motion at θ = 0 and a given radial position r, if its initial parallel and perpendicular velocities satisfy v /v ⊥ < √ 2ε, the particle will not complete a poloidal turn. Instead, it will bounce back at a poloidal angle θ 0 where its parallel velocity vanishes and then changes its sign. This particle is called trapped particle and its trajectory is usually referred to as banana, since its projection onto the poloidal cross-section has a banana shape. On the other hand, if v /v ⊥ > √ 2ε, the particle will complete a poloidal turn and its parallel velocity will never vanish. This particle is called passing or circulating particle.

Let us now use the invariance of the toroidal canonical momentum for a trapped particle to determine the width of the banana. For this purpose, we apply its conservation between the point (r + δ b , θ = 0) and the point (r -δ b , θ = 0). This invariance reads

-eZψ (r + δ b ) + mR (r + δ b ) v ϕ = -eZψ (r -δ b ) -mR (r -δ b ) v ϕ (2.13)
where we have used the fact that v ϕ (r + δ b , θ = 0) = -v ϕ (r -δ b , θ = 0). Taylorexpanding the poloidal flux and using the relation between parallel and perpendicular velocities at the boundary trapped-passing, one obtains the ordering for the width of the banana

δ b ∼ q √ ǫ ρ (2.14)
Similarly, one can calculate the deviation of a passing particle orbit with respect to the flux surface, resulting in the same ordering δ c ∼ qρ/ √ ǫ.

Examples of the 3D view of the trajectories for passing and trapped particles are given in figures 2.2a and 2.2b, respectively. Their projections onto the poloidal cross-section are provided in figures 2.2c and 2.2d, respectively. The top panels represent the 3D views, whereas the bottom panels represent the projection onto the poloidal cross-section.

PARTICLE TRAJECTORIES

Relationship between the angles α and the tokamak coordinates

What is important to realise is that the frequencies Ω 2 and Ω 3 have the same physical meaning for both passing and trapped particles: they represent the frequencies in the poloidal and toroidal directions, respectively. Nevertheless, their names are different depending on the trajectory they describe. The frequency Ω 2 is called transit frequency for a passing particle and bounce frequency for a trapped particle. The frequency Ω 3 is the toroidal frequency for passing particles, whereas for trapped particles it is divided into to frequencies: one which is of the same order as the bounce frequency and one which is one order less than the bounce frequency. The latter one is called precession frequency. In general, we have the ordering

Ω 1 ≫ Ω 2 ≫ Ω 3 (2.15)
Strictly speaking, these frequencies correspond to the periodic motion of par-

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

ticles in the directions α 2 and α 3 , but these angles do not correspond to the angles θ and ϕ that we use as periodic coordinates for the tokamak. In other words, the way particles see the periodicity is not the same as the way we see the periodicity.

Trapped particles

It is clear that the trapped particle illustrated in figure 2.2d does not have a periodic motion in θ, because -θ 0 ≤ θ ≤ θ 0 . The periodicity is reflected by the angle α 2 that is given in figure, although this representation is not exact either. This implies that there is the relation θ = f (α 2 ), with f a periodic function. For deeply trapped particles, it can be shown that θ = θ 0 sin α 2 . Similarly, the toroidal angle ϕ is equivalent to the third angle α 3 , but modulated by a periodic function of α 2 , i.e. ϕ = α 3 + qf (α 2 ). Another variable that is essential to express in terms of α is the parallel velocity, for it is not a motion invariant and therefore depends on the position of the particle. For trapped particles, we can approximate the parallel velocity by v = v ,0 cos α 2 , where v ,0 ≈ (2/m) µB 0 ǫ (1 -cos θ 0 ). Finally, the radial position of a trapped particle is expressed in terms of the radial position of the reference magnetic surface r s and the angle α 2 as r = r s + δ b cos α 2 .

Passing particles

The fact that θ = α 2 is not so straightforward for passing particles. Indeed, by having a look at figure 2.2c, it might be tempting to say that θ is the periodicity angle. In fact, this is almost true for passing particles with low energy (thermal particles), because the deviation of their trajectories with respect to the magnetic flux surfaces is not significant. To understand this one can simple realise the small width of the banana in figure 2.2d. A simple analysis can reveal that for passing particles θ = α 2 -σ (δ c /r s ) sin α 2 , where σ = v /|v | and r s is the radial position of the reference magnetic surface. It is to be noted that this relation is valid only if δ c < r s . Similarly, ϕ = α 3 -qσ (δ c /r s ) sin α 2 . The parallel velocity reads v = v ,0 -µB 0 /v ,0 ǫ sin (α 2 /2). Finally, the radial position is expressed as

r = r s + δ c cos α 2 .
For energetic passing particles, the difference between the tokamak angles and the periodicity angles is more dramatic and significant than for thermal passing particles. Also, for energetic particles (both trapped and passing) the deviation of their orbits with respect to the magnetic surface is more important.

Exact linear kinetic response of particles to the electro-magnetic field of a tokamak

The Vlasov equation can be solved analytically when linearizing around an equilibrium. For this purpose, we decompose the total Hamiltonian into an equilibrium and a perturbation H = H eq + δH (2.16)

EXACT LINEAR KINETIC RESPONSE OF PARTICLES TO THE ELECTRO-MAGNETIC FIELD OF A TOKAMAK

In the same way, we can decompose the distribution function into an equilibrium and a perturbation F = F eq + δF (2.17)

with the equilibrium distribution function defined as

F eq ≡ dαF (2.18)
In practice, particles are embedded in an equilibrium magnetic field B eq , with toroidal and poloidal components, and their response provides an electro-magnetic field composed of an electrostatic potential φ and a perturbed magnetic field δB. In this manuscript, we will assume that the perturbed magnetic field is perpendicular to the equilibrium magnetic field, in other words, the perturbed magnetic field derives from a parallel vector potential δA .

Owing to the periodicity of the motion, we can project the perturbed Hamiltonian and distribution function onto a Fourier basis

δH = n,ω δH nω e i(n•α-ωt) (2.19a) δF = n,ω δF nω e i(n•α-ωt) (2.19b)
Linearizing the Vlasov equation results in the exact linear solution

δF = n,ω δF nω e i(n•α-ωt) = - n,ω n • ∂ J F eq ω -n • Ω δH nω e i(n•α-ωt) (2.20)
Note that this has been obtained independently of the Maxwell's equations and for a general equilibrium. We have so far made the assumptions that the equilibrium does not depend on the angular variables and that a kinetic description of a plasma is accurate enough.

Gyro-kinetic approach: when the fastest motion can be averaged out

In order to perform analytic calculations, we need to substitute the expressions of the motion invariants in the numerator of the exact solution. This is almost straightforward for the first and the third invariants. Indeed, J 1 is trivially shown to be proportional to the magnetic moment, by means of the perpendicular energy µB and the cyclotron frequency Ω 1 , so that J 1 = -m s (eZ s ) -1 µ. The third motion invariant is simply the toroidal canonical momentum J 3 = P ϕ . However, J 2 does not have a trivial expression. Therefore, we make use of another motion invariant in the absence of any perturbation: the equilibrium Hamiltonian H eq , which is also called the kinetic energy. Therefore, we perform a change of variable J → I = (J 1 , H eq , J 3 ). The numerator can therefore be written as

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES n • ∂ J F eq = n • ∂ J I • ∂ I F eq .
By doing this, it is easy to rewrite the exact solution as follows

δF nω = - n 1 ∂ J 1 F eq + n • Ω∂ Heq F eq + n 3 ∂ J 3 F eq ω -n • Ω δH nω (2.21)
One can now sum and substract the quantity ω∂ Heq F eq to decompose the exact linear solution into

δF nω = ∂ Heq F eq δH nω - n 1 ∂ J 1 F eq + n 3 ∂ J 3 F eq + ω∂ Heq F eq ω -n • Ω δH nω (2.22)
The inverse Fourier transform gives the particle distribution function

δF = ∂ Heq F eq δH - n,ω n 1 ∂ J 1 F eq + n 3 ∂ J 3 F eq + ω∂ Heq F eq ω -n • Ω δH nω e i(n•α-ωt) (2.23)
We can now make use of the ordering Ω 1 ≫ Ω 2 , Ω 3 to further decompose the third term into two terms, one containing the derivative of the equilibrium distribution function with respect to the magnetic moment, and one containing the resonance

-B -1 ∂ µ F eq n,ω,n 1 =0 δH nω e i(n•α-ωt) - n,ω,n 1 =0 n 3 ∂ J 3 F eq + ω∂ Heq F eq ω -n • Ω δH nω e i(n•α-ωt)
(2.24) The last term contains the sum over all the wavenumber satisfying n 1 = 0. Mathematically, this is equivalent to performing an average over the angle α 1 , which is called a gyro-average. Physically, this corresponds to averaging out the fastest motion of the particle, owing to the time-scale separation between the cyclotron motion and the other two motions. The term proportional to ∂ µ F eq contains all the n 1 = 0 components, which can be obtained by taking all the n 1 components and substracting the n 1 = 0 ones. We can decompose the perturbed distribution function into adiabatic and non adiabatic responses, respectively

δF = δF ad + G (2.25)
where

δF ad = ∂ Heq F eq -B -1 ∂ µ F eq (1 -J 0 •) δH (2.26a) G = - n,ω n 3 ∂ J 3 F eq + ω∂ Heq F eq ω -n • Ω (J 0 • δH) n,ω e i(n•α-ωt) (2.26b)
with J 0 • the gyro-average operator. The adiabatic part corresponds to the response of particles proportional to the perturbed Hamiltonian. The non adiabatic part contains all the resonances. It is to be noted that expression 2.25 is still the exact linear solution of the Vlasov equation, within the ordering Ω 1 ≫ Ω 2 , Ω 3 .

EXACT LINEAR KINETIC RESPONSE OF PARTICLES TO THE ELECTRO-MAGNETIC FIELD OF A TOKAMAK

Existence of multiple resonances in the reduced 5D phasespace

The linear interaction between particles and one single (m, n, ω) mode is not necessarily reduced to one single resonance. In this subsection and the next one we show how a single helicity mode can lead to an ensemble of resonances in phasespace which can eventually overlap. For this purpose, we assume the ordering

Ω 1 ≫ Ω 2 , Ω 3
, ω so that we can reduce the resonance to the simple condition

ω -n 2 Ω 2 -n 3 Ω 3 = 0 (2.27)
It can be shown [START_REF] Zarzoso | Gyro-kinetic theory and global simulations of the collisionless tearing instability: the impact of trapped particles through the magnetic field curvature[END_REF][START_REF] Nasr | Fluid and kinetic descriptions of the mutual interaction between tearing modes and thermal and energetic particles in tokamak plasmas[END_REF] that the second and third frequencies are written, respectively, as

Ω 2 = Ω c,th ρ * ǫ a σ 1 q E -(1 -ǫ) µB 0 π √ 2 K -1 (κ) for κ < 1 κK -1 (1/κ) for κ > 1 (2.28a
)

Ω 3 = Ω c,th ρ * ǫ a Ω t Ω c,th qE ǫ Ωd + δ passing q (r) Ω 2 (2.28b)
where Ω c,th is the cyclotron frequency of thermal particles, Ω t = v th /R 0 , ρ * is the thermal Larmor radius normalised to the minor radius of the tokamak (ρ * = ρ/a), ǫ a is the inverse aspect ratio (ǫ a = a/R 0 , with R 0 the major radius measured at the magnetic axis), q is the safety factor, B 0 is the magnetic field at the magnetic axis, E = v2 /2 with v the velocity normalised to the thermal velocity v th , µ is the magnetic moment, σ is the sign of the parallel velocity, K is the complete elliptic integral of the first kind, ǫ = r/R 0 , which can be expressed in terms of the poloidal flux ψ of the reference magnetic surface as

ψ = r 0 r ′ B 0 q (r ′ ) dr ′ (2.29)
The symbol δ passing is 1 for passing particles and 0 otherwise and κ is the trapping parameter, defined as

κ 2 = 2ǫΛ 1 -(1 -ǫ) Λ (2.30)
with Λ = µB 0 /E the invariant pitch angle. The particles are classified in terms of κ as passing (κ < 1) or trapped (κ > 1). In Eq. 2.28b, the normalised frequency Ωd is expressed in terms of the complete elliptic integrals of the first and second kind as follows [START_REF] Zarzoso | Gyro-kinetic theory and global simulations of the collisionless tearing instability: the impact of trapped particles through the magnetic field curvature[END_REF][START_REF] Nasr | Fluid and kinetic descriptions of the mutual interaction between tearing modes and thermal and energetic particles in tokamak plasmas[END_REF]

Ωd = 4ǫ + κ 2 (1 -2ǫ) 2ǫ + (1 -ǫ) κ 2 1 + 2 κ 2 E K -1 + 4s κ 2 E K -π 2 √ 1-κ 2 K , for κ < 1 2 K E -1 + 4s K E + 1 κ 2 -1 for κ > 1 (2.31)
where the argument of the elliptic integrals for passing particles is κ 2 , whereas the argument is 1/κ 2 for trapped particles. Therefore, the second and third frequencies normalised to the thermal cyclotron frequency are functions of E, ψ, µ . The resonances can be obtained in this 3D space solving the equation for n 2 , n 3 ∈ Z, with the ˆsymbol representing the normalised frequencies to the thermal cyclotron frequency. This is nothing else but the equation of 2D surfaces parametrised by n 2 and n 3 in a 3D space. Therefore, for a given µ, the surfaces are reduced to lines in ψ, E space. In the present dissertation, we focus on n 3 = 0 mode in chapters 3 and ?? and low finite n 3 = 0 in chapter 6. The lines where resonances occur are represented for the case n 3 = 0 (corresponding to EGAMs) in figure 2.4 for two values of magnetic moment: µ = 8 (left) and µ = 16 (right). These two values correspond to particles further from and closer to the trapping-passing boundary, respectively. The resonances corresponding to different values of n 2 are represented in colours and conveniently labeled in the figure. The dashed black line represents the separation between the passing (right region) and the trapped (left region) particles. To produce these figures and only as an example we have imposed the frequency obtained from simulations in chapter 3 [START_REF] Zarzoso | Particle transport due to energetic-particle-driven geodesic acoustic modes[END_REF] and a parabolic q-profile with q(0) = 1.1, q(r = a) = 3.5 and a derivative q ′ (r = 0) = 5 • 10 -4 ρ -1 th (as plotted in figure 2.3). It is observed that the resonances tend to accumulate when approaching the passing-trapping boundary. Since µ is conserved in the adiabatic limit, the motion of particles occurs in the plane ψ, E . However, for the axi-symmetric case of EGAMs (n 3 = 0) the particle is not free to move everywhere in that plane, since the restriction to axi-symmetric modes obliges particles to conserve the toroidal canonical momentum P ϕ .

ω -n 2 Ω2 -n 3 Ω3 = 0 (2.
For the non axi-symmetric case of low n 3 wave numbers (like the tearing mode analysed in chapter 6) we use the same q-profile and a frequency of the order of the one obtained from simulations [START_REF] William A Hornsby | The non-linear evolution of the tearing mode in electromagnetic turbulence using gyrokinetic simulations[END_REF]. Figure ?? depicts this situation for n 3 = -1 and different values of n 2 . These figures show the fact that a single helicity mode will result in several resonances. Seminal work for the tearing mode was done in [START_REF] Harry E Mynick | Transport of energetic ions by low-n magnetic perturbations[END_REF] proposing the existence of these resonances in real space (not in phase-space), leading potentially to the so-called orbit stochastization in the presence of significant resonance overlap in the radial direction satisfying the Chirikov criterion [START_REF] He Mynick | Stochastic transport of mev ions by low-n magnetic perturbations[END_REF]. In this subsection, we have shown that these resonances exist (for both non anisymmetric and axi-symmetric modes) in phase-space and we will provide numerical evidences of the induced chaotic transport in chapters 4 (for the EGAMs) and in chapter 6 for the tearing mode. In order to determine whether an overlap in phase-space is possible or not we must calculate the Fourier mode in phase space, i.e. the (n 2 , n 3 , ω) mode, for both electrostatic and magnetic perturbations.

EXACT LINEAR KINETIC RESPONSE OF PARTICLES TO THE ELECTRO-MAGNETIC FIELD OF A TOKAMAK

Determining the (n, ω) Fourier mode

In practice, the analysis of the instabilities will be done by introducing a perturbation φ or δA . These perturbations are periodic in the tokamak geometry, i.e. their periodic angles are θ and ϕ. Therefore, it is possible to introduce a (m, n) mode. However, as explained earlier, particles do not see the perturbation in the same way as we do. We must find a way to express the perturbation as a funcion of (α 2 , α 3 ), which are the coordinates of periodicity seen by particles. The analyses for an electrostatic and a purely magnetic perturbation are very similar, but an essential difference appears in principle because of the introduction of the parallel velocity in the Hamiltonian in the presence of a magnetic prturbation δA (δH = -eZv A ). This difference is due to the fact that strictly speaking the parallel velocity is not a motion invariant and therefore depends on the angles, which introduces an additional dependence on α 2 . In this section, we calculate the Fourier modes for passing particles, anticipating the fact that we want to study the transport of passing particles well confined in the core that will be magnetically trapped and eventually lost. For the magnetic case, strictly speaking we must introduce also the angular dependence of the parallel velocity. However, this introduces higher order dependences in the amplitude of the perturbation. We can therefore assume a constant parallel velocity and perform the calcula- tion for both electrostatic and magnetic perturbations. We therefore calculate the (n, ω) Fourier mode for H. The poloidal and toroidal angles for passing particles are written, respectively, as

θ = α 2 + Θ (α 2 ) (2.33a) ϕ ≈ α 3 + q (r) Θ (α 2 ) (2.33b)
where Θ is a periodic function of α 2 . We can approximate, for the sake of simplicity, Θ (α 2 ) ≈ θ b sin α 2 , with θ b the poloidal angle corresponding to the position where the particle is at the maximum Z (vertical coordinate). This angle is, for passing particles

θ b ≈ tan -1 r qρ th (2.34)
Since the poloidal flux is not a conserved quantity, a convenient way to represent the vector potential is in terms of the poloidal flux of the reference magnetic surface ψ, which is an invariant of the unperturbed trajectories. The relation between ψ and ψ is done through a periodic function of α 2

ψ = ψ + ψ (J, α 2 ) (2.35)
where we approximate again ψ (J, α 2 ) ≈ ψ0 cos α 2 , with ψ0 the deviation at θ = 0 from the reference magnetic surface, which is a function of the magnetic moment and the safety factor for passing particles, ψ0 ∼ q √ µ. Therefore, the potential can be written as

H = e i(mα 2 +nα 3 -ωt) e (m+nq)θ b sin α 2 + ψ0 cos α 2 ∂ ψ H ψ (2.36)
For the electrostatic modes that we consider in this dissertation, the main component comes from the (m, n) = (0, 0) mode. Therefore, we have φ = e ψ0 cos α 2 ∂ ψ φ 00 ψ (2.37)

COUPLING VLASOV AND MAXWELL'S EQUATIONS

Using the identity

e ix cos α = p i p J p (x) e ipα (2.38)
we can write the mode (n 2 , n 3 = 0, ω) as

φ n 2 ,n 3 =0,ω = i n 2 J n 2 ψ0 k ψ φ 00 (ψ) (2.39)
where the argument of the Bessel function must be understood as a differential operator. For magnetic perturbations with finite (m, n), the ψ0 term can be neglected compared to the term (m + nq) θ b , which leads to

A ,n 2 ,n 3 ,ω = i n 2 J n 2 ((m + nq) θ b ) A ,0 (2.40)
with A ,0 the amplitude of the magnetic potential. Therefore, not only have we shown that several n 2 = 0 resonances exist in phase-space, but also we have proved that the Fourier mode of a single helicity mode is not limited to one couple (n 2 , n 3 ), which can potentially lead to resonance overlap.

Coupling Vlasov and Maxwell's equations

In order to provide the response of particles in the presence of a self-consistent electro-magnetic field, the expression 2.25 has to be injected into the Maxwell's equations. In this thesis, we consider three types of modes: ion temperature gradient modes, energetic geodesic acoustic modes and tearing modes. Although strictly speaking all these isntabilities have to be described within an electromagnetic framework, the two first can be understood using an electrostatic approach and the third one can be as a first step analysed using a purely magnetic formalism.

The main ingredients of the instabilities: resonance and gradient of the equilibrium

The Maxwell's equations introduce integrals in velocity space of the response of particles. The third term in the right-hand side of expression 2.25 leads therefore to a resonance in phase-space. Although very commonly used, the concept of resonance is sometimes misunderstood. Mathematically, the resonance provides the condition for particles to have the same velocity as the wave they interact with. But it is important to stress here that the instability (or the damping) is not due to those resonating particles, but to the nearly resonating ones. In other words, the instability or the damping occurs because particles are in the neighbourhood of the resonance and they synchronise with the wave. Particles having more energy than the wave will be slowed-down and they will transfer part of their energy to the wave. Conversely, particles having less energy than the wave will be accelerated and they will take energy from the wave. Our intuition makes us believe

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

that when there are more particles with more energy than the wave than particles with less energy than the wave, there will be a net transfer of energy from the particles to the wave and an instability will be triggered. On the other hand, when there are less particles having more energy than the wave, the wave will be damped by the plasma. This translates into the criterion for the instability: if the derivative of the distribution function evaluated at the resonance is positive, the mode is excited and damped otherwise. It is to be noted that the gradient of the distribution function appears in the numerator of the expression 2.20, which comes from the linearisation of the Poisson bracket. Therefore, the stability criterion finds its origin in the Poisson bracket of the Vlasov equation. The numerator of 2.25 indicates that the instability comes actually from the gradients of the distribution function with respect to J 3 and H eq . In particular, for an instability with frequency ω and wavenumber (n 2 , n 3 ) to develop, the following condition must hold: n 3 ∂ J 3 F eq + ω∂ Heq F eq > 0 evaluated at the position in phase-space given by the equation ω -n 2 Ω 2 -n 3 Ω 3 = 0. In Appendix the expressions for the eigenfrequencies are given and they depend on the class of particles, i.e. on whether they are passing or trapped. It is important to realise here that the drive can come from the derivative of the equilibrium distribution function with respect to either J 3 or H eq or both. In the present manuscript we will consider that each instability can come from one partial derivative of the equilibrium. Strictly speaking, the equilibrium must depend only on J, i.e. F eq = F eq (J 1 , J 2 , J 3 ). However, in order to perform analytic calculations, this dependence will sometimes be relaxed and we will choose an analytic equilibrium which does not depend on the motion invariants.

For the case of EGAMs, which are axisymmetric modes, the only remaining term in the numerator is ω∂ Heq F eq . The partial derivative must be understood as keeping the other motion invariants constant. Therefore, assuming a constant magnetic field, it is straightforward to write ω∂ Heq F eq = ω mv ∂F eq ∂v (2.41)

The drive comes therefore from the derivative of the equilibrium distribution function with respect to the parallel velocity, evaluated at ω -n 2 Ω 2 = 0.

For the case of ITG modes, we will assume that the equilibrium depends only on H eq and J 3 as

F eq = n eq (J 3 ) (2πT eq (J 3 ) /m) 3/2 e - Heq Teq(J 3 ) (2.42)
what is usually referred to as canonical Maxwellian. The dependence on J 3 introduces the radial dependence, since J 3 ∝ -eZψ for low parallel velocity. The numerator of the resonant response therefore reads

n 3 ∂ J 3 F eq + ω∂ Heq F eq = - F eq T eq (ω -ω ⋆g ) (2.43)
where ω ⋆g = n 3 T eq ∂ J 3 log F eq is the generalised diamagnetic frequency. In the case of the equilibrium 2.42, it contains the gradients of the equilibrium density and temperature profiles

ω ⋆g = ω ⋆ ≡ ω ⋆n + H eq T eq - 3 2 ω ⋆T (2.44)
where ω ⋆n = n 3 T eq d log n eq dJ 3 (2.45a)

ω ⋆T = n 3 T eq d log T eq dJ 3 (2.45b)
In this case it is clear that the drive for the instability comes from the gradients in the radial direction of the equilibrium distribution function.

Finally, for the case of tearing modes (TM), we will consider only electrons, characterised by an equilibrium with an electron parallel flow u e (J 3 ), such that

F eq = n eq (J 3 ) (2πT eq (J 3 ) /m) 3/2 e
-Heq-mv ue(J 3 )

Teq(J 3 ) (2.46)
where the parallel velocity v is built in such a way that it is a motion invariant. In this case, the numerator of the resonant response is given by the expression 2.43, where the generalised diamagnetic frequency reads

ω ⋆g = ω ⋆ + ω ⋆u (2.47)
where ω ⋆u = n 3 mv du e dJ 3 (2.48) represents the gradient of the electron flow, which, as we will see later in this manuscript, results in the drive of the tearing mode.

Ampère's law

Ampère's law is written as

∇ × B = µ 0 j (2.49)
where µ 0 is the vacuum permeability, we assume that all the current comes from electrons, for their mass is much less than that of ions. After decomposing the magnetic field into equilibrium and perturbed parts, writing B = ∇ × A and projecting onto the parallel direction, we can write Ampère's law as

∇ 2 ⊥ A = -µ 0 j (2.

50)

where A = A • b eq and j = j • b eq , with b eq = B eq / B eq . The perturbed parallel current is calculated as j = eZ e d 3 p v ,eq δF + δv F eq (2.51)

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

where the perturbed parallel velocity is expressed in terms of the parallel vector potential δv = -eZ e m e δA (2.52)

Therefore, the perturbed parallel current reads j = eZ e d 3 p v ,eq δF -eZ e m e δA F eq (2.53)

Passing from (x, p) to action-angle variables (α, J) at this point requires to pass from a three-dimensional to a six-dimensional space. For this purpose we use the δ-Dirac function such that j for electrons is calculated as a 6D integration.

j (x, t) = eZ e d 3 pd 3 x ′ δ (x ′ -x) v ,eq δF - eZ e m e δA F eq (2.54)
where x ′ is a spatial position for the integration. We can now make the canonical change of variable (x, p) → (α, J) which gives

j (x, t) = eZ e dτ ⋆ dα 1 δ(x ′ -x) v ,eq δF - eZ e m e δA F eq (2.55)
where we have separated the integration variables into gyrophase dα 1 ≡ dϕ c and dτ ⋆ = dα 2 dα 3 dJ 1 dJ 2 dJ 3 . We make use of Parseval's identity for Fourier series,

dα 1 2π δ(x ′ -x)δF = n 1 [δ(x ′ -x)] n 1 δF † n 1 (2.56)
where the † designates the complex conjugate of the perturbed quantity. The n 1 mode of the Dirac function can be calculated using the definition of the Dirac distribution

[δ(x ′ -x)] n 1 = dα 1 2π δ(x ′ -x)e -in 1 α 1 = dα 1 2π d 3 k (2π) 3 e ik•x ′ -ik•x e -in 1 α 1 (2.57)
Therefore, the parallel current density in Eq. (2.55) can be expressed using the Bessel function of the first kind, defined as

dα 1 /(2π)e i(k ⊥ ρc cos α 1 -n 1 α 1 ) = (-i) n 1 J n 1 (k ⊥ ρ c ) j (x, t) = eZ e dτ ⋆ n 1 (-i) n 1 d 3 k (2π) 3 e ik•(x G -x) J n 1 (k ⊥ ρ c ) v ,eq δF † n 1 - eZ e m e δA † n 1 F eq (2.
58) Here the perturbed quantities are functions of (x G , H eq , P ϕ ). Note that through this Hamiltonian formalism one can obtain the parallel current associated to the particle motion, in terms of particle coordinates x, using physical quantities describing the guiding-center motion, which depend on guiding-center coordinates

COUPLING VLASOV AND MAXWELL'S EQUATIONS

x G . Similarly, we calculate the n 1 -mode of each component of δF and we plug it into the expression of the parallel current in Eq. (2.58). Anticipating the fact that we will later use the deeply passing particle approximations, we pass to an integration in guiding center coordinates dτ ⋆ = dα 2 dα 3 dJ 1 dJ 2 dJ 3 → J dx G dv ,eq v ⊥ dv ⊥ , where J is the Jacobian of the transformation. We now make use again of the ordering Ω 1 ≫ ω, Ω 2 , Ω 3 , which implies that the third term in the right-hand side of Eq. (??) is negligible except for n 1 = 0. Therefore, we keep only the n 1 = 0 terms, which strictly speaking is the exactly the same as performing a gyro-average. Using the Bessel function property +∞ n=-∞ J 2 n (x) = 1 [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF] and the inverse Fourier transform

d 3 k (2π) 3 δ Ĥ † (k, t) e -ik•x = δH(x, t) (2.59)
we obtain the following general expression for the parallel current in the gyrokinetic approach using the magnetic limit, i.e. neglecting the electrostatic potential j (x, t) = -e 2 J dv ,eq dµ v 2

,eq ∂ Heq F eq δA (x, t)

+ v 2 ,eq 1 B ∂F eq ∂µ 1 -J 2 0 δA (x, t) + n 2 ,n 3 v 2 ,eq ω∂ Heq F eq + n 3 ∂ J 3 F eq ω -n 2 Ω 2 -n 3 Ω 3 J 2 0 δA n 2 ,n 3 ,ω (x) ×e i(n 2 α 2 +n 3 α 3 -ωt) + 1 m e F eq δA (x, t) (2.60)
where J 0 ≡ J 0 (k ⊥ ρ c ), J 2 0 represents the gyro-average operator applied twice and the squared electron charge Z 2 e = 1 has been removed for the sake of clarity. We choose the equilibrium distribution function

F eq = n eq
(2πT eq /m e ) 3/2 e -Heq Teq

(2.61)

where n eq and T eq are the equilibrium density and temperature, respectively. After integrating the first and last terms of the right-hand side in expression (2.60), they cancel out. For electrons, one can take the small orbit limit and approximate the Bessel function as

J 2 0 (k ⊥ ρ c ) ≈ 1 for k ⊥ ρ c ≪ 1.
Therefore, the second term with the ∂F eq /∂µ vanishes. Finally, j (x, t) is written as

j (x, t) = e 2 T eq n 2 ,n 3 v 2 ,eq ω -n 3 T eq ∂ J 3 log F eq ω -n 2 Ω 2 -n 3 Ω 3 δA n 2 ,n 3 ,ω (x) e i(n 2 α 2 +n 3 α 3 -ωt) (2.62)
where the notation • • • has been used for simplicity to represent an average over gyro-centre equilibrium velocity space weighted by the equilibrium distribution function

• • • = J dv ,eq dµ • • • F eq (2.63)

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

If we restrict our analysis to deeply passing electrons, we can write α 2 = θ and α 3 = ϕ, where θ and ϕ are the poloidal and toroidal angles, respectively. Consequently n 2 = m and n 3 = n, with m and n the poloidal and toroidal mode numbers, respectively. The frequencies Ω 2 and Ω 3 are calculated in detail in Appendix ??. For deeply passing particles, these frequencies can be approximated as

Ω 2 ≈ v ,eq qR 0 (2.64
)

Ω 3 ≈ 2q r m e v 2 ,eq + µB 0 eZ e B 0 R 0 + v ,eq R 0 (2.65)
where q is the safety factor, which represents the helicity of the magnetic field lines, and r is the minor radius. B 0 and R 0 are respectively, the modulus of the magnetic field and the major radius, measured at the magnetic axis. The perturbed parallel current will be written in normalized units. For this purpose we normalize the velocities to the thermal velocity of electrons v th , the distances to R 0 , the frequencies to the transit frequency ω t = v th /R 0 , the equilibrium distribution function to n 0 /v 3 th , with n 0 some normalizing density, the parallel vector potential to B 0 R 0 ρ 2 ⋆ and the temperature to a normalizing temperature defined as

T 0 = m e v 2
th /2. The parallel current can then be rewritten as

j (x, t) = ρ ⋆ en 0 v th 2 Teq m,n,ω v2 ω -ω⋆g ω -k v -ωD δ Â m,n,ω (x) e i(mθ+nϕ-ωt) (2.66)
with • representing normalized quantities. The parallel wave vector k is given by k = (m/q -n) /R 0 and the magnetic drift frequency ω D is given by ω D = 2qn/r m e v 2 ,eq + µB 0 / (eZ e B 0 R 0 ). Notice that k vanishes on the rational surface defined as q = m/n. In the remainder of this paper, the eq subscript for the parallel velocity will be dropped for the sake of simplicity. Note that ρ ⋆ gives the typical ordering between the equilibrium and the perturbed distribution function and en 0 v th is the normalization for the current. In this expression we have introduced the generalized diamagnetic frequency ω ⋆g = nT eq ∂ J 3 log F eq , which includes all the spatial dependence of the equilibrium. For thermal particles, at lowest order in ρ ⋆ , one can write J 3 ≈ -eZ e ψ. Therefore, the derivative with respect to J 3 can be reduced to a derivative with respect to the radial position.

The quasi-neutrality condition

In a plasma, the charge density beyond the Debye sphere vanishes. This leads to the so-called quasi-neutrality condition, which replaces the Maxwell-Poisson equation

s Z s δn s = 0 (2.67)
where the sum is performed over the different species. In this manuscript, when solving analytically the quasi-neutrality condition, we will assume that electrons response adiabatically to the potential perturbations. In other words, their perturbed density is proportional to the perturbed potential. In addition, we will assume that the electron equilibrium distribution function is a canonical Maxwellian. Therefore, the adiabatic response of electrons leads to the following perturbed density

δn e = eδφ T e d 3 p e F eq,e = eδφ T e n e (2.68)
where n e and T e are the equilibrium electron density and temperature, respectively. The perturbed potential is obtained from the definition of the equilibrium by the relation δφ = φ -(2π) -3 φd 3 α. We can consider that the deviation of the electron orbits from the magnetic surfaces is negligible. This means that J 0 • = 1, α 2 = θ and α 3 = ϕ. Therefore the average over the angles is trivially replaced by the average over θ and ϕ. For the case of adiabatic electrons, the quasi-neutrality condition therefore reads

e T e (φ -φ ) = 1 n e s Z s d 3 p s δF s (2.69)
where the sum is performed over the ion species and the brackets • represent the flux-surface average

• = 1 4π 2 •dθdϕ (2.70)
The quasi-neutrality condition 2.69 can be written for each mode (m, n) by projecting onto the corresponding Fourier mode e T e (φ mn -

φ 00 ) = 1 4π 2 1 n e s Z s d 3 p s δF s e -i(mθ+nϕ) dθdϕ (2.71)
The integral with respect to p s can be written as an integral over the whole phase-space as follows

d 3 p s = δ (x -x ′ ) d 3 p s d 3 x ′ = δ (x -x ′ ) dα 1 dτ ⋆ (2.72)
where dτ ⋆ = d 3 Jdα 2 dα 3 . We can therefore write

d 3 p s δF s = dτ ⋆ n 1 [δ (x -x ′ )] n 1 δF † sn 1 (2.73)
where † represents the complex conjugate. The n 1 Fourier mode of the Dirac distribution is calculated as

[δ (x -x ′ )] n 1 = 1 2π d 3 kdα 1 e ik•(x-x ′ ) e -in 1 α 1 = d 3 ke ik•(x-x G ) i n 1 J n 1 (k ⊥ ρ s ) (2.74)

CHAPTER 2. GYRO-KINETIC DESCRIPTION OF INSTABILITIES

where we have decomposed the position x ′ into the guiding-centre position and the Larmor radius, i.e. x ′ = x G +ρ s , and used the definition of the Bessel function.

We can therefore write

d 3 p s δF s = dτ ⋆ d 3 ke ik•(x-x G ) n 1 i n 1 J n 1 (k ⊥ ρ s ) δF † sn 1 (2.75)
Dividing the linear response into adiabatic and non adiabatic parts and performing similar calculations to those presented for Ampère's law, one finds the quasi-neutrality equation e T e (φ mn -φ 00 ) -

1 n e s Z s ∇ ⊥ • n e B 0 ω c ∇ ⊥ φ = 1 n e s Z s δn Gs (2.76)
where ∇ ⊥ = ∂ r , 1 r ∂ θ is the perpendicular gradient and δn Gs is the perturbed density of guiding-centres, expressed as

δn Gs = J dv dµJ 0 • Fs -Fs,eq
(2.77)

THE LINEAR ORIGIN OF THE EGAM

The linear origin of the EGAM

EGAMs were predicted theoretically [START_REF] Fu | Energetic-particle-induced geodesic acoustic mode[END_REF]50,[START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF]42,44,1], observed experimentally [39,38,[START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF][START_REF] Hl Berk | Fast excitation of egam by nbi[END_REF] and analysed in numerical simulations [42,43,46,45,47]. As the standard GAMs [START_REF] Winsor | Geodesic acoustic waves in hydromagnetic systems[END_REF], they are characterised by axisymmetric components of electrostatic potential, density and pressure. The electrostatic potential is dominated by a zonal structure, whereas density and pressure exhibit a poloidally up-down asymmetric structure, which is also visible in the electrostatic potential. The EGAMs are axisymmetric modes and are therefore driven mainly by the positiveness of the slope of the distribution function of ions at a certain position in velocity space, namely at v = qRω EGAM [42]. This can be understood by looking at Eq. 2.26b. The denominator provides the resonance condition, whereas the numerator provides the drive of the instability. For axi-symmetric modes n 3 = 0 and therefore the only drive comes from the gradient of the equilibrium distribution function in velocity space. This is the reason why EGAMs must be described kinetically. The existence of these modes is based on the geodesic curvature of the magnetic field lines, as occurred with standard GAMs [START_REF] Winsor | Geodesic acoustic waves in hydromagnetic systems[END_REF]. When coupling Vlasov and Poisson equations, one can derive a dispersion relation which take the general form

D (Ω) = 0 (3.1)
where Ω = ℜ (Ω) + ℑ (Ω) is the complex frequency. We have made the choice in this dissertation to normalise the frequency to v th /R 0 , where v th is the thermal velocity and R 0 is the major radius measured at the magnetic axis. Such dispersion relation has been derived in different works using different approaches [START_REF] Fu | Energetic-particle-induced geodesic acoustic mode[END_REF]50,[START_REF] Qu | Energetic geodesic acoustic modes associated with two-stream-like instabilities in tokamak plasmas[END_REF]. In essence, one has to define an equilibrium distribution function for both thermal and energetic particles. Solving this dispersion relation can be a difficult task, but simplifications can be introduced in order to allow the analytical or numerical calculations of the solution. For instance, we can consider that the energetic particle population is modelled by a Maxwellian distribution function shifted in parallel velocity. Whereas such assumption is valid for NBI dominated plasma with tangential injection, it might revealed incomplete in the case of perpendiculat NBI injection or ICRH discharges. Under this assumption of Maxwellian distribution function, an analytic expression for the dispersion relation can be derived as in [44,1]. Its solution is possible only numerically though. Figure 3.1 shows the position of the zeros of D (Ω) for q = 2 when the fraction of energetic particles, defined as the ratio of energetic particle density n EP to total particle density n eq , is increased little by little. It can be observed that the GAM, which is located initially at a frequency Ω = 2.6 and labeled as Mode 1 in figure 3.1a, is first damped when introducing a few energetic particles, as observed in figure 3.1c. At the same time, the mode labelled as Mode 2 in figure 3.1a exhibits an increase of frequency and an increase of growth rate (which remains always negative). The evolution of Mode 2 is accompanied by a decrease of the GAM frequency (Mode 1). This Mode 1 will eventually reach the real axis, as observed in figure 3.1e, and be excited with a further increase of energetic particle concentration as observed in figure 3.1f . Therefore, the mode EGAM is born from the standard GAM whose frequency has been decreased with respect to the situation where only thermal particles exist. This picture changes when increasing the safety factor, for instance q = 3. In that case we have observed that the Mode 1, which is initially the GAM, is not excited when introducing energetic particles and its frequency is increased. In that case, it is the Mode 2 which is excited, travelling on the complex plane towards decreasing frequencies, positive growth rates and therefore becoming the so-called EGAM. Note that this EGAM in the absence of energetic particles does exist, but is highly damped. It is contained in the branch of highly damped poles. This evolution is illustrated in figure 3.2.

Therefore, the dispersion relation facilitates the understanding of the evolution of the modes in the complex plane and is especially useful to understand what the origin of the EGAM is. We see that modifying the value of the safety factor leads to the excitation of the EGAM from a different zero of D (Ω). This is easily understood since the value of the safety factor plays an essential role in the resonance condition derived from the denominator of Eq. 2.20, when assuming passing particles, i.e. ω -(qR 0 ) -1 v = 0. In conclusion, the mode called EGAM is not really a new mode, but an already existing mode excited in the presence of energetic particles. Depending on the values of q, it is born either from the initial GAM or from an initially highly damped mode.

LINEAR GYROKINETIC SIMULATIONS OF EGAMS

Linear gyrokinetic simulations of EGAMs: comparison with analytic theory and quantification of the wave-particle energy exchange

The solution of the dispersion relation can be compared to the results of gyrokinetic simulations where the Vlasov equation is coupled to the quasi-neutrality condition in the electrostatic limit with adiabatic electrons. A scan on the fraction, mass and charge of energetic particles using deuterium as thermal population can reveal the good agreement between theory [44,1] and simulations [START_REF] Zarzoso | Particle transport due to energetic-particle-driven geodesic acoustic modes[END_REF] using for example the GYSELA code [82]. This is shown in figure 3.3, where lines represent theoretical predictions from the dispersion relation derived in [1] and symbols represent values measured from GYSELA simulations. For the simulations reported in that figure, all the radial profiles are flat, the safety factor is set to q = 1.8, the temperatures are chosen so that T th = T EP = T e and the energetic particle mean velocity is set to v EP = 4v th,EP , where v th,EP is the thermal velocity of energetic particles defined as v th,EP = T EP /m EP .

More in-depth analysis can be done by examining the energy which is directly transferred from particles to the mode. In that respect, a multi-species code like represents an optimal tool to determine the energy which is transferred from each (lines) and as measured from GYSELA simulations (symbols). The scan is performed using the same species for the thermal population and modifying the mass and the charge for the energetic particle population.

species to the mode. The energy transfer is quantified by means of the expression

dE kin,s dt = d 3 x J s • E (3.2)
which represents the increase or decrease of kinetic energy of species s interacting with the electric field E. In that expression, J s is the current carried out by species s. The equation 3.2 is calculated in a gyro-kinetic code as follows

dE kin,s dt = -eZ s dx GC dt • ∇ (J 0 • φ) δF s dτ (3.3)
where x GC is the position of the guiding-centre, J 0 • is the gyro-average operator, δF s is the perturbed distribution function of species s and dτ = B * d 3 x GC dv dµdα 1 is the phase-space volume element. The perturbed distribution function is defined as δF s = F s -F s,eq , with F s,eq the equilibrium distribution function defined by equation 2.18. It is straightforward to see that eZ s ẋGC δF s is the current density due to species s, namely J s , and -∇ (J 0 • φ) is the electric field as seen by the guiding-centre. In that sense, equation 3.3 represents indeed the wave-particle exchange of energy integrated over the whole phase-space. The expression 3.3, in the absence of collisions between species, must equal the decrease or decrease of energy of the electric field due to each species E f,s , i.e. [83]

dE kin,s dt = - dE f,s dt (3.4)
where the field energy is defined to the lowest order in perturbations as [84,85]

E f,s (t) = eZ s 2 J 0 • [n s (x GC , t) -n 0,s (x GC )] φ (x GC , t) d 3 x GC (3.5)
The total energy of the field can then be obtained simply by summation over all the species, i.e.

E f = s E f,s (3.6)
Neglecting second order terms due to the E × B drift, the time derivative of the kinetic energy can be decomposed into different contributions of the guidingcentre velocity

dE kin,s dt = eZ s δF s v b + v µ∇B + v curv • ∇ (J 0 • φ) dτ (3.7)
where

v µ∇B = B eZ s B 2 × µ s ∇B (3.8)
is the µ∇B drift, due to the gradient of the magnetic field, and

v curv = m s v 2 eZ s B b (×b • ∇b) (3.9)
is the curvature drift.

Using this decomposition, we can analyse separately which terms in the J • E are significant for the growth of the mode, especially if one wants to determine the nature of the excited mode. For instance, is the mode driven mainly by trapped or passing particles? In the mode driven by the dynamics in the parallel direction or by the curvature of the magnetic field? To shed light on these questions, the different components of the time derivative of the energetic particle kinetic energy as written in Eq. 3.7 are plotted in figure 3.4. The thin solid line corresponds to the total time derivative. The thin dashed line represents the contribution of the curvature term. The thick dashed line corresponds to the contribution of the µ∇B term and the thick solid line gives the time evolution of the parallel velocity contribution. We can observe the double nature of the EGAM. On one hand, the main contribution comes from the curvature, which is in agreement with the fact that the curvature is responsible for the existence of the geodesic acoustic mode. On the other hand, the parallel velocity term becomes important when the curvature term exhibits a local minimum, which is in agreement with the fact that mainly passing particles excite the mode through the resonance in the parallel direction. The existence of an energy exchange via the curvature drift is in agreement with a recently reported EGAM branch due to the magnetic drift [54].

Using the energy exchange diagnostic, it is also possible to analyse the energy exchanged between thermal/energetic particles and the mode in order to provide further understanding on how the mode is excited. In figure 3 the time evolution of the energy transferred from the energetic population to the mode (solid red curve), from the thermal particles to the mode (dashed blue line) and the contribution of all particles (dotted-dashed black line), which is the sum of energetic and thermal particle contributions. The left and right panels represent, respectively, the time trace at and beyond the threshold of the instability. In particular, at the threshold of the instability, the energy transferred to the mode remains constant and oscillates around zero, as expected. We can observe that, beyond the threshold, the energetic particles give, on average, energy to the mode, whereas the thermal particles receive energy from the mode. This is of course in agreement with the standard picture of direct and inverse Landau damping, where the exciting population transfers energy to the mode and the mode in its turn transfers energy to the stabilizing population. However, it can be observed that the transferred energies oscillate, exhibiting maxima and minima, which is in agreement with the fact that both the current and the electric field are oscillating quantities at a frequency which is twice the EGAM frequency. This means in particular that thermal particles are sometimes more stabilizing and sometimes less stabilizing. Similarly, energetic particles are sometimes more destabilizing and sometimes less destabilizing. A magnified view for a given time window is provided, showing that thermal particles can actually transfer energy to the mode when the energy transferred from energetic particles to the mode exhibits a local minimum value. This occurs during the growth of the total energy transferred from particles (both thermal and energetic) to the mode. Therefore, though the global energy exchange between particles and mode could be understood as a transfer of energy from the exciting particles to the mode and from the mode to the stabilizing particles, in-depth analysis of the energy channels could pro-

NONLINEAR SATURATION OF EGAMS

vide additional information to understand more in detail the contribution of each species during the growth phase of the mode. For this purpose, we analyse in the next section the interaction of particles and the mode in the nonlinear phase of the EGAM, where the formation of islands in phase-space reveals the existence of a wave-particle interaction.

Nonlinear saturation of EGAMs

As occurs with any instability, the EGAM exhibits a linear growth, where the amplitude of the mode increases exponentially in time, and nonlinear phase, where the amplitude of the mode saturates and exhibits eventually oscillations related to the trapping and detrapping of particles inside the potential well. A standard picture of the linear growth and nonlinear saturation of the EGAM is given in figure 5.8a, obtained from a gyro-kinetic simulation using the multi-species version of GYSELA and representing the time trace of the electrostatic potential. Figure 3.6b gives the Fourier transform in time, which evidences the presence of a dominant frequency.

The nonlinear saturation of EGAMs was analysed for the first time in [42], where the flattening of the equilibrium distribution function was evidenced by conveniently averaging the full distribution function over flux surfaces. The flattening of the distribution function can be expected, since the saturation occurs with the formation of structures that tend to cancel the drive of the instability. If the instability is driven by the gradient of the equilibrium distribution function in velocity space, the nonlinear saturation is expected to occur with a flattening of that equilibrium distribution function. Moreover, this flattening indicates the position in phase-space of particles resonating with the mode. Therefore, using different distribution functions for thermal and energetic particles allows one to determine to what extent the excitation of EGAMs by energetic particles affects the thermal population. In particular, the existence of higher order resonances can lead to a significant interaction between energetic and thermal particles through the excitation and damping of EGAMs. For instance, let us assume that energetic particles excite the main component m = 1 of EGAM by resonating around the velocity v = v res leading to the EGAM frequency ω EGAM = v res / (qR). Then, the appearance of higher order resonances m > 1 leads to resonant velocities v (m) res = qRω EGAM /m with smaller absolute values when m increases and therefore approaching the thermal velocity. In addition, higher order resonances in frequency may occur, leading to frequencies lω EGAM , for l > 1. Therefore, for a given poloidal mode number m, this will result in more general expressions for resonant velocities v (m,l) res = qRlω EGAM /m. For the main poloidal harmonic m = 1, we have v (m=1,l) res = qRlω EGAM , representing an interaction between the EGAM and suprathermal particles. The modification of the equilibrium distribution function is quantified by the distance between the equilibrium distribution function and the initial equilibrium distribution function, averaged over time:

δF eq v = 1 T T 0 dt F v , µ = 0, θ, ϕ, t FS -F v , µ = 0, θ, ϕ, t = 0 (3.10)
where • • • FS represents the flux-surface average and T is a sufficiently large time. It is to be noted that if a flattening occurs during the nonlinear saturation, this is indicated by a δF eq v exhibiting an odd parity around the resonant velocity. An example of δF eq v is plotted for both thermal and energetic passing particles in figure 3.7 and represented by thick lines (solid blue line for thermal particles on the left axis and dashed red line for energetic particles on is given by dashed black lines. The higher order resonances m = 2 and m = 3 for l = 1 are given by dashed red and blue lines, respectively. The higher order resonance in l for m = 1 is represented by the dashed-dotted black line. Note that only l = 2 is visible in the figure. For l > 2 and m = 1, the resonant velocity is outside the parallel velocity grid used for the simulation. Finally the higher order resonance in both m and l is shown by the solid magenta line for (l, m) = (2, 3). For the sake of clarity all the discussed resonances are labelled in the figure and a magnified view of the (l, m) = (2, 3) resonance is provided on the top left panel. In figure 3.7 it can be observed that the distribution function of both thermal and energetic particles is perturbed around the various resonances discussed above, showing that although energetic particles are mainly responsible for the excitation of the mode, thermal particles can also exchange energy. This provides an efficient mechanism for exchange of energy between energetic particles and thermal plasma in the absence of collisions through the excitation of EGAMs.

The regions where the distribution function is perturbed have a width that is linked to the width of the island in phase space given by the amplitude of the mode. Since the EGAM is given by φ ∼ sin θ, observing the EGAM island requires to keep the dependence on the poloidal angle, hence not to perform any flux surface average. Note that avoiding the flux-surface average will not introduce any toroidal dependence, since the EGAM is an axi-symmetric mode. Therefore, the analysis can be performed at a given toroidal angle. ing particles, at ω c,s0 t = 4 • 10 4 , which corresponds to a time after nonlinear saturation. The EGAM island at the main resonance (l, m) = (1, 1) is clearly visible for the energetic particle population, whereas the thermal population is mostly affected by the higher order resonance (l, m) = (1, 2) EGAM island. Figure 3.9 shows the same quantities at the same time but when µ is increased, so that the presence of the trapping cone can be observed, mainly in the thermal population, which is closer to the trapping cone than energetic population. Moreover further higher order resonances are evidenced for thermal particles, getting closer to the trapping cone and being therefore distorted. The EGAM island that is formed in phase-space can interact with the island related to the magnetically trapped particles. In the next chapter of this dissertation we will see how both island can interact, leading to a chaotic region which results in losses of particles.

Finally, it is to be noted that a spectrogram of the time trace shown in figure 5.8a would reveal that the mode does not exhibit any chirping. The chirping of EGAM was reported in previous publications, both experimentally [START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF] and numerically [46]. Very recently, in global gyrokinetic simulations using the ORB5 code [START_REF] Jolliet | A global collisionless pic code in magnetic coordinates[END_REF], a chirping of EGAM was also observed [START_REF] Biancalani | Saturation of energetic-particle-driven geodesic acoustic modes due to wave-particle nonlinearity[END_REF], when only wave-particle nonlinearities are kept in the gyro-kinetic equations. It was later shown in [START_REF] Biancalani | Nonlinear gyrokinetic investigation of energetic-particle-driven geodesic acoustic modes[END_REF] that adding wave-wave nonlinearities in that case may result in the strong reduction of the chirping. In our simulations, both wave-particle and wave-wave nonlinearities are kept. This might explain why no chirping is observed for our particular case. In addition, as reported in [START_REF] Vann | Fully nonlinear phenomenology of the berk-breizman augmentation of the vlasov-maxwell system[END_REF][START_REF] Lesur | Nonlinear categorization of the energetic-beamdriven instability with drag and diffusion[END_REF], a steady-state solution is not completely excluded, which indicates that a set of parameters exist for which GYSELA should not predict any chirping. Nevertheless, further parametric studies and quantitative comparisons between GYSELA and other gyro-kinetic codes such as ORB5 in nonlinear regime are required in order to conveniently catego- rize the behaviour of the EGAM depending on the different parameters that characterize the gyro-kinetic simulations, which is beyond the scope of the present dissertation.

BRIEF INSIGHT INTO THE CONTEXT OF EGAM-INDUCED TRANSPORT

Brief insight into the context of EGAM-induced transport

Since EGAMs are axisymmetric modes, the toroidal canonical momentum, P ϕ , is conserved, as explained in chapter 2. Having a look at Eq. 2.11, it is clear that P ϕ can be approximated (up to a multiplying factor) by the poloidal flux for low toroidal velocities, which is the case of thermal particles. Therefore, for thermal particles, the conservation of P ϕ is equivalent to the conservation of the poloidal flux, which is a function of the radial position. This means that axisymmetric modes (and in particular EGAMs) are not expected to induce a significant radial transport of thermal particles. However, for energetic particles, the contribution of the toroidal velocity can no longer be neglected in the toroidal canonical momentum. In this case, any modification of the toroidal velocity results in a substantial variation of the radial position, in virtue of the conservation of P ϕ . This was first observed in DIII-D discharges [START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF], where 10% -15% drops in the neutron emission occurred, indicative of strong beam ion redistribution and even losses. It was speculated that EGAMs might modify the pitch angle of counterpassing particles, resulting in unconfined trapped particles. This was later corroborated in [41] by comparison between a DIII-D discharge, where fast-ion loss detector (FILD) [START_REF] Garcia-Munoz | ASDEX Upgrade Team[END_REF] data were available, and numerical results from full orbit simulation code SPIRAL [92]. Also, a recent work provides analytical calculation of the effect of the EGAM on well passing particles that are trapped in the EGAM island and subsequent impact on the neutron emission [93].

Numerical observation of the EGAM-induced losses of particles

Within the electrostatic limit, the evolution of the gyro-centre coordinates of species s is given by the equations [82]

dx i dt = v b * • ∇x i + v E • ∇x i + v D • ∇x i (4.1a) m s dv dt = -µb * • ∇B -eZ s b * • ∇J 0 φ + m s v B v E • ∇B (4.1b)
where x i is the ith contravariant component of the coordinate x (x 1 ≡ r in the radial direction, x 2 ≡ θ in the poloidal direction and x 3 ≡ ϕ in the toroidal direction), v the parallel component of the velocity along the magnetic field lines, v E is the E × B drift, v D is the magnetic drift, µ is the magnetic moment, which is an invariant within the present model, m s is the mass of particles, e is the elementary charge, Z s is the atomic number, φ is the electrostatic potential, B is the magnitude of the magnetic field, J 0 is the gyro-average operator and b * is defined as

b * = B B * + m s v eZ s B * B ∇ × B (4.2) CHAPTER 4. TRANSPORT DUE TO EGAMS with B * = B + m s eZ s v b • ∇ × b (4.3)
the volume element in guiding-centre velocity space, where we have introduced the magnetic field B and the unit vector along the magnetic field lines b. The basics of tensor calculus and toroidal geometry used in this dissertation are summarised in appendix A. In particular, we use the simple model of concentric circular magnetic flux surfaces, with the magnetic field defined by equation (A.18).

The ith contravariant components of the drifts are calculated as follows

v E • ∇x i = 1 B * b • ∇J 0 φ × ∇x i (4.4) v D • ∇x i = m s v 2 + µB eZ s B * B b • ∇B × ∇x i (4.5)
where the triple product is calculated as

b • (∇F × ∇G) = √ gǫ ijk b k ∂ i F ∂ j G (4.6)
with g = det {g ij } the determinant of the metric tensor and ǫ ijk the contravariant components of the Levi-Civita tensor. An example of the numerical integration of these equations is given in figure 2.2, which was obtained using a 4 th order Runge-Kutta method in the absence of any perturbation and with an axisymmetric equilibrium magnetic field. These equations can also be solved in the presence of an electrostatic potential obtained from a gyro-kinetic code, for instance from GYSELA, ORB55 or any other code. In that case, the equations of motions are integrated for passive tracers, i.e. for particles which do not have any impact on the fields. Of course, these equations must satisfy the conservation of the motion invariants, i.e. in the absence of any perturbation both P ϕ and E must be conserved and in the presence of an axi-symmetric time-dependent perturbation only P ϕ must be conserved (the total energy is conserved, but not the kinetic energy E). This is illustrated in figure 4.1. In the left panel, the surfaces E = cst and P ϕ = cst intersect each other forming a 1D sub-space which corresponds exactly with the trajectory of the particle. When an EGAM is introduced (obtained from GYSELA), it is observed that the particle, whose trajectory is highlighted by a dotted yellow line, remains on the P ϕ = cst surface, whereas it departs from the E = cst one.

The transport of passive tracers can be evaluated as a function of their initial positions for an ensemble of particles, instead of doing it for a single tracer. This is illustrated in figure 4.2. This figure has been obtained by tracking an ensemble of passive particles intialised with 0 < µ < 16 and -7 < v < 7 in the presence of an electrostatic potential obtained from GYSELA with 0.02 fraction of energetic particles. Whenever a particle intercepts the wall of the tokamak (assumed a torus with circular cross-section) a point is marked in the figure. Blue points represent prompt losses. These are losses that exist even in the absence of EGAMs due to the fact that the initialisation of particles is such that their orbits intercept the boundary and are lost from the very beginning of the simulation. Red points represent losses that exist only in the presence of EGAMs and are therefore called EGAM losses. The left panel represents losses of particles initialised at an inner radial position (r/a = 0.2) and the right panel indicates the losses of particles initialised in an outer region (r/a = 0.6). As already conjectured in [START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF] and observed numerically in [41], mainly counter-passing particles experience EGAM losses, independently of the radial position. The first interpretation proposed in [START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF] was that the transport of particles in phase space due to the EGAM might modify the orbits in the sense that passing particles might be trapped and trapped particles might be de-trapped. This is due to the conservation of the toroidal momentum P ϕ = -eZ s ψ + m s Bϕ B v , which makes the link between the radial position given by the poloidal flux ψ and the parallel velocity v . For ions (Z s > 0), conserving P ϕ implies that an acceleration of particles leads to a radial transport towards the outer region of the tokamak. The acceleration in the parallel direction makes counter-passing particles approach the trapping-cone and become trapped. Their orbits can subsequently intercept the wall of the tokamak. Note that for co-passing particles that approach the trapping-cone, their parallel velocity is decreased. Therefore the radial position is also reduced owing to the conservation of P ϕ . They are therefore transported to a region where the trapping-cone is narrower. As a result co-passing particles are unlikely trapped by the magnetic field. This explains why neither the trapping/de-trapping mechanism due to EGAM nor the losses shown in figure 4.2 are symmetric with respect to v .

Chaotic transport in the presence of a non-chaotic electrostatic potential

The previous section proves that, in agreement with experimental observations [START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF] and full-orbit simulations [41], a transport due to EGAMs is possible. However, further understanding of the underlying mechanism is needed. For this purpose, one can analyse the trajectories of individual particles whose losses are identified as EGAM-induced losses. The trajectories of three of thoses particles are plotted in the top panel of figure 4.3. The left panel shows the projection onto the poloidal cross-section of one of the particle trajectories represented by solid blue lines (note that for the sake of clarity, only one trajectory can be plotted on the cross-section). It is actually crearly observed the transition from the counter-passing trajectory to the trapped and lost orbit, as explained in the previous section. The middle panel shows the time trace of the parallel velocity for the three particles (solid blue line, dashed black line and dotted red line), and the radial position of the particles as a function of time is given in the right panel. For reference, the bottom panel shows the same quantities for the unperturbed trajectory of the blue particle, i.e. the trajectory of the blue particle without EGAM.

For both perturbed and unperturbed trajectories, the boundary domain in real space, corresponding to a circular flux-surface, is given in the left panel and represented by a dotted black line. The three particles are initialised with slightly different parallel velocities in the neighborhood of v (t = 0) = -5.5v th , which is close to the main EGAM resonance, and in figure 3.7 it is located on the left side of the resonance indicated by the black dashed line. The three particles are also initialised with slightly different magnetic moments in the neighborhood of µ = 13T 0 /B 0 . When the EGAM is excited, the particle parallel velocities start exhibiting low frequency oscillations during a synchronisation process [94] with the m = 1 component of the mode. This phase is represented by the magnified view only for the blue particle. The fact that the trajectories of the three particles, while initialised very close to each other in phase-space, exhibit such different behaviour is an indication of a possible chaotic transport. More detailed analyses can be done to study the dynamics of particles in the vicinity of the main EGAM resonance. Moreover, for the sake of interest for the fusion community, it is essential to analyse what happens to energetic particles initialised in the inner region of the tokamak. Indeed, this is crucial for alpha particles born from nuclear fusion reactions as well as energetic ions generated by external heating are mainly localised in the core of the plasma. Since the main dynamics of the EGAM occurs in the (θ, v ) sub-space, higher resolutions in v than in µ are preferred. Also, in order to improve the statistics of the analyses of the particles close to the main EGAM resonance, the total number of integrated trajectories must be as high as possible. The domain selected for this detailed analysis is indicated in the left-hand side top panel of the trapping cone. This indicates that the interaction between the EGAM and the trapping cone might indeed be a key point for the dynamics of the physical system. Nevertheless, no clear pattern can be identified within the cloud of points, suggesting that a deterministic process is unlikely to result in the observed losses. The number of particles that are lost is plotted as a function of time in the middle panel of figure 4.4, represented by a black line. The time-trace of the flux-surface averaged electrostatic potential is over-plotted and represented by a magenta line. This figure indicates that there is a maximum in the number of losses right after the EGAM saturates and the losses occur during the whole nonlinear phase of the EGAM, suggesting that the losses are due to nonlinear effects. Finally, the time-Fourier transform of the losses of particles has been performed and plotted in the right panel of figure 4.4. The amplitude of the Fourier modes exhibit dominant components at the EGAM frequency and its harmonics, in agreement with experimental results [START_REF] Nazikian | Intense geodesic acousticlike modes driven by suprathermal ions in a tokamak plasma[END_REF].

Finally, it is especially enlightening to show the distribution of particles in phase-space. This is represented in figure 4.5 at different instants, indicated by dashed vertical lines (blue, red and green, respectively) in the middle panel of figure 4.4: (1) before, (2) just after and (3) well after the nonlinear saturation of the EGAM. For the sake of clarity, the radial position of particles is indicated by colours, being blue and red the two limits representing r/a → 0 and r/a → 1, respectively. The trapping cone is represented by a dashed black line. In the left panel of figure 4.5 it is observed that some particles remain passing (in terms of EGAM island) and other particles are trapped inside the EGAM. The trapping of some of the particles allows one to observe clearly the formation of the EGAM island in phase space. When the EGAM saturates, the EGAM island strongly interacts with the trapping-cone, as observed in the middle panel of figure 4.5. This panel shows the instant where the first particles of the ensemble are lost and indicates that the losses are produced through a chaotic channel established between the counter-passing particles and the lost particles. The right panel of figure 4.5 shows the same structure later on in the simulation. This interaction occurs during the whole nonlinear phase. Note that this picture is slightly different from the general understanding à la Chirikov [95]. It is usually believed that chaos originates from the overlap of two contiguous islands, making it possible for a particle to explore regions inside the two islands. This overlap is likely to be produced in the region where the separatrices are closer, i.e. when the two O-points are aligned with each other. Nevertheless, one has to consider that the separatrix of one island is nothing else but the stable and unstable manifolds of its X-point, which merge into one single manifold when the Hamiltonian system is integrable in the sense of Poincaré. When the system is not integrable (which is our case here), the unstable and stable manifolds of the X-point do not merge into one single surface and they start being folded, leading to the homoclinic tangle. Now, if one considers the X-point of the contiguous island, the same picture occurs, resulting in the heteroclinic intersection of the stable manifold of one Xpoint (that of the trapping cone) and the unstable manifold of the other X-point (that of the EGAM) [96,[START_REF] Allan | Regular and stochastic motion[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF].

Statistical description of the particle losses using a reduced model

The previous results indicate that EGAM can induce losses of particles due to a stochastization of the passing-trapping boundary layer. The question arises about how fast energetic particles are transported and lost and about the nature of the transport. Is it diffusive, super-or sub-diffusive? Will energetic particles stay sufficiently well confined to transfer their energy to the thermal population? We will answer these questions through statistical analysis. However, in order to have meaningful statistics, we need to simulate a huge number of test particles during sufficiently long gyro-kinetic simulations. Due to computational restrictions, we avoid this approach by replacing the EGAM potential obtained from the expensive direct gyro-kinetic simulations by an analytical model containing the main physics of the EGAM. This strategy is the most favorable in terms of CPU time, since no interpolation of the field is required. In addition, it allows us to simulate trajectories on time scales comparable with the experimental measurements. The main characteristics of a zonal (n = 0) structure such as an EGAM are: (1) its frequency, (2) its spatial (radial and poloidal) structure and (3) its amplitude. A zonal structure can therefore be modeled as

φ (r, θ, t) ≈ [φ 00 (r) + φ 10 (r) sin θ] cos (ωt) (4.7)
Based on the ordering φ 10 ∼ 10 -1 φ 00 , we neglect in the following the poloidal dependence and focus only on the dominant component. Following gyro-kinetic simulations [42,43,48] we can model the radial dependence as

φ 00 (r) = φ 00 1 -tanh r -r 0 δr (4.8)
where φ 00 is the value of the potential at r = r 0 and δr controls the width of the mode. This gives a radial electric field of amplitude E r,0 = φ 00 /δr at r = r 0 and localised in a region r 0 -δ r /2 < r < r 0 + δ r /2.

Description of particle losses in terms of motion invariants

The description of particle dynamics in terms of motion invariants is especially suitable when one wants to describe how the trajectories are modified in the presence of a perturbation. This is due to the fact that a particle initialised with a set of three motion invariants in the absence of any perturbation will explore the phase-space keeping the motion invariants constant, which translates into a 1D curve in the 3D real space. We know that this trajectory will correspond to the one of any other particle initialised in such a way that at t = 0 it has the same motion invariants. The unperturbed trajectory of a particle can for instance be described by the kinetic energy E and the pitch angle Λ = µB 0 /E, providing the radial position of the injection is known. An example is illustrated in figure 4.6, where the left panel represents the projection onto the poloidal cross-section of the trajectories of two particles, one deeply counter-passing and the other barely counter-passing, both injected at the radial position r/a = 0.1. Of course, with a real heating, the particles are injected within a certain range of energies and pitch angles. If all the particles are injected roughly at the same position, the resulting trajectories will cover the area represented by the blue region in the right panel of figure 4.6. One can use then the analytic expression for the EGAM potential given by Eq. 4.7 and study the losses of all those particles injected at r/a ≈ 0.1 within a certain range of energies and pitch-angle. This is given in figure 4.7. The left panel represents the fraction of lost particles as a function of their initial kinetic energy and pitch angle. Overlaid are the iso-contours of the toroidal canonical momentum, for the sake of completeness. The right panel represents the Poincaré map of a lost particle. The blue circle indicates its unperturbed trajectory. This is a clear example of how an initially confined counter-passing particle can be lost and the area where it is deconfined is actually given by the intersection of the chaotic region (the red one) with the tokamak wall (the black circle). In the right panel it is also possible to observed two islands: the larger one corresponds to the island of the magnetically trapped particles and the smaller one corresponds to the island of the electrostatically (EGAM) trapped particles. Of course, the region where the chaotic sea intercepts the wall of the tokamak might be different for each particle depending on the initial energy and pitch-angle. Also, analysis of the characteristics of the particles when they are lost reveals especially interesting.

For example, we can determine how the poloidal angle at which the particles are lost correlates with the velocity of the particle. This can be quantified by using the cosine of the pitch angle, which is defined by λ = v /v. It is to be noted that this quantity is not a motion invariant. that sense, it is intuitive to think that the exit time of the lost particles will also depend on the initial conditions. Analysing the dependence of the exit time on the initial conditions is especially relevant, since it allows us to identify the existence of patterns or structures. This cannot be done with the diagram of lost particles represented in figure 4.2, because the quantity plotted there was binary (either the particle is lost or it remains confined). Of course, since the exit time of particles that are never lost is infinite, the best way to represent the exit time is by plotting its inverse, i.e. t -1 exit . Figure 4.8 provides the inverse of the exit time as a function of the initial parallel velocity and magnetic moment particles. The top panel indicates the inverse of the exit time for all the simulated particles. A clear pattern of structures aligned with the trapping cone is observed. Half of the diagram represents particles that are never lost. We can focus on a smaller region containing lost particles. This region is represented by the yellow rectangle. A new simulation with increased resolution in time and velocity space in order to detect the small structures was performed. The exit time is plotted in the middle panel of figure 4.8, which reveals smaller and more intricate structures that were not visible in the top panel. A new simulation focusing on an even smaller region provides smaller structures, as indicated by the bottom panel. This property where the structure does not seem to depend on the scale at which the problem is analysed is called fractal behaviour [START_REF] Benoit | The fractal geometry of nature[END_REF], and is found in many other situations, such as in biological systems [START_REF] Nonnenmacher | Fractals in biology and medicine[END_REF], in the financial market [START_REF] Walter | Les structures du hasard en économie: efficience des marchés, lois stables et processus fractals[END_REF] in magnetized plasmas [102] and in the time-space diagram of some linear cellular automata [START_REF] Stephen | Cellular automata can generate fractals[END_REF][START_REF] Stephen | Computing fractal dimensions for additive cellular automata[END_REF].

Scale invariance: quasi-fractal behaviour and anomalous diffusion

We can now focus the analysis on a more restricted region in velocity space, selecting almost mono-energetic EP injected in a localised region of the tokamak and determining the probability distribution function (PDF) of their exit time.

STATISTICAL DESCRIPTION OF THE PARTICLE LOSSES USING A REDUCED MODEL

ping cone. This cone is a well defined surface in phase-space, also called separatrix, since it represents the separation between the two classes of particles. The black lines in figure 4.10 represent the Poincaré map of the unperturbed trajectories for passing and trapped particles. The left panel represents the projection onto the poloidal cross-section, i.e. onto the (R, Z) sub-space, and the right panel represents the projection onto the (r 2 /2, θ) sub-space. The dashed blue lines with arrows in the left panel indicate the direction of the trajectories of the particles contained in each region, as considered in this chapter. We assume that counter-passing particles are injected in the inner part of the tokamak. Therefore, those particles rotate clock-wise. When they become trapped and eventually co-passing, they rotate in the anti-clock-wise direction. This occurs in the outer region of the tokamak, where particles can intercept the wall and be lost.

The red region represents the Poincaré map of particles located on the separatrix in the presence of an EGAM. It is clearly observed that the separatrix is transformed into a chaotic area connecting the inner and outer parts of the tokamak. More interestingly, it is to be noted that the separation between inner and outer regions is done strictly speaking in the radial direction. Since the radial region that the particle explores when going from one region to another is necessarily upper bounded by the minor radius of the tokamak, the statistics might be meaningless when focusing on the radial excursion of particles. However, due to the conservation of P ϕ , the radial position is intrinsically linked to the parallel velocity, which is in its turn linked to the time derivative of the poloidal angle through the Eq. 4.1a applied to i = 2, corresponding to x 2 ≡ θ. In other words, when the counter-passing particle is confined in the core of the tokamak the poloidal angle decreases, and when the particle becomes co-passing in the outer region the poloidal angle increases. Contrary to the behaviour of the radial position, the poloidal angle that the particle explores can be as large as possible. Therefore, the statistical analysis of the poloidal excursion can be easily done with the possibility to be connected to the radial excursion. In the following, we study only the statistics in the poloidal angle. Our analysis follows closely the one reported in [START_REF] Del Castillo-Negrete | Asymmetric transport and non-gaussian statistics of passive scalars in vortices in shear[END_REF] for the transport of passive scalars in vortices in the presence of a shear layer. For this purpose, we initialise an ensemble of ∼ 10 5 energetic particles with E ≈ 20E th in the chaotic region. The time traces over 10 6 cyclotronic periods of the poloidal displacement defined as

∆θ (t) = θ (t) -θ (0) (4.9)
for a subset of these particles are plotted as illustration in the left panel of figure 4.11, where a clear spreading is observed. The question arises whether this spreading results from a diffusion in the poloidal direction or not. This can be analysed with the variance of the poloidal displacement, defined as

σ 2 θ (t) = (∆θ -∆θ ) 2 (4.10)
where the time dependence of the poloidal displacement has be dropped for the sake of clarity and the brackets • • • represent an ensemble average. This vari- ance can be expressed as σ 2 θ (t) ∼ t γ . The classification of the transport is done in terms of the exponent γ as follows γ < 1:

sub-diffusive γ = 1: diffusive 1 < γ < 2: super-diffusive γ = 2: ballistic Both super-and sub-diffusion enter into the category of anomalous diffusion. The case of a diffusion is usually referred to as normal diffusion to differentiate it from the anomalous diffusion. Another classification is possible based on the coherent/incoherent behaviour of the transport [START_REF] Lesieur | Turbulence in Fluids[END_REF]: the normal diffusion is an incoherent diffusion whereas the ballistic motion is a coherent diffusion. Such classification was used in magnetized plasmas for the analysis of the transport of energetic particles in the presence of Alfvén waves [START_REF] Asdex Upgrade Team | Convective and diffusive energetic particle losses induced by shear alfvén waves in the asdex upgrade tokamak[END_REF], for instance. It is of course obvious that in between both coherent and incoherent diffusions we can observe the already defined anomalous diffusion. In the right panel of figure 4.11 we represent by open red symbols the time trace (log-log) of the variance of the poloidal displacement as measured from our simulations using the expressions 4.9 and 4.10. The solid red line represents the linear fit in scale log-log. For comparison, we provide the result that we should obtain if the transport was produced by ballistic (dotted-dashed black line) or diffusive (dotted blue line) processes. It is clear that our simulations are bounded by both processes, meaning that the spreading observed in the left panel is due to a super-diffusion, with an exponent γ = 1.64.

The existence of anomalous diffusion implies that the motion of the guidingcentres cannot be modelled using a diffusion equation, which has important consequences when trying to predict the transport of energetic particles by means of

STATISTICAL DESCRIPTION OF THE PARTICLE LOSSES USING A REDUCED MODEL

(a) reduced fluid models. As explained in [START_REF] Del Castillo-Negrete | Asymmetric transport and non-gaussian statistics of passive scalars in vortices in shear[END_REF], the anomalous diffusion is understood as follows [START_REF] Lesieur | Turbulence in Fluids[END_REF]. Let us consider the Lagrangian velocity

d∆θ dt = v θ (t) (4.11)
and the Lagrangian diffusion coefficient

K (t) = 1 2 d dt (∆θ -∆θ ) 2 (4.12)
which can simply be expressed as

K (t) = t 0 v θ (τ ) v θ (0) dτ = t 0 C (τ ) dτ .
Therefore the time derivative of the variance is expressed in terms of the integral of the Lagrangian velocity auto-correlation function as follows

dσ 2 θ dt = 2 t 0 C (τ ) dτ (4.13)
Therefore, the scaling of the variance depends on how the Lagrangian velocity auto-correlation function decays in time. If the auto-correlation function decays fast enough in time, the integral 4.13 exists in the limit t → ∞, meaning that σ 2 θ ∼ t and defining the diffusion coefficient K. If the auto-correlation function exhibits an algebraic decay (for γ = 1) as C ∼ t γ-2 , then σ 2 θ ∼ t γ . The anomalous diffusion is therefore related to the slow decay of the auto-correlation function. This behaviour is understood in terms of the physics at play. Indeed, an energetic counter-passing particle is injected in the inner region of the tokamak and will remain rotating in the clock-wise direction unless something (the chaotic separatrix) makes it change its radial position until it becomes magnetically trapped. Once it is trapped, the poloidal displacement vanishes on average (there is no transport). The particle will remain trapped unless the chaotic separatrix makes it change again its radial position. It will become either co-passing, evolving as if the particles was flying towards positive poloidal angles, or counter-passing, evolving as if the particle was flying towards negative poloidal angles. A superdiffusion can therefore be understood as a compromise between trapping periods and short and rare events called that we call flights, which tend to de-trapped the particles. Basically, it is caused by the presence of coherent structures (magnetically trapped and passing regions) which make particles spend an anomalous amount of time moving slowly (trapping region) or fast (passing region) [START_REF] Del Castillo-Negrete | Asymmetric transport and non-gaussian statistics of passive scalars in vortices in shear[END_REF], the bridge between both being ensured by the chaotic separatrix. The existence of the flights is clearly visible in the left panel of figure 4.12 we plot the poloidal displacement of two tracers during the first 10 6 cyclotron periods. It can be observed that sometimes the particle is magnetically trapped and therefore the poloidal displacement does not evolve on average. Sometimes, there are either positive or negative flights which de-trappe the particles. As a comparison, we give in the right panel of the figure a time trace assuming an asymmetric standard random walk, in the absence of any flights.

Coming back to figure 4.10, when an EP becomes trapped it is lost if the wall of the tokamak intercepts the chaotic region, which is the case here since a = 150ρ th . The exit time of a counter-passing particle is related to the time a particle spends moving towards negative poloidal angles, since the particle remains in the inner region of the tokamak. Accordingly, the probability distribution function (PDF) of the exit time, P exit , corresponds to the PDF of the negative flights of duration t, P - flight , i.e. To verify this connection, figure 4.13 shows the probability distribution function (PDF) of negative (left panel) and positive (right panel) flights of duration t. As expected, the PDF of negative flights exhibits an algebraic decay, P - flight ∼ t -µ f , with µ f ≈ µ e , where as shown in Fig. 4.9 P exit ∼ t -µe . The PDF of positive flights decays faster, following an exponential scaling (P - flight ∼ e -λt ). This finding has a physical impact in terms of radial transport: a counter-passing particle spends, probabilistically speaking, more time in the inner region than a co-passing particle in the outer region. This leads to an asymmetrically poloidal (and therefore radial) transport in the presence of the chaotic separatrix. Note that, because µ f < 3, the second moment of the PDF of negative flights diverges, ∞ 0 t 2 P - flight dt → ∞. That is, the negative flights are Lévy flights which invalidates the use of the central limit theorem (CLT) as it is customary done in the Brownian random walk model of diffusive transport. On the other hand, within the CTRW, superdiffusive behavior, γ > 1, is a natural consequence of the existence of Lévy flights, µ f < 3. In particular, according to CTRW theory γ = 2/(µ f -1) which for the numerically determined exponent µ f ≈ 2.2 predicts γ = 1.66, a value very close to the numerically observed γ ≈ 1.64. Let us remind that the theory of the Brownian motion relies upon the application of the CLT, which states that the sum of of N i.i.d. random variables {x i } 1≤i≤N is described by a Gaussian distribution in the limit N → ∞, as long as the first and second moments exist, i.e. x i < ∞ and x 2 i < ∞. One can naturally ask what happens in the case where one of the moments (or even both) does not exist, which is our case. Fortunately, there is a generalization of the CLT for this kind of situations, which was formulated by P. Lévy in the 1930s. The Gaussian distribution function as limit of the sum of i.i.d. variables is replaced by the so-called Lévy or α-stable distribution, characterised by long heavy tails and diverging moments. The definition of this variable is not a coincidence. Indeed, the collapse of the rescaled PDFs at different times provides numerical evidence that poloidal transport exhibits self-similar dynamics with anomalous exponent γ which, consistent with the CTRW model, is equal to the numerically determined super-diffusive exponent γ ≈ 1.64. Formally, the observed self-similar evolution implies the existence of a scaling function F satisfying P ∆θ = t -γ/2 F (χ).

It is observed that the scaling function departs significantly from a Gaussian distribution (represented by a dashed grey line for comparison, which is representative of diffusive processes). Also, the asymmetry in the flights is reflected in the asymmetry of the scaling function. Moreover, according to the CTRW, it should exhibit an algebraic decay of the χ < 0 tail of the form F ∼ χ -(α -+1) with α -= µ f -1, a results fully consistent with the numerically obtained values µ f ≈ 2.2 and α -≈ 1.2. For the χ > 0 tail, within the CTRW, the exponential decay of P + f light implies, consistent with the numerical results, α + > 2. This is represented in the insets of figure 4.14 showing the log-log plots of the tails of the PDF. In the previous chapter we have seen how EGAMs can impact directly the transport of particles by modifying the trajectories through a chaotic regime. That study was done in the absence of any turbulent mode. In the present chapter, we focus our analysis on the possible impact of EGAMs on turbulent transport, which can have an effect on the overall confinement. Turbulent transport is one of the main issues to be understood, predicted and controlled, to achieve steadystate operation in future nuclear fusion devices such as ITER, for it is believed to limit the performance and the confinement time. Basically, a particle embedded in a turbulent flow is transported following an eddy, providing an efficient way of transferring energy from one region to another region (non local transport). The eddies are nothing else but the contours of the E × B flow. Therefore, particles will drift following the electrostatic potential iso-contours, with a net radial transport if the eddies break apart after a given time. Otherwise, the particle comes back to the initial position. Assuming that the transport is only due to turbulence and that turbulent transport is diffusive, we can give an ordering for the diffuion coefficient D ≈ v 2 E,r τ c , with τ c a correlation time. Assuming also that the radial drift is dominated by the radial component of the E × B, we can write v E,r ≈ E θ /B 0 . The correlation time of fluctuations is approximated as τ c ≈ L ,c /v th , where L ,c is a parallel correlation length. Therefore, the diffusion coefficient yields D ≈ k 2 θ ρ 2 th (eZ s φ/T eq ) 2 v th L ,c . Taking the typical orderings of a tokamak (k θ ρ th ≈ 0.1, L ,c ≈ 10m, v th ≈ 10 5 ms -1 and eZ s φ/T eq ≈ 0.01, we have D ≈ 1m 2 s -1 . This value is very close to the measured experimental ones, which indicates that the confinement is indeed limited by turbulence.

Chapter 5

Nonlinear interaction between energetic particle modes and turbulence

Because of the important role that turbulence plays in the performance of tokamaks and because energetic particles are also an essential key for the success of fusion, in this chapter we are interested in analysing how turbulence and energetic particles interact with each other. We will show that this interaction can be direct, in the sense that injecting energetic particles can modify the stability of turbulent modes even without any energetic particle mode excitation, but also indirect, in the sense that the excitation of any energetic-particle-driven instability can lead to an exchange of energy between energetic particles and turbulent modes. In this chapter, we focus on the EGAMs, for they are considered as the simplest paradigm of energetic particle modes, but the analysis can be extended to the electromagnetic regime. Since the present chapter is particularly dense, we anticipate its structure. In section 5.1 we give the fundamental physics behind the linear excitation of ITG modes from a fluid approach (subsection 5.1.1) to both improved fluid and gyro-kinetic descriptions (subsection 5.1.2). Section 5.3 is devoted to a brief discussion on how to inject energetic particles in a gyrokinetic code. The remaining sections of the chapter are devoted to explaining how energetic particles can interact with turbulence. For this purpose, in section 5.4 we characterise the reference turbulent state that we have chosen in which energetic particles will be embedded. Afterwards, we proceed with the following reasoning. In section 5.5 we show that energetic particles can transiently reduce, but not suppress, turbulent transport. In section 5.6 we illustrate the excitation and saturation of EGAMs in the presence of radial profiles and turbulent modes.

CHAPTER 5. INTERACTION BETWEEN EPM AND TURBULENCE

In section 5.7 we provide evidences that EGAMs result in an increase of turbulent transport, which is not due to any modification of the radial profiles, but likely to the interaction between EGAMs and ITG modes. After showing the radial structure of the EGAM in section 5.8, detailed analysis based on wavelet transform is presented in section 5.9, proving the nonlinear coupling between EGAMs and ITG modes, responsible for a local exchange of energy.

Destabilization of ion temperature gradient modes

Before analysing the interaction between EGAMs and turbulence, we have considered necessary for the sake of completeness to provide an interpretation of the destabilization of non axi-symmetric modes by an ion temperature gradient (ITG). We will first of all give a heuristic interpretation of the mechanism of the destbilization through a fluid approach and then we will explain the limitations and present the results using a gyro-kinetic approach.

Heuristic interpretation of the the marginal stability in the fluid limit

The quasi-neutrality condition can be derived by means of the variational principle, which states that the action

Ξ = -d 3 xδnφ † + τ d 3 x n eq |φ -φ 0 | 2 (5.1)
is extremum for any variation of the electrostatic potential φ † , where the second term represents the adiabatic response of electrons. Here, τ is the ratio between electron and ion temperatures, namely T e and T i respectively, the electrostatic potential φ is normalized to T i /e, with e the electron charge and n eq is an equilibrium density. The flux surface average of φ is noted φ 0 . The perturbed density δn is obtained from the identity δn = d 3 vδF . Using the adiabatic and non adiabatic decomposition introduced in chapter 2, the Lagrangian for a (n = 0, ω) mode is

L = L ad + L non ad (5.2)
where

L ad = |φ n =0,ω | 2 τ n eq -T i d 3 v∂ H F eq (5.3a) L non ad = |φ n =0,ω | 2 T i d 3 v ω∂ H F eq + n 3 ∂ Pϕ F eq ω -n • Ω (5.3b)
The denominator of the Lagrangian due to the non adiabatic response for the resonant modes (k = 0) may be expressed as ωn • Ω = ω -n 3 Ω d E, where, in circular unshifted magnetic flux surfaces, Ω d = q erB 0 R 0 Ωd and Ωd is a normalized frequency [START_REF] Zonca | Electron fishbones: theory and experimental evidence[END_REF]. In the following, Ωd will be assumed constant, which means that

DESTABILIZATION OF ION TEMPERATURE GRADIENT MODES

no difference is made between trapped and passing particles. Therefore, the non adiabatic term reads (up to a normalization factor)

L non ad = -T i d 3 v E ω ∂ H F eq + Ω -1 d ∂ Pϕ F eq E -E ω (5.4)
where E ω = ω n 3 Ω d . We now assume that the equilibrium distribution function is given by

F eq = n i (2πT i /m) 3/2 e -E
(5.5)

with n i = n i (P ϕ ), T i = T i (P ϕ ) and E = E T i .
The diamagnetic frequency is generalized as follows

Ω -1 dT ∂ Pϕ F eq = Ω * n + E - 3 2 Ω * T F eq (5.6)
where

Ω * n = Ω -1 d 1 n i dn i Pϕ and Ω * T = Ω -1 d 1 T i dT i
dPϕ . Note that the derivative with respect to P ϕ of a quantity decreasing with the minor radius is positive. Therefore, for profiles decreasing with the minor radius we have

Ω * n , Ω * T , Ω * T ⊥ > 0.
The integral is performed with respect to the energy and the adiabatic invariant in order to conveniently analyse the resonance. Therefore, the integral in velocity space is transformed as follows

• • • d 3 v = T 3/2 i 4π m 2 m 2 ∞ 0 dE E 0 • • • du √ E -u
where we have assumed that everything inside the integral is even in v .

The fluid limit can be obtained by Taylor-expanding the denominator up to the second order in E/E ω . The dispersion relation reads (with τ = 1 and n eq = n i )

L (E ω ) ≈ L (E ω ) = E 2 ω + A 1 E ω + A 0 = 0 (5.7)
where

A 1 = Ω ⋆ n - 3 2 Ω ⋆ T + (Ω ⋆ T -1) E (5.8a) A 0 = Ω ⋆ n - 3 2 Ω ⋆ T E + (Ω ⋆ T -1) E 2 (5.8b) and E = 2 √ π R x 3 g 1 (x) e -x 2 dx
(5.9)

E 2 = 2 √ π R x 5 g 1 (x) e -x 2 dx (5.10)
where we have introduced the function The fluid dispersion relation finally reads

g 1 (x) = x 0 e -β(
E 2 ω + Ω ⋆ n - 3 2 E ω + 3 2 Ω ⋆ n + Ω ⋆ T - 5 2 = 0 (5.12)
and the stability curve is given by the equation

Ω ⋆ T = 5 2 -Ω ⋆ n + 1 6 Ω ⋆ n - 3 2 2
(5.13)

From an improved fluid approximation to the gyro-kinetic dispersion relation

Of course, the previous analysis is based on several approximations which facilitate the analytical developments. The first one is that the solution we are looking for is the marginally stable one, i.e. the condition L (E ω ) = 0 provides the solution without any growth rate, satisfying therefore ℑ (ω) = 0. This is an important assumption, since it means that we simply need to look for the condition to have real solutions, and therefore we do not face the problem of looking for the zeros in the complex plane. The second approximation is that we can Taylor-expand the resonance in the fluid limit up to the second order in frequency. This expansion and the subsequent truncation can be done assuming that the mode frequency is great than the characteristic frequency of the particle orbits, namely

|E ω | > E.
This assumption is quite strong and, in addition, when Taylor-expanding one has to take into account that the integrals 5.9 and 5.10 should not be performed in the

DESTABILIZATION OF ION TEMPERATURE GRADIENT MODES

whole velocity space, for the Taylor expansion is valid only within the convergence radius of the truncated series. This problematic has been recently invoked in [START_REF] Migliano | An improved approximation for the analytical treatment of the local linear gyrokinetic plasma dispersion relation in toroidal geometry[END_REF], where the name High Frequency Approximation (HFA) was used to refer to the following Taylor-expansion on the parameter

E/E ω 1 E -E ω = 1 E ω lim n→∞ n l=0 E E ω l (5.14)
Note that in the previous subsection we have truncated this series at n = 2. In Ref. [START_REF] Migliano | An improved approximation for the analytical treatment of the local linear gyrokinetic plasma dispersion relation in toroidal geometry[END_REF], a detailed comparison between the solution obtained from a gyrokinetic code [110] against the one obtained using the HFA is presented and it is explained why this limit is not applicable to relevant tokamak scenarios [START_REF] Romanelli | Ion temperature-gradient-driven modes and anomalous ion transport in tokamaks[END_REF][START_REF] Romanelli | Kinetic theory of the iontemperature-gradient-driven mode in the long wavelength limit[END_REF], in terms of the region of validity for the convergence of the series 5.14. The problematic is schematically represented in figure 5.2a, where we plot in velocity space (v , µ) the limit inside which the expansion 5.14 converges. Basically, if the integrals 5.9 and 5.10 can be performed only inside the red region, we are losing information about the interaction between waves and particles. In Ref. [START_REF] Migliano | An improved approximation for the analytical treatment of the local linear gyrokinetic plasma dispersion relation in toroidal geometry[END_REF] we proposed a different Taylor expansion, which leads to a more accurate description of the wave-particle interaction. The Taylor expansion is actually performed after having applied some elementary algebra to the resonance

1 E -E ω = E ω,r -E -iE ω,i E 2 ω,r + E 2 ω,i 1 -2EEω,r-E 2 E 2 ω,r +E 2 ω,i (5.15)
where E ω,r = ℜ (E ω ) and E ω,i = ℑ (E ω ). The new resonance is slightly modified and we can Taylor-expand

1 1 -2EEω,r-E 2 E 2 ω,r +E 2 ω,i = lim n→∞ n l=0 2EE ω,r -E 2 E 2 ω,r + E 2 ω,i l (5.16) provided that |2EE ω,r -E 2 | < E 2 ω,r + E 2 ω,i .
We refer to this expansion as Improved High Frequency Approximation (IHFA). Figure 5.2b shows the comparison of the two convergence domains for both HFA and IHFA. The colours indicate the exchange of energy between particles and waves, given by the J • E diagnostic as obtained with a local gyro-kinetic code [110]. Negative values indicate stabilising regions (particles take energy from the mode) and positive values indicate destabilising regions (particles give energy to the mode). The thin dashed line is the boundary between stabilising and destabilising regions. It is observed that the IHFA allows one to capture in a more accurate way the exchange of energy. Therefore, even though we employ an expansion of the resonance, we expect to recover a solution which is closer to the gyro-kinetic one than that obtained using the HFA.

The dispersion relation is now obtained by coupling the solution of the gyrokinetic Vlasov equation to the quasi-neutrality condition. Wave-particle energy exchange (colour) together with the convergence domains of HFA and IHFA. For the sake of completeness, the boundary between stabilising and destabilising regions is given by the thin dashed line.

frequency (left panel) and growth rate (right panel) using the gyro-kinetic dispersion relation (dashed line with open circles), the IHFA restricted to its convergence domain using n → ∞ (solid line with circles) and the truncated IHFA with n = 2 (dashed line). For comparison, we plot also the solution using the HFA with n = 1. Further analyses can be done to study the solution of each method based on different parameters. Detailed discussion is reported in [START_REF] Migliano | An improved approximation for the analytical treatment of the local linear gyrokinetic plasma dispersion relation in toroidal geometry[END_REF] and the main idea is that the HFA (the one used usually for the ITG instability) cannot be applied to obtain physical results for the relevant tokamak parameters. In essence, there is a destabilisation due to the ion temperature gradient and the real frequency exhibits different behaviours depending on the approximation employed, with the IHFA very close to the exact gyro-kinetic solution, showing that the mode rotates in the ion diamagnetic direction, which consistent with the fact that it is excited by ions.

The nonlinear coupling between modes

It is to be noted that the previous discussion applies to linear excitation of one single (m, n, ω) mode. This is quite straightforward to understand in the case of analytic calculations where we impose the mode for which we want to analyse the stability and we linearise the gyro-kinetic equations. In gyro-kinetic simulations, we can either perform a filter in Fourier space to select one single mode, or give a whole set of modes. The first method is required if one wants to analyse the linear excitation of one single mode with a nonlinear gyro-kinetic code (similarly to the analysis of chapter 3). GYSELA belongs to this class of codes, intrinsically written to solve the nonlinear equations. The second method is used when the gyro-kinetic code exhibits the possibility to run in linear regime. Codes like ORB5 [83], GKW [110] or GENE [START_REF] Jenko | Nonlinear electromagnetic gyrokinetic simulations of tokamak plasmas[END_REF] belong to this class. In that case, all the unstable modes will be excited, but only the most unstable one will be observed, since it dominates the system. Within the aforementioned linear regime, it is essential to realise that no turbulence exists. Turbulence occurs only in nonlinear regimes. Basically, a wide spectrum of (m, n, ω) modes will be linearly excited. As occurred for the EGAM, each of these modes enters the nonlinear regime, characterising the saturation phase where all the modes are able to interact with each other through the nonlinear terms in the Vlasov-Poisson system. These terms are basically the E × B drift and the nonlinearity in velocity space, namely ∂F ∂t

THE NONLINEAR COUPLING BETWEEN MODES

(N.L.) v E = v E • ∇F = √ g 1 B * ǫ jik b k ∂ j (J 0 • φ) ∂ i F (5.17a) ∂F ∂t (N.L.) v = dv dt ∂F ∂v (5.17b)
One can easily make the nonlinearity appear by decomposing the distribution function into an equilibrium and a perturbation, as we did in chapter 2, and replacing in 5.17a and 5.17b the perturbed distribution function by its expression in terms of the electrostatic potential using the quasi-neutrality condition. For the sake of clarity and for pedagogical purposes, let us consider only the v E nonlinearity with J 0 = 1 together with a simplified expression of the quasineutrality condition using adiabatic electrons and J 0 = 1

1 n eq dvδF = φ -∇ 2 ⊥ φ (5.18)
where we assume that φ is already of first order in perturbation, i.e. there is not equilibrium potential. Integrating over velocity space the Vlasov equation, considering only the nonlinear terms in v E , we have an nonlinear equation for the evolution of the perturbed density ∂n ∂t

(N.L) v E = √ g n eq B * ǫ jik b k ∂ j φ∂ i ∇ 2 ⊥ φ (5.19)
The (k, ω) mode of the perturbed density reads

-iωn kω = √ g n eq B * k ′ ω ′ ,k ′′ ω ′′ ǫ jik b k k ′ i k ′′ j k ′ ⊥ 2 φk ′ ω ′ φk ′′ ω ′′ (5.20)
where

k = k ′ + k ′′ (5.21a) ω = ω ′ + ω ′′ (5.21b)
and where the sum is explicitly written only for the wave-vector and the frequency. For the repeated indexes (i, j, k), the Einstein notation is used. The Fourier mode of the perturbed density is expressed in terms of the Fourier mode of the electrostatic potential

nkω = n eq 1 -k ⊥ 2 φkω (5.22)
Introducing the last expression into equation 5.20 yields

φkω ∝ i ω (1 -k ⊥ 2 ) k ′ ω ′ ,k ′′ ω ′′ ǫ jik b k k ′ i k ′′ j k ′ ⊥ 2 φk ′ ω ′ φk ′′ ω ′′ (5.23)
with the constraints given by equations 5.21. This shows how the amplitude of each Fourier mode of the potential depends on the amplitudes of all the other Fourier modes whose wave-vectors and frequencies satisfy the equations 5.21. This regime where all the existing modes interact nonlinearly with each other is characteristic of the turbulent simulations we present in this chapter. Of particular interest is the nonlinear coupling between (k ′ , ω ′ ) and (-k ′ , ω ′′ ), with |ω ′ω ′′ | ≪ 1. This situation represents the most efficient nonlinear transfer of energy from non axi-symmetric modes to large scale low frequency (m, n) = (0, 0) mode, which is usually referred to as zonal flow [START_REF] Diamond | Zonal flows in plasmas -a review[END_REF]. This nonlinear coupling will be evidenced later on in this chapter regarding the interaction between ITG modes and EGAMs, which is essential to understand the transfer of energy from the EGAMs to the turbulent modes.

Introducing energetic particles in a gyro-kinetic code

In addition to turbulence, in present devices additional heating such as Neutral Beam Injection (NBI) and Ion Cyclotron Resonance Heating (ICRH) leads to the

INTRODUCING ENERGETIC PARTICLES IN A GYRO-KINETIC CODE

formation of long tails in the distribution function, which indicates the presence of energetic particles (EP), having energies larger than those of thermal particles. EP will also be present in future nuclear reactors, where an important population of α particles will be produced during the fusion reactions. References [35,36] are recent reviews on EP and EP modes with appropriate references to previous work on this topic. Energetic particles must be sufficiently well confined in order to transfer their energy, via Coulomb collisions, to thermal particles and sustain this way the nuclear fusion reactions. However, energetic particles excite EP modes. This excitation is due to the inversion of the slope of the distribution function in phase space, i.e. either in real space or in energy space (parallel and/or perpendicular velocity). Intense research has been done in the direction of EGAMs for they have been believed to play a role in the stabilisation of turbulence, as was demonstrated for the standard GAM [START_REF] Hallatschek | Transport control by coherent zonal flows in the core/edge transitional regime[END_REF] and observed in experiments in the AUG tokamak [START_REF] Gd Conway | Mean and oscillating plasma flows and turbulence interactions across the lh confinement transition[END_REF], similar to what is believed for the zonal flow [START_REF] Diamond | Zonal flows in plasmas -a review[END_REF]. However, they have been found to modify nonlocally turbulent transport, for instance when fronts propagate outwards (as observed in ITG turbulence) and couple to static oscillations of radial profiles (characteristic of EGAMs) [43]. Because of the importance of both turbulence and EP and because they naturally co-exist in present and future fusion devices, the present chapter of this dissertation is devoted to characterising the way EP modes and turbulence interact. Because EGAMs constitute the simplest example of EP modes, we focus our effort on the investigation of the fundamental underlying mechanism of the interaction between turbulence and EGAMs.

All the simulations presented in this chapter were performed with the full-F , global, flux-driven gyro-kinetic code GYSELA [?], in the presence of a collision operator as a simplified Lorentz-type operator in parallel velocity space. To avoid strong gradients at the boundaries, radial diffusion is applied in the inner and outer buffer regions by means of a radial diffusion operator. It is to be noted that so far, no other gyro-kinetic code has been used for the analysis of the interaction between energetic particle modes and turbulence.

There are typically two ways to introduce energetic particles in a gyro-kinetic code: either by means of an equilibrium distribution function which depends only on the motion invariants or by means of an external source which builds the distribution function. The first method uses the gyro-kinetic code as an initialvalue solver. Basically, at t = 0 we fix the distribution function of energetic particles. If it does depend only on the motion invariants, it will not evolve until an instability occurs. However, there are situations where defining a distribution function which depends only on the motion invariants is rather complicated. For instance, if there is a radial dependence of temperature and density or if there is a collision operator, defining such distribution functions is quite a challenging task. These kind of scenarios are the most relevant ones though. In essence, the first method where the distribution function depends only on the motion invariants is useful for academic purposes, in the sense that it allows us to compare theoretical predictions with gyro-kinetic simulations, as we did in chapter 3. The second method provides a way to introduce energetic particles independently of the presence of collisions and/or radial profiles which will tend to modify the distribution function towards a canonical or a local Maxwellian.

For our purposes we make use of two sources. A source of thermal energy (S th ) to maintain the radial profiles well above the threshold for the excitation of ITG modes and an energetic particle source (S EP ) that accelerates particles in the parallel direction eventually providing an inversion of the slope of the distribution function in parallel velocity space. The mathematical expressions of these two sources can be found in [48]. The radial envelope of the thermal and energetic particle sources together with the radial envelope of the buffer region are given in figure 5.4. The blue line represents the thermal source and the red line the energetic particle source. The shadowed area indicates the buffer region. Indicated is also the radial position where the analysis of the equilibrium distribution function is done in the following. An example of the dependence of the sources on the parallel velocity is given in figure 5.5. It can be observed that the thermal source accelerates particles around a velocity v ≈ 0, whereas the energetic particle source accelerates and decelerates particles around a velocity v 0 ≈ 2.

The location and the extension of the regions in velocity space where particles are accelerated are critical for the inversion of the slope of the equilibrium distribution function. The EP source used in this simulation, although it represents a model, tries to mimic the generation of energetic particles, making the system evolve towards a state where the slope of the distribution function is positive in a region of phase space. Therefore, the EP source cannot be used to predict the direct impact of external heating on confinement or turbulence properties, but only as a means to introduce energetic particles in the system and analyse their nonlinear interaction with background turbulent modes. Note also that in the simulations the presence of the collision operator is essential to balance the source of energetic particles. The system must therefore be understood as a very complex physical system in five dimensions with injection of energy that is dis- tributed in phase space via collisions and transported in real space via neoclassical and turbulent transport. Part of this energy will be used to excite EGAMs, that can eventually interact with turbulence, modifying this way the transport. The system is closed by the presence of sinks of energy localised at the buffer regions.

DESCRIPTION OF A TYPICAL INITIAL TURBULENT STATE

Description of a typical initial turbulent state

The simulations are run first with the thermal source S th but without the EP source S EP , to reach a statistically turbulent steady state, up to ω c t 0 ≈ 9 • 10 5 . Once this state is achieved, the EP source is switched on to evidence the impact of energetic particles on ITG modes. On the left panel of figure 5.6 we give the radial profiles of the guiding-centre temperature and density, together with the radial envelope of the buffer region. Both temperature and density are averaged over a time window ω c ∆t ≈ 3 • 10 4 . On the right panel of figure 5.6 we plot the radial profile of the turbulent intensity in the stastistically turbulent steady-state, calculated as follows

I turb (r) = 1 ∆t ∆t dt m>1 φm (r, ϕ = 0, t) 2 1/2
(5.24

)
where ∆t is the width of the time window (ω c ∆t ≈ 3 • 10 4 ) where turbulent intensity is calculated, and φm is the Fourier mode in the poloidal direction evaluated at the fixed toroidal position ϕ = 0. Since we are considering a fixed toroidal position, we are considering also the axisymmetric modes. Ideally one should use the Fourier mode φmn for n > 0 in both poloidal and toroidal directions. However, due to numerical and memory constraints, one cannot have a sufficiently well time resolved 3D diagnostic and must use the electrostatic potential at a given toroidal position for radial profiles. Although expression 5.24 is an approximation for the turbulent intensity, the resonant condition of ITG modes q ≈ m/n, ensures that the axisymmetric modes n = 0 do not contribute significantly. The only axisymmetric components that could enter in the expression are related to higher poloidal harmonic GAMs. Since the GAM amplitude decreases with the poloidal mode number m, the main components of GAMs are m = 0 and m = 1 and therefore GAMs do not contribute significantly in expression 5. [START_REF] Lee | Theory of ion-temperature-gradient-driven turbulence in tokamaks[END_REF].

The right panel of figure 5.6 is of particular interest, since the radial pattern of the turbulent intensity is reminiscent of self-organised flow structure called E × B staircase, already observed in previous GYSELA simulations [115] and more recently identified in Tore Supra tokamak by means of high-resolution fastsweeping X-mode reflectometry [116]. This staircase structure is characterised by regions with increased turbulence, the separation between those regions being of the order of 20 -30ρ i [116].

Onset of the energetic particle source: initial reduction of turbulence

Once the statistically turbulent steady-state is reached, the EP source is switched on. After switching on the EP source, the simulations are run long enough (up to ω c t ≈ 1.3 • 10 6 ) to reach a new statistically steady state, where both ITG modes and EGAMs coexist in nonlinear regime.

As explained in Ref. [43], the first effect of the EP source is to accelerate par- Figure 5.7: (Left) Radial dependence of the turbulent intensity as given by expression 5.24 for the three simulations used in this chapter. (Right) Time evolution of the temperature gradient after switching on the energetic particle source, averaged over the radial region 50ρ i < r < 100ρ i . ticles from thermal to suprathermal parallel velocities. This leads to a depletion of the distribution function around the thermal velocity, where particles resonate with ITG modes. The depletion results in a stabilisation of ITG modes and therefore in a reduction of turbulent transport. In the same direction, it is important to notice that recent works have been conducted to analyse more in detail the reduction of turbulence by injection of energetic particles in electrostatic regimes. In essence, the mechanism is reduced to the modification of the resonance between particles and the ITG modes, as established in [START_REF] Di Siena | Fast-ion stabilization of tokamak plasma turbulence[END_REF].

ONSET OF THE ENERGETIC PARTICLE SOURCE: INITIAL REDUCTION OF TURBULENCE

On the left panel of figure 5.7 we plot the radial dependence of the turbulent intensity for the three amplitudes of the EP source that we use. For comparison, we give also the radial dependence of the turbulent intensity before switching on the EP source, whose radial profile is given by the shadowed region. It is clearly observed the overall decrease of turbulent intensity where the EP source is applied. On the right panel of figure 5.7 the time evolution of the temperature gradient averaged over the radial region 50ρ i < r < 100ρ i is plotted. It can be observed that the reduction of turbulence is accompanied by an increase of temperature gradient within the time interval 10 6 < ω c t < 1.1 • 10 6 . The observed reduction is only temporary since the turbulent transport is afterwards increased from a time ω c t inc that depends on the amplitude of the EP source. This can be to either a re-organisation of the system (in both real and velocity space) in the presence of the EP source, or a modification of turbulent modes when EGAMs are excited. In the remainder of this chapter we investigate the possibility that EGAMs contribute to the increase of turbulent transport.

Excitation and nonlinear saturation of EGAMs in the presence of turbulence

The EP source implemented in GYSELA and detailed in the previous section is able to invert the slope of the distribution function. The larger the injection of parallel energy is, the larger the slope of the distribution function is. The inversion of the slope of the distribution function leads to the excitation of the EGAM [42], which can be observed when plotting the time evolution of the flux-surface averaged electrostatic potential at a given radial position (in our case we consider the mid radial position). This is shown in figure 5.8a for the three simulations used in this chapter. It can be observed that the excitation occurs earlier when increasing the amplitude of the EP source, i.e. when increasing the drive of the EGAM. This is due to the fact the evolution of the slope of the equilibrium distribution function is accelerated and crosses the zero value earlier. Also, the amplitude of the electrostatic potential during the nonlinear saturation is increased with the amplitude of the EP source. Further analysis will be done in the future to analyse in detail the dependence of the amplitude of the nonlinearly saturated potential in terms of the drive in the absence of turbulence. After filtering out the time averaged signal, time Fourier transform of the flux-surface averaged electrotatic potential reveals the existence of a dominant frequency, as shown in figure 5.8b, which corresponds to the EGAM frequency ω EGAM . This frequency can be obtained from analytic theory, as reported in [44,1], which has been verified against previous gyro-kinetic simulations [44]. Also benchmarking between gyro-kinetic codes [45] has revealed that this frequency is robust and depends only on the fraction of energetic particle pressure, decreasing down to roughly half the standard GAM frequency. The EGAM frequency can also be observed in figures ?? and 5.7b from ω c t 1.1 • 10 6 . Note that although the slope of the distribution function is also inverted (as we will see later) for the case S EP /S th = 0.312, the excitation of the mode is not visible in figure 5.8a. This indicates the presence of additional damping mechanisms different from the Landau damping.

Increase of turbulent transport in the presence of EGAMs

In figure 5.9a we plot the slope of the equilibrium distribution function, ∂ E F eq , as a function of time for the three amplitudes reported in the paper and evaluated at the resonant velocity, where ω is the measured frequency of the EGAM, i.e. v G = qR 0 ω EGAM . Time-averaged quantities over an EGAM cycle are represented by thick lines. Also, the time-averaged over an EGAM cycle of the heat flux due to the E × B drift velocity is plotted for the sake of completeness, represented by the shadowed area and calculated as follows where F is the guiding-centre distribution function and v E×B,r is the radial component of the E × B drift velocity. For the sake of completeness, the equilibrium distribution functions without energetic particles (thin dashed line), with S EP /S th = 0.437 (thick dashed line) and with S EP /S th = 0.625 (solid line) just before the excitation of EGAMs are given in figure 5.9b, where the inversion of the slope is clearly visible for co-and counter-passing particles, leading to a bumpon-tail distribution. Note that the distribution function for S EP /S th = 0.312 is not given in the figure for clarity. It is observed that increasing the amplitude of the EP source increases also the slope of the equilibrium distribution function and therefore the drive. Note that as mentioned earlier, the slope is inverted for the three amplitudes, even if the EGAM is only visible for two of them. Time Fourier transform of the slope of the distribution function (not shown here for the sake of clarity) reveals oscillations at the EGAM frequency. The oscillations are larger when the amplitude of the EP source is increased, in other words when the EP drive is stronger. The most interesting feature related to the evolution of the distribution function that we observe in these simulations is the evolution of the slope of the distribution function during the nonlinear saturation of the mode. The saturation of EGAM in the absence of turbulence and without radial profiles of density and pressure is due to a total flattening of the distribution function in parallel velocity space, as observed in previous gyro-kinetic simulations [42]. However, in the simulations presented in this chapter, with turbulence and with radial profiles of density and pressure, we observe that the flattening is only partial. This indicates that either the presence of turbulence or the finite radial gradients result in an additional sink of energy. To analyse the effects of radial profiles on the EGAM excitation and saturation we have performed a simulation similar to the one where S EP /S th = 0.625, but without turbulence. This is done by filtering out every n = 0 mode after solving the quasi-neutrality equation. We have selected an amplitude of EP source so that the maximum slope of the distribution function is close to the one obtained in the presence of turbulence, while keeping the same radial gradients. This way the simulation is equivalent in terms of EP drive and diamagnetic damping. The time trace of the flux-surface averaged electrostatic potential is plotted in figure 5.10a in logarithmic scale. The solid red line corresponds to the simulation with turbulence, whereas the thin dashed line corresponds to the simulation without turbulence. Both signals have been plotted using the same time vector, for comparison. It can be observed that in both cases EGAMs are excited exhibiting similar linear growth rates. However, the absence of turbulence leads to an increased nonlinearly saturated amplitude. In addition, we plot the time evolution of the slope of the equilibrium distribution function at the resonant parallel velocity in figure 5.10b and represented by the magenta curve. As in figure 5.9a, time-averaged quantities over an EGAM cycle are represented by thick lines of the same colour. Several comments can be made based on these figures. First, the reduction of the slope during the nonlinear saturation is more dramatic in the absence of turbulence. Second, the time evolution of the slope evolves smoothly in the presence of turbulence, compared to the situation where turbulence is suppressed, characterised by a more pronounced bursty behaviour. Finally, the amplitude of the oscillations of the distribution function is increased when filtering out the non axisymmetric modes. These observations in the electrostatic potential and the distribution function indicate that in the absence of turbulence kinetic effects play a more pronounced role in the saturation of EGAM, whereas in the presence of turbulence a channel of energy from EGAMs to turbulence is possible.

Q E×B = 1 4π 2 dθdϕ J v dv G dµ 1 2 mv 2 v E×B,r (J 0 • F ) (5.
A second indication of the interaction between ITG modes and EGAMs, is that the nonlinear saturation of the EGAM, visible mainly because of the partial flattening of the distribution function, is accompanied by a sharp increase of the E × B heat flux (around ω c t = 1.1 • 10 6 ), which occurs in the same time scale as the nonlinear saturation of the EGAM. This time evolution of the E × B heat flux indicates a possible nonlinear interaction between ITG modes and EGAMs, which can occur nonlocally by means of coupling between outwards propagating fronts linked to ITG modes and static oscillations of radial profiles linked to EGAMs, as invoked in [43]. A local mechanism might also exist by means of a parametric decay, as in the case of standard GAM and ITG modes [START_REF] Zonca | Radial structures and nonlinear excitation of geodesic acoustic modes[END_REF], represented by a resonant three-wave interaction dominating over the wave-particle interaction.

Note that two additional mechanisms could be invoked for the modification and increase of turbulence. First, since the source is depleting the population resonating with ITG modes, one could expect that the nonlinear saturation of EGAMs would redistribute particles in velocity space and therefore feed again the ITG instability. However, as can be observed from figure ??, after the EGAM saturation, the value of the equilibrium distribution function remains below 50% the initial value (i.e. dµF eq / dµF eq (t = 0) < 0.5). Since the reduction of turbulence occured with a reduction of more than 30% (i.e. dµF eq / dµF eq (t = 0) < 0.7) one can conclude that the redistribution of particles in velocity space during

CHAPTER 5. INTERACTION BETWEEN EPM AND TURBULENCE

the nonlinear saturation of EGAM is not enough to explain the increase of turbulent intensity. Second, in the calculation of the E × B heat flux we are considering also the axisymmetric modes. Therefore, one could think that those modes also contribute to the increase of the heat flux. Nevertheless, it was reported in [43] that axisymmetric modes excited by energetic particles do not contribute directly to the turbulent transport estimated by the E × B heat flux. For these reasons, in the remainder of this chapter we look for evidences of a local nonlinear coupling between EGAMs and non axisymmetric modes that could explain the impact of EGAM saturation on turbulence.

Radial structure of EGAM and modification of turbulence intensity and transport

The structure of EGAM has been little explored in the previous litterature. EGAM is believed to have a global structure [START_REF] Fu | Energetic-particle-induced geodesic acoustic mode[END_REF] and has been found in simulations with imposed background distribution function to exhibit a structure constrained by the Dirichlet boundary conditions [42], with a frequency that is constant in the whole radial domain, as in experiments [39, 38, ?]. In addition, in the presence of a localised source, the EGAM has been found to be trapped where the drive is strongest with an exponential damping due to the continuous GAM [50]. In our simulations, the source cannot be strongly localised, in order to minimise the direct impact on the radial gradients and therefore on the ITG modes. However, due to the way the EP source is implemented, the EP drive, i.e. the positiveness of ∂F eq /∂E, does not follow the same radial dependence as the EP source [43].

In that way, we can infere that the effective EP source is localised in a region 0.4 < r/a < 0.8 [43], which is equivalent to 60ρ i < r < 120ρ i . The radial structure of the EGAM is analysed in this chapter using the electrostatic potential φ. We determine the frequency of the mode from the time dependent flux-surface averaged electrostatic potential φ 00 (r mid , t) evaluated at the mid radial position r mid .

Owing to the poloidal and toroidal structures of the EGAM-like potential, we filter out all the non-axisymmetric toroidal components and keep only the m = 1, 2 poloidal components that oscillate at the EGAM frequency. Once the filter applied, the sine and cosine components of the m = 1 and m = 2 modes can be obtained. These components are plotted in figure 5.11 for the simulation with an EP source S EP /S th = 0.625. It can be observed that the sin θ component dominates over the cos θ component and is mainly localised in the region 60ρ i < r < 120ρ i , where the EP drive is strongest. It is also observed that the cos θ, sin 2θ and cos 2θ are of the same order, extending slightly radially further in. The localisation of the sin θ component results in the localisation of the radial electric field calculated as E r ≈ -∂φ 00 /∂r at the EGAM frequency. The oscillating radial electric field can therefore be able to interact with ITG modes, modulating and modifying this way the amplitude of the turbulent intensity and the E × B heat flux. The modification of the radial structure of turbulent intensity is analysed by using the expression 5.24 evaluated at the end of the simulations within a suf- .12: For the three simulations used in this chapter, plotted is the turbulent intensity at the end of the simulations as measured by 5.24 for a sufficiently long time window where the system exhibits a statistically steady state. For reference, the case without EP source is given (dashed thin line). The shadowed region represents the radial localisation of the EGAM. ficiently long time window to ensure that the system has reached a statistically steady state. The radial dependence of the turbulent intensity is plotted in figure 5.12, where also the turbulent intensity in the absence of EP is given for reference (thin dashed line), as was done in figure 5.7a. In addition, we plot the radial dependence of the EGAM (shadowed region). Several comments can be made based on this figure. First, the staircase structure observed in the absence of EP source disappears in the presence of EP. The plasma staircase plays an important role in organising the energy transport, as the long-range radial propagation of energy by means of avalanches is interrupted statistically by avalanche-staircase interplay [115,116]. When the staircase structure disappears, energy transport that would be interrupted can expand over a larger radial domain, as observed in [43]. Second, in the region where the EGAM is not excited, the turbulent intensity is close to the one in the absence of EP, for any amplitude of the EP source. Third, the main modification occurs in the region where the EGAM is localised. It is to be noted that the turbulent intensity remains lower than the reference intensity for the two cases where the EGAM amplitude is weaker (blue lines). When the EGAM amplitude is increased (red line), the turbulent intensity is also increased, but remains of the same order as the intensity in the absence of EP. The main quantity to be analysed within the framework of the transport in a tokamak is the radial E × B heat flux, given by equation 5.25. In figure 5.13a, the time-averaged E × B heat flux over a time window ω c ∆t > 10 5 is plotted as a function of the minor radius. It can be observed that overall the heat flux when the EP drive is stronger is increased with respect to the simulations where the amplitude of the EP source is reduced.

RADIAL STRUCTURE OF EGAM AND MODIFICATION OF TURBULENCE INTENSITY AND TRANSPORT

NONLINEAR COUPLING BETWEEN EGAMS AND ITG MODES

Nonlinear coupling between EGAMs and ITG modes

Note that in figure 5.13a we plot time-averaged quantities. However, Fourier analysis shows that the E × B heat flux oscillates at the EGAM frequency. The amplitude of the Fourier mode is represented in figure 5.13b as a function of the frequency. The dominant peak corresponds to the EGAM frequency. For completeness, the first and second harmonics are also included in the figure.

Since the E × B heat flux is a nonlinear quantity, those oscillations must be due to the nonlinear interaction between modes satisfying the matching conditions

ω 1 + ω 2 = ω EGAM and k 1 + k 2 = k EGAM .
The EGAM is mainly represented by the mode φ 00 oscillating at the EGAM frequency. Therefore, we need to look for the possibility to have non axisymmetric modes with opposite poloidal and toroidal mode numbers and whose frequencies sum up the EGAM frequency. This has been by means of a bispectral analysis of the toroidal and poloidal Fourier modes of the electrostatic potential at the mid radial position, called in the following φm,n . The bispectral analysis used in this chapter is based on a wavelet decomposition in time. The term wavelet was introduced by the French geophysicist Jean Morlet [START_REF] Morlet | Wave propagation and sampling theory. part i: Complex signal and scattering in multilayered media[END_REF] and it is based on a generalized filter of a given signal f (x), as explained in [START_REF] Frisch | Turbulence[END_REF]. The given signal can be decomposed, without any loss of generality, into two components [START_REF] Obukhov | Spectral energy distribution in a turbulent flow[END_REF] f

(x) = f > K (x) + f < K (x) (5.26)
where

f < K (x) =
k≤K fk e ikx (5.27a)

f > K (x) =
k>K fk e ikx (5.27b) which are not Fourier transforms, since they are still functions of x, but depend also on another variable K, which is what we call scale variable. The generalization of the decomposition 5.26 to filters of arbitrary shape is called the wavelet transform. Here we make use of the wavelet transform as was done in [START_REF] Van Milligen | Nonlinear phenomena and intermittency in plasma turbulence[END_REF] and explained in previous works [START_REF] Meneveau | Analysis of turbulence in the orthonormal wavelet representation[END_REF][START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF][START_REF] Hudgins | Wavelet transforms and atmopsheric turbulence[END_REF]. For this purpose a Morlet wavelet function is chosen, having the following expression

ψ ω (t) = ω 2π e iωt e -1 2 ( ωt 2π ) 2 (5.28)
where ω ∈ R is the frequency representing here the wavelet scale. The wavelet transform of the electrostatic potential φm,n , similarly to the Fourier transform, is defined as

W φm,n (ω, t) = R φm,n (τ ) ψ ω (τ -t) dτ (5.29)
The wavelet transform is calculated considering that it is a convolution product. Therefore, the Fourier transform of the wavelet transform is related to the

w (ω 1 , ω 2 ) = dt W φ0,0 (ω, t) † m,n W φm,n (ω 1 , t) W φ-m,-n (ω 2 , t) (5.30)
where the integral is in practice performed over a finite time interval, W † represents the complex conjugate of the quantity W and the frequencies ω 1 and ω 2 satisfy the matching condition ω 1 + ω 2 = ω. If a nonlinear interaction exists between two waves of frequencies ω 1 and ω 2 and a third one of frequency ω this must be observed by the presence of straight lines intersecting the vertical axis at the frequency ω.

The analysis has been performed only for the simulation with the largest amplitude of energetic particle source. Four time windows are analysed. They correspond to (1) the phase where turbulence is not reduced yet, (2) the linear excitation of EGAMs until the end of the nonlinear saturation, (3) the phase just after the nonlinear saturation of EGAMs and the increase of turbulence and (4) the steadystate phase where EGAMs and turbulence coexist. A structure characterised by a straight line crossing the point (ω 1 = 0, ω 2 = 0) is clearly observed before the onset of EGAMs (figure 5.14a), meaning that the main interaction occurs between non axisymmetric modes and the zonal structure oscillating at ω ≈ 0, which is the well-known mechanism of self-regulation of turbulence by stationary zonal flows [START_REF] Miyato | Global structure of zonal flow and electromagnetic ion temperature gradient driven turbulence in tokamak plasmas[END_REF][START_REF] Miyato | Study of a drift wave-zonal mode system based on global electromagnetic landau-fluid itg simulation in toroidal plasmas[END_REF][START_REF] Diamond | Zonal flows in plasmas -a review[END_REF]. During the linear excitation and nonlinear saturation of EGAMs (figure 5.14b), the same structure characterising the nonlinear interaction between ITG modes and zonal flows is still present, together with a second one characterising the nonlinear interaction between non axisymmetric modes oscillating at ω 1 and ω 2 and a zonal component oscillating at a frequency ω = 0. This non-zero frequency is evidenced by the analysis of the wavelet bispectrum at later times, as shown in the next figures. In particular, just after the saturation of EGAMs (figure 5.14c), two nonlinear interactions are clearly observed: one is the standard interaction between non axisymmetric modes and zero frequency zonal modes and the other one is the nonlinear interaction between non axisymmetric modes and zonal components at EGAM frequency (highlighted by two horizontal yellow dashed lines). When analysing the nonlinear interaction during the steadystate (figure 5.14d), the standard self-regulation of turbulence by zonal flows that should be visible by the straight line crossing the point (ω 1 = 0, ω 2 = 0) disappears and the dynamics is completely dominated by the nonlinear interaction between three waves satisfying ω 1 + ω 2 = ω EGAM . It is to be noted that analysis of different energy channels during the nonlinear interaction between turbulence and zonal flows was reported in the past using a global electromagnetic Landau

CHAPTER 5. INTERACTION BETWEEN EPM AND TURBULENCE

fluid code [START_REF] Miyato | Global structure of zonal flow and electromagnetic ion temperature gradient driven turbulence in tokamak plasmas[END_REF][START_REF] Miyato | Study of a drift wave-zonal mode system based on global electromagnetic landau-fluid itg simulation in toroidal plasmas[END_REF], evidencing that a part of the energy of the oscillatory zonal flows (GAMs) is transferred back to the turbulence via the poloidally asymmetric component of the pressure perturbation. This was the case for nonlinearly self-generated zonal flows in the absence of energetic particles. However, in the present paper, the generation of oscillatory zonal flows is done by means of energetic particles. Moreover, all the modification of the distribution function by the energetic particle source is essential for the initial reduction of turbulence, as reported in [43]. Therefore, although the results presented in this section are based on a three-wave interaction of the electrostatic potential as could be done in fluid codes, the gyro-kinetic approach reveals essential to evidence the energy transfer from energetic particles to turbulence via EGAMs.

Chapter 6

Towards the low-frequency electro-magnetic regime: the tearing mode

THE TEARING MODE IN TOKAMAK PLASMAS

The results shown in the previous chapters are restricted to the so-called electrostatic regime, where no magnetic perturbations are allowed to grow. The natural extension of this work is the analysis in the electro-magnetic regime, where other classes of instabilities can develop or where energetic particles can exhibit different behaviour with respect to the electrostatic one. The purpose of this chapter is to provide the reader with the fundamental understanding of some of the electro-magnetic instabilities observed in tokamaks, study how they are excited, how they interact with turbulence and with particles. It is to be noted that some electro-magnetic modes present in a tokamak are not necessarily excited by the particles, but are introduced externally to control the stability of the plasma. This is the case of the Resonant Magnetic Perturbations (RMP) [], which are used to control the Edge Localized Modes (ELMs) []. We will devote a small section in this chapter to explain how these RMP can induce a transport of energetic particles.

The tearing mode in tokamak plasmas

In tokamaks, an ideal magnetic equilibrium consisting of axisymmetric nested flux surfaces can be perturbed on a rational flux-surface by a radial magnetic component δB r (θ, ϕ), generated by the so-called unstable tearing mode [59,60,61]. The tearing mode is linearly driven by the gradient of the equilibrium parallel current around the rational surface. The growth of this linear instability is clearly unphysical when the amplitude of δB r becomes significant. At that point, the situation must be regularized and the instability saturated by non-ideal effects leading to magnetic re-connection [128], characterized by the destruction of the resonant surfaces into magnetic islands. These instabilities are found in both astrophysical and laboratory plasmas [62]. In the case of laboratory plasmas (like tokamak plasmas), the presence of magnetic islands has strong consequences for the core confinement [129] This is due to the fact that the radial perturbation connects different equilibrium radial flux surfaces, resulting in an increase of transport which subsequently leads to the flattening of the radial profiles and therefore to an overall loss of core confinement and eventually a disruption [63]. Because of the impact that magnetic islands have on the tokamak performance, it is essential to analyse their stability and their interaction with thermal and energetic particles as well as with turbulence.

In order to better understand why this class of instability is driven by the radial gradient of the parallel current around the resonant surface, it is useful to have a look back to expression 2.48 in section 2.4.1, where it is shown tha the numerator of the resonant response is written in terms of the radial gradient of the electron mean velocity. As explained therein, this leads to an instability driven by the radial gradient of the equilibrium parallel current. We can provide a further deep insight into the excitation of the tearing mode by deriving a simplified gyro-kinetic dispersion relation, which captures the main physics behind the interaction between the tearing mode and the particles (electrons in this case). In essence, we assume that the current is carried out only by deeply passing elec-CHAPTER 6. ELECTRO-MAGNETIC REGIME trons. This means that we can make use of the parallel current given by the expression 2.66 and inject it into Ampère's law 2.50. Using normalized quantities and projecting onto one single (m, n, ω) Fourier mode, we have

δ2 e ∇2 ⊥ δ Â m,n,ω (x) = 2 Teq v2 ω -ω⋆g ω -k v -ωD δ Â m,n,ω (x) (6.1)
where δe = δ e /R 0 and δ e = m e /µ 0 e 2 n 0 is the electron skin depth. In the following, the ˆsymbol will be dropped for the sake of simplicity and all quantities are assumed to be normalized. The dispersion relation of the tearing mode is obtained integrating equation (6.1), but for this purpose we need to give convenient equilibria that will allow us to perform analytic calculations. For this purpose, an intermediate step is the explicit calculation of ω ⋆g following the discussion of chapter 2. for this purpose, we need to express the electron mean velocity in terms of the equilibrium. Indeed, the equilibrium parallel current is directly related to the electron mean velocity through the equation J eq = eZ e u e , where Z e = -1. At the equilibrium, Ampère's law must also be satisfied, which leads to the relation between the safety factor and the electron mean velocity

q (r) = δ 2 e ρ e R 0 Z e n eq r 2 r 0 dr ′ r ′ u e (r ′ ) (6.2)
where all quantities are conveniently normalised in order to bring all the dimensions to the front of the right-hand side. Therefore, the differential equation (6.1) can be written using only the radial profile of the safety factor, which yields the linear dispersion relation of the tearing mode. Nevertheless, solving analytically this dispersion relation in the whole radial domain for general profiles of density, temperature and electron velocity is rather arduous. Therefore the tearing mode dispersion relation is solved by splitting it into two linear differential equations that will be solved in two different regions of the radial domain: an ideal outer region and a narrow resonant layer. For each region we will make some assumptions. In the resonant layer, |k | ≪ 1 and therefore the whole resonance ω -k v must be kept leading to the ordering ω ≈ k v th , however in the ideal region, |ω| ≪ |k v |, |ω ⋆g |. The two solutions are matched by considering the so-called constant-ψ approximation [60] that assumes that the perturbed parallel scalar potential of the magnetic field is constant inside the non-ideal layer.

The tearing mode equation in the outer (ideal MHD) region

In the ideal region where ω ⋆g ≫ ω and k v ≫ ω, ω D the expression in between brackets of equation (6.1) reads

v 2 ω -ω ⋆g ω -k v -ω D = v 2 ω ⋆ k v (6.3)

THE TEARING MODE IN TOKAMAK PLASMAS

Since ω ⋆ (resp. ω ⋆u ) is an even (resp. odd) function in v , the only remaining term in the numerator of expression (6.3) is ω ⋆u . This leads to the differential equation

δ 2 e ∇ 2
⊥ δA m,n,ω (x) = -ρ ⋆ nn eq q rk du e dr δA m,n,ω (x) (6.4) Using the relation between the safety factor and the mean electron parallel velocity given by equation (6.2) and assuming again that the density gradient does not affect significantly the parallel current gradient we can write

d 2 A (r) dr 2 + 1 r dA (r) dr - m 2 r 2 A (r) = Λ (r) A (r) (6.5)
where we have noted A ≡ δA m,n,ω for simplicity and

Λ (r) = -3q ′ (r) + 2rq ′2 (r) /q (r) -rq ′′ (r) r m n -q (r) (6.6)
with q ′ and q ′′ representing respectively, the first and second derivatives with respect to r of the safety factor. We consider that q > 1, q ′ > 0 and q ′′ > 0 for all r.

We will also suppose that there is a radial position r s such that q (r s ) = m/n. This means that for r < r s ⇒ q (r) < m/n and for r > r s ⇒ q (r) > m/n. Then, for r < r s (resp. r > r s ) we have the inequality Λ (r) < 0 (resp. > 0), which means that for A > 0 the second derivative must be negative (resp. positive) and therefore the solution is concave (resp. convex). In addition, note that lim r→r ± s Λ (r) = ±∞. This means that the first derivative is discontinuous at the position r = r s , which allows us to define the parameter ∆ ′ [60]

∆ ′ = lim ε→0 A ′ (r s + ε) -A ′ (r s -ε) A (r s ) (6.7)
where A ′ denotes the derivative with respect to r of A . This parameter represents the jump of the solution accross the resonant surface and it will be used as a matching parameter between the solution in the inner region and the solution in the outer region. Equation (6.5) can be solved using the shooting method [START_REF] Daintith | A Dictionary of Computing[END_REF] for a general q profile, but in the following we use a q profile of the Wesson-type [START_REF] Wesson | [END_REF]

q(r) = q a r 2 /ǫ 2 1 -(1 -r 2 /ǫ 2 ) ν+1 (6.8)
where q a is the value of the safety factor at the position r = a and ν is a parameter that controls the current density peaking. The solution of equation 6.1 is given in figure 6.1 for three values of ν. It is observed how the position of the inflection point, which corresponds to the position of the resonant surface, is modified when modifying the safety factor profile. This is essential when solving Ampère's law in the inner region, because the position of the resonant layer determines the fraction of trapped particles and has therefore an impact on the stability, as we will see in the following.

The tearing mode equation in the inner region

In the region around the resonant surface k ≪ 1 and therefore the resonance in the expression in between brackets of Eq. (6.1) is kept as it is. Due to parity reasons, the contribution from ω ⋆u in the numerator can be neglected, so that the tearing mode equation reads

δ 2 e ∇ 2 ⊥ δA m,n,ω (x) = 2 T eq v 2 ω -ω ⋆ ω -k v -ω D δA m,n,ω (x) (6.9)
We assume that the width of the resonant layer is 2∆. Integrating the equation over this region yields

δ 2 e ∆ -∆ d 2 A dx 2 dx = A (0) ∆ -∆ Λdx (6.10)
where we have again noted A ≡ δA m,n,ω and the so-called constant-ψ approximation has been made, assuming that the solution is constant in the resonant layer. In addition, to perform the integration in the resonant layer the assumption of slab geometry is made, such that ∇ 2 ⊥ ≡ d 2 /dx 2 , where x = (r -r s ) /R 0 represents the distance to the rational surface at the radial position r s , normalized to R 0 . This assumption is justified due to the small width of the resonant layer. Moreover, the coefficient Λ is now given by the more general expression

Λ = 2 T eq v 2 ω -ω ⋆ ω -k v -ω D (6.11)
The integral on the left hand-side of Eq. (6.10) gives ∆ ′ A (0). For the integral on the right-hand side we perform a change of variable v ′ = -v and x ′ = -x

THE TEARING MODE IN TOKAMAK PLASMAS

to the negative domain of the space integral. In addition, since the perturbed parallel current is localised in a very narrow region around the resonant surface, we can extent the integral with respect to the radial distance to infinity without introducing significant errors. Under these assumption, the collisionless tearing mode dispersion relation in the magnetic limit is expressed in normalized units as

δ 2 e ∆ ′ = 8n eq T 5/2 eq √ π +∞ 0 dx ∞ -∞ v 2 dv ∞ 0 dv ⊥ v ⊥ e - v 2 +v 2 ⊥ Teq ω -ω ⋆ ω -k (x)v -ω D (6.12)
In this expression, the parallel wave vector can be expressed as k (x) = k y /L s x = k ′ x, with L s the magnetic shear scale length and k y = m/r, and the diamagnetic frequency is expressed as

ω ⋆ = ω ⋆n 1 + η e v 2 - 3 2 (6.13)
where v 2 = v 2 + v 2 ⊥ , ω ⋆n = ρ ⋆ nT eq q/2r n ′ eq /n eq and η e = n ′ eq /n eq / T ′ eq /T eq . After some algebra, decomposing the frequency into real and imaginary parts ω = ω r + iγ, the real frequency and growth rate read

ω = ω ⋆n 1 + 2T eq - 3 2
n eq (6.14a)

γ = 1 4 √ π T eq n eq δ 2 e |k ′ |∆ ′ (6.14b)
In addition, it can be shown that the result does not depend on ω D as long as only deeply passing electrons are considered [START_REF] Zarzoso | Gyro-kinetic theory and global simulations of the collisionless tearing instability: the impact of trapped particles through the magnetic field curvature[END_REF]. When considering trapped particles, one has to keep in mind the ordering (for deeply trapped particles) Ω 2 ≫ ω, Ω 3 , where ω is the tearing mode frequency. In order to satisfy the resonance condition, the only possibility is n 2 = 0, which reduces to a bounceaverage. Then, if a resonant effect is due to trapped electrons, it must satisfy ω = n 3 Ω 3 , where Ω 3 is proportional to the magnetic drift frequency ω D . If we neglect ω D for trapped electrons, their effect is not resonant anymore and therefore no transfer of energy between tearing mode and trapped electrons occur. This implies that the integral 6.12 should be performed in the passing domain. This integration leads to a modified growth rate γ trapping γ no trapping = 1 + 2ǫ (6.15)

Hence, trapping electrons tend to increase the growth rate. Physically, this can be understood as follows [START_REF] Peter | Trapped electron modifications to tearing modes in the low collision frequency limit[END_REF]. The free energy for the tearing mode to be excited comes from the radial gradient of the equilibrium parallel current. This free energy is encapsulated in the stability parameter ∆ ′ . Therefore, the free energy comes from the ideal region outside the resonant layer. This energy is subsequently transferred to particles within the resonant layer, where kinetic effects 100 CHAPTER 6. ELECTRO-MAGNETIC REGIME must be retained. The difference between the free energy that the mode takes from the ideal region and the energy that the mode transfers to the particles in the resonant layer is the energy available for the mode to grow. Since deeply trapped electrons do not contribute to this transfer of energy, there is an increase of energy available for the mode to grow, resulting in an increased growth rate. A comparison between this analytic prediction and the one obtained from global gyrokinetic simulations for different values of the stability parameter is presented in figure 6.2 using the code GKW [110], showing a good quantitative agreement. The fact that trapped electrons do not contribute to the tearing instability can be verified by means of the energy exchange diagnostic. This is illustrated in figure 6.3, where the J • E at the resonant position is plotted in velocity space using flat profiles and no magnetic drift in the case where all electrons are passing (top panel) and in the case where trapping effects are included (bottom panel). It is clearly observed that trapped electros do not exchange energy with the tearing mode, which contributes to the increase of the growth rate of the instability.

Impact of the magnetic field curvature on the tearing instability

The magnetic field of a tokamak is characterised by a curvature κ = (b • ∇) b responsible for a drift proportional to v 2 . This drift introduces toroidal couplings between different poloidal harmonics [START_REF] Connor | Resonant magnetohydrodynamic modes with toroidal coupling. part i: Tearing modes[END_REF][START_REF] Fitzpatrick | Stability of coupled tearing modes in tokamaks[END_REF]. The effect of the curvature of the magnetic field has also been analysed using an effective or averaged curvature which oscillates in the poloidal direction [START_REF] Kotschenreuther | Nonlinear dynamics of magnetic islands with curvature and pressure[END_REF], whose effect was found stabiliz- ing for the large scale tearing mode [START_REF] Glasser | Resistive instabilities in a tokamak[END_REF] In addition, this curvature is responsible for the toroidal branch of the ITG instability or for that of the electron temperature gradient (ETG) instability [START_REF] Zielinski | Electromagnetic electron temperature gradient driven instability in toroidal plasmas[END_REF]. In the case of the tearing mode, it has been demonstrated using a fluid approach that similar interchange-like mechanisms can exist [START_REF] Nasr | Interchange destabilization of collisionless tearing modes by temperature gradient[END_REF]. In particular, small scale (micro-) tearing modes that are stable in the absence of curvature can be destabilised by the electron temperature gradient through the magnetic curvature. Also, large scale tearing modes have been found to be further destabilised when the magnetic curvature is taken into account. The impact of the curvature can be analysed using global gyro-kinetic simulations. This is shown in figures 6.4a and 6.4b with flat profiles and with temperature gradient, respectively. The first observation is the fact that the linear griwth rate increases with ∆ ′ , as expected. The second observation is the good agreement with analytic theory showing that the magnetic drift does not play any role when the whole electron population is passing. This is evidenced in the figures by the two black curves lying on top of each other. As predicted analytically, inclusing trapped electrons destabilizes the tearing mode. Finally, the impact of the magnetic drift on the tearing instability depends on the electron temperature gradient.

The impact of the electro temperature gradient can be further analysed by calculating the growth rate for a given ∆ ′ > 0 as a function of R/L T . This is plotted in figure 6.5. It can be again observed that when only passing electrons exist in the system, the magnetic drift does not play any role, independently of the temperature gradient. Regarding the situation where trapped electrons are considered, an overall increased of the growth rate is evidenced. In addition, the temperature gradient tends to destabilize the tearing mode only when considering trapped particles, as predicted in [START_REF] Peter | Trapped electron modifications to tearing modes in the low collision frequency limit[END_REF]. Finally, it is clearly observed that the temperature 102 gradient combined with the magnetic drift further destabilizes the tearing mode, which goes in the direction of an interchange-like destabilization predicted by the fluid theory [START_REF] Nasr | Interchange destabilization of collisionless tearing modes by temperature gradient[END_REF].

Impact of the tearing mode on the ITG-dominated turbulence

The previous sections are based on the linear excitation of the tearing mode. The nonlinear theory of the tearing instability [139] predicts that the growth of the mode is algebraic in time, i.e. there is a transition between the linear growth and the nonlinear saturation during which the amplitude grows proportionally to t.

In nonlinear regime and with strong temperature gradients, not only is the tearing mode further destabilized until a saturated island is obtained, but also other modes such as the ITG or the ETG can be destabilized, entering into the nonlinear phase and providing a turbulent state characterized by an increased transport. When both tearing mode and turbulence exist, their mutual interaction cannot be neglected. The mutual interaction between turbulence and magnetic islands was analysed in slab geometry where, in particular, the impact of magnetic islands on linearly stable ITG modes was reported [START_REF] Hu | Nonlinear mutual destabilization of the tearing mode and ion temperature gradient mode[END_REF]. Here, we consider only linearly unstable ITG modes bringing the system towards a strong turbulent state. Very recent gyro-kinetic simulations have provided for the first time the self-consistent growth of magnetic islands in the presence of electromagnetic turbulence [START_REF] William A Hornsby | The non-linear evolution of the tearing mode in electromagnetic turbulence using gyrokinetic simulations[END_REF][START_REF] Wa Hornsby | On seed island generation and the nonlinear self-consistent interaction of the tearing mode with electromagnetic gyro-kinetic turbulence[END_REF]. Nevertheless, these simulations are extremely expensive and analysis of the interaction between islands and turbulence remains difficult when both evolve self-consistently. An alternative to this is the implementation of an externally imposed magnetic island. This simplification is based on the fact that the growth of the tearing mode is typically slow, of the order of the resistive time τ R , compared to the growth of the small scale perturbations like the ITG modes. This technique has been used in the literature to determine the electrostatic potential in rotating islands and analyse the adiabatic response of ions [142], to determine the impact of small rotating islands on the bootstrap current [143], to study the impact of a static island on ITG turbulence [144,145,146], to prove the existence of a vortex mode related to the presence of an island [147] and to analyse how the breaking of toroidal symmetry due to the presence of a nonrotating magnetic island induces a nonambipolar current and a plasma torque [START_REF] Waltz | Gyrokinetic simulations with external resonant magnetic perturbations: Island torque and nonambipolar transport with plasma rotation[END_REF].

Here we focus on the impact of small and intermediate rotating magnetic islands on electrostatic turbulence, performing a scan on both the island width and the rotation frequency. This analysis is useful to understand what is observed in self-consistent simulations where the rotation frequency and the island width do evolve, but also to shed some light on experimental observations where the modification of the transport in the presence of a magnetic island might have an impact on the stability of the island itself and therefore on the measured island width. This approach where the magnetic island is externally imposed is justified by the time scale separation between ITG and tearing modes.

For the quantitative analysis of the impact on turbulence, we use from the electrostatic potential only the modes N ζ ≥ 3, since the mode N ζ = 1 represents the zonal flow and the mode N ζ = 2 represents the island mode. The filtered electrostatic potential is called δφ (x, t). We then define the turbulent intensity for

I turb (x) = 1 T T 0 |δφ (x, t) | 2 dt 1/2 (6.16)
where T is a sufficiently large number representing the length of the nonlinear saturated phase. Typically, we use T ≈ 800. In addition, we have a reference value of turbulent intensity in the absence of magnetic island, I turb,ref , calculated in the same way as in the simulations with magnetic island. The relative modification of turbulent intensity for each position x is defined as

δI turb (x) = I turb (x) -I turb,ref (x) I turb,ref (x) (6.17)
so that δI turb (x) < 0 (resp. > 0) indicates a decrease (resp. increase) of turbulent intensity. Finally, since we are interested in the impact of the island, we use a numerical filter to separate the region inside the separatrix from the region outside. This filter is based on the flux-surface label Ω = 2x 2 W 2 -cos π ζmax ζ , where Ω = 1 defines the separatrix and Ω > 1 represents the region outside the island, characterized by open flux surfaces. The filter is defined as follows

δ region (x) = H (Ω region -Ω) (6.18)
where H is the Heaviside distribution and Ω region corresponds to the flux-surface label which limits the region we are interested in. For example, if we are interested in analysing the turbulent intensity inside the whole island, Ω region = 1. If we are interested in the turbulent intensity only around the O-point, we have to define an arbitrary value to limit the region of the so-called O-point, i.e. it does not make any sense to analyse what is happening only at one precise location. In our case, we have decided to define the O-point as the region inside the Figure 6.6 illustrates an example of a ITG simulation, in particular the one for W = 6 rotating at the ion diamagnetic frequency. The quantity we plot is the turbulent intensity I turb (x). The region around the O-point used for the numerical filter is represented by the area inside the yellow line. The separatrix is given by the black line.

We apply the technique described to every simulation of the scan in island width and rotation frequency and we plot the modification of the turbulent intensity inside the island.

In figure 6.7 we give the relative modification of the turbulent intensity with respect to the reference simulations as a function of the island width and rotation frequency. The ratio ω isl /ω * i > 0 (resp. ω isl /ω * i < 0) corresponds to the ion (resp. electron) diamagnetic direction. This is explicitly indicated in the figure for the sake of clarity. Some observations can be made. First, the presence of an island tends to stabilize turbulence around the O-point. This general result is in agreement with previous gyrokinetic results [146,144,145], where the island was large enough to lead to a rather significant flattening of radial profiles and in addition it was a static island, i.e. there was no rotation frequency. The reduction of turbulence was due to the flattening of radial profiles leading to a reduction of the turbulence linear drive. However, in our case, we consider smaller islands and the flattening is not necessarily as significant as in the case of larger islands (see next section for further details). Second, we observe that islands rotating in CHAPTER 6. ELECTRO-MAGNETIC REGIME the ion diamagnetic direction are more efficient to reduce ITG turbulence than islands rotating in the electron diamagnetic direction. Third, very small islands (W ≈ 2) seem to impact turbulence in a less significant way than larger islands, characterized by W ≥ 18. Finally, the impact on turbulence due to intermediate islands (W ≈ 6 -9) depends more strongly on the island rotation frequency than the impact of larger islands. We observe that the width W = 6 -9 seems to be the optimal range for the rotation dependence. Beyond those island widths, the dependence on the island rotation frequency tends to disappear and the island seems to reduce the turbulence independently of the direction of the rotation.

Profile flattening in the presence of tearing modes and ITG-dominated turbulence

The stability of the magnetic island strongly depends on the non-inductive parallel current. One of the components of this current is the bootstrap current [149,150,151], which comes from the density and temperature gradients in the presence of trapped particles and low collisionality. A decrease of the bootstrap current within the magnetic island due to the flattening of the radial profiles produces a current perturbation with the same helicity as the magnetic island, exhibiting a destabilizing effect. The impact on the bootstrap current depends on the island width and rotation frequency. Nevertheless, the natural rotation frequency of an island is poorly known in experiments and simulations. In addition, the island evolves until a saturated width is reached. This is the reason why it is essential to predict the level of bootstrap current for different island widths and rotation frequencies. In Ref. [63] it is predicted that the contribution of the density gradient is about 65%, whereas the contributions of electron and ion temperatures are 20% and 15%, respectively. This implies that the modification of the density profile would have a larger impact on the bootstrap current than the modification of the temperature profile [?]. Predicting the level of density flattening is therefore essential to analyse the impact of an island on the bootstrap current, which in turn affects the stability of the island itself. In the absence of turbulent transport, it has been observed that small magnetic islands rotating in the ion diamagnetic direction flatten the density due to the adiabatic response of trapped ions, reducing in this way the stabilizing effect of the bootstrap current [143]. We quantify the impact on the density and temperature profiles by performing a time average over a sufficiently long nonlinear saturated phase. This aims at removing all fluctuating parts in order to have only equilibrium profiles. Therefore, the quantities that we use in this subsection are the time-averaged density and electron and ion pressures, defined respectively as In figure 6.8 we plot the normalized density gradient averaged around the Opoint by means of the same numerical filter. We observe that the island w = 6 rotating in the ion diamagnetic direction is able to flatten the density profile and this flattening level is reduced when increasing the island width. This is due to the adiabatic response of trapped ions [143,142]. The fact that the ion diamagnetic response dominates the flattening might be due to the reduction of turbulence when the island is rotating in the ion diamagnetic direction.

n i T (x) = 1 T T 0 n i (x, t) dt (
In figure 6.9a we plot the quantity R/L T,i for the analysed magnetic islands as a function of the island rotation frequency and the island width. Also in this figure, for reference, we give the value of R/L T,i for a static island W = 12 when the turbulent modes are filtered out in the simulation. This means that this simulation has been run with only a maximum value of κ ζ ρ ref = 0.2 and N ζ = 5. Note that the necessity of using a sufficiently high number of binormal modes to correctly describe the separatrix was pointed out in a separate publication [START_REF] Zarzoso | Verification of a magnetic island in gyro-kinetics by comparison with analytic theory[END_REF]. Therefore, the value that we give with the blue dashed line in figure 6.9a must be considered as a slightly underestimated reference value, since the flattening cannot be correctly described due to the artificially wide separatrix in that simu- lation. Some observations can be made from this figure. First, for islands W = 9 and W = 12, the ion temperature profile tends to be only slightly flattened, i.e. R/L T,i remains close to the equilibrium value. Whether this absence of flattening is due or not to the presence of turbulence can be elucidated with the reference simulation represented by the blue dashed line. We observe that the reduction of turbulent transport leads to a significant decrease of R/L T . Therefore, the island W = 12 would be in principle large enough to provide a flattening due to the parallel streaming, but the level of turbulence is strong enough to keep the gradient well beyond the linear threshold. Note that when the island becomes larger (W = 18) the degree of flattening is significantly increased. Second, for the island W = 6, the impact on the ion temperature profile is different. In particular, we observe an increase of ion temperature gradient when the island rotates in the ion diamagnetic direction. This increase is due to the fact that the ion pressure is almost unperturbed. However, the ion temperature gradient is determined from the pressure and density profiles. Since the density is flattened, the ion temperature must be steeper. This means that the linear drive is not suppressed and therefore the reduction of turbulence for the intermediate island rotating in the ion diamagnetic direction is likely due to a mechanism different from the standard suppression of linear drive. Finally, when the island is large enough (W ≈ 18), the observed flattening results in the flatter ion temperature profiles, leading to a significant reduction of the linear drive and therefore providing a stabilization of ITG turbulence [144,145], though the same mechanism as in the island with W = 6 might co-exist.

For completeness, we also show in figure 6.9b, the normalized electron tem- perature gradient. It can be observed that this quantity follows a trend understood in terms of width of the island, independently of the rotation frequency. This is expected from the time scales that we gave at the beginning of this section. Let us note that in particular, for W = 12, the flattening is almost the same as the one obtained without turbulence and when the island is increased no significant extra-flattening is observed. Owing to the flattening of electron density and temperature profiles, one could think that the reduction of turbulence observed with an intermediate island rotating in the ion diamagnetic direction might be due to a linear stabilization of sub-dominant TEM turbulence. To clarify the validity of this conjecture, we have run a simulation without magnetic island equivalent to the ITG simulation for the case W = 6, ω isl = ω * i using the normalized profile gradients given by the simulation with magnetic island, i.e. R/L T,i ≈ 9, R/L T,e ≈ 3 and R/L n ≈ -1. A quantification of the level of fluctuations with respect to the reference simulation gives an increase of ∼ 40%. This is observed in figure 6.10, where we plot the time evolution of the turbulent intensity of the CBC simulation (solid line) and the turbulent intensity of the simulation without the magnetic island, but using the background gradients given by the simulation with magnetic island (dashed line). Therefore, the stabilization of turbulence for that island must be due to nonlinear effects.

Possible mechanism of turbulence reduction

From the previous sections we can conclude that the impact of magnetic islands on turbulence is more subtle than the reduction of linear drive. The suppression of the linear drive seems to be the dominating mechanism for turbulence reduction when the island is large enough. However, for intermediate magnetic islands, another mechanism should be considered. We focus here on a possible nonlinear mechanism which might be responsible for the reduction of fluctua-CHAPTER 6. ELECTRO-MAGNETIC REGIME tions inside and around the island. This nonlinear mechanism can be related to the increased of a radial electric shear around the separatrix when the island is rotating with respect to the plasma [?]. The transfer of energy is expected to be more efficient when island and turbulence structures are co-rotating because the time of interaction between large-scale (zonal flows and island) and small-scale (turbulence) structures is maximized. In addition to the time interaction, one must consider also the scale lengths. In particular when the shear of the radial electric field exhibits scale lengths close to the decorrelation length of turbulence, the transfer of energy can lead to a reduction of turbulence. To analyse the radial correlation length of turbulence we define the autocorrelation function as For the parameters we are considering in this paper, the turbulence in the absence of magnetic island is characterized by a radial autocorrelation function given in figure 6.11. The normalized integral of this autocorrelation function is used to define the so-called radial correlation length λ c = 1 2 rmax r min R (∆r) d∆r (6.25) This radial correlation length can be linked to the half width of the radial autocorrelation function evaluated at half maximum, which is noted W 1/2 in figure 6.11, as W 1/2 = λ c log 2, assuming an exponential decrease. For the ITG parameters used in this paper the radial correlation length as defined by equation 6.25 is λ c = 5.2 in ion Larmor radius units. If we assume an exponential decrease, from figure 6.11 we obtain W 1/2 = 4.2, which is in agreement with global flux-driven gyrokinetic simulations of ITG turbulence [START_REF] Sarazin | Predictions on heat transport and plasma rotation from global gyrokinetic simulations[END_REF], and the radial correlation length assuming an exponential decrease is λ c,exp = 6.09. The radial auto-correlation function is characterized by a long tail. The presence of this tail was reported in a previous publication [?], where the existence of two scales was suggested. The first scale is the mean microscopic eddy size, which is of the order of a few ion Larmor radii. The second scale is due to the mesoscale eddies.

An interaction via an energy transfer between large and small binormal scales is possible within the macro and meso radial scales. As explained earlier in this section the transfer is expected to be more efficient when both island structure and eddies are co-rotating. In that case, we expect a dramatic reduction of turbulence: stronger at the microscropic eddy size scale, where the autocorrelation is statistically larger, and less pronounced at a scale beyond the microeddies up to some tens of ion Larmor radii. In the case of ITG turbulence, the eddies rotate in the ion diamagnetic direction. Therefore, for a magnetic island as wide as the statistically relevant eddies, the impact on turbulence depends on the rotation frequency of the island and, in particular, turbulence is expected to be reduced when the island rotates in the ion diamagnetic direction, due to a maximization of the interaction time between turbulence and large-scale structures. This implies that the strongest dependence on the island rotation frequency occurs when W ≈ λ c . This condition is closer to the simulation W = 6 in the scan presented in this paper. In that case, energy is transferred from ITG turbulence to island when the island rotates in the ion diamagnetic frequency. If a transfer of energy exists from turbulence to large-scale structures, this must lead to an increase of energy of the modes k ζ = 0 and k ζ,isl . This transfer of energy comes from the second derivative of the electrostatic potential, which can be increased in the presence of an island around the separatrix. In figure 6.12 we plot for the island width W = 6 the radial profile of the second derivative of the time-averaged electrostatic potential through the O-point. This is done for two rotation frequencies: ion (solid line) and electron (dashed line) diamagnetic frequencies. The island separatrix is highlighted by two vertical lines. It can be observed that the shear is clearly increased when the island is rotating in the ion direction.

The mechanism of nonlinear coupling between island and turbulent modes is rarely isolated and in a general scenario where the island exhibits an intermediate size, two mechanisms co-exist: (1) the reduction of the linear drive due to either the parallel streaming of particles or the adiabatic response of trapped ions and (2) the transfer of energy from the turbulent modes to the island structure at two radial scales related to the nature of the turbulence. .12: Radial profile through the O-point of the shear of the radial electric field for two island rotation directions: electron diamagnetic direction (dashed line) and ion diamagnetic direction (solid line). The case of an island width W = 6 is considered.

A few words on the transport of particles in the electro-magnetic regime

As we saw in chapter 4, one single electrostatic wave can lead to a chaotic transport of particles, characterised by anomalous diffusion. This transport was interpreted as the result of a stochastization of the separatrix between magnetically passing and trapping particles in the presence of an oscillating radial electric field. This phenomenon could be used to understand the losses induced by EGAMs in current fusion devices. Within the framework of the present chapter, one can be naturally tempted to extend those results to the case with magnetic perturbations. In this chapter we are interested in particular in tearing modes, which are characterised by low frequency and low poloidal and toroidal wave numbers. As in the case of a single electrostatic wave, one could intuitively believe that a single helicity mode (m, n) has an impact on particle transport reduced to the mere radial motion with respect to the unperturbed magnetic surfaces. This is basically the main idea behind the presence of a magnetic island: the particle follows the magnetic field lines, given by B = B eq + δB, where δB is oriented mainly in the radial direction for a tearing mode. The radial component of the total magnetic field makes particles explore a radial region whose extension is given by the island width, therefore proportional to the square root of the magnetic perturbation. In the presence of a single helicity mode, this implies that the maximum radial excursion of particles would be given by the island width. However, coherent MHD instabilities can have an impact on the confinement of energetic particles [START_REF] Heidbrink | The behaviour of fast ions in tokamak experiments[END_REF]. If MHD modes are driven by energetic particles, they naturally resonate, stay in phase and particles can be transported at large scale. For low frequency (like tearing modes) or almost static MHD instabilities, resonances might not play a significant role and transport could come from the fact that the motion of particles (characterised by magnetic drifts) can couple to the single helicity perturbation. This can lead to the presence of multiple resonances in phase space and eventually to a chaotic transport [START_REF] Harry E Mynick | Transport of energetic ions by low-n magnetic perturbations[END_REF][START_REF] He Mynick | Stochastic transport of mev ions by low-n magnetic perturbations[END_REF]. As discussed in [START_REF] Em Carolipio | Simulations of beam ion transport during tearing modes in the diii-d tokamak[END_REF], the transport induced by stochastization in phase space is intrinsically different from parallel transport along ergodic magnetic field lines. Indeed, the former one comes from the finite Larmor radius effects (FLR), which introduces orbit drifts. Therefore, FLR tend to increase chaotic transport in phase space. On the other hand, parallel transport in ergodic magnetic field lines is reduced due to gyro-average over the perturbations. The transport of particles in the presence of low-frequency MHD instabilities is especially important because of experimental evidences in ASDEX Upgrade [START_REF] García-Mu Ñoz | Ntm induced fast ion losses in asdex upgrade[END_REF] and DIII-D [START_REF] Em Carolipio | Simulations of beam ion transport during tearing modes in the diii-d tokamak[END_REF] tokamaks, where NBI born energetic ions are transported and subsequently lost at the frequency of the mode [START_REF] García-Mu Ñoz | Ntm induced fast ion losses in asdex upgrade[END_REF], leading to a ∼ 65% reduction of the beam current drive efficiency [START_REF] Em Carolipio | Simulations of beam ion transport during tearing modes in the diii-d tokamak[END_REF]. A generalization of the code used in chapter 4 to integrate the guiding-centre equations can be developed in the electro-magnetic regime to analyse the transport of guiding-centres in the presence of a tearing mode. Let us assume a purely magnetic regime (φ = 0) with a parallel vector potential given by [START_REF] García-Mu Ñoz | Ntm induced fast ion losses in asdex upgrade[END_REF] A = A ,0 r m (ψ -ψ w ) cos (mθ + nϕ + ωt) (6.26)

where A ,0 is the amplitude and ψ w is the poloidal flux at the outer boundary. This parallel potential is used to calculate the total magnetic field as B = B eq + ∇A × b (6.27)

Note that in the presence of the parallel vector potential 6.26, neither the kinetic energy nor the toroidal canonical momentum are conserved. In [START_REF] Em Carolipio | Simulations of beam ion transport during tearing modes in the diii-d tokamak[END_REF][START_REF] García-Mu Ñoz | Ntm induced fast ion losses in asdex upgrade[END_REF], a static perturbation was used, in which case the kinetic energy and the magnetic moment are conserved, but not P ϕ . Knowing that in a three-degrees-of-freedom system with two motion invariants the trajectories are not integrable in the sense of Poincaré, we can already presume that a transport will occur, even in the presence of a single helicity perturbation. Including a time dependence reduces the number of motion invariants to the magnetic moment. Figure 6.13a shows the Poincaré map of the magnetic field perturbed by a parallel potential whose amplitude has been chosen so that the width of the island is ∼ 20 thermal Larmor radii on the low field side, which is equivalent to ∼ 20 cm for a 3 keV plasma of Deuterium, which is close to the island widths observed in experiments. Figure 6.13a shows the presence of m = 2 and m = 3 islands, which appear at the positions q = 2 and q = 3 (cf. q-profile shown in figure 6.13b), as expected. One would expect that the main island (m = 2) would trapped particles inside. However, as shown in chapter 2 several resonances at different n 2 = 0 values exsit in phase space, what can lead potentially to stochastization of the orbits. An example of these orbits is shown in figure 6.13a, where is depicted the projection onto the poloidal cross-section of the trajectory of an initially counter-passing particle (red line). For comparison, the blue line represents its unperturbed trajectory (in the absence of any perturbation). Also, for the sake of completeness, it is given the Poincar ḿap of the equilibrium and perturbed magnetic field, as in figure 6.13a. To show how the particle is able to cross different n 2 resonances, its motion in phase-space is plotted in figure 6.14b by the magenta line. The other lines in that figure where described in chapter 2.

Chapter 7

Conclusions and future directions

The previous chapters have shown to what extent the kinetic description of particles in the presence of instabilities is essential to understand their confinement and losses. We have also shown how the mutual interaction between energetic particle driven modes and turbulence cannot be neglected, for a nonlinear coupling can exist between the interacting modes. The work presented in this dissertation focuses on the special class of electrostatic modes called EGAMs, on the electrostatic ITG-dominated turbulence and the purely magnetic tearing modes. Restricting the analysis to these modes is certainly easier than analysing the complete electro-magnetic problem in fully kinetic simulations with experimental equilibria in the presence of instabilities covering wider ranges of spatio-temporal scales. Nevertheless, the choice of presenting the work this way is not arbitrary. Indeed, the EGAMs are taken as the paradigm of energetic particle driven mode. Because it is mainly electrostatic and axi-symmetric, it must be considered as the one of the simplest instabilities in tokamaks and without any doubt the first step when trying to understand the fundamental mechanisms behind the excitation of energetic particle driven modes. And yet, these simple modes reveal in this dissertation mechanisms that were unimaginable so far, like their negative impact on turbulence due to a nonlinear coupling of modes or their impact on particle trajectories leading to a chaotic transport, characterised by anomalous diffusion. Similarly, tearing modes constitute one of the most studied magnetic cases in nuclear fusion using fluid approaches, but not much work has been done within a gyro-kinetic framework and even less using global gyro-kinetic simulations, where only few works exist in the literature. We have shown how a gyro-kinetic description of the tearing mode can lead to important consequences for the tearing stability, linked to the trapping of particles in the equilibrium magnetic field. Also the impact of the kinetic effects on the interaction between magnetic islands and turbulence and eventual impact on radial profiles has been evidenced. We have taken advantage of the apparent simplicity of these modes to offer a rigorous theoretical framework based on a Hamiltonian approach built upon the existence of action-angle variables. This is essential to understand the main physical mechanisms responsible for the results shown in this dissertation. Of course, once the framework has been set up, it is natural to extent this work to other regimes and scenarios, which constitutes the current and future directions of my career as researcher and supervisor. Three directions are envisioned.

• Electro-magnetic interaction between turbulence and energetic particles in the Alfvén regime.

The results I have found in the electrostatic limit should now be applied to the electro-magnetic regime. It is important to realise that the alpha particles born from nuclear reactions will excite non axi-symmetric in the Alfvén range of frequencies. On the route towards this electromagnetic regime, it is essential to understand the physics of these instabilities excited by energetic particles present in current fusion devices and how they interact with turbulence. Recent research on this topic has been published [START_REF] Citrin | Nonlinear stabilization of tokamak microturbulence by fast ions[END_REF][START_REF] Garcia | Key impact of finite-beta and fast ions in core and edge tokamak regions for the transition to advanced scenarios[END_REF], but much work remains to be done, especially to understand the impact of energetic ions on a class of turbulence different from that driven by ITG modes. This is the purpose of a PhD I supervise, where experimental profiles from the JET and JT-60U tokamaks and the gyro-kinetic code GENE are used to determine the impact of energetic ions on trapped-electron driven turbulence. This is important because it is a class of turbulence which is expected to dominate the scenarios of the JT-60U tokamak, built within the framework of the broader approach [START_REF] Kikuchi | Overview of modification of jt-60u for the satellite tokamak program as one of the broader approach projects and national program[END_REF][START_REF] Matsuda | The eu/ja broader approach activities[END_REF][START_REF] Tsunematsu | Broader approach to fusion energy[END_REF].

• Interaction between energetic particles and energetic particles.

This project is currently being developed within the framework of a EU-ROfusion project. As explained in the previous chapter, the interaction between energetic particles and tearing modes is still far from being understood. The goal of this project is to assess at the same time the impact of tearing modes on the transport of energetic particles and the impact of the energetic particles on the stability of the tearing mode. The selected approach is a reduced model where the guiding-centre tracking code that I have developed will be coupled to the equations governing the evolution of the nonlinear stability of the tearing mode, improving the approach used in [START_REF] Chris | Suppression of the tearing mode by energetic ions[END_REF].

• Investigating new numerical techniques to model fusion plasmas. Using analytical calculations to describe the behaviour of complex systems like fusion plasmas reveals sometimes a challenging task. A possible solution consists in solving numerically the equations by means of massively parallelised codes. Nonetheless, this approach remains nowadays quite demanding in terms of required CPU time, even though the performance of supercomputers has considerably increased in the past ten years. Also, the accelerated production of high-dimensionality data forces us to find new ways to analyse the data. In essence, having huge amounts of experimental and numerical data does not necessarily mean increasing our knowledge.

We have reached a critical point in Plasma Physics and Nuclear Fusion research where a new way of thinking and analyzing the data is mandatory to make further progress. This is the reason why part of my future research

1 . 1

 11 Graphical representation of a system composed of N=20 entities. Their mutual interactions are represented by red arrows. . . . . . . 1.2 (Left) Speed of each car at the last instant of the simulation. (Middle) Speed of the first car as a function of time. (Right) Speed of all the cars as a function of time. . . . . . . . . . . . . . . . . . . . . . . 1.3 Number of normal (solid blue line) and cancerous (dashed red line) for mutation probabilities p mut = 10 -4 (left panel) and p mut = 10 -2 (right panel). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Tokamak geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Examples of passing (left) and trapped (right) particle trajectories. The top panels represent the 3D views, whereas the bottom panels represent the projection onto the poloidal cross-section. . . . . . . . 2.3 Radial profile of the safety factor as a function of the minor radius (left) and the normalised poloidal flux (right). . . . . . . . . . . . . . 2.4 Resonances in the ψ, E space for two values of magnetic moment and n 3 = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Resonances in the ψ, E space for two values of magnetic moment and n 3 = -1. The colors for the values of n 2 are the same as in figure 2.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Evolution of the zeros of the dispersion relation with increasing energetic particle concentration with q = 2. . . . . . . . . . . . . . . . . . . . . . .

OF FIGURES 3 . 6 48 4. 1 OF FIGURES 5 . 5 80 5. 6 ( 81 5. 7 ( 84 5. 9 (

 3648155806817849 . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Comparison of two different terms to the total energy transfer from energetic particles to the mode. The fact that the main contribution comes from the curvature term indicates that the mode has indeed a geodesic nature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Time evolution of the change of the kinetic energy of the different species at q = 2 at the threshold (left) and beyond the threshold (right). . . . . . . xi LIST For the simulation n EP = 0.02n eq . (Left panel) Time trace of the flux-surface averaged electrostatic potential evaluated at the mid radial position. (Right panel) Fourier transform of the signal plotted on the left panel over the entire time window. . . . . . . . . . . 45 3.7 For the simulation with n EP = 0.02n eq , plotted is the perturbed distribution function as defined by expression (3.10) for thermal (thick solid blue line) and energetic (thick dashed red line). The vertical thin lines indicate the different resonances (l, m), explained in the text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.8 Snapshots of the perturbed distribution function for the simulation n EP = 0.02n eq for thermal (left panel) and energetic (right panel) particles as a function of the parallel velocity and the poloidal angle, for µ ≈ 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.9 Snapshots of the perturbed distribution function for the simulation n EP = 0.02n eq for thermal (left panel) and energetic (right panel) particles as a function of the parallel velocity and the poloidal angle, for µ = 6, normalised to B 0 /T 0 . . . . . . . . . . . . . . . . . . . . Intersection of the surfaces E = cst and P ϕ = cst and trajectory of a particle in the absence (left panel) and in the presence (right panel) of an EGAM. The trajectory clearly departs for the E = cst surfaces, whereas it remains on the P ϕ = cst in the case of an axisymmetric time-dependent electrostatic potential. . . . . . . . . . . 52 4.2 Losses of particles in the presence of EGAMs excited by a n EP = 0.02n eq fraction of energetic particles. The left (resp. right) panel represents the losses for particles initialised at r/a = 0.2 (resp. r/a = 0.6). The black line represents the trapping region. Blue (resp. red) points represent prompt (resp. EGAM) losses. . . . . . . 53 4.3 Perturbed (top) and unperturbed (bottom) trajectories of test gyrocentres. The middle and right panels on the top present trajectories of three particles that are lost at different times during the simulation. The left panel on the top presents only the projection onto the poloidal cross section for one of the particles (dotted blue line). For both perturbed and unperturbed trajectories, the boundary domain in real space, corresponding to a circular flux-surface, is given in the left panel and represented by a dotted black line. . . 54 4.4 (Left) Losses of particles for the zoomed area shown in figure 4.2. (Middle) Time trace of the losses of particles (black line) together with the time-trace of the electrostatic potential (blue line). (Right) Fourier transform of the time-trace of the losses of particles, exhibiting dominant components at the EGAM frequency and its harmonics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.5 Position of particles in (v , θ) space at three different instants: before (left), just after (middle) and well after (right) the nonlinear saturation of the EGAM. . . . . . . . . . . . . . . . . . . . . . . . . . 56 xii LIST OF FIGURES 4.6 (Left) Trajectories of two counter-passing particles: one deeply counterpassing with (Λ = 0.4, E = 25E th ), represented by the almost circular projection, and one barely counter-passing with (Λ = 0.8, E = 43E th ). (Right) Ensemble of all the possible trajectories with energies within the range 25E th ≤ E ≤ 43E th and pitch angle 0.4 ≤ Λ ≤ 0.8. . . . . . 4.7 (Left) Fraction of lost particles as a function of their initial energy and pitch angle. The contours P ϕ = cst are overlaid. (Right) Poincaré map of the trajectory of a lost particles with initial Λ = 0.54 and E = 31.5E th . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Inverse of the exit time as a function of the initial parallel velocity and magnetic moment. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Probability distribution function (PDF) of the exit time for counterpassing EP initialised at r = 0.4a, θ = 0, ϕ = 0, with E = 20E th and µB 0 /T i = 14. The dashed line shows an algebraic fit, with µ e = 2.5. 4.10 Poincaré map of unperturbed trajectories (black lines) and particles initialised on the separatrix in the presence of an EGAM (red region). The direction of rotation of particles in the inner and outer regions of the tokamak is represented by dashed blue lines in the left panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 (Left) Poloidal displacement of passive tracers showing a spreading in the poloidal direction. (Right) Time dependence of the variance of the poloidal displacement. . . . . . . . . . . . . . . . . . . . 4.12 (Left) Poloidal displacement of two passive tracers, showing the existence of positive and negative flights. (Right) Poloidal displacement assuming a standard random walk. . . . . . . . . . . . . 4.13 PDF of negative (left) and positive (right) flight events of duration t. 4.14 Rescaled PDF of poloidal displacements at different times: ω c t = 9 • 10 6 (dashed magenta), ω c t = 9.5 • 10 6 (dotted red) and ω c t = 10 7 (solid black). The dotted-dashed grey line corresponds to a Gaussian PDF. The insets represent the log-log plots of the tails, showing the asymmetric algebraic decays. . . . . . . . . . . . . . . . 5.1 Diagram showing the marginal stability line separating the stable (below the line) and unstable (above the line) regions. . . . . . . . . 5.2 (Left) Schematic view of the convergence domain of the HFA. (Right) Wave-particle energy exchange (colour) together with the convergence domains of HFA and IHFA. For the sake of completeness, the boundary between stabilising and destabilising regions is given by the thin dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Real frequency (left) and growth rate (right) using different methods to write to the dispersion relation. . . . . . . . . . . . . . . . . . 5.4 Radial envelope of thermal source (dashed blue line), energetic particle source (red line) and buffer region (shadowed area). The amplitudes are arbitrary. . . . . . . . . . . . . . . . . . . . . . . . . . xiii LIST Schematic view of the thermal (left) and energetic particle (right) sources as a function of the velocity and parallel velocity, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Left) Radial profiles of ion guiding-centre temperature (blue line) and density (red line). The shadowed region represents the radial envelope of the buffer region in arbitrary units. (Right) Radial profile of the turbulent intensity as calculated by expression 5.24, averaged over ω c ∆t ≈ 3 • 10 4 . . . . . . . . . . . . . . . . . . . . . . . . . Left) Radial dependence of the turbulent intensity as given by expression 5.24 for the three simulations used in this chapter. (Right) Time evolution of the temperature gradient after switching on the energetic particle source, averaged over the radial region 50ρ i < r < 100ρ i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.8 (Left) Time evolution of the flux-surface averaged electrostatic potential for the three simulations. (Right) For the simulation with S EP /S th = 0.625, time Fourier transform of the electrostatic potential after filtering out the time averaged signal. . . . . . . . . . . . . Left) Time evolution of the slope of the equilibrium distribution function (thin lines) evaluated at v = qω EGAM R 0 . Time-averages over an EGAM cycle are plotted by thick lines of the same colour. The shadowed region represents the time evolution of the heat flux due to the E × B drift velocity. (Right) Equilibrium distribution functions without energetic particles (thin dashed line) and with energetic particles just before the excitation of EGAMs using S EP /S th = 0.437 (thick dashed line) and S EP /S th = 0.625 (solid line). 85 5.10 (Left) Time evolution of the flux-surface averaged electrostatic potential with (solid red line) and without (dashed magenta line). (Right) Comparison of the slope of the equilibrium distribution function at the resonant velocity with (red line) and without (magenta line) turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.11 Radial dependence of the m = 1 (thick lines) and m = 2 (thin lines) poloidal components of the electrostatic potential oscillating at the EGAM frequency. The sine and cosine components are evidenced by solid and dashed lines, respectively. . . . . . . . . . . . . . . . . . 88

  . . . . . . . . . . . . . . . . . . . 88 xiv LIST OF FIGURES 5.13 (Left) Radial dependence of the time-averaged E × B heat flux for the simulations of this paper. For comparison, the heat flux in the absence of EP is also plotted and represented by the thin dashed black line. (Right) Time Fourier transform of the E × B heat flux evaluated at the mid radial position for the simulation with S EP /S th = 0.625. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 Wavelet bispectrum for the four phases indicated in the main text of the paper and on top of each of the subfigures. For clarity and only in the last two figures, the EGAM frequency has been indicated by horizontal yellow dashed lines. . . . . . . . . . . . . . . . . 6.1 Solution of equation 6.5 for different values of ν, q a = 3.5 and ǫ = 1/3.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Ratio between growth rates with and without trapping electrons in the absence of magnetic drift and comparison with global gyrokinetic simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Energy exchange between the mode and particles with flat profiles and ω D = 0 in the presence of only passing particles (top panel) and when the trapping effects are taken into account (bottom panel).102 6.4 Growth rate of the tearing mode as a function of the stability parameter with flat profiles (left) and with electron temperature gradient (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Growth rate of the tearing instability as a function of the electron temperature gradient. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 For the simulation W = 6, ω ils /ω * i = 1, plotted is the turbulent intensity I turb (x) as a function of radial and binormal positions.

1 The multi-scale problem in complex systems Contents 1 . 1 2 1. 2 4 1. 2 . 1 4 1. 2 . 2 6 1. 3 8 1. 3 . 1 9 1. 3 . 2 9 1. 3 . 3 10 1. 4 11 1 1 . 1 .

 11122421422638319329331041111 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.14 (Left) Trajectory of a lost particle in the poloidal cross-section together with the Poincaré map of the magnetic field perturbed by a tearing-like instability. Blue (resp. red) line represents the unperturbed (resp. perturbed) trajectory. (Right) Perturbed trajectory in the space (ψ, E) diagram as shown in chapter 2 (right). . . . . . . . 115 xvi Chapter Complex systems: from the smallest entity to the global behaviour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Some examples of complex systems . . . . . . . . . . . . . . . . . Modelling the traffic in a round-about or how to get a jam without any bottleneck . . . . . . . . . . . . . . . . . . . . . Modelling a tumour growth . . . . . . . . . . . . . . . . . . Nuclear fusion plasmas as complex systems with multi-scale instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ion-temperature-gradient (ITG) driven instability . . . . . Energetic particle driven modes (EPMs) . . . . . . . . . . . The tearing mode (TM) instability . . . . . . . . . . . . . . The multi-scale problem in my research . . . . . . . . . . . . . . COMPLEX SYSTEMS: FROM THE SMALLEST ENTITY TO THE GLOBAL BEHAVIOUR

Figure 1 . 1 :

 11 Figure 1.1: Graphical representation of a system composed of N=20 entities. Their mutual interactions are represented by red arrows.

Figure 1 . 2 :

 12 Figure 1.2: (Left) Speed of each car at the last instant of the simulation. (Middle) Speed of the first car as a function of time. (Right) Speed of all the cars as a function of time.

Figure 1 . 3 :

 13 Figure 1.3: Number of normal (solid blue line) and cancerous (dashed red line) for mutation probabilities p mut = 10 -4 (left panel) and p mut = 10 -2 (right panel).

Contents 2 . 1 2 . 2

 2122 The tokamak and the toroidal geometry . . . . . . . . . . . . . . Particle trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Tokamak geometry

Figure 2 . 2 :

 22 Figure 2.2: Examples of passing (left) and trapped (right) particle trajectories. The top panels represent the 3D views, whereas the bottom panels represent the projection onto the poloidal cross-section.

32 )Figure 2 . 3 :

 3223 Figure 2.3: Radial profile of the safety factor as a function of the minor radius (left) and the normalised poloidal flux (right).

Figure 2 . 4 :

 24 Figure 2.4: Resonances in the ψ, E space for two values of magnetic moment and n 3 = 0.

Figure 2 . 5 :

 25 Figure 2.5: Resonances in the ψ, E space for two values of magnetic moment and n 3 = -1. The colors for the values of n 2 are the same as in figure 2.4.

CHAPTER 3 . 1 ReFigure 3 . 1 :

 3131 Figure 3.1: Evolution of the zeros of the dispersion relation with increasing energetic particle concentration with q = 2.

Figure 3 . 2 :

 32 Figure 3.2: Evolution of the zeros of the dispersion relation with increasing energetic particle concentration with q = 3.

Figure 3 . 3 :

 33 Figure 3.3: Frequency (left) and growth rate (right) of EGAM as predicted in[1] (lines) and as measured from GYSELA simulations (symbols). The scan is performed using the same species for the thermal population and modifying the mass and the charge for the energetic particle population.

Figure 3 . 4 :

 34 Figure 3.4: Comparison of two different terms to the total energy transfer from energetic particles to the mode. The fact that the main contribution comes from the curvature term indicates that the mode has indeed a geodesic nature.

Figure 3 . 5 :

 35 Figure 3.5: Time evolution of the change of the kinetic energy of the different species at q = 2 at the threshold (left) and beyond the threshold (right).

Figure 3 . 6 :

 36 Figure 3.6: For the simulation n EP = 0.02n eq . (Left panel) Time trace of the flux-surface averaged electrostatic potential evaluated at the mid radial position. (Right panel) Fourier transform of the signal plotted on the left panel over the entire time window.

Figure 3 . 7 :

 37 Figure 3.7: For the simulation with n EP = 0.02n eq , plotted is the perturbed distribution function as defined by expression (3.10) for thermal (thick solid blue line) and energetic (thick dashed red line). The vertical thin lines indicate the different resonances (l, m), explained in the text.

Figure 3 . 8 :

 38 Figure 3.8: Snapshots of the perturbed distribution function for the simulation n EP = 0.02n eq for thermal (left panel) and energetic (right panel) particles as a function of the parallel velocity and the poloidal angle, for µ ≈ 0.

Figure 3 . 9 :

 39 Figure 3.9: Snapshots of the perturbed distribution function for the simulation n EP = 0.02n eq for thermal (left panel) and energetic (right panel) particles as a function of the parallel velocity and the poloidal angle, for µ = 6, normalised to B 0 /T 0 .

Figure 4 . 1 :

 41 Figure 4.1: Intersection of the surfaces E = cst and P ϕ = cst and trajectory of a particle in the absence (left panel) and in the presence (right panel) of an EGAM. The trajectory clearly departs for the E = cst surfaces, whereas it remains on the P ϕ = cst in the case of an axi-symmetric time-dependent electrostatic potential.

52

 52

Figure 4 . 2 :

 42 Figure 4.2: Losses of particles in the presence of EGAMs excited by a n EP = 0.02n eq fraction of energetic particles. The left (resp. right) panel represents the losses for particles initialised at r/a = 0.2 (resp. r/a = 0.6). The black line represents the trapping region. Blue (resp. red) points represent prompt (resp. EGAM) losses.

4. 3 .

 3 Figure 4.3: Perturbed (top) and unperturbed (bottom) trajectories of test gyrocentres. The middle and right panels on the top present trajectories of three particles that are lost at different times during the simulation. The left panel on the top presents only the projection onto the poloidal cross section for one of the particles (dotted blue line). For both perturbed and unperturbed trajectories, the boundary domain in real space, corresponding to a circular flux-surface, is given in the left panel and represented by a dotted black line.

Figure 4 . 4 :

 44 Figure 4.4: (Left) Losses of particles for the zoomed area shown in figure 4.2. (Middle) Time trace of the losses of particles (black line) together with the timetrace of the electrostatic potential (blue line). (Right) Fourier transform of the time-trace of the losses of particles, exhibiting dominant components at the EGAM frequency and its harmonics.

Figure 4 . 5 :

 45 Figure 4.5: Position of particles in (v , θ) space at three different instants: before (left), just after (middle) and well after (right) the nonlinear saturation of the EGAM.

Figure 4 . 6 :

 46 Figure 4.6: (Left) Trajectories of two counter-passing particles: one deeply counter-passing with (Λ = 0.4, E = 25E th ), represented by the almost circular projection, and one barely counter-passing with (Λ = 0.8, E = 43E th ). (Right) Ensemble of all the possible trajectories with energies within the range 25E th ≤ E ≤ 43E th and pitch angle 0.4 ≤ Λ ≤ 0.8.

Figure 4 . 7 :

 47 Figure 4.7: (Left) Fraction of lost particles as a function of their initial energy and pitch angle. The contours P ϕ = cst are overlaid. (Right) Poincaré map of the trajectory of a lost particles with initial Λ = 0.54 and E = 31.5E th .

Figure 4 . 10 :

 410 Figure 4.10: Poincaré map of unperturbed trajectories (black lines) and particles initialised on the separatrix in the presence of an EGAM (red region). The direction of rotation of particles in the inner and outer regions of the tokamak is represented by dashed blue lines in the left panel.

Figure 4 .

 4 Figure 4.11: (Left) Poloidal displacement of passive tracers showing a spreading in the poloidal direction. (Right) Time dependence of the variance of the poloidal displacement.

Figure 4 .

 4 Figure 4.12: (Left) Poloidal displacement of two passive tracers, showing the existence of positive and negative flights. (Right) Poloidal displacement assuming a standard random walk.

PFigure 4 . 13 :

 413 Figure 4.13: PDF of negative (left) and positive (right) flight events of duration t.

Figure 4 . 14 :

 414 Figure 4.14: Rescaled PDF of poloidal displacements at different times: ω c t = 9 • 10 6 (dashed magenta), ω c t = 9.5 • 10 6 (dotted red) and ω c t = 10 7 (solid black). The dotted-dashed grey line corresponds to a Gaussian PDF. The insets represent the log-log plots of the tails, showing the asymmetric algebraic decays.

Contents 5 . 1 5 . 2 5 . 3 5 . 4 5 . 5 5 . 6 5 . 7 5 . 8 5 . 9

 515253545556575859 Destabilization of ion temperature gradient modes . . . . . . . . 5.1.1 Heuristic interpretation of the the marginal stability in the fluid limit . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 From an improved fluid approximation to the gyro-kinetic dispersion relation . . . . . . . . . . . . . . . . . . . . . . . The nonlinear coupling between modes . . . . . . . . . . . . . . Introducing energetic particles in a gyro-kinetic code . . . . . . Description of a typical initial turbulent state . . . . . . . . . . . Onset of the energetic particle source: initial reduction of turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Excitation and nonlinear saturation of EGAMs in the presence of turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Increase of turbulent transport in the presence of EGAMs . . . Radial structure of EGAM and modification of turbulence intensity and transport . . . . . . . . . . . . . . . . . . . . . . . . . . Nonlinear coupling between EGAMs and ITG modes . . . . . .

Figure 5 . 1 :

 51 Figure 5.1: Diagram showing the marginal stability line separating the stable (below the line) and unstable (above the line) regions.

Figure 5 .Figure 5 . 2 :

 552 Figure 5.2: (Left) Schematic view of the convergence domain of the HFA. (Right)Wave-particle energy exchange (colour) together with the convergence domains of HFA and IHFA. For the sake of completeness, the boundary between stabilising and destabilising regions is given by the thin dashed line.

Figure 5 . 3 :

 53 Figure 5.3: Real frequency (left) and growth rate (right) using different methods to write to the dispersion relation.

Figure 5 . 4 :

 54 Figure 5.4: Radial envelope of thermal source (dashed blue line), energetic particle source (red line) and buffer region (shadowed area). The amplitudes are arbitrary.

Figure 5 . 5 :

 55 Figure 5.5: Schematic view of the thermal (left) and energetic particle (right) sources as a function of the velocity and parallel velocity, respectively.

Figure 5 . 6 :

 56 Figure 5.6: (Left) Radial profiles of ion guiding-centre temperature (blue line) and density (red line). The shadowed region represents the radial envelope of the buffer region in arbitrary units. (Right) Radial profile of the turbulent intensity as calculated by expression 5.24, averaged over ω c ∆t ≈ 3 • 10 4 .

Figure 5 . 8 :

 58 Figure 5.8: (Left) Time evolution of the flux-surface averaged electrostatic potential for the three simulations. (Right) For the simulation with S EP /S th = 0.625, time Fourier transform of the electrostatic potential after filtering out the time averaged signal.

Figure 5 . 9 :

 59 Figure 5.9: (Left) Time evolution of the slope of the equilibrium distribution function (thin lines) evaluated at v = qω EGAM R 0 . Time-averages over an EGAM cycle are plotted by thick lines of the same colour. The shadowed region represents the time evolution of the heat flux due to the E × B drift velocity. (Right) Equilibrium distribution functions without energetic particles (thin dashed line) and with energetic particles just before the excitation of EGAMs using S EP /S th = 0.437 (thick dashed line) and S EP /S th = 0.625 (solid line).

Figure 5 .

 5 Figure 5.10: (Left) Time evolution of the flux-surface averaged electrostatic potential with (solid red line) and without (dashed magenta line). (Right) Comparison of the slope of the equilibrium distribution function at the resonant velocity with (red line) and without (magenta line) turbulence.

Figure 5 . 11 : 2 ]Figure 5

 51125 Figure 5.11: Radial dependence of the m = 1 (thick lines) and m = 2 (thin lines) poloidal components of the electrostatic potential oscillating at the EGAM frequency. The sine and cosine components are evidenced by solid and dashed lines, respectively.

Figure 5 .

 5 Figure 5.13: (Left) Radial dependence of the time-averaged E × B heat flux for the simulations of this paper. For comparison, the heat flux in the absence of EP is also plotted and represented by the thin dashed black line. (Right) Time Fourier transform of the E × B heat flux evaluated at the mid radial position for the simulation with S EP /S th = 0.625.

CHAPTER 5 .

 5 INTERACTION BETWEEN EPM AND TURBULENCEFourier transforms of the electrostatic potential and the Morlet wavelet function as F W φm,n = F φm,n F [ψ ω ], with F [•] the Fourier transform. Then, the wavelet transform simply reads W φm,n = F -1 F φm,n F [ψ ω ] , where F -1 [•] is the inverse Fourier transform. In this chapter we restrict the analysis to the interaction between φ0,0 and all the possible poloidal and toroidal modes (m 1 , n 1 ), (m 2 , n 2 ) satisfying m 1 + m 2 = 0 and n 1 + n 2 = 0. This is reflected by the definition of the wavelet bispectrum b

Figure 5 . 14 :

 514 Figure 5.14: Wavelet bispectrum for the four phases indicated in the main text of the paper and on top of each of the subfigures. For clarity and only in the last two figures, the EGAM frequency has been indicated by horizontal yellow dashed lines.

Contents 6 . 1 6 . 2 6 . 3 6 . 4

 61626364 The tearing mode in tokamak plasmas . . . . . . . . . . . . . . . 6.1.1 The tearing mode equation in the outer (ideal MHD) region 6.1.2 The tearing mode equation in the inner region . . . . . . . 6.1.3 Impact of the magnetic field curvature on the tearing instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Impact of the tearing mode on the ITG-dominated turbulence . Profile flattening in the presence of tearing modes and ITGdominated turbulence . . . . . . . . . . . . . . . . . . . . . . . . . Possible mechanism of turbulence reduction . . . . . . . . . . . 6.5 A few words on the transport of particles in the electro-magnetic regime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6 . 1 :

 61 Figure 6.1: Solution of equation 6.5 for different values of ν, q a = 3.5 and ǫ = 1/3.3.

Figure 6 . 2 :

 62 Figure 6.2: Ratio between growth rates with and without trapping electrons in the absence of magnetic drift and comparison with global gyro-kinetic simulations.

Figure 6 . 3 :

 63 Figure 6.3: Energy exchange between the mode and particles with flat profiles and ω D = 0 in the presence of only passing particles (top panel) and when the trapping effects are taken into account (bottom panel).

Figure 6 . 4 :

 64 Figure 6.4: Growth rate of the tearing mode as a function of the stability parameter with flat profiles (left) and with electron temperature gradient (right).

Figure 6 . 5 :

 65 Figure 6.5: Growth rate of the tearing instability as a function of the electron temperature gradient.

Figure 6 . 6 :

 66 Figure 6.6: For the simulation W = 6, ω ils /ω * i = 1, plotted is the turbulent intensity I turb (x) as a function of radial and binormal positions. The region inside the yellow line represents the area around the O-point used to calculate the averaged quantity δI turb . For reference, the separatrix is given by a black line.

Figure 6 . 7 :

 67 Figure 6.7: Scan on island rotation frequency (ω isl ) and island width (W ). The vertical axis represents the modification of turbulence around the O-point of the island, normalized to the reference turbulence value of a simulation without island. A negative (resp. positive) value implies reduction (resp. increase) of turbulence.

Figure 6 . 8 :

 68 Figure 6.8: Normalized density gradient R/L n averaged around the O-point region as a function of the island rotation frequency and the island width. For the sake of clarity only one error bar is provided, measured as the standard deviation in the radial direction of the normalized density gradient.

108 CHAPTER 6 Figure 6 . 9 :

 108669 Figure 6.9: Normalized ion (left) and electron (right) temperature gradients R/L T,i , R/L T,e averaged around the O-point region as a function of the island rotation frequency and the island width. For the sake of clarity only one error bar is provided, measured as the standard deviation in the radial direction of the normalized ion temperature gradient.

Figure 6 . 10 :

 610 Figure 6.10: Time evolution of the turbulent intensity for the CBC simulation (solid line) and the simulation without island but taking into account the profiles given by the presence of the island (dashed line).

  [?] R (∆r, ∆ζ, t) = δφ (r + ∆r, ζ + ∆ζ, t) δφ (r, ζ, t) δφ 2 (r + ∆r, ζ + ∆ζ, t) δφ 2 (r, ζ, t) (6.22) where the brackets are defined as an integral over the whole domain in r and ζ as follows time allows us to define a radial autocorrelation function, which depends only on the radial lag ∆r [?]

Figure 6 . 11 :

 611 Figure 6.11: Radial autocorrelation functions for the ITG parameters used in this paper in the absence of magnetic island.

Figure 6

 6 Figure 6.12: Radial profile through the O-point of the shear of the radial electric field for two island rotation directions: electron diamagnetic direction (dashed line) and ion diamagnetic direction (solid line). The case of an island width W = 6 is considered.

Figure 6 . 13 :

 613 Figure 6.13: Poincaré map of the magnetic field perturbed by a tearing-like instability, leading to the presence of m = 2 and m = 3 islands (left panel), which are located on the resonant surfaces q = 2 and q = 3 (right panel).

Figure 6 .

 6 Figure 6.14: (Left) Trajectory of a lost particle in the poloidal cross-section together with the Poincaré map of the magnetic field perturbed by a tearing-like instability. Blue (resp. red) line represents the unperturbed (resp. perturbed) trajectory. (Right) Perturbed trajectory in the space (ψ, E) diagram as shown in chapter 2 (right).

2.3 Exact linear kinetic response of particles to the electro-magnetic field of a tokamak . . . . . . . . . . . . . . . . . . . . . . . . . . .

	2.2.2 Action-angle variables and motion invariants . . . . . . .
	2.2.3 Passing and trapped particles: orbit width . . . . . . . . .
	2.2.4 Relationship between the angles α and the tokamak co-
	ordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.3.1 Gyro-kinetic approach: when the fastest motion can be
	averaged out . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.3.2 Existence of multiple resonances in the reduced 5D phase-
	space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.3.3 Determining the (n, ω) Fourier mode . . . . . . . . . . . .
	2

.4 Coupling Vlasov and Maxwell's equations . . . . . . . . . . . . .

	2.4.1 The main ingredients of the instabilities: resonance and
	gradient of the equilibrium . . . . . . . . . . . . . . . . . .
	2.4.2 Ampère's law . . . . . . . . . . . . . . . . . . . . . . . . . .
	2.4.3 The quasi-neutrality condition . . . . . . . . . . . . . . . .

As explained earlier, the transport of particles exhibits a chaotic behaviour. This implies that the trajectories of particles are very sensitive to initial conditions. In

We assume experiment-relevant parameters, taking the minor radius of the tokamak (a) and the thermal ion Larmor radius (ρ th ) such that ρ ⋆ = 1/150, with ρ ⋆ = ρ th /a, and we calculate the exit time of an ensemble of counter-passing EP. Such EP are characteristic of NBI heating in medium-size tokamaks like DIII-D. For this purpose, we follow ∼ 4 • 10 5 deuterium tracers initialised at the position r = 0.4a, θ = 0, ϕ = 0, with energy E ≈ 20E th and magnetic moment such that µB 0 /T i = 14. These particles are confined in the absence of any perturbation. We use GYSELA normalizations [82], but one can recover the units by choosing parameters for standard tokamaks. For instance, with T i ≈ 4 keV and B 0 ≈ 2 T, one gets a ≈ 0.67 m, which is typical of medium-size tokamaks. Using the amplitude of the EGAM in nonlinear GYSELA simulations (φ 00 = 1.5), and using the electron temperature T e ≈ 3 keV, the amplitude of the radial electric field is E r,0 ≈ 14 kV • m -1 , which is of the same order as the one obtained in [41]. Despite the simple structure of φ, its time dependence leads to radial transport and losses. The PDF of the exit time, P exit , is plotted in Fig. 4.9 in log -log scale, showing an algebraic decay P exit ∼ t -µe , in contrast with the exponential decay one would expect in the case of a diffusive transport [?]. For the parameters chosen here, the tail of the PDF is developed from 1 -100 ms, which is of the order of, or greater than, the slowing-down time. It is to be noted that a long decay time was mentioned in [41], although no scaling was provided.

To understand why this algebraic decay occurs and therefore why the radial transport is not diffusive in the presence of an oscillating radial electric field, we focus on the region responsible for the chaotic transport of particles, i.e. the stochastic layer separating the passing and trapped particles. Let us remind that, in the absence of any perturbation, the particles in a tokamak are divided into trapped and passing and the boundary between these two classes is called trap-CHAPTER 7. CONCLUSIONS AND FUTURE DIRECTIONS will be devoted to the use of the most advanced computing techniques from Artificial Intelligence to improve our understanding, prediction and control of energetic particles and instabilities in plasmas.

Appendix A

Geometry basics: tensor calculus and toroidal geometry

Note that we give here only a brief introduction to the geometry basics. This is done for the sake of completeness. Nevertheless, the reader can find more material in textbooks such as Ref. [START_REF] Richard | Plasma confinement[END_REF].

A vector field A can be expressed by its components in terms of a basis set, which is a linearly independent set of vectors spanning the set of directions. In Cartesian coordinates, the basis set is composed of a set of mutually orthogonal unit and constant vectors so that the vector field A is decomposed as

In a more general way, we can use a basis set built upon a set of three vector fields {x i } i=1,2,3 that span the set of directions at each point (x, y, z). The covariant representation of any vector field A is written as

where the covariant component is given by the product

the Jacobian in space and ǫ ijk the covariant components of the Levi-Civita tensor. The term covariant comes from the transformation properties of the components. To express the vector in a different covariant system xi we use the chain rule for the derivatives as follows

which gives straightforwardly the covariant components in the basis ∇x i

Since the new components vary with the basis vectors in the sense that x i is in the numerator, this representations is called covariant. However, the contravariant components vary against the basis vectors in the sense that x i appears in the denominator. The contravariant components are obtained by scalar product

Applying the chain rule for the derivatives, we have

The contravariant vectors can be built using expression (A.6), which gives

The relationship between covariant and contravariant components is written

where g ij are the contravariant components of the metric tensor

and the covariant tensor is given by the inverse of the contravariant tensor and satisfies det (g ij ) = J 2

x (A.12)

The gradient of a scalar field A can be written as

and the contravariant components of the rotational of a vector field A are expressed as

We will use the toroidal coordinate system, i.e. the set of coordinates (x 1 , x 2 , x 3 ) is reduced to the radial position and the poloidal and toroidal angles

where r is the radial position, θ is the poloidal angle and ϕ is the toroidal angle. The metric tensor is given by its contravariant components as

APPENDIX A. GEOMETRY BASICS: TENSOR CALCULUS AND TOROIDAL GEOMETRY are then g 11 ≡ g rr = 1, g 22 ≡ g θθ = r -2 , g 33 ≡ g ϕϕ = R -2 and g ij = 0 for all i = j. The major radius of the tokamak is given by R (r, θ) = R 0 + r cos θ (A.17)

We use simplified magnetic geometry, with poloidal cross-sections of the magnetic surfaces as circular and concentric. In this case ∇r •∇θ = ∇θ •∇r = 0 and the metric tensor is diagonal with g 11 ≡ g rr = 1, g 22 ≡ g θθ = r -2 and g 33 ≡ g ϕϕ = R -2 . Therefore, the Jacobian in space is simply J x = rR. The magnetic field is written in the simplified form

with B 0 the magnetic field at the magnetic axis, R 0 the major radius at the magnetic axis, R the major radius at a given position R = R 0 + r cos θ, e θ = r∇θ, e ϕ = R∇ϕ are the unit vectors in the poloidal and toroidal directions, respectively, and ζ = r qR 0 (A. [START_REF] Santos | Study of cancer hallmarks relevance using a cellular automaton tumor growth model[END_REF] where q is the safety factor, which measures the helicity of the magnetic field lines. The covariant and contravariant components of the magnetic field are calculated using expressions (A.3) and (A.6), respectively, which leads to