Analysis of the nonlinear dynamics of a chirping-frequency Alfv\'en mode in a Tokamak equilibrium

Chirping Alfv\'{e}n modes are considered as potentially harmful in burning Tokamak plasmas. In this paper, the nonlinear evolution of a single-toroidal-number chirping mode is analysed by numerical particle simulation. This analysis can be simplified if the different resonant phase-space structures can be investigated as isolated ones. This can be done adopting a coordinate system that includes two constants of motion. In our simulations, we adopt as constants of motion, the magnetic momentum and the initial particle coordinates. For each resonant structure, a density-flattening region is formed around the respective resonance radius, with radial width that increases as the mode amplitude grows. It is delimited by two large negative density gradients, drifting inward and outward. With constant mode frequency, this density flattening would be responsible for the exhausting of the drive when large negative density gradients leave the resonance region. The frequency chirping, however, causes the resonance radius and the resonance region to drift inward. This drift delays the moment in which the inner density gradient reaches the inner boundary of the resonance region. On the other side, the island reconstitutes around the new resonance radius; as a consequence, the large negative density gradient further moves inward. This process continues as long as it allows to keep the large gradient within the resonance region. When this is no longer possible, the resonant structure ceases to be effective in driving the mode. To further grow, the mode has to tap a different resonant structure, possibly making use of additional frequency variations.


I. INTRODUCTION
Alfvén modes can be driven unstable, in Tokamak plasmas, by the resonant interaction with alpha particles produced by fusion reactions and/or energetic ions produced by auxiliary heating methods, characterised by speeds of the order of the Alfvén velocity [1][2][3][4][5][6][7][8][9].Modeparticle interaction can in turn deteriorate the confinement of these particles, preventing their thermalisation in the core plasma and eventually damaging the first wall.Some of the Alfvén modes, like the Toroidal Alfvén Eigenmodes [10] (TAEs) have a MHD counterpart; that is, they exist as marginally stable or quasi marginally stable modes even in the absence of energetic-particle drive.Other modes, like the Energetic Particle Modes [11] (EPMs) have no such counterpart: they are oscillations of the Alfvén continuum driven unstable only when the energetic-particle drive exceeds the continuum damping; in the absence of this drive, they are strongly damped.While the former modes are typically characterised by an approximately constant frequency, even in the nonlinear stage, essentially constrained to keep within the frequency gaps opened in the Alfvén continuum by the coupling between different poloidal harmonics, the latter are able to vary their frequency (they have been named chirping modes), during the nonlinear evolution, in order to extract as much power as possible from the energetic particles.Both because of this feature and the strong dependence on the energetic-particle drive they usually present above the instability threshold, such chirping modes are considered as potentially harmful for the confinement of energetic particles in burning Tokamak plasmas and have attracted much interest.
In this paper we want to investigate, by numerical particle simulation, the nonlinear dynamics of chirping modes.For the sake of simplicity, we consider the case of an Alfvén spectrum characterised by a single toroidal mode number n.We have already studied the case of a single-n mode, with constant frequency ω in Refs.[12,13].Our approach consisted in adopting a coordinate system that includes two invariants: namely, the magnetic moment M and the quantity C ≡ ωP φ − nE, with P φ and E being, respectively, the toroidal angular momentum and the energy of the particle [14].The phase space can then be seen as a set of slices characterised by given values of M and C. Because of the invariance of these two coordinates, particles cannot cross their birth slice and the gradients of the distribution function orthogonal to the slices do not play any direct role in mode-particle power exchange: the slices can then be treated as isolated ones.Once the most relevant resonances (that is, the slices yielding the largest contribution to the mode drive or damping) are identified, saturation dynamics can be analysed by focusing on the evolution of each of them.To this aim, in Refs.[12,13] we have sampled the selected slices by a large number of test particles moving in the fields computed in a full-population, self-consistent simulation, and its interaction with the mode has been analysed in detail.We have shown that an island-like structure, enclosing the bounded orbits of particles instantaneously trapped in the potential well of the wave, forms around the resonance radius in the plane (Θ, r eq ); here, Θ corresponds to the wave-phase seen by the particle and r eq is the "equatorial radius" (see Sec. VI and, respectively, Sec.III for their definition).Mixing of trapped particles originating from the opposite sides of the resonance radius occurs.This gives rise to a flattened-density region delimited by large negative density gradients moving inward and outward as the island width grows with the increasing mode amplitude.For fixed mode frequency saturation is reached when the density flattening completely covers the resonant-interaction region [12,13].The width of this region is typically limited by the smallest of the resonance and mode widths: in the former case, the saturation mechanism has been named resonance detuning; in the latter, radial decoupling.
In the case of chirping modes, we can expect a more complex saturation dynamics.Indeed, the capability of the mode of changing its frequency implies both that the resonant interaction of a given slice can be prolonged by counteracting the saturation mechanism just described, and that other phase-space regions can take over from the initially dominant one in destabilising the mode.Unfortunately, a slavish application of the above procedure is not possible, as the quantity C is no longer an exact invariant of the particle motion when ω is not a constant.The invariance of C can be, of course, approximately satisfied if the rate of variation of ω is small enough.Moreover, the possibility of generalising the definition of C in such a way to get an isolating invariant cannot be ruled out [15][16][17].In the present paper, however, we adopt a different choice: including in the phase-space coordinate system a constant of motion rather than a dynamical invariant; namely, a suited function of the initial values of the particle coordinates.Here, we use the words constant of motion with the trivial meaning of a quantity which remains constant along the particle trajectory: the initial values of particle coordinates have, of course, this feature, then labelling each particle in a permanent way.It is worth observing that resorting to the initial values of particle coordinates (or to a suited function of them) is a viable choice only in the frame of numerical simulation, where the transformation that brings from the actual coordinates to the initial ones is readily available.
We will show that our approach allows us identifying the different resonant structures contributing to mode drive during the evolution of the system.In the case analysed in this paper, it is possible to take two of these structures as representative of the mode-particle interaction during the linear stage and, respectively, the nonlinear one.In both cases, the relatively small growth rate and the quite large mode structure would make the resonance detuning the dominant saturation mechanism in the constant-frequency case.We will see, however, that the capability of the mode to change its frequency (a downward chirping, in the considered case) alters this mechanism.Indeed, it causes the resonance radius and resonance region to move inward.This drift allows to delay the moment in which the inner density gradient reaches the inner boundary of the resonance region and leaves it.On the other side, the island reconstitutes around the new resonance radius, drifting inward too.
As a consequence, the large negative density gradient further moves inward.This process goes on as long as the frequency can decrease and the resonance region can move inward in such a way to keep the large gradient within the resonance region.When this is no longer possible, either because a further change in frequency is disadvantageous in terms of drive/damping balance, or because such a change does not result in a significant inward shift of the resonance radius and resonance region, the gradient ceases to be effective in driving the mode.To further grow, the mode has to tap a different resonant structure, possibly making use of additional frequency variations.We will analyse the analogies and the differences between the behaviours of the two resonant structures, identified as the representative ones for the linear and, respectively, the nonlinear stages.
The paper is structured as follows: Sec.II shows how to compute the power transfer in a coordinate system different from that adopted to push particles in the phase space.The only requirement is that the alternative coordinates of each particle are known in terms of the standard ones.Section III presents some useful choices for such alternative coordinates.
In particular, some advantages in adopting the so-called exact-invariant coordinates in the constant-frequency single-toroidal-number case are discussed.For more general cases, the possibility of resorting to constants of motion instead of invariants is proposed, and it is observed that, in the view of a numerical simulation approach, a suited set of constants of motion can be immediately recognised in the initial coordinates of the particle.The numerical experiment performed to analyse the dynamics of a chirping mode is described in Sec.IV.Section V presents the search for the most relevant phase-space resonant structures.In Sec.VI, the Hamiltonian-mapping test-particle approach is adopted to investigate some aspect of the nonlinear dynamics of a resonant structure.Section VIII analyses the phenomenon of particle trapping and de-trapping in detail.The relationship between modeparticle power transfer and fulfilment of the resonance condition is examined in Sec.VII.
Section IX shows the evolution of density-flattening and resonance regions accompanying the mode chirping and describes the succession of different resonant structures in driving the mode during the nonlinear stage.Finally, a short summary of the paper and a discussion of the critical points of our approach are presented in Sec.X.

NATES
In this section we want to show how the calculation of the power exchange between mode and particles can be performed in the framework of a gyrokinetic simulation, and how it is possible, from a numerical point of view, to arbitrarily choose the set of phase space coordinates to represent this exchange.
Gyrokinetic simulations consist in solving the Vlasov equation for the particle distribution function for each relevant species, coupled with the equations for electromagnetic fields, in the low-frequency limit, that is for phenomena characterised by frequencies much smaller than the Larmor frequency of each species.The equations for the fields can be in the form of Poisson-Ampere equations (fully gyrokinetic codes [18,19]) or MHD equations containing suited momenta of the particle distribution functions (hybrid MHD-gyrokinetic codes [20])1 .
Here, we consider for simplicity the case of a hybrid MHD-gyrokinetic code and concentrate on the energetic-particle population, though all the following treatment can be immediately generalised to any code and species.
The simulation is performed by discretising the phase space into microscopic volumes, representing each volume by a marker (macro-particle) that brings the whole electric charge contained in the volume, computing the electromagnetic fields in terms of the momenta of the marker distribution function and updating the phase-space coordinates of each marker according to the fields it experiences.
Let us adopt, for the phase space, the gyrocenter coordinate system Z ≡ (r, θ, φ, M, U, ϑ), where r is the radial coordinate, θ and φ are the poloidal and toroidal angle, respectively (cf.Fig. 2 below), M is the conserved magnetic momentum, U is the parallel (to the equilibrium magnetic field) velocity and ϑ is the gyrophase.Discretising the phase space allows us to represent the particle distribution function in terms of markers in the following way: Here, F H is the energetic ("hot") particles, D zc→Z is the Jacobian of the transformation from canonical to gyrocenter coordinates, Z l = Z l (t) are the gyrocenter coordinates of the l-th macroparticle, which evolve according to the equations of motion.Moreover, and Note that the integration in Eq. 1 extends to a 5-D space, as D zc→Z F H does not depend on the gyrophase ϑ.Note also that ∆ l , the phase-space volume element corresponding to the l-th macroparticle, is a constant of motion of the l−th macro-particle (Liouville's theorem).
The power transfer from energetic particles to the wave is given by −dE H /dt, where is the total energy of fast particles, is the single particle energy, m H and e H are the energetic-particle mass and electric charge, respectively, c is the speed of light, and δϕ and δA are the gyro-average of the fluctuating scalar potential and, respectively, the parallel (to the equilibrium magnetic field) component of the fluctuating vector potential.We can then write After integrating over the gyrophase ϑ and replacing the quantity D zc→Z F H by its discrete form, Eq. 1, we get In the following, we will indicate by Z the coordinates (r, θ, φ, M, U ), neglecting the gyrophase.
In order to identify the wave-particle resonances responsible for the destabilisation of the mode, it is worth reducing the dimensionality of the phase space by averaging the contribution of each marker over the poloidal and toroidal angles.In other words, we can define a power-transfer density in the 3-D space (r, M, U ) in the following way with P (r, M, U, t) obtained from Eqs. 1 and 7 in the form Let us now consider a different coordinate system, Z.We can write the power transfer rate as Note that the transformation Z → Z could be very complicate and possibly a time dependent one.Equation 11 however shows that, if we choose these coordinates to perform the computation instead of the Z ones, the only further quantities we need to know are the Z coordinates of each marker, while we do not need to compute the Jacobian of the transformation (which could be, as we will see in the following, a difficult or even impossible task).
If the coordinate system Z includes, as Z does, the angles θ and φ, we can define a power transfer density P in the reduced space (r, M , Ũ ) as or Here, we have used the notation (r, M , Ũ ) for the other Z coordinates only conventionally, as there is no constraint on the meaning of any of these coordinates.

III. USEFUL CHOICES OF PHASE SPACE COORDINATES
Among all the possible choices of phase-space coordinate system, some are particularly suited to yield pregnant information about particle dynamics.So, it can be worth replacing r and U by the angular momentum and the kinetic energy E; in Eq. 15 R is the major-radius coordinate, ψ eq (r, θ) is the poloidal flux of the equilibrium magnetic field B eq ≡ R 0 B φ0 ∇φ + R 0 ∇ψ eq × ∇φ and the subscript 0 indicates the quantities computed at the equilibrium magnetic axis.P φ and E are invariants of the unperturbed motion, and this allows to immediately obtain evidence, in the evolution of particle coordinates, of the effects of the mode-particle interaction.
An alternative and useful choice is that of the equatorial coordinates r eq and U eq , here defined as the value that the radial coordinate and, respectively, the parallel velocity of the particle with actual coordinates (r, θ, φ, M, U ) would assume at the next crossing of the equatorial plane at θ = 0 if its motion were unperturbed.In particular, for trapped 2 particles, we can convene to refer to the outermost equatorial crossing.These coordinates can be numerically computed, as functions of the Z coordinates, from the conservation, in the unperturbed motion, of P φ , M and E. One of the merits of this choice is that the equatorial coordinates, along with the magnetic momentum M , are able to immediately identify unperturbed particle orbits, clearly separating passing particle contributions from trapped particle ones 3 .This is only an approximation, of course, as particle motion is perturbed by fluctuating electromagnetic fields, but it is a satisfactory one until the perturbation is weak; that is, in the linear stage.During the nonlinear stage the description offered by this coordinate choice is still able to enlighten mode-particle dynamics, provided that the nonlinear perturbation of particle orbits keeps relatively small during a transit or bounce period.
A third possible choice resorts to exact invariants of the (perturbed) motion, if we are able to identify them.There is an immediate advantage in adopting coordinate systems that include such invariants: namely, the unambiguous identification of the most active particles in driving or damping the mode at each time.It can happen, during the nonlinear stage, that the mode is driven/stabilised by particles whose Z coordinates are different from those of the particles driving/stabilising the mode during the linear stage.An interesting issue is distinguishing whether the nonlinearly driving/stabilising particles are actually the same as in the linear stage, with nonlinearly modified Z coordinates, or different particles.In the former case, we could describe the mode-particle power exchange as a process in which the mode causes resonant particles to modify their orbits but, at the same time, adapts itself to further extract energy from them.In the latter case, the interaction between mode and resonant particles would apparently be described as a continuous search, by the mode, for the particles that can transfer power more efficiently.
2 Note that here the word trapped means "trapped in the magnetic well".When the word means "trapped in the potential well of the wave", it is explicitly indicated. 3Trapped-particle coordinates satisfy the following condition: , π − ψ eq (r eq , 0) − m H c eR 0 (R 0 + r eq )U eq > 0.
Let us assume that the coordinates (r, M , Ũ ), or, more generally, some of them, are exact invariants of the perturbed motion.If we define (r max (t), Mmax (t), Ũmax (t)) the coordinates of the maximum power transfer at time t, any variation of the value corresponding to an exact invariant would indicate that new particles have replaced the previous ones in destabilising or stabilising the mode.
In practice, M can be considered as a conserved quantity up to the required asymptotic expansion order in the collisionless gyrokinetic limit.So, a succession of different values of M would readily be recognised as succession of different relevant particles taking part in the considered physical process.If, instead, the coordinates of the maximum power transfer corresponds to r max (t), M max , U max (t), with constant M max (this is typically the case of modes driven unstable by transit resonances), the Z coordinates do not allow us to distinguish between the two above situations (same particles versus succession of different particles).In order to succeed in this task, we have to identify more invariants, if they exist.
In particular, if two more invariants can be selected, we get a complete identification, at any time, of the most driving or stabilising particles.
Another advantage in adopting a phase-space coordinate system including exact invariants of motion is that we can reduce the dimensionality of the free-energy source.This can be seen in the following way.The Vlasov equation for the energetic-particle distribution function can be written in the form with FH (t, Z) = F H (t, Z( Z)).If Zj is an invariant of the motion, we will have d Zj /dt = 0.
Then, the gradient ∂ FH /∂ Zj will not take part in the mode-particle dynamics.We can express the same concept in a different way: if we cut the phase space into slices orthogonal to the invariant coordinate, there will be no particle flux from one slice to the other; that is, the nonlinear particle evolution will not yield any change of the distribution function along that coordinate.The dynamics of a given slice will then depend on other slice dynamics only through the fields, not because of mixing of the respective populations.This feature allows us to investigate separately each slice.Once the slices providing the most relevant resonances are identified, the details of mode-particle interaction can be analysed, e.g., by the so-called Hamiltonian mapping technique, which greatly enhance the resolution with which the resonance is examined by sampling the corresponding slice with a large number of test particles evolving in the fields self consistently computed by the considered simulation.
The advantage of this approach is twofold: on one side, the dynamics of resonant particles is not obscured by the behaviour of non resonant ones; in particular, the local modifications of the distribution function can be highlighted even in the case of weak resonances, in spite of a general insensitivity of the overall population to the presence of the mode.On the other side, examining the evolution of a single-slice dynamics can make describing and identifying the saturation mechanisms easier.
In previous works [12,13], we have investigated the nonlinear dynamics of Alfvén modes characterised by a single toroidal number n and a constant frequency ω.In such cases, the quantity C = ωP φ −nE is an exact invariant of the perturbed motion.Indeed, the equations of motion in the Hamiltonian form read and where H is the single particle Hamiltonian.The conservation of C immediately follows [14] from the dependence on time and toroidal angle of the Hamiltonian in the form H = H(ωt− nφ).In those cases, we adopted a coordinate system Z C ≡ (r eq , θ, φ, M, C), performing the Hamiltonian mapping analysis of phase-space slices characterised by fixed values of M and C. As explained in Ref. [12], once a relevant slice has been fixed, the analysis proceeds storing the information concerning each test-particle orbit at its equatorial plane crossing; in particular, its coordinates r eq and φ (θ = 0 by definition) and the power transferred from/to the mode during the last poloidal orbit.In this way, the distribution function for particles belonging to the considered slice can be computed, along with its integral over φ, which depends only on r eq and represents the free-energy source with which the slice can contribute to the destabilisation of the mode.Time evolution of such integral can then be related to that of the power transfer, getting insight in the saturation dynamics.
If multiple toroidal numbers are present and/or the frequency is not constant, C is no longer an exact invariant of the perturbed motion.This is the general case, of course, both because a large n spectrum of Alfvén modes can be destabilised in tokamak equilibria and chirping modes (that is, modes with time-varying frequency) are often observed in experiments.In these cases, slices orthogonal to the C direction would not be isolated and the dynamics would be influenced also by the gradient of the distribution function along C and the corresponding fluxes crossing each slice.In other words, the evolution of the freeenergy source associated to a certain slice could not be adequately reproduced by simply looking at the gradient of the distribution function along r eq .
It is however possible to envisage other, alternative, coordinate transformations, capable to yield the same advantages as the "exact-invariant" one.In this paper, we propose to look at simple constants of motion (that is, quantities that are conserved during particle motion), without addressing the question whether an isolating integral or global invariant exists [16,17] besides M (and can be considered instead of C) and which form it assumes.
The reason is that any constant of motion could be used instead of the additional global invariant (isolating integral) of motion, as it will be shortly discussed at the end of this Section.Identifying constants of motion is an elementary task, as the initial coordinates of a particle are, by definition, conserved along its motion.The transformation that links a certain coordinate system, Z, to the system represented by its initial values Z0 can be obtained by the equations of motion in the following form Zi More generally, we can consider a coordinate transformation in which only some of the Zi are replaced by the corresponding Zi 0 , or by any function of them; in the following, any of these coordinate system will be named an initial-value coordinate system.Note that the choice of an initial-value coordinate system is, in general, impractical in the context of an analytical treatment, as it requires complicated backward transformations from the current particle coordinates to the initial ones.In the frame of a numerical approach, however, it is sufficient to memorise, in addition to the current coordinates of each marker, its initial coordinates.
On the basis of the arguments exposed in this Section, we find worth considering a threeconstant system, Z3 ≡ (r eq0 , θ, φ, M, U eq0 ) for the identification of instantaneous maxima of the power transfer, and a two-constant system, Z2 ≡ (r eq , θ, φ, M, U eq0 ), for investigating the nonlinear saturation dynamics.Here, r eq0 and U eq0 are the initial values of the equatorial coordinates r eq and U eq defined above.Note that in the latter system the first coordinate is the current value of the equatorial radial coordinate: it is not a constant of motion.The corresponding gradient of the distribution function plays the role of a free-energy source for the mode destabilisation, and the associated fluxes contribute to mode saturation.
In concluding this section, we would like to emphasise the following points.For the purpose of a complete and unambiguous identification of instantaneous power transfer maxima, the use of any triad of invariant coordinates is equally effective, regardless of whether they are global invariants or not.As far as the goal of simplifying the description of nonlinear dynamics is concerned, the conclusion is similar, but requires some clarification.First of all, it is apparent how the choice of any invariant coordinate (besides M ) allows a subdivision of the phase space into slices that do not exchange particles with each other.The fact, however, that these slices are treated as isolated is a mere artifice, consisting in treating the fluctuating field as an exogenous field: this is what is done in sampling each of the considered slices by a large number of test particles that passively evolve according to the field calculated, previously, in a self-consistent simulation.In the evolution of the real system, however, each slice communicates with the others, though not exchanging particles, precisely through the field.Moreover, one slice can take over from another in destabilising the mode, due to the modification of the frequency and/or the structure of the mode.One cannot therefore expect a priori that the relevant dynamics will remain confined to a particular slice.In principle it is also possible that such a confinement occurs in correspondence with the choice of a given invariant, but, so far, no recipe is known to guide this choice in the general case: indeed, having such a recipe would be equivalent to a full ability to predict the nonlinear evolution of the system, even in very complex situations.And all this is irrespective of the particular invariant (global or not) included in the coordinate system.From this point of view there is, then, in general, no qualitative difference between the adoption of a coordinate system including a global invariant and that of a system containing any other quantity conserved during the motion of the particle: in both cases, only a posteriori one can establish whether the choice made is effective for the purpose of a simplified description: namely, whether it allows to restrict the analysis to a small number of relevant slices.

IV. SETTING OF THE NUMERICAL EXPERIMENT
In this Section, the numerical experiment envisaged for analysing the dynamics of a chirping Alfvén mode is presented.
We have simulated, by the hybrid MHD-gyrokinetic code XHMGC [25][26][27], the evolution of an n = 3 Energetic Particle Mode [11], with poloidal harmonics m = 4 ÷ 9, in a large aspect ratio Tokamak equilibrium (R 0 /a = 10, with a being the minor radius of the torus), characterised by a safety factor q monotonically increasing from 1.6 to 2.8 (Fig. 1).The spatial coordinate system is shown in Fig. 2. A deuterium plasma is considered, along with an energetic deuteron population.Both bulk-ion and energetic-ion populations are treated kinetically 4 .The relevant dimensionless parameters are the following: n H0 /n i0 = 0.01, v H0 /v A0 = 0.3, ρ H0 /a = 0.01, v i0 /v A0 = 0.06, ρ i0 /a = 0.002.Here, n H and n i are the densities of the two species, v H and v i their thermal velocities, ρ H and ρ i their thermal Larmor radii, v A is the Alfvén speed and subscripts 0 denote quantities computed at the magnetic axis.The density profile of energetic particles is shown in Fig. 3, while their temperature profile, as well as bulk-ion density and temperature profiles are assumed flat.
The initial distribution function is Maxwellian for both species.For the energetic particles, it is multiplied by an anisotropy factor Ξ(α; α 0 , κ), with cos α = U/(U 2 + 2M Ω H /m H ), α 0 = 0 and κ = 0.1.The small initial relaxation due to the fact that such functions are not equilibrium ones is inhibited.
Figure 4 shows the time evolution of mode energy and frequency.The first saturation of the mode is followed by a stage in which the mode gets further drive and saturates again, several times.A significant decrease of frequency is observed in the nonlinear phase (see also Fig. 5, which shows the power spectrum of the scalar potential in the space (r, ω) at two times: t = 240τ A0 (with τ A0 ≡ R 0 /v A0 ), during the linear phase, and t = 600R 0 /v A0 , in the nonlinear one).In spite of the significant nonlinear variation of the frequency, the radial structure of the mode is scarcely modified, as it can be seen from Fig. 6.

V. SEARCH FOR THE MOST RELEVANT PHASE-SPACE STRUCTURES
Our aim is yielding a detailed description of the nonlinear evolution of the mode.As stated before, the frequency variation causes C to lose its exact invariance.Although, in the presente case, the rate of variation of C keeps quite small, instead of controlling the amount of non conservation and evaluating its relevance, we prefer to include in the coordinate system different conserved quantities, as proposed in Sec.III: the initial values of minor radius and/or parallel velocity; more precisely, those of the equatorial minor radius and/or the equatorial parallel velocity.Then, while pushing particles in the usual gyrocenter coordinates Z, we adopt two alternative systems: a three-constants system, Z3 ≡ (r eq0 , θ, φ, M, U eq0 ) for the identification of instantaneous maxima of the power transfer, and a two-constant one, Z2 ≡ (r eq , θ, φ, M, U eq0 ), for investigating the nonlinear saturation dynamics, with r eq0 and U eq0 being the initial values of the equatorial coordinates r eq and U eq .
The first issue we want to face is the identification of the most relevant resonant structures; that is, the phase-space regions playing the most important role in destabilising or stabilising the mode.To this aim, we consider the power transfer from energetic particles to the mode it can be more interesting referring the power transfer to an equatorial-coordinate system, Z ≡ (r eq , θ, φ, M, U eq ).As the triad (r eq , M, U eq ) fully identify a poloidal particle orbit, this choice is justified by the interest in identifying the contribution of each particle along its whole orbit, rather than the instantaneous one.To state this in a slightly different way, until a certain particle maintains its unperturbed orbit, it will contribute to the power transfer to the same point in the space (r eq , M, U eq ), while it would spread its contribution on a curve U (θ) = U (r(θ), M ) in the space (r, M, U ).
Figure 7 shows the time evolution of the coordinates of the maximum of P (r eq , M, U eq , t) for the considered case.We see that the drive is always related to co-passing particles (M = 0, U eq > 0).During the nonlinear stage the maximum power transfer coordinates r eq max (t) and U eq max (t) change, while M max (t) keeps constant.This time dependence is consistent with both a nonlinear drive yielded by the same particles responsible for the linear destabilisation, with coordinates progressively modified by the interaction with the mode; and with a drive yielded by a succession of different particles.As discussed in Sec.III, the fact that there is no time dependence of M max (with M being the only invariant coordinate in the considered system) makes it impossible to readily show which of these two interpretations is correct.To do so unambiguously, we adopt the Z3 coordinate system and, after averaging the power transfer density over the toroidal and poloidal angles, look at its maxima in the reduced space (r eq0 , M, U eq0 ) (Fig. 8).In this case, the time variation of r eq0 max (t) and U eq0 max (t) clearly indicates that the maximum drive is generally due, at different times, to different particles (characterised by different initial coordinates).
Figure 9 shows the quantity P versus time at the grid points, in the discretised space (r eq0 , M, U eq0 ), that result to host the maximum value for some time, in the time interval [0, 800τ A0 ]5 .We see that grid points dominating the drive during the linear stage, progressively lose importance, being replaced by grid points less and less contributing to linear destabilisation.
Computing power transfer in the coordinate system Z3 is worth to label particles providing the maximum-drive at different times, but it is not suited for investigating the saturation dynamics.To address this issue, we have to adopt a coordinate system with only two constants of motion, in such a way that particle fluxes, responsible for the relaxation of the FIG.9: Time evolution of the power transfer (averaged over θ and φ) at the most relevant grid points in the discretised space (r eq0 , M, U eq0 ).
Only those grid points that yield a maximum of the power transfer at some time in the interval [0, 800τ A0 ] have been considered.Grid points are labelled by three indices.The first one, relative to the radial coordinate, ranges from 46 (r eq0 0.48a) to 56 (r eq0 0.58a).The spacing between adjacent points is ∆r eq0 0.01a.The second one, relative to M , is 0 for all grid points (M = 0).The third one, relative to the parallel velocity, ranges from 10 (U eq0 1.48v H0 ) to 14 (U eq0 2.07v H0 ), with a grid spacing ∆U eq0 0.1475v H0 .
free-energy source and the related saturation mechanism can be observed.Then, we resort to the Z2 = (r eq , θ, φ, M, U eq0 ) coordinate system, containing the two constants of motion M and U eq0 .In this way, the gradient of the distribution function along the equatorial radial coordinate r eq represents the free-energy source, and the corresponding flow contributes to mode saturation.yielded by the most relevant phase-space slices.Each slice corresponds to a grid point in (M, U eq0 ).As in Fig. 9, only grid points yielding a maximum power transfer at some time in the considered interval are plotted.Grid point labels, relative to M and U eq0 , have the same meaning as in Fig. 9.
In principle, we can expect that in the present case, differently from the cases characterised by single toroidal number and constant frequency [12,13], in addition to the linearstage resonance, further resonances may contribute the mode drive during the nonlinear stage.Then, we have first to identify the most relevant structures in the 2-D reduced phase space (M, U eq0 ).Once these structures are identified, we can analyse each of them independently of the others.Indeed, the phase space can be cut into distinct slices, orthogonally to the subspace (M, U eq0 ), with flows entering or leaving the slices being forbidden by the conservation of M and U eq0 .
Figure 10 show the time behaviour of the power yielded by the most relevant slices.As in Fig. 9, only those grid points that yield a maximum power transfer at some time in the considered interval are plotted.Consistently with the previous observation, we find that all the relevant slices are characterised by M = 0. Figure 11 shows the time evolution of the 2-D maximum coordinates obtained by approximating, as done for Fig. 8, the average (integrated over all the other coordinates) at t = 300τ A0 and t = 564τ A0 .We observe that, although qualitatively similar, Fig. 10 presents a simpler situation than Fig. 9.The reason is, of course, that the former is obtained by integration over the radial coordinate, and it is insensitive to redistribution of the power transfer density within the same (M, U eq0 ) slice.In the following, we will take advantage from this simplification, concentrating our analysis on two slices only: namely, the slice dominating during the linear stage (M = 0, U eq0 1.92) and that dominating around t = 564τ A0 (M = 0, U eq0 1.63), after the contribution of the former one has strongly decreased.We will conventionally refer to these two slices as the linear slice and the nonlinear one.

VI. HAMILTONIAN MAPPING ANALYSIS
In order to analyse the behaviour of each of the two phase-space structures identified in the previous Section, we apply the Hamiltonian mapping approach, described in previous papers [12,13].The slice is sampled by a large number of test particles, initialised with the following coordinates: r = r min ÷ r max , with [r min , r max ] being a radial interval covering the whole radial domain of interest (here, r min = 0.4a, r max = 0.65a), θ = 0, φ = 0 ÷ 2π, M = 0, U = Ūeq0 , with Ūeq0 being the value characterising the slice.As all particles are initialised at θ = 0, the initial values of r and U assume the meaning of the equatorial coordinates r eq and U eq .Test particles are pushed in the fields computed self-consistently, at each time step, in the considered simulation.Their coordinates are detected at each crossing of the equatorial plane, and the wave-particle phase Θ ≡ t 0 dt ω(t ) + mθ − nφ is calculated.Note that we have defined the phase taking into account that the mode frequency is not a constant.The wave-phase at the j−th crossing will then be given by Θ j = t j 0 dt ω(t ) + 2πjmσ − nφ j , where σ ≡ sign(U ), t j is the time at which the crossing occurs and φ j the value assumed by the toroidal angle at that time.The resonance condition can be written as with k = 0, ±1, ±2, ... being the bounce harmonic.Here, ∆φ j ≡ φ j − φ j−1 .The average power exchanged with the mode along the poloidal orbit is also computed.
Figure 13 shows, for each test particle, a coloured marker in the plane (Θ, r eq ), where Θ ≡ Θ module 2π.From Eq. 22, we see that the resonance conditions reads, in terms of Θ, To get a clearer view of particle evolution, a companion marker is also drawn at Θ + 2π.
Marker colour is chosen according to the birth r eq value of the particle.Three times are considered, relative, respectively, to linear (t = 100.8τA0 ), early-nonlinear (t = 351.0τA0 ) and fully-nonlinear stage of mode evolution (t = 480.τA0 ).In the unperturbed motion, FIG.13: From left to right: test-particle markers in the (Θ, r eq ) plane at three successive times of the considered simulation: t = 100.8τA0 (left), t = 351.0τA0 (center) and t = 480.6τA0 (right).Each marker is coloured according to the birth r eq value of the particle.
r eq is constant.Then, during the linear phase of the mode evolution (Fig. 13-left), particle trajectories in the (Θ, r eq ) plane essentially reduce to fixed points for r eq = r eq res (in this case, r eq res 0.54a), while they correspond to drift along the Θ axis in the negative (positive) direction, for r eq greater (less) than r eq res .As the mode amplitude grows (Fig. 13 center), r eq varies because of the mode-particle interaction (e.g., radial E × B drift).Even particles that were initially resonant are brought out of resonance, drifting in phase and eventually undergoing an inversion of the drift in r eq .Particles that cross the r eq = r eq res line revert the phase drift as well.Thus, their trajectories are bounded, giving rise to the formation of an island-like structure in the (Θ, r eq ) plane (Fig. 13 right).Its radial extension grows with mode amplitude, consistent with equations of motions accounting for radial E × B drift and particle motion due to radial magnetic field perturbation.Particles outside the island undergo a secular drift in phase, as the E × B drift is not able to cause them crossing the resonance radius.
The formation of the island mixes particles born on both sides of the resonance radius, causing a density flattening around that radius.This can be seen from Fig. 14, where the initial density profile and the flattened one at t = 480.6τA0 are compared.At the same time, large negative density gradients emerge at the boundaries of this flattening region6 .This flattening process is the basis for mode saturation, or, to be more precise, for the exhausting of the destabilising contribution of the considered phase-space slice.If the mode were constrained to keep constant its frequency and radial structure, saturation would be reached as soon as the flattened-density region extends over the whole resonant interaction region.This, in turn, would be limited by the intersection between the resonance region (defined as the region, around the resonance, in which the rate of phase drift is smaller than a certain amount proportional to the linear growth rate) and the radial extension of the mode.Typically, this intersection coincides with the narrowest of these two regions.If it is the resonance region, the saturation mechanism has been dubbed resonance detuning; if it is the radial extension of the mode, radial decoupling [9,12,13,15].If the equilibrium allows the mode to modify its frequency and/or structure, several elements can make it able to further extract energy from particles, prolonging its growth.such that the exhaustion of the contribution of one resonance does not imply a definitive saturation of the mode.We will see in Sec.IX that, in the present case, the latter two elements are essential in determining the nonlinear evolution of the mode.

VII. POWER TRANSFER RATE AND RESONANCE CONDITION
Before addressing the problem of mode saturation in greater detail, we want to investigate the behaviour of the most destabilising particles, within the same slice, and its connection with their capability of fulfilling the resonance condition.
Figure 15 shows analogous plots to those presented in Fig. 13.In this case, however, markers are coloured according their own power-transfer rate (averaged over the last poloidal orbit).Two times are considered: t = 248.4τA0 and t = 419.4τA0 .We see that, consistently with island formation, density flattening and large negative gradients appearance, the power transfer maxima, originally located around the resonance radius, separate from each other, following, respectively, the inner and the outer gradient.We wonder whether the most destabilising particles at t = 248.4τA0 (represented by the red markers) are the same that drive the mode at t = 419.4τA07 .To answer this question, we plot, in Fig. 16  The right frame makes apparent the position, at t = 419.4τA0 , of the t = 248.4τA0 most destabilising particles.Comparing it with the right frame of Fig. 15, it is evident that the mode is driven by not the same particle clusters at the two considered times.
Given that a particle does not permanently belong to the group of the most destabilising particles, we want to check now whether the power exchange between the mode and a certain particle is limited to only a short time interval or it can recur later.To this aim, we follow (Fig. 17 colours to markers in Fig. 19, we conventionally consider a band as "wave trapped" if at least one of its particles has its last phase variation smaller than a certain value (0.015, in this case), corresponding to the resonance condition fairly well satisfied by at least one particle.Fig. 20 reports, at each time, the indices of actually wave-trapped bands.We observe a progressive, asymmetric (mainly inward) trapping of bands.From this analysis, no de-trapping clearly emerges.The reason is that a band is considered "wave trapped" if at least one of the particles it includes is fairly resonant.A different picture can be obtained by looking at the individual behaviour of test particles.In Fig. 21 only markers corresponding to particles that have already fairly satisfied the resonance condition are plotted ("already-trapped" particles).Colour is red for particles that appear still trapped in the wave ("actually-trapped" particles); it switches to blue as the particle cumulates a phase drift greater than a certain conventional threshold (2.5π, in this case).The same two times considered in Fig. 19 are examined.We see that, at t = 726τ A0 , a significant amount of particle de-trapping takes place on the outer side of the island.Here, the de-trapping phenomenon becomes apparent because each particle is considered according to its specific dynamics, which can be different from that of closely born particles.
The continuous trapping/de-trapping phenomenon can also be observed by resorting to the Lagrangian Coherent Structure (LCS) technique [29][30][31][32][33].It is well known that these to diverge from each other, as the first particle maintains an unbounded orbit, while the second gets trapped.The second and third particles tend to move together (both come from unbounded-orbit region; both get trapped).The third and fourth particles come from well separated positions (the fourth particle comes from the region of trapped particles) and converge towards close positions (both in the trapped-particle region).This results, that is substantially independent on the specific particles examined, show that particles that enter the lower channel become trapped.Figure 24 shows the same LCSs, along with the trajectories of four different particles, located, at the reference time, just below (the first particle), within (second and third particles) and just above (the fourth particle) the upper channel.
It can be seen that the second and third particles become de-trapped, while the first particle remains trapped and the fourth one maintains an unbounded orbit.These findings can be easily justified repeating, mutatis mutandis, the above considerations concerning the relative positions of each couple with respect to repulsive and attractive lines.We then arrive at the specular conclusion: particles that enter the upper channel get de-trapped.
The possibility of distinguishing actually-trapped particles from already-trapped ones (Fig. 21) allows us to evidence a subtle distinction between the island structure and the density-flattening region.Indeed, the island only includes, by definition, actually-trapped particles.However, particles that get de-trapped continue to take part to the density flattening, as they contribute to the mixing phenomenon around the resonance radius.Then, the density flattening region extends over the whole region including the set of alreadytrapped particles (both actually-trapped and de-trapped particles).This can be seen from domain, with a trade off between the need of accuracy and that of containing noise.It can be seen that identifying the boundaries of the already-trapped set with the radial position of the gradients is a fair approximation.Then, in the following, we will assume that such boundaries adequately represent the large-gradient positions and, then, the boundaries of the flattened-density region.On the other side, we have observed that, due to the inward island drift (consequent to the downward frequency chirping), trapping of new particles takes place mainly on the inner side of the island structure8 ; particle de-trapping, on the outer side only.Therefore, the inner boundary is the same for the region including the actually trapped particles (the island) and that including the already trapped ones (the flattening As last remark, we want to emphasize that this trapping and detrapping process that accompanies frequency chirping fluctuations has been recently pointed out by hybrid simulations [36] and corresponding theoretical framework [37,38] of "chorus emission" in the Earth's magnetosphere, suggesting the universal nature of the underlying nonlinear dynamics.

IX. FREQUENCY CHIRPING AND RESONANCE EVOLUTION
We want to investigate how the collective power transfer is related to the the fulfilment of the resonance condition and the evolution of the free-energy source (namely, the radial We have already observed that, during the linear stage (Fig. 15-left), the maximum power transfer is yielded by resonant particles, around r eq 0.54a.In the early nonlinear stage, the island formation and growth and the corresponding flattening of the density profile around the resonance radius is accompanied by the splitting of the power transfer maximum into two separate maxima, each of them following one of the two large negative density gradients delimiting the flattening region and moving, respectively, inward and outward (Fig. 15 right).
Figure 26 compares the radial boundaries of the island and the density-flattening regions (shown in Fig. 25) with those of the resonance region, conventionally defined as the region containing particles that satisfy the condition |∆Θ j | ≤ 2γ L (t j − t j−1 ) ∼ 4πγ L qR 0 /U eq0 ∼ 0.153 π, with γ L being the linear growth rate, and ∆Θ j , t j and t j−1 being defined in Eq. 22; this choice corresponds to include in the resonance width up to one fifth of the resonance peak.Some broadening of the resonance region is observed.Such broadening keeps the drifting gradients in the resonance region for a slightly longer time than the "linear" resonance width would be able to do, allowing for further driving the mode.This would however go quite soon to an end if the mode were not able to vary its frequency.As soon as the density-flattening region completely covered the broadened resonance region (whose width, in the present case, is smaller than the mode width), the drive would be exhausted.In the present case, however, the frequency can chirp down.This moves the resonance radius and the entire resonance region inward.The symmetry of the two large density gradients is, thus, broken: the inner one can contribute to resonant drive much more efficiently than the outer one.On the other side, wave-trapped trajectories (Fig. 17-right and Fig. 19) and the wave-trapping process (Fig. 20) will move their center around the new resonance radius, along with the island structure (cf. the black lines in Fig. 25), and the inner gradient will further drift inward, requesting further frequency variation.
We can expect that this process continues until the frequency change becomes either detrimental to the growth of the mode from the point of view of the drive/damping balance, or unable to cause a significant inward shift of the resonance radius and/or the resonance region.The latter fact effectively occurs, as seen from Fig. 27, in which all markers are reported, at four different times, in the plane (r eq , ∆Θ j ), with j referring, for each marker, at the last poloidal orbit completed.We observe that the resonance broadening (corresponding to the broadening of the red region around ∆Θ = 0) can be considered as the effect of two competing factors: the spindle-shaping of the resonant structure and its migration towards the left-bottom corner of the plot (due to the downward frequency chirping); resonance broadening is favoured by the former, contrasted by the latter, as it moves the resonance condition towards the spindle end.At a certain time the spindle-end factor prevails and the inward drift of the resonance region suddenly slows down.This is clearly shown in Fig. 28, which compares the trajectories, in the space (ω, r eq ) of the inner boundary of the density-flattening region and that of the resonance region.The time evolution of the mode frequency is also reported to make it clear how these trajectories are travelled over time.
It is possible to see that as the frequency falls below values of the order ωτ A0 0.18 (that is, at t 500τ A0 ), the inner gradient remains irretrievably outside the resonance region and any further frequency decrease is not able to prolong the relevance of the destabilising contribution yielded by the linear slice.Correspondingly, the inner power-transfer maximum loses its prevalence, as shown in Fig. 29, which plots the same boundaries as in Fig. 26, along with the radial coordinates of the time dependent maximum of the power transfer and the time evolution of the slice-integrated power.No further frequency decrease is then able to prolong the relevance of the destabilising contribution yielded by the linear slice, and the mode has to tap to a different slice in order to extract more power from the energetic particles.In the following, we examine the behaviour of the nonlinear slice; that is, the slice The time evolution of the mode frequency is also reported (green dashed line) to make it clear how these trajectories are travelled over time.As the frequency falls below values of the order ωτ A0 0.18, at t 500τ A0 , the inner gradient remains irretrievably outside the resonance region.
stage the maximum power transfer of the nonlinear slice neither occurs around the resonance radius, nor is directly related to the density gradient: it is influenced by the fluctuating field localisation, mainly determined by the interaction with the linear slice.Only later, because of resonance broadening and frequency chirping, the maximum power transfer matches a full resonance condition and appears dominated by the inner density gradient.The third relevant element shown by Fig. 33 is that the de-trapping is negligible, in the considered time interval, for such slice, so that even the outer actually-trapped particle boundary and the already-trapped particle one essentially coincide.Finally, and more important, the decoupling between density-flattening and resonance-region boundaries is not only delayed, but also less pronounced than in the linear-slice case, as shown in Fig. 32.This is consistent with the fact that the inner power-transfer maximum maintains its prevalence, as shown in Fig. 33, for a longer time and that the the power transfer yielded by the slice, after reaching its maximum, falls down at t 600.0τA0 .Marker colours follow the recipe adopted in Fig. 13.The delay in the island formation for the nonlinear slice is evident.

X. SUMMARY AND DISCUSSION
The dynamics of a single-n chirping Alfvén mode has been investigated.
We have discussed the relevance of adopting a coordinate system containing global invariants of motion in order to unambiguously identify resonant particles destabilising/stabilising the mode.In particular, this allows us to distinguish whether, in the non-linear evolution of the mode, it is driven by the same set of particles, with phase-space coordinates modified by the interaction with the fluctuating fields, or by a succession of different sets.Moreover, it makes us able to cut the phase space into slices characterised by fixed values of the global invariants; as there is no particle flux from one slice to another, the only gradients that can play the role of free-energy source for the mode destabilisation are those along the coordinates that vary inside the slice, while the gradients orthogonal to the slice do not play any role.The analysis of the nonlinear dynamics and, in particular, the saturation mechanism can be simplified by investigating separately the most relevant slices (i.e., the most effective in exchanging power with the mode).
In the gyrokinetic collisionless limit, for Alfvén modes characterised by a single toroidal mode number and a constant frequency, both the magnetic moment M and the linear combination of energy and angular momentum, defined as C = ωP φ − nE, are global invariants of  motion.In the general case and, in particular, for chirping modes (i.e., modes with varying frequency), C is no longer a constant.It is, however, possible to choose a coordinate system including other conserved quantities, with the same advantages offered by global invariants of motion.Here, in the view of a numerical approach to the investigation of mode-particle dynamics, we have proposed to adopt as constant-of-motion coordinates M and the initial values of the parallel velocity (more precisely, the initial values of the corresponding equatorial coordinate; that is, the value assumed when the particle crosses the equatorial plane at θ = 0).These coordinates are readily connected with the actual particle coordinates within a numerical approach by simply storing, in addition to the actual particle coordinates, their initial values.Such choice allows us to analyse the nonlinear evolution of the system in terms of isolated resonances.
The evolution of the particle-mode power transfer shows that a succession of different resonance takes place in driving the mode.We have identified two dominant driving phasespace slices during the linear and nonlinear phases.The analysis of the dynamics of these slices has been performed by sampling each of them by sets of test particles evolving in the fluctuating fields computed by a full-population, self-consistent simulation.We have shown that, within each set, the same bunch of particles drives or damps the mode, depending on their inward or, respectively, outward radial drift during their bounce in the potential well of the wave.Moreover, efficient power transfer to or from the mode always corresponds to the instantaneous fulfilment of the resonance condition.
The nonlinear motion of particles gives rise to two different kind of trajectories in the space (Θ, r eq ); namely, bounded orbits for particles instantaneously trapped in the wave, and unbounded ones for streaming particles.The former yield the formation of a relatively closed island-like structure around the resonance radius, characterised by mixing of particles originating from the different sides of that radius.A density-flattening region is then produced, delimited by large (negative) density gradients.In the absence of other effects, the mode would saturate as the flattening region covers the whole resonant-interaction region; that is, the intersection between the resonance region and the mode structure.In the considered case, because of the relatively small value of the growth rate and the relatively large mode structure, saturation would occur as the density gradients reach the limit of the resonance region.Frequency decreasing causes, however, the resonance region to drift inward and keeps the inner gradient within the resonant-interaction region, helped by a significant resonance broadening.The island now grows around the new resonance radius, by trapping new particles on its inner boundary, while other particles get de-trapped from the outer one.De-trapped particles still continue to contribute to the density flattening, as they have previously taken part it the mixing phenomenon.Then, while the inner density gradient is still associated to the inner boundary of the island, the position of the outer gradient is mainly related to the outer boundary of the larger region including particles that have been wave-trapped, even if they are no longer contained in the island.The inward drift of the island and its further growth cause the inner gradient to drift inward as well; to further maintain the drive capability of the slice, further frequency decreasing is needed.
The process continues as long as the frequency change is effective in causing a inward shift of the resonance region able to recapture the density gradient.
At a certain time, the combined effect of frequency chirping and resonance broadening weakens the sensitivity of the inner boundary of the resonance region to the frequency decreasing; this boundary is no longer able to lock to the inner density gradient, and the process goes to an end.After that time, the power transfer of the considered slice decreases.
To further grow, the mode has to resort to the interaction with a different resonant structure (the "nonlinear slice"), using, if needed, additional frequency variations.This slice is poorly resonant during the linear phase.It becomes fully resonant because of frequency chirping and resonance broadening.The corresponding island formation is delayed with respect to the linear slice.Moreover, its resonance broadening is more effective than that of the linear set, and it is able to prolong the permanence of the large density gradient in the resonance region.This gives the mode further drive and makes it reaching larger amplitudes than that obtained because of the linear-slice drive.
In this paper the analysis has been focused on the drive evolution during the "early" nonlinear stage, where by early we mean the stage that brings to the first fall of the destabilising contribution of the linearly dominant resonant structure (the linear slice) and to the growth and successive fall of the nonlinear-slice contribution.We have not addressed the subsequent rich and complicate evolution.Moreover, we have especially looked at the drive evolution, without investigating the nonlinear damping mechanisms.In this respect, our comprehension of the nonlinear dynamics of a chirping mode is still far from being complete, although we hope to have successfully enlightened some relevant aspects.
We have also to stress that including a constant of motion (in our case, the initial value of

FIG. 7 :FIG. 8 :
FIG.7: Coordinates, in the reduced phase space (r eq , M, U eq ), of the maximum of the power transfer density averaged over θ and φ.

FIG. 11 :FIG. 12 :
FIG. 11: Coordinates, in the reduced phase space (M, U eq0 ), of the maximum of the power transfer density averaged over θ, φ and the radial coordinate.

FIG. 15 :
FIG.15: Power transfer structures in the plane (Θ, r eq ).Each marker is coloured according to its own power-transfer rate (averaged over the last poloidal orbit).Two times are considered: t = 248.4τA0 and t = 419.4τA0 .
FIG. 16: Same as Fig. 15, but with markers coloured according to their power-transfer rate at t = 248.4τA0 .Two times are considered: t = 248.4τA0 (left) and t = 419.4τA0 (right).The left frame is, by definition, coincident with the left frame of Fig. 15.Comparing the right frame with that of the latter Figure shows that the mode is driven, at different times, by not the same particle clusters.

FIG. 20 :
FIG. 20: At each time, the indices of the radial bands of test particles owning at least one particle able to satisfy fairly well the resonance condition (last phase variation less than 0.015).A progressive asymmetric (mainly inward) trapping of new bands is observed, while no bands get fully de-trapped in this case.

FIG. 21 :Fig. 19 :
FIG.21:Only markers corresponding to particles that have already fairly satisfied the resonance condition are plotted ("already-trapped" particles).Particles that appear to be still trapped in the wave ("actually-trapped" particles) are coloured in red, while those that have cumulated a phase drift greater than a conventional threshold value (2.5π) are coloured in blu.Two different times are considered, as in Fig.19: t = 480.6τA0 (left) and t = 726τ A0 (right.).A significant amount of particle de-trapping is observed, on the outer side of the island, at t = 726τ A0 .

Fig. 25 ,
Fig.25, which shows the time evolution of the inner and outer radial boundaries of the regions including the set of actually-trapped particles and, respectively, that of the alreadytrapped particles (both actually-trapped and de-trapped particles).The radial position of the maximum negative density gradients, delimiting the density-flattening region, are also shown.The latter are reported only after a certain time, and their calculation appears to be affected by noise; the reason is that computing them requires a discretisation of the radial

FIG. 24 :
FIG. 24: Analogous to Fig. 23, but with particles, at the reference time, located just below (the blue big dot), within (violet and light blue big dots) and just above (the black big dot) the upper channel.The arrow indicates the verse of motion.Particles that enter the upper channel get de-trapped.

FIG. 25 :
FIG. 25: Radial boundaries of the already-trapped region, for the linear slice (thick red line), compared with those of the actually-trapped region (black lines).As trapping of new particles takes place on the inner side of the island structure, and de-trapping on the outer side, the inner boundary is the same for the two regions; the outer boundaries are instead different.The radial positions of the maximum negative density gradients (green circles) delimiting the density-flattening region are also shown.The boundaries of the flattened-density region are then adequately represented by those of the region including the already-trapped particles; as far as the inner boundary is concerned, however, it is well represented by that of the actually-trapped region as well; that is, by the inner boundary of the island.

FIG. 28 :
FIG.28:Trajectories, in the space (ω, r eq ) of the inner boundary of the density-flattening region (red) and that of the resonance region (blue).

FIG. 29 :
FIG.29: Same boundaries as in Fig.26, along with the radial coordinates of the time dependent maximum of the power transfer (orange dots) and the time evolution of the slice-integrated power (green line; cf Fig.10).

FIG. 31 :
FIG.31: Time evolution of the island radial width for the two slices.

FIG. 32 :
FIG.32: Same as Fig.28, but for the nonlinear slice.In this case, the decoupling between density-flattening and resonance-region boundaries is both delayed and less pronounced than in the linear-slice case.