Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

On the Eﬀect of Beating during Nonlinear Frequency Chirping∗)
Andreas BIERWAGE, Roscoe B. WHITE1) and Vinícius N. DUARTE1)
National Institutes for Quantum and Radiological Science and Technology (QST),
Rokkasho Fusion Institute, Aomori 039-3212, and Naka Fusion Institute, Ibaraki 311-0193, Japan1)
Princeton Plasma Physics Laboratory, Princeton University, Princeton, NJ 08543, USA(Received 22 March 2021 / Accepted 19 May 2021)
Spectroscopic analyses of energetic particle (EP) driven bursts of MHD fluctuations in magnetically confined plasmas often exhibit chirps that occur simultaneously in groups of two or more. While the superpositionof oscillations at multiple frequencies necessarily causes beating in the signal acquired by a localized external
probe, self-consistent hybrid simulations of chirping EP modes in a JT-60U tokamak plasma have demonstrated
the possibility of global beating, where the mode’s electromagnetic field vanishes globally between beats andreappears with opposite phase [Bierwage et al., Nucl. Fusion 57, 016036 (2017)]. This implies that there can bea single coherent field mode that oscillates at multiple frequencies simultaneously when it is resonantly driven bymultiple density waves in EP phase space. Conversely, this means that the EP density waves are mutually coupled
and interfere with each other via the jointly driven field, a mechanism ignored in some theories of chirping. In
this thesis-style treatise, we study the role of field pulsations in general and beating in particular using the Hamiltonian guiding center orbit-following code ORBIT with a reduced wave-particle interaction model in realistic
geometry. Beating is found to drive the evolution of EP phase space structures. A key mechanism is the pulsationof eﬀective phase space islands combined with the alternation of their eﬀective O- and X-points due to phasejumps between each beat. Observations: (1) Beating causes density wave fronts to advance radially in a pulsedmanner and the resulting chirps become staircase-like. (2) The pulsations facilitate convective transfer of material
between neighboring layers of phase space density waves. On the one hand, this may inhibit the early detachmentof solitary phase space vortices. On the other hand, it facilitates the accumulation of hole and clump fragments
into larger structures. (3) Long-range chirping is observed when massive holes or clumps detach and drift awayfrom the turbulent belt around the seed resonance. It is remarkable that the detached vortices remain robust and,on average, maintain their concentric nested layers while being visibly perturbed by the field’s continued beating.
c 2021 The Japan Society of Plasma Science and Nuclear Fusion ResearchKeywords: tokamak, fast particle, Hamiltonian guiding center simulation, phase space dynamics, chirping
DOI: 10.1585/pfr.16.1403087

Contents

3.3

1

Introduction

2

Physical Concepts and Open Questions
4
2.1 Phenomenology of a chirping system . . . . 4
2.2 δ f structures = ensemble of pump waves . 6
2.3 Phase slippage between pump and field . . 7
2.4 Field amplitude pulsations . . . . . . . . . 9
2.5 Characterization of chirping systems . . . . 10
2.6 Discussion of the adiabatic limit . . . . . . 11
2.7 Some open questions: Formation of solitary vortices and the role of beating . . . . . 1232Model and Methods
13
3.1 Global beating as a rationale for a semiperturbative model . . . . . . . . . . . . . 13
3.2 Equations and parameters . . . . . . . . . . 14

author’s e-mail: bierwage.andreas@qst.go.jp∗) This article is based on an invited talk at the 37th JSPF Annual Meeting(2020, Matsuyama, online meeting).

Initialization and stepping . . . . . . . . . . 16

4 Structure Formation and Saturation
4.1 Linear response, phase mixing and halos . .
4.2 First saturation with and without damping .
4.3 Potentially misleading density profiles andthe cause of pulsations and phase jumps . .
4.4 Primordial hole-clump pairs, frequency
splitting and the onset of beating . . . . . .

17
17
20

5 Analysis of the Role of Beating
5.1 Phase jumps and halo dynamics during
symmetric and asymmetric beats . . . . . .
5.2 Feedback loop for pulsed front propagation and nonadiabatic onset of chirping . . .
5.3 Growth, detachment and interference ofsolitary vortices . . . . . . . . . . . . . . .
5.4 Particle trapping and transport . . . . . . .

24

30
33

6 Summary, Discussion and Outlook

37

1403087-1

22
23

24
27

c 2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles
6.1
6.2
6.3
6.4
6.5
6.6

Volume 16, 1403087 (2021)

Halos reveal phase relations and may fuelthe growth of holes and clumps . . . . . . . 38
Feedback loop for staircase-like nonadiabatic onset of chirping . . . . . . . . . . . 38
Zonal and vortical structures . . . . . . . . 39
Robust nested vortices in the presence ofstrong beating . . . . . . . . . . . . . . . . 40
Chirp-induced resonance overlap . . . . . . 41
Concluding remarks on chirping and beating 41

A Theory of Two-Wave Beating

42

B Pulsations of Undamped Resonances
44
B.1 Spectral fluctuations . . . . . . . . . . . . . 45
B.2 Nonadiabatic boundary layer . . . . . . . . 45
C Other Interesting Observations
C.1 Splitting of spectral lines (in DMUSIC) . .
C.2 Anomalies in bounce frequency spectra . .
C.3 EP distribution during long-range chirps . .
C.4 Ghost chirps and resonance overlap . . . .

47
47
47
47
48

D Characterization and Verification of the Model
D.1 Background plasma response . . . . . . . .
D.2 Case (A’): Nonzero loading width ΔE0 . . .
D.3 Case (A”): Relaxed initial EP distribution .
D.4 Auto-optimization of resonant drive . . . .
D.5 Quiet start and sensitivity . . . . . . . . . .
D.6 Numerical convergence . . . . . . . . . . .
D.7 Calibration . . . . . . . . . . . . . . . . .

49
49
50
51
52
54
55
57

E Characterization of Diagnostics
E.1 Fourier spectra . . . . . . . . . . . . . . .
E.2 DMUSIC spectra . . . . . . . . . . . . . .
E.3 Animations . . . . . . . . . . . . . . . . .

57
57
58
59

References59

1. Introduction
In magnetically confined fusion (MCF) research,
bursting and chirping Alfvénic fluctuations have attracted
much interest since they may aﬀect the confinement offusion-born alpha particles, which are expected to provide most of the plasma heating in a power plant. Much
work has been done on this subject during the last several
decades, as can be appreciated from recent reviews and tutorials, such as Refs. [1–3]. However, there is more to belearned about these important and fascinating phenomena,
and about the underlying nonlinear self-organization processes that occur when the plasma is driven away fromequilibrium under the influence of strong sources, sinks
and dissipation. A thorough understanding can help withthe construction of reduced models that yield quantitative
predictions. It is possible to adopt analytic methods thathave been successfully applied in other research fields,such as particle acceleration and quantum field science,with adaptations made to account for the specific conditions of MCF plasmas, in particular their strong magnetization and field geometry [1, 4] (more details will follow
shortly).
A key issue in many systems is whether or not there
exists a separation of time scales between the dynamics ofthe field (growth, saturation, pulsations, chirping) and thetime scales of the dynamics of collective structures in theparticle distribution, because this determines whether adiabatic invariants can utilized. In typical cases of practicalinterest, time scale separation may be vague and exist only
in a portion of the system. The problem can be illustratedusing the cartoon in Fig. 1 (a), which shows an island inaction-angle phase space, with “radius” representing theaction variable. The island and surrounding contours appear when taking Poincaré maps of the perturbed orbitsof particles (or guiding centers) in the frame of reference
moving with the field wave. Poincaré maps of particles1
that are trapped inside the eﬀective potential well of the
field wave circulate around the island center with a characteristic bounce frequency ωb , which peaks in the center ofthe island and drops to zero at the separatrix between thetrapped and untrapped domains [5]. No matter how slowly
the amplitude A(t) and phase φ(t) of the field vary, there
will always be a boundary layer, where the value of ωb issmaller than the rate at which A(t) or φ(t) evolve.
The nonadiabatic boundary layer provides a path forparticles outside the remaining island to pass from one sideof the eﬀective resonance to the other. In the case of frequency chirping at constant amplitude, a well-defined separatrix still exists in suitable coordinates, but its domain
is reduced so that it no longer isolates the regions aboveand below the resonance [6]. As illustrated in Fig. 1 (c),this allows the remaining adiabatic island core to propagate through the surrounding phase space fluid.
The situation becomes more complicated when thefield amplitude fluctuates. For instance, consider the process of saturation of an undamped resonant instability as
shown in Fig. 1 (b). Deriving an analytical formula thatpredicts the saturation level of such an instability turned
out to be a nontrivial task. It was solved by Dewar [7], who
obtained a fairly accurate estimate by examining two idealized extremes as lower and upper limits. The problem that
Dewar encountered was that the system is actually fully
nonadiabatic during the saturation process.
The same is true for the spontaneous onset of nonlinear frequency chirping as observed in simulations of instabilities that contain a coherent energy loss channel. Theenergy loss may occur in the field (damping) or in the particles (drag), its role being to permit phase slippage between
field waves and density waves, as we will discuss in moredetail later. An analytic formula for the nonadiabatic onset of chirping in tokamak geometry has been derived by
1 Hereafter, when we speak loosely of “particle motion in/around aphase space island”, we mean the motion of the particle’s Poincaré map.

1403087-2

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 1 Schematic illustration of processes in resonantly interacting systems; here, energetic particles (EP) with shear Alfvén waves. The
upper diagrams show time traces of the amplitude A(t) of a symmetry-breaking field perturbation, whose time-dependence has theform A(t) sin(−ω0 t − φ(t)), with initial frequency ω0 and dynamic phase φ(t). In columns (a)-(c), only a single mode is considered,giving rise to resonant phase space islands whose internal contours are drawn as gray ellipses. The background EP density
decreases with increasing radius. Red and blue arrows indicate uphill and downhill convection, respectively. Here, “radius” is a
proxy for a generalized action variable in suitable action-angle coordinates. Column (a) illustrates a steady state; i.e., a mode with
fixed frequency A and phase φ. The closed separatrix of the perfectly adiabatic steady state (a) becomes a nonadiabatic boundary
layer (orange band) in the semi-adiabatic cases (b) and (c) due to a pulsating amplitude A(t) and/or an accelerating phase φ(t) (i.e.,
chirping). Faster phase adjustments and larger fluctuation amplitudes yield a wider boundary layer around the adiabatic core, until
the entire system becomes nonadiabatic. Column (d) shows an example of convective amplification, where three radially adjacentmodes labeled k = 1, 2, 3 successively rise to large amplitudes and “relay” a resonant particle bunch radially outward [8].

Zonca & Chen [8] based on the concept of a convective
instability as illustrated in Fig. 1 (d). The chirping rate ispredicted to follow the curve of a hyperbolic tangent; i.e.,
exponential acceleration followed by exponential deceleration. We will confirm that this prediction compares wellwith simulations during the stage of spontaneous formationand separation of a nonadiabatic hole-clump wave pair.
Rapid short-lived chirps that last only a millisecondor less are often observed experimentally in well-resolved
spectrograms (e.g., Fig. 6 of Ref. [9]). These chirps may
remain fully nonadiabatic until their end, but they may also
spawn vortices with an adiabatic core. Suﬃciently massive
vortices may gain a large degree of control over the field
and become self-sustained for several milliseconds. Theirpropagation can be observed as long-range chirps and they
may cause ballistic long-range transport. We will analyze some concrete examples where semi-adiabatic vortices form after a few milliseconds of total nonadiabaticity.
Berk, Breizman and Petviashvili [10] have derived
a widely used formula that predicts the rate of vortexmediated chirps in regimes near marginal stability; i.e.,
when strong damping is overcome by a slightly stronger
drive. The predicted overall range and rate of chirping often lies in the right “ball park”, even in cases with fairly
rapid chirps and even relatively far from marginality. However, the actual form of long-range chirps in a realistic setting often deviates from the predicted square-root time dependence. This is to be expected because the assumptions

made are, at best, satisfied only at intermediate times, after the nonadiabatic layer has narrowed and before background nonuniformities and collisions become important.
A rigorous treatment of the nonadiabatic processes
remains challenging and usually requires numerical simulations, but analytical theory can help to gain insightand minimize the computational eﬀort by reducing theMaxwell-Vlasov system. A comprehensive theoreticalframework has been formulated by Chen & Zonca [1] on
the basis of concepts from quantum field theory, invoking a Schrödinger-like equation for the intensity evolution (bursts or solitons) and a Dyson-like equation for(driven or spontaneous) emission of (zonal or vortical)
phase space structures that break the symmetry of the reference state [4]. For instance, this theory captures nonadiabatic avalanches in a process called convective amplification, as illustrated in Fig. 1 (d). The Schrödinger-Dyson
system of equations is less complex than the original gyrokinetic Maxwell-Vlasov system, but handy formulas forthe quantitative prediction of chirping dynamics are not yet
easy to come by, except in certain limits, such as shortwavelength modes interacting with magnetically deeply
trapped particles [8], or in certain space plasmas [11].
It is particularly diﬃcult to reduce analytically the dynamics of long-wavelength modes in a MCF plasma, because the plasma profiles and the global geometry of the
field can have a significant influence on the dynamics. In
this case, it can be diﬃcult to justify the assumption of

1403087-3

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

the existence of two disparate spatial scales that underlies
some of the reductions used by Chen & Zonca. However,the long-wavelength regime is of practical interest, becausesuch modes are common in present-day MCF experiments
and the associated transport tends to be global.
The present work was motivated by a study of longwavelength chirping modes with toroidal mode numbers
n = 1, 2, 3 driven by energetic particles (EP) in JT-60U
tokamak plasmas [12], where experimental and numericalevidence for beating in the field signal due to interference
between multiple chirps was reported. Using the ORBIT
code with a reduced model for the fluctuating field [13], we
investigate here the physical implications of beating. We
examine the wave-particle interactions in realistic tokamak
geometry, but use a somewhat simpler scenario than theoriginal JT-60U case that motivated this work, where theprimary resonances were located near the magnetic axis,
touching the domain of stagnation orbits. Here, we avoid
such boundary eﬀects by working with a resonance locatedat mid-radius in the domain of circulating particles. Our
plasma has the dimensions of the conceptual fusion reactor
FIRE, but a weaker magnetic field with a flux density of
B0 = 0.49 T at the axis located at R0 = 2.15 m.
The eﬀect of field pulsations in general and beatingin particular is relevant for several aspects of nonlinear frequency chirping, including the nonadiabatic onset of chirping, hole-clump pair formation, and the applicability of reduced transport models, such as the so-called bucket transport [6] and waterbag [14–16]. In Sec. 2, we introduce essential physical concepts, review methods used to characterize chirping systems, and discuss open questions that
may be addressed by elucidating the eﬀects of field pulsations and beating. In Sec. 3, we describe the physicalmodel and numerical methods used in this study. The processes of phase space structure formation, nonlinear saturation2 and the onset of beating are revisited in Sec. 4. The
eﬀects of beating during the nonadiabatic onset of chirping, the formation and propagation of massive vortex structures and the particle transport they cause are analyzed inSec. 5.
This study makes extensive use of comparisons between cases near and far from marginal stability. Comparisons are also made with simulations where no sustained
chirping occurs due to the absence of field damping. The
contrast between these examples will help us to highlightthe trends. By presenting known and less known eﬀects
in a unified systematic treatise, we hope to contribute to abetter understanding of nonlinear frequency chirping.
Sections 2, 4 and 5 can be regarded as three shorter
articles that have been compiled here into a unified thesisstyle treatise. We conclude with a summary, discussion
2 “Saturation” is used here as a generic term, referring to a stagnationin the growth of the field amplitude, such as a (quasi-)steady state due toa permanent exhaustion of resonant drive, or temporary peaking due to atransient dynamic balance between resonant drive and damping comingfrom diﬀerent regions of phase space that we integrate over in Eq. (2).

and outlook in Sec. 6. The comprehensive Appendicescontain a mathematical theory of two-wave beating, additional simulation results, and a thorough characterizationof the model and methods used, including tests for numerical convergence and sensitivity. At the end of the paper,
there are links to animated movies for several simulations.

2. Physical Concepts and Open Questions
In this section, we establish the conceptual frameworkthat we will use in our discussion and interpretation of thesimulation results. The rationales behind the chosen physical pictures are explained and useful quantities characterizing chirping systems are defined. Finally, we discuss why
we expect that studying the eﬀects of beating may improve
our understanding of chirping systems. Before we discuss
the theory, let us show a concrete example of nonlinear
chirps as seen in our simulations.

2.1 Phenomenology of a chirping system
Figure 2 shows results from one of the simulations that
we will analyze in detail in this paper. Appendix E.3 contains a link to the movie file. The simulations use a reduced
model of EP-driven ideal magnetohydrodynamic (MHD)
fluctuations in the frequency band of shear Alfvén waves
in a tokamak plasma. Our EPs are fast deuterons with ki-

Fig. 2 Evolution of ideal MHD fluctuations in one of our numerical simulations of nonlinear frequency chirping using the
ORBIT code [13, 17, 18] that we analyze in detail in thispaper. Panel (a) shows the raw signal s(t) of Eq. (1) incode units, and panel (b) shows the Fourier spectrogram
obtained with a sliding time window of size 0.47 ms. Thiscase is marginally unstable and develops semi-adiabatic
vortices. The dotted parabola is the theoretically predicted chirp front assuming marginality, adiabaticity, anduniform background [10]. The “ghost chirp” seems to be
an artifact and is discussed in Appendix C.4.

1403087-4

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

netic energies around 80 keV. Figure 2 (a) shows the raw
signals(t) = A(t) sin(−ω0 t − φ(t)),(1)

which represents ideal MHD displacements ξ. In thepresent case, the peak values of this signal fluctuate atlevels below 10−3 , which corresponds to radial displacements of resonantly trapped particles by δrres /R0  6×10−3
and ideal MHD bulk displacements δrmhd /R0  7 × 10−4
(cf. Appendix D.7). Panel (b) shows the corresponding
Fourier spectrogram obtained with a sliding time windowof size Δtwin = 0.47 ms. The chirps in Fig. 2 (b) proceed at rates on the order of φ̈ = δν̇  10 kHz/1 ms(φ̇ ≡ dφ/dt) around the seed resonance with initial frequency ν0 = ω0 /(2π) = 100 kHz. In the present example,
the chirping dynamics enhance the displacement of resonant and near-resonant particles by an order of magnitude
to δrchirp /R0 ≈ 4 × 10−2 , which illustrates why this phenomenon is of practical interest.
The example in Fig. 2 was simulated in the perfectly
collisionless limit and without external EP sources, so thechirps tend to continue and dominate the signal for a longtime, whereas the original mode is not revived. In reality,
rapid nonlinear chirps may be visible for less than a millisecond, and even long-range nonlinear chirps usually last
no longer than 10 ms, so this phenomenon can be taken to
be insensitive to changes in the background plasma conditions (magnetic geometry, kinetic profiles, rotation) on
global scales, which typically evolve on the time scale of
100 ms or more under the influence of heating, fueling, collisional diﬀusion and turbulent transport.
The δ f -weighted kinetic Poincaré plots in Fig. 3 area way to visualize the structure and motion of the incompressible Vlasov fluid that makes up the EP phase space.
Density waves in this incompressible fluid are realized by
the convective interchange of denser and lighter fluid elements that exist due to an initially imposed density gradient. Here, the gradient is along the canonical toroidal
angular momentum Pζ , which represents the radial direction along the vertical axis in Fig. 3 (a). EP phase spacestructures causing upward chirps travel upward in Fig. 3and vice versa. The vertical axis of Fig. 3 (b) has been converted to frequency units for easier comparability with thespectrogram in Fig. 2 (b).
The Poincaré section in Fig. 3 moves with the phasevelocity of the seed wave along the toroidal angle ζ, so δ fstructures at the seed resonance appear stationary. Shadesof blue indicate a reduction in EP phase space density(δ f < 0, “hole”) where lighter fluid has moved into regions
previously occupied by denser fluid, and shades of red indicate increased density (δ f > 0, “clump”). In the present
scenario, four cross-sections of the same phase space structure appear at diﬀerent poloidal angles ϑ since the underlying resonance has p = 4 elliptic points along ϑ. Notethat we use the terms “hole” and “clump” in a more general manner than some readers may be used to. Namely,

Fig. 3 Two snapshots of EP phase space density waves around aseed resonance with p = 4 poloidal elliptic points foundin the case of Fig. 2. The structures are visualized using
δ f -weighted kinetic Poincaré plots, with shades of blue
representing decreased density (δ f < 0, “hole”) and red
increased density (δ f > 0, “clump”). The Poincaré maps
are accumulated during a short interval of one toroidaltransit, so, strictly speaking, we plot a time-integrated
phase space density perturbation δ f transit . The canonicaltoroidal angular momentum Pζ serves as a radial coordinate, here normalized by the poloidal magnetic flux at theplasma edge ΨP,edge . The initial EP density (not shown)decreases from the bottom (inner plasma) to the top (outerplasma). In the rotating frame of reference chosen here,
fluctuations at the seed resonance (P̂ζ,res = 0.719) appear stationary, and the poloidal phase velocities increaseaway from the resonance as indicated by the horizontal
arrows. Using the calibration in Eq. (D6), the vertical axisof panel (b) has been converted to frequency ν for easiercomparability with the spectrogram in Fig. 2 (b).

we use them not only to refer to solitary vortices but also
for short-lived convective perturbations in EP density.
Some notable dynamic structures seen in the system
studied here have been labeled in Figs. 2 and 3. Chirping
starts with the radial propagation of convective EP density“wave fronts”. In the wake of these wave fronts, a turbulent
belt forms, which consists of more or less stratified layersof diﬀerentially rotating hole and clump waves. The interactions between these layers are reminiscent of sheared“convective interchange” and “wave breaking”: wave tops
rise, overtake wave bottoms, then collapse back down.
Sometimes, we observe the complete detachment of a massive hole or clump from the turbulent belt. Such structures
can be seen in the upper and lower parts of Fig. 3 (b) andare referred to here as “solitary vortices”. The massive solitary clump vortex, whose interior looks like a spiral galaxy,is responsible for the long-range upward chirp in Fig. 2 (b)that spurs away from the turbulent domain up to 180 kHz.
The propagating phase space structures and associated chirps may be aﬀected by but are not necessarily constrained to the linear MHD response spectra of the background plasma, which (in the limit t → ∞) consist of singular continua or discrete global eigenmodes. In fact, theMHD spectra are entirely ignored in the reduced model

1403087-5

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

that we employ here — boldly assuming that the plasma
responds equally well at any frequency, being subject only
to an ad hoc damping rate γd that is spatially uniform,
independent of frequency and constant in time. Our reduced model captures the essential physical mechanisms
for chirping with full complexity in the particle dynamics
and reduced complexity in the field dynamics, which enhances physical transparency and computational eﬃciency.

2.2

δ f structures = ensemble of pump waves

In the previous section, we have, in a casual way, related the coherent structures appearing in Poincaré plotsof δ f in Fig. 3 to the chirps seen in the spectrogram in
Fig. 2 (b). Let us review the underlying theoretical foundation.
When applying a symmetry-breaking electromagnetic
perturbation that displaces our energetic charged particles
in the direction of a preexisting gradient in their initial
guiding center distribution F0 = F(t = 0), one obtains
a perturbed distribution F(t) = F0 + δ f (t), whose fluctuating part δ f (t) consists of wave-like structures. Shortwavelength portions of δ f that cease to contribute to the
dynamics of interest act like a zonal (angle-independent)
structure, which we denote by δ f¯ and interpret as a modification of the axisymmetric reference state F0 . During
the ensuing mutual interactions between the phase spacedensity waves and the electromagnetic field, remaining or
newly formed gradients in F̄(t) = F0 + δ f¯(t) may facilitate
the excitation of additional waves or continue to feed existing ones. Although conventionally referred to as “waveparticle interactions” — which is, of course, correct fromthe energetic point of view — we interpret these dynamics
in terms of interactions between two types of waves:• electromagnetic field waves, here represented by E⊥ ,and• collective density waves represented by δ f (t).
Speaking of density waves instead of particles is not just asemantic matter. It is motivated by Maxwell’s equations,which state that the field responds to collective structuresin δ f , whose group velocities generally diﬀer from particle
velocities. Let us elucidate this point within the scope ofthe reduced field model that we will be using in this work.
We express the combined time-dependence of theideal MHD field waves as E⊥ (t) ∝ A(t) cos(−ω0 t − φ(t)),
with an amplitude factor A(t) and nonlinear phase shift φ(t)relative to a reference wave with frequency ω0 .3 Our reduced equations governing the field response have the form
ν2dA
= −Ωc A02
dt
ω0
A



ν2
dφ
= −Ωc A03dt
ω0

d5 Zgc δ f
gc · E⊥− γd A,A(2a)

gc · ∂t E⊥
;
A

(2b)

d5 Zgc δ f

3 A(t) cos(Θ(t)) corresponds to the real part of the analytic signal s̃(t) =s(t) + i ŝ(t) = A(t)eiΘ(t) , where ŝ(t) is the Hilbert transform of s(t).

(see Sec. 6.9 of Ref. [19]) where νA0 = A0 /R0 is the onaxis Alfvén frequency, Ωc the cyclotron frequency, d5 Zgc
is a volume element in guiding center (GC) phase space,gc is the GC velocity, and ∂t ≡ ∂/∂t in Eq. (2b) acts
only on the “fast” oscillations ω0 . All damping mechanisms unrelated to the particle distribution F0 are captured by γd , which we take to be constant. The Jacobian
has been absorbed in the distribution function F, so
d5 Zgc F = N is the number of particles (whose conservation implies d5 Zgc δ f = 0). The set of GC coordinates Z gc = (Pζ , K, μ, ϑ, ζ) consists of the canonical toroidal angular momentum Pζ , kinetic energy K, themagnetic moment μ, and poloidal and toroidal angles ϑand ζ. For particles with mass M, electric charge Qeand velocity , in a toroidally axisymmetric magnetic field
B = ∇ζ ×∇ΨP +Bζ ∇ζ, these GC coordinates have the formPζ = −Ψp + Bζ


,
Ωc

μ=

M2⊥
,
2B

K=

M2
,
2

(3)

where 2πΨP is the poloidal magnetic flux,  is the GC
velocity parallel to the magnetic field, 2⊥ ≡ 2 − 2 , andΩc = QeB/M with e the electric charge of a positron.
Note that the phase space density gradients ∂ Z gc F do
not appear in the field Eq. (2) explicitly; only the density
perturbation δ f (Z gc ) appears. The Vlasov equation∂t δ f = − Ż gc · ∂ Z gc F
= − (ϑ̇∂ϑ δ f + ζ̇∂ζ δ f ) − (Ṗζ ∂Pζ F + K̇∂K F) , 
free streaming

acceleration(4)
states that the instantaneous rate of change of δ f is determined by phase space flows Ż gc ≡ dZ gc /dt across phasespace density gradients ∂ Z gc F. Meanwhile, the instantaneous structure of δ f reflects the time history of these dynamics (until erased by external sources, sinks and collisions). Thus, like the fields in Maxwell’s equations, ourfield governed by Eq. (2) sees only the result of that history in the form of an evolving density perturbation δ f (t).
For instance, let us assume that the zonal gradients
have been flattened, so that F̄ ≈ const., and what remains
is a single density wave with frequency ωpump . The absence of zonal gradients means that there is no net guiding center current δ f gc along the transverse electric fieldE⊥ , so the total energy transfer d5 Zgc δ f gc · E⊥ vanishes, because there is no “DC” component and the “AC”
components δ f ∝ sin(−ωpump t) and E⊥ ∝ cos(−ω0 t − φ)are out of phase. Ignoring damping, Eq. (2a) becomesȦ = 0. Meanwhile, the time derivative on the righthand side of the phase Eq. (2b) aligns the phases of density and field fluctuations. Using the fact that the integral
over ∂t (δ f gc · E⊥ ) vanishes, we can let δ f gc · ∂t E⊥ →12 (δ f gc · ∂t E⊥ − gc · E⊥ ∂t δ f ), so that Eq. (2b) becomesφ̇ = ωpump − ω0 . This means that the field wave respondswith oscillations E⊥ ∝ cos(−ωpump t) of the same frequency
as the density wave, which acts as a pump wave as it propagates subject to the free streaming terms in Eq. (4). In the1403087-6

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

presence of zonal gradients, the pump waves propagate up
or down those gradients and grow or decay in magnitude.
Therefore, Eq. (2) states that the field responds to thepumping action exerted by an ensemble of density waves,
so it behaves as a driven oscillator; albeit an active one,that is capable of shaping the δ f landscape4 . One may saythat the integral on the right-hand side of Eq. (2a) measuresthe correlation between the structure of the density waves’
current density Qeδ f gc and the field wave E⊥ , so the response is stronger the better the phases match. Conversely,
the coherence of density waves is maintained by resonantinteractions with the field, where “resonance” means thatphase-matching is maintained over many field oscillation
cycles. Thus, the eﬀect of resonance emerges on intermediate time scales situated between rapid field oscillations
and transport, allowing coherent δ f structures to form andto remain robust.
The total field signal E⊥ (t) cannot be phase-locked tomultiple pump waves simultaneously, and as the relative
phases shift, the sign of the phase correlation integral in
Eq. (2a) can change5 . The beating phenomenon, whose effect we study in this work, can be interpreted as the temporal interference pattern of two or more pump waves inphase space density, which propagate with diﬀerent phasevelocities around the plasma torus. The math for the special case of two-wave beating is presented in Appendix A.
Variations in the pump wave distribution are observable in spectrograms of the field signal s(t). The spectral
patterns one sees depend on the time window used for the
spectral analysis. Short time windows can reveal details
about rapid phase adjustments performed by the field, but
may not reflect the existence of multiple pump waves, only
their combined instantaneous action. Long time windows
show only the overall trend, but allow to discern a larger
number of spectral harmonics that can then be attributed tomultiple pump waves as in Figs. 2 and 3. We will analyze
dynamics in our simulations on short and long time scales.
When the EP density waves propagate radially in atokamak plasma, their transit frequencies change, giving
rise to the chirps in the spectrogram of the field waves they
drive. We will consider scenarios with a destabilizing “radial” gradient ∂Pζ F̄ > 0 and will observe “clump waves”(δ f > 0) that tend to propagate downhill and “hole waves”(δ f < 0) that tend to propagate uphill. The radial propagation of these phase space structures is facilitated by fielddamping γd and the reason for this is considered next.
4 Parallels with the concept of parametric instabilities are evident.
5 For instance, even when the density waves δ f are still growing, thefield may saturate or even vanish if the phase relation is such that thecontributions of pump waves with positive and negative δ f cancel in thecorrelation integral of Eq. (2a). This illustrates the diﬀerence between
the linear regime — where phases are spontaneously aligned such that∂t ln A = ∂t ln δ f regardless of amplitude and position — and the nonlinear regime, where the phase relations can become arbitrarily complicatedand one has to consider the flows Ż gc of GC phase space fluid acrossphase space density gradients ∂ Z gc F. The local flow contours, in turn, aredetermined by the field’s instantaneous phase and amplitude.

2.3 Phase slippage between pump and field
Interactions between pump waves in EP phase spacedensity and electromagnetic field waves do not necessarily
lead to spontaneous chirping. There needs to be a mechanism that causes phase slippage between these density
waves and field waves, which we review in this section.
When the field wave is unaﬀected by the presenceof the EP density perturbation δ f (e.g., when it is externally enforced), one obtains the situation shown in Fig. 4.
Around the resonance (vertical dash-dotted line), an island structure forms in phase space, which can be visualized using a kinetic Poincaré plot as in panel (a), using
action-angle coordinates like (Pζ , ϑ) in the frame of reference moving with the phase velocity of the wave field.
The circulation of particles that are trapped inside the island causes a perturbation in phase space density δ f with acoarse-scale structure as shown in (b), with a hole (δ f < 0,blue) on one side and a clump (δ f > 0, red) on the other.
The hole and clump in this configuration are stationary collective modulations of the particle density, not tied
to the motion of individual particles. Asymptotically, thesystem develops a non-axisymmetric EP distribution that

Fig. 4 Perturbation of the EP trajectories and EP density in thepresence of an imposed oscillating field with fixed amplitude and phase velocity. The kinetic Poincaré plot inpanel (a) shows a phase space island around the resonance at P̂ζ,res = 0.719. The colors in (a) represent the instantaneous total (kinetic + potential) energy E of a particle. Panel (b) shows the relaxed phase space density perturbation δ f , averaged over the last 50 transits. Panel(c) shows the radial profile of the axisymmetric EP reference distribution F0 (shaded triangle) and three views of
the perturbed EP distribution F = F0 + δ f : poloidally
averaged, O-point profile and X-point profile.

1403087-7

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 5 Schematic illustration of the nonlinear frequency chirping process in a monoenergetic beam that is abruptly perturbed out of
equilibrium. Note that, in reality, stages 1-4 are not as distinct as they are shown here. For instance, phase slippage (item 3) occurs
before the wave amplitude saturates and decays; in fact, we will show that, in a marginally unstable system, this phase slippage
and the resulting beats (item 5a) contribute actively to an early temporary saturation and decay, and subsequent revival of the fieldwave. In the case of EPs interacting with shear Alfvén waves in a toroidal plasma, acceleration (deceleration) of particles andresonant damping (drive) of the field is connected with radially inward (outward) transport.

is consistent with the perturbed field6 . If the amplitude of
the symmetry-breaking field perturbation is reduced, someparticles are released from the island and carry with them
the previously established density modulation. One willthen observe a counter-propagating hole-clump wave pair.
When the dynamics of the field’s amplitude A(t) and phaseφ(t) are coupled to the dynamics of the density waves δ f asin Eq. (2), the propagating hole and clump structures will
act as pump waves that drive field fluctuations, which respond auto-resonantly to the pump frequencies.
From this we infer that instability is not essential for
triggering a frequency shift. Although a gradient in the reference state F0 is required to produce density waves in an
incompressible fluid, that initial gradient by itself does notneed to be destabilizing. The energy required to launch oﬀresonant pump waves can also be provided by an abruptly
imposed perturbation that breaks the symmetry of the par6 It must be noted that Fig. 4 does not show an equilibrium state. We
have fixed the amplitude of the field perturbation at an arbitrary value
A0 = 10−6 that is far below the level of about 3 × 10−5 required forthe formation of a (quasi-)steady state that is consistent with the initial
EP density gradient (cf. Fig. 14). The non-equilibrium character of thisrelaxed EP distribution can be inferred from the lack of mirror symmetryof δ f in Fig. 4 (b), which implies that there is still potential for net energy
transfer (which we have inhibited here by freezing the amplitude). Thus,
the flat O-point profile in (c) is misleading! There are still active gradients
nearby. The X-point profile is also misleading; the actual direction of the
X-point gradient is diagonal in (P̂ζ , ϑ) and the energy stored there is noteasily tapped (cf. Sec. 4.3). The result of launching a simulation from thisrelaxed state is shown in Fig. D4 of Appendix D.3. Since a primordial
hole-clump pair has already been set up, this setup eﬀectively skips the
exponential growth phase and first beat.

ticle distribution. What is essential for nonlinear chirpingis a mechanism that facilitates phase slippage by preventing complete trapping of displaced particles; for instance,a mechanism that damps the field wave.
The process leading to the onset of chirping is illustrated schematically by the cartoon in Fig. 5. As shownon the left-hand side, the process begins with symmetry
breaking and (partial) particle trapping in the potential wellof the field wave:
0. Consider a monoenergetic bunch of particles circulating around the plasma torus at a characteristic transit frequency Ω0 . The particle bunch is abruptly perturbed by a resonant field wave that causes an abrupt
finite displacement and oscillates at the same frequency ω0 = Ω0 . This perturbation breaks the symmetry of the system and takes it out of equilibrium.
1. The collective response of the guiding centers produces a modulation δ f in their distribution function
F = F0 + δ f in position and velocity space. Depending on the local phase of the field wave, some particles
are accelerated, others decelerated.
When the field and the particles constitute an energyconserving closed system, the particles oscillate indefinitely inside the field’s potential well. A high-resolution
spectrogram may show transient frequency splitting, but
no sustained chirping occurs.
The situation changes when this field-particle system
is not isolated and exchanges energy with other compo-

1403087-8

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

nents of the plasma or its surroundings. On the right-handside of Fig. 5, we illustrate how sustained frequency chirping can occur when
2. the field wave experiences some form of energy loss
that is not related to the resonance at hand and exceeds
a certain threshold. Possible damping mechanisms
include phase mixing in the wave-carrying medium(a.k.a. Landau damping), nonlinear wave-wave coupling, or viscous dissipation. When this happens, theparticles are detrapped, and . . .
3. . . . the density waves disperse. For simplicity, thiscartoon ignores the continuity of the velocity distribution and shows only the density wave fronts constituted by maximally decelerated and maximally accelerated particles, whose characteristic transit frequencies satisfy Ω− < Ω0 < Ω+ .7 The consequence of detrapping is phase slippage between the density waves
and the (weakened) field wave.
4. The released density waves excite auto-resonantly
phase-matched field waves. In the present symmetriccase, the result is a pair of new field waves, one downshifted (ω− ) and one up-shifted (ω+ ) in frequency.
As indicated by the triangle of green arrows in the centerof Fig. 5, this cycle of trapping, detrapping, phase slippage
and auto-resonance can repeat, spawning multiple generations of chirps that successively appear as time passes by.
In reality, steps 1-4 illustrated simplistically in Fig. 5
are not so distinct and may overlap in time. For instance,
the field is usually damped continuously, so the phase slippage described in step 3 may already begin before the field
amplitude saturates and decays. The essential point is that,
because of damping, the field wave has irreversibly lost apart of the energy that it had received when disturbing thegradients in the particle’s phase space density, so it does
not possess enough energy to reverse all the displacementsthat have occurred in the particle distribution. In otherwords, the (weakened) field is able to trap only a fraction of
the particles that have been displaced, and that fraction decreases with increasing damping rate. The particles that thefield fails to pull back across the resonance will carry along
density modulations that act as pump waves at shifted frequencies ω+ and ω− .
The same principles apply in more complex systems.
The consequences of irreversible field energy loss is particularly manifested in (but not limited to) marginally unstable scenarios subject to strong drive and nearly equallystrong damping, like our examples in Figs. 2 and 3. In
such cases, the onset of chirping occurs early and proceeds
gradually. In contrast, in cases with strong drive but weakdamping, chirping tends to be retarded and begin abruptly.

7 Like a two-stream approximation [20] of the continuous system [21].

2.4 Field amplitude pulsations
In order to characterize the evolution of the field amplitude, we define the instantaneous growth rate γ(t) as
A(t) = A0 eγt ,γ = d ln A/dt = γk − γd ;(5)where γd represents some form of damping (here taken
to be constant), and γk (t) represents the instantaneous kinetic drive; namely, the integral on the right-hand side ofEq. (2a) divided by A(t). At the beginning of a simulation
with a suﬃciently small initial perturbation A0 and an axisymmetric EP distribution F0 ,8 the kinetic drive has an approximately constant value that corresponds to the growthrate of the linearized system: γk ≈ γL = const. The linear
growth rate γL (∂Pζ F0 , ∂K F0 ) is determined by the gradientsof the reference state F0 (Pζ , K, μ) at the seed resonance.
When γL > γd , the field amplitude A(t) begins to grow(nearly) exponentially with a rate γ ≈ γL − γd .
Figure 6 shows two examples of the amplitude evolution in ORBIT simulations in the scenario that we study inthis paper. Panel (a) shows an undamped case (γd = 0)
with relatively weak drive γL /ω0 = 0.5%. The rapid
oscillation of the field signal s(t) with frequency ω0 =2π × 100 kHz are not visible, since we plot only its temporal envelope A(t). One can see that the amplitude grows
to about A ≈ 2.7 × 10−5 . The period τpulse  1 ms of thesubsequent pulsations seen in Fig. 6 (a) is more or less con-Fig. 6 (a) Weak pulsations and (b) strong beating in the amplitude A(t) associated with (a) an undamped and (b) astrongly damped electromagnetic field wave destabilized
by resonant EPs. Both cases have the same small initial
growth rate γ/ω0 = (γL − γd )/ω0 ≈ 0.5% (γ/ν0 ≈ 0.03)
with ν0 = ω0 /(2π) = 100 kHz. The damped case in (b)exhibits strong chirping as was shown in Fig. 2. Panel (b)is split in two parts with diﬀerent scales (factor 5) in orderto show the small amplitudes during the first 2 ms.

8 As shown in Apendix D.3, exponential growth and the first beat areeﬀectively skipped if one starts with a relaxed non-axisymmeteric EP distribution that has adapted to the initial perturbation of the field as in Fig. 4.

1403087-9

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

sistent with Dewar’s theoretical prediction [7]previous works dealing with nonlinear frequency chirping— theoretical, numerical and experimental. Clarifying therole of these pulsations in diﬀerent stages of a chirping
simulation is the main subject of the present study.

ωDewar=b

256
γL ≈ 2.88 × γL ,
9π2

(6)

which gives a bounce period τDewar= 2π/ωDewar≈ 0.7 msbb
for our parameters.
Pulsations like those appearing in Fig. 6 (a) are sometimes said to be caused by the continued over- and undershoots of the gradients around the resonance that, in turn,are caused by the bouncing motion of near-resonant particles that are trapped inside the eﬀective potential wellof the field wave. Although this interpretation may appeal to intuition, it is important to keep in mind that pulsations can be caused not only by the acceleration terms(Ṗζ ∂Pζ + K̇∂K )F of the Vlasov Eq. (4), but also by the freestreaming terms (ϑ̇∂ϑ + ζ̇∂ζ )δ f . The details of these interactions are complicated, which is why we use numericalsimulations. As suggested in Sec. 2.2, it may be more instructive to think of these amplitude pulsations as being
a consequence of the motion of density waves, because
it is δ f (not Ż gc · ∂ Z gc F) that appears in the field Eq. (2).
The field pulsations cause a repeated trapping and detrapping of density waves in the boundary layer of the phasespace island as illustrated schematically in Fig. 1 (b). The
resulting variations in the relative phase between densityand field waves, in turn, feed back on the evolution of the
field via the phase correlation integral between EP current density and electric field in Eq. (2). We suspect thatthese boundary layer dynamics, in combination with radial asymmetries in our simulation setup, are the primarycause for the sustained and somewhat irregular oscillations
seen in Fig. 6 (a). The fact that the bounce frequencies varyacross the island radius [5] may also play a role. Limited
numerical accuracy can cause spurious growth of the field
amplitude (see Appendix D.6). We will analyze these data
in more detail in Appendix B, examining the eﬀect of thepulsations and demonstrating that, although there is no sustained chirping, the field’s oscillation frequency can varymeasurably even in the absence of field damping.
Figure 6 (b) shows that a very diﬀerent evolution withstrong pulsations of the field amplitude A(t) is seen in thepresence of strong damping, which is here nearly equal tothe drive. This is the same case as in Figs. 2 and 3, with
γd /ω0 = 7.5%. The damping is marginally overcome bya steep gradient ∂F0 /∂Pζ > 0 that gives a linear kinetic
drive γL /ω0 ≈ 8%. Thus, the overall growth rate during
the exponential phase is small, γ/ω0 ≈ 0.5%, and equal
to that in the undamped example in Fig. 6 (a). We observein Fig. 6 (b) an early saturation of the field amplitude during the first millisecond, at a still tiny amplitude of onlyA ≈ 6 × 10−6 , not far from the initial value A0 = 10−6 .
Subsequently, the amplitude grows by nearly two orders ofmagnitude to about 3 × 10−3 , while performing pulsations
with a magnitude of 100% on time scales of 0.1 ms or less,not much longer than the 10 μs seed wave period.
These amplitude pulsations have been shown in many

2.5 Characterization of chirping systems
Aside from geometric eﬀects and other nonuniformities of the plasma, the behavior of our reduced system islargely determined by its initial proximity to the state of
marginal stability, where drive and damping approximatelybalance,|γL − γd |/γL  1 (marginal (in)stability),(7)

and by the amount of free energy stored in the initial gradient, which can be measured by the ratio of γL to the initialwave frequency ν0 at the seed resonance9 ,γL /ν0  1(weak drive),(8)

γL /ν0  1 (strong drive).

(9)

The cases discussed in this paper are situated in thestrongly driven regime, γL /ν0 ≈ 0.5 (γL /ω0 ≈ 8%). We
consider one case close to marginal stability, γL /γd ≈ 1.07,
and one case with relatively rapid growth, γL /γd ≈ 1.88.
An important measure that characterizes the phasespace dynamics underlying the chirps is the degree of adiabaticity of the dynamics, which is determined by the ratio of the mean and maximal displacement of a simulation
particle (representing an EP Vlasov fluid element) during
a field pulsation period τpulse ,10δPζ pulseδEpulse(adiabatic∼< 1 regime); (10)max|δPζ |pulse max|δE|pulse
where δPζ
= Pζ − Pζ pulse and . . .pulse
≡
dτ
(.
.
.)
denotesthepulse
average.
E
=
K
+QeΦ
τ−1
pulse
pulse
is the total particle energy, consisting of kinetic energy K
and electrostatic potential energy QeΦ. In simple terms,
Eq. (10) means that, in the adiabatic regime, the coherent
phase space structures (if any) rotate faster than they advance, so the chirps are associated with the drift of vortexlike structures in phase space, which we have labeled “solitary hole/clump vortices” in Fig. 3 (b). In this case, thepulse period τpulse may be identified with the bounce period
τb , and the pulse average becomes the bounce average. In
contrast, chirps in the nonadiabatic regime are associatedwith the propagation of convective plume-like structuresin phase space density as in Fig. 3 (a). These structures advance radially faster than they rotate, taking the form ofa convective instability [4]. In this case, a bounce time τbcannot be easily defined, and this is the reason why we
9 In Eq. (9), we write γ /ν instead of γ /ω in order to be able to
L 0
L0make a clearer distinction between strong ( 1) and weak ( 1) drive.
Apart from this, we usually work with γL /ω0 .
10 As indicated in Eq. (10), the variations of a particle’s canonicaltoroidal momentum Pζ and energy E = K + QeΦ are closely related
in the system studied here. We will use this relation later (cf., Eq. (26)).

1403087-10

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

have chosen the amplitude pulsation time τpulse to measure(non)adiabaticity in Eq. (10).
Nonlinear frequency chirping always begins nonadiabatically, with the propagation of convective wave fronts
in whose wake dynamics reminiscent of convective interchange and wave breaking can be observed in plots ofδ f (Pζ , ϑ), with corresponding chirping patterns appearing
in the spectrograms. In Figs. 2 and 3, this turbulent belt can
be clearly seen in the form of complicated patterns aroundthe seed resonance. The region around the seed resonancetends to remain nonadiabatic, but solitary vortices with aneﬀectively adiabatic core can be emitted from its boundary. Thus, the same system may simultaneously contain
adiabatic and nonadiabatic portions. If the evolution of the
wave field is dominated by one of the adiabatic vorticalstructures in phase space, Eq. (10) can be written as112πδω̇τb
= τ2b δν̇ =<12
τchirpωb(single smooth
adiabatic chirp); (11)where τ−1chirp ≡ δν(τb ) ≈ τb δν̇ is the frequency shift that thewave experiences during one bounce time.
In summary, we will refer to a subdomain of phasespace as being “adiabatic” when the dynamics in that region exhibit a distinct separation of time scales, where aparticle’s mean nonlinear displacement δPζ during onebbounce cycle τb is smaller than the maximal displacementmax|δPζ |b during that period. The physical relevance ofsuch a time scale separation is that it facilitates the formation of robust vortex structures in phase space, that canbe viewed as a generalization of Bernstein-Greene-Kruskal(BGK) modes formed by electron Langmuir waves [23].
We will show that the large vortex that is visible in the upper part of Fig. 3 satisfies the adiabaticity condition only
marginally (τb  τchirp /2), but it nevertheless resembles a
BGK mode since its interior consists of concentric nested
layers as illustrated in Fig. 1 (c), with the remarkable ability to withstand the strong and incoherent beats of the field.

2.6

Discussion of the adiabatic limit

The scale separation in the adiabatic limit is useful for
analytical treatments and has allowed Berk, Breizman etal. (BB) to develop a model that yields a quantitative prediction of the chirping rate if the values of γL and γd are
known [10, 24]. With some approximations, their integral
relation between δω and ωb reduces to the analytical formδωBB = ωBB
b
with

2γd
t
3

1/2

ωBB
b =

≈ 0.44 × γL (γd t)1/2 ,16
γL ≈ 0.54 × γL .
3π2Note that since A ∝ ω2b , Eq. (12b) predicts the field amplitude associated with a chirping BGK mode. Comparison
Dewar≈ 0.18, so the field amwith Eq. (6) shows that ωBBb /ωbplitude of the BGK wave is predicted to be 30 times smaller
than the saturation amplitude of the undamped wave driven
by gradients that give the same value of γL . Conversely,
this means that a BGK wave is predicted to require a
1/0.18 = 5.6 times larger value of γL in order to yield thesame bounce frequency (and amplitude). We suspect that
this is connected with the assumption of marginal stability
Dewar. Consequently, we expectthat underlies ωBBb but not ωbto find reasonable agreement with ωb in simulations whereDewar∼ O(0.1 . . . 0.2).
(γL − γd )/γL is comparable to ωBBb /ωbIn our simulations, this marginality parameter will have thevalues 0.06 and 0.47 and we will find bounce periods thatlie slightly above and below the prediction τBBb ≈ 0.24 ms.
The dotted green curve in Fig. 2 (b) indicates the theo√
retical prediction δνBB ∼ 24 kHz × t[ms] from Eq. (12a)for our parameters. One can see that the curve lies in theright “ball park”, especially for the upward chirp12 . Similar agreement in the overall extent of the chirps was seen
in many other studies, even when the parameters clearly
stretch the range of validity of the theoretical model, as inour marginally unstable but strongly driven example.
This remarkable robustness of the model may be attributed to the fact that Eq. (12) is the result of balancing
the rate at which the fluctuating field receives energy fromthe destabilizing gradients (γL ) and loses energy via thead hoc damping γd . The theory simply postulates that the
BGK wave is somehow able feed on the gradients at thesame rate as the field loses energy, implicitly relying on theexistence of a nonadiabatic layer as illustrated in Fig. 1 (c),without modeling the detailed processes inside that layer.
Therefore, the overall agreement means only that both theory and simulation satisfy the energy balance constraint. It
does not necessarily mean that the theoretical picture on
which the theory is built captures the processes that actually occur in the simulation (or in a real plasma)13 .
A generalized model that accounts for nonlinear distortions and shrinking of the field’s eﬀective potential wellhas been derived by Breizman [25] and further extended to
expanding potential wells by Hezaveh et al. [16]. A key
point that is captured in these and related studies of longrange chirping [26–29], and which is closely related to the
topic of the present paper, is that the generalized BGK-likemodes have an active boundary layer that can dynamically
expand or shrink as the field amplitude grows or decays14 .

(12a)
(12b)

11 The quantity α = δω̇/ω2 has been used as an expansion parameter inbtheoretical treatments. Here, we included a factor 2π since we feel that
this quantifies the upper limit of adiabaticity in a more precise and intuitive way: 2πα = τb /τchirp < 1 instead of α  1. An indication for thiscan be seen in Fig. 13 of Ref. [22] and we will oﬀer further justificationin Sec. 5.4 when inspecting the motion of tracer particles in a vortex.

12 One reason for the smaller extent of the downward chirp is that it
propagates away from the peak of the field mode in our setup.
13 As an analogy, consider the phenomenon of so-called magnetic reconnection, which can be realized by diﬀerent physical mechanisms. All
mechanisms break the frozen-in-flux condition of ideal MHD, but the energy conversion rates and side eﬀects (such as wave excitation) diﬀer.
14 Our interpretation of the results presented in Ref. [25] is that theleading-order eﬀect is the variation of the potential well’s depth. The distortions of its shape seem to be corrections of higher order, comparable tovarious other eﬀects (e.g., geometric nonuniformity) that were ignored.

1403087-11

Plasma and Fusion Research: Regular Articles

2.7

Volume 16, 1403087 (2021)

Some open questions: Formation of solitary vortices and the role of beating

Having reviewed the concepts that will be needed forthe discussion of our results, we proceed now with an outline of some open questions that we wish to tackle.
One question that, to our knowledge, has not been satisfactorily answered yet is how solitary vortices are formed
and emitted from the turbulent belt around the seed resonance. We believe that some of the results reported here
can contribute to a better understanding of the nonadiabatic
onset of chirping and how it can lead to the detachment ofsolitary vortices. We will present evidence indicating thatthese processes are facilitated by pulsations of the field amplitude A(t) like those in Fig. 6 (b).
Such pulsations and accompanying phase jumps havebeen reported in a previous analysis of moderate burstsof chirping shear Alfvén waves seen in JT-60U tokamak
experiments [12]. The patterns were seen both in the experimental data and in self-consistent hybrid simulations
with realistic EP sources and collisions, and they were explained (loosely speaking) as the beating between multiple
chirps that occur simultaneously at diﬀerent frequencies.
Some examples were shown in Fig. 13 of Ref. [12].
Indeed, the chain of events sketched in Fig. 5 anddiscussed in Sec. 2.3 implies that one chirp rarely comesalone. As is illustrated schematically in the bottom part of
Fig. 5, we may distinguish two ways for beating to occur:
5a. interference of wave pairs ω± produced by frequency
splitting (described in Sec. 2.3);5b. interference between diﬀerent generations of chirps
ω(1) , ω(2) , . . . (like those seen in Fig. 2).
The physical picture we will use may be described as follows. The convective interchange of EP Vlasov fluid produces multiple more or less coherent structures in the perturbed EP phase space density δ f that are more or lessstratified radially along the vertical (Pζ ) axis in Fig. 3. As
indicated by the horizontal arrows in Fig. 3, these structures propagate at diﬀerent phase velocities around theplasma torus (horizontal axis and out of plane), depending on the local value of Pζ . These EP density waves
drive the same field mode at diﬀerent frequencies simultaneously. The field mode responds by beating as shownin Fig. 7, as if it consisted of multiple field waves, eachlocked to a diﬀerent pump wave. Indeed, the total signals(t) = A(t) sin(−ω0 t−φ(t)) in Eq. (1) represents an arbitrary
number of harmonic waves, whose superposition causes
the combined amplitude A(t) to pulsate and the combinedphase φ(t) to jump by ±π between pulses as in Fig. 7. In
Appendix A, we show that this physical picture is mathematically consistent with the field equations we solve.
When we represent the time-dependence of the combined signal in the form s(t) = A(t) sin(−ω0 t −φ(t)), we demand that its harmonic components must all have the same
spatial structure. One may argue that representing EPdriven shear Alfvén waves in this form is an oversimpli-

Fig. 7 Example of global beating in a short time window aroundt ≈ 2 s of Fig. 6 (b). Panel (a) shows the raw signal s(t)(gray) and the envelope A(t) (red). The equidistant vertical yellow stripes help to discern changes in the signal’s
phase. Panels (b) and (c) show, respectively, the global
mode structure of the ideal MHD displacement in code
units at times of maximal constructive and destructive interference. While global beating is inevitable in the semiperturbative model used in the ORBIT code with which
this signal was computed, it is remarkable that such behavior was also observed in self-consistent hybrid simulations (cf. Fig. 12 of Ref. [12]).

fication, because a self-consistent (nonperturbative) fieldmode may simply decompose into several spatially separated nonlinear modelets, whose amplitudes and phases
can evolve more or less independently, each driven by its
own resonant phase space structure. This is certainly a possible scenario, but it is not a necessity, as one can infer froma remarkable observation made in the above-mentioned
self-consistent hybrid simulation of a JT-60U plasma [12]:Fig. 12 of that paper shows that the energetic particle mode
(EPM) [30] with toroidal mode number n = 1 that was responsible for the chirps in that system sometimes vanishes
completely between beats and reappears with perfectly opposite phase. In principle, the field in that hybrid simulation could have decomposed into multiple independent
EPMs, but this did not occur, at least not at a readily visible level (i.e., not at leading order). The same behavior isreproduced here in idealized form in Fig. 7.
We interpret the possibility of global beating as a consequence of a certain degree of rigidity that is inherent
to MHD modes with long wavelengths and which is further enhanced by large magnetic drifts. For instance, in
the above-mentioned JT-60U scenario, EP orbits spanned
20. . . 30% of the minor radius (cf. Fig. 3 of Ref. [31]). Thisrigidity preserves the coherence of the wave field across asignificant portion of the plasma radius and one might say

1403087-12

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

that this is the reason why we call it a mode. All chirps
studied in the present paper are caused by radially propagating phase space structures that stay within the range ofan unperturbed drift orbit, which justifies a posteriori theassumption of a rigid mode whose evolution is captured by
the time-dependence of two scalars, A(t) and φ(t).
One physical consequence that is highly relevant forthe present study is that the resonant phase space density
waves are mutually coupled and interfere with each othervia the jointly driven MHD field mode. Such an interaction
has been explicitly ignored in the model of chirping BGK
modes formulated by Berk et al. (Eq. (23) of Ref. [24]).
We believe that this is a point where interesting unexploredphysics may still be lurking and we are aware of one moregroup of researchers who are tackling this topic by considering the eﬀect of beating on chirps driven by a pairof electron resonances [16]. It is diﬃcult to foresee all
the implication of this complex feedback; apart from theperhaps obvious fact that all particles are temporarily detrapped at the times of destructive interference, and that
the EP density waves find themselves out-of-phase when
the field wave regrows after performing a phase jump as inFig. 7. Thus, instead of independent BGK waves, we havea complex multi-body system consisting of several more orless strongly coupled components.
Such multi-component interactions can become moreregular in cases where chirping is suppressed, e.g., by processes that scatter the resonant particles. In such cases, thespectrum of a resonantly driven mode can still broaden, but
it does so in a well-structured manner by splitting into a sequence of discrete lines with equal frequency spacing Δν
as reported by Fasoli et al. [32]. According to the theory
presented by Fasoli et al. [33] based on Berk et al. [10],such phenomena can be explained by a period doubling bifurcation of a near-threshold kinetic instability. The modulation of the field amplitude (beating) through the interference between the waves associated with such a bunch ofspectral lines is an integral part of the process, but in thatcollision-dominated (non-chirping) regime the phase space
dynamics as well as the beats of the field become regular,
settling into a cyclic pattern. The feature that makes our
chirping scenarios (with weak or no collisions) more complicated is that irregular convective motion can persist, especially near the seed resonance. Fascinatingly, however,coherent structures in the form of long-lived phase spacevortices do form in this turbulent environment. One purpose of our study is to throw light on the role that beating
may play in the formation of such robust vortices.
One may argue that beating has a zeroth-order eﬀect
only on the motion of individual particles, while its eﬀect
on the coherent phase space structures is of higher order,
because the interaction is nonresonant15 . However, chirping itself is often a higher-order eﬀect, especially in the
15 The beats may resonate with other regions in phase space, which isan interesting story on its own, and not pursued here.

limit of slow (adiabatic) chirps. Therefore, we think that it
is worthwhile to consider the eﬀects of beating on the samefooting. In this paper, we take a close look at how beating aﬀects the phenomenon of nonlinear frequency chirping and the underlying motion of Vlasov fluid in EP phasespace, which we represent here using discrete particles.

3. Model and Methods
3.1 Global beating as a rationale for a semiperturbative model
MHD-PIC hybrid simulations often compare wellwith experiments, at least for EP-driven Alfvén modes
with long wavelengths (n < 10) [12, 34–36]. However,such simulations require high-performance supercomputers, primarily because of the expensive field solver thathas to evolve the electromagnetic fields, plasma densityand pressure on a dense three-dimensional (3D) mesh consisting of tens of millions of grid points. The globality ofthe beating observed in Ref. [12] oﬀers a justification for asemi-perturbative approach, where the spatial structure of
the Alfvénic field modes is prescribed, reducing the MHD
solver’s large number of degrees of freedom to the evolution of merely two scalars; namely, the amplitude A(t)and phase φ(t) appearing in Eq. (1). Numerical studies ofnonlinear frequency chirping in realistic tokamak geometry can then be performed on a laptop PC.
In the presence of multiple modes k = 1, . . . Nmode ,with toroidal mode numbers nk and initial seed frequencies
ω0k , the field model in Eq. (1) can be written asξ̃(ψP , ϑ, ζ, t) =

N
mode
l=1

Ak (t)

ξ̂k,m (ψP )eiΘk,m (t) , (13a)m

Θk,m (t) = nk ζ − mϑ − ω0k t − φk (t);(13b)

where ξ̃ is the complex-valued ideal MHD displacement
vector and ξ̂k,m (ψP ) is the radial profile of each poloidal
Fourier harmonic m, here written as a function of normalized poloidal flux 0 ≤ ψP ≡ ΨP /ΨP,edge ≤ 1 that increases
monotonically from the center (0) to the edge (1) of theplasma. The imaginary unit is denoted by i, and Θk,m is thecomplex phase at a given point in time and space, here expressed in toroidal coordinates (ψP , ϑ, ζ) with poloidal andtoroidal angles ϑ and ζ. Our convention is that the physicalelectric and magnetic field perturbations E and δB are thereal (cosine) component of the complex signal, so that thepotentials Φ and δ A, and, thus, the ideal MHD displacement ξ = {ξ̃}, are the imaginary (sine) components.
Reduced field models based on Eq. (13) have been
implemented in codes such as ORBIT [17, 18, 37], HAGIS[38, 39] and MEGA [40, 41], which have then been used to
study nonlinear EP-Alfvén wave interactions and chirpingin realistic geometries. The reduced number of degrees
of freedom in the field dynamics enhances physical transparency, and the low computational cost combined with aset of free parameters (Table 1), permits detailed physics

1403087-13

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

studies through extensive parameter scans that are not feasible with self-consistent simulations. Meanwhile, the fullgeometric complexity of the EP orbits is retained, allowing
us to explore practically relevant conditions with respect tospatio-temporal scales and geometric nonuniformity.
While global beating gave us an excuse to ignore perturbations of the mode structure in Eq. (13), the implicit
assumption made in the field Eq. (2) of the existence ofa global “mode” that can oscillate well in a wide rangeof frequencies is a rather strong and disputable simplification of the Alfvénic plasma response. The dynamics inthe reduced system (2) do show some preference for resonant frequencies at the radius of the mode’s peak: We
observe that an oﬀ-peak seed resonance tends to be produce stronger and more rapid chirps in the direction of themode’s peak (hence the asymmetry in Fig. 2 (b)). However,
we are missing entirely the eﬀect of continuous spectra ofshear Alfvén waves, and their band gaps that may containdiscrete eigenmodes. The radial structure ωA (ψP ) of these
spectra is known to aﬀect chirping by providing a preferredtrajectory for the chirps (e.g., see Refs. [42–44]) and bycontributing to the field damping [45]. For the time being,
we choose to ignore these eﬀects without further justification other than saying that this simplification appears to be
acceptable for our purposes. Other factors that can play arole but are ignored here include fluid nonlinearities, collisions and sources of EPs. Due to these limitations, thescope of our study is constrained to the qualitative featuresof chirping in the semi-perturbative limit, with no intention to make quantitative predictions for experiments. Further considerations related to nonperturbative eﬀects canbe found in Appendix D.1.

3.2

Equations and parameters

In the present work, we use the ORBIT code [17] withthe reduced δ f field-particle interaction model derived by
Pinches et al. [38], extended by Chen et al. [18], and reformulated by White et al. [13,19]. First simulations of chirping using ORBIT were reported in [13,37]. We use the same
MHD equilibrium as in those recent papers, which resembles a plasma of the conceptual reactor FIRE with reduced
field strength B0 = 0.49 T at the magnetic axis located at
R0 = 2.15 m. The profiles characterizing the magnetic
geometry are shown in Fig. 8. The following Eqs. (14)(22) are normalized by the deuteron cyclotron frequency
Ωc0 = QeB0 /M (inverse time unit) and the magnetic flux
density B0 (field strength unit) at the magnetic axis.
The time-independent axisymmetric magnetic field Bof the torus and time-dependent perturbations that resemble ideal incompressible electromagnetic flute modes are
represented in Boozer coordinates as
B = g∇ζ + I∇ϑ + BΨP ∇ΨP (equilibrium),(14a)δB = ∇ × αB

(perturbation),(14b)(nq − m) Φ̂k,m
α̂k,m =
(gq + I) ω0k(ideal MHD),

Fig. 8 Equilibrium parameters and profiles for the magnetic
field in Eq. (14a). This is the same case as in Ref. [13].

Ψ=(m − nq) ξ̂k,m(mg + nI) q(radialdisplacement),(14c)

where α = δ A · B/B2 = δA /B represents the parallel component of the perturbed vector potential, Φ is the electrostatic potential, q(ψP ) = dΨ/dΨP is the field line helicity(tokamak safety factor), and Ψ the toroidal flux16 . The coordinate surfaces ψP = const. and ϑ = const. of the equilibrium field B can be seen in Figs. 7 (b,c) (orange lines).
The fluctuating potentials α and Φ are decomposed like ξ
in Eq. (13a). The radial profile of the contravariant radialΨ= ξ̂k,m · ∇Ψ/B0 of the ideal MHD displacecomponent ξ̂k,mment vector defined in Eq. (13) is given as input (specified
below).
The stepping equations for the amplitude and phaseappearing in the Fourier representation of Eq. (13) are17
Np −2A0 dAk=
wk, j {S k,m, j } − γd Ak ,dt
ω0k Dk j=1 m

(15a)

Np 
−2A0 dφk=
wk, j {S k,m, j },dt
ω0k Ak Dk j=1 m(15b)with the on-axis Alfvén velocity A0 , particle index j, and
S k,m, j = ρ B2 α̂k,m (ψP ) − Φ̂k,m (ψP ) eiΘk,m , (16)
j

2
Ψ
2
Dk = 4π(17)dψP [ξ̂k,m (ψP )] ,m

where ρ =  /B = (gζ̇ + I ϑ̇)/B2 with the parallel guiding
center (GC) velocity  = gc · B/B. The GC phase space
is sampled by j = 1 . . . Np particles, whose distributionG(ψp , ϑ, ζ, K, μ, t) satisfies dG/dt = 0 and whose weights
16 Strictly speaking, our Ψ and Ψ

P are the fluxes divided by 2π.
17 Our Eqs. (15)-(17) correspond to those in Ref. [13], except that here

the lengths are not normalized by R0 , so the velocity A0 (normalized byΩc0 ) appears here instead of the frequency νA0 = A0 /R0 . Moreover, we
let Dk = 1 since only a single Fourier harmonic m/n = 6/5 is used.

1403087-14

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Table 1 Model parameters used in our ORBIT simulations. Cases (A) and (B) will be examined in detail. The undamped cases (A0) and(B0) are used for some comparisons in Sec. 4, and Appendix B contains further details. Case (C) is used in Appendix D.4 to
examine the eﬀect of changing the mode peak location relative to the resonance, which can aﬀect the eﬀective seed frequency aswell as the direction and rate of chirping. The values of the gradients GF00 ∂Pζ F0 (dominant) and GF00 ∂K F0 (negligible) are those of
the input variables df0dpz and df0de in ORBIT when A0 = falf = 0.07 and (for historical reasons) F0 /G0 = 1 + df0dpz. The
values of γL = γ + γd were measured in the simulations during the phase of exponential growth using Eq. (5).

Case
(A)
(B)
(A0)
(B0)
(C)

Harmonic
ν0 [kHz] m/n
100
6/5
100
6/5
100
6/5
100
6/5
100
6/5

Mode profile
r0 /a rw /a
0.65 0.15
0.65 0.15
0.65 0.15
0.65 0.15
0.7
0.15

Drive
G0
F0 ∂Pζ F 0

2.801642.801640.46227
1.79305
5.53325

G0
F0 ∂K F 0−0.084898−0.084898−0.014008−0.054335−0.16767w j = δ f j /G are evolved as
F0d ln F0
1 dF0
dw
=
=−
,
−w
dt
G dt
G0
dt

(18)

beginning with δ f (t = 0) = 0. Here, δ f represents theperturbation around the initial reference distribution F0 of
physical guiding centers, and G0 = G(t = 0) is the initial
distribution of simulation particles (= phase space markers) that will be specified below. The particle trajectories
are determined by solving equations for ρ (t), ΨP (t), ϑ(t)and ζ(t), whose explicit form can be found in Ref. [17] orSec. 3.9.2 of Ref. [19]. These equations of motion can becompactly summarized in Hamiltonian form,ϑ̇ = ∂Pϑ H,
ζ̇ = ∂Pζ H,

P̃˙ ϑ = −∂ϑ H,
P̃˙ = −∂ H,ζ

(19)

H = ρ2 B2 /2 + μB + Φ,

(20)

ζ

with the Hamiltonian

and the canonical angular momentaP̃ζ = (ρ + α)g − ΨP ,P̃ϑ = (ρ + α)I + Ψ.

(21)

For marker particle loading and phase space diagnostics,
we use the canonical toroidal angular momentum of theunperturbed system,Pζ = gρ − ΨP = P̃ζ − αg.

(22)Equations (19) and (20) for the particle motion can be
regarded as exact within the realm of drift-kinetic theory18 .
The field-particle coupling is overestimated by the neglectof the gyroaveraging eﬀect, but otherwise the assumptions
18 The polarization drift and the ponderomotive force [46] are ignored.

All other aspects of GC motion are retained to the extent permitted by a
Hamiltonian formulation in Boozer coordinates, which requires neglecting the field component BΨP (denoted by “δ” in Eq. (2.26) of Ref. [19])that arises from the coordinates’ nonorthogonality. The model is appropriate for fast ions interacting with shear Alfvén modes.

γL /ω0⇒ 0.080⇒ 0.075⇒ 0.005⇒ 0.035⇒ 0.081

Damping
γd /ω0
0.075
0.040
0
0
0.075

Remarks
Marginally unstable
Strongly unstable(A) without damping(B) without damping(A) with peak shifted outward

made in the derivation hold in the parameter regime considered here. Meanwhile, Eq. (15) for the field should be
regarded as a toy model because the derivation involves theneglect of terms containing products of more than one time
derivative, such as ∂2t Ak , ∂2t φk , ∂t Ak ∂t φk . These terms are
not strictly negligible by mere ordering arguments in a system exhibiting rapid amplitude pulsations and phase jumps
that are typical for cases with nonlinear chirping (Fig. 7).
In addition, the factorization (13) is based on a neglectof the background plasma response (continuous spectra,etc.) as discussed in Sec. 3.1 above. Here, we accept thesecaveats and simply postulate that the field mode evolves in
the manner prescribed by Eq. (15).
The model parameters are listed in Table 1. The fluctuating field in our simulations is initiated with a seed frequency of ν0 = ω0 /(2π) = 100 kHz and contains only asingle toroidal harmonic n = 5 and a single poloidal harmonic m = 6, so we omit the mode label k. Its radial profile
has a Gaussian shape,(23)ξ̂Ψ (r) = exp (r − r0 )2 /rw2 ,as shown in Fig. 9 (a). Contour plots of the mode structure in the toroidally rotating poloidal plane at ζ  = ζ −ω0n t = 2πl with integer l are shown in Fig. 7 (b) in cylinder coordinates (R, Z, ζ  ) and Fig. 9 (b) in polar coordinates(ψP , ϑ, ζ  ). Due to the large magnetic drifts performedby the circulating energetic deuterons in our simulations,which have kinetic energies K ≥ K0 = 85 keV and magnetic moment μB0 = 2 keV, this mode has a rich harmonic
content from the point of view of the particles [31], so thateﬃcient transit resonances [47]ω = nωζ − pωpol ,(24)with various values of p are possible. As one can inferfrom the kinetic Poincaré plot in Fig. 9, our seed resonancehas p = 4 poloidal elliptic points. In the (P̂ζ , ζ) plane (notshown), there are n = 5 elliptic points, so  = 1. The
quantities ωpol = 2π/τpol and ωζ = Δζ(τpol )/τpol = 2π/τζ

1403087-15

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 9 Structure of the field mode and seed resonance. The panels on the left show (a) the radial profile ξ̂Ψ (ψP ) and (b) gray-scale contoursof the mode structure in the (ψP , ϑ)-plane in the toroidally rotating frame of reference where ζ  = ζ − ωn0 t = 2πl with integer l. The
dots in (c) represent a uniformly randomized distribution of 100 particles loaded along the line E  = E− ωn0 Pζ = const. (≈ 75.6 keV)in the range 0.1 ≤ ψP ≤ 0.5 on the outer midplane ϑ = 0 (horizontal dashed line in (b)). Kinetic Poincaré plots of the particle
trajectories in the perturbed system with fixed amplitude A0 = 10−3 and phase φ = 0 are shown in (b) and (d), with colors
representing energy E = K + QeΦ. The seed resonance for this mode with ω0 = 2π × 100 kHz and m/n = 6/5 is located atP̂ζ,res = 0.719. It has p = 4 elliptic points in (P̂ζ , ϑ) and n = 5 elliptic point in (P̂ζ , ζ) (not shown). In the simulations, particles
are loaded uniformly randomized in ψ2P , whereas here they are uniformly randomized in ψP and followed for 500 transit times τζ0 .

Table 2 Default parameters for particle loading and pushing.

Num. of marker particles
Loading range @ ϑ = 0
Kinetic energy @ ψP,max
Magnetic moment
Loading duration
Pushing time stepNp[ψP,min , ψP,max ]K0 [keV]μB0 [keV]
T load /τζ0Δtstep /τζ0106[0.1, 0.5]
85
2
10
1
1
400 , 1000

in Eq. (24) are the angular frequencies for an unperturbed
transit in the poloidal and toroidal directions, respectively.
The toroidal transit period on axis,
q0 g0 + I0
2π
≈ 4.7 μs,τζ0 ≡ √
q0 B0
2K0 /M 

(25)≈R0 =2.15 m

will sometimes be used here as a unit of time: tˆ ≡ t/τζ0 .

3.3

Initialization and stepping

Table 2 contains the default parameters used for loading and advancing the simulation particles, which is doneusing a standard 4th-order Runge-Kutta scheme. With time
steps of size Δtstep /τζ0 = 1/400 (damped cases) or 1/1000(undamped cases), we obtain an accuracy that is suﬃcientfor our purposes (see Appendix D.6).
The default number of particles in our simulations
is Np = 106 . An approximate quiet start is obtainedby loading only a portion of the simulation particles at
each time step Δtstep during an interval of 10 transit timesτζ0 (see Appendix D.5). The particles are injected at
the outer midplane at ϑ = 0 (horizontal dashed line in
Fig. 9 (b)), where we distribute them along the toroidal
angle and normalized poloidal flux as ζ j = 2πx j andwith uniform randomψP j = ψP,min + (ψP,max − ψP,min )x1/2
jnumbers 0 ≤ x j ≤ 1. With uniform particle weighting(F0 /G0 = const.), this choice yields an approximately uniform gradient ∂Pζ F0 ≈ −∂ΨP F0 = const.19
The magnetic moment is fixed (μB0 = 2 keV) and the
kinetic energies are chosen to satisfy [48]E  = E − ω0 Pζ /n = const.,(26)

as shown in Fig. 9 (c). The simulations reported here wereobtained with E(t = 0) = K0 = 85 keV at (ψP , ϑ) =(ψP,max , 0), so that E0 ≈ 75.6 keV. During the course of thesimulations, deviations around this value will be on the order of 3%, so we can assume E  ≈ E0 = const. Eﬀectively,
we are simulating a thin slice of the phase space structures,
which have a cylindrical shape in the (P̂ζ , E  , ϑ)-space.
Figure 10 shows a 2-D projection in the (P̂ζ , E  )-plane.
The black box indicates a slice of width ΔE0 = 7.5 keV
as used in a simulation performed for validation purposes(see Appendix D.2). All simulations discussed in the mainpart of this paper were performed with ΔE0 = 0.
19 In the radial direction, the initial marker distribution G , which mea0sures the number of simulation particles dN per cell dψP is chosen tosatisfy G0 ∝ dN/dψP ∝ ψP , which is realized by loading particles with
dN/dx = const. and a uniformly randomized variable x ∝ ψ2p . For uniform particle weighting F0 /G0 = const., this yields ∂ψP F0 = ∂2ψP N ∝γL = const. With ψP = (r/a)2 this corresponds to an EP density profiledN
dN
of the form nEP (r) = rdr∝ dψ∝ 1 − ψP = 1 − (r/a)2 .
P

1403087-16

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 10 Phase space structures (e.g., holes and clumps) have theform of bent cylinders in (P̂ζ , E  , ϑ)-space. This plotshows a projection of the marker distribution into the(P̂ζ , E  )-plane, with colors representing δ f . The dynamics are primarily in the vertical direction (E  = const.),so we simulate only a thin slice of width ΔE0 as indicated by the black rectangle. In most cases, only thecenter line was simulated; i.e., ΔE0 ≈ 0. The only exception is case (A’) in Figs. D1 and D2 of Appendix D.2,
where ΔE0 = 7.5 keV. In that region, the phase spacestructures are approximately perpendicular to the axis of
the modified radial coordinate P̂ζ = P̂ζ − (E  − E0 ) × C
with C ≈ n/ω0 , which is indicated by the long diagonalarrow.

as large as γd /ω0 = 4%. Both cases exhibit strong chirping and an overview of the observations made in diﬀerent
stages of these simulations is given in Figs. 11 and 12.
Panels (a) and (c) show the time trace of the amplitude A(t) obtained by solving Eq. (15a). The kinetic drive
plotted in panel (b) was measured by computing the expoΔA(t)+0.99×γ(t−Δt)nential running average γ(t) = 0.01× A(t)Δtof the increment ΔA on the right-hand side of Eq. (15a) andadding the damping rate γd . The above values of the linear drive γL correspond to γk = γ + γd during the earlyphase of (nearly) exponential growth. Panel (d) shows
the spectrogram obtained by Fourier analyzing the signals(t) = A(t) sin(−ωt − φ(t)) with a sliding time window ofsize Δtwin = 100 τζ0 ≈ 0.47 ms.
The lower parts of Figs. 11 and 12 show five snapshotsof the perturbed EP phase space density as δ f -weighted kinetic Poincaré plots. Our Poincaré section is the torodiallyrotating plane ζ  = ζ− ωn0 t = 2πl with integer l, where particle weights w = δ f /G are accumulated in a histogram-likefashion on a mesh consisting of 200 cells in the poloidal
angle −π ≤ ϑ ≤ π, and 200 . . . 1000 cells in the radial
range 0.5 ≤ P̂ζ ≤ 0.9. The accumulation was performed
for a short time interval of one toroidal transit, τζ0 ≈ 4.7 μs,
during which each particle crosses the ζ  = 2πl plane approximately p = 4 times. The resulting quantityt+τζ0 /2dt δ f (t |ζ  = 2πl),δ f transit (P̂ζ , ϑ, t) =(27)

t−τζ0 /2

After loading the markers during the interval −10 ≤tˆ ≤ −Δtˆ, with tˆ ≡ t/τζ0 , the kinetic instability is kickstarted at t = 0 with a small perturbation A0 . We use
A0 = 10−6 in the marginally unstable case (A). In thestrongly unstable case (B), we sometimes use a smaller
A0 = 10−8 with the purpose of extending the duration of
exponential growth enough to capture the onset of chirping in the Fourier spectrogram.
A comprehensive characterization of our simulation
model, including convergence and sensitivity tests, is given
in Appendix D. Note that if one starts the instability simulation with a relaxed non-axisymmetric EP distributionthat has adapted to the initial perturbation of the field, as
in Fig. 4, the system immediately enters the stage of nonlinear beating, eﬀectively skipping the exponential growth
and the first beat. This is demonstrated in Apendix D.3.

4. Structure Formation and Saturation
As a preparation for studying the eﬀects of beating, we
examine in this section the processes of structure formation, nonlinear saturation and the onset of beating in cases(A) and (B), whose parameters are given in Table 1. Case(A) is marginally unstable, with a linear drive γL /ω0 ≈ 8%
only slightly larger than the damping γd /ω0 = 7.5%. Case(B) is strongly unstable, with γL /ω0 ≈ 7.5% nearly twice

will be referred to simply as “δ f ”.
In the following subsections, we examine the dynamics occurring during the early part of the simulations, until
the instabilities saturate; namely, the first three snapshotsin Figs. 11 and 12. This includes the birth of hole-clump
pairs and the onset of beating.

4.1 Linear response, phase mixing and halos
In the linearized system, the magnitudes of field perturbations, density perturbations and particle displacements bear no meaning. There exists only a trend that isexpressed by the complex frequency ω̃ = ω + iγ whose
value is directly determined by the gradients of the chosen reference state. Strictly speaking, resonant instabilities(also known as “nonlinear Landau damping”) are neverlinear because they always involve finite displacements.
However, nearly exponential growth does occur initiallybecause the imposed finite-amplitude perturbation spontaneously induces flows with a phase relation such that∂t ln A ≈ ∂t ln δ f (independent of amplitude and position).
This phase relation changes when nonlinearities begin todominate over the linear response. This transition can beclearly seen at the beginning of our simulations in the formof a marked change in the structure of δ f . In this subsection, we will take a closer look at this process, becausesome of the insights gained will be useful later.

1403087-17

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 11 Overview of the dynamic stages in the marginally unstable case (A) (same as in Fig. 2). Panels (a) and (b) show time traces of the
field amplitude A(t) and the instantaneous kinetic drive (growth rate γ(t) plus damping γd ) during the first 210 transits. Panel (c)shows A(t) for a larger portion of the simulation (1800 τζ0 ≈ 8.5 ms) and panel (d) shows the corresponding Fourier spectrogram
obtained with a 0.47 ms time window. The vertical dotted lines in (a)-(d) indicate the times at which the five snapshots of theδ f -weighted kinetic Poincaré plots were taken that are shown in the lower part of the figure. The last two snapshots are the sameas those shown in Fig. 3, except that the contrast at small amplitudes is enhanced here by using a nonlinear color scale (squareroot). As in all contour plots of δ f that appear in this work, we have inverted the direction of the P̂ζ axis for easier comparison
with the spectrograms, so that structures with higher frequencies (smaller P̂ζ , larger radii) appear at the top (cf. Fig. D14).

Fig. 12 Overview of the dynamic stages in the strongly unstable case (B) with reduced initial perturbation amplitude, A0 = 10−8 . Results
are arranged as in Fig. 11. The first two snapshots of Fig. 11, where the resonance gradually emerges, are very similar in thepresent case, so they are not shown again. Instead, we show additional snapshots of δ f during the first two beats (100  tˆ  140).

The structure of δ f in the first snapshot (tˆ = 0.5) inFig. 11 shows the linear short-time response of the EP distribution caused by the E × B drift associated with theabruptly applied oscillating field perturbation. The attribute “linear” refers to the fact that the structure of δ freflects only the instantaneous trend of the motion of EP

Vlasov fluid. No appreciable particle displacement has
occurred yet, since only a very short time (few μs) has
passed and the field amplitude is still very small (≈ 10−6 ).
A zoom-up is shown in Fig. 13 (a), where black contour
lines indicate schematically the structure of the resonanceand the bold black arrows indicate the trend of fluid mo-

1403087-18

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

tion. Along the red stripes (δ f > 0), there is a net outward
motion of denser fluid into regions previously occupied bylighter fluid, causing a local increase of the EP phase spacedensity. In the blue stripes (δ f < 0), the density drops due
to a net inward flow of lighter fluid.
Returning to the full-size snapshot in Fig. 11 (left):
The periodicity of δ f (ϑ) indicates that a p = 4 resonancedominates in the region P̂ζ  0.62 and a p = 5 resonancedominates in the region P̂ζ  0.62. The p = 4 seed resonance will appear at P̂ζ,res = 0.718 as indicated by a horizontal dashed-dotted line. However, in the first snapshot attˆ = 0.5 in Fig. 11 this resonance is not visible in the structure of δ f yet. This may be interpreted as a consequenceof the fact that the short time interval that has passed corresponds to a large uncertainty Δν ∼ 1/τζ0 ≈ 200 kHz infrequency. Another (closely related) interpretation is thatphase mixing (discussed below) has not had enough time
to become eﬀective in the frequency band we are looking
at; it requires dozens of transits to become noticeable.
Of course, the most eﬃcient coherent accumulation ofphase space density perturbations does occur near the resonance, where the particles move more or less in phase withthe field mode. For instance, after 10 transit times, the resonance width becomes approximately Δν ∼ 1/(10τζ0 ) ≈20 kHz, which corresponds to 0.7  P̂ζ  0.74 (cf.
Eq. (D6)) and is precisely the region where δ f is large inthe second snapshot at tˆ = 10.5 in Fig. 11.
Outside that (near-)resonant layer, the magnitude ofδ f is smaller by an order of magnitude and represents thenonresonant ideal MHD displacement in Eq. (D4). The
corrugated structure of δ f in the oﬀ-resonant domain as
seen in the second snapshot at tˆ = 10.5 in Fig. 11 may beinterpreted as a superposition of patterns produced by two
eﬀects that reflect the emergence of nonlinearities:• One eﬀect is phase mixing of the remnants from theinitial perturbation in δ f that was produced by theabruptly imposed field at t = 0. Phase mixing isrealized by the radially sheared poloidal propagation(horizontal arrows in Fig. 11) of the phase space structures due to the convective nonlinearity ϑ̇∂ϑ δ f .
• The second eﬀect is that δ f in our Poincaré plots increasingly shows the cumulative response; namely,
the time integral of the weight Eq. (18).
Phase mixing, also known as “Landau damping”, causes
a progressive tilting and thinning of the oﬀ-resonant red
and blue stripes in our snapshots of δ f . However, this isvisible only for a short time because the second eﬀect —the cumulative δ f response — eventually dominates.
Consequently, the meaning of the patterns we seein δ f -weighted Poincaré plots changes when the system
leaves the linear and enters the nonlinear regime. Let
us elucidate this change of meaning using Fig. 13, whichshows zoom-ups of (a) the first snapshot (tˆ = 0.5) and (b)the third snapshot (tˆ = 150.5) of Fig. 11, highlighting thediﬀerences between the linear and nonlinear δ f structures.

Fig. 13 Illustration of the meaning of patterns in kinetic
Poincaré contour plots of the EP phase space density
perturbation δ f . The data are taken from (a) the firstand (b) the third snapshot of Fig. 11 for case (A). Here,
only three of the resonance’s four periods are visible.
Black lines indicate schematically the eﬀective positions of phase space islands and neighboring flow contours that one would obtain after waiting infinitely longwith the field amplitude A and phase φ fixed. Bold arrows indicate the flow direction. The island width issort of realistic in (b) but strongly exaggerated in (a). In
fact, the island does not really exist in this fully nonadiabatic stage of the simulation. Here and in the following, we sketch “eﬀective islands” only to indicate
the phase of the field and the instantaneous directionof the associated flow lines. The phases in (a) and (b)have been aligned by using the modified poloidal angleϑ = ϑ − ϑ0 + λϑ Δν0 t, where λϑ is the wavelength of thecentral island (2-headed arrow), Δν0 is the prompt frequency shift (cf. Appendix D.4), and ϑ0 = 0.105π aligns
the O-point with ϑ = 0. Halos in (b) are indicated byalternating “+” and “–” signs.

When compared to the linear short-time response in
Fig. 13 (a), one can see that the vertical stripes outside theresonant layer in Fig. 13 (b) are shifted to the left by −λϑ /4above the resonance (P̂ζ < 0.718) and to the right by +λϑ /4
below the resonance (P̂ζ > 0.718), so that the sign of δ finside those stripes flips across the resonance. Here, λϑ
is the poloidal wavelength of the p = 4 resonance, whose
value is about λϑ ≈ 0.365π in the center of Fig. 13 (b) as
indicated by the two-headed arrow. (For our GC drift orbits, the value of λϑ in Boozer coordinates increases fromthe outer midplane ϑ = 0 towards ϑ = ±π.)
While the structure of the δ f pattern has changed in
the way we have just described, the structure and phase of
the resonance, which is schematically illustrated by black
contour lines in Fig. 13, remains the same in panels (a) and(b) in the frame where the (promptly shifted) resonanceis stationary. Thus, the δ f -weighted Poincaré plots in the
nonlinear regime (b) do not represent the instantaneous
trend as in (a), but the accumulated modulations in the EP

1403087-19

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

phase space density relative to the background. In otherwords, the nonlinear δ f patterns in (b) tell us where dense
fluid (red) and light fluid (blue) has been accumulated untilthat time; i.e., the history of transport. The red arrows in
Fig. 13 (b) indicate where dense fluid has recently moved
radially outward (resulting in a local density increase), andthe blue arrows indicate where light fluid has moved radially inward (resulting in a local density reduction).
Another consequence of this cumulative nature of thenonlinear δ f patterns is that they are fairly robust. A particularly impressive example of this robustness is the fact thatthe oﬀ-resonant stripes indicated by alternating “+” and“–” signs in the third snapshot (tˆ = 150.5) in Fig. 11 andin Fig. 13 (b) are almost perfectly vertical. These robust
nonlinear density modulations, which seem to emanate radially from the nearest hole and clump wave fronts, will in
the following be referred to as halos.
The values of δ f inside the halos is small, so we use anonlinear color scale (square root) when we want to visualize them clearly. Their amplitude tends to decrease withincreasing distance from the resonant δ f structures, which
can be attributed to phase mixing. However, the shearing
eﬀect of phase mixing is almost invisible. This implies thatthe density modulations that constitute these halos are decoupled from the particle motion along ϑ. Thus, the halos
are of an entirely collective nature.
The halos exhibit a large degree of phase-locking withrespect to the field oscillations that, in turn, are controlled
by the dominant resonant δ f structures. This implies thatthe halos can tell us the eﬀective locations of O- and Xpoints as in Fig. 13 (b), so they are a welcome substitutefor kinetic Poincaré plots like that in Fig. 4 (a), which are
not applicable when the field varies rapidly in time. For instance, the halos will help us to analyze the eﬀect of phasejumps associated with beating (Sec. 5.1).
The processes that we have described in this subsection are universal for resonant field-particle interactions(regardless of growth) and, therefore, are qualitatively thesame in all cases studied here. Due to this reason, we haveomitted the phase mixing stage in Fig. 12 and, instead, included additional snapshots of the saturation stage that will
be discussed next. The subsequent dynamics depend on theparameters, so we choose to proceed with a comparison between the saturation processes with and without damping
in the marginally unstable and strongly unstable case.

4.2

First saturation with and without damping

Beating and nonadiabatic chirping begin during thefirst peak of a resonant instability, so it is instructive to examine the saturation process in some detail. We will do sohere by comparing damped cases (A) and (B) that do chirpwith undamped cases (A0) and (B0) that do not. The latterare obtained by setting γd = 0 and reducing the gradientF0 ∝ γL such that the growth rate γ = γL − γd remains thesame.
For the marginally unstable cases (A) and (A0), withγ/ω0 ≈ 0.5%, Fig. 14 (a) shows the evolution of the field
amplitude A(t) during the first 2 ms. After growing at thesame exponential rate, the instability in case (A) saturates
at a low amplitude of A ≈ 6 × 10−6 at time tˆ = 154.5(0.73 ms). Panels (b) and (c) show the structure of thephase space density perturbation δ f in a portion of the(P̂ζ , ϑ) plane when A ≈ 6 × 10−6 in both cases (A) and(A0). The corresponding perturbations of the radial profileof the EP density distribution, F(P̂ζ ) = F0 + δ f , measured
across the O-point (red) and across the X-point (blue) are
shown in panels (f) and (g). A ϑ-averaged density profile(green) is also plotted.
One can see that the O-point profiles in Fig. 14 (f) for
the damped case (A) and in panel (g) for the undamped
case (A0) have both significantly reduced gradients at theseed resonance, with the profile of (A) appearing somewhat flatter than in (A0). Note that the vertical axes ofpanels (f) and (g) have diﬀerent scales, which diﬀer by afactor F0 (A)/F0 (A0) ≈ 16. Flattening a steeper gradientrequires a larger perturbation δ f , which is indeed an order of magnitude larger in case (A) than in case (A0) on
the left-hand side of Fig. 14, where the field amplitudes areequal. Recalling that the energy released during the interchange of denser and ligher EP Vlasov fluid that is represented by δ f has all been transferred to the field wave, thefact that the values of δ f in cases (A) and (A0) diﬀer by anorder of magnitude for equal field amplitudes, reflects thefact that the field in the damped case (A) has lost nearly all
the energy that it had previously received from the initial
EP density gradient.
In contrast to the undamped case, the field in thedamped case does not possess enough energy to reverse
the radial displacement of the particles. Equivalently, one
may say that only a small amount of trapping occurs before the field’s energy is exhausted and its amplitude drops
to a small value. This explains the diﬀerences seen in the
structure of δ f in panels (d) and (e) on the right-hand sideof Fig. 14: While the EP phase space fluid revolves aroundthe resonance in the undamped case (A0) in panel (e), thefirst-generation hole and clump waves in the damped case(A) in panel (d) has departed radially from the resonance,
closely followed by a second-generation hole-clump pairthat has formed during the second pulse of the field.
This multi-layered structure of the near-resonant δ fpatterns in panel (d) makes it diﬃcult to determine the effective locations of O- and X-points, and this is where thehalos introduced in Sec. 4.1 above come in handy. The
halo-based O- and X-points are indicated in panel (d) bydashed and dotted lines, respectively. If one takes into
account that the eﬀective resonance has shifted by Δϑ =−λϑ Δν0 Δt = −0.365π × 0.4 kHz × 0.43 ms ≈ −0.06π (cf.
Fig. 13) during the time between snapshots (b) and (d), it
becomes clear that the O- and X-point locations are eﬀectively reversed. This is the consequence of a phase jump

1403087-20

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 14 EP phase space density perturbation during the saturation of (A) damped and (A0) undamped instabilities in the marginallyunstable cases with γ/ω0 = 0.5%. Panel (a) in the center shows the time traces of the field amplitude A(t). The situation at low
amplitudes A ≈ 6×10−6 , where case (A) peaks for the first time, is shown on the left-hand side. The situation at higher amplitudes
around 3 × 10−5 , where case (A) peaks for the second and case (A0) for the first time, is shown on the right-hand side. Panels(b)-(e) show contour plots of δ f (P̂ζ , ϑ), and panels (f)-(i) show the corresponding perturbed EP density profiles F(P̂ζ ) = F0 + δ f :one measured at the poloidal angle of the O-point, one at the X-point, and one poloidal average. The eﬀective locations of O- andX-points were inferred from the halos (cf. Sec. 4.1) and are indicated, respectively, by dashed and dotted vertical lines in panels(b)-(e).

Fig. 15 EP phase space density perturbation during the saturation of (B) damped and (B0) undamped instabilities in the strongly unstablecases with γ/ω0 = 3.5%. The results are arranged as in Fig. 14. The initial perturbation amplitude in these simulations was
A0 = 10−6 , so the first saturation in case (B) occurs earlier here than in Figs. 12, 20, 24 and 27, where we used A0 = 10−8 .

Δφ = π that was performed by the field during the amplitude minimum between snapshots (b) and (d), aroundt ≈ 0.9 ms, which will be examined in detail in Sec. 5.1.
The dynamics in the strongly unstable cases (B) and(B0) in Fig. 15, with γ/ω0 ≈ 3.5%, seem to be qualitatively
similar to the marginally unstable cases (A) and (A0), buteverything is faster and larger in both magnitude and spatial extent. Here, the saturation amplitudes diﬀer only by
about a factor of 2, and the magnitudes of δ f diﬀer by afactor of about 2. . . 3, consistently with the ratio of the initial gradients, F0 (B)/F0 (B0) ≈ 2.
After the undamped instabilities in cases (A0) and(B0) saturate, their amplitudes remain near the level ofthe first peak. At the instant of saturation, Figs. 14 (i) and
15 (i), show an O-point profile that is flat around the resonance. In fact, the profile near the O-point can be expected
to flatten for the first time after half a bounce cycle. Indeed,
the first flattening of the O-point profile near the resonancecan be seen long before the amplitude A(t) reaches its first
peak. In case (A0), with the resolution of our diagnostics,
this occurs at about 30% of the peak amplitude (shortly after the profile snapshot in Fig. 14 (g)), and in case (B0) thisoccurs at about 18% (before the snapshot in Fig. 15 (g)).
At first glance, this might look like an overshoot, but thatis not necessarily the case because, in order for the net outward flow of EPs (= net energy transfer to the mode) to1403087-21

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

cease, it does not suﬃce to flatten the gradient only alonga thin line across the O-point. A suﬃcient amount of profile flattening has to occur in a suﬃciently wide portion of
the eﬀective island’s poloidal length λϑ , which is the distance between X-points along the resonance line as indicated in Fig. 13 (b). The structure of δ f inside the phasespace islands in Figs. 14 (e) and 15 (e) implies that A(t)stops growing only after the deeply trapped particles havecompleted about 1 to 11/2 revolutions20 around the nearby
O-point (see also Fig. B1). In other words, the steady-statelevel around which the field amplitude A(t) subsequently
oscillates is reached when a majority of the particles havemore or less returned to or passed their initial positions.
This observation is corroborated by the time traces of individual particles as shown in Figs. B2 (a,b) of Appendix B
and is consistent with Dewar’s theoretical prediction [7].
Another noteworthy observation is that the saturationin the damped cases (A) and (B) occurs only after the density profile across the O-point has been flattened substantially. This is especially remarkable in the strongly damped
marginally unstable case (A) in Fig. 14 (f): since that casewould be linearly stable if the initial density gradient had
been just a little less steep, one might naïvely think thatthe instability would saturate as soon as the gradient at theresonance has been reduced below that threshold, whichshould require only a small perturbation in δ f . This, however, does not happen and it reflects the fact that, evenduring the phase of approximately exponential growth, theinstantaneous stability of the perturbed system is a more
complex problem than the linear stability of the initial reference state, whose gradient is uniform along the resonance. In the perturbed state, the gradient ∂F/∂Pζ withrespect to the action variable P̂ζ varies along the angle variable ϑ, which is precisely the feature we discussed in theprevious paragraph, where we were saying that saturationin the undamped cases required the gradient to flatten in acertain range of ϑ around the O-point. In the damped cases,that range of ϑ is shorter, requiring only a small amount ofrevolution (trapping) around the resonance.

4.3

Potentially misleading density profiles
and the cause of pulsations and phasejumps

As mentioned in Sec. 2.2, we emphasize again thatdensity gradients ∂Pζ F can have a stabilizing or destabilizing eﬀect only when there are flow lines in phase space thatare directed uphill or downhill, respectively. Those lines offlow are not always easy to determine when the field’s amplitude and phase evolve rapidly, and without that information, it can be diﬃcult to interpret the role of the gradients
correctly. For instance, the gradients seen in a radial crosssection at an O-point are relatively ineﬀective for resonant
20 The first revolution is part of the resonant trapping process, but its
duration is longer than the bounce period τb because the first revolution
includes the formation process of the field’s eﬀective potential well.

drive or damping, because they are oriented transversely tothe lines of flow at that location (even after phase jumps by±π). Instead, it is primarily the degree of profile flattening in the region between the O- and X-points that determines at what stage a resonant instability ceases to grow(saturates), as mentioned in the previous section. A simple
direct link between destabilizing gradients and the growthof the field perturbation exists only in the linear regime,
where the gradient is still uniform and the relative phases
are such that ∂t ln A ≈ ∂t ln δ f regardless of magnitude andposition. The relative phases between the field and the density waves in δ f at diﬀerent radii begin to change well before the instabilities saturate, and the evolution of the fieldis then determined by the field’s spatial correlation with
δ f ; namely, the integral in Eq. (2), where gradients do notappear explicitly. The gradients enter the field equations
only in a complicated implicit way through the time integration of the equations of motion — or, equivalently, theVlasov Eq. (4) — that yield δ f .
These considerations show that radial profiles likethose shown in the lower parts of Figs. 14 and 15 should
be interpreted with great caution. For instance, it might
be tempting to hypothesize that the steep gradient at the
X-point in Figs. 14 (f) and 15 (f) is the force that causes
the ±π jump in the field’s phase φ(t) during the subsequent
amplitude minimum and allows its amplitude A(t) to growagain to a large value during the next pulse, where the Oand X-point positions are switched [13]. However, uponcloser inspection this intuitively appealing hypothesis turns
out to be self-contradictory: As illustrated in Fig. 16 (a),the form of δ f around the X-point is such that the gradients would actually be stabilizing when the X-point would
be replaced by an O-point. This is because the counterclockwise flows would have to move dense fluid (red) uphill and light fluid (blue) downhill. Even if one considers
a gradual phase shift, a contradiction arises: Energy may
be extracted from the X-point gradient if the flow contours
would shift to the right (ϑ̇ > 0) as illustrated in Fig. 16 (b),which corresponds to a negative phase shift in the field(φ̇ < 0). However, as we will show later, the first phasejump in our simulations is towards the left (ϑ̇ < 0, φ̇ > 0;
cf. Fig. 17 (b)), which can be attributed to an asymmetryin our simulation setup that strengthens the leftward moving pump waves; namely, the clumps traveling downhill
towards the mode’s peak.
In other words, from the energetic point of view, one
may say that phase jumps and subsequent pulses occur not
because but in spite of the steep X-point gradient. This
statement holds not only at the time of the amplitude peakshown in Fig. 16, where we illustrated the situation for ahypothetical phase shift. The X-point gradient continues
to play a stabilizing role all the way to the amplitude minimum, where the actual phase jump occurs (cf. Fig. 19).
This indicates that something else must be driving theamplitude pulsations and phase jumps, and our suspects
are the counter-propagating EP fluid packets of increased

1403087-22

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)
Table 3 Estimates based on Eq. (28) of the frequency shift andbeat pulse length corresponding the maximal radial separation ΔP̂ζ of the primordial hole-clump pair during
the first pulse in cases (A) and (B).

Case
(A)
(B)
Hole-clump separation ΔP̂ζ
 0.005
 0.07Frequency split Δν (upper bound)  2.5 kHz  35 kHz
Pulse length τpulse (lower bound)  0.4 ms  0.03 ms
Fig. 16 Schematic illustration of how the perturbed EP distribution δ f near the eﬀective X-point in Fig. 15 (b) would be
aﬀected by the hypothetical phase space flows (arrows)that one would obtain (a) by exchanging O- and Xpoints, or (b) by slightly shifting the flow contours rightward (ϑ̇ > 0, φ̇ < 0). The situation in (a) corresponds toa phase jump by ±π and the situation in (b) corresponds
to a negative phase shift by about Δφ ≈ −π/3.

(δ f > 0) and decreased (δ f < 0) density that can be seenin panel (b) of Figs. 14 and 15. In this work, we choose
to interpret the rapid field pulsations and phase jumps as aconsequence of the beats that result from the superpositionof multiple density waves in δ f that simultaneously drive
the field wave. This is more easily said than proven ingeneral, but at least in the limit of two-wave beating, one
can demonstrate that our ansatz is consistent with Eq. (15),as shown in Appendix A.

4.4

Primordial hole-clump pairs, frequency
splitting and the onset of beating

The packets of dense and light fluid in Figs. 14 (b) and(c) may be regarded as a primordial hole-clump pair. The
hole wave propagates faster than the seed wave, and theclump wave lags behind, which means that there are now
two pump waves in EP phase space density δ f with elevated and reduced frequencies ω+ and ω− as was illustrated schematically in Fig. 5. The similarity between panels (b) and (c) in Fig. 14 implies that such a pair of pumpwaves exists transiently also in the undamped case during
the first cycle of hole-clump revolution around the O-point.
Their eﬀect is analyzed in Appendix B, where we present
evidence for transient frequency splitting. In cases withstrong damping, the frequency splitting is sustained andthe superposition of the two pump waves produces a beating pump wave signal, to which the electromagnetic field
wave responds likewise.
In the limit of symmetric beating between a pair ofpump waves with equal amplitudes but opposite signs (cf.
Eq. (A5)), their combined signal is sin(ω+ t) − sin(ω− t) =2 cos(ωt) sin( 21 Δωt) with base frequency ω = 12 (ω+ +ω− ) ≈ω0 and pulse frequency Δω = ω+ − ω− . Using Eq. (D6)in Appendix D.7, we can estimate the magnitude of frequency splitting and the initial pulse length of the beat,τpulse = 1/Δν = 2π/Δω ≈ 2 × 10−3 ms/ΔP̂ζ ,(28)from the hole-clump separation in phase space, ΔP̂ζ . The
results are summarized in Table 3.
In the marginal case (A), we have ΔP̂ζ  0.005 at
the time of amplitude peaking in Fig. 14 (b), which corresponds to a frequency split of Δν  2.5 kHz and a pulse
length of τpulse  0.4 ms. The interval between the firstpeak and the first amplitude minimum in Fig. 14 (a) is
about (0.91 − 0.73) ms = 0.18 ms, which corresponds toa beat pulse length of 0.36 ms close to that estimated fromthe hole-clump separation. This observation supports ourinterpretation of the first amplitude minimum as being a
consequence of interference between the two pump waves
constituted by the primordial hole-clump pair.
In the strongly unstable case (B), we have ΔP̂ζ  0.07
in Fig. 15 (b), which corresponds to an upper bound ofΔν  35 kHz for the frequency splitting, and a lower
bound of τpulse  0.03 ms on the pulse length. A glanceat the spectrogram in Fig. 12 (d) shows that the actual frequency splitting at the time where the instability in case(B) saturates is significantly smaller than the upper limitof 35 kHz. A detailed measurement shows that Δν at thattime is only about 10 kHz (cf. Fig. 20 (f)), which corresponds to a pulse length of τpulse ≈ 0.1 ms. This value
is consistent with the interval between the first peak andthe first amplitude minimum in Fig. 15 (a), which is about(0.58 − 0.54) ms = 0.04 ms and corresponds to a beat pulse
length of ≈ 0.08 ms. As in case (A), this supports ourinterpretation of the first amplitude minimum as being a
consequence of pump wave beating. The fact that the actual values of Δν and τpulse in case (B) diﬀer significantly
from their respective upper and lower bounds in Table 3
also makes sense. The upper and lower bounds should
yield accurate estimates only in the limit of marginal instability, which was indeed true for case (A) discussed in
the previous paragraph. With increasing departure frommarginality due to reduced damping as in case (B), the initial frequency splitting should become smaller as particle
trapping becomes stronger and inhibits the phase slippageof the hole and clump waves with respect to the seed resonance. Eventually, in the absence of damping, sustained
frequency splitting does not occur (Δν → 0), whereas ΔP̂ζ(= island width) can be large.
The small (retarded) frequency splitting at early times

1403087-23

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

means that weakly damped scenarios like our case (B)should undergo a relatively abrupt large frequency splitting
at later times, when the large primordial hole-clump pairescapes the grasp of the seed resonance. In contrast, theonset of chirping should be a more gradual process in themarginally unstable case (A) with strong damping, where
nearly no trapping occurs and phases slip when the holeclump pair is still small. These expectations are confirmed
numerically and supported by theory in Sec. 5.2 below.

(i). The ensuing wave front propagation and nonadiabatic
onset of chirping (ii) is analyzed in Sec. 5.2, followed by ananalysis of the formation and detachment of solitary vortices (iii) in Sec. 5.3. While the first three subsections deal
with the collective structures in EP phase space density,
the last Sec. 5.4 is dedicated to an analysis of the motionof individual simulation particles, which yields information about global transport and the internal dynamics ofthe large vortical structures that manage to survive in thepresence of a beating field.

5. Analysis of the Role of Beating
In the above overview Figs. 11 and 12, the snapshotsof the EP phase space density perturbation δ f taken after
the saturation of the instabilities show the following processes:(i) the formation of hole-clump wave pairs (until thethird snapshot),(ii) the radial propagation of hole and clump wave fronts,
in whose wake a turbulent belt forms (fourth snapshot), and(iii) solitary hole and clump vortices (fifth snapshot).
Here, the word “front” will be used to refer to the steep δ fgradients at the inner and outer edges of the turbulent belt
around the seed resonance. These fronts are indicated byhorizontal dashed lines in the fourth snapshot in Figs. 11
and 12. The distinction between interchanging hole &clump waves on the one hand and “turbulence” on theother is somewhat arbitrary, based mainly on appearancein scale and complexity. In both cases, there is a convective interchange of light and dense packets of EP Vlasov
fluid, with plumes of dense fluid advancing outward andvice versa. These phase space structures are nonadiabaticand may be regarded as a semi-perturbative manifestationof what Wang Tao et al. [44] call “convective branch”21 .
All these structures propagate poloidally at characteristic frequencies that depend on their instantaneous radial
location measured by P̂ζ . Their (increasingly uncoordinated) pumping action drives the oscillations of the field.
As soon as two or more pump waves are present in the EP
phase space density, their interference can cause the field
to beat. Among other interesting features, we will show
shortly that, at least during the first dozen or so beats, thewave fronts advance radially in pulses that correlate wellwith the beating of the field wave.
In Sec. 5.1, we characterize the beats and their feedback on the phase space dynamics at the beginning of thesimulations, where the first hole-clump wave pair forms
21 In Ref. [44], the “convective branch” was only responsible for relatively minor transient chirps, while the dominant signal was associatedwith an MHD eigenmode that was referred to as “relaxation branch”. In
our semi-perturbative model, there are no discrete eigenmodes, but theseed wave may be regarded as the counterpart of the “relaxation branch”,
and the early nonadiabatic chirps may be regarded as manifestations of
the “convective branch”. The (semi-)adiabatic phase space structures thatare responsible for long-range chirping may then be referred to as “vortex
branch”. Here, all three branches have the same field mode structure.

5.1 Phase jumps and halo dynamics during
symmetric and asymmetric beats
Amplitude pulsations and phase jumps are ubiquitousin experimental measurements of chirping modes and insimulated chirps. The first explicit experimental reportsthat we are aware of came from JET [49], where these dynamics were interpreted as a signature of chaos. Here,
we do not go to such lengths and interpret the pulsations
and phase jumps simply as a manifestation of beating (asin Ref. [12] and Appendix A). In this section, we perform a detailed analysis of the phase jumps in our simulations, with focus on the associated phase space dynamics and feedback loops. Our results are summarized
in Fig. 17, which shows the dynamics of the field (left)
and EP phase space density perturbations (right) during
the first few beats. The beats consist of amplitude pulsations as shown in panel (a) and phase jumps between thepulses as shown in panel (b). The situation is representative for all our simulations with strong damping, both near
marginal stability and strongly unstable. Here, we discuss
the marginally unstable case (A).
Before we examine the results, we should point out
that, in most simulations, we observe a prompt frequency
shift during the first few transits. This phenomenon is attributed to a self-optimization process and is discussed inAppendix D.4. The oscillation frequency in case (A) becomes ν0 + Δν0 ≈ (100 + 0.4) kHz and the resonance location P̂ζ,res shifts from 0.719 to ≈ 0.718. Similarly to theadjustment made to the poloidal angle ϑ in Fig. 13 above,
the phase φ(t) in Fig. 17 (b) has also been adjusted by acorresponding oﬀset φoﬀs = 2πΔν0 t. The location of thepromptly shifted eﬀective resonance is indicated by dashdotted lines in panels (c)-(i) of Fig. 17.
During pulse #1, the field amplitude in Fig. 17 (a) saturates around tˆ = 150 (0.7 ms), and then it rapidly drops,
here by about a factor of 3. Around the instant of the first
amplitude minimum, at about tˆ = 194.5 (0.9 ms), we observe in panel (b) a arctan-like phase jump Δφ = +π thattakes about 20 transits (0.1 ms) to complete. Around thesame time, the normalized Fourier spectrum in panel (c)shows transient frequency splitting into ν+ ≈ 102.4 kHz
and ν− ≈ 98.1 kHz, which corresponds to a radial separation of |ΔP̂ζ |  5 × 10−3 around the resonance and iscomparable to the extent of the phase space density pertur-

1403087-24

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 17 Onset of beating due to the interference of hole & clump (H&C) pump waves in the marginally unstable case (A). Panels (a) and(b) show, respectively, the time trace of the field amplitude A(t) (log scale) and phase φ(t) − φoﬀs (t), where the oﬀset φoﬀs = 2πΔν0 t
cancels the phase drift associated with the small prompt frequency shift Δν0 ≈ 0.4 kHz (cf. Appendix D.4). Panel (c) shows theFourier spectrogram whose amplitude has been normalized at each time step. Panels (d)-(f) and (g)-(i) show snapshots of thephase space density perturbation δ f around the first two minima of the amplitude in (a), where phase jumps by +π and −π canbe observed in (b). The dash-dotted horizontal lines in (c)-(i) indicate the eﬀective resonance location (ν0 + Δν0 = 100.4 kHz,P̂ζ,res = 0.718). A nonlinear color scale (square root) is used to make the halos visible in (d)-(i). It is important to note that thearrows in (d)-(i) indicate halo motion, which are collective structures, not tied to particle motion. Arrows are drawn with diﬀerent
thickness for easier identification. A detailed view of the third amplitude minimum at tˆ = 312 (1.467 ms) between pulses #3 and#4 can be gained from Fig. 18. The phase jump in the raw signal between pulses #10 and #11 (out of range) was shown in Fig. 7.

bations seen in the contour plots of δ f in panels (d)-(f). A
qualitatively similar behavior can be seen between pulses#2 and #3, during the second amplitude minimum aroundtˆ = 270.5 (1.3 ms), where the spectral peaks are located at
ν+ ≈ 103.0 and ν− ≈ 98.2.
The upper branch ν+ is pumped by a clump wave inthe region P̂ζ < P̂ζ,res (upper part of panels (d)-(i)) and thelower branch ν− is pumped by a hole wave in the regionP̂ζ > P̂ζ,res (lower part of panels (d)-(i)). According to
the theory in Appendix A, the interference between thesetwo pump waves causes beating in the kinetic drive on theright-hand side of the field Eq. (15), to which the field responds with identical beats. Thus, the amplitude minima
highlighted in Fig. 17 correspond to instances of destructive interference that occurs when pump waves with opposite signs of δ f (namely, a hole-clump pair) are aligned.
Consistently with Eq. (13b), which yields the conditionpϑ + φ = const. in the Poincaré section nζ − ωt = 2πl withinteger l, we observe positive phase jumps Δφ = +π when
the clump wave dominates (ϑ̇ < 0), whereas Δφ = −πwhen the hole wave dominates (ϑ̇ > 0).
There are quantitative diﬀerences between the firstand second phase jumps in Fig. 17. During the second
event, panel (a) shows a deeper and sharper amplitude
minimum, and panel (b) shows a phase jump that is more
abrupt and in the opposite direction (Δφ = −π). The spectrogram in panel (c) shows us the reason for these diﬀer-

ences between the first and second amplitude minima: thefrequency splitting is asymmetric during the first event,with a dominant upper branch ν+ , whereas relatively symmetric splitting is seen during the second event, with an
only slightly (2%) stronger lower branch ν− . This meansthat the interfering pump waves have diﬀerent strengths
during the (asymmetric) first beat, while they are comparable during the (symmetric) second beat.
The beating dynamics of the field wave become increasingly complicated as time advances. For instance,in the presence of more than two pump waves — pairsof which may dominate and control the field oscillations
at diﬀerent times — pulses may overlap and slow phasejumps may be cut short (|Δφ| < π) by faster ones.
So far, we have discussed the beating phenomenon
only from the point of view of the field behaving as adriven oscillator. However, in our system there is a feedback loop, by which the pulsations and phase jumps of the
field wave acts back on the motion of individual particles aswell as on their collective density waves. One manifestation of this feedback are the oﬀ-resonant halos introducedin Sec. 4.1, which appear in the form of vertical stripes of
red and blue shade in the snapshots of δ f on the right-handside of Fig. 17. We emphasize again that these halos are
nonresonant collective modulations of EP density, which
are not tied to the motion of individual particles. Our observations indicate that the halo patterns are largely phase-

1403087-25

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

locked to the field wave (that, in turn, is controlled by theresonant phase space structures), which makes these halos
useful for the present analysis of phase jumps.
We distinguish “near halos” (meaning: near-resonant)and “far halos” (meaning: far oﬀ-resonant). Diﬀerent halo
dynamics are observed during phase jumps associated withsymmetric and asymmetric beats:(i) Symmetric beat: Positive (red) and negative (blue)
halos abruptly switch places globally during instantsof perfect or nearly perfect destructive interference oftwo pump waves. This is realized in Fig. 17 (g)-(i).
(ii) Asymmetric beat: Far halos remain locked to thedominant pump wave, while the subdominant pumpwave controls only the nearest halo. Phase slippage
is realized here by halos being torn and reconnectedwith their neighbors as seen in Fig. 17 (d)-(f).
Scenario (i) is simpler, so we start with Figs. 17 (g)-(i).
While the hole and clump waves counter-propagate alongϑ at a more or less steady pace (with only minor jolting),
the far halos tend to stay locked in place during the preceding pulse #2, so that they stay behind the resonant structures that they are attached to. This results in some bending of the near halos. Towards the end of pulse #2, thefar halos begin to accelerate as indicated by the arrows inpanel (g): The halos above the resonance are attached tothe clump wave and tend to shift to the left as indicated by
the pair of red arrows. Similarly, the halos below the resonance are attached to the hole wave and tend to shift tothe right as indicated by the pair of blue arrows. Here, anarrow’s thickness is used merely to identify it in diﬀerent
snapshots. When the field disappears at the time of snapshot (h), the halos jump ahead of the respective clump orhole wave. Then they decelerate and more or less freeze
at the new positions seen in snapshot (f). The motion of
the halos corresponds closely to the time trace of the field
phase φ(t) in panel (b), whose sign is determined by that ofthe dominant phase space structure. Here, the hole wave
is only about 2% stronger, so the pump waves are almost
equally strong, causing a very abrupt and radially global
phase jump of the halos, which remain attached to the nearest hole or clump wave front.
More or less symmetric beats do occur fairly often
in our simulations, but more frequently we observe asymmetric beats, where one of the pump waves dominates, sothat destructive interference is not complete, and amplitude minima are smoother and shallower. This situation,
namely scenario (ii), is realized during the first amplitude
minimum in all of our chirping simulations. Panels (d)-(f)of Fig. 17 show the corresponding phase space dynamics.
In the present case, the clump wave dominates and we observe that not only the halos above the resonance follow
the clump wave, but even the far halos below the resonance do so, as indicated by the red arrows in panel (d).
The weaker hole wave controls only the near halo as indicated by the blue arrow in panel (d). During pulse #1,

Fig. 18 Evolution of δ f during the third amplitude minimum
at tˆ = 312 (1.467 ms) in Fig. 17. This event is similar
to that in Fig. 17 (e), but here the signal-to-noise ratio
is larger. As in Fig. 17, the dotted rectangles indicate
the region where halo tearing and reconnection occurs.
The arrows indicate halo motion (decoupled from particle motion).

all structures are more or less locked in phase, so the arrows in panel (d) indicate only the trend of their relative
motion. When the field is weakened due to partial destructive interference at the time of snapshot (e), these relative
drifts in ϑ become significant, and a phase jump of magnitude π occurs, causing the red and blue halos in panel (f)to exchange places compared to panel (d). However, thistime, all far halos have moved to the left with the clump,as one can infer from the positions of individual arrows in
the snapshots. This means that around snapshot (e) the halos below the resonance have been torn and reconnected in
the region that is roughly indicated by the dotted rectangle.
Unfortunately, the signal-to-noise ratio is somewhat low in
Fig. 17 (e), so we show another example of an asymmetricbeat in Fig. 18. Here one can see snapshots taken just before and just after the instant of halo reconnection during
the third amplitude minimum in Fig. 17 (a).
In summary, we have shown in this subsection howphase jumps occur during amplitude minima. See also themovie in Fig. E3 (b). Hole and clump waves and nearby
halos propagate poloidally at a relatively uniform pace that
depends on their distance from the resonance. In contrast,
the motion of the far halos is nonuniform, with a arctanlike time trace that resembles the evolution of the phaseφ(t) in Fig. 17 (b). In the case of symmetric beating, withtwo pump waves of similar intensity, the phase jump is
global along P̂ζ . In the case of asymmetric beating, in
our example dominated by the clump wave above the resonance, the phase jump is realized by tearing and reconnection of near and far halos on the side of the subdominant
pump wave, in our example on the hole side22 . The phasejumps have important consequences that we will examine
next.

22 Halo reconnection always occurs somewhere, but during symmetricbeats it occurs far away from the resonance, outside our field of vision.

1403087-26

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 19 Illustration of the eﬀect of a phase jump, using the snapshots of δ f in Fig. 17 (g)-(i) around the second amplitude minimum in case (A). The solid lines in (a) and (c)indicate schematically the positions of eﬀective phasespace islands (similarly to Fig. 13). The island widthis somewhat exaggerated. Islands are absent during theamplitude minimum in (b), so the dashed lines indicate
their past and future forms. Bold arrows indicate theflow direction. Dotted lines highlight the inclinationof the halos in the vicinity of the resonance (far halos
above and below the visible P̂ζ range are almost vertically straight).

5.2

Feedback loop for pulsed front propagation and nonadiabatic onset of chirping

In Sec. 4.1, we have asserted that the halos that appear in our δ f -weighed kinetic Poincaré plots can be used
to identify the eﬀective instantaneous location of O- andX-points. Here, this procedure is applied to the three snapshots from Figs. 17 (g)-(i) that we analyzed in the previous
subsection. Figure 19 shows schematically the eﬀective
phase space islands (black sinusoidal curves) plotted on
top of the phase space density perturbation δ f at the timesof those snapshots: (a) at the end of pulse #2, (b) at the
instant of the phase jump by Δφ = −π, and (c) at the beginning of pulse #3.
The arrows in Fig. 19 (a) show how, at the end of apulse in the field amplitude A(t), the fronts of hole & clump(H&C) waves are approaching the eﬀective X-point, whichmeans that they are pulled “inward” (towards the resonance line). However, this process is terminated by the disappearance of the field during the amplitude minimum at
the time of snapshot (b), before any significant revolution
around the eﬀective O-point can occur. Around that time,the H&C waves are free to drift poloidally at their characteristic speeds, as indicated by the arrows in Fig. 19 (b).
When the field is revived in the next pulse, it does so with

opposite phase. This means that the new O-points appear
at (or near) the locations of the former X-points and viceversa, yielding a situation as illustrated schematically in
Fig. 19 (c). The H&C wave fronts that were previously
drawn towards the resonance have now begun to travel“outward” (away from the resonance).
This constitutes an important feedback mechanism.
The new location of the eﬀective island contours (representing flow lines) with respect to the location of the H&C
pair at the time of the phase jump facilitates a further departure of the H&C structures from the seed resonance.
The further outward motion of the clump and inward motion of the hole transfers energy to the field. The growing
field amplitude makes the eﬀective island larger and theflow contours steeper, allowing the H&C structures to depart even further. This positive feedback continues until
the poloidal drift of the H&C takes them across the eﬀective O-point, reversing the trend: The flow contours at thehead of the H&C structures point towards the resonance,
making clumps move inward and holes outward, so that themode transfers energy to the particles and is, therefore, resonantly damped. As mentioned previously in Sec. 4.2, theadditional damping γd makes the process irreversible since
the mode’s remaining energy is already exhausted after apartial reversal of radial H&C displacement.
This process is repeated during each beat cycle. The
resulting evolution of the system during the first dozen or
so beats of the marginally unstable case (A) is shown in
the left half of Fig. 20. Panel (a) shows the evolution of
the field amplitude A(t), with the first 9 pulses labeled (thesmall peak of pulse #1 is hardly visible). Panel (b) showsthat the fronts of the hole and clump waves advance radially in a pulsed manner that is perfectly correlated with thepulsations of the field seen in panel (a). One can see thatthe radial expulsion of the wave fronts is partially reversed
at the end of most pulses, but the net displacement clearly
increases with time, which is consistent with the discussion in the previous paragraph. Note that there is often asignificant diﬀerence in the magnitude of the displacementof the hole and clump wave fronts during a given pulse, especially when the amplitude does not drop very far due toasymmetric beats or due to an overlap between two pulses.
Qualitatively similar behavior is seen in the strongly
unstable case (B) shown in the right half of Fig. 20. The
main diﬀerences are the shorter time scale (factor 2. . . 3)
and the much larger amplitude of the first pulse (factor
100). In the marginal case (A), the magnitude of the displacement of the wave fronts during pulse #1 is below
the accuracy of our diagnostics, whereas pulse #1 in thestrongly unstable case (B) does more displacement workthan all the subsequent pulses that are visible in Fig. 20.
The high-resolution DMUSIC spectrograms in panels(c) and (f) of Fig. 20 show the frequencies that are inducedby those phase space structures. A striking feature of these
spectrograms is that the evolution of the field’s oscillation
frequencies is far from smooth. The interpretation of these1403087-27

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 20 Propagation of EP phase space density wave fronts and nonadiabatic onset of chirping during the first dozen pulses of themarginally unstable case (A) (left) and first half dozen pulses in the strongly unstable case (B) (right). Panels (a) and (d) show
the time traces of the field amplitudes A(t). Panels (b) and (e) show the distance |ΔP̂ζ | of the hole wave fronts (blue) and clump
wave fronts (red) from the seed resonance P̂ζ,res = 0.718. The wave front positions in the δ f -weighed Poincaré data — which
can be seen indicated by dashed horizontal lines in the fourth snapshot of δ f (P̂ζ , ϑ) in Figs. 11 and 12 — were determinedusing an automatic threshold-based algorithm. Panels (c) and (f) show high-resolution spectrograms obtained with the DMUSIC
algorithm (see Appendix E.2), using a time window of size Δtwin = 0.094 ms in case (A) and Δtwin = 0.047 ms in the more
rapidly evolving case (B). The horizontal green bars indicate Δtwin for easy comparison with the spectral patterns. The ΔP̂ζ (t)curves from panels (b) and (e) have also been converted to frequency using Eq. (D6) as |Δν| ≈ |ΔP̂ζ | × 486.8 kHz, and then plottedas νfront ≈ 100.4 kHz ± |Δν| together with the spectrograms in (c) and (f). Vertical orange lines indicate the times of amplitude
minima, where phase jumps ±π occur. The dashed black curves in (c) and (f) represent fits of the ZC model in our parametrization(30) using the values in Table 4.

spectrograms requires caution (cf. Appendices B.1 and E),
but at least some features seem to be correlated with thewave fronts’ radial propagation in a pulsed manner. In
many instances, the frequency varies rapidly during an amplitude minimum between successive pulses, so that theearly chirping has a staircase-like appearance23 .
The spectrograms in Figs. 20 (c) and (f) also show
copies of the ΔP̂ζ (t) curves from panels (b) and (e), thathave been converted to frequency using Eq. (D6) as |Δν| ≈|ΔP̂ζ | × 486.8 kHz. One can see that the signals in the spectrogram — and, thus, the eﬀective pump wave frequencies — lie well within the range νfront = 100.4 kHz ± |Δν|
spanned by the wave fronts. The separation between |νfront |and the spectral peaks is smaller in the marginally unstablecase (A) than in the strongly unstable case (B). One reason
is simply that the radial width of the primordial holes andclumps in case (B) is much larger than in case (A) during
23 In case (B), the first few frequency steps in Fig. 20 (f) during the in-terval 0.50 ms  t  0.55 ms occur during pulse #1, so on a time scalethat is shorter by about a factor of 3 compared to beating and particle
bouncing. So far, we have not been able to attribute these steps to any
physical feature of δ f , so their cause remains unknown. Statistical noise
associated with the imperfect quiet start may play a role.

the onset of chirping.
In addition, the poloidal drift speed of the pump waves
relative to the seed wave is initially reduced due to stronger
resonant trapping in case (B), whose damping rate γd /ω0 =
4% is smaller than the 7.5% in case (A). Eventually, thepump waves in (B) escape from the seed resonance andchirping proceeds similarly to case (A). Thus, the onset of
chirping in case (B) is first retarded with respect to mode
growth and then occurs in a more abrupt manner than incase (A). This can be seen clearly seen in Figs. 20 (c) and(f): while upward chirp takes about 1 ms to reach 110 kHz
in case (A), it takes only about 0.1 ms in case (B). This
trend should continue with decreasing γd until, in the absence of damping (γd = 0), the spectral peaks remain at orvery near the seed frequency (i.e., no sustained chirping),
whereas ΔP̂ζ remains large and represents the half-widthof the phase space island as in Fig. B124 .
The pulsed radial propagation of the wave fronts andthe more or less discrete stepping in frequency in Fig. 20
during the nonadiabatic phase of the chirp may be regarded
24 Here, we have expressed once more in diﬀerent words what has already been mentioned near the end of Sec. 4.4 above.

1403087-28

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

as one possible realization of the relay runner paradigm,a phenomenological model proposed by Zonca & Chen
[8, 50] on the basis of the mode pumping mechanism discussed by White et al. [51]. In the original model, therelay runners were identified with the poloidal harmonics
and EP-driven distortions of a shear Alfvén wave packet intoroidal geometry, as illustrated schematically in Fig. 1 (d).
The physics of evolving mode structures are absent inour semi-perturbative model, but this caveat can be tolerated in the present cases, since we have long-wavelength
modes and large magnetic drifts. The successive (andphase shifted) pulses in our simulations can play the roleof the relay runners who are taking turns to carry the batons in the form of phase space density waves. Moreover,in the marginally unstable case (A), the first dozen or sopulses have gradually increasing amplitudes, which corresponds to the process of convective amplification [4], except that here the radial convection occurs only in phasespace, not in the mode structure. In addition, recall that
we have placed the peak of our field mode at a larger radius than the resonance (cf. Fig. 9). From the perspectiveof the clumps, which propagate radially outward towards
the mode’s peak, the field appears to get stronger with increasing radius, even when A(t) fluctuates around the samelevel. This process may be referred to as “preconvective
amplification” as discussed in Appendix D.1. Preliminary
results of parameter scans indicate that this arrangement
significantly aﬀects the preferred direction of chirping. In
the present setup, we observe stronger (more intense andmore rapid) upward chirping, which means radially outward, towards the peak of the mode.
The Zonca-Chen (ZC) model in the form presented inRef. [8] predicts that the chirping rate should first accelerate exponentially, then pass through a phase of linear-intime chirping, and eventually decelerate exponentially; inother words, it should evolve like a hyperbolic tangent: γ t 1exp(γ0 t)
0=
1 + tanh, (29)
δνZC (t) ∝1 + exp(γ0 t) 22
where γ0 is the linear growth rate, and the time variable t issuch that the chirping begins at t → −∞, and the transition
from accelerating to decelerating chirping occurs at t = 0.
The ZC relay runner model was formulated for magnetically deeply trapped toroidally precessing particles
in the short-wavelength limit with small magnetic drifts,
where the particles are passed on from one poloidal harmonic to the next. However, here we deal with a singleharmonic long-wavelength mode and EPs with large magnetic drifts, so we cannot simply adopt the coeﬃcientsfrom Ref. [8]. Instead, we use the parametrization
γΔν±∞ 
0(t − t0 ) ; (30)
1 + tanhδνZC = ν−∞ +22where the height Δν±∞ of the hyperbolic tangent is treatedas an unknown parameter that is fitted separately to theup- and downward chirping branches in our spectrograms.

Fig. 21 DMUSIC spectra for the first 5 ms in cases (A) and (B).
The spectral analysis was performed with a longer time
window than in Fig. 20; here Δtwin = 0.47 ms. Blue andred lines represent the time traces of the hole and clump
wave front frequencies νfront = 100.4 kHz ± |Δν|, as inFigs. 20 (c) and (f). The dotted black curves represent
the prediction of the BB model (12) with the actual values of γL and γd for the respective case, but a somewhat
arbitrarily chosen pole location: 1.0 ms in case (A) and
0.3 ms in case (B). The dashed black curves represent
fits of our parametrization of the ZC relay runner model(30) using the values in Table 4.

Table 4 Parameters for our adaptation of the Zonca-Chen (ZC)
relay runner model in Eq. (30).

Parameter
γ0 = γL − γd [kHz]t0 [ms]ν−∞ [kHz]Δν±∞ [kHz]Case (A)
Case (B)
3.14
21.99
1.65
0.55
100.4
100.2
+13 (clump), −10 (hole)

The location t0 of the inflection point is also fitted, whereas
the promptly shifted seed frequency ν−∞ = ν0 + Δν0 (Appendix D.4) and γ0 = γL − γd are measured in the simulations. The parameter values are given in Table 4 and thefitted curves are shown as dashed black lines in Figs. 20 (c)and (f). One can see that both in the marginally unstable case (A) and in the strongly unstable case (B), the timescale of the chirps is captured well by the tanh(γ0 t/2) function with the nominal growth rate values. The fitted curves
are shown once more in Fig. 21, which covers the first 5 ms

1403087-29

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)of the simulations and shows DMUSIC spectrograms obtained with a larger time window, Δtwin = 0.47 ms. For
comparison, the result of the BB model (12) is also shown.
After the hyperbolic tangent levels oﬀ, it eﬀectively
measures the width of the turbulent belt that has formed
around the resonance (although that width continues to
evolve). Here, we simply matched the amplitude Δν±∞of the hyperbolic tangent function to the observed spectrum, but a quantitative model based on readily available
parameters such as γL , γd , drift orbit width, field geometry
and mode structure would be desirable. DMUSIC spectrograms with high temporal resolution as in Fig. 20 show
only the dominant chirps and tend to follow the vortices. In
contrast, spectrograms obtained with relatively long time
windows as in Fig. 21 show the turbulent belt as a region
characterized by multiple fluctuating spectral peaks that
may be attributed to the nonadiabatic interchange dynamics of hole and clump waves of various sizes as seen in theregion around P̂ζ ≈ 0.7 of the last two snapshots in Fig. 11.
The relay runner model (29) does not capture the longrange chirps produced by the emission of solitary vortices.
In our simulations, chirping continues beyond the turbulent belt as massive solitary vortices detach and propagate
away from the seed resonance. The detachment process
and the role that beating may play there are examined in
the following subsection.

5.3

Growth, detachment and interference ofsolitary vortices

When we discussed phase jumps in Sec. 5.1 above, we
noted that, at the end of a pulse, the far halos tend to lag
behind their H&C wave fronts, so that the near halos tend
to be tilted to the right as indicated by the dotted lines inFig. 19 (a). At the beginning of the next pulse, the halos arelocated ahead of their H&C waves fronts, so that the nearhalos are tilted to the left as indicated by the dotted lines
in Fig. 19 (c). This has an important implication: Since theclump in Fig. 19 (c) moves into a halo region where the EP
density has already increased above the reference (δ f >0, represented by the halo’s red shade) and since the field
amplitude is growing at that time, the clump grows in size
as additional dense fluid gets trapped in its neighborhood.
The same holds in reverse for the hole. This leads us to the
topic of this subsection: the formation of massive holes andclumps and their detachment from the turbulent belt. The
detachment process is gradual and diﬃcult to quantify, so
this will be a mostly observation-based discussion of theresults summarized in Figs. 22-24.
We will be telling the story primarily from the pointof view of clumps because, in our setup, they move into
regions of higher field amplitude, so they tend to dominate the picture. Apart from this imposed asymmetry, thestory may be told in the same way for the holes because
here the only diﬀerence between holes and clumps is thesign of δ f . It must be kept in mind that packets of VlasovFig. 22 Schematic illustration of the formation of EP clump
fragments in the inner plasma (bottom) and their accumulation into a massive clump in the outer plasma (top).
The horizontal axis corresponds to the poloidal angle ϑ.
As in the phase space plots shown in Figs. 23 and 24,
we are in the frame of the electromagnetic (EM) seedwave, so structures above the seed resonance drift to theleft and those below drift to the right. The yellow-greenshaded arrows in the text boxes on the left indicate theaction of the EM field wave. The red arrows indicate themotion of clump fragments. The black arrows indicate
flow patterns in and around the reference clump.

fluid that are part of a clump at one time may be part of ahole at another time, especially when they are located inside the turbulent belt, where they may be advected up anddown between the belt’s boundaries and become clumplike or hole-like depending on whether the density at thatlocation was lower or higher initially, at t = 0. After a vortex has detached from the turbulent belt, its interior may
become a “pure” clump or hole, but its surface (boundary
layer) remains transformative when the field fluctuates (cf.
Fig. 1 (c)).
With the conceptual framework put in place, let us
now discuss the vortex formation and detachment process
in some detail. Figure 22 illustrates schematically how weunderstand the growth of a clump while it is still embedded in the turbulent belt. The turbulent belt in the center ofFig. 22 can be thought of as consisting of clump and holefragments of various sizes. Here, only the clump fragments
are indicated, which are gradually advected downhill; i.e.,
upward in Fig. 22, towards lower EP density and higher
frequency. The clump fragments originate from the lower
part of Fig. 22, where new clump fragments are flaked oﬀ
the EP density profile with every new pulse of the field (recall also Fig. 19). Whenever a clump fragment is flaked
oﬀ, its place is taken by a hole, which is shown in blue inthe lower part of Fig. 22. Animations of this process (Appendix E.3) make it appear as if a hole wave front is digging into the EP profile, advancing uphill with each new
pulse of the field as indicated by the wavy cyan-colored

1403087-30

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 23 Detachment of solitary vortices in the marginally unstable case (A). Panel (a) shows the Fourier spectrogram obtained withlarge time window Δtwin = 0.47 ms. The dotted green parabola is the prediction of the BB model (12). Panel (b) shows theevolution of the field amplitude A(t) during the interval 2.35 ms ≤ t ≤ 3.9 ms. Four sets of snapshots labeled (I)-(IV) show thestructure of the phase space density perturbation δ f before, during and after the detachment of a massive solitary clump vortex.
The δ f snapshots in set (II) have been aligned with the center of the uppermost clump, so the ϑ-window shown varies, moving
leftward from positive to negative values. Horizontal dashed lines indicate the locations of the hole and clump wave fronts, whose
evolution was shown in Figs. 20 and 21. The white and black circles and black cross symbols in the δ f snapshots indicate thelocations of particles whose trajectories we analyze in Figs. 27 and 28 below. Each particle appears approximately 4 times during
the τζ0 = 4.7 μs interval where the δ f Poincaré maps were accumulated.

arrow drawn towards the center of the hole structure.
While being advected downhill, the clump fragments
may merge or decompose as they are being continuouslyinterchanged inside the turbulent belt. Large clumps may
already form along the way, but in Fig. 22 we assume that
this occurs at the upper rim of the turbulent belt. The wavy
dark red arrow inside this massive clump structure indicates how it advances outward with each new pulse. At thesame time, neighboring clump fragments may be absorbed
into the larger structure.
Whether or not our reference clump repels or absorbs
a nearby fragment depends on the instantaneous relative
phases of three players: (i) the reference clump, (ii) the
fragment, and (iii) the field wave. The phase of the field
wave is indicated at the top of Fig. 22, where the black ellipses indicate the separatrix of the eﬀective phase spaceisland. O- and X-point locations are also marked. The
bold black arrows depict the direction of the vertical component of phase space flows in diﬀerent parts of the fieldwave. The corresponding pattern of phase space flow inthe boundary layer of our reference clump are indicatedby thin black arrows. The green-outlined clump fragmentsthat happen to be located just behind the center of our reference clump are attracted and can be absorbed, if they areable to reach the same altitude on time. Meanwhile, theyellow-outlined fragments that lie ahead (in the directionof the poloidal drift) are repelled by the downward flows at
this time. Fragments that are repelled at one time may beattracted and absorbed later. Overall the reference clump
grows due to a net upward flow of denser EP Vlasov fluid.
Snapshots of this accumulation phase as it occurs in
the actual simulation can be seen in Fig. 23 for case (A)and Fig. 24 for case (B), where they are labeled as stage(I). The accumulation process still continues for some time
during stage (II), which we consider to be the detachment
phase. The detachment of a clump or hole vortex fromthe turbulent belt is a gradual process, during which interactions with neighboring fragments become less intenseand strong interactions become less frequent. Note thatthe δ f snapshots taken during stage (II) still show a clear
correlation between the radial location of the wave fronts(horizontal dashed lines) and the field pulsations in panel(b). Thus, the fronts’ radial advancement still occurs in apulsed manner (recall Figs. 20 (b) and (e)), suggesting thatbeating plays a role in this process.
Although we cannot identify a precise timing, we consider the detachment to be more or less complete in stage(III). In stage (III) of the strongly unstable case (B), the1403087-31

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 24 Detachment of solitary vortices in the strongly unstable case (B). The results are arranged as in Fig. 23. Here, panel (b) shows
the time interval 0.71 ms ≤ t ≤ 2.35 ms.

large first-generation clump is subject to strong distortions, which can be attributed to the presence of a secondgeneration clump vortex that has grown large near the first.
These two structures continue to strongly interfere witheach other for some time. Finally, the snapshots for stage(IV) show the regime where the outermost clump and innermost hole can be regarded as completely detached vortices as they have separated far from and are hardly distorted by other clumps or holes.
The last snapshots in Figs. 11 and 12 (Sec. 4) werealso taken during the fully detached stage (IV), and thepoloidally averaged radial profiles of δ f for those snapshots are shown in Fig. 25. Here, the detached H&C vortices have the appearance of radially propagating solitary
waves. An interesting observation can be made in themarginally unstable case (A): In the wake of the solitary
clump vortex on the right-hand side of Fig. 25 (a), there isa depression in the EP density profile that we have labeled
zonal hole. The attribute “zonal” refers to the nearly featureless uniformity of that region, which appears in the lastsnapshot of Fig. 11 as a cyanish band that contains no obvious macroscopically coherent structures and separates thesolitary clump from the turbulent belt. The same is true forthe detached solitary hole vortex on the other side, which isseparated from the turbulent belt by a zonal clump that appears as a yellowish band in the lower part of the last snapshot of Fig. 11. The presence of a zonal hole implies thatthe solitary clump vortex has not merely plowed throughthe surrounding EP Vlasov fluid but also absorbed some ofit until the time of this snapshot (t ≈ 6 ms). Conversely, thesolitary hole vortex has deepened. In the next subsection,
we will see that the clump begins to disintegrate thereafter.

Fig. 25 Poloidally averaged EP density profiles for the last two
snapshots of Figs. 11 and 12 for cases (A) and (B). The
snapshot at t = 6.01 ms (tˆ = 1277.5) was taken after
solitary hole and clump vortices had detached and became well-separated from the turbulent belt. In case(A), a band of reduced density (zonal hole) lies behind the solitary clump, and a band of increased density(zonal clump) lies behind the solitary hole vortex (seealso Fig. 11). Note that the labels “solitary clump/holevortices” are not quite accurate: Since we are plotting the poloidally averaged density, the mean value ofF(Pζ ) = F0 + δ f ϑ at the radius of a hole or clump
vortex includes fluid elements at neighboring ϑ-values
that originate from diﬀerent radii. Of course, the local value of F(Pζ , ϑ) of each EP Vlasov fluid element
is conserved, as required by the Liouville theorem, but
this is not visible in the profiles shown here.

1403087-32

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)
shot in Fig. 11)25 .
The robustness of the solitary vortices that survive inspite of the persistent beating of the field wave (which is
caused by the interference between all the pump waves inthe system) indicates that, after detaching from the turbulent belt, their interiors may have attained a certain degree of adiabaticity, with more or less conserved actions
as assumed in so-called waterbag models (e.g., see Refs.
[14–16]). This can be confirmed by examining the motionof individual simulation particles, which constitutes the final part of this study.

5.4 Particle trapping and transport

Fig. 26 Instances of constructive interference causing largeamplitude spikes through spontaneous alignment ofphase space structures in the advanced stages of thestrongly unstable case (B), where multiple detachedsolitary hole and clump vortices are present. Panel (a)shows the Fourier spectrogram (Δtwin = 0.47 ms), withvertical dotted lines indicating the times where relatively large spikes in the amplitude A(t) are observed.
The spikes are shown in panels (b) and (c), and coincidewith the spontaneous phase alignment of multiple pumpwaves as shown in the snapshots of δ f in the lower partof the figure.

It must be noted that even after detachment, the interference between the spatially separated phase space structures has noticeable eﬀects. Besides ubiquitous beating,one manifestation of such an interference is the observation of relatively large spikes in the field amplitude that
occur when multiple phase space structures happen to havetheir phases aligned. According to Appendix A, we expect
large spikes when multiple clumps or multiple holes arealigned, and it is beneficial when the clumps are out-ofphase with respect to the holes, or simply more intense.
Two instances of such transient phase alignments can beobserved in Fig. 26. Although transient, such large spikes
may have an influence on the rate at which the hole andclump vortices advance radially and accumulate (or lose)material. Similar processes may be relevant for triggering
Abrupt Large-amplitude Events (ALE) [36] whose precise
timing, according to our current understanding, depends
sensitively on the trajectories of multiple players, just likethe events of spontaneous phase alignment in Fig. 26.
It is remarkable that the outermost solitary clump vortex in Fig. 26 retains its structure with p = 4 elliptic points
for a long time even after entering the domain of the p = 5
resonance; namely, the region P̂ζ  0.6 (cf. the first snap-

The δ f snapshots shown in Figs. 23 and 24 contain
black and white circles and black crosses, which indicate
the locations of three selected simulation particles. During the τ0ζ = 4.7 μs time interval used to accumulate
δ f Poincaré maps, each particle appears approximately 4times at diﬀerent ϑ. In Fig. 23 for the marginally unstable case (A), one can see that the particles represented by
the black and white circles are part of the massive clump
structure at the time of snapshot (I), whereas the particle
represented by a black cross is absorbed by that clump atsome later time.
The time traces P̂ζ (t) for these three tracer particles
in case (A) are plotted in Fig. 27 (a), where one can seethat the black cross, here labeled as tracer #3, is trappednear the clump’s boundary at t ≈ 3 ms, then transported
outward with the clump for several milliseconds, before
being released and deposited at a larger radius at t ≈ 7 ms,
indicating that the clump has begun to shrink. The timet ≈ 6 ms where the clump begins to shrink coincides withthe time where the clump enters the domain P̂ζ  0.6,
where the presence of the next resonance with p/n = 5/5
becomes noticeable in the δ f halos (see the first snapshot
in Fig. 11). That timing also coincides with a reductionin the clump’s radial propagation speed (= chirping rate),
which had been remarkably linear at earlier times, as canbe seen more clearly in the inset panel in Fig. 27 (a)26 .
The time traces P̂ζ (t) for three tracer particles in case(B) are shown in Fig. 27 (c). Here the rate of radial propa25 The vortex is even able to travel beyond P̂ = 0.55, where no simuζ
lation particles have been loaded. Although not shown here, we find that
this initially empty region of phase space is subsequently filled by particles that leak from the massive clump vortex. However, the δ f simulation
becomes invalid in that region because outward propagating Vlasov fluid
elements can no longer be replaced with inward propagating ones. When
this boundary eﬀect feeds back on the phase space marker trajectories, we
expect that it breaks the Liouville theorem. Our sensitivity tests indicatethat this eﬀect remains tolerable for the duration of our simulations.
26 The linearity of the first half of the long-range chirp in Fig. 27 (a)does not seem to be a robust feature of the simulation. The occurrence,
form and number of long-range chirps is sensitive to details. They canvary significantly when the physical parameters (such as damping anddrive) are changed slightly, and they are even aﬀected by numerical parameters. This can be seen by comparing the spectrogram of case (A) in
Fig. 11 (b) with those in Figs. D1, D10 and D12 (a) in Appendix D. All
these simulations were performed with the same values of γL and γd as
in case (A).

1403087-33

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 27 Time traces of the radial position P̂ζ (t) and rotating frame energy E  (t) of three tracer particles labeled #1, #2 and #3 in case (A)(top) and case (B) (bottom). One tracer is trapped in the center, one in the bulk, and one near the boundary of the first-generation
solitary clump vortex. The inset in panel (a) shows a portion of the data tilted such that the dashed line becomes horizontal tohighlight the constancy of the clump’s radial propagation in that time window. The zoomed-up inset in panel (c) is meant to show
more clearly the rapid modulation caused by the field’s beating. For the reader’s convenience, the approximate bounce periods
τb during some stages of the evolution are shown as text labels. See also Fig. 28 for the full time traces of νb (t) for each tracer
particle. The vertical axis on the right-hand side of panels (a) and (c) have been converted to frequency ν using Eq. (D6).

gation (and chirping) can be seen to decrease continuously.
As in case (A), the clump begins to lose particles as soon as
it enters the domain P̂ζ < 0.6. First, the outermost tracer#3 (orange curve) is detrapped at t ≈ 8 ms, followed bytracer #2 (red curve) at t ≈ 10.5 ms. The innermost tracer#1 (blue curve) remains trapped until the end of the simulation. Note how tracer #3 gets temporarily trapped by

the second-generation clump around t ≈ 10 ms and by thethird-generation clump after t ≈ 13 ms, as these follow
closely behind the first clump. As it happens, the fate oftracer #3 was that the secondary trapping event moved thisparticle radially back inward.
These observations show the extent of global transportcaused by the propagating solitary phase space vortices,

1403087-34

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. 28 Evolution of the bounce frequencies νb of tracers #1, #2 and #3 whose time traces P̂ζ (t) were shown in Fig. 27. For orientation andcomparison, panels (a) and (e) show again the chirps observed in case (A) (left) and case (B) (right) using a Fourier time windowof size Δtwin = 0.47 ms. Panels (b)-(d) and (f)-(h) show Fourier spectrograms of P̂ζ (t) using Δtwin = 3 × 0.47 ms = 1.41 ms.

which bears similarity to the physical picture of so-called“bucket transport”, with the addition that — besides the
chirping φ̈ [6] — the pulsations of the field amplitude A(t)assist with the transmission of untrapped particles aroundthe bucket by widening the nonadiabatic boundary layerthat would constitute only a singular separatrix in the limitof a field wave with fixed amplitude and phase, and an isolated single resonance (see Fig. 1).
Another closely related and important piece of information that we can gather from Fig. 27 is the following.
In both cases (A) and (B), one can see that particles thatare located near the vortex center (namely, tracer #1) remain there for the entire duration of the simulation. Similarly, particles in the bulk of the vortex and particles in theboundary layer remain there, with the peripheral particles
being the first to escape when the vortex begins to shrink.
This means that these vortices maintain a nested structure.
What we find particularly remarkable about this observation is that the nested structure of the vortices is maintained in spite of the persistent beating of the field wave,
and in spite of the fact that the dynamics in case (A) areonly marginally adiabatic. Let us elaborate these important points in some more detail.
The eﬀect of the persistent violent fluctuations of the
field in the form of amplitude pulsations and phase jumps
can be clearly seen in Fig. 27 in the form of high-frequency
modulations of all time traces P̂ζ (t). Typically, we see between 5 and 7 oscillations per bounce period τb , so thetime scale separation is not all that large27 . Moreover, inthe case of tracer #1 (blue curve), which is located near theclump center in both cases (A) and (B), the perturbations
caused by the beating have amplitudes that are comparable
to the radial excursion of the particle due to its revolution
around the vortex center. In Fig. 27 (a) it is sometimes evendiﬃcult to see the bouncing signal of tracer #1 within the“beat noise”. Although we are looking only at very few
particles, the overall integrity of the vortices suggest that
other particles behave likewise. The possibility of maintaining such a nested structure in the presence of strong
perturbations is fascinating and strengthens the stance of
reduced adiabatic models such as the “waterbag” [14, 15];
in particular, when cast in a generalized form, allowing thestructure to shrink or grow [16, 25]. This implies the presence of a nonadiabatic boundary layer, even when it does
not appear explicitly in the idealized theoretical model.
Thus motivated, let us examine the degree of adiabaticity
in our vortices.
Although it is possible to determine the characteristic bounce periods τb by-eye when looking closely at thetime traces in Fig. 27, the evolution of the bounce frequencies νb of our tracers can be conveniently seen at a glancefrom the Fourier spectra that are summarized in Fig. 28.
With Δtwin = 1.41 ms, we have used a fairly large sliding
time window for the Fourier transforms, but one can never27 The inset in panel (h) of Fig. 28 shows the spectral peak of the beatsignal in case (B) at about 20 . . . 25 kHz, whereas νb ≈ 3 . . . 7 kHz.

1403087-35

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

theless see significant variations of the bounce frequencies
over time. Various observations that are interesting but not
immediately relevant here are discussed in Appendix C.2.
What we would like to emphasize here is how the bounce
frequency compares to the chirping rate: The bounce frequency near the clump center, namely for tracer #1 inFigs. 28 (b) and (f), is approximately νb ≈ (4 . . . 6) kHz
during the first few milliseconds after the clump detaches.
The chirping rates in that part of the simulation can bereadily inferred from Fig. 27, where the vertical axis on
the right-hand side shows the radial position P̂ζ converted
to frequency ν using Eq. (D6). In case (A), the slope is1
ms−1 (as shown explicitly in the inset panel),δP̂ζ /δt = 39kHzso we have δν̇ = 486.8= 12.5 kHz/ms. In case (B), we
39 ms
have roughly δν̇ ≈ 20 kHz/2.5 ms = 8 kHz/ms. With this,the adiabaticity parameter δν̇/νb2 in Eq. (11) becomes⎧ 12.5 kHz/ms⎪⎪⎪⎨ (4.5 kHz)2 ≈ 0.6 : Case (A),(31)δν̇/νb2 = ⎪⎪⎪⎩ 8(6kHz/ms≈ 0.2: Case (B);kHz)2where the overbar (. . .) means that we are looking onlyat the overall trend of the chirp, ignoring rapid fluctuations in δν̇, such as phase jumps. Equation (31) highlights in a quantitative manner what may already havebeen obvious to experienced eyes from the time traces inFig. 27: Our strongly unstable case is more adiabatic thanthe marginally unstable case (A), with the latter being only
marginally adiabatic since its δν̇/νb2 value is not far fromunity. In fact, if it had not been for the observed nested
structure, one may have been tempted to classify case (A)as nonadiabatic. Instead, it seems that an adiabatic treatment of the detached solitary clump vortices (and presumably also the solitary hole vortices) can be justified in bothcases that were analyzed here, near and far from marginal
stability, and with fairly strong drive γL /ω0 ≈ 8%.
Finally, we examine to what extent the simulation particles depart from their initial position E  = 75.6 keV in
the rotating frame energy. Figure 29 shows that the particle distribution becomes slightly crescent-shaped but remains close to the initial E  = 75.6 keV line during the first≈ 4 ms in case (A) and during the first ≈ 2 ms in case (B).
After that, the distribution broadens and reaches a widthof about δE  ≈ 2.5 keV in both cases. The particle clusters in these broadened energy bands can be clearly associated with the upward propagating clumps. A more detailed
view of the distribution’s 4-D structure in (δ f, P̂ζ , E  , ϑ)and (P̂ζ0 , P̂ζ , E  , ϑ) space, and how it evolves in time isgiven in Figs. C1 and C2 of Appendix C.3.
Time traces of E  (t) for a few tracer particles areplotted in panels (b) and (d) of Fig. 27. For case (A),Fig. 27 (b) shows that E  is well conserved during the first
4 or 5 bounces. Oscillations, whose frequency seems tomatch the bounce frequency νb , become clearly noticeable at about t ≈ 3 ms, which coincides with the stage
where the first massive clump starts to detach from the turbulent belt (cf. Fig. 23). Starting around t ≈ 5 . . . 6 ms,

Fig. 29 EP transport in the E  -P̂ζ plane, where E  = E − ω0 P̂ζ /nis the particle energy in the frame moving with the seedwave. Several snapshots of the particle distribution are
shown for case (A) (top) and case (B) (bottom). Moredetails can be found in Appendix C.3, Figs. C1 and C2.

where the clump begins to shrink and interact with the
p/n = 5/5 resonance, we observe nonsinusoidal oscillations with a magnitude on the order of 1 keV and long periods of about 3 . . . 4 ms, the cause of which is still unknown.
In Figs. 27 (b) and (d), there seems to be an overall trend
for E  to decrease in time, and the crescent-shaped bendingof the full distribution as seen in Fig. 29 is probably related
to this trend.
Note that the magnitude of high-frequency beat noise
is smaller in the time traces of E  (t) than in the time tracesof P̂ζ (t), so that the curves in panels (b) and (d) of Fig. 27
appear significantly smoother than in panels (a) and (c).
This implies that there is a significant cancellation between
beat-induced fluctuations in E and P̂ζ . This observation
leads us to conclude that beating may not be the primarymechanism that leads to the decoupling of E and P̂ζ and,thus, the breaking of the E  = E − ωPζ /n = const. condition. Therefore, we suspect that the solitary vortices themselves may be the primary cause for the perturbations in E  ;in particular, their radial propagation (chirping and plowing through surrounding material) and direct mutual interactions (e.g., when their nonadiabatic boundary layers in-

1403087-36

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

terfere or overlap).
The E  values of the three tracers in Figs. 27 (b) and(d) vary less in case (B) than in case (A), but this trend is
not seen in Fig. 29, where the full particle distribution isscattered more broadly in case (B) with weaker damping
than in case (A) with stronger damping. Presumably thisis due to the larger number of more closely spaced vortices
in case (B). If conditions were equal, we expect that thevariations in E  should decrease with decreasing damping
rate, because γd = 0 yields eﬀectively constant E  as canbe confirmed in Fig. B3 of Appendix B for case (A0).
The E  line broadening seen in Fig. 29 has been reproduced with 2.5 times shorter time steps, which gives better
energy conservation, so we may assume that its cause is nota numerical one. In particular, the crescent-shaped distortion that dominates during the first few milliseconds of thesimulations seems to be a robust physical feature. We suspect that this (small) modification of E  is a consequenceof chirping, as it aﬀects the domains of holes and clumps
alike. However, this hypothesis remains to be tested.
At present, a systematic limitation of our simulations
prohibits a meaningful study of the observed E  line broadening. Recall from Fig. 10 that the phase space structures
have the form of elongated cylinders in (Pζ , E  , ϑ)-space,which we have eﬀectively truncated by loading particles
only along a line with E  ≈ E0 = 75.6 keV = const. Thismeans that our simulations using the δ f method do not accurately represent dynamics in the E  direction, where nomarker particles have been loaded, because the continuityof GC phase space and, thus, the condition dG/dt = 0 underlying Eq. (18) is violated. This does not immediately
make the entire simulation unphysical, but it can produce
artifacts, and the “ghost chirp” in Fig. 2 (b) may be one
such artifact (see Appendix C.4). For the scope of thepresent work, this does not seem to be a problem. Theline broadening reaches only δE  /E  ≈ 3% in 10 ms, andwe show in Appendix D.2 that similar dynamics are obtained when marker particles are loaded in a band of width
ΔE0 = 7.5 keV that fully encompasses the broadened distribution in Fig. 29.

6. Summary, Discussion and Outlook
When multiple field fluctuations with the same polarization are simultaneously present in the same spatial domain of a plasma, their linear superposition causes beating in the combined signal s(t) = A(t) sin(−ω0 t − φ(t)),
here measured relative to a reference wave with frequency
ω0 . The beats consist of pulsations in the amplitude A(t)and jumps by ±π in the phase φ(t). The pulse length of
the beats, their magnitude, and the time scale of the phasejumps depend on the separation between the frequencies
and the relative amplitudes of the interfering waves.
The amplitude and phase modulations associated withbeating may contain useful information about possible
nonlinear couplings between the interfering waves. Thishas been exploited, for instance, in studies of coexistingEP-driven Alfvén modes and rotating kink/tearing modes[52–54]28 . Since diﬀerent modes generally have diﬀerent spatial structures, such cases exhibit not only temporalmodulation but also spatial modulation. The spatial interference may add to the complexity of the dynamics.
Beating can occur as a passive side-eﬀect, merely producing spectral sidebands. However, the beats of the field
can also play an active role in the dynamics through direct
or indirect feedback on the particles and/or field perturbations. One prominent example is spectral line splitting
in regimes where chirping is suppressed by strong scattering [32, 33]. Beating is an integral part of that process,where spectral broadening occurs in a well-organized fashion through limit cycle bifurcations (see also the introductory discussion in the last paragraph of Sec. 2.7).
In the present study, we examined the role of beatingin cases with strong chirping in a “collisionless” plasma;i.e., in the absence of drag or scattering. In contrast
to the discrete line splitting seen in the above-mentioned
scattering-dominated cases, scenarios with strong chirping
exhibit a high degree of temporal incoherence, as a large
portion of the phase space tends to be dominated by transient convective dynamics. Coherent vortex structures tend
to form only after a significant delay (millisecond scale).
Our beats can be viewed as the linear modulation of
one field mode that is driven by multiple sources (density
waves) in the phase space of one particle species (here,energetic deuterons). Although the fundamental process
underlying this phenomenon is the exchange of energy between particles and field waves, the observations are perhaps most easily interpreted as a collective eﬀect; namely,as the result of beating between multiple pump waves inphase space density that jointly drive the field oscillation.
This physical picture was introduced in Sec. 2.2 and supported by the mathematical analysis in Appendix A. The
purpose of this study was to clarify how the beating of the
field acts back on the pump waves, and what role beating
plays in the formation and evolution of solitary vortices
in phase space, which are responsible for long-range frequency chirping and associated particle transport.
The system considered was an ideal incompressible
electromagnetic flute mode in realistic tokamak geometry and with an initial frequency in the shear Alfvén band
(ν0 = ω0 /(2π) ∼ 100 kHz) that interacts with circulatingenergetic deuterons with kinetic energies around 80 keV.
The use of a semi-perturbative model, namely a fixed spatial mode structure, constrains the validity of our numericalanalysis to the long-wavelength limit, namely low toroidal
mode numbers n ∼ O(1), whose poloidal harmonics tend tohave a large radial width and would undergo only relatively
little distortion. Such low-n modes are typically driven by
28 When the frequency diﬀerence is large, as for high-frequency Alfvén
and low-frequency tearing modes, the slowly evolving modes may also
be treated as a stationary non-axisymmetric perturbation of the reference
state (or equilibrium) around which the high-frequency modes oscillate.

1403087-37

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

energetic particles (EP) with large magnetic drifts. Although the phase space structures responsible for chirping remain within the range of these magnetic drifts, thedrifts’ large radial extent makes the chirp-induced transport global and, therefore, practically relevant.
In the following subsections, we summarize the main
insights that we have gained in this work and discuss possible implications and directions for further study. Concluding remarks are made in the last Sec. 6.6.

6.1

Halos reveal phase relations and may
fuel the growth of holes and clumps

We demonstrated the existence of robust nonresonant
modulations in EP phase space density δ f that we call halos. They are visible as radially elongated stripes of alternating sign (δ f ≷ 0) in δ f -weighted kinetic Poincaré plots
accumulated over a short interval of one toroidal transit
time, here τζ0 = 4.7 μs. These halos are interesting forseveral reasons, which we summarize here.
In order to correctly interpret the phase space dynamics, it is often essential to compare the phase space structures seen in contour plots of δ f with the instantaneous
flow lines of the EP Vlasov fluid. In principle, the flows are
determined by the instantaneous amplitude A(t) and phaseφ(t) of the field, but reconstructing the flow lines from that
information is diﬃcult in realistic geometry and with large
magnetic drifts. For instance, when trying to estimate locations of X- and O-points along ϑ at a certain time t, wecannot simply shift the initial field pattern forward in timeusing the relation Δϑ = φ(t)/m, because the magnetic drifts
turn that relation into something more complicated. Conventional (long-time) kinetic Poincaré plots are not applicable in a rapidly varying field, and inferring the locationsof eﬀective O- and X-points through visual inspection ofnear-resonant δ f patterns also becomes diﬃcult after the
first beat, because those patterns become complicated.
Our observations indicate that far oﬀ-resonant halos(Fig. 13 (b)) are strongly correlated with the instantaneous
phase of the field, so they can be used to track the eﬀective
locations of O- and X-points. On that basis, it is possible todetermine at least the radial direction of the flows at eachpoloidal angle. Here, this technique has been used mainly
during the early stages of a chirping simulation, where solitary vortices have not formed yet. It helped us to estimate
flow lines in phase space and interpret the processes associated with the onset of beating and chirping. The eﬀect ofmore or less abrupt phase jumps between successive pulsesof the beating field can be observed clearly, and the halos
have revealed details about that process (Fig. 17).
Besides being an indicator of the phase relation between field and density waves, we suspect that the halos
may also play an active role in the dynamics: The phase relation between halos and density wave fronts is such that,during the beginning of a pulse, a halo of increased density (δ f > 0) lies ahead of the clump wave, and a halo of

reduced density (δ f < 0) lies ahead of a hole. We haveargued based on Fig. 19 that this may facilitate the inflation of these structures, making them grow larger with eachpulse of the beating field.
Although we can explain halos as collective density
modulations, the reasons for why they exist in this particular form and why they evolve in the way they do remains
to be understood. On the one hand, the halos are attachedto the hole and clump wave fronts, so we usually see twocounter-propagating halos. On the other hand, the halos areclosely linked to the phase φ of the field. In that respect,
the halos are similar to the equilibrated stationary modulations of phase space density in Fig. 4 (b), which are also
phase-locked collective structures, not tied to the motionof individual particles. The close connection with the field
phase is also manifested in the fact that halos can moveahead of or lag behind the hole and clump wave fronts thatthey are more or less attached to. Overall, the halos seemto be a manifestation of the feedback between collective
phase space structures and the field. It would be interesting to take this “explanation” from the philosophical to thescientific level.

6.2 Feedback loop for staircase-like nonadiabatic onset of chirping
One of the most pronounced eﬀects of beating wasobserved during the nonadiabatic onset of chirping during the first few milliseconds of our simulations (Fig. 20).
Based on our interpretation of the field dynamics as being a
consequence of a superposition of multiple pump waves in
EP phase space, we have identified the following feedback
mechanism (Sec. 5.2, Fig. 19):
• The emergence of two pump waves with frequencies
ω+ and ω− causes the field to beat, with amplitude
pulsations of order 100% and phase jumps by ±π.
• These beats drive the pump waves further away fromthe seed resonance by making the eﬀective phasespace islands (= flow contours) pulsate with alternating eﬀective O- and X-point positions.
• The resulting increase in the frequency diﬀerenceΔω = ω+ − ω− causes more beats per time. Moreover, new pump waves (hole-clump pairs) are generated with each beat.
• The increased radial transport also pumps more energy into the field perturbation, whose growth increases the size of the islands (= steeper flow contours) and the step size of the radial displacements.
Initially, this feedback causes an exponential accelerationof the up- and downward chirps ω+ (t) and ω− (t). Highresolution spectrograms obtained with the DMUSIC algorithm indicate that the frequency advances in a somewhat
discontinuous fashion (Fig. 20)29 . The staircase-like on29 Qualitatively similar behavior known as “parametric ladder climbing” has been reported in other systems [55].

1403087-38

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

set of chirping is consistent and often correlated with thepulsed radial expansion of the pump wave fronts. However, the increasing frequency diﬀerence Δω also leads toa shorter pulse length τpulse = 2π/Δω and reduced energy transfer to the field. The shorter pulses and saturating
field amplitude yield smaller steps in frequency, the consequence being a deceleration and saturation of the chirps.
Exponential acceleration followed by exponential deceleration is the nature of a hyperbolic tangent. In a simplified limit of their general theory [1], Zonca & Chen haveconstructed a phenomenological “relay runner” model [8],that predicts nonadiabatic chirping of the formδω(t) ∝ tanh(γ0 (t − t0 )/2);(32)

where t0 is the point of inflection and γ0 = γL −γd the initial
growth rate (linear drive γL minus field damping γd ). This
formula matches the overall trend and time scale of theearly staircase-like nonadiabatic chirps in our simulations,both in the marginally unstable and strongly unstable cases(Figs. 20 and 21). The magnitude of the chirping range had
to be fitted with a free parameter, because the original relay runner model was formulated for a diﬀerent scenario
than ours: The model coeﬃcients in Ref. [8] were derived
for magnetically deeply trapped EPs and short wavelength
modes, where the batons in the form of resonant pumpwaves in phase space density are relayed by neighboring
poloidal harmonics of a toroidal Alfvén wave packet ina process of convective amplification [4], as illustrated inFig. 1 (d). In the future, it would be beneficial to formulate
reduced quantitative models of nonadiabatic chirping forvarious scenarios.
In our case, there was only a single long-wavelength
mode (m/n = 6/5) interacting with circulating EPs that
perform large magnetic drifts. The relay runners were realized here by successive pulses (beats) of the same modewith alternating phase φ. In close proximity to marginal
stability (γL − γd  γL ), namely in case (A), the first saturation of the resonant instability occurred at a very small
amplitude (Fig. 14), and successive pulses gradually rose
to higher levels (Figs. 17 and 20). This can be regarded
as a proxy for the above-mentioned convective amplification process, except that in our semi-perturbative model
the convection occurred only in EP phase space, without
shifting the peak of the field mode. However, from the perspective of the clumps, the local field did appear to be “preconvectively” amplified, because we placed the seed resonance oﬀ-peak, so that the clumps’ radial propagation tookthem into regions with a stronger field. See Appendix D.1
for further discussion on nonperturbative eﬀects.

6.3

various scales. Although complicated in their detailed appearance, the phase space structures exhibit a certain degree of radial stratification (e.g., see the last snapshots in
Figs. 11 and 12). Considering that the overall dynamics
are fed by background gradients (generally, in combination with sources and sinks), one may describe the overall
transport of material and energy with equations that govern
only the zonal component and ignore the angular dependence. That idea has been realized in the theory by Zonca& Chen [1, 4], a specific form of which describes the evolution of so-called “phase space zonal structures (PSZS)”.
A quantitative comparison with the transport occurring in
the turbulent belt in our simulations would be interesting.
One complication that might cause discrepancies is thatthe assumption of two disparate spatial scales (backgroundplasma and mode structure) may be violated in our case,since we deal with long-wavelength instabilities comparable to system size.
Generally speaking, a representation in terms of PSZS
only as illustrated in Fig. 1 (d), without any angular dependencies, would yield reliable quantitative predictions
for transport of particles and energy only if suitable transport coeﬃcients can be calculated. Obviously, these coefficients will be very diﬀerent in strongly turbulent regimes
involving many modes, in cases with intermittent solitonlike avalanches, and in scenarios with propagating vortices subject to a beating field mode as we have simulated.
In fact, the eﬀective transport coeﬃcients may even vary
within one scenario. For instance in our case, transport
coeﬃcients are likely to diﬀer between the near-resonant
turbulent belt and the inter-resonant region occupied by detached solitary vortices. The results reported in this work
may contain information that can be used in the construction of suitable transport coeﬃcients for reduced models,such as the PSZS paradigm.
From an observational perspective, the formation ofmassive vortical clump structures involves a gradual aggregation of smaller fragments that flake oﬀ on the highdensity side of the turbulent belt and are advected downhill
towards the clump wave front (see the cartoon in Fig. 22),with the reverse being true for holes. Beating plays a keyrole in this process and a feedback loop similar to that identified during the nonadiabatic onset of chirping (summarized in Sec. 6.2 above) is present here as well:

Zonal and vortical structures

In the wake of the pump wave fronts that are observed during the nonadiabatic onset of chirping, a turbulent belt forms, whose dynamics are reminiscent of radially sheared convective interchange and wave-breaking on
1403087-39• Radially sheared poloidal motion of pump waves(holes and clumps) causes the field to beat.
• The beating field causes the hole and clump wave
fronts to advance radially in a pulsed manner andspawns new clump and hole fragments with eachpulse. The broadening of the spectrum and increasing number of pump waves causes more (and more
complicated) beats.
• By driving radial convection of the fragments, thebeats also facilitate their accumulation into “massive” (wide and deep) holes and “massive” (wide andPlasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

dense) clumps.
In other words, the beating enhances mixing between
radially stratified layers of pump waves (hole and clump
fragments), with net transport occurring due to the initial
density gradient. One the one hand, the stronger mixing
tends to inhibit the early detachment of vortical structures.
On the other hand, it can fuel the growth of massive holes
and clumps, especially at the boundaries of the turbulent
belt. The larger value of |δ f | inside a massive structure
makes it a more powerful pump wave (cf. Eq. (2a)). When
driven by a few dominant massive structures, the beats ofthe field may become more regular and more intense. Ultimately, this may facilitate the detachment of a massive
clump or hole, because it is likely that some material will
have to stay behind during the detachment. A massive vortex structure is more likely to remain intact and detachedafter loosing some material than a smaller one, because it
has more influence over the auto-resonant field oscillations(due to larger integrated |δ f | in Eq. (2a)), thus, perpetuating
its own stability by prolonging phase locking.
Although a higher degree of on-average phase locking
to the field makes a vortex more robust, too much phaselocking may lead to enhanced resonant particle trapping
and, thus, inhibit chirping. Here, phase slippage is facilitated by strong field damping. In addition, beating (due tothe presence of other pump waves) and chirping [6] bothcontribute to the maintenance of a relatively wide nonadiabatic boundary layer (Fig. 1 (c)), which prevents the vortices from falling into a state of adiabatic stagnation.
We also observed the formation of truly zonal structures (indicated in Figs. 11 and 25). After the first solitary
clump vortex in the marginally unstable case (A) had departed from the turbulent belt, it left behind a region of reduced density that we have called “zonal hole”. Similarly,
a “zonal clump” was observed behind the first solitary holevortex on the other side of the turbulent belt around theseed resonance. No such structures were visible in the
more strongly unstable (less damped) case (B). It may be
worthwhile to investigate in a future case study under what
conditions such zonal holes and clumps form and what role
they play. For instance, one may speculate that such zonal
modulations of the EP density profile play a role in the formation and radial propulsion of subsequent generations ofhole and clump vortices. Meanwhile, the leading solitaryvortices that have produced those zonal perturbations will
experience a reduced EP density gradient across their radial width. This may contribute to the reduction of theirradial propagation speed and, thus, their chirping rates.

6.4

Robust nested vortices in the presence ofstrong beating

Even before complete detachment, the internal structure of hole and clump waves was found to resemble, on
average, concentric nested layers of circulating flows. Especially in the marginally unstable case (A), this obser-

vation was somewhat surprising, because the adiabaticity
parameter defined in Eq. (31) has a value close to unity
in this case: δν̇/νb2 ≈ 0.6. It would not have been surprising to see this case dominated by nonadiabatic dynamics, where long-lived vortices do not form and where phasespace structures tend to be of collective nature, in the sensethat diﬀerent particles would constitute these structures at
diﬀerent times. In such a scenario, there would be no ballistic “bucket” transport.
Moreover, we expected that the amplitude pulsations
and phase jumps associated with ubiquitous beating would
scramble the interior of any apparently coherent phasespace structures and, if not prevent their formation, at least
reduce their life time. Therefore, we were fascinated by
the fact that prior simulations of such chirping systems(with varying degrees of complexity) routinely show thepresence of long-lived phase space vortices. In particular, our curiosity was raised by the fact that the observed
beats of the field are characterized by instantaneous growth
rates γ = d ln A/dt = Ȧ/A and phase shift rates φ̇ that canbe comparable to and even larger than the oscillation frequency of the seed wave,
 
max |γ| ∼ max φ̇  ν,

(33)

which seems to contradict assumptions like |γ| ∼ |φ̇|  νthat are attached to theories of adiabatic chirping; namely,
the separation of time scales needed to construct adiabatic
invariants. In the theoretical derivations, this problem had
been sidestepped by considering only one solitary vortex
structure at a time (see Eq. (23) in Ref. [24]). By ignoring
explicitly the interference between multiple pump waves,beating is taken out of the picture, leaving only slowly
evolving amplitudes and phases.
One of our main results, reported in Sec. 5.4, is that we(re)confirmed the robustness of the phase space vortices inour simulations and showed that, on average, their interiors
maintain a nested structure (i.e., good adiabatic constantsof motion), in spite of persistent strong beating.
One may argue that this is not surprising since thelargest and fastest changes in A(t) and φ(t) occur at timesof destructive interference, when A(t) is small and, therefore, may hardly aﬀect the particle motion. In addition, untrapped phase space structures in a collisionless plasma canbe destroyed only by phase mixing, and one might thinkthat the phase mixing times would be much longer than the
typical beat periods. However, we did not find order-ofmagnitude estimates suﬃciently convincing, because forrealistic tokamak geometry and with realistic parameters
for the particle motion and field oscillations, the separationof time scales is often vague, and it was certainly vague inthe cases we have simulated here. For instance, Fig. 7 (a)shows that the pulse period τpulse associated with the beatsis only a couple of wave periods T 0 = 2π/ω0 long, and thesame counts for the ratio of τpulse and the nonlinear bounce

1403087-40

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

times τb in Fig. 27: in our cases, we typically haveτpulse /T 0 ∼ τb /τpulse ≈ 5 . . . 7,(34)

which is less than an order of magnitude and numericalfactors of “order unity” (like π) begin to matter. Moreover,for our largest solitary clump vortices (Fig. 11), whose
width ΔP̂ζ ≈ 0.25 corresponds to a frequency diﬀerenceΔν = 12 kHz ≈ ν0 /8 (Eq. D6) between the vortice’s innerand outer edge, has a phase mixing time (cf. Sec. 4.1) thatis even comparable to the beat pulse period,
τpulse Δν  1.

(35)

On top of that, the time traces of the canonical toroidal angular momentum (radial position) P̂ζ (t) for individual particles in Fig. 27 demonstrated that the “beat noise” can becomparable to and even exceed the magnitude of the “signal” of the nonlinear bouncing motion, especially for particles that are (on average) trapped deep inside the field’s
eﬀective potential well; i.e., near the center of the vortex.
These qualitative and quantitative observations givethe impression that the phase space vortices in our simulations live in a fairly hostile environment. Therefore, we
think that it is remarkable that, even under these conditions with (at best) vague scale separation, the vortices remain robust and (on average) maintain their nested structure over the course of many milliseconds (a hundred beats
and more). This strengthens the stance of “semi-adiabatic”
models. Our results confirm the applicability of the idea
underlying the “bucket transport” model, with the addition
that, besides the chirping φ̈ [6], the beating of the field assists with the transmission of particles through the nonadiabatic boundary layer. Similarly, our results support extended “waterbag” models that include an active boundary
layer as in Ref. [16], where the eﬀect of multiple interfering long-range chirps is also considered.
If one of the solitary clump vortices would be able
to gain full control over the field, the resonance and theassociated eﬀective phase space island would be centeredaround that clump permanently. We expect that, in sucha case, the field would stop beating and we should obtain a single smooth long-range chirp. Such non-beating
chirps were also observed in JT-60U experiments as shownin Fig. 14 of Ref. [12]. With δν̇ ≈ 2 kHz/ms, those lone
chirps are 4 . . . 6 times slower than those studied in thepresent work. Eﬀorts are currently underway to reproduce
such a single (non-beating) long-range chirp using ORBIT
with JT-60U equilibria. Those plasmas were driven by onaxis negative-ion-based neutral beams, so it is likely thatthe seed resonances are located near the axis and nonstandard (stagnation) orbits contribute. Such boundary eﬀects
are often important in real-world problems, but have been
avoided in the present study.

6.5

Chirp-induced resonance overlap

The solitary clump vortices in our simulations wereseen to travel far enough to reach the domain of a neigh-boring resonance. This can be interpreted as a manifestation of “nonlinear resonance overlap” that does not requirea direct overlap of neighboring phase space islands. The
chirp-induced (or vortex-mediated) resonance overlap allows transfer of material between neighboring resonanceseven at fairly low amplitudes of the field perturbation.
The chirp-induced resonance overlap process requires
further study. We observed that the chirping rate is reduced
and the p/n = 4/5 clump vortex begins to disintegrate
when it enters the domain of the neighboring resonance.
It would be interesting to see whether the structure evaporates completely, or whether a new p/n = 5/5 structure
forms spontaneously. Such a study may have to deal more
carefully with the boundary eﬀects discussed at the end ofSec. 5.4 and in Appendix D.2; namely, the E  line broadening and associated problems with the δ f method.

6.6 Concluding remarks on chirping andbeating
Amplitude pulsations in general and beats in particular are part of and can have an impact on various processes
in systems governed by resonant interactions and nonlinear
frequency chirping. In this treatise, we have begun to elucidate the influence that beating has on the dynamics (formation and propagation) of collective structures and on the
motion of individual particles in numerical simulations ofstrong chirping with diﬀerent strengths of field damping.
In the presence of a symmetry-breaking perturbationin the field (a seed wave) and a gradient in the EP distribution, the EP distribution spontaneously develops collective density perturbations that we refer to as a primordial(nonadiabatic) hole-clump wave pair. Sustained chirping
occurs when there is a mechanism that facilitates phaseslippage between the hole-clump wave pair and the resonant seed wave. In our simulations, phase slippage resulted
from field damping: When the field constantly loses energy, it cannot fully reverse the particle displacements, sothat some particles remain untrapped and carry with them
the above-mentioned collective density modulations.
Once phase slippage has set in, beating is both a consequence and driver of chirping, where chirping reflects
the evolution of EP phase space structures. The feedback
mechanism by which beating drives the radial advancement of hole and clump waves consists of• pulsating eﬀective island widths, and• alternating O- and X-points locations,with field damping making the processes irreversible.
A manifestation of this feedback on a longer timescale (millisecond), comprising many beats and associatedfrequency steps (“micro-chirps”), is that the nonadiabatic
onset of chirping takes the form of a hyperbolic tangent:
first accelerating, then decelerating. This behavior has, albeit in a diﬀerent setting, been predicted theoretically [8]and we find qualitative agreement with our simulations.
Although the beats become increasingly irregular as

1403087-41

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

the number of pump waves increases, each beat constitutes
a global perturbation of the field that moves particles in a
globally coherent manner (unlike a stochastic process thatscatters individual particles independently). Perhaps it is
this large degree of spatial coherence (here exaggerated by
the use of a semi-perturbative model) and small degree of
temporal coherence of the beating field that makes it aneﬀective mechanism for driving the evolution of multiple
coherent phase space structures, from the first primordial(nonadiabatic) hole-clump pair, all the way to the semiadiabatic solitary vortices of the advanced stages of oursimulations. The dynamics are self-sustained as the pumping action of the phase space structures with their diﬀerent
characteristic frequencies cause the field to beat.
It appears that the overall integrity of vortices andtheir radial propagation is maintained not only in spiteof but also owing to ongoing kicks associated with therapid amplitude pulsations and phase jumps. Together withmechanisms that damp the field fluctuations, chirping andbeating prevent the system from falling into a nonlinear
equilibrium state of adiabatic stagnation by maintaining a
nonadiabatic boundary layer around quasi-adiabatic vortex structures. The dynamics in that boundary layer erode
nearby density gradients and facilitate the radial propagation of the vortices and, with that, long-range chirping andlong-range ballistic transport. Like biological lifeforms —which tend to die or hibernate in the absence of external
stimuli and without repeated throughput of material — coherent resonant structures become dynamic when receiving
stimuli that are partially incoherent and oﬀ-resonant. And
like social animals, the resonant structures can stimulate
each other by jointly causing the field to beat.
The system we have studied here was simplified inmany respects — ignoring dynamic adaptations of the fieldmode structure, external sources, sinks and collisions —which can often be important in real MCF plasmas on thelong time scales considered (up to ∼ 10 ms). However,these simplifications have made our simulations physicallymore transparent and yielded insights that would be diﬃcult to obtain with more realistic but numerically expensive
codes. Even the simplified system poses challenges, someof which have been addressed (e.g., with our approximate
quiet start technique), while others remain to be tackled(e.g., dealing with phase space loading boundaries). On themodeling side, there seems to be potential for a more realistic description of the plasma response (MHD spectra),field damping, and dynamic changes in the mode structure.
A few ideas were discussed in this paper.
There are also open questions on the physics side,
such as the following. A prompt frequency shift was observed and characterized, but a full understanding remains
to be attained. The role of beating in maintaining a nonadiabatic boundary layer around solitary vortices was emphasized, but there may be more to be uncovered. For instance, the role of the boundary layer dynamics for the radial propagation of a vortex may deserve further attention.

It is conceivable that the rate at which gradients are eroded
depends on the rate and magnitude at which the boundary
layer around a phase space vortex pulsates due to the beating field. One may then ask whether semi-adiabatic vortices have energetic or evolutionary (dis)advantages compared to nonadiabatic spreading of phase space turbulence,elaborating on related thoughts presented in this work.
Our results are relevant for the understanding and prediction of particle transport and confinement in MCF plasmas. We have presented evidence showing that the use ofquasi-adiabatic models of ballistic “bucket” transport [6]may be justified even in parameter regimes where the formal adiabatic ordering begins to fail. This seems to suggest that the eﬀect of beating may be accounted for in a
simplified manner in the transport coeﬃcients of reduced
models.

Acknowledgments
A.B. thanks Xin Wang for sharing her Python scripts
for DMUSIC analysis. Enlightening discussions with Fulvio Zonca are gratefully acknowledged, and we also thankthe members of the PPPL EP Group as well as Herbert
L. Berk and Yasushi Todo for stimulating communications during the early stages of this work (noting that ourviews expressed here may diﬀer from theirs). The work by
R.B.W. and V.N.D. was supported by the US Departmentof Energy (DOE) under contract DE-AC02-09CH11466.
The workstation used for many of the simulations reported here was funded by QST President’s Strategic Grant(Creative Research). Some of the simulations were carriedout using the JFRS-1 supercomputer system at Computational Simulation Centre of International Fusion Energy
Research Centre (IFERC-CSC) in Rokkasho Fusion Institute of QST (Aomori, Japan) (Project ID: EPTRANS).

Appendix A. Theory of Two-Wave
Beating
Consider the integral on the right-hand side of Eq. (2a)or, equivalently, the sums in Eq. (15). Express α in termsof Φ using Eq. (14c) as ω0 αB = k Φ with k ∝ nq − m, andwrite the electrostatic potential in complex form as Φ̃ =A(t)Φ̂(ψP )eiΘ(t) with Θ(t) ≡ nζ − mϑ − ω0 t − φ(t). Ignoring
damping, the equations for the amplitude (15a) and phase(15b) can be combined as
(A1)
Ȧ + iAφ̇ =d5 ZgcC0 (k  − ω0 )Φ̂eiΘ δ f ;where C0 is a real-valued function of ψp . The Jacobianhas been absorbed in the distribution F = F0 + δ f , sod2 Zgc F = N is the number of particles ( d5 Zgc δ f =0). Taking the physical density wave signal δ f to be
the imaginary (sine) component of the complex functioniΘorb, we have δ f ≡ {δ f˜} = − 2i (δ f˜ − δ f˜∗ ), whereδ f˜ = 2CC0 e
the asterisk indicates a complex conjugate. The form of
the phase Θorb of the density wave in GC phase space will

1403087-42

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

be specified later. Substitution into Eq. (A1) yields
Ȧ + iAφ̇ = − d5 Zgc (ik  − iω0 − iφ̇)C Φ̂ei(Θ+Θorb )+ d5 Zgc (ik  − iω0 − iφ̇)C Φ̂ei(Θ−Θorb )− iφ̇ d5 ZgcC Φ̂ei(Θ+Θorb )+ iφ̇ d5 ZgcC Φ̂ei(Θ−Θorb ) .
(A2)
The integrals are nonzero only when the complex exponential is a constant; namely, for phase-matched waves. Let us
assume that this is the case for the diﬀerence between Θand Θorb , so that Θ − Θorb = β0 + 2πl with β0 = const. andinteger l. The factor (k  − ω0 − φ̇) is readily recognized
as the time derivative of the phase diﬀerence between thefield and density waves, so that we can writedd5 ZgcC Φ̂ei(Θ−Θorb )Ȧ + iAφ̇ =dt
+ iφ̇ d5 ZgcC Φ̂ei(Θ−Θorb ) .
(A3)
Evidently, the integral expression corresponds to A(t).
For simplicity, we ignore all the radial dependencies,
so that Φ̂ = 1. Moreover, the distribution is assumed tocontain particles with only one value of the magnetic moment μ and to be monoenergetic in the rotating frame withenergy E  = E − ωn Pζ . We can then integrate the distribution of phase space density waves δ f over Pζ along theline E  = const. In this way, the time-dependencies of allpump waves are merged into a single signal and we obtaina combined complex pump wave functiondPζ dKdμ δ f˜(Pζ , μ, K, ϑ, ζ, t);δg̃(t)einζ−ipϑ ≡(A4)
where p is the number of elliptic points of the resonancealong ϑ.
To be more concrete, suppose that the perturbed
component of the distribution function δ f consists of
two sinusoidal pump waves located at diﬀerent radiiPζ : a clump wave with frequency ω+ and a hole
wave with frequency ω− , which can be written as± ±
A0 w± 1±
2 1 + sin(nζ − pϑ − ω t) ≥ 0 with weights w =±1. The parameter −1 ≤  ≤ 1 measures the diﬀerence inthe pump wave amplitudes. Here, the clump will dominate
for  > 0 and the hole will dominate for  < 0. The resulting combined complex pump wave function δg̃ is then


−
A0
+
−iω+ t+ w− (1 − )e−iω tδg̃(t) = 2(1+2 ) w (1 + )e
A0 −iω̄t Δωe
 cos( Δω=2 t) − i sin( 2 t) , (A5)21+

base waebeat=b(t) exp(−iβ(t))

with mean frequency ω̄ = 12 (ω+ + ω− ), and frequency difference Δω = ω+ − ω− . The amplitude is normalized toFig. A1 Example for asymmetric beating of two pump waves;namely a hole and a clump separated by Δν = 10 kHz
and with relative amplitudes diﬀering by 2 = 0.4.
The plot shows the time traces of b(t) and β(t) from
Eq. (A6).

A0 and constant oﬀsets are ignored. The complex exponential factor in Eq. (A5) represents an eﬀective base wave
oscillating at the mean frequency ω̄, which is not present
in the original signal30 . The beat consists of a modulation
amplitude b(t) and phase β(t), given by1/2
,b(t) = sin2 (Δωt/2) +  2 cos2 (Δωt/2)
β(t) = arctan  −1 tan(Δωt/2) ,

(A6a)
(A6b)

in the case of a superposition of a hole-clump pair. The
frequency of the beats is Δω = ω+ − ω− (without the factor
1/2), so the pulse length is τpulse = 2π/Δω. An example
with  = 0.2 and Δω = 2π × 10 kHz is shown in Fig. A1.
Times t0 satisfying Δωt0 /2 = πl with integer l correspond to instants of maximal destructive interference,
where the modulation factor b(t) has a minimum and therate of phase shift |β̇(t)| is maximal. Maximal constructive
interference occurs at times t1 satisfying Δωt1 /2 = π/4+πl.
To avoid confusion, one should keep in mind that destructive interference between a hole wave and a clump wave
occurs when their phases are aligned at t = t0 , because theirweights w± have opposite signs. If we superimpose twoclump waves, the functions i sin(Δωt/2) and cos(Δωt/2)
would switch places in Eq. (A5), and times t = t0 would be
the instants of constructive interference.
Equation (A6b) shows that the total phase shift between successive pulses of two-wave beating is always ±π
and that the direction of a phase shift is determined by thedominant pump wave at that time; i.e., the sign of . In our
setup, clump and hole waves yield positive and negative
phase shifts, respectively. In the symmetric case  = 0,
with two pump waves of identical amplitude, the phaseshift becomes instantaneous and has no specific direction(degenerate sign); the sign of the modulation factor simply
flips abruptly when b = 0 and remains constant at othertimes. In the present paper, the term “phase jump” refers
collectively to both smooth phase shifts and sign flips.
30 In signal processing, Eq. (A5) is known as “double-sideband modulation” and the eﬀective base wave is called “suppressed carrier”.

1403087-43

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Finally, we determine the response of the combined
field amplitude that, according to Eq. (A3), obeys
A(t) =  d5 ZgcC Φ̂ei(Θ−Θorb ) .
(A7)
Based on Eq. (A5), the phase of the combined pump waves
with resonance number p can be written
Θorb = nζ − pϑ − ω̄t − β(t).

(A8)

Note that the integration along GC drift orbits turns
the single poloidal harmonic e−imϑ into a broader spectrum morb e−imorb ϑ , which includes the resonance morb =p [31].
Letting again Φ̂ = 1 and substitutingdPζ dKdμ Ce−iΘorb = δg̃∗ eipϑ−inζ , we obtainA=dϑdζ δg̃∗ ei[(p−morb )ϑ−ω0 t−φ(t)] ,

(A9)

morb

=

A0 b(t)  i(ω̄−ω0 )t i(β(t)−φ(t))
 ee.
1 + 2(A10)

When the hole and clump wave frequencies are situatedsymmetrically around the seed frequency, we have ω̄ =
ω0 . Assuming phase matching, φ(t) = β(t) + β0 + 2πl, andignoring β0 ,31 we obtain the result
A0b(t),1 + 2φ(t) = β(t) + 2πl.
A(t) =(A11a)(A11b)

This shows that the field wave beats in precisely the same
way as the combined pump wave function δg. Since the latter was a superposition of multiple pump waves, we mayview the combined field wave also as a linear superposition of multiple harmonic waves with amplitudes Ak (t) andphases φk (t), each driven by a diﬀerent pump wave:A(t)eiΘ(t) =Ak (t)einζ−imϑ−iω0 t−iφk (t) .
(A12)k

The response of the field to the concert of pumpwaves, as in the process of two-wave beating analyzed
above, may be viewed as an auto-resonance phenomenon;
albeit with one additional twist: Usually, the term “autoresonance” is used to refer to the automatic phase lockingof particle bunches or waves in response to (chirped) external drive, as in a synchrotron [56, 57]. However, in thesystem we are considering, the prefix “auto” can also be interpreted as “self” because the combined field-particle system contains feedback loops that make it self-resonant andself-propelled in the presence of a source of free energy inthe form of destabilizing gradients. In other words, we are
dealing with auto-resonance in an active medium. When
this field-particle system is “opened”, in our case through31 The value of β

0 determines the phase lag between the field and themean density wave: ei(ω̄−ω0 )t+iβ0 . Presumably, one would have to invoke

the Vlasov equation in order to determine the proper value of β0 , but it isirrelevant for the amplitude modulations and phase jumps.

the introduction of field damping, a directional flow of energy is established, and nonlinear frequency chirping is one
manifestation of that.
We would like to point out that there exists an alternative, albeit purely phenomenological, interpretation of thebeats: Loosely speaking, one may also interpret each pulseof A(t) as a new instability whose phase is shifted by ±πwith respect to its predecessor and successor. This picture
exemplifies the connection between the beats in our simulations and the so-called “relay runner” model of nonadiabatic chirping proposed by Zonca & Chen [8], which we
discussed in Secs. 5.2 and 6.2. However, the purely phenomenological nature of this picture in the present context is apparent from the time trace of A(t) in Fig. 7, whichdrops sharply before and recovers equally sharply after anear-zero minimum. This corresponds to a hyperbolic time
trace of the growth rate γ(t) (cf. Fig. D7), and we see noway to reconcile this with any known form of resonant
drive or damping, where dA/dt should decrease (not increase) when the value of A goes to zero, since the amountof resonantly transferred energy per unit time should decrease with decreasing amplitude A as the interaction becomes weaker. At present, the only physically viable interpretation for the observed signals that we can oﬀer isthe beating picture, where the pulses and phase jumps are
the consequence of a superposition of multiple coexisting
pump waves whose amplitudes do not vanish even when
the overall field amplitude does so temporarily during instants of destructive interference. Supported by the above
derivations for the case of two-wave interference, the physical picture of beating is the foundation for this work.
An important implication of the beating picture is thatthe instantaneous growth rate γ(t) of the beating signal
does not represent resonance width. While the growth or
damping rate |γ(t)| associated with the combined amplitude A(t) can become very large (e.g., see Fig. D1 (c)), theamplitudes of individual components Ak (t) in Eq. (A12)usually vary relatively slowly. Thus, instead of one γbroadened resonance, we have multiple resonances whose
widths correspond to the size of a quasi-adiabatic islandcore (if any) plus the width of a nonadiabatic boundary
layer, as illustrated in Figs. 1 (b,c). The beats (or any
pulsations) of the field contribute to this resonance width
through their influence on the boundary layer, where particles are detrapped and retrapped repeatedly, as shown for
the undamped case in Appendix B.2 and in the movie inFig. E3 (d) below. Clearly, this resonance width depends
on the magnitude of the pulsations.

Appendix B. Pulsations of Undamped
Resonances
Figure B1 shows an overview of the field and phasespace dynamics in the absence of damping. The weakly
driven case (A0) on the left and the strongly driven case(B0) on the right are qualitatively similar, so snapshots of

1403087-44

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. B1 Overview of the field and phase space dynamics during and after saturation of a weakly driven [case (A0), left] and a strongly
driven resonant instability [case (B0), right] in the absence of damping. These simulations were performed with shorter time stepΔtˆ = 10−3 (default: 1/400) to reduce spurious growth. Panels (a) and (c) show the respective time traces of the field amplitude
A(t), and panels (b) and (d) show the corresponding high-resolution DMUSIC spectrograms (twin = 0.047 ms). For each case,several snapshots of the phase space density perturbation δ f are shown as colored contour plots. For case (B0), the three arrows
labeled “1”, “2” and “3” roughly indicate the observed flow of EP Vlasov fluid in the boundary layer during about half of abounce period. A link to the movie for case (B0) is provided in Fig. E3 (d).

δ f are shown mostly for case (B), where the fluctuations
are larger, both in magnitude and spatial extent, so they arealso better resolved with our diagnostic mesh. Figure B2
shows the motion of a few tracer particles in the domain of
these resonant structures.
In Sec. 2.4, we have already given a brief review ofthe saturation process of such undamped resonant instabilities and demonstrated that the subsequent pulsations of
the mode amplitude correlate well with the bouncing motion of particles trapped in the eﬀective potential well of
the field wave. Here, we analyze these dynamics in somemore detail and highlight two interesting features that are
related to the topic of the present work; namely:(i) We find evidence for variations in the field’s oscillation frequency and even transient frequency splitting.
(ii) The continuing pulsations of the field amplitude
A(t) maintain a nonadiabatic boundary layer, insidewhich marginally resonant particles are repeatedly
detrapped and retrapped.

B.1

Spectral fluctuations

Transient frequency splitting can be seen in the highresolution DMUSIC spectrogram in Fig. B1 (d), aroundt = 0.3 . . . 0.4 ms, just before the field amplitude reaches its
peak value. During that short interval, the initial 100 kHz
signal seems to split temporarily into ν− ≈ 98 kHz andν+ ≈ 103 kHz. The magnitude of the observed frequency
shifts by Δν± ≈ ±2 . . . 3 kHz relative to the seed wave frequency ν0 = 100 kHz, is consistent with the poloidal groupvelocity ±ϑ of the primordial hole and clump wave fronts

that can be seen to pass above and below the resonance’s
O-point around the time of the first δ f snapshot for case(B) in Fig. B1, at tˆ = 70.5 (0.33 ms): our measurement
gives
±ϑ
∓0.045πΔϑ±≈= ∓2.5 kHz; (B1)=λϑ λϑ Δt 0.365π × 0.05 ms
where λϑ ≈ 0.365π is the poloidal island length in the region around ϑ = 0 (cf. Fig. 13).
Thus, one may speculate that the pumping action ofthese phase space structures is the cause of the observed
transient frequency splitting. The subsequent persisting
frequency fluctuations seen in the two spectrograms in panels (b) and (d) of Fig. B1 may also be associated with the
dynamics of density waves, especially those in the island’s
nonadiabatic boundary layer.
Although the spectral fluctuations appear to be consistent with the dynamics of phase space structures, the spectral data have to be interpreted with care. We were not able
to confirm the frequency splitting phenomenon seen in the
DMUSIC analysis of case (B0) using conventional Fourier
analysis. At present, we cannot rule out the possibility thatfrequency splitting as seen in Fig. B1 (d) has origins otherthan the motion of the primordial hole and clump waves.

B.2 Nonadiabatic boundary layer
Consider the second and third snapshots of δ f for case(B0) in Fig. B1 at tˆ = 81.5 and 98.5, where we have drawn
a set of three short arrows that roughly indicate the flowof EP Vlasov fluid during about half of a bounce period.

1403087-45

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. B2 Motion of well-trapped (“center”, “bulk”) and temporarily trapped (“boundary”) particles in the undamped cases (A0) (left)and (B0) (right). Panels (a) and (f) show the time traces of the radial position P̂ζ (t). Panels (b) and (g) show the locations of
the particles inside the δ f phase space structures during an advanced stage of the simulations. These δ f -weighted Poincaré
plots were accumulated during one transit period τζ0 = 4.7 μs, during which each tracer particle appears in the Poincaré section
approximately 4 times. Panels (c)-(e) and (h)-(j) show Fourier spectra of the P̂ζ (t) signals obtained with sliding time windows ofsize Δtwin = 1.41 ms, where one can infer the bounce frequencies νb of resonantly trapped particles.

The first arrow, labeled “1”, passes near an eﬀective Xpoint at (P̂ζ , ϑ) ≈ (0.719, −0.1π) while the field amplitudeis large. The pair of arrows labeled “2” and “3” approaches
the same X-point during an amplitude minimum. The innerarrow “2” circulates back around the eﬀective O-point of
the resonant phase space structure, while the outer arrow“3” passes above the eﬀective X-point and continues its
poloidal drift towards the left. This process repeats during
each pulse and can also be seen in the motion of individual
particles, which we discuss in the following.
For both cases (A0) and (B0), we have chosen threeparticles, whose trajectories P̂ζ (t) are plotted in panels (a)and (f) of Fig. B2. Fourier spectrograms of these time
traces are plotted in panels (c)-(e) and (h)-(j), which effectively show the bounce frequencies (if any). In panels(b) and (g), we show snapshots of the phase space structures in δ f after several bounces. The circles and crosses
indicate the positions where the three tracer particles haveappeared in that Poincaré section (which rotates toroidallywith the seed wave) during one transit time τζ0 ≈ 4.7 μs.
• Tracer #1 labeled “center” circulates not far from theO-point. It has the largest bounce frequency, which isparticularly clear in panel (j) for case (B0).
• Tracer #2 labeled “bulk” is also well-trapped, but circulates closer to the boundary, with a smaller bounce
frequency than the “center” particle.
• Tracer #3 labeled “boundary” lives in the boundary
layer and undergoes repeated detrapping and retrapping. Depending on whether it approaches an X-point
during an amplitude maximum or minimum, it is reflected or passes. The Fourier spectra in this case areunintelligible, but the aﬀective bounce time can be inferred from the raw P̂ζ (t) data in (a) and (f).

Fig. B3 Time trace of the rotating frame energy E  (t) for thesame three tracers as in Fig. B2 (a) for case (A0).

The behavior of the “boundary” particle was illustrated
schematically in Fig. 1 (b). The nonadiabatic boundary
layer is maintained by the amplitude pulsations seem inpanels (a) and (c) of Fig. B1. Judging by the size of thestructures seen around the X-points in panels (b) and (g)of Fig. B2, the radial width of this layer seems to be about
1/5. . . 1/3 of the width of the eﬀective phase space island.
Dewar’s theoretical prediction for the bounce frequency of deeply trapped particles [7], νbDewar = 2.88 ×
γL /(2π) in Eq. (6), gives 1.44 kHz for case (A0) and10.1 kHz for case (B0). Given the approximations madein the theory, we consider these estimates to be in good
agreement with the simulation results in Figs. B1 and B2.
Finally, Fig. B3 demonstrates the high degree of energy conservation in the undamped simulation: |δE  |/E  ∼10−5 , which is two orders of magnitude better than in thesimulations of damped cases (cf. Fig. 29).

1403087-46

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Appendix C. Other Interesting Observations
C.1

Splitting of spectral lines (in DMUSIC)

Panels (a) and (b) of Fig. 21 both show a splitting ofsome spectral lines, especially those associated with the
upward chirping solitary clump vortices. Another view isgiven in the second zoom-up of Fig. E1 (b).
So far, we have seen this phenomenon only in DMUSIC spectrograms exploiting the complex frequency plane
as described in Appendix E.2, not in FFT spectrograms orwhen DMUSIC is applied only along the real frequency
axis. Since the split strands here are separated by only
about 2 kHz, the splitting is seen only with time windows
Δtwin that have width of about half a millisecond or longer.
More strands tend to appear with larger time windows.
The line splitting seems to become stronger with increasing chirping rate, so it is most easily seen along the spectral lines of massive clumps when they chirp rapidly upward, while traveling radially outward, towards the peakof the mode. But with suﬃciently large time windows,
line splitting can also observed on downward chirps, especially around times where the chirping rate changes relatively abruptly. As the chirping rate or intensity of thesignal varies, the line splitting pattern may change. For instance, it may alternate between odd and even numbers of
strands, or exhibit an oscillatory pattern.
It remains to be clarified what factors (physical, numerical, post-processing) cause these spectral patterns. For
instance, one may check whether there is a relation to Fresnel ripples known to aﬀect chirp spectra32 . The patterns
vary depending on DMUSIC parameters, such as the time
window size and the number of damped sinusoids. However, that does not necessarily mean that it is entirely apost-processing artifact. The spectrum in complex frequency space may truly (for physical reasons) be so complicated that its appearance changes drastically depending
on how it is being processed and visualized.
If this kind of splitting has physical reasons, a better understanding of such spectral patterns may enable us
to extract information about processes such as the accumulation and detachment of massive hole and clumps, orthe internal dynamics of structures such as the solitary
clump vortex that resembles a spiral galaxy in Fig. 3 (b).
At present, however, this is mere speculation.

C.2

Anomalies in bounce frequency spectra

The spectrograms of the bouncing motion of resonantly trapped particles that were shown in Fig. 28 exhibit
rich dynamics, and it may be interesting to analyze the underlying reasons. Here, we merely speculate.
In the stationary case, with fixed amplitude A andphase φ, it is known that the bounce frequency decreases
monotonically from the O-point to the boundary (separa32 When we apply a Hann window before the DMUSIC algorithm, thepattern becomes smoother, but the line splitting as such remains.

trix) of a resonant phase space structure. This has recently
been reconfirmed using ORBIT [5]. However, our simulations show interesting anomalies that violate this rule inthe case of chirping solitary vortices. On the left-hand sideof Fig. 28 for case (A), one can see that, during the period
4 ms  t  6 ms, tracer #3 near the boundary bounces approximately as rapidly as tracer #1 in the core of the clump
with νb ≈ 4 . . . 5 kHz; i.e., faster than tracer #2 in the bulk
with νb ≈ 3 . . . 4 kHz. Similarly, on the right-hand side ofFig. 28 for case (B), one can see that, during and around theinterval t = 8 . . . 10 ms, the tracer #2 in the bulk bounces
with a significantly higher frequency νb ≈ 5 kHz than thecentralized tracer #1, which has νb ≈ 3 kHz.
The reason for this behavior is unclear. We suspecta connection to the processes that make these vortices
grow, detach from the turbulent belt and propagate radially (chirp) as discussed in Sec. 5.3. This speculation ismotivated by the fact that the bounce frequencies vary significantly during the course of a simulation. The largest
bounce frequencies can be observed in Fig. 28 during thedetachment process and when the chirping rate is large.
The overall trend is that the bounce frequencies decrease in time. This observation is in itself surprising,
because these clumps propagate radially outward towards
the peak of the mode, so they are actually experiencing
stronger field oscillations of progressively larger amplitude when A(t) fluctuates more or less at the same level.
This contradicts the expectation that ω2b ∝ A(t) × ξ̂(ψP ).
However, this prediction is based on linear drive, ωb ∝ γL
(cf. Sec. 2.6), so it should be expected to hold only during the first bounces after the resonant instability saturates.
Further study is necessary for clarification. One may also
investigate the cause of the relatively rapid oscillation of
the bounce frequencies νb (t) in Fig. 28, which have a timescale comparable to the bounce period itself.

C.3 EP distribution
chirpsduringlong-range

In Sec. 5.4, we have shown that particles in our chirping simulations depart from the line E  = 75.6 keV inthe (Pζ , E)-plane where they have initially been loaded(Fig. 29). To motivate further study, we include here amore detailed view of the distribution’s 4-D structure in(δ f, P̂ζ , E  , ϑ) and (P̂ζ0 , P̂ζ , E  , ϑ) space, four snapshots ofwhich are shown in Figs. C1 and C2 for cases (A) and (B),respectively. The left half of these figures shows the domains occupied by structures with diﬀerent values of δ f ; inparticular, holes (δ f > 0, red) and clumps (δ f < 0, blue).
The right half of these figures shows the radial transportand mixing of phase space.
The crescent-shaped distortion that is centered aroundthe seed resonance P̂ζ,res ≈ 0.72 seems to be physical, since
it is also present in simulations with particles loaded in aband of width ΔE0 = 7.5 keV (see Appendix D.2 below).
We suspect that this crescent is due to chirping, but this1403087-47

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. C1 Particle transport in (P̂ζ , E  , ϑ)-space in the marginally unstable case (A). This figure provides a more detailed view of theparticle distribution data in Fig. 29 (a) at four snapshot times: t = 2.1, 4.0, 4.5 and 5.5 ms (top to bottom). The colors in columns(a) and (b) on the left-hand side represent marker weights w j = δ f j /G, so red means that a marker represents an increased phasespace density (δ f > 0) and blue means reduced density (δ f < 0) at that instant of time. The colors in columns (c) and (d) on theright-hand side represent the initial radial position P̂ζ (t = 0). The dashed lines in columns (b) and (c) represent the inclinationof E  − (ω − ω0 )P̂ζ /n = const. lines for frequencies ω = 2πν at the edges of the crescent-shaped distortion of the distribution
or the tip of the outermost massive clump structure. The values of ν associated with these dashed lines are shown in column(c). Note that column (b) has a fixed scale in both E  and P̂ζ , so that one can clearly see the broadening of the distribution fromsnapshot (b1) to (b4). In columns (a), (c) and (d), the limits of the E  -axis vary between snapshots, so that one can clearly see
the structures.

hypothesis remains to be confirmed. A physics study ofthis and other features seen in Figs. C1 and C2 requires
more careful considerations of boundary eﬀects associatedwith the δ f method, as discussed at the end of Sec. 5.4.

C.4

Ghost chirps and resonance overlap

The spectrogram in Fig. 2 (b) contains a downward
chirp labeled “ghost”, which diﬀers from the other downchirps in that it does not seem to have a corresponding inward propagating hole structure in EP phase space. Further examples can be seen in Fig. D12. These ghost chirps
seem to be somewhat correlated with the long-range upward chirps. According to our observations, a downward
chirping ghost appears whenever an upward chirp exceeds
150 kHz. This coincides with the time after which we observe enhanced E  line broadening in Figs. 29, C1 and C2
discussed in the previous Sec. C.3. Indeed, we do not see
any ghost chirps in a simulation with ΔE0 = 7.5 keV that
we performed for validation purposes (see Fig. D1 (d) be-low). Therefore, we suspect the ghosts to be a consequenceof boundary eﬀects caused by letting ΔE0 = 0.
Although the cause of the ghosts seems to be an artifact of sampling only a part of the GC phase space,
the mechanism of their generation in that particular form
may be physical. The first snapshot of Fig. 11 shows thata neighboring resonance with p/n = 5/5 elliptic points
dominates the structure of the δ f landscape in the region P̂ζ  0.62, which corresponds to ν  150 kHz.
When the upward-chirping solitary clump vortices with
p/n = 4/5 enter that domain, we expect that these structures are modulated by the p/n = 5/5 resonance, and that
this modulation produces an image of the original oscillation at a diﬀerent frequency. Poloidal transit frequency estimates suggest that the resonances are separated by aboutωϑ /(2π) ≈ 210 kHz. This is consistent with the appearance of ghosts around 60 kHz when a clump hits the markof 150 kHz. Thus, the ghosts may also be an indication of
the chirp-mediated nonlinear resonance overlap discussed
in Sec. 6.5. In fact, even the enhanced spreading of parti-

1403087-48

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. C2 Particle transport in (P̂ζ , E  , ϑ)-space in the strongly unstable case (B). Arranged in the same way as Fig. C1, this figure provides
a more detailed view of the particle distribution data in Fig. 29 (b) at four snapshot times: t = 0.7, 2.0, 4.2 and 14 ms (top tobottom).

cles away from the initial line E0 = 75.6 keV in Figs. 29,
C1 and C2 may at least partially be a consequence of thiskind of resonance overlap.
Note that ghost chirps have also been reported in otherstudies. For instance, see the recent analysis of ghost features in ion cyclotron emission (ICE) spectra [58].

diﬀerent poloidal harmonics m for a given toroidal harmonic n. This procedure may allow to simulate resonant
instabilities with relatively short wavelengths (which usually interact best with EPs that have relatively small drift
orbits). Nevertheless, as discussed in Sec. 6.2, even ourpresent simulations with a single harmonic (m, n) = (6, 5)
can realize a form of convective amplification:• Firstly, although the peak of the mode does not shift,
the phase space structures are convected radially witheach pulse of the beating field. Successive pulses canbe amplified, especially in cases close to marginal stability (Fig. 20, left).
• Secondly, by placing the seed resonance oﬀ-peak (cf.
Fig. 9), some phase space structures (in our case theclumps) travel into regions with stronger fluctuations
even if A(t) remains around the same level.

Appendix D. Characterization andVerification of the Model
D.1

Background plasma response

The semi-perturbative approach used in this work,where we prescribe the radial profile of the field mode,
may be regarded as an application of the Rayleigh-Ritz
method: We are probing the system’s frequency responsewith an approximate (eigen)mode, thus, reducing the infinite degrees of freedom of the real system to a few variables; here, the field amplitude A(t), the phase φ(t), and abunch of marker particles that sample the EP guiding center phase space.
The process of convective amplification [4] in theform illustrated in Fig. 1 (d) could also be simulated in ourmodel using Eq. (13) with the mode index k identifying independent pairs of Fourier harmonics (m, n); i.e., ignoring
the geometric (toroidal, elliptic, etc.) coupling between

In some sense, the use of an oﬀ-peak resonance can be said
to exchange the cause and eﬀect of convective amplification: instead of a new peak growing at a new location after
the resonant drive has propagated radially (post-convective
amplification), here the resonant phase space structures experience a larger field amplitude on one side than on theother and their resonant drive eﬀectively pre-amplifies thefield further in the direction of the peak (pre-convective
amplification). This may accelerate the growth and radial

1403087-49

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

drift of the phase space structures that are predisposed topropagate towards the peak; in our case the clumps, because our peak is located in a region of lower EP density(at a larger radius) than the resonance. We suspect thatthe process of pre-convective amplification is not merely
a curiosity of our simulation setup, but may bear practicalrelevance, because oﬀ-peak resonances are common andsince Alfvén mode profiles seem to have a certain degreeof rigidity, at least in the case of long wavelength modes.
Motivated by the observation of global beating in hybrid simulations [12] as discussed in Sec. 2.7, we assume
that, for our current purposes, it is acceptable to ignore
perturbations of the amplitude profile for long-wavelength
modes. However, it is not evident to what extent our ignoring of the dynamic adjustments of the phase profile isjustified, since this directly touches the quantity φ(t) responsible for chirping. This may require further consideration, which we choose to postpone. Discussions concerning the role of the phase profile during interaction between
EPs and shear Alfvén waves can be found, for instance, inSec. III D of Ref. [42] and in Refs. [59, 60].
While perturbations of the field mode structure are ignored, the frequency is entirely determined by the resonant
drive in our simulations. We even observe a prompt frequency shift that we attribute to an auto-optimization process, which will be discussed in Sec. D.4 below. Although
it is possible to include the MHD plasma response without
directly solving MHD equations, the required extensions
are not trivial. Wang Ge et al. [22] have developed andsuccessfully applied a reduced model for Alfvénic chirping that includes the eﬀect of continuous spectra. We havemade no attempt to adopt this model yet, because it requires sophisticated numerical methods to deal with thesmall scale structures that inevitably develop in an ideal
MHD formulation. In the end, however, these troublesome
micro-structures seem to have little eﬀect on the chirping
dynamics when the EP’s magnetic drift is significant. Perhaps, it is worthwhile to develop a reduced kinetic description of the continua, where the ion Larmor radius would
impose a lower bound on the spatial scales. This may reduce the numerical obstacles and give an even more realistic representation of the continua in the form of kinetic
Alfvén waves, whose wave vector contains also a radial
component (radiative damping) [61, 62]. Another possibility one may consider is to use a precomputed damping rate
γd (r, ω) that depends on radius and frequency, for instanceusing a gyrokinetic solver for the plasma dispersion relation, such as LIGKA [63, 64]. Moreover, it has been shownthat radial electric fields can also have an impact on the
structure of Alfvén continua in the frequency range of interest [65] and, thus, may aﬀect chirping dynamics.

D.2

Case (A’): Nonzero loading width ΔE0

The condition dG/dt = 0 underlying our implementation of the δ f method, namely Eq. (18), does not strictlyFig. D1 Evolution of case (A’), where particles are loaded in
an energy band of width ΔE0 = 7.5 keV around E0 =75.6 keV as shown in Fig. D2. All other parameters are
the same as in the marginally unstable case (A), where
ΔE0 = 0. Panels (a) and (c) show, respectively, thetime traces of the amplitude A(t) and growth rate γ(t)plus damping rate γd /ω0 = 7.5%. The data of case(A) is also plotted for comparison. Panel (b) shows along-time Fourier spectrogram obtained with large sliding time window Δtwin = 0.47 ms, and the dotted green
curves represent the prediction of the BB model (12).
Panel (d) shows a snapshot of δ f (P̂ζ , ϑ) taken at timet = 7 ms, using the modified radial coordinate P̂ζ asshown in Fig. 10. The full 3-D distribution δ f (P̂ζ , E  , ϑ)is shown in Fig. D3.

hold in our simulations, since particles are found to travel
outside the region that was initially filled with phase spacemarkers. This was shown in Figs. 29, C1 and C2. Nevertheless, our choice to load marker particles only along theline E  = E0 = 75.588 kHz = const. in the present study
is justified because the errors remain small and because weare able to reproduce the essential features of earlier simulations, where markers were loaded in an energy band ofnonzero width ΔE0 [13, 37].
This is demonstrated here using the example inFig. D1, which was obtained with the same parameters as
our marginally unstable case (A), except that the initial
distribution of marker particles covered a band of width
ΔE0 = 7.5 keV around the E  = 75.588 kHz line. This
modified case is called (A’). Figure D1 shows that theearly evolution of the field is nearly identical in cases (A)and (A’). Diﬀerences in the long-term evolution can be
attributed to the high sensitivity of these dynamics withrespect to physical and numerical parameters (e.g., see
Figs. D10 and D12 below).

1403087-50

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. D2 Comparison between the marker particle distributions
in our default case (A), where ΔE0 = 0, and the modified case (A’), where ΔE0 = 7.5 kHz. The number ofparticles is Np = 106 in both cases.

fact of the δ f method as discussed at the end of Sec. 5.4.
Finally, returning once more to the overview in
Fig. D1, note
√ how the up-chirping clump wave front isclose to the t curve (green dots) of the BB model (12) for
several milliseconds. We would like to emphasize that one
should not interpret too much into this “agreement”, which
seems to be largely accidental, because we have repeatedly
observed that a small change in the parameters can significantly alter the chirping pattern; for instance, giving rise
to rapid chirps as in Fig. 2 (another example will appear in
Ref. [66]). Disagreements between the BB model and oursimulation results (here the downward chirps, but generally also the upward chirps) may be easily explained: toostrong drive, too rapid chirping, nonuniform mode structure and other asymmetries, plasma geometry, large magnetic drifts. The point is that an (apparent) agreement in asingle case can be misleading. Sensitivity checks are important and constitute a large portion of this Appendix.

D.3

Case (A”): Relaxed initial EP distribution

Our standard initialization procedure was to use anunperturbed axisymmetric particle distribution and a fieldthat contains a non-axisymmmetric perturbation in theform of a traveling wave with initial frequency ν0 =100 kHz, poloidal mode number m = 6 and toroidal mode
number n = 5. This means that, when a simulation isstarted with a destabilizing gradient, two things happen simultaneously:• The perturbation amplitude grows and, with that, thesize of the resonantly perturbed phase space domain.
• The particle distribution adapts to the growing nonaxisymmetric perturbation.

Fig. D3 Snapshot of the 3-D particle distribution δ f (P̂ζ , E  , ϑ)at t = 7 ms in case (A’), where ΔE0 = 7.5 kHz.

Figure D2 shows that the broadened distribution ofcase (A) at t = 10 ms (light blue) lies well within the particle loading range of case (A’) (gray). The 3-D view ofδ f (P̂ζ , E  , ϑ) in Fig. D3 shows the cylinder-like structureof the clumps and (albeit partially hidden) the holes. Taking advantage of this quasi-2-D structure, the simulations
analyzed in this paper were performed for only a vanishingly thin slice of this distribution (cf. Fig. 10).
Figure D2 also shows that the small crescent-shaped
distortion of the distribution that we observed in case (A)extends all the way to the edges of the 7.5 keV wide ΔE0band in case (A’). This supports our assertion made inSec. 5.4, that this crescent-shaped distortion breaking theE  = E − ωPζ /n = const. condition may be a consequenceof chirping; i.e., variations in the frequency ω. Meanwhile,enhanced scattering of particles seen in Figs. 29, C1 andC2, especially in the domain of the clumps, is not visible
in Figs. D2 and D3. Such a diﬀuse spreading may be occurring inside the ΔE0 = 7.5 keV window, but we have notchecked this yet. The absence of such diﬀuse spreading incase (A’) would suggest that it is a (small) boundary arti-

We were curious to see what happens if we let the particle distribution adapt before we let the instability grow,so we have launched a simulation from the state shownin Fig. 433 . That state was produced by fixing the field
amplitude at a small value of A = 10−6 and letting the
EP distribution adapt to the perturbed field for 500 transits(2.35 ms). The subsequent simulation with evolving amplitude A(t) and phase φ(t) used the same parameters as
the marginally unstable case (A), and we call this modified
case (A”).
The results are shown in Fig. D4. In panel (a), one cansee that the instability in case (A”) is significantly enhanced
compared to the default case (A). Exponential growth is
not observed, as can be confirmed in panel (c). Instead, theamplitude grows almost in a straight line. The first pulse,which saturated at a tiny amplitude of A ≈ 6 × 10−6 in case(A), is eﬀectively skipped. Instead, case (A”) saturates forthe first time at the level A ≈ 3 × 10−5 , similar to the peakvalue 3.5 × 10−5 of the second pulse in case (A). The high33 Here, we have strong symmetry breaking in the relaxed initial EP
distribution. A related problem, albeit with an unmodulated phase spacedensity plateau perturbed only by noise, has been considered in Ref. [67].

1403087-51

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. D4 Evolution of case (A”), where the amplitude and phasevelocity of field were fixed for 500 transits (2.35 ms)
before launching the simulations with the same parameters as the marginally unstable case (A). The EP distribution at the beginning of the simulation was shownin Fig. 4. Panels (a) and (c) show, respectively, the time
traces of the amplitude A(t) and growth rate γ(t) plus
damping rate γd /ω0 = 7.5%. The data of case (A) arealso plotted for comparison. Panel (b) shows a longtime Fourier spectrogram obtained with sliding time
window Δtwin = 0.47 ms. For the first 0.5 ms, panel(d) shows a high-resolution DMUSIC spectrogram obtained with Δtwin = 0.047 ms.

resolution DMUSIC spectrogram in panel (d) of Fig. D4shows eﬀectively immediate frequency splitting, as may
be expected from the preestablished hole-clump pair seenin Fig. 4. The long-time Fourier spectrogram in panel (d)of Fig. D4 shows long-range chirping dynamics similar to
those observed in case (A).
In summary, the main consequence of initializing asimulation of nonlinear frequency chirping with a relaxed
non-axisymmetric EP distribution is that one eﬀectively
skips the exponential growth phase and the first pulse.
Since the primordial (nonadiabatic) hole-clump wave pair
has already been established in case (A”), the simulation
goes straight into beating and chirping. The two initialization methods used in cases (A) and (A”) — axisymmetricvs. relaxed EP distribution — are perhaps both unrealistic, and the real situation may lie somewhere in between.
The choice of initialization method does not seem to play asignificant role for most of the dynamics examined in thiswork, except for the processes of early structure formationand first saturation described in Sec. 4.

D.4

Auto-optimization of resonant drive
At first glance, Eqs. (15) and (18) suggest that theFig. D5 Dependence of the linear kinetic drive γL on (a,b) theradial density gradient F0 ≡ ∂Pζ F0 and (c) the damping rate γd . The results in panels (a) and (b) are obtained without damping (γd = 0) and the axes are normalized by the values of γL and F0 in cases (A0) and(B0), respectively. The dashed lines indicate perfectlinear scaling γL ∝ F0 . Note that the γd scan in (c)is performed with fixed steep gradients F0 , so the data
points for γd = 0 do not correspond to cases (A0) and(B0), which have shallower gradients. Panel (d) shows
how the destabilizing gradient F0 has to be increased
exponentially in order to maintain the same linear drive
γL /ω0 = 8% when the mode’s peak is shifted radiallyoutward, away from the resonance (cf. Fig. 9 (a)).

linear drive γL = d ln A/dt + γd should vary proportionally with the dominant destabilizing gradient, in our caseF0 ≡ ∂Pζ F0 . Figures D5 (a,b) confirms that this expectation is approximately satisfied in our simulations without
damping. However, there is a small but noticeable deviation of about 10 . . . 13%. Furthermore, Fig. D5 (c) showsthat γL increases nearly linearly by about 15% when thedamping rate is raised to γd /ω0 = 7.5%, the value of case(A).
Another notable feature of our simulations is that thegradient F0 needed to obtain a certain value of γL depends
on the location of the resonance relative to the mode’s
peak. As one can see in Fig. D5 (d), a 3.5 times steeper
gradient is needed in case (C) in order to obtain the samevalue of γL /ω0 = 8% as in case (A) after moving the peakof the mode from ψP,peak = 0.475 to 0.54 (cf. Fig. 9 (a)).
Since its vertical axis has a logarithmic scale, Fig. D5 (d)shows a nearly exponential dependence of F0 on ψP,peak ,which may be due to the Gaussian shape of the mode’s radial profile.
Although we do not fully understand this behavior,
we suspect that it is the result of a self-optimization process, where the particle-to-field energy transfer rate is max-

1403087-52

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)
47 μs. Note that the frequency shift does no longer change
after the instability has begun to grow exponentially. The
small erratic fluctuations that can be seen in the spectrograms may be taken as a measure for the eﬀect of statisticalnoise and an imperfect quiet start (see Appendix D.5).
The measured frequency shifts in the marginally unstable case (A) and in the strongly unstable case (B) are!
+0.45 kHz : Case (A),
(D1)
Δν0 ≈
+0.15 kHz : Case (B).

Fig. D6 A prompt relaxation occurs during the first transit(4.7 μs), where the amplitude A(t) and phase φ(t) droptransiently as shown in panels (a) and (b) for case (A).
The relaxation is associated with a prompt frequency
shift Δν0 = ν − ν0 , which depends on the destabilizing gradient F0 , the damping rate γd and the location of the resonance relative to the mode peak ψP,peakas shown in panels (g) and (h). The frequency shifts
are measured using high-resolution DMUSIC spectrograms, some examples of which are shown in (c)-(f),where the initial frequency ν0 = 100 kHz is indicatedby a dashed horizontal line.

imized. In a self-consistent simulation as in Ref. [12], boththe spatial structure and the frequency of the MHD modeare optimized during the initial structure formation process. In the semi-perturbative model used here, the mode
structure is held fixed, so only the frequency can selfoptimize, depending on the local gradient F0 , on the rateof damping γd , and on the radial excursions of the GC drift
orbits across the mode’s profile.
This self-optimization hypothesis is supported by theobservation of a prompt relaxation at the beginning ofmost simulations that we have performed. As an example,
Figs. D6 (a,b) shows how the field amplitude A(t) drops
and the phase φ(t) adjusts transiently in the marginally unstable case (A). This relaxation lasts only about one transit
time τζ0 ≈ 4.7 μs and causes a prompt shift Δν0 = ν − ν0 in
the frequency of the field oscillations, away from the prescribed value of ν0 = 100 kHz. The shifted frequencies canbe seen in the high-resolution spectrograms plotted in panels (c)-(f) of Fig. D6, which were obtained with the DMUSIC algorithm using a short time window Δtwin = 10τζ0 =As shown in Fig. D6 (g), the frequency shift is reduced tozero and becomes slightly negative when the damping rate
γd is reduced to zero while keeping the values of the radial gradients F0 fixed. A reduction of the gradients to
those in cases (A0) and (B0) increases the frequency shift
to Δν0 ≈ +0.05 kHz and +0.1 kHz, respectively, as indicated by the small symbols at γd = 0 in Fig. D6 (g). Animations of the evolution of the phase space density perturbations δ f (P̂ζ , ϑ) in the undamped cases (A0) and (B0) exhibit a corresponding poloidal drift of the eﬀective O- andX-points of the resonance in the frame rotating toroidallywith the initial frequency ω0 /n. The poloidal drift velocity
is consistent with the values of Δν0 obtained using DMUSIC, so we may rule out the possibility that this is merely asignal processing artifact. Although limited numerical accuracy does aﬀect the value Δν0 , the phenomenon does notseem to be entirely numerical, at least for |Δν0 |  0.1 kHz.
As was noted at the beginning of this section, the location of the resonance relative to the mode peak matters. As
one can see in Fig. 9 (a) above, the resonance in our simulation setup is located in the region ψP,res ≈ 0.13 . . . 0.33,
whereas the peak of the mode in cases (A) and (B) is located at ψP,peak = 0.475. When we move the resonancefurther inward (away from the peak) or, equivalently, when
we move the mode’s peak outward (away from the resonance) while increasing F0 such that γL remains constant as in Fig. D5 (d), we observe that the positive frequency shift increases significantly. This is demonstrated
in Fig. D6 (h), where we scan ψP,peak up to the value of 0.54
of case (C) and obtain a prompt frequency shift as large asΔν ≈ +3 kHz, much larger than in case (A), which has thesame values of γL and γd .
The above-mentioned drift of the O- and X-points inthe undamped cases suggests that the location of the resonance has shifted. In the cases studied here, we load
particles only along a line in phase space, where E  =E − ωn0 Pζ = const., and the numerical spread around thisline is only about ΔEnum  3 × 10−3 keV (cf. Fig. D11).
This gives an upper bound on the possible resonance shift
along E  . With Eq. (26) written more conveniently asE  [kHz] ≈ E[kHz] −ν[kHz]
P̂ζ × 0.85,n(D2)

and noting that the particle energy E is eﬀectively constantduring the prompt relaxation phase due to the small field
amplitude A  10−6 , we can estimate the radial shift of the1403087-53

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

resonance using the formula
ΔP̂ζ ≈

n × ΔEnum [kHz] Δν0−P̂ζ,res .
0.85 × ν0 [kHz]
ν0

(D3)

For case (A) with Δν0 ≈ 0.4 kHz and P̂ζ,res = 0.719, we
obtain ΔP̂ζ ≈ 2 × 10−4 − 3 × 10−3 ≈ −3 × 10−3 . This islarger by a factor 3 than the observed shift ΔP̂ζ ≈ −10−3 toP̂ζ ≈ 0.718 of the resonant phase space structures in case(A). This means that the resonance does no longer lie in
the region populated by marker particles and would imply
that the field-particle interactions are occurring slightly oﬀresonantly. More likely, however, is the possibility that ourformula (D3) may not be suﬃciently accurate. We suspectthat the discrepancy may be due to the large magnetic drifts
and the radial variation of the Gaussian mode structure in
the domain occupied by the drift orbits (cf. Fig. 9).
Indeed, panels (c) and (d) in Fig. D6 show that thefrequency shift seen in case (A) with ΔE0 = 0 remains
the same in case (A’), where we have loaded marker particles in an energy band with a relatively large width ΔE0 =7.5 keV. In that case, the first term in Eq. (D3) is approximately 0.4, which would allow huge frequency shifts, but
the system does not take advantage of his possibility. This
insensitivity of Δν0 with respect to ΔE0 and the abovementioned sensitivity with respect to the relative locationsof the resonance and the mode’s peak (Fig. D6 (f)) supports
our hypothesis that the prompt relaxation is determined by
the slope of the mode structure at the resonance (possibly
in combination with drift eﬀects [42, 43]). This is further
corroborated by observations made in other case studies(not presented here) showing that (at least with suﬃciently
large γd ) the prompt frequency shift Δν0 seems to be directed towards the peak of the mode when the resonance islocated oﬀ-peak.
This subject requires further investigation. In themeantime, we choose to speak of “promptly shifted eﬀective resonances” in the present paper.

D.5

Quiet start and sensitivity

While loading particles only at the outer midplane(ϑ = 0) suﬃces for kinetic Poincaré plots, this procedure introduces spurious oscillations in the mode evolution when applied in a simulation where Eq. (15) is solved.
Spreading the simulation particles uniformly along thepoloidal angle −π ≤ ϑ ≤ π does not suﬃce either: dueto the mirror force and the large magnetic drifts that canbe seen in Fig. 9 (b), the spurious oscillations remain noticeable throughout the simulation. Their eﬀect is readily
seen in the growth rate, which becomes extremely noisy as
shown in Fig. D7. Long-term eﬀects can be observed in the
EP distribution in the form of spurious radially propagating
micro-structures as shown in Fig. D8 (b).
Although the evolution of the field remains more orless the same, such spurious oscillations are undesirable
in detailed analyses of phase space dynamics as were performed in this work. Perhaps the most accurate way toFig. D7 Noise reduction with quiet start (QS). Time traces of(a) the field amplitude A(t) and (b) the instantaneous
kinetic drive γk (t) = γ(t) + γd in the marginally unstablecase (A) simulated without and with QS. Two non-QS
simulations were performed: one with particles loaded
only on the outer mid-plane (ϑ = 0) and one with particles distributed uniformly randomized along ϑ (i.e., notconsistent with the magnetic mirror force and drifts).
The growth rates γ(t) were measured using an exponenΔA(t)+0.99×γ(t−Δt).
tial running average: γ(t) = 0.01× A(t)ΔtFig. D8 Noise reduction with quiet start (QS). Simulation of thecase studied in Ref. [13], without (left) and with quiet
start (right). Panels (a) and (c) show time traces of theamplitude A(t). The contour plots in (b) and (d) showδ f − δ f transit ϑ (P̂ζ , t) where δ f transit is δ f averaged
over one transit time τζ0 , and . . .ϑ is an average overthe poloidal angle. The quiet start smoothes the amplitude and eliminates the propagating micro-structures
seen in (b) that were reported in Fig. 11 of Ref. [13].

deal with this problem is to load particles uniformly intime along precomputed unperturbed GC orbits [68]. For
the present study, the ORBIT code has been extended witha routine that approximates this method in a simpler way
as described in Sec. 3.3. As one can see in Figs. D7 andD8 (d), the growth rates become much smoother and the
spurious micro-structures in EP phase space disappear.
Our quiet start procedure is not perfect and still gen-1403087-54

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. D9 Sensitivity test with respect to the initial perturbation
amplitude A0 . The diagrams show the early evolution of the amplitude A(t) in the marginally unstablecase (A), starting from amplitudes in the range A0 =10−6 . . . 10−9 . The log-scale panel (a) shows the nearly
exponential growth. The circles in the linearly scaled
panel (b) highlight the points of first saturation. The
dashed blue curve shows the result for A0 = 10−7 with
Np = 107 marker particles (otherwise Np = 106 ).

erates a certain amount of “statistical noise” (including unphysical modulations). We believe that this is one reason
for the results in Fig. D9, where the amplitude at the time offirst saturation in the marginally unstable case (A) changes
in a somewhat irregular manner when we reduce the initial
perturbation amplitude A0 from 10−6 to 10−9 . This case isvery sensitive to statistical noise because its proximity tomarginal stability causes the first saturation to occur at anextremely small amplitude, which corresponds to a narrow
resonant structure in EP phase space (see Figs. 4 and 11).
Statistical noise also aﬀects the chirping patterns as
can be seen in Fig. D10, where we show results of case(A) simulated with diﬀerent realizations of the quiet start(a)-(c) and with uniform particle loading (d). Our default
procedure (a) was to load a fraction of the particles at every
time step during the first 10 transits (10τζ0 = 47 μs). WithΔt/τζ0 = 1/400, this makes 4000 loading steps. In (b), we
have loaded particles only every 4th step. In (c), we loaded
particles at each step, but only for 1 transit time instead of
10. The onset of chirping is very similar in all cases, but
after about 2 ms, the patterns vary substantially, although
the overall extent of the chirps remains similar.

D.6

Fig. D10 Sensitivity test with respect to quiet start (QS) parameters. Fourier spectrograms are shown for four simulations performed with the parameters of the marginallyunstable case (A). Panels (a)-(c) show results of diﬀerent realizations of the QS, and panel (d) shows the result with particles uniformly randomized along ϑ (nonQS). The size of the sliding FFT time window was
Δtwin = 0.47 ms. The green dotted parabola represents
the BB model (12) and is shown here only as a reference for orientation and easier comparison of patterns
in diﬀerent simulations.

Numerical convergence

The ORBIT code uses a 4th-order Runge-Kutta algorithm to push particles and evolve the fields, so there isa numerical error associated with the time step Δtstep thatcauses the energy conservation condition to be violated bya small amount. Figure D11 shows snapshots of the particle distribution in the (P̂ζ , E  )-plane taken at a very earlystage of a simulation of case (A) with quiet start, after onlytˆ = 5 transit times τζ0 . Results are shown for simulations
performed with τζ0 /Δtstep = 100, 400 and 1000 steps pertransit, where 400 is the default value used in this paper.
One can see that the spread in E  is δE  /E0 ≈ 4 × 10−4for τζ0 /Δtstep = 100, which decreases by an order of mag-Fig. D11 Numerical convergence with respect to time step Δtstep .
The particle distribution in the (P̂ζ , E  )-plane obtained
after 5 transits is shown for simulations of case (A)with τζ0 /Δtstep = 100, 400, 1000 steps per transit. The
upper diagram shows an enlarged view of the latter twocases.

nitude to δE  /E0 ≈ 4 × 10−5 for τζ0 /Δtstep = 400, and isnearly invisible on the scales shown for τζ0 /Δtstep = 1000.
Note that the spread occurs only along the E  axis ofFig. D11. It does not aﬀect P̂ζ , so it is entirely limited to1403087-55

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Fig. D12 Numerical convergence with respect to time step Δtstep .
Fourier spectrograms are shown for cases (A) and (B)obtained in simulations with τζ0 /Δtstep = 1000 steps
per transit. The results are similar to those in Fig. 28
obtained with 400 steps per transit. The green dottedparabola represents the BB model (12) and is shown
here only as a reference for easier comparison between
diﬀerent cases.

a (small) variation in the particle energy E. Moreover, theenergy spread appears during the first few transit times ofa simulation and remains more or less at the same level
thereafter. This can be seen in the lower part P̂ > 0.65
of Fig. 29 (a). (The broadening of the E  distribution in
the upper part P̂  0.65 of Fig. 29 (a) has other reasons
discussed in Sec. 5.4 and Appendix D.2.) In other words,
the energy conservation error is bounded and the bound is
set by the size of the time step Δtstep .
Figure D12 shows Fourier spectrograms of the longtime evolution in cases (A) and (B) simulated withτζ0 /Δtstep = 1000 steps per transit instead of the default
400, spectrograms for which were shown in Fig. 28. One
can see that the results are very similar; especially, in thestrongly driven case (B), which is less sensitive to statistical noise than the marginally unstable case (A). The rapid
downward chirps in the lower part of both spectrograms in
Fig. D12 are the spurious “ghosts” discussed in Sec. C.4.
A higher (or more clearly visible) sensitivity with respect to the time step Δtstep is found in simulations of undamped instabilities (γd = 0), where a quasi-steady state is
expected to form. Figure D13 (a) for the weakly driven
case (A0) and Fig. D13 (c) for the strongly driven case(B0) show that simulations with too large time steps failto reproduce the saturation of the kinetic instabilities. In
fact, this is a well-known problem and a standard method
to test codes and choose suitable time steps. In addition,
Fig. D13 (a) shows that unphysical fluctuations in simulations without quiet start can enhance the saturation problem.
In case (A0), taking τζ0 /Δtstep = 1000 steps per transit

Fig. D13 Numerical convergence test with respect to time stepΔtstep and particle number Np . The diagrams show time
traces of the field amplitude A(t) in simulations without damping (γd = 0). Panels (a) and (b) show results
for the weakly driven case (A0) simulated with diﬀerent time steps Δtstep and diﬀerent numbers of marker
particles Np , respectively. Panel (a) also shows results
obtained with and without quiet start (QS). Panel (c)show results for the strongly driven case (B0) simulated with Np = 106 markers and diﬀerent time steps
Δtstep .

Fig. D14 Calibration: Relation between the oscillation frequency ν and the radial location P̂ζ in EP phase space.
Formulas based on a linear fit are given for the absolute value ν(P̂ζ ) and for the distance δν(δP̂ζ ) from theseed resonance.

1403087-56

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

yields good saturation as one can see in Fig. D13 (a). Even
400 steps per transit yield an acceptable result, with a very
similar oscillation pattern. Figure D13 (b) shows that the
number of simulation particles can influence the oscillation
pattern but seems to have no significant influence on thesaturation (or lack thereof).
In case (B0), simulations with τζ0 /Δtstep = 1000 steps
per transit still suﬀer from spurious growth as one can see
in Fig. D13 (c). However, the results for the first millisecond or so, which we have analyzed in Sec. 4.2 and Appendix B, seem to be acceptable.

D.7

Calibration

Ideal MHD displacement. The magnitude of the ideal
MHD displacement δr = ξ · ∇r ≈ ξΨ /Ψ = ξΨ /(qΨP ) withΨP ≡ ∂r ΨP of the (nonresonant) bulk plasma can be estimated as follows. Using the relation ψP = ΨP /ΨP,edge ≈(r/a)2 , making explicit the normalization of Ψ̂P,edge =ΨP,edge /B0 , and recalling that ξ̂Ψ = |ξΨ |/A is normalized
to unity (so that A has units of length squared), we obtain
Aξ̂Ψ
Aξ̂Ψ
a/R0
δrmhd
,
=
≈
×

R0
R0 qΨP /B0 qΨP,edge /B0 2r0 /a


1/3.7
10−3
×
= 7 × 10−4 ,1.15 × 0.276
1.3

(D4)

for A  10−3 . Meanwhile, the kinetic Poincaré plot in
Fig. 9 (b) shows an island with a half-width of about δψP ∼±0.03 for A = 10−3 . This means that the resonant particles
experience displacements
0.03
δrres
δψP
≈
≈ 6 × 10−3 , (D5)≈R0
2r0 /a×R0 /a 1.3×3.7
that are one order of magnitude larger than δrmhd /R0 for
the nonresonant bulk.

Profile of EP reference state F0 . The slope of the reference EP density profile F0 (Pζ ) is not known precisely,
since our marker loading method gives only an approximately uniform density gradient and, thus, we do not know
the exact form of F0 /G0 ≈ const. The calibration of the
EP profile was performed using Fig. 4 (c), which shows
the first snapshot of case (A’) described in Appendix D.3,
where we have allowed the EP distribution to adapt to thefield perturbation with fixed amplitude A = A0 = 10−6 andphase φ = 0. The O-point profile of such a relaxed distribution F̄(Pζ ) = F0 (Pζ ) + δ f¯(Pζ ) is expected to be flat and wehave calibrated the reference profile F0 (Pζ ) accordingly.
The resulting profile F0 (Pζ ) is shown as a yellow shaded
triangle in Fig. 4 (c). The calibration was also confirmed
independently using the undamped cases (A0) and (B0),
where the perturbed profile F̄(Pζ ) is also expected to become flat across the O-points (see Figs. 14 and 15).
Pump frequency. The relation between the radial location P̂ζ of a phase space structure and the corresponding
oscillation frequency ν in a spectrogram of the field signals(t) is shown in Fig. D14. The data points in Fig. D14 wereroughly taken at the approximate centers of large hole andclump vortices in the advanced stages of the simulation,
which is why there are no data points near the resonanceP̂ζ,res ≈ 0.72. A linear fit givesδP̂ζ ≈ δν/486.8 kHz.

(D6)

Appendix E. Characterization of Diagnostics
E.1 Fourier spectra
The Fourier analysis was performed using a Hann-

Fig. E1 Comparison between FFT and DMUSIC spectrograms for case (A). Window sizes Δtwin = 0.47 ms (left) and 0.094 ms (right)
were used. For better comparability, the FFT spectra were normalized to their maximal value at each time, so intensity information is lost (as in DMUSIC). The panels at the bottom show enlarged portions of the full spectrograms in panels (a)-(d).

1403087-57

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)
weighted sliding time window. The default window size
was Δtwin = 100 τζ0 ≈ 0.47 ms and the spectrum has beensmoothed by extending the array length to T pad = 16×Δtwin
and padding with zeros.
T
pad /2Y(t, ω) =dτ y(t − τ)H(τ)e−iωτ .

(E1)−T pad /2

The Hann window function H(τ) eliminates artifacts fromaperiodicity.

E.2 DMUSIC spectra

Fig. E2 Characterization of DMUSIC spectrograms. Using thesignal s(t) from case (A), whose envelope A(t) is plotted in panel (a), we show the form of the spectrogram
during the onset of chirping for short time windows
Δtwin = 0.096 ms or 0.047 ms, and for diﬀerent numbers K = 3, 5 and 21 of complex sinusoids used in thedecomposition.

DMUSIC is a spectral peak-finding algorithm originally developed for nuclear magnetoresonance (NMR) signal analysis by Li et al. [69]. Our application of the method
to fusion plasmas was inspired by Slaby et al. [70], whoused it to visualize continuous spectra of shear Alfvén
waves in initial value simulations.
The strength of the algorithm lies in its applicability tonoisy signals with phases of rapid growth or strong damping. Its weaknesses are that signal intensity information iseﬀectively lost and that it can be diﬃcult to attribute physical meaning to all the features appearing in a DMUSIC
spectrogram. The spectral analyses in this study were performed with sliding time windows with sizes in the range
Δtwin = (10 . . . 300)τζ0 = (0.047 . . . 1.410) ms. For eachvalue of the real frequency ν, we scanned the imaginary
frequency axis γ and selected the maximal value of thetransform. At each time, the peak value of the spectrogram
was normalized to unity.
Figure E1 shows a comparison between FFT spec-

Fig. E3 Click on panels (a)-(d) to view animated movies for cases (A), (B) and (B0). Playing the animations directly in the PDF file may
require Adobe Acrobat Reader and VLC to be installed. On some systems, a separate media player may be launched.

1403087-58

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

trograms and DMUSIC spectrograms obtained with longand short time windows. With the long time window
Δtwin = 0.47 ms, the spectra have a similar structure, although the DMUSIC spectra are much sharper and seem toreveal more detail. It is not clear how much of that detailbears physical meaning. For instance, the splitting of spectral lines seen in the upper part of Fig. E1 (b) (enlarged in
the third panel in the bottom row) remains to be understoodas discussed in Appendix C.1. We use DMUSIC primarily
for the purpose of spectral analysis with high temporal resolution (d), where FFT spectra become unintelligible (c).
Our default parameters for DMUSIC are as follows.
We take a portion of N = Δtwin /Δtsample data points froma uniformly sampled signal, where Δtsample is the samplingstep. We choose N to be even, and the parameter J of Ref.
[69] is chosen to be J = N/2 as recommended by the authors. Typically, we use K = 21 complex sinusoids for thedecomposition. Typically, a real frequency window covering the range 40 kHz ≤ ν ≤ 180 kHz is sampled by at leastNν = 201 points. The window size for the imaginary component of the frequency is usually −4 kHz ≤ γ ≤ 4 kHz
and sampled by Nγ = 31 points (tests with Nγ = 121 gave
essentially identical results).
The role of the parameters K and Δtwin can be seen in
Fig. E2, which shows spectrograms of the onset of chirping during the first 2 ms of case (A). The spectrograms arecomputed with short time windows Δtwin = 0.094 ms and0.047 ms. When using a small number of only K = 3 complex sinusoids (b), the spectrogram eﬀectively shows theinstantaneous oscillation frequency. Phase jumps by ±πthat occur around the amplitude minima and have the formof an arctan function (see Eq. (A6b)) appear here in panel(b) as rapid up- and downward chirps. Frequency splitting
becomes visible with K = 5 in panel (c), and becomes only
a little richer in detail with K = 21 in (d). With a smaller
time window as in panel (e) and (f), the phase jumps start tostrongly aﬀect the signals even for K = 5 and 21. Further
reduction of the window size would yield a result similar
to that in panel (b) since the number of harmonics that canbe captured decreases, so choosing a large value of K has
no eﬀect for short time windows.

E.3

Animations

Figure E3 contains links to animations for cases (A),(B) and (B0). Each movie window consists of three panels: one showing a δ f -weighted Poincaré plot of the phasespace structures (left), one for the field amplitude A(t) (topright), and a Fourier spectrogram (bottom right).
[1] L. Chen and F. Zonca, Physics of Alfvén waves and energetic particles in burning plasmas, Rev. Mod. Phys. 88,
15008 (2016).
[2] Y. Todo, Introduction to the interaction between energetic
particles and Alfvén eigenmodes in toroidal plasmas, Rev.
Mod. Plasma Phys. 3, 1 (2019).
[3] W.W. Heidbrink and R.B. White, Mechanisms of energetic-

particle transport in magnetically confined plasmas, Phys.
Plasmas 27, 030901 (2020).
[4] F. Zonca, L. Chen, S. Briguglio, G. Fogaccia, G. Vlad andX. Wang, Nonlinear dynamics of phase space zonal structures and energetic particle physics in fusion plasmas, New
J. Phys. 17(1), 013052 (2015).
[5] G. Meng, N.N. Gorelenkov, V.N. Duarte, H.L. Berk, R.B.
White and X.G. Wang, Resonance frequency broadeningof wave-particle interaction in tokamaks due to Alfvénic
eigenmode, Nucl. Fusion 58(8), 082017 (2018).
[6] C.T. Hsu, C.Z. Cheng, P. Helander, D.J. Sigmar and R.
White, Particle dynamics in chirped-frequency fluctuations, Phys. Rev. Lett. 72(16), 122506 (1994).
[7] R.L. Dewar, Saturation of kinetic plasma instabilities byparticle trapping, Phys. Fluids 16(3), 431 (1973).
[8] F. Zonca and L. Chen, Destabilization of energetic particle
modes by ICRF-induced fast minority ion tails on TFTR, In
JAERI-Conf 2000-004 (6th IAEA TCM on EP in Magnetic
Confinement Systems, 1999, Naka, Japan), pages 52-56.
JAERI, 2000.
[9] W.W. Heidbrink, M.A Van Zeeland, M.E. Austin, N.A.
Crocker, X. Du, G.R. McKee and D.A. Spong, Stability ofbeta-induced Alfvén eigenmodes (BAE) in DIII-D, Nucl.
Fusion 61(6), 066031 (2021).
[10] H.L. Berk, B.N. Breizman and N.V. Petviashvili, Spontaneous hole-clump pair creation in weakly unstable plasmas,
Phys. Lett. A 234, 213 (1997), Erratum: Phys. Lett. A 238,
408 (1998).
[11] F. Zonca, Private communication concerning a theory ofauroral chorus generation, Article in preparation.
[12] A. Bierwage, K. Shinohara, Y. Todo, N. Aiba, M. Ishikawa,
G. Matsunaga, M. Takechi and M. Yagi, Self-consistentlong-time simulation of chirping and beating energetic particle modes in JT-60U plasmas, Nucl. Fusion 57(1), 016036(2017).
[13] R.B. White, V.N. Duarte, N.N. Gorelenkov, E.D.
Fredrickson and M. Podesta, Phase-space dynamics ofAlfvén mode chirping, Phys. Plasmas 27, 052108 (2020).
[14] H.K. Berk and K.V. Roberts, Nonlinear study of Vlasov’s
equation for a special class of distribution functions, Phys.
Fluids 10, 1595 (1967).
[15] P. Khain and L. Friedland, A waterbag theory of autoresonant Bernstein-Greene-Kruskal modes, Phys. Plasmas 14(8), 082110 (2007).
[16] H. Hezaveh, Z.S. Qu, M.J. Hole and R.L. Dewar, Theoretical description of chirping waves using phase-spacewaterbags, Plasma Phys. Control. Fusion. 63(6), 065008(2021).
[17] R.B. White and M.S. Chance, Hamiltonian guiding center
drift orbit calculation for plasmas of arbitrary cross section,
Phys. Fluids 27(10), 2455 (1984).
[18] Y. Chen, R.B. White, G.-Y. Fu and R. Nazikian, Numerical study of the nonlinear evolution of toroidicity-induced
Alfvén eigenmodes, Phys. Plasmas 6, 226 (1999).
[19] R.B. White, The Theory of Toroidally Confined Plasmas,(Imperial College Press, London, 3rd edition, 2014).
[20] K.V. Roberts and H.L. Berk, Nonlinear evolution of a twostream instability, Phys. Rev. Lett. 19(6), 297 (1967).
[21] T.P. Armstrong, Numerical studies of the nonlinear Vlasov
equation, Phys. Fluids 10(6), 1269 (1967).
[22] G. Wang, H.L. Berk, B.N. Breizman and L.-J. Zheng, Frequency chirping in the Alfvén continuum, Nucl. Fusion
58(8), 082014 (2018).
[23] I.B. Bernstein, J.M. Greene and M.D. Kruskal, Exact non-

1403087-59

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

linear plasma oscillations, Phys. Rev. 108(3), 546 (1957).
[24] H.L. Berk, B.N. Breizman, J. Candy, M. Pekker and N.V.
Petviashvili, Spontaneous hole-clump pair creation, Phys.
Plasmas 6, 3102 (1999).
[25] B.N. Breizman, Nonlinear travelling waves in energetic
particle phase space, Nucl. Fusion 50(8), 084014 (2010).
[26] R.M. Nyqvist, M.K. Lilley and B.N. Breizman, Adiabatic
description of long range frequency sweeping, Nucl. Fusion 52(9), 094020 (2012).
[27] R.M. Nyqvist and B.N. Breizman, Modeling of long range
frequency sweeping for energetic particle modes, Phys.
Plasmas 20(4), 042106 (2013).
[28] H. Hezaveh, Z.S. Qu, B. Layden and M.J. Hole, Impact ofenergetic particle orbits on long range frequency chirpingof BGK modes, Nucl. Fusion 57(12), 126010 (2017).
[29] H. Hezaveh, Z.S. Qu, B.N. Breizman and M.J. Hole, Long
range frequency chirping of Alfvén eigenmodes, Nucl. Fusion 60(5), 056014 (2020).
[30] L. Chen, Theory of magnetohydrodynamic instabilities excited by energetic particles in tokamaks, Phys. Plasmas
1(5), 1519 (1994).
[31] A. Bierwage and K. Shinohara, Orbit-based analysis of resonant excitations of Alfvén waves in tokamaks, Phys. Plasmas 21(11), 112116 (2014).
[32] A. Fasoli, D. Borba, C. Gormezano, R. Heeter, A. Jaun,
J. Jacquinot, W. Kerner, Q. King, J.B. Lister, S. Sharapov,
D. Start and L. Villard, Alfvén eigenmode experiments intokamaks and stellarators, Plasma Phys. Control. Fusion,
39, B287 (1997).
[33] A. Fasoli, B.N. Breizman, D. Borba, R.F. Heeter, M.S.
Pekker and S.E. Sharapov, Nonlinear splitting of fast particle driven waves in a plasma: Observation and theory,
Phys. Rev. Lett. 81(25), 5564 (1998).
[34] Y. Todo, M.A. Van Zeeland, A. Bierwage, W.W. Heidbrink
and M.E. Austin, Validation of comprehensive magnetohydrodynamic hybrid simulations for Alfvén eigenmode induced energetic particle transport in DIII-D plasmas, Nucl.
Fusion 55(7), 073020 (2015).
[35] Y. Todo, M.A. Van Zeeland and W.W. Heidbrink, Fast ion
profile stiﬀness due to the resonance overlap of multiple
Alfvén eigenmodes, Nucl. Fusion 56(11), 112008 (2016).
[36] A. Bierwage, K. Shinohara, Y. Todo, N. Aiba, M. Ishikawa,
G. Matsunaga, M. Takechi and M. Yagi, Simulations tackle
abrupt massive migrations of energetic beam ions in a tokamak plasma, Nature Comms. 9, 3282 (2018).
[37] R.B. White, V.N. Duarte, N.N. Gorelenkov, E.D.
Fredrickson, M. Podesta and H.L. Berk, Modeling of chirping toroidal Alfvén eigenmodes in NSTX, Phys. Plasmas
26, 092103 (2019).
[38] S.D. Pinches, L.C. Appel, J. Candy, S.E. Sharapov,
H.L. Berk, D. Borba, B.N. Breizman, T.C. Hender, K.I.
Hopcraft, G.T.A. Huysmans and W. Kerner, The HAGIS
self-consistent nonlinear wave-particle interaction model,
Comp. Phys. Comm. 111, 133 (1998).
[39] S.D. Pinches, H.L. Berk, M.P. Gryaznevich, S.E. Sharapov
and JET-EFDA Contributors, Spectroscopic determinationof the internal amplitude of frequency sweeping TAE,
Plasma Phys. Control. Fusion 46(7), S47 (2004).
[40] Y. Todo, H.L. Berk and B.N. Breizman, Simulation of intermittent beam ion loss in a Tokamak Fusion Test Reactor
experiment, Phys. Plasmas 10(7), 2888 (2003).
[41] S. Nishimura, Y. Todo, D.A. Spong, Y. Suzuki and N.
Nakajima, Simulation study of Alfvén-eigenmode-induced
energetic ion transport in LHD, Plasma Fusion Res. 8(7),2403090 (2013).
[42] A. Bierwage and K. Shinohara, Orbit-based analysis ofnonlinear energetic ion dynamics in tokamaks. I. Eﬀective mode number profile and resonant frequency tracking,Phys. Plasmas 23(4), 042511 (2016).
[43] A. Bierwage and K. Shinohara, Orbit-based analysis ofnonlinear energetic ion dynamics in tokamaks. II. Mechanisms for rapid chirping and convective amplification,
Phys. Plasmas 23(4), 042512 (2016).
[44] T. Wang, Z. Qiu, F. Zonca, S. Briguglio and G. Vlad, Dynamics of reversed shear Alfvén eigenmode and energetic
particles during current ramp-up, Nucl. Fusion 60, 126032(2020).
[45] F. Zonca and L. Chen, Resonant damping of toroidicityinduced shear-Alfvén eigenmodes in tokamaks, Phys. Rev.
Lett. 68(5), 592 (1992).
[46] J.R. Cary and A.J. Brizard, Hamiltonian theory of guidingcenter motion, Rev. Mod. Phys. 81, 693 (2009).
[47] R.B. White and A. Bierwage, Particle resonances intoroidal fusion devices, Phys. Plasmas 28(3), 032507(2021).
[48] C.T. Hsu and D.J. Sigmar, Alpha-particle losses fromtoroidicity-induced Alfvén eigenmodes. Part I: Phase-space
topology of energetic particle orbits in tokamak plasma,
Phys. Fluids B: Plasma Physics 4, 1492 (1992).
[49] R.F. Heeter, A.F. Fasoli and S.E. Sharapov, Chaotic regimeof Alfvén eigenmode wave-particle interaction, Phys. Rev.
Lett. 85(15), 3177 (2000).
[50] F. Zonca, S. Briguglio, L. Chen, G. Fogaccia and G. Vlad,
Collective eﬀects and self-consistent energetic particle dynamics in advanced tokamaks, In Proc. 19th Fusion Energy
Conference, Lyon, France, 14-19 Oct. 2002, pages IAEACSP-19/CD/TH/4-4, Vienna, 2003. IAEA.
[51] R.B. White, R.J. Goldston, K. McGuire, A.H. Boozer, D.A.
Monticello and W. Park, Theory of mode-induced beam
particle loss in tokamaks, Phys. Fluids 26, 2958 (1983).
[52] M.J. Hole and L.C. Appel, A modulation model for mode
splitting of magnetic perturbations in the Mega Ampere
Spherical Tokamak, Plasma Phys. Control. Fusion 51,
045002 (2009).
[53] W. Chen, M. Jiang, Y. Xu, P.W. Shi, L.M. Yu, X.T. Ding,
Z.B. Shi, X.Q. Ji, D.L. Yu, Y.G. Li, Z.C. Yang, W.L.
Zhong, Z.Y. Qiu, J.Q. Li, J.Q. Dong, Q.W. Yang, Yi. Liu,
L.W. Yan, M. Xu and X.R. Duan, Experimental observationof multi-scale interactions among kink/tearing modes andhigh-frequency fluctuations in the HL-2A core NBI plasmas, Nucl. Fusion 57(11), 114003 (2017).
[54] J. Zhu, Z.W. Maa, S. Wang and W. Zhang, Nonlinear dynamics of toroidal Alfvén eigenmodes in the presence oftearing modes, Nucl. Fusion 58, 046019 (2018).
[55] I. Barth and L. Friedland, Quantum phenomena in a chirped
parametric anharmonic oscillator, Phys. Rev. Lett. 113(4),040403 (2014).
[56] V.I. Veksler, A new method of acceleration of relativistic
particles, J. Phys. Acad. Sci. U.S.S.R. 9, 153 (1945).
[57] E.M. McMillan, The Synchrotron—a proposed high energy
particle accelerator, Phys. Rev. 68(1), 143 (1945).
[58] B. Chapman, R.O. Dendy, S.C. Chapman, K.G.
McClements, G.S. Yun, S.G. Thatipamula and M.H.
Kim, Nonlinear wave interactions generate high-harmonic
cyclotron emission from fusion-born protons during a
KSTAR ELM crash, Nucl. Fusion 58(9), 096027 (2018).
[59] G.J. Kramer, B.J. Tobias, A. Turnbull and E.M. Bass, Enhanced radial energy transport induced by radially curved

1403087-60

Plasma and Fusion Research: Regular Articles

Volume 16, 1403087 (2021)

Alfvén eigenmode wavefronts, Nucl. Fusion 59(9), 094001(2019).
[60] G. Meng, Ph. Lauber, Z. Lu and X. Wang, Eﬀects of the
non-perturbative mode structure on energetic particle transport, Nucl. Fusion 60, 056017 (2020).
[61] F. Zonca, L. Chen and R.A. Santoro, Kinetic theory of lowfrequency Alfvén modes in tokamaks, Plasma Phys. Control. Fusion 38, 2011 (1996).
[62] F. Zonca, L. Chen, R.A. Santoro and J.Q. Dong, Existenceof discrete modes in an unstable shear Alfvén continuous
spectrum, Plasma Phys. Control. Fusion 40, 2009 (1998).
[63] Ph. Lauber, S. Günter, A. Könies and S.D. Pinches,
LIGKA: A linear gyrokinetic code for the description of
background kinetic and fast particle eﬀects on the MHD
stability in tokamaks, J. Comp. Phys. 226, 447 (2007).
[64] A. Bierwage and Ph. Lauber, Gyrokinetic analysis of low-n
shear Alfvén and ion sound waves spectra in a high-betatokamak plasma, Nucl. Fusion 57(11), 116063 (2017).
[65] A. Könies, C. Slaby, R. Kleiber, T. Fehér, M. Borchardt andA. Mishchenko, The MHD continuum with a radial electricfield, Phys. Plasmas 27(12), 122511 (2020).
[66] R.B. White, V.N. Duarte, N.N. Gorelenkov, A. Bierwage,

E.D. Fredrickson, M. Podesta and H.L. Berk, Alfvén mode
chirping, In Proc. 28th Fusion Energy Conference, Virtual
Event, 10-15 May 2021, pages TH/P1-2, Vienna, 2021.
IAEA.
[67] M.K. Lilley and R.M. Nyquist, Formation of phase spaceholes and clumps, Phys. Rev. Lett. 112(15), 155002 (2014).
[68] A. Bierwage, C. Di Troia, S. Briguglio and G. Vlad, Orbitbased representation of equilibrium distribution functions for low-noise initialization of kinetic simulations oftoroidal plasmas, Comp. Phys. Comm. 183, 1107 (2012).
[69] Y. Li, J. Razavilar and K.J.R. Liu, A high-resolution
technique for multidimensional NMR spectroscopy, IEEE
Trans. Biomed. Eng. 45(1), 78 (1998).
[70] C. Slaby, S. Äkäslompolo, M. Borchardt, J. Geiger, R.
Kleiber, A. Könies, S. Bozhenkov, C. Brandt, A. Dinklage,
M. Dreval, O. Ford, G. Fuchert, D. Hartmann, M. Hirsch,
U. Höfel, Z. Huang, P. McNeely, N. Pablant, K. Rahbarnia,
N. Rust, J. Schilling, A. von Stechow, H. Thomsen andthe Wendelstein 7-X team, Investigation of mode activity
in NBI-heated experiments of Wendelstein 7-X, Nucl. Fusion 60(11), 112004 (2020).

1403087-61