\thechapter Introduction

Ecole Doctorale de Japan Atomic IRFM, Association l’Ecole Polytechnique Energy Agency Euratom-CEA Route de Saclay 6-9-3 Higashi-Ueno Taito-ku Centre de Cadarache F-91128 Palaiseau, France. 110-0015 Tokyo, Japan F-13108 St Paul-Lez-D., France.

The Berk-Breizman Model as a Paradigm for

Energetic Particle-driven Alfvén Eigenmodes

by

Maxime Lesur

A dissertation submitted for the degree of

Doctor of Philosophy

in

Plasma Physics

Ecole Doctorale de l’Ecole Polytechnique

Dr. Saddrudin Benkadda CNRS Research Director Dr. Xavier Garbet CEA Research Director - Advisor Dr. Virginie Grandgirard CEA Research Scientist Dr. Yasuhiro Idomura JAEA Research Scientist - Advisor Dr. Jean-Marcel Rax Professor at Ecole Polytechnique - PhD Advisor Dr. Laurent Villard Professor at EPFL - Referee Dr. Fulvio Zonca Senior Researcher at ENEA - Referee

Sep 2007 - Dec 2010

## Acknowledgements

A PhD thesis is like a marathon in some sense. More than reaching the finish line, the point is what is learned in the process, from the first training session to the last spark. But, unlike a marathon, which, with enough inner strength, can be achieved without external help, the success of a PhD hinges on the qualities of surrounding people. A few lines in this manuscript is the least I can do to thank those who contributed to the quality of this work, or the quality of my PhD student’s life. And since it is human nature to read the Acknowledgment section only, to check if one’s name is included, I hope I do not forget anybody…

My first thanks naturally go to my two supervisors, Yasuhiro Idomura and Xavier Garbet, who both showed incredible pedagogy and patience, while keeping relationships simple and friendly.

Y. I. has been a supervisor like every PhD student should dream of. He formulated many ideas that are developed in this manuscript, helped me structure my thoughts and writings, instilling his analytic way of thinking, and proofread most of my output with the utmost care, while giving me complete freedom in my final choices. I am very honored to have been his first PhD student, and hope this was not my last opportunity to work with him. I am also grateful to him for helping me set up in Japan and explaining parts of Japanese culture.

X. G. always had time to spare for me, despite incredibly tight agenda. While we had less occasions to interact because of the distance, these interactions were extremely precious to the present work and to me. Since he magically has the answer for any problem, he was like a trump card when I was stuck. I am grateful for his teaching the importance of advertisement and diplomacy in research.

Jean-Marcel Rax kindly accepted to take the formal role of thesis director, and was the eye of the doctoral school of Ecole Polytechnique while I was away.

I sincerely thank Shinji Tokuda for profound lessons on Landau damping and MHD, which were given with extreme pedagogy. Koji Shinohara gave me invaluable help in understanding experimental problematics and diagnostic issues. The quality of this work was greatly improved by eye-opening discussions with Yasushi Todo, Boris Breizman, Matthew Lilley, Simon Pinches, Fulvio Zonca, and Patrick Diamond. I am grateful for fruitful collaborations and nice memories with Christine Nguyen and Nicolas Dubuit.

I want to thank the friends I found in Japan, who made me feel like I was leaving home after three years, while I was coming back to my country. Warm thanks go to my sempai Sébastien Jolliet, for his friendship. I am indebted to him for answering my questions on many topics, and for sharing some of my best memories of Japan. Thanks to my kouhai Miho Janvier for listening to my patronizing speeches and for her words of motivation. Thanks to Yann Camenen for his advices on planning the after-thesis and for his relaxed attitude to which I take example. I also benefited from interactions with Ken Kawaoto and Remy Cardinet, since teaching is often even more formative for the teacher than for the student. I thank all the people of CCSE in Inarichô, especially Masahiko Machida for always having something to chat about. I thank the theory group and the JT-60 team in Naka, Nobuyuki Aiba for frequent rides from the station, Masatoshi Yagi, Yasutomo Ishii, Naoaki Miyato, Makoto Hirota, Junya Shiraishi and Mitsuru Honda for memorable welcoming and farewell parties. I thank Toshio Hirayama and the late Toshihide Tsunematsu for approving my stay.

In France, I owe many thanks to my long-time friends, Audrey, Julie, Christophe, Edouard, and Mathieu, whose trust made me keep looking ahead and reaching forward. Witnessing their friendship, untouched despite the distance, was a real comfort. I also thank my mother who dealt with personal issues in France while I was in Japan, and my whole family for their constant support. I thank Virginie Grandgirard for insightful discussions on numerical issues and for helping me set up in the legendary Hameau ; Guillaume Latu for introducing me to OpenMP and helping me set up on Cadarache super-calculators. I thank my office-mates, Farah Hariri who really brightened my working days and weekends in the Hameau, and David Zarzoso who was so nice with me even when I was using all our computing resources. Thanks to Jérémie Abiteboul for helping me send the manuscript in time, and for his communicative cheerfulness. Thanks to Shimpei Futatani for visiting our office many mornings. Thanks to the other students or ex-students for removing any residual stress, Guilhem Dif-Pradalier, Patrick Tamain, Eric Nardon, Stanislas Pamela, Antoine Merle, Nicolas Mellet, Thimothée Nicolas, Didier Vezinet, Dmytro Meshcheriakov, Hugo Bufferand, Grégoire Hornung, Francois-Xavier Duthoit, Sylvain Rauch, and Jonathan Jacquot. Thanks to the regulars of the 945 AM coffee break for making me laugh every day, Gloria Falchetto, Marina Bécoulet, Chantal Passeron, Yanick Sarazin, Patrick Maget, Joan Decker, Rémi Dumont, Philippe Ghendrih, Simon Allfrey, Guido Huysmans, and Maurizio Ottaviani. I thank Tuong Hoang, Alain Bécoulet, and Xavier Litaudon for their support.

I thank the people of NIFS, CRPP, and NFRI who hosted my visits. I also thank people I met in conferences or seminars, with whom I had insightful discussions, Guo-Yong Fu, Kazuo Toi, Rowdy Vann, Choong-Seock Chang, Chang-Mo Ryu, Raffi Nazikian, Kenji Imadera, and many others. Thanks to the organizers and participants of Festival de Théorie 2007 and 2009, where I learned so much.

Warm thanks to Anthony Genot, who just obtained a PhD from Merton college, and remained an invaluable friend throughout my thesis. He put me back on track when metaphysical doubts were seizing me. He also experimented and reported the various phases of a thesis a few months before me so that I would always be prepared…

For taking care of administrative procedures, I thank Kazutoyo Onozaki, Stéphanie Villechavrolle, and Valerie Icard, who also organized the defense.

Finally I thank the members of my jury, for taking the time to read my manuscript and to attend my defense.

This work was performed within the frame of a co-thesis agreement between the doctoral school of the Ecole Polytechnique, the Japan Atomic Energy Agency, and the Commissariat à l’Energie Atomique. It was supported by both the JAEA Foreign Researcher Inviting Program and the European Communities under the contract of Association between EURATOM and CEA, and by the MEXT, Grant No. 22866086. Financial support to attend the IAEA-FEC 2010 conference in Daejeon was kindly provided by NFRI-WCI Center for Fusion Theory. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Computations were performed on Altix3700 and BX900 systems at JAEA, and on Norma system at CEA.

## Abstract

The achievement of sustained nuclear fusion in magnetically confined plasma relies on efficient confinement of alpha particles, which are high-energy ions produced by the fusion reaction. Such particles can excite instabilities in the frequency range of Alfvén Eigenmodes (AEs), which significantly degrade their confinement and threatens the vacuum vessel of future reactors. In order to develop diagnostics and control schemes, a better understanding of linear and nonlinear features of resonant interactions between plasma waves and high-energy particles, which is the aim of this thesis, is required. In the case of an isolated single resonance, the description of AE destabilization by high-energy ions is homothetic to the so-called Berk-Breizman (BB) problem, which is an extension of the classic bump-on-tail electrostatic problem, including external damping to a thermal plasma, and collisions. A semi-Lagrangian simulation code, COBBLES, is developed to solve the initial-value BB problem in both perturbative () and self-consistent (full-) approaches. Two collision models are considered, namely a Krook model, and a model that includes dynamical friction (drag) and velocity-space diffusion. The nonlinear behavior of instabilities in experimentally-relevant conditions is categorized into steady-state, periodic, chaotic, and frequency-sweeping (chirping) regimes, depending on external damping rate and collision frequency. The chaotic regime is shown to extend into a linearly stable region, and a mechanism that solves the paradox formed by the existence of such subcritical instabilities is proposed. Analytic and semi-empirical laws for nonlinear chirping characteristics, such as sweeping-rate, lifetime, and asymmetry, are developed and validated. Long-time simulations demonstrate the existence of a quasi-periodic chirping regime. Although the existence of such regime stands for both collision models, drag and diffusion are essential to reproduce the alternation between major chirping events and quiescent phases, which is observed in experiments. Based on these findings, a new method for analyzing fundamental kinetic plasma parameters, such as linear drive and external damping rate, is developed. The method, which consists of fitting procedures between COBBLES simulations and quasi-periodic chirping AE experiments, does not require any internal diagnostics. This approach is applied to Toroidicity-induced AEs (TAEs) on JT-60 Upgrade and Mega-Amp Spherical Tokamak devices, which yields estimations of local kinetic parameters and suggests the existence of TAEs relatively far from marginal stability. The results are validated by recovering measured growth and decay of perturbation amplitude, and by estimating collision frequencies from experimental equilibrium data.

## Résumé

Le succès de la fusion nucléaire par confinement magnétique repose sur un confinement efficace des particules alpha, qui sont des ions hautement énergétiques produits par les réactions de fusion. De telles particules peuvent exciter des instabilités dans le domaine de fréquence des modes d’Alfvén (AEs) qui dégradent leur confinement et risquent d’endommager l’enceinte à vide de réacteurs futurs. Afin de développer des diagnostiques et moyens de contrôle, une meilleure compréhension des comportements linéaire et non-linéaire des interactions résonantes entre ondes plasma et particules énergétiques, qui constitue le but de cette thèse, est requise. Dans le cas d’une résonance unique et isolée, la description de la déstabilisation des AEs par des ions énergétiques est homothétique au problème de Berk-Breizman (BB), qui est une extension du problème classique de l’instabilité faisceau, incluant un amortissement externe vers un plasma thermique, et des collisions. Un code semi-Lagrangien, COBBLES, est développé pour résoudre le problème aux valeurs initiales de BB selon deux approches, perturbative () et auto-cohérente (full-). Deux modèles de collisions sont considérés, à savoir un modèle de Krook, et un modèle qui inclue la friction dynamique et la diffusion dans l’espace des vitesses. Le comportement non-linéaire de ces instabilités dans des conditions correspondantes aux expériences est catégorisé en régimes stable, périodique, chaotique, et de balayage en fréquence (sifflet), selon le taux d’amortissement externe et la fréquence de collision. On montre que le régime chaotique déborde dans une région linéairement stable, et l’on propose un mécanisme qui résoud le paradoxe que constitue l’existence de telles instabilités sous-critiques. On développe et valide des lois analytiques et semi-empiriques régissant les caractéristiques non-linéaires de sifflet, telles que la vitesse de balayage, la durée de vie, et l’asymétrie. Des simulations de longue durée démontrent l’existence d’un régime de sifflets quasi-périodiques. Bien que ce régime existe quel que soit l’un des deux modèles de collision, la friction et la diffusion sont essentielles pour reproduire l’alternance entre sifflets et périodes de repos, telle qu’observée expérimentalement. Grâce à ces découvertes, on développe une nouvelle méthode pour analyser des paramètres cinétiques fondamentaux du plasma, tels que le taux de croissance linéaire et le taux d’amortissement externe. Cette méthode, qui consiste à faire correspondre les simulations de COBBLES avec des expériences d’AEs qui présentent des sifflets quasi-périodiques, ne requiert aucun diagnostique interne. Cette approche est appliquée à des AEs induits par la toroidicité (TAEs) sur les machines JT-60 Upgrade et Mega-Amp Spherical Tokamak. On obtient des estimations de paramètres cinétiques locaux qui suggèrent l’existence de TAEs relativement loin de la stabilité marginale. Les résultats sont validés en recouvrant la croissance et décroissance de l’amplitude des perturbations mesurées, et en estimant les fréquences de collision à partir des données expérimentales d’équilibre.

## Abbreviations

BB: Berk-Breizman.
EP: Energetic Particles.
MHD: Magneto-HydroDynamics.
SAW: Shear-Alfvén Wave.
AE: Alfvén Eigenmode.
TAE: Toroidicity-induced Alfvén Eigenmode.
BAE: Beta-induced Alfvén Eigenmode.
KAW: Kinetic Alfvén Wave.
EPM: Energetic Particle Mode.
1D 6D: one dimension six dimensions. Real space dimensions, unless explicitly specified to be phase-space dimensions.
NBI: Neutral Beam Injection.
ICRH: Ion Cyclotron Resonant Heating.
DCE: Displacement Current Equation.
BGK: Bernstein-Greene-Kruskal.
CIP: Constrained Interpolation Profile.
R-CIP: Rational interpolation version of CIP.
CIP-CSL: Conservative Semi-Lagrangian version of CIP.
ITER: International Thermonuclear Experimental Reactor (Cadarache, France).
DIII-D: Doublet III-D (San Diego, USA).
MAST: Mega Amp Spherical Tokamak (Culham, UK).
START: Small Tight Aspect Ratio Tokamak (formerly Culham, UK).
NSTX: National Spherical Torus Experiment (Princeton, USA).
LHD: Large Helical Device (Toki, Japan).
CHS: Compact Helical System (Toki, Japan).
COBBLES: COnservative Berk-Breizman semi-Lagrangian Extended Solver.
HAGIS: Hamiltonian Guiding Center System.

## Chapter \thechapter Introduction

This manuscript has been written in an effort toward pedagogy, privileging clarity over details, physical pictures over mathematical developments. We hope it can be used as an introduction to the Berk-Breizman problem and its numerical computation. However we cannot pretend being able to introduce notions such as nuclear fusion, magnetic plasma confinement, Magneto-HydroDynamics (MHD), kinetic theory, tokamak geometry, or particle orbits, with as much clarity as in well-known textbooks. For this reason, the reader is assumed to be familiar with the basis of the above fields. In this introduction we expose the background and motivation for our work.

### 1 Energetic particles-driven AEs

In an ignited tokamak operating with a deuterium-tritium mix, the confinement of -particles is critical to prevent damages on the first-wall and to achieve break-even. The reason is that these high energy particles must be kept long enough in the plasma core to allow enough of their energy to heat thermal populations by slowing-down processes. A major concern is that high energy ions can excite plasma instabilities in the frequency range of Alfvén Eigenmodes (AEs), which significantly enhance their transport. Ever since the recognition of this issue in the 1970’s, considerable progress has been made in the theoretical understanding of the principal Alfvénic instabilities. However, the estimation of the mode growth rate is complex, and the question of their stability in ITER [ACH01] remains to be clarified. In addition, estimation of the effect on transport and development of diagnostics and control schemes require significant progress on our understanding of nonlinear behaviors.

We limit our framework to the tokamak configuration. In this work, we focus our interest on the Toroidicity induced Alfvén Eigenmode (TAE) [CCC85], which is a shear-Alfvén wave (SAW) located in a gap of the SAW continuum that is created by toroidal coupling of two successive poloidal modes, and which can be destabilized by resonant interactions with high-energy ions. TAEs have been observed to be driven by -particles in burning plamas [WSB96, NFB97], Ion-Cyclotron-Resonance Heating (ICRH) [WWC94], and Neutral Beam Injection (NBI) [WFP91, WSD91]. In this work we consider only the latter situation (NBI-driven TAEs).

### 2 The BB problem as a paradigm for TAEs

In general, TAEs are described in a three-dimensional (3D) configuration space. However, near a phase-space surface where the resonance condition is satisfied (resonant phase-space surface), it is possible to obtain a new set of variables in which the perturbation is described by a one-dimensional (1D) effective Hamiltonian in 2 conjugated variables [BBP97a, BBP97c, WB98, GDPN08], if we assume an isolated single mode. In this sense, the problem of kinetic interactions between TAEs and fast-particles is homothetic to a simple 1D single mode bump-on-tail instability. The so-called Berk-Breizman (BB) problem [BBY93, BBP97a, BBP97c, BBP96] is a generalization of the bump-on-tail problem, where we take into account an external wave damping accounting for background dissipative mechanisms, and a collision operator. Observed qualitative and quantitative similarities between BB nonlinear theory and both global TAE simulations [WB98, PBG04] and experiments [FBB98, HFS00] are an indication of the validity of this reduction of dimensionality. Similar arguments exist for other Alfvén wave instabilities such as the geodesic acoustic mode (GAM) [BBB06], and the beta Alfvén eigenmode (BAE) [Ngu09]. These analogies enable more understanding of fully nonlinear problems in complex geometries by using a model that is analytically and numerically tractable. This approach is complementary to heavier 3D analysis.

### 3 Frequency sweeping

A feature of the nonlinear evolution of AEs, the frequency sweeping (chirping) of the resonant frequency by 10-30 on a timescale much faster than the equilibrium evolution, has been observed in the plasma core region of tokamaks JT-60 Upgrade (JT-60U) [KKK99], DIII-D [Hei95], the Small Tight Aspect Ratio Tokamak (START) [MGS99], the Mega Amp Spherical Tokamak (MAST) [PBG04], the National Spherical Torus Experiment (NSTX) [FGB06], and in stellerators such as the Large Helical Device (LHD) [OYT06], and the Compact Helical Stellerator (CHS) [TTT99]. In general, two branches coexist, with their frequency sweeping downwardly (down-chirping) for one, upwardly (up-chirping) for the other. In many experiments, asymmetric chirping has been observed, with the amplitude of down-chirping branches significantly dominating up-chirping ones. Chirping TAEs have been reproduced in 3D simulations with a hybrid MHD/drift-kinetic code [TSTI03], and with a drift-kinetic perturbative code [PBG04]. Qualitatively similar chirping modes are spontaneously generated by the BB model, and have been shown to correspond to the evolution of holes and clumps in the velocity distribution. Theory relates the time evolution of the frequency shift with the linear drive and the external damping rate [BBP97b], when the evolution of holes and clumps is adiabatic. The idea of using chirping velocity (or sweeping rate) as a diagnostic for growth rates looks promising. If we assume that the mode is close to marginal stability, , and if we further assume that holes and clumps are in the adiabatic regime, then the growth rates are simply obtained by fitting analytic prediction to experimental measurement of chirping velocity. However, these two assumptions are not verified in general, which limits the applicability of this simple approach. In most experiments, chirping events are quasi-periodic, with a period in the order of the millisecond. On the one hand, the long-time simulation of repetitive chirping with 3D codes is a heavy task, which has not been undertaken yet. On the other hand, simulations of the BB model with many chirping events have been performed [VBSC07, LBS10], but such solutions do not feature any quasi-periodicity. In this sense, repetitive chirping that qualitatively agrees with an experiment has not been reproduced, neither with 1D nor 3D simulations.

### 4 Aim of this thesis

Our main goals are to improve our understanding of wave-particle nonlinear resonant interactions, develop diagnostics and identify control parameters for AEs. From these backgrounds, a straight-forward plan is to:

• clarify the link between BB model and AEs,

• extend BB theory where it is insufficient to explain experimental observations,

• analyze experimental data by applying the BB model to AEs.

### 5 Outline

In brief, we reduce the problem of TAEs to a numerically and analytically tractable paradigm, the BB problem. We make a survey of linear and nonlinear physics of the BB model, and work on improving our understanding, by extending theory and by using systematic numerical observations, focusing on the frequency sweeping regime. Armed with new findings and robust numerical tools, we finally go back to the original problem of TAEs, explaining quantitative features of experimental observations, and developing a new method for accurate linear predictions.

In Chap. \thechapter, we review the physics of the TAE, and simplify the description of the problem from 6D to 2D in phase-space, around a single resonant phase-space surface. The first step in this procedure is to isolate the gyromotion, which is decoupled from the wave for typical TAEs with kHz, by changing variables to the guiding-center coordinates. This change of variables is based on the so-called Lie transform perturbation theory, which we review in order to get a better grasp of involved hypothesis. Then we show how to reduce the guiding-center Hamiltonian, as well as collision operator and background mechanisms, and discuss limitations of this simplified description. In Chap. \thechapter, we review basic nonlinear physics of the BB problem. We develop and verify numerical tools that we use in following chapters. In particular, we develop and verify a kinetic code to solve the initial-value problem. We show that numerical simulations in experimentally relevant conditions, with cold-bulk and weak, warm-beam distributions, require a careful choice of advection scheme. Among the family of so-called Constrained Interpolation Profile (CIP) schemes, a locally conservative version (CIP-CSL) displays the best convergence properties. As an intermediate summary, we cast the analogies between BB model and 1D description of TAEs. In Chap. \thechapter, we investigate kinetic features of self-consistent interactions between an energetic particle beam and a weakly unstable electrostatic wave in 1D plasma, within the BB framework, with an emphasis on chirping regime and subcritical instabilities. We show that the nonlinear time-evolution of chirping in 1D simulations can be used to retrieve information about linear input parameters with good precision. We identify a regime where chirping events are quasi-periodical. This regime exists whether the collision model is annihilation/creation type, or takes into account dynamical friction and velocity-space diffusion. Based on these findings, in Chap. \thechapter, we propose a new method to estimate local linear drive, external damping, and collision frequencies from the spectrogram of magnetic field variations measured by a Mirnov coil at the edge of the plasma. This method, which relies on a fitting of normalized chirping characteristics between the experiment and BB simulations, is described and applied to TAE experiments in MAST and JT-60U. We show that the BB model can successfully reproduce features observed in the experiment only if the collision operator includes drag and diffusion terms. We find a quantitative agreement for the diffusion frequency, and a qualitative agreement for the drag frequency, between the values obtained with our fitting procedure, and an independent estimation obtained from experimental equilibrium measurements. In the conclusion (Chap. \thechapter), we summarize our findings, and we suggest numerical and experimental investigations these findings enable.

## Chapter \thechapter Energetic particle-induced TAEs

The TAE is a shear-Alfvén wave located in a gap of the SAW continuum that is created by toroidal coupling of two successive poloidal modes, which can be destabilized by energetic particles. This chapter deals with the modelization at several levels of complexity of an isolated single-mode TAE. In Sec. 6, we give a short review of the basic TAE physics, and of the state-of-the-art of its linear stability analysis, emphasizing a need for improved accuracy. When the TAE is an isolated single-mode, it is possible to reduce the problem to a simple harmonic oscillator, by expanding the perturbed Hamiltonian around a resonant phase-space surface. This reduction from 3D to 1D becomes particularly transparent in angle-action variables (, ), once the perturbation has been put in the form , where are integers. In Sec. 7, we show how to obtain the latter form. The first step is to separate the gyromotion, which does not resonate with typical TAEs, by changing variables to guiding-center coordinates. This procedure can be done while conserving the Hamiltonian properties by making use of Lie transform perturbation theory, which we review thoroughly. In Sec. 8, we show how to reduce the Hamiltonian, collision operator, and background damping mechanisms from their 3D description to a 1D model.

### 6 Review

This short review is intended to introduce notions and notations that we use when applying the BB model to the TAE, and to motivate our linear analysis of TAEs. For a more comprehensive review of energetic-particle driven AEs, see Ref. [Hei08].

#### 6.1 Toroidal geometry

Hereafter, we make use of magnetic flux coordinates [RDH03], also called as straight field-line coordinates, (, , ), where , and are minor radius, poloidal and toroidal angle, respectively, to describe the nested poloidal flux surfaces of the equilibrium magnetic field . In these coordinates, takes the following form,

 B0=∇χ×∇(qθ−ζ), (1)

where is the poloidal flux (divided by ), and the safety factor, defined by

 q(r)≡B0\boldmath⋅∇ζB0\boldmath⋅∇θ, (2)

is the flux surface-averaged number of toroidal rotation that a field line undergoes in one poloidal rotation.

#### 6.2 Physics of the TAE

In an homogeneous magnetized plasma, linear ideal-MHD arguments show the existence of a shear-Alfvén wave of frequency with the dispersion relation

 ω2A=k2∥v2A, (3)

where

 vA=B0√μ0∑inimi (4)

is the Alfvén velocity, and is the wave number in the direction of the equilibrium magnetic field . Let us consider axisymmetric toroidal plasmas. In the cylindrical limit, the periodicities of the system require that there exists two integers, a toroidal mode number and a poloidal mode number , such that

 k∥=n−m/q(r)R0, (5)

where is the distance from the symmetry axis of the tokamak to the magnetic axis. In a non-homogeneous plasma in a sheared magnetic field, both and are functions of . The simple dispersion relation Eq. (3) is still valid in this configuration [CC86], and it is called the Alfvén continuum. Since phase velocity is a function of radius, a wave packet with finite radial extent would suffer from phase-mixing, the so-called continuum damping. Except for energetic particle modes, which are outside of the scope of this work, resonant drive by fast particles is not enough to overcome this damping. However, a toroidal coupling of two successive poloidal modes and breaks up the continuous spectrum. This is illustrated in Fig. 1, which shows the Alfvén continuum for , and , in cylindrical geometry, where two poloidal continuum are decoupled, and in toroidal geometry, with a two-mode coupling model [CC86, FD89]. The latter is obtained with equilibrium plasma parameters corresponding to JT-60U shot E32359 at s, which we analyze in Chap. \thechapter, assuming concentric circular magnetic flux surfaces, retaining toroidicity effects in the first order in inverse aspect ratio. Though we show only the half-plane, the continuum is symmetric with respect to . Coupled modes are (, ) and (, ) for , and (, ) and (, ) for . The gap is centered at a radius such that , where the two continuous spectra would cross in the absence of coupling, and where . The resulting discrete eigenmode is a TAE, at a frequency .

For a deuterium plasma with typical magnetic field and density , the Alfvénic energy is , which is in the range of passing particles induced by NBI. For ITER parameters, , which is in the range of passing -particles born from fusion reactions. In both cases, TAEs can be driven unstable by resonance with energetic particles. For far-passing particles, the resonance condition is , where

 Ω=nωζ+lωθ, (6)

where and are frequencies of toroidal motion and poloidal motion, respectively, and for co-passing particles, for counter-passing particles [TS98]. Since we analyze TAEs driven by co-injected ions, we can simplify following discussions by considering only co-passing particles. Then, the resonance condition is

 ωA−nv∥R0+mv∥qR0=0. (7)

#### 6.3 Linear stability

In theory, the linear growth rate is positive when the drive by fast particles overcomes damping processes to background plasma. The growth rate can be estimated either by linear stability codes, such as PENN [JAVV95], TASK/WM [FA03], NOVA-K [Che92], or CASTOR-K [BBB02] ; or by gyro- or drift- kinetic perturbative nonlinear initial value codes, such as FAC [CBB97] or HAGIS [PAC98]. The analysis requires internal diagnostics that are not always available.

The linear drive depends on several factors such as spatial and energy gradients of EP distribution, and alignment between particle orbit and the eigenmode. In the limit of large aspect ratio, analytic theory [FD89] yields successful estimations of for well-defined numerical models [TSW95]. However, in the general case, complicated factors cited above forbid accurate analytic estimation, as yet. Moreover, improvements in measurement of EP distributions are needed to allow estimation of for experimental TAEs.

The external damping rate involves complicated mechanisms, which include continuum damping (since parts of the mode extend spatially into the continuum) [ZC92], radiative damping [MM92], Landau damping with thermal species [BF92, ZCS96], and collisional damping by trapped electrons [GS92]. Most of these mechanisms are still under investigation, and cannot be quantified by existing theory. Experimentally, can be estimated by active measurements of externally injected perturbations [FBB95, FBB00]. However, the applicability of this technique is limited to dedicated experiments.

Moreover, if the system is close to marginal stability, is sensitive to small variations of driving and damping terms, and is very sensitive to plasma parameters such as the profile. In addition, the existence of unstable AEs in a regime where linear theory predicts , or subcritical AEs, has not been ruled out. To access the subcritical regime, nonlinear analysis is necessary. Therefore, accurate linear stability analysis requires significant theoretical and experimental improvements, or a new approach.

### 7 Angle-action description

#### 7.1 Kinetic description

Since the kinetic equation is at the center of interest of this thesis, it is useful to get a deep understanding of its meaning, by reviewing its derivation from fundamental principles. The steps we summarize here are detailed in Ref. [KT86] for example.

A many-body system is completely described by the microscopic distribution , where , and are the position and velocity of a particle indexed as , and the sum is taken over all particles. To simplify the present discussion, we consider a single-species system of particles, without external forces, and normalize the total phase-space volume, assumed constant, to . Substituting Newton equations of motion into the partial time derivative of yields the so-called Klimontovich-Dupree equation,

 ∂N∂t+v\boldmath⋅∂N∂x+amicro\boldmath⋅∂N∂v=0, (8)

where is the acceleration due to microscopic electromagnetic fields, at the exception of those due to a particle located at , if such a particle exists. The microscopic electromagnetic fields obey Maxwell equations, where density and current are velocity moments of .

Since it is impossible to reproduce any many-body experiment at the microscopic level, it is much more efficient to take an ensemble point-of-view, where distributions and fields are smooth functions of phase-space. The statistical properties are completely determined by the distribution , where is the probability of finding, at a time , particle 1 within , particle 2 within , and particle within , and is an infinitesimal phase-space volume at the neighborhood of . By integrating over all but , we can define the one-particle distribution function , where is the probability of finding one particle within at . is related to the microscopic distribution by

 f1(z1,t)=⟨N(z1,t)⟩Np, (9)

where is the ensemble average . Similarly, by integrating over all but and , we can define the two-body distribution function , which is related to the microscopic distribution by

 f2(z1,z2,t)=⟨N(z1,t)N(z2,t)⟩N2p−δ(z1−z2)Npf1(z1,t), (10)

where we used a large particle number approximation, . In the absence of atomic and nuclear processes, is a constant, then the ensemble average of Eq. (8) yields

 ∂f1∂t+v\boldmath⋅∂f1∂x+a\boldmath⋅∂f1∂v=df1dt∣∣∣coll., (11)

where the average acceleration is given by electromagnetic fields that obey Maxwell equations, where density and current are velocity moments of . The collision term,

 df1dt∣∣∣coll.≡a\boldmath⋅∂f1∂v−1Np⟨amicro\boldmath⋅∂N∂v⟩, (12)

is shown to vanish in the absence of particle interactions, yielding Vlasov equation, which is valid on a time-scale much shorter than a collisional time-scale. Eq. (11) is not a closed equation, because the second part of the collision term involves expressions of the form . Under the Coulomb approximation, which forbids any retardation effect, and which is valid if the thermal velocity is much slower than the speed of light, we can reduce the latter term as a function of . Similarly, the equation which gives the evolution of involves terms of the form , and so on. Altogether, we have a chain of equations for which we need a closure. A collision operator is an approximative statistical operator that accounts for particle interactions, which provides such closure. A clear review of collision operators is given in Ref. [HS02]. In this thesis, we focus on a Fokker-Planck collision operator, which is based on the fact that, in a Tokamak plasma, collisions are dominated by binary Coulomb collisions with small-angle deflection.

In the following, we write simply as . The kinetic equation, Eq. (11), can be put in Hamiltonian form,

 Extra open brace or missing close brace (13)

where is the Hamiltonian, and are Poisson brackets. The l.h.s. is the variation of following particle orbits, or so-called Lagrangian derivative of , noted .

The fact that, in the absence of collisions, is conserved along particle orbits, can be seen as a direct consequence of Liouville’s theorem, which states that the density in phase-space is constant along particle orbits. Let us consider an infinitesimal volume of phase-space surrounding a particle at at , as illustrated in Fig. 2. The particle and all points of evolve following the equations of motion until a time where the particle is at , and is changed to a volume . Liouville’s theorem arises from the following properties,

• Since the points of the boundary of follow the equations of motion, in the absence of collisions, the number of particles inside is constant;

• Time-evolution can be seen as a series of infinitesimal canonical transformations generated by the Hamiltonian;

• Poincaré’s invariant: the volume element in phase-space is a canonical invariant.

#### 7.2 Review of Lie transform theory

When working in canonical variables, it is easy to exploit the properties of a Hamiltonian system. However, for some systems it is difficult to find canonical variables, and the most convenient variables may not be canonical. Moreover, if some variables are canonical for the unperturbed system, they may be non-canonical for the perturbed one. By applying the Lie near-identity transformation theory to the phase-space Lagrangian, one can change any Lagrangian into a simpler form in coordinates that reveal the symmetries of the system, and use this formulation to study the properties of a perturbed Hamiltonian system in any arbitrary noncanonical variables [CL83].

Noncanonical Hamiltonian Mechanics Consider a Hamiltonian system with degrees of freedom. Hamilton equations are given by the application of a variational principle to the scalar Lagrangian , in some arbitrary coordinates on the -dimensional extended phase space, (, ),

 δ∫L(zμ(λ),λ)dλ=0, (14)

where is an arbitrary parameter. Hereafter, Greek indices , and run from to , whereas Latin indices and run from to , and , , from to . Thus, any Hamiltonian system is characterized by its Lagrangian

 L=γμdzμdλ, (15)

or equivalently by its fundamental one-form, . In the canonical extended phase-space , when is the time coordinate,

 γ0 = −h(zμ,t), (16) γa = pa, (17) γa+N = 0, (18)

where is the Hamiltonian. However, in noncanonical variables, all the may be nonzero.

From Eq. (14), which has the same form in any new coordinate system with , Euler-Lagrange equations are obtained as

 ωμνdzνdλ=0, (19)

where is a tensor defined by

 ωμν=∂γν∂zμ−∂γμ∂zν, (20)

and its restriction to the symplectic coordinates is the Lagrange tensor. Thus the flow is an eigenvector of with eigenvalue . The solution is unique only after a normalization. In physical terms, we can take . Since the Jacobian of the coordinate transformation is nonsingular, we can invert the Lagrange tensor. The Poisson tensor then yields the equations of motion,

 dzidz0=Jij(∂γj∂z0−∂γ0∂zi). (21)

Lie transform technique We consider a fundamental one-form which consists of an equilibrium part for which the flow is well-known, and a small perturbation . We want to study the effect of the perturbation on the flow. The strategy is to search for a near-identity transformation that will reveal the symmetries of the perturbed system. Indeed, if the fundamental one-form is independent of some coordinate , then, as an application of Noether’s theorem,

 dzαdt=∂γα∂zμdzμdt=−ωαμdzμdt=0, (22)

revealing an exact invariant. The general form of a near-identity transformation with a small parameter is

 Zμ=zμ+ϵZμ1f(z)+ϵ2Zμ2f(z)+… (23)

Rather than an expansion in which is difficult to invert, we express the transformation in operator form. We denote a forward transformation , and a backward transformation . In the Lie transform technique, the coordinate transformation is specified by a generator such that

 ∂Zμf∂ϵ(z,ϵ)=gμ(Zf(z,ϵ)), (24)

and .

The forward transformation of a scalar into is given by

 F = e−ϵLgf, (25) f = e+ϵLgF, (26)

where is defined by its action on a scalar , and its action on a one-form . The transformation of coordinates is just a special case of scalar transformation,

 Zαb=e−ϵLgIα, (27)

where is the coordinate function. The transformation of a one-form into is given by the functional relationships,

 Γ = e−ϵLgγ+dS, (28) γ = e+ϵLgΓ+ds, (29)

where and are scalar functions.

Higher order perturbation theory We now consider a fundamental one-form in the form of an expansion . We introduce a push-forward transformation operator , where each is a Lie transform operator, and is a short for . Substituting these new definitions into Eq. (28), we have , which yields for each successive order, , and

 Γ1 = dS1−L1γ0+γ1 Γ2 = dS2−L2γ0+γ2−L1γ1+12L21γ0 … Γn = dSn−Lnγ0+Cn,forn≠0. (30)

Let us recall that our aim is to simplify the fundamental one-form in the new coordinates . To this aim we have to solve successively the latter equations for the gauge scalar and the generating vector (Each is given by and the result of the preceding order ). In these equations, the generating vector appears only in the one-form . We already discussed that has a unique null eigenvector. Then we can add any multiple of this eigenvector to without changing the equations we are now trying to solve. Let us suppose this eigenvector has a nonzero component in the time direction, then we can set

 g0n=0, (31)

so that the time-coordinate doesn’t change (). At this point, we are left with unknowns, namely generating functions and one scalar gauge . A priori we should be able to bring the components of into a simpler form. A good strategy is to make all its symplectic components vanish by choosing

 gjn=(∂iSn+Cni)Jij0. (32)

Let us now focus on the time component of the new one-form,

 Γn0=−Hn=∂tSn−gjnˆω0j0+Cn0. (33)

It is convenient to introduce the lowest order velocity vector, defined as the time derivative along the unperturbed orbits of the coordinates :

 Vμ0≡dzμdt∣∣∣0. (34)

Substituting the unperturbed equations of motion (21), we find that the scalar gauge is given by its total time derivative over the unperturbed orbit,

 DSnDt∣∣∣0≡Vμ0∂Sn∂zμ=Γn0−Vμ0Cnμ. (35)

In integrating the latter equation, we want to avoid any secularity effect. Then we should remove any secular perturbation by taking

 Γn0=⟨Vμ0Cnμ⟩0, (36)

where the average is to be taken over the unperturbed orbits.

Finally, in the new coordinates , the new one-form is given by

 Γn=⟨Vμ0Cnμ⟩0dt, (37)

if we choose the generating functions

 gjn=(∂iSn+Cni)Jij0, (38)

where the gauge scalar is given by

 Sn=−∮0˜Vμ0Cnμ, (39)

and a tilde represents the oscillatory part.

To illustrate the benefits of Lie transform theory compared to classical perturbation theory, and to quantify its validity limit, we apply it to a simple mathematical problem in App. \thechapter.

#### 7.3 Guiding-center Lagrangian

Guiding-center theory provides reduced equations of motion of a charged particle in a slowly-varying (in time and space) electromagnetic field, where the fast gyromotion is decoupled from the relatively slow drifting motion of the guiding-center. Guiding-center theory is based on the closeness to the limit of a fixed and uniform magnetic field, where the orbit of a particle in a frame following its gyration center is a circle. The perturbation from this situation is quantified by a small parameter , if we assume the following ordering,

 ϵ∼ω/Ωc∼ρ/L≪1, (40)

where is the cyclotronic frequency (or gyrofrequency), is the Larmor radius (or gyroradius), is a characteristic frequency of interest, and is the scalelength of field variation.

On the one hand, the standard derivation of guiding-center equations of motion is based on an averaging procedure [Nor63], which removes important properties of the Hamiltonian formulation. On the other hand, the modern derivation [Lit83] is based on Lie-transform perturbation theory. This approach preserves the validity of Noether’s theorem and the validity of Liouville’s theorem, to each order in . In addition, expansion to arbitrary order is straightforward. Moreover, since we keep the Hamiltonian formulation, further reductions of dimensionality are still possible with Lie transforms.

The starting point is the one-form in canonical cartesian phase-space (, ),

 γ=p\boldmath⋅dx−h(x,p,t)dt, (41)

where , and is the vector potential. The one-form can be expressed in noncanonical phase-space , and written as an expansion in ,

 γ=γ0+ϵγ1+…, (42)

where

 γ0 = eA(x,t)\boldmath⋅dx−eφ0(x,t), (43) γ1 = mv\boldmath⋅dx−(eφ1(x,t)+m2v2)dt, (44)

and where is the scalar potential.

In Ref. [Lit83], Lie-transform theory is applied to change variables to new coordinates, where the new one-form does not depend on the gyroangle. When this procedure is carried up to the second order, we obtain the guiding-center one-form,

 Γ=(mUb(R,t)+eA(R,t))% \boldmath⋅dR+Mdξ−Hdt, (45)

where

 H=eφ(R,t)+m2U2+(e/m)MB(R,t), (46)

and the guiding-center coordinates as

 R ≡ x−ϵv⊥Ωcia, (47) U ≡ v\boldmath⋅b=v∥, (48) M ≡ mv2⊥2Ωci, (49) ξ ≡ tan−1(v\boldmath⋅e1v\boldmath⋅e2), (50)

where is an orthogonal unit vectors system, , and .

In Eq. (45), the angle-action variables of the gyromotion, (, ), appear in canonical form in the symplectic part, and does not depend on . Thus, since the gyromotion is irrelevant to fast-ion/TAE interactions, we can drop the term in the guiding-center one-form.

#### 7.4 Guiding-center action-angle variables

In Ref. [MH90], the zeroth order guiding-center one-form is put in the form,

 Γ0=Pθdθ+Pζdζ+Mdξ−H(Pθ,Pζ,θ,M)dt, (51)

where is the toroidal canonical angular momentum,

 Pζ≡−eχ(r)+mbR0v∥, (52)

and is the poloidal canonical angular momentum. The latter expression can be used as a starting point to develop an action-angle formalism where the perturbed Hamiltonian takes a standard form. In Ref. [BBP95b], a canonical transformation is performed to obtain

 Γ0=Jθd˜θ+Jζd˜ζ+Jξd˜ξ−H(Jθ,Jζ,Jξ)dt, (53)

where , , and are the unperturbed frequencies of the gyromotion, poloidal motion, and toroidal motion, respectively. and reduce to the geometric angles and if we neglect finite aspect ratio effects (But we do not neglect them, since we would remove the toroidicity from which the TAE originates).

#### 7.5 Application to TAE

Although TAEs resonate mainly with passing particles, when the source of high-energy ions is isotropic, a large fraction of the energy transfer may be accounted by resonance with the bounce-motion (or banana motion) of toroidally trapped particles [TS98]. However, for tangential NBI ions, to which we confine our analysis, it is sufficient to describe resonance with far passing particles.

In a gauge where the perturbed scalar potential is zero, the TAE can be described by a perturbation Hamiltonian,

 H1=−eA1⊥\boldmath⋅vgc, (54)

where is the guiding-center velocity. Here, is the perpendicular part of the perturbed vector potential, and we have neglected a second order term . In a small plasma pressure (small ) limit, we can neglect parallel gradients, then can be split into magnetic compression,

 Ac1⊥=b0×∇Π, (55)

and magnetic shear,

 As1⊥=∇Φ−(b0% \boldmath⋅∇Φ)b0, (56)

where . For the TAE, which is a shear Alfvén wave, the latter part only is relevant, hence the excitation is described by a single scalar function .

In Ref. [BBP95b] the perturbed Hamiltonian is put in the form

 H1=V(J)cos(l\boldmath⋅α−ωt), (57)

in arbitrary tokamak geometry, where (, ) are angle-action variables for the solvable unperturbed Hamiltonian , and is a triplet of three integers. Substituting Eq. (56) into Eq. (54) yields

 H1 = −evgc\boldmath⋅∇Φ+e(b0\boldmath⋅∇Φ)b0\boldmath⋅vgc (58) = −edΦdt∣∣∣0+e∂Φ∂t+e(b0\boldmath⋅∇Φ)v∥, (59)

where is the Lagrangian derivative of along unperturbed particle orbit, which can be removed from the Lagrangian without altering Euler equations.

For a single mode of frequency , the eigenfunction takes the following form,

 Φ(r,θ,ζ)=C(t)e−ıωt−ıφ(t)m+1∑l=mΦl(r)eınζ−ılθ+c.c., (60)

where and are the amplitude and phase of the wave. Substituting the eigenfunction and the expression into Eq. (59), the TAE excitation is reduced to

 H1=−ıeC(t)e−ıωt−ıφ(t)m+1∑l=mΦl(r)eınζ−ılθ[ω−v∥b0\boldmath⋅(n∇ζ−m∇θ)]+c.c., (61)

where we have neglected time-derivation of slowly varying phase and amplitude of the wave.

Then we change the variables to the canonical angle-actions (, , , ), in order to express the perturbed Hamiltonian as a Fourier series in ,

 H1=−ıeC(t)e−ıωt−ıφ(t)eın˜ζ+∞∑p=−∞Vp(Jθ,Jζ)eıp˜θ+c.c., (62)

where

 Vp=∫d˜θ2πe−ıp˜θ−ın˜ζm+1∑l=mΦl(r)eınζ−ılθ[ω−v∥b0\boldmath⋅(n∇ζ−m∇θ)]. (63)

Formally, the problem is expressed in the desired form of Eq. (57), but the numerical computation of the Fourier components , which is needed for quantitative comparison of absolute physical quantities between 3D and 1D model, requires to relate geometric angles and canonical angles, which may be complicated depending on the equilibrium configuration. However, since each particle of the resonant phase-space surface interact with the wave in the same way at the same frequency (though with different strength), comparison of quantities that are normalized to the mode frequency is possible, even without evaluating coefficients.

### 8 Reduction to a one dimensional problem

If we consider only small toroidal mode numbers , and modes are isolated. Let us consider a single toroidal mode number. On the one hand, since on a flux surface where a poloidal mode number is centered, the safety factor is , then we can estimate the distance between two neighbouring modes by writing , as

 Δr≈1nq′. (64)

On the other hand, the characteristic width of TAE modes is of the order of [CCC85]. Hence, for typical parameters, , where is the magnetic shear. Thus, TAEs have a two-scale radial structure, the larger scale corresponding to the enveloppe of the TAE. In our analysis, we assume that the number of harmonics involved is small enough to consider resonances one by one, as isolated , mode. The latter hypothesis is reasonable for sufficiently core-localised, low- TAEs. We must keep in mind, though, that high- TAEs are likely to be destabilized in future devices such as ITER, in which case it may be necessary to take into account multiple-wave resonances.

The evolution of the distribution of energetic ions is described in 3D configuration space by the kinetic equation (13), with the perturbed Hamiltonian Eq. (57). In the following, we reduce the problem to a 1D Hamiltonian system, by considering a single , mode.

#### 8.1 Reduction of the Hamiltonian

Formally, the resonant phase-space surface, , is defined by a function . The resonance condition,

 l\boldmath⋅Ω(JR)=ω, (65)

where , is satisfied on the resonant phase-space surface.

Once the perturbed Hamiltonian has been put in the form of Eq. (57), we can reduce the problem to one action and one angle [Lic69, GDPN08], by performing a canonical transformation with the generating function

 S=−I\boldmath⋅ψ+I3(l%\boldmath$⋅$α−ωt)+I1α1+I2α2+F(I1,I2)α3. (66)

This procedure yields,

 J1= I1+l1I3 ψ1= α1+α3∂I1F J2= I2+l2I3 ψ2= α2+α3∂I2F J3= F(I1,I2)+l3I3 ψ3= l\boldmath⋅α−ωt,

and . Thus, near the resonant phase-space surface, , and we can expand the new Hamiltonian around this surface,

 H′(ψ,I) = H0(JR+I3l)+V(JR+I3l)cosψ3−I3ω (68) = H0(JR)+I3(l\boldmath⋅Ω(JR)−ω)+12DI23 +V(JR+I3l)cosψ3,

with .

If the variations of are small around , we can replace by in the latter expression, and obtain the new Hamiltonian , with

 H1,JR(ψ,I)≡12DI2+Vcosψ. (69)

Thus, the problem has been reduced to a 1D Hamiltonian problem for the angle-action variables (, )(, ).

Substituting the expression of the TAE perturbation, Eq. (62) into Eq. (57), we obtain

 ψ = p˜θ+n˜ζ−ωt, (70) I = Jζ−JζRn≈−eχ−χRn, (71) D ≈ n2∂2H0∂J2ζ(JR), (72) V = −ıeC(t)Vp(JθR,JζR), (73)

where the subscript means that the quantity is evaluated at the resonance, and with , , and