# Nonlinear equilibria and transport processes in burning plasmas

Matteo Valerio Falessi1,2, Liu Chen3,4,1, Zhiyong Qiu3,1 and Fulvio Zonca1,3
1Center for Nonlinear Plasma Science and C.R. ENEA Frascati, C.P. 65, 00044
Frascati, Italy
2Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Roma, Piazzale Aldo
Moro 2, Roma, Italy
3Institute of Fusion Theory and Simulation, School of Physics, Zhejiang
University, Hangzhou, China
4Department of Physics and Astronomy, University of California, Irvine CA
92697-4575, USA<EMAIL_ADDRESS>
###### Abstract

In this work, we put forward a general phase-space transport theory in
axisymmetric tokamak plasmas based upon the concept of zonal state (ZS).
Within this theoretical framework, the ZS corresponds to a renormalized plasma
nonlinear equilibrium consisting of phase-space zonal structures (PSZS) and
zonal electromagnetic fields (ZFs) which evolve self-consistently with
symmetry breaking fluctuations and sources/collisions. More specifically, our
approach involves deriving governing equations for the evolution of particle
distribution functions (i.e, PSZS), which can be used to compute the
corresponding macro-/meso-scale evolving magnetized plasma equilibrium
adopting the Chew Goldberger Low (CGL) description, separating the
spatiotemporal microscale structures. The nonlinear physics of ZFs and of
geodesic acoustic modes/energetic particle driven geodesic acoustic modes is
then analyzed to illustrate the implications of our theory.

##  1 Introduction

One fundamental aspect of magnetized fusion plasmas is the study of physics
processes that are occurring in burning conditions, where $\alpha$ particles
are produced, e.g., by deuterium-tritium fusion reactions. Understanding the
behavior of $\alpha$ particles and, more generally, of energetic particles
(EPs)111In a broader sense, when we refer to energetic particles (EPs), we
mean not only charged fusion products, but also suprathermal ions and
electrons generated by external power sources for plasma heating and current
drive. in fusion plasmas is crucial, as they play a critical role in mediating
couplings between microscopic and macroscopic plasma dynamics significantly
influencing its behavior as a complex system [1, 2, 3]. Reactor relevant
plasmas will be fundamentally different with respect to those existing in
present day experimental devices due to EP characteristic orbit size and
predominant contribution to reactor power balance.

A thorough description of EP transport processes is essential since they
involve resonantly excited fluctuations, which have different time scales and
structures with respect to thermal plasma instabilities. EPs may induce non-
local (global) behaviors as well as excite singular radial mode structures at
shear Alfvén wave continuum resonances which, through mode conversion,
generate radially propagating microscopic fluctuations that may be absorbed at
different radial locations. For these reasons, a global approach is necessary,
where the entire plasma volume is described taking into account realistic
magnetic geometry and non-uniformities. Due to the characteristic features of
EP distribution functions in the velocity space, the excitation of collective
fluctuations around the cyclotron frequency is usually of minor importance [2]
and, therefore, the relevant resonant frequencies are those characterizing the
particle guiding center motion in the equilibrium magnetic fields.
Consequently, current research on EP-driven instabilities and transport is
generally pursued using nonlinear gyrokinetic theory [4, 5, 6]. Numerical
simulations of EP physics usually require costly and time-consuming global
gyrokinetic codes. While these simulations provide valuable insight into
fundamental physics processes, they typically cover only a relatively limited
time range of the dynamics. Thus, they typically do not properly address the
coupling of different spatio-temporal scales and, hence, have limited
predictive capability of transport processes.

To overcome some of the challenges faced in the description of multi-scale
burning plasmas dynamics, we have developed the transport theory of phase-
space zonal structures (PSZS) [7, 8]. PSZS represent slowly evolving (on the
transport time scale) structures in the phase-space that are not affected by
fast collisionless dissipation and provide a proper definition of plasma
nonlinear equilibrium distribution function extending the concept of plasma
transport processes to the phase-space [7]. This is particularly relevant for
weakly collisional plasmas that exhibit significant deviations from local
thermodynamic equilibrium, often described by Maxwellian distribution
functions. Notably, using the PSZS theory, the usual plasma transport
equations can be obtained as a particular limiting case where the deviation
from the local Maxwellian is small. This was demonstrated in a previous work
[7] for the energy and density transport equations. However, in the most
general case, PSZS will not be associated with a reference Maxwellian since
they will be the result of the competition between resonantly induced
nonlinear transport, sources and only weakly collisional effects, thus
requiring a phase-space description. The consistency of the PSZS theory with
“gyrokinetic transport theory” [9, 10, 11] and global gyrokinetic codes stems
from its foundation in the well-established nonlinear gyrokinetics [4, 5], as
emphasized in [7]. The novelty, however, stands in explicitly identifying the
part of the toroidally symmetric distribution function that must be
incorporated in the non-Maxwellian reference state. Over time, in fact, this
contribution may increasingly deviate from the reference equilibrium due to
nonlinear processes eventually invalidating the usual transport analyses that
rely on local Maxwellian equilibria, for instance. This is evidently the case
of transport processes of weakly collisional EPs. For example, the role of
PSZS in EP transport due to energetic particle modes (EPM) [12] was recently
investigated by means of gyrokinetic particle-in-cell simulations [13]. There,
the linear scaling of the chirping rate with mode amplitude of nonlinear
coherent EPM fluctuations is numerically demonstrated, consistent with
theoretical analyses, which predict the ballistic propagation of the PSZS [1,
2, 14]. However, one may argue that even the thermal component of magnetized
plasmas, for which Coulomb collision tend to continuously restore the nearly
Maxwellian local thermodynamic equilibrium, may have significantly different
evolutions on the long time scales if the small but finite deviation from
local Maxwellian is accounted for. In other words, the thermal plasma may
evolve into nearly Maxwellian equilibria with completely different radial
profiles (cf., e.g., the recent work in Ref. [15]).

In this work, we first derive the PSZS evolution equation in conservative form
using the equilibrium constants of motion as phase-space coordinates. The
orbit averaging approach, adopted here, has analogies with the methodologies
that are used for neoclassical transport studies in stellarators in the weakly
collisional regimes [16]; and in Hamiltonian formulations of quasi-linear
transport [17, 18]. However, our present approach takes into account that EP
induced transport processes may be induced by a quasi-coherent fluctuation
spectrum of non-perturbative nature and characterized by ${\cal O}(1)$ Kubo
numbers [1, 2, 3], invalidating fundamental assumptions of quasi-linear
theory. Meanwhile, non-local processes that question the Ansatz of local
diffusive transport are also accounted for in the present analysis. In fact,
by solving the PSZS transport equations through a hierarchy of transport
models [8] ranging from global gyrokinetics [19] to quasilinear theories, we
can develop and validate advanced reduced EP transport models capable of
capturing the long timescale evolution of burning plasmas, and provide
insights into the non-locality of the underlying transport processes [8].

In order to further articulate and discuss the present gyrokinetic theory of
transport in phase-space, we also represent the PSZS evolution equation in the
magnetic-drift/banana center frame using standard flux coordinates and the
relative shift operator accounting for the guiding center magnetic drift
motion in the slowly evolving equilibrium. The corresponding self-consistent
modifications to the reference magnetic equilibrium can be obtained applying
the Chew Goldberger Low (CGL) description [20] by means of the macro-/meso-
scopic component of the PSZS moments [21]. The renormalized nonlinear
equilibrium evolving on the transport time scale due to self-consistent
fluctuations, sources and collisions is described by the zonal state (ZS) [7,
8], consistent with the PSZS transport theory. The ZS, thus, consists of the
PSZS and its counterpart, i.e., the zonal electromagnetic fields (ZFs), which
represent the long-lived component of toroidally symmetric electromagnetic
fields. In fact, the ZS does not evolve in the absence of fluctuations and/or
sources and collisions, which is consistent with its definition as a proper
nonlinear equilibrium. A more rigorous definition of the ZS is given below in
Section 2.1. The evolution of the ZS discussed in this work expands upon the
results in Ref. [7] and is predominantly due to toroidal symmetry breaking
fluctuations. Here, as a further step, we derive an expression for the plasma
polarizability that generalizes the expressions derived recently to arbitrary
geometry and equilibrium distribution functions, i.e., PSZS [22, 23, 24, 25].
We also show that transport equation can be cast in the form of a flux surface
averaged continuity equation including neoclassical transport in the banana
regime as well as sources/sinks. An in-depth discussion of phase-space
transport processes due to toroidally symmetry breaking perturbations as well
as sources/sinks will be reported in a separate work, where we will also
address the possibility of constructing reduced transport models within a
unified theoretical framework.

The article is structured as follows. In Section 2, we introduce the concept
of ZS based on the notion of PSZS and ZFs, which is explored in more detail in
Section 2.1. Next, in Section 2.2, we demonstrate how PSZS can be interpreted
as a renormalization of the reference distribution function in the presence of
a finite level of fluctuations. In Section 2.3, we explore the self-consistent
modification of the reference magnetic equilibrium due to the PSZS. Section 3
focuses on the self-consistent evolution of the ZS, showing how a gyrokinetic
transport theory on long time scales can be consistently developed within the
present theoretical framework adopting the conservative form of nonlinear
gyrokinetic equations and reconnecting to the previous work discussed in Ref.
[7]. Finally, we summarize our findings and discuss future directions in
Section 4. As further illustrative example of applications of the present
theoretical framework, A presents to interested readers a detailed discussion
of the physics of Geodesic acoustic mode (GAM)/ Energetic particle driven
geodesic acoustic mode (EGAM) in general geometry. Although GAM/EGAM do not
belong to the ZS as nonlinear equilibrium due to their finite frequency and
collisionless damping/drive, their nonlinear dynamics can affect the ZS
nonlinear evolution in a peculiar fashion. Thus, they are presented here as
particular yet paradigmatic case.

##  2 Phase-space zonal structures and the zonal state

As already mentioned in the Introduction, PSZS are characterized by being
“slowly evolving” which means that they are not affected by collisionless
dissipation, e.g., Landau damping [1, 2, 3, 7, 8]. To satisfy this criterion,
PSZS must be calculated by adopting a two-step averaging procedure. More
precisely: first, an average along guiding center equilibrium orbits is
applied; second, a filter is used, on the axisymmetric particle response, to
remove the fast spatiotemporal variations on the characteristic particle orbit
length-scale and/or the hydrodynamic time-scale. Consequently, PSZS depend
solely on the equilibrium invariants of motion, such as the particle energy
(per unit mass) $\mathcal{E}_{0}=v^{2}/2$222For simplicity of the present
analysis, we assume that equilibrium radial electric field, if it exists,
corresponds to sufficiently slow $\bm{E}\times\bm{B}$ flow that is consistent
with the gyrokinetic ordering and, thus, can be incorporated within the
perturbed radial electric field. If needed, this assumption could be readily
dropped., the magnetic moment $\mu=v_{\perp}^{2}/2B_{0}+\ldots$ and the
toroidal angular momentum $P_{\phi}$ [4]. It is worth noting that any other
combination of three invariants of motion can be used, for example involving
the ‘second invariant’ $J=m\oint v_{\parallel}dl=J(P_{\phi},{\cal E},\mu)$
[1]. In this section, we apply this approach to derive the governing equation
for PSZS in conservative form. Furthermore, we derive the equations providing
the deviation of the axisymmetric particle response from the PSZS and its
dynamic evolution. Then, the notion of PSZS is used to introduce the concept
of ZS, which, together with the ZFs defined below in Sec. 2.1, provides a
proper definition of plasma nonlinear equilibrium [7, 8, 26] that evolves
consistently with the (toroidal) symmetry breaking fluctuation spectrum as
well as with sources and collisions.

###  2.1 Orbit averaged particle response: PSZS and zonal state

Since PSZS depend only on the equilibrium constants of motion, their evolution
equation can be readily cast using these as coordinates in the phase-space.
Proceeding along these lines, we write the phase-space velocity appearing in
the gyrokinetic equation [4, 27, 5] as the sum of two contributions, i.e.
$\dot{\bm{Z}}=\dot{\bm{Z}}_{0}+\delta\dot{\bm{Z}}$, representing,
respectively, the integrable particle motion in the reference magnetic field
and the remaining particle response that we generically attribute to the
effect of fluctuations. This decomposition is general and could be applied to
any nearly integrable (Hamiltonian) system. It assumes that the reference
equilibrium, defined by the reference magnetic field and by the plasma
profiles that are consistent with it and with the PSZS, varies on a
sufficiently slow time scale. A more rigorous definition of reference or
“equilibrium” magnetic field, is given in Sec. 2.3. The self-consistency of
this description and approach can be rigorously checked a posteriori.
Consequently, the gyrokinetic equation in conservative form reads:

$\frac{\partial}{\partial
t}(DF)+\frac{\partial}{\partial\mathbf{Z}}\cdot\left(D\bm{\dot{}}{\bm{Z}}_{0}F\right)+\frac{\partial}{\partial\mathbf{Z}}\cdot\left(D\delta\bm{\dot{}}{\bm{Z}}F\right)=0\,,$
(1)

where $D$ is the velocity space Jacobian and $F$ the gyro-center distribution
function [27, 5]. Here, for the sake of simplicity, we have temporarily
suppressed collisions and source terms. We introduce
$(\theta,\zeta,P_{\phi},\mathcal{E}_{0},\mu)$ as phase space coordinates333For
circulating particles, the sign of particle motion parallel or anti-parallel
to the equilibrium magnetic field must be specified as well, and will be
implicitly assumed. where $\theta$ and $\zeta$ are, respectively, poloidal and
toroidal magnetic flux coordinates. Our focus now turns to the zonal
distribution function that is the toroidally symmetric part of $F$. This can
be obtained by extracting the $n=0$ component of its Fourier expansion, where
$n$ is the toroidal mode number. For symmetry reasons, this is the obvious
starting point for the definition of an equilibrium distribution function in
axisymmetric Tokamak plasmas. Without loss of generality, we assume an
axisymmetric equilibrium magnetic field, i.e.,
$\bm{B}_{0}=\hat{F}\bm{\nabla}\phi+\bm{\nabla}\phi\times\nabla\psi$ where
$\hat{F}=RB_{\phi}$ and $\phi$ is the toroidal angle, which is connected to
$\zeta$ as $\zeta=\phi-\nu(\psi,\theta)$, with $\nu(\psi,\theta)$ chosen such
that magnetic flux coordinates are characterized by straight magnetic field
lines. We now note that, in the zonal component of Eq. (1), we can re-write
the term describing the equilibrium motion as:

$\frac{\partial}{\partial\bm{Z}}\cdot\left(D\bm{\dot{}}{\bm{Z}}_{0}F\right)_{z}=\nabla\cdot\left(D\bm{\dot{}}{\bm{X}}_{0}F\right)_{z}=\frac{1}{{\cal
J}_{P_{\phi}}}\frac{\partial}{\partial\theta}(D{\cal
J}_{P_{\phi}}F\bm{\dot{}}{\bm{X}}_{0}\cdot\bm{\nabla}\theta)_{z}\,,$ (2)

where ${\cal J}_{P_{\phi}}={\cal J}(\partial P_{\phi}/\partial\psi)^{-1}$,
${\cal J}=(\bm{\nabla}\zeta\cdot\bm{\nabla}\psi\times\bm{\nabla}\theta)^{-1}$
is the Jacobian in flux coordinates, $\psi$ is the poloidal magnetic flux, and
the toroidal symmetry of the reference state has been used along with and the
conservation of $P_{\phi}$ and energy characterizing particle motion in the
equilibrium magnetic field, i.e., respectively
$\dot{\bm{X}}_{0}\cdot\bm{\nabla}P_{\phi}=0$ and $\mathcal{\dot{E}}_{0}=0$. We
can now orbit average the zonal component of Eq. (1) in the reference
equilibrium (slowly evolving in time); that is, an averaging operator along
$\theta$ on the gyrokinetic equation while using ${\cal J}_{P_{\phi}}$ as
weight annihilating, as expected, the term described in Eq. (2). Assuming that
the reference magnetic equilibrium is slowly evolving; e.g., on the resistive
current diffusion time, we finally obtain:

$\partial_{t}\oint d\theta{\cal J}_{P_{\phi}}DF_{z}+\oint d\theta{\cal
J}_{P_{\phi}}\frac{\partial}{\partial\bm{Z}}\cdot(D\delta\dot{\bm{Z}}F)_{z}=0\,.$
(3)

Recalling the governing equation in the absence of fluctuations [28, 21];
i.e., $\dot{\psi}=-v_{\|}\partial_{\theta}\bar{\psi}/({\cal J}B_{\|}^{*})$,
where $\bar{\psi}=\psi-RB_{\phi}v_{\|}/\Omega\simeq-(c/e)P_{\phi}$444Here, by
the $\simeq$ notation, we mean that $\bar{\psi}$ contains the leading order
expression of $P_{\phi}$. This is not a limitation of the present approach,
which can be carried out at the desired order of accuracy consistent with the
adopted gyrokinetic formulation [27, 5].,
$\dot{\theta}=v_{\parallel}\partial_{\psi}\bar{\psi}/({\cal
J}B_{\parallel}^{*})$ and $D=B_{\parallel}^{*}/|v_{\parallel}|$, with
$B_{\|}^{*}\equiv\mathbf{B}^{*}\cdot\mathbf{b}$,
$\bm{b}\equiv\bm{B}_{0}/B_{0}$,
$\mathbf{B}^{*}\equiv\nabla\times\mathbf{A}^{*}$,
$(e/c)\mathbf{A}^{*}\equiv(e/c)\mathbf{A}_{0}+m\left(v_{\|}\mathbf{b}\right)$,
$\mathbf{B}_{0}\equiv\nabla\times\mathbf{A}_{0}$, we can recognize that the
averaging used to derive Eq. (3) is indeed a time averaging along the
integrable particle orbit; denoted as

$\overline{(\ldots)}^{\rm{(O)}}=\tau_{b}^{-1}\oint
d\theta(...)/\dot{\theta}\;,$ (4)

with $\tau_{b}=\oint d\theta/\dot{\theta}$. Thus, we obtain the following
(equilibrium) orbit averaged kinetic equation:

$\frac{\partial}{\partial
t}\overline{F_{z}}^{\rm{(O)}}+\frac{1}{\tau_{b}}\left[\frac{\partial}{\partial
P_{\phi}}\overline{\left(\tau_{b}\delta\dot{P}_{\phi}F\right)_{z}}^{\rm{(O)}}+\frac{\partial}{\partial\mathcal{E}}\overline{\left(\tau_{b}\delta\dot{\mathcal{E}}F\right)_{z}}^{\rm{(O)}}\right]=\overline{\left(\sum_{s}C_{s}^{g}\left[F_{,}F_{s}\right]+\mathcal{S}\right)_{z}}^{\rm{(O)}}\,,$
(5)

where $\delta\dot{P}_{\phi}=\delta\dot{\mathbf{X}}\cdot\bm{\nabla}P_{\phi}$,
$\delta\dot{\mathcal{E}}$ is defined in Eq.(14), and we have restored
collisions and source terms on the RHS. The expressions of gyrokinetic
collision operators of the considered test particles with the field particle
species $s$, as denoted by the subscript in $C_{s}^{g}$, are given in Refs.
[27, 5]. Meanwhile, for the sake of notation simplicity, the summation over
different field particle species in the collisions term will be omitted from
now on. Denoting the spatial-temporal slowly evolving component of
$\overline{F_{z}}^{\rm{(O)}}$, i.e. PSZS, as
$\overline{F_{0}}^{\rm{(O)}}\equiv[\overline{F_{z}}^{\rm{(O)}}]_{S}$, the
relevant evolution equation is obtained by additionally extracting the
macro-/meso- scopic component of Eq. (5); i.e.:

$\frac{\partial}{\partial
t}\overline{F_{0}}^{\rm{(O)}}+\frac{1}{\tau_{b}}\left[\frac{\partial}{\partial
P_{\phi}}\overline{\left(\tau_{b}\delta\dot{P}_{\phi}\delta
F\right)_{z}}^{\rm{(O)}}+\frac{\partial}{\partial\mathcal{E}}\overline{\left(\tau_{b}\delta\dot{\mathcal{E}}\delta
F\right)_{z}}^{\rm{(O)}}\right]_{S}=\overline{C^{g}}^{\rm{(O)}}_{S}+\overline{S}^{\rm{(O)}}_{S}\,,$
(6)

where $[\ldots]_{S}$ denotes an appropriate (ad hoc) spatio-temporal averaging
procedure to be illustrated; and where we have postulated a bi-linear
collision term such that:

$\overline{C^{g}}^{\rm{(O)}}_{S}=\overline{C^{g}\left[\overline{F_{0}}^{\rm{(O)}},\overline{F_{0}}^{\rm{(O)}}\right]}^{\rm{(O)}}+\left[\overline{C^{g}\left[\overline{F_{0}}^{\rm{(O)}},\delta
F\right]}^{\rm{(O)}}+\overline{C^{g}\left[\delta
F,\overline{F_{0}}^{\rm{(O)}}\right]}^{\rm{(O)}}+\overline{C^{g}\left[\delta
F,\delta F\right]}^{\rm{(O)}}\right]_{S}.$ (7)

In the previous expression we have introduced the following decomposition at
each instant of time:

$F=\overline{F_{0}}^{\rm{(O)}}+\delta F\,.$ (8)

The aforementioned spatio-temporal averaging over the micro-scales is what
allows us to separate $\overline{F_{0}}^{\rm{(O)}}$ from
$\overline{F_{z}}^{\rm{(O)}}$, given by Eq. (5). It is arbitrary, to some
extent, and can be specified for convenience according to the problem of
interest. Our choice, here, is to write Eq. (6) as definition of
$\overline{F_{0}}^{\rm{(O)}}$, and based on Eq. (5), include all residual
response in $\overline{\delta
F_{z}}^{\rm{(O)}}\equiv\overline{F_{z}}^{\rm{(O)}}-\overline{F_{0}}^{\rm{(O)}}$
consistent with Eqs. (9) and (10) below. This point will be further discussed
in Sec. 2.2. The approach proposed in the present analysis could be considered
a “full-F” description [29, 30] of the nonlinear evolving equilibrium, and a
“delta-F” approach [31] for the (toroidal) symmetry breaking perturbations
(cf. [32] for a general review for gyrokinetic simulations of turbulent
transport).

Here, it is also worthwhile to briefly remark that the ratio between the third
and the second terms of LHS in Eq. (5) can be estimated as
$\delta\,\mathcal{E}/\Delta\mathcal{E}$ with $\delta P_{\phi}/\Delta
P_{\phi}\sim{\cal O}(1)$, where $\Delta\mathcal{E}$ and $\Delta P_{\phi}$ are,
respectively, PSZS characteristic scales in energy and toroidal angular
momentum space; and $\delta\mathcal{E}$ and $\delta P_{\phi}$ are the
corresponding nonlinear distortions due to the considered fluctuation
spectrum. Using the typical nonlinear gyrokinetic ordering, this is consistent
with the relatively small effect of the so-called parallel nonlinearity on
fluctuation induced phase-space transport. This is not longer the case for Eq.
(6), where the two terms are generally of the same order. However, once the
effect of the third term is integrated over in energy space, the corresponding
overall effect can, again, be shown to be negligible at the leading order.
Consistently, in Ref. [7] a PSZS transport theory has been formulated omitting
the parallel non-linearity term and adopting the classical Frieman-Chen
formulation of the nonlinear gyrokinetic equation, which is appropriate up to
leading order in the multi-spatiotemporal scale asymptotic expansion [4]. In
the present work, the parallel nonlinearity is retained; consistent with the
conservative formulation of nonlinear gyrokinetics [27, 5].

Having introduced the concept of PSZS, we can decompose the whole gyrocenter
particle response and, consequently, the zonal component of the gyrokinetic
distribution $F_{z}$, into the sum of different terms accounting for the
various relevant physics processes, i.e.:

$F_{z}=\overline{F_{z}}^{\rm{(O)}}+\delta\tilde{F}_{z}=\overline{F_{0}}^{\rm{(O)}}+\overline{\delta
F_{z}}^{\rm{(O)}}+\delta\tilde{F}_{z}.$ (9)

In particular, as already stated, the PSZS, $\overline{F_{0}}^{\rm{(O)}}$,
describes the evolution of the macro-/meso-scopic equilibrium. The phase-space
transport theory, derived in this work is primarily motivated by the fact that
this contribution may increasingly deviate in time from the reference
thermodynamic equilibrium due to nonlinear processes; eventually invalidating
the usual transport analyses that rely, e.g., on local Maxwellian equilibria.
Notable applications are analyses of EP transport in fusion plasmas [1, 2, 3],
but deviation of (bulk) particle equilibria from local Maxwellian can be also
important to explain, e.g., the nonlinear up-shift (the so-called “Dimits
shift” [33]) of the threshold for ion temperature gradient driven turbulence
[26]. In the following, we will show that a multipole expansion can be applied
to the PSZS fluid moments [5, 21] yielding an anisotropic CGL pressure tensor
[20, 34] and a self-consistently evolving nonlinear magnetic equilibrium.
Further to $\overline{F_{0}}^{\rm{(O)}}$, the residual components of $F_{z}$
either describe the micro-scale spatio-temporal variation of the orbit
averaged distribution function or have zero average along equilibrium orbits.
More precisely, the residual orbit averaged particle response,
$\overline{\delta F_{z}}^{\rm{(O)}}$, characterizes the transition between
neighboring nonlinear equilibria, which are all undamped by collisionless
dissipation [7, 26] and slightly deviate from the reference
$\overline{F_{0}}^{\rm{(O)}}$ as schematically described in Fig. 1.

Figure 1: Schematic picture describing the equilibrium orbit averaged
distribution function $\overline{F_{z}}^{\rm{(O)}}$. The solid line represents
the slowly varying component of the orbit averaged distribution function,
while the oscillation around it corresponds to $\overline{\delta
F_{z}}^{\rm{(O)}}$.

These neighboring nonlinear equilibria [26] should be understood as the
ensemble of different realizations of the system, thereby representing the
connection between time average, introduced above in the definition of
$\overline{F_{0}}^{\rm{(O)}}$ by means of Eq. (6), and “ensemble average” in a
statistical sense [7]. The distribution function
$\overline{F_{z}}^{\rm{(O)}}{=\overline{F_{0}}^{\rm{(O)}}+\overline{\delta
F_{z}}^{\rm{(O)}}}$, together with the undamped components of the
electromagnetic potentials (ZFs), constitute the ZS, i.e., “the collisionless
undamped (long-lived) nonlinear deviation of the plasma from the reference
state as a consequence of fluctuation-induced transport processes, due to
emission and reabsorption of (toroidal equilibrium) symmetry-breaking
perturbations” [7]. We note that the symmetry breaking fluctuations (with
$n\neq 0$) are not explicitly mentioned as elements of the ZS, but they are
self-consistently accounted for in its definition, as they have determined
$\overline{F_{z}}^{\rm{(O)}}$ during the dynamic evolution of the system. We
also note that, although PSZS can be obtained separating out the macro-/meso-
scopic component of $\overline{F_{z}}^{\rm{(O)}}$ consistently with the usual
definition of equilibrium and transport, this separation is somehow arbitrary
and depends on the specific problem of interest. For example, when describing
EPM [35],
$\overline{F_{z}}^{\rm{(O)}}=\overline{F_{0}}^{\rm{(O)}}+\overline{\delta
F_{z}}^{\rm{(O)}}$ is best considered as a whole since in this case phase-
space transport occurs on the same time scale of the nonlinear dynamic
evolution of the spectrum of fluctuations [1, 2]. Following the previous
derivation, we obtain the governing equation for the spatio-temporal micro-
scale component of the equilibrium orbit averaged distribution function:

$\frac{\partial}{\partial t}\overline{\delta
F_{z}}^{\rm{(O)}}+\frac{1}{\tau_{b}}\left[\frac{\partial}{\partial
P_{\phi}}\overline{\left(\tau_{b}\delta\dot{P}_{\phi}\overline{F_{0}}^{\rm{(O)}}\right)_{z}}^{\rm{(O)}}+\frac{\partial}{\partial\mathcal{E}}\overline{\left(\tau_{b}\delta\dot{\mathcal{E}}\overline{F_{0}}^{\rm{(O)}}\right)_{z}}^{\rm{(O)}}\right]\\\
+\frac{1}{\tau_{b}}\left[\frac{\partial}{\partial
P_{\phi}}\overline{\left(\tau_{b}\delta\dot{P}_{\phi}\delta
F\right)_{z}}^{\rm{(O)}}+\frac{\partial}{\partial\mathcal{E}}\overline{\left(\tau_{b}\delta\dot{\mathcal{E}}\delta
F\right)_{z}}^{\rm{(O)}}\right]_{F}=[\overline{C^{g}}^{\rm{(O)}}_{z}+\overline{S}^{\rm{(O)}}]_{F}.$
(10)

Note that, here, $[\ldots]_{F}$ denotes the spatio-temporal micro-scale
component of the argument such that $[\ldots]\equiv[\ldots]_{S}+[\ldots]_{F}$;
that is, the spatial variation on Larmor radius and finite magnetic drift
orbit width length scale, and the temporal variation on the hydrodynamic time
scale. It should be pointed out that, as opposed to the governing equation for
PSZS, there is a formally linear term in the orbit averaged response on the
LHS. This term may have fast as well as slow spatio-temporal variations and,
thus, the subscript $F$ is omitted there. Furthermore, this same term is also
responsible for the high frequency oscillation characterizing the geodesic
acoustic mode (GAM) [36] and, therefore, cannot be included into the
definition of a macroscopic equilibrium consistent with usual transport time
scale orderings [7]. Interested readers can find linear as well as nonlinear
GAM/EGAM physics discussed in detail in A.

###  2.2 Nonlinear equilibrium as renormalized particle response

In the previous subsection, we showed that the concept of PSZS is
intrinsically related to the integrable “equilibrium” guiding center motion
and, thus, it is naturally described using $P_{\phi}$ as phase-space
coordinate. In particular, the governing equations for the different
components of the zonal distribution function take very compact expressions.
However, when describing the self-consistent evolution of ZFs, we must adopt
standard flux coordinates $(\psi,\theta,\zeta)$. We can define the associated
change of coordinates between these two representations noting that
$P_{\phi}=P_{\phi}(\psi,\theta,\mathcal{E},\mu)=-(e/c)(\psi-\delta\tilde{\psi}(\psi,\theta,\mathcal{E},\mu))$.
Thus, consistently with previous works [37, 1, 7], one can introduce a shift
operator, formally represented as $e^{iQ}$, accounting for the guiding center
equilibrium magnetic drift motion, which therefore provides the push-forward
transformation from gyrocenter to magnetic drift/oscillating-centers. Then,
the (equilibrium) orbit average of a scalar function
$\hat{H}(P_{\phi},\theta)=H(\bar{\psi}+\delta\tilde{\psi},\theta)$ reads:

$\oint\frac{d\theta}{\dot{\theta}}\hat{H}(P_{\phi},\theta)=\oint\frac{d\theta}{\dot{\theta}}H(\bar{\psi}+\delta\tilde{\psi},\theta)=\oint\frac{d\theta}{\dot{\theta}}e^{iQ}H(\bar{\psi},\theta)$
(11)

where $\bar{\psi}=-(c/e)P_{\phi}$ was defined in Sec. 2.1 above and
$\delta\tilde{\psi}=RB_{\phi}v_{\parallel}/\Omega$ at the leading order, with
its dependence on $\theta$ and the other phase-space variables left implicit.
It follows by direct inspection that, as expected, (equilibrium) orbit
averaging is equivalent to a bounce/transit average combined with the action
of the shift operator $e^{iQ}$. As a further remark, we recall that
bounce/transit averaging is also connected with flux surface averaging of
velocity space integrals and, consequently, with the standard representation
of plasma (radial) transport equations. For this reason, the PSZS governing
equation is particularly relevant for describing plasma transport and allows
recovering well known results [10], and further generalizing them [7]. In
order to see the equivalence between orbit averaging and “shifted bounce
averaging” more clearly, let us define $\overline{(\ldots)}=\tau_{b}^{-1}\oint
d\theta(...)/\dot{\theta}$, where, now, the closed poloidal orbit integral
follows the constant-$\theta$ projection of the actual guiding center orbit on
the $\bar{\psi}$ flux surface. This definition of $\tau_{b}$ with respect to
the orbit averaging approach, is unambiguous since it is uniquely defined
being $\theta$ a dummy integration variable. Then,

$\left.\overline{(\ldots)}^{\rm{(O)}}\right|_{P_{\phi}}=\left.\overline{e^{iQ}(\ldots)}\right|_{\bar{\psi}}\,,$
(12)

where, for further clarity, we have explicitly denoted by the additional
subscripts the reference value of $P_{\phi}$ on the LHS, and of $\bar{\psi}$
on the RHS. Rephrasing this concept, Eq. (12) states that orbit average for
given $P_{\phi}$ (implicitly assuming given ${\cal E}$ and $\mu$) corresponds
to a proper shifted bounce averaging on the flux surface labeled by
$\bar{\psi}$.

Due to the equivalence between orbit and bounce/transit averaging, the
governing equation for PSZS introduced in the previous subsection are
consistent with those obtained in Refs. [1, 7], with the further inclusion of
the effects of the so called parallel nonlinearity. In what follows, we re-
derive the governing equations for the different components of $F_{z}$ using
$(\psi,\theta,\zeta,\mathcal{E},\mu)$ as phase-space coordinates for two main
reasons: first, in order to establish a close contact with the representation
that is most conveniently adopted for writing equations describing mode
structures, either ZFs or symmetry breaking perturbations; and second, for
demonstrating that the PSZS definition adopted in this and earlier works [1,
2, 3, 7] corresponds to a renormalization of equilibrium particle response.

As a first step towards the second goal, we note that the decomposition of Eq.
(9) can be written by introducing the drift/banana center pull-back operator
$e^{-iQ_{z}}$ where
$Q_{z}=RB_{\phi}\left(v_{\|}/\Omega\right)k_{z}/(d\psi/dr)$ [7, 8, 26, 37].
This is the explicit expression for the shift operator introduced previously
and the subscript $z$ to the radial wave number $k_{z}\equiv-i\partial_{r}$
reminds that it is acting on toroidally symmetric response. As noted above,
$Q_{z}$ is the leading order expression of $Q$; and more accurate expressions
for $Q_{z}$ could be given by means of corresponding more accurate expressions
of $P_{\phi}$, consistent with the adopted gyrokinetic description [27, 5]. In
particular, $F_{z}$ can be written as:

$F_{z}\equiv\overline{F}_{0}+e^{-iQ_{z}}\left(\overline{\delta
F_{Bz}}\Big{|}_{F}+\delta\tilde{F}_{Bz}\right)=\overline{F}_{0}+e^{-iQ_{z}}\left(\left.\overline{e^{iQ_{z}}\delta
F_{z}}\right|_{F}+\delta\tilde{F}_{Bz}\right)$ (13)

where the function $\delta F_{Bz}$ is the drift/banana center particle
response, the bar stands for bounce/transit averaging, the tilde denotes the
vanishing bounce/transit average response and $\bar{F}_{0}$ is, from now on, a
short notation for $\overline{F}^{\rm{(O)}}_{0}$ and, thus, describes the PSZS
component. Note the one-to-one correspondence of Eq. (13) and Eq. (9), which
also illuminates the notation. Following the previous subsection, we now
proceed in deriving the governing equations for the different terms of this
decomposition. In particular, we recall the gyrokinetic expression for the
time variation of the energy per unit mass, i.e.:

$\delta\dot{\mathcal{E}}=-\frac{e}{m}\dot{\bm{X}}_{0}\cdot\nabla\left\langle\delta
L_{g}\right\rangle_{z}\,,$ (14)

where angular brackets denote gyro-phase averaging, $\delta
L=\delta\phi-\bm{v}\cdot\delta\bm{A}/c$, $\delta\phi$ is the scalar potential,
$\delta\bm{A}$ is the vector potential with $\delta
L_{g}=e^{\bm{\rho}\cdot\bm{\nabla}}\delta L(\bm{X})=\delta
L(\bm{X}+\bm{\rho})$, $\bm{\rho}=\bm{b}_{0}\times\bm{v}/\Omega$ and
$\Omega=eB_{0}/(mc)$. Note that, $\left\langle\delta
L_{g}\right\rangle_{z}=J_{0}(\lambda)(\delta\phi_{z}-v_{\parallel}\delta
A_{\parallel z})+(2/\lambda)(m/e)\mu J_{1}(\lambda)\delta B_{\parallel z}$,
$\lambda^{2}=2(\mu B_{0}/\Omega^{2})k_{\perp}^{2}$ and $J_{0,1}$ are Bessel
functions. We also recall the conservation of the toroidal component of the
canonical angolar momentum in the presence of toroidally symmetric
perturbations:

$\delta\dot{\theta}\frac{\partial
P_{\phi}}{\partial\theta}+\delta\dot{\psi}\frac{\partial
P_{\phi}}{\partial\psi}+\delta\dot{\mathcal{E}}\frac{\partial
P_{\phi}}{\partial\mathcal{E}}=-\frac{e}{c}\left(\partial_{t}+\dot{\bm{X}}_{0}\cdot\nabla\right)\frac{RB_{\phi}\langle\delta
A_{\parallel g}\rangle_{z}}{B_{0}}\;.$ (15)

Thus, we can re-write the toroidally symmetric component of Eq. (1) as
follows:

$\displaystyle D($
$\displaystyle\partial_{t}+\dot{\bm{X}}_{0}\cdot\bm{\nabla})\left(F_{z}-\frac{e}{m}\langle\delta
L_{g}\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\Big{|}_{\bar{\psi}}+\frac{RB_{\phi}}{B_{0}}\frac{\partial\bar{F}_{0}}{\partial\bar{\psi}}\langle\delta
A_{\parallel g}\rangle_{z}\right)+$
$\displaystyle{-D\frac{RB_{\phi}}{B_{0}}\langle\delta A_{\parallel
g}\rangle_{z}\frac{\partial}{\partial\bar{\psi}}\partial_{t}\bar{F}_{0}}+D\frac{e}{m}{\frac{\partial}{\partial
t}\left(\left.\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\right|_{\bar{\psi}}\langle\delta
L_{g}\rangle_{z}\right)}+$ $\displaystyle+\frac{1}{{\cal
J}}\frac{\partial}{\partial\theta}({\cal J}D\delta\dot{\theta}\delta
F)+\frac{1}{{\cal J}}\frac{\partial}{\partial\psi}({\cal
J}D\delta\dot{\psi}\delta
F)+\frac{\partial}{\partial\mathcal{E}}(D\delta\dot{\mathcal{E}}\delta
F)=D(C^{g}+\mathcal{S}).$ (16)

This equation, consistently with Ref. [8, 4], suggests to introduce the
following definition:

$G_{z}\equiv F_{z}-\left.\frac{e}{m}\left\langle\delta
L_{g}\right\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\right|_{\bar{\psi}}+\frac{RB_{\phi}}{B_{0}}\left\langle\delta
A_{\|g}\right\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\bar{\psi}}$ (17)

where, as radial coordinate, we are using $\bar{\psi}\equiv-(c/e)P_{\phi}$
introduced earlier. From the definition above, the role of $\bar{F}_{0}$ as
renormalized reference distribution function taking into account nonlinear
plasma behaviors (self-interactions) consistently with the theoretical
framework introduced in [1, 7] is made clear. In fact, consistently with Eqs.
(9) and Eq. (13), no distinction is made in Eqs. (2.2) and (17) between the
$\overline{\delta F_{z}}^{\rm{(O)}}=\overline{e^{iQ_{z}}\delta F_{z}}$
contribution that should be kept distinct from $\bar{F}_{0}$ and the one that
can be reabsorbed in it. Thus, the distinction can be made for convenience of
identification of a reference magnetic equilibrium involving macro- and meso-
scale kinetic profiles only (cf. next subsection); but, the physics analysis
of phase-space structures that are undamped by linear collisionless processes
is “full-$F$” by construction. We may also note that, consistently with Eqs.
(9) and Eq. (13), the reference state appearing in Eq. (17) generally includes
a spatio-temporal micro-scale contribution. However, consistently with the
gyrokinetic ordering [4, 5], this term can be neglected at the relevant
leading order. We can now write $G_{z}$ in terms of the drift/banana shift
operator, i.e. $G_{z}=e^{-iQ_{z}}G_{Bz}$, substitute this expression in Eq.
(2.2) and apply $e^{iQ_{z}}$ on both sides. We find:

$e^{iQ_{z}}D\dot{\bm{X}}_{0}\cdot\bm{\nabla}G_{z}=e^{iQ_{z}}\frac{v_{\parallel}}{{\cal
J}|v_{\parallel}|}\left[1-\frac{\partial}{\partial\psi}\frac{RB_{\phi}v_{\parallel}}{\Omega}\right]e^{-iQ_{z}}\frac{\partial
G_{Bz}}{\partial\theta}=\frac{v_{\parallel}}{{\cal
J_{\bar{\psi}}}|v_{\parallel}|}\frac{\partial}{\partial\theta}G_{Bz}\,,$ (18)

where ${\cal J}_{\bar{\psi}}=-(e/c){\cal J_{P_{\phi}}}$ is computed at the
actual gyrocenter particle position; and we have noted
$iQ_{z}=(\psi-\bar{\psi})\partial_{\psi}$. Considering the effect of the shift
operator on $D$, we obtain the following kinetic equation:

${{\cal
J}_{\bar{\psi}}D}\partial_{t}(G_{Bz})+\frac{v_{\parallel}}{|v_{\parallel}|}\partial_{\theta}G_{Bz}=\\\
e^{iQ_{z}}\Big{[}-\frac{e}{m}{\cal J}_{\bar{\psi}}D\frac{\partial}{\partial
t}\left(\langle\delta
L_{g}\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\Big{|}_{\bar{\psi}}\right)+{\cal
J}_{\bar{\psi}}D\frac{RB_{\phi}}{B_{0}}\langle\delta A_{\parallel
g}\rangle_{z}\frac{\partial}{\partial\bar{\psi}}\partial_{t}\bar{F}_{0}\Big{]}\\\
-\frac{\partial}{\partial\theta}(e^{iQ_{z}}{\cal
J}_{\bar{\psi}}D\delta\dot{\theta}\delta
F)-\frac{\partial}{\partial\psi}(e^{iQ_{z}}{\cal
J}_{\bar{\psi}}D\delta\dot{\bar{\psi}}\delta
F)-\frac{\partial}{\partial\mathcal{E}}(e^{iQ_{z}}{\cal
J}_{\bar{\psi}}D\delta\dot{\mathcal{E}}\delta F)+e^{iQ_{z}}{\cal
J}_{\bar{\psi}}D(C_{g}+\mathcal{S}).$ (19)

In fact, it can be shown that the shift operator commutes with the partial
derivatives of nonlinear terms in Eq. (2.2) taking into account formal
simplifications among commutators. Note that $\delta\dot{\bar{\psi}}$ is
typically dominated by $\delta\dot{\psi}$ for $n\neq 0$ symmetry breaking
perturbations. The nonlinearities caused by ZFs, meanwhile, need a special
attention, since $\delta\dot{\psi}_{z}$ vanishes at the leading order. We will
come back to this technical but important point while describing some
applications of this theory in A. Integrating over $\theta$ on a closed
trajectory, the $\theta$ derivatives can be annihilated. Recalling the bounce
average definition introduced previously, the following expression is finally
obtained:

$\displaystyle\partial_{t}\overline{G_{Bz}}=$
$\displaystyle-\overline{e^{iQ_{z}}\frac{e}{m}\partial_{t}\left[\left.\left\langle\delta
L_{g}\right\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\right|_{\bar{\psi}}\right]}+\overline{e^{iQ_{z}}\frac{RB_{\phi}}{B_{0}}\left\langle\delta
A_{\|g}\right\rangle_{z}\frac{\partial}{\partial\bar{\psi}}\partial_{t}\bar{F}_{0}}+\left.\overline{e^{iQ_{z}}\left[C_{g}+\mathcal{S}\right]}\right|_{z}$
(20)
$\displaystyle-\frac{1}{\tau_{b}}\frac{\partial}{\partial\psi}\left[\tau_{b}\overline{e^{iQ_{z}}\delta\dot{\bar{\psi}}\delta
F}\right]_{z}-\frac{1}{\tau_{b}}\frac{\partial}{\partial\mathcal{E}}\left[\tau_{b}\overline{e^{iQ_{z}}\delta\dot{\mathcal{E}}\delta
F}\right]_{z}$

This is a generalization of Eq. (25) in Ref. [7] written in conservative form
and retaining the role of parallel nonlinearity, collisions and source terms.
Consequently, recalling the relationship between bounce/transit and flux
surface averaging, from this expression it is possible to derive all the usual
flux surface averaged transport equations [7]. The governing equation for
$\bar{F}_{0}$ follows directly from Eq. (20) and is consistent with Eq. (6):

$\partial_{t}\overline{e^{iQ_{z}}\bar{F}_{0}}=-\left.\overline{e^{iQ_{z}}\frac{RB_{\phi}}{B_{0}}\partial_{t}\left\langle\delta
A_{\parallel
g}\right\rangle_{z}\frac{\partial}{\partial\bar{\psi}}\bar{F}_{0}}\right|_{S}+\left.\overline{e^{iQ_{z}}\left[C_{g}+{\cal
S}\right]}\right|_{zS}+\\\
-\frac{1}{\tau_{b}}\frac{\partial}{\partial\psi}\left[\tau_{b}\overline{e^{iQ_{z}}\delta\dot{\bar{\psi}}\delta
F}\right]_{zS}-\frac{1}{\tau_{b}}\frac{\partial}{\partial{\cal
E}}\left[\tau_{b}\overline{e^{iQ_{z}}\delta\dot{\cal E}\delta
F}\right]_{zS}\;.$ (21)

Here, it is worthwhile noting that, except for the first (ZFs inductive) term
on the RHS, Eq. (21) shows that PSZS evolution is either caused by nonlinear
interactions or by sources/collisions. Physically, the ZFs inductive term is
due to the externally or nonlinearly generated perturbation of the magnetic
flux function. Thus, the present definition of PSZS, consistent with earlier
works [1, 2, 3, 7], describes the renormalized reference distribution function
taking into account nonlinear plasma behaviors (self-interactions). The orbit
averaged fast spatiotemporal deviation of the plasma response about the PSZS
is given by:

$\displaystyle\partial_{t}\left.\overline{\delta g_{Bz}}\right|_{F}$
$\displaystyle=-\left.\overline{e^{iQ_{z}}\frac{e}{m}\partial_{t}\left[\left\langle\delta
L_{g}\right\rangle_{z}\left.\frac{\partial}{\partial{\cal
E}}\right|_{\bar{\psi}}\bar{F}_{0}\right]}\right|_{F}+\left.\overline{e^{iQ_{z}}\frac{RB_{\phi}}{B_{0}}\left\langle\delta
A_{\parallel
g}\right\rangle_{z}\frac{\partial}{\partial\bar{\psi}}\partial_{t}\bar{F}_{0}}\right|_{F}$
$\displaystyle+\left.\overline{e^{iQ_{z}}\left[C_{g}+{\cal
S}\right]}\right|_{zF}-\frac{1}{\tau_{b}}\frac{\partial}{\partial\psi}\left[\tau_{b}\overline{e^{iQ_{z}}\delta\dot{\bar{\psi}}\delta
F}\right]_{zF}-\frac{1}{\tau_{b}}\frac{\partial}{\partial{\cal
E}}\left[\tau_{b}\overline{e^{iQ_{z}}\delta\dot{\cal E}\delta
F}\right]_{zF}\;.$ (22)

where $\delta g_{z}=e^{-iQ_{z}}\delta g_{Bz}$, consistent with Eq. (17), is
the nonadiabatic particle response that is connected with $\delta F_{z}$ by:

$\delta g_{z}\equiv\delta F_{z}-\left.\frac{e}{m}\left\langle\delta
L_{g}\right\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\right|_{\bar{\psi}}+\frac{RB_{\phi}}{B_{0}}\left\langle\delta
A_{\|g}\right\rangle_{z}\frac{\partial\bar{F}_{0}}{\partial\bar{\psi}}\;.$
(23)

Similarly, after some lengthy but straightforward algebra, one can obtain the
governing equation for $\delta\tilde{g}_{Bz}=\delta g_{Bz}-\overline{\delta
g_{Bz}}$ [8]. Equation (21), or the equivalent Eq. (6), and Eq. (2.2)
completely describe the ZS, introduced and defined in the previous subsection,
once the reference magnetic equilibrium and the evolution equations for the
ZFs are given along with the symmetry breaking fluctuation spectrum. This is
done in the next subsection. In particular, the ZS, consisting of neighboring
nonlinear equilibria [26] which, can be thought of as ensemble of different
realizations of the system [7], can be written as [8]

$F_{0*}\equiv\bar{F}_{0}+e^{-iQ_{z}}\left.\overline{e^{iQ_{z}}\delta
F_{z}}\right|_{F}.$ (24)

###  2.3 Chew Goldberger Low reference equilibrium and zonal fields

The motivation for analyzing phase-space features of transport processes in
low collisionality burning plasmas is that PSZS could significantly deviate
from a given model plasma equilibrium, e.g. Maxwellian, over long time scales
and, thus, the usual transport description as evolution of macroscopic radial
profiles may become inadequate [7]. Consistently with these motivations, it is
necessary to describe the modification of the reference magnetic equilibrium
self-consistently with PSZS. Following Refs. [2, 21, 5], we recall that the
transformation from the gyrocenter to the particle distribution function can
be cast as:

$\displaystyle f=$ $\displaystyle
e^{-\rho\cdot\nabla}\left[{F}-\frac{e}{m}\left(\frac{\partial{\overline{F_{0}}}}{\partial\mathcal{E}}+\frac{1}{B_{0}}\frac{\partial{\overline{F_{0}}}}{\partial\mu}\right)\left\langle\delta
L_{g}\right\rangle\right]+\frac{e}{m}\left[\frac{\partial{\overline{F_{0}}}}{\partial\mathcal{E}}\delta\phi+\frac{1}{B_{0}}\frac{\partial{\overline{F_{0}}}}{\partial\mu}\delta
L\right].$ (25)

Using this expression, we can write every fluid moment in terms of its so
called push-forward representation [5]. In the following, differently from the
usual approach, we will describe the guiding center transformation using the
Dirac delta formalism instead of the $e^{-\bm{\rho}\cdot\bm{\nabla}}$
operator, extending to phase-space the velocity space integrals on the
particle distribution function, expressed as in Eq. (25). As a simple example,
the toroidally symmetric plasma current density $\mathbf{J}_{z}$ reads:

$\displaystyle\mathbf{J}_{z}(\mathbf{r})=e\int d\mathcal{E}d\mu d\alpha
d^{3}\mathbf{X}\,D\left(T^{-1}_{gc}\mathbf{v}\right)\delta(\mathbf{X}+\bm{\rho}-\mathbf{r})\Big{[}\bar{F}_{0}+\delta
F_{z}-\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\langle\delta
L_{g}\rangle_{z}+$
$\displaystyle-\frac{e}{m}\frac{1}{B_{0}}\frac{\partial\bar{F}_{0}}{\partial\mu}\langle\delta
L_{g}\rangle_{z}\Big{]}+\frac{e^{2}}{m}\int d\mathcal{E}d\mu d\alpha
d^{3}\mathbf{X}\,\mathbf{v}D\Big{[}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}\delta\phi_{z}+\frac{1}{B_{0}}\frac{\partial\bar{F}_{0}}{\partial\mu}\delta
L_{z}\Big{]}\;,$ (26)

where $\alpha$ is the gyrophase, $T^{-1}_{gc}\mathbf{v}$ represents the
guiding-center transformation of the velocity $\mathbf{v}$ and the argument of
the delta function accounts for the relation between the particle position
$\mathbf{r}$ and the guiding center position $\mathbf{X}$ [21]. The pressure
tensor can be derived analogously. In the present approach, ZFs are considered
explicitly as a distortion of the nonlinear equilibrium, that is of the zonal
state. Thus, the reference magnetic equilibrium must be computed assuming only
the PZSZ as describing the reference state; i.e., only the
$\propto\bar{F}_{0}$ term in the push forward representation of the fluid
moments such as Eq. (26). Due to the macro-/meso-scopic nature of PSZS, and
applying the usual multipole expansion in the push-forward representation of
the fluid moments [5, 21], one can obtain a CGL pressure tensor and a
toroidally symmetric current satisfying the following force balance equation:

$\sigma\frac{\mathbf{J}_{z}\times\mathbf{B}_{0}}{c}=\nabla
P_{\|}+(\sigma-1)\nabla\left(\frac{B_{0}^{2}}{8\pi}\right)+\frac{B_{0}^{2}}{4\pi}\nabla_{\perp}\sigma$
(27)

where $\perp$ and $\parallel$ denote the components perpendicular and parallel
to $\bm{B}_{0}$ and

$\sigma=1+\frac{4\pi}{B_{0}^{2}}\left(P_{\perp}-P_{\|}\right).$ (28)

It is well-known that, assuming
$\bm{B}_{0}=\hat{F}\nabla\phi+\nabla\phi\times\nabla\psi$, the radial
component of this expression reads:

$\Delta^{*}\psi+\bm{\nabla}\ln\sigma\cdot\bm{\nabla}\psi=-\frac{4\pi
R^{2}}{\sigma}\frac{\partial
P_{\parallel}}{\partial\psi}-\frac{1}{\sigma^{2}}\frac{\partial
G}{\partial\psi}$ (29)

where $\Delta^{*}$ is the usual Grad-Shafranov operator and
$G(\psi)=(\sigma\hat{F})^{2}/2$ is a flux function. Meanwhile, pressure
components and $\hat{F}(\psi,B_{0})$ function are connected by the parallel

$\frac{\partial P_{\parallel}}{\partial
B_{0}}=\frac{P_{\parallel}-P_{\perp}}{B_{0}}$ (30)

and bi-normal

$\frac{\partial P_{\perp}}{\partial
B_{0}}=\frac{P_{\perp}-P_{\parallel}}{B_{0}}-\sigma\frac{B_{0}^{2}}{4\pi}\frac{\partial\ln\hat{F}}{\partial
B_{0}}$ (31)

components of Eq. (27). The resulting solution of Eqs. (27) to (31) defines
the magnetic equilibrium that is consistent with the presence of PSZS, which
$\bm{J}_{z}$ as well as $P_{\perp}$ and $P_{\parallel}$ have been computed
from. More precisely, $P_{\perp}$ and $P_{\parallel}$ can be calculated
integrating the PSZS distribution function and, then, $\hat{F}(\psi,B_{0})$ is
obtained from the expression of the poloidal plasma current, which is also
computed from $\bar{F}_{0}$, and from Eq.(31). Finally a standard Grad
Shafranov problem must be solved . This produces, as expected [21, 5, 38], an
anisotropic MHD equilibrium. It is worth noting that this result holds at the
leading order in the multipole expansion. At higher order, we could compute
the macro-/meso-scopic deviations from the CGL pressure tensor, which is
expected to become relevant when the multipole expansion does not hold; e.g.,
when steep gradient regions are encountered, viz., near the last closed
magnetic surface. More generally, whenever the length scale of the gradients
becomes comparable with the characteristic length of particle orbits, the
proposed separation of scales to isolate the PSZS is no longer applicable and
a “full-$F$” approach is mandatory [39].

The neighboring nonlinear equilibria [8, 26, 7]; that is, the micro
spatiotemporal deviation from the reference state given by PSZS and the just
constructed anisotropic (CGL) reference magnetic equilibrium, can also be
self-consistently determined along with the ZFs; i.e., $\delta\phi_{z}$,
$\delta A_{\parallel z}$ and $\delta B_{\parallel z}$ are obtained by means of
quasineutrality and Ampère equations following the well known theoretical
framework described in [2]:

$\displaystyle\sum_{s}\left\langle\frac{e^{2}}{m}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}\right\rangle_{v}\delta\phi_{z}+\nabla\cdot\sum_{s}\left\langle\frac{e^{2}}{m}\frac{2\mu}{\Omega^{2}}\frac{\partial\bar{F}_{0s}}{\partial\mu}\left(\frac{J_{0}^{2}-1}{\lambda^{2}}\right)\right\rangle_{v}\nabla_{\perp}\delta\phi_{z}$
$\displaystyle\hskip 20.00003pt+\sum_{s}\left\langle eJ_{0}(\lambda)\delta
g_{z}\right\rangle_{v}+\sum_{s}\left\langle
eJ_{0}(\lambda)\bar{F}_{0s}\right\rangle_{v}=0\;,$ (32)
$\displaystyle\frac{\partial}{\partial t}\delta
A_{\|z}=-\left[\frac{1}{B_{0}}\bm{b}_{0}\times\bm{\nabla}\delta
A_{\|}\cdot\bm{\nabla}\left(\nabla_{\parallel}^{-1}\partial_{t}\delta
A_{\parallel}\right)\right]_{z}\;,$ (33) $\displaystyle\nabla_{\perp}\delta
B_{\|z}=\bm{\kappa}_{0}\delta B_{\parallel
z}+\nabla_{\parallel}\delta\mathbf{B}_{\perp
z}+\nabla\mathbf{b}_{0}\cdot\delta\mathbf{B}_{\perp
z}+\frac{4\pi}{c}\delta\mathbf{J}_{\perp z}.$ (34)

Here, $\sum_{s}$ denotes summation on all particle species,
$\left\langle...\right\rangle_{v}$ stands for velocity space integration,
$\nabla_{\parallel}^{-1}$ is the inverse operator of $\nabla_{\parallel}$,
$\bm{\kappa}_{0}\equiv\bm{b}_{0}\cdot\bm{\nabla}\bm{b}_{0}$ is the magnetic
field curvature, $\delta\mathbf{J}_{\perp z}$ is readily obtained from Eq.
(26), $\delta\bm{B}_{\perp z}$ and $\delta B_{\parallel z}$ are expressed in
terms of the fluctuating vector potential as in Ref. [2], and the Coulomb
gauge $\bm{\nabla}\cdot\delta\bm{A}=0$ is assumed. Thus,

$\displaystyle\delta\bm{B}_{\perp z}$ $\displaystyle=$
$\displaystyle\bm{\nabla}_{\perp}\delta A_{\parallel
z}\times\bm{b}_{0}+\bm{b}_{0}\times\bm{\kappa}_{0}\delta A_{\parallel z}$ (35)
$\displaystyle+\bm{b}_{0}\times\nabla_{\parallel}\delta\bm{A}_{\perp
z}+(\bm{b}_{0}\times\bm{\nabla}\bm{b}_{0})\cdot\delta\bm{A}_{\perp z}\;,$

and, therefore, Eqs. (32) to (34) can be solved for $\delta\phi_{z}$, $\delta
A_{\parallel z}$ and $\delta B_{\parallel z}$ as independent field variables
uniquely defining ZFs [8]. Note that we have maintained the last term on the
RHS of Eq. (32) despite it usually vanishes assuming equilibrium
quasineutrality. However, in the present theoretical framework where
$\bar{F}_{0}$ is assumed to vary consistently with Eq. (21), equilibrium
quasineutrality is not imposed separately, while plasma quasineutrality is
satisfied overall. This means that ZFs are allowed to develop slow
spatiotemporal mean field structures due to fluctuation induced transport.

PSZS, with their micro spatiotemporal scale counterpart, and the ZFs
constitute the zonal state, introduced in Sec. 2.1, which is consistent with
the finite level of $n\neq 0$ symmetry breaking fluctuations. Here, the $n\neq
0$ fluctuation spectrum is assumed as given, but can generally be computed by
means of standard nonlinear gyrokinetic theory.

##  3 Zonal state self-consistent evolution

In order to illuminate the evolution of the zonal state, let us focus on the
case where a given $n\neq 0$ spectrum is assumed.

Following [2], we adopt low-$\beta$ ordering with good separation of SAW and
compressional Alfvén wave frequencies and we calculate $\delta B_{\parallel}$
from the perpendicular pressure balance equation:

$\nabla_{\perp}\left(B_{0}\delta B_{\|z}+4\pi\delta P_{\perp z}\right)\simeq
0,$ (36)

where $\delta P_{\perp z}$ represents the perpendicular pressure perturbation.
Having solved for $\delta B_{\|z}$ explicitly, the fluctuation spectrum is
entirely described by the scalar potential $\delta\phi_{z}$ and the parallel
vector potential $\delta A_{\parallel z}$. Furthermore, for the sake of
simplicity, we also assume the ZS is predominantly characterized by
$\delta\phi_{z}$. Consequently, in what follows, we will describe the zonal
state by means of the scalar potential ZFs only. Thus, we are left with the
solution of the zonal component of the quasi neutrality condition, i.e., Eq.
(32) and of the corresponding particle responses. The equations obtained
below, may thus be readily adopted for discussing electrostatic turbulence
and, in particular, can be used to describe GAM/EGAM (energetic particle
induced GAM [40, 41]) physics (cf. A). To this aim, we rewrite the gyrokinetic
equation for the non adiabatic drift/banana center distribution function,
i.e.:

$(\partial_{t}+v_{\parallel}\nabla_{\parallel})\delta
g_{Bz}=-e^{iQ_{z}}\left[\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\partial_{t}\delta\phi_{z}+N.L.\right]$
(37)

where, for the sake of brevity, the nonlinear terms
$\delta\dot{\bm{X}}\cdot\bm{\nabla}\delta
F+\delta\dot{\mathcal{E}}\partial_{\mathcal{E}}\delta F$ have been indicated
as $N.L.$. Introducing the lifting of a generic scalar field to the particle
phase-space [1] and the action angle coordinates $\vartheta_{c}$ and
$\zeta_{c}$, i.e., such that $\omega_{b}=\dot{\vartheta}_{c}$ and
$\dot{\zeta}_{c}=\bar{\omega}_{d}$ where $\omega_{b}$ and $\bar{\omega}_{d}$
are, respectively, the bounce/transit frequency and the precession drift
frequency, we obtain:

$(\partial_{t}+\omega_{b}\partial_{\vartheta c})\delta
g_{Bz}=-e^{iQ_{z}}\left[\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\partial_{t}\delta\phi_{z}+N.L.\right]_{z}\;.$
(38)

For now, we neglect the nonlinear term and Fourier decompose the RHS with
respect to the $\vartheta_{c}$ coordinate. Meanwhile, we introduce the
$\delta\hat{G}_{l}$ function which is connected to the Fourier series of the
scalar potential by the following definition:

$-e^{iQ_{z}}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\partial_{t}\delta\phi_{z}\equiv\sum_{l}e^{il\vartheta_{c}}\partial_{t}\delta\hat{G}_{l}.$
(39)

Here, the coefficients in the Fourier series, $\delta\hat{G}_{l}$, can be
calculated as:

$\delta\hat{G}_{l}\equiv\frac{1}{2\pi}\oint
d\vartheta_{c}e^{-il\vartheta_{c}}\left[-\frac{e}{m}e^{iQ_{z}}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\phi_{z}\right]=-\overline{e^{-il\vartheta_{c}+iQ_{z}}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\phi_{z}}.$
(40)

It can be readily shown that the spectral representation of the linear
solution of Eq. (38) reads:

$\delta
g_{Bz}=\sum_{l}\frac{\omega_{z}}{\omega_{z}-l\omega_{b}}e^{il\vartheta_{c}}\delta\hat{G}_{l}\;,$
(41)

where $\omega_{z}\equiv i\partial_{t}$ is the ZFs characteristic frequency
that must be intended as an operator. Thus, $(\omega_{z}-l\omega_{b})^{-1}$ in
Eq. (41) must be intended as the inverse operator of
$(\omega_{z}-l\omega_{b})$. We emphasize that this is a formal solution since
it requires integration along characteristics and, therefore, it involves an
integral equation. The same procedure can be straightforwardly applied to
solve the equation including the nonlinear term, and the corresponding
solution can be substituted into Eq. (32). Thus, restoring the species index
$s$ and explicitly denoting summation on particle species as well as summation
on $\hat{\sigma}=\pm$, where $\hat{\sigma}=v_{\parallel}/|v_{\parallel}|$ for
circulating particles, while, for magnetically trapped particles,
$\hat{\sigma}=\pm$ represents the right-/left-handed rotation of the particles
on the outer leg of their poloidal orbit, the magnetic flux surface averaged
Eq. (32) reads555By magnetic flux surface average we mean
$[...]_{\psi}=(2\pi/V_{\psi}^{\prime})\int_{0}^{2\pi}{\cal J}(...)d\theta$
with $V_{\psi}^{\prime}=2\pi\int_{0}^{2\pi}{\cal J}d\theta$. The following
equation, for simplicity, does not report a factor
$(4\pi^{2}/V_{\psi}^{\prime})$ that should appear in front of the double sum,
$\sum_{s}\sum_{\hat{\sigma}}$, after flux surface averaging of Eq. (32).:

$\displaystyle\sum_{s}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}\frac{e_{s}^{2}}{m_{s}}\left(\overline{\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}\delta\phi_{z}}-\sum_{l}\frac{\omega_{z}}{\omega_{z}-l\omega_{bs}}\overline{e^{il\vartheta_{c}-iQ_{zs}}J_{0}}\overline{e^{-il\vartheta_{c}+iQ_{zs}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\delta\phi_{z}}\right)$
$\displaystyle+\sum_{s}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\frac{1}{d\psi/dr}\frac{\partial}{\partial
r}\overline{\left[\frac{e_{s}^{2}}{m_{s}}\frac{2\mu}{\Omega_{s}^{2}}\frac{\partial\bar{F}_{0s}}{\partial\mu}\left(\frac{J_{0}^{2}-1}{\lambda_{s}^{2}}\right)\tau_{bs}\frac{d\psi}{dr}\frac{\partial}{\partial
r}\delta\phi_{z}\right]}=$ $\displaystyle\sum_{s}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}e_{s}\sum_{l}\frac{i}{\omega_{z}-l\omega_{bs}}\overline{e^{il\vartheta_{c}-iQ_{zs}}J_{0}}\,\overline{e^{-il\vartheta_{c}+iQ_{zs}}N.L.}.$
(42)

Note that, here, $\sum_{\hat{\sigma}}$ applies to circulating particles only,
since it is reabsorbed by the bounce averaging for trapped particles.
Furthermore, $\tau_{bs}=2\pi/\omega_{bs}$ and, for simplicity, we have ignored
the possible contribution due to breaking the PSZS quasineutrality, discussed
above. That contribution can be easily restored, if needed, along with the
contribution of sources and collisions by letting $N.L.\rightarrow
N.L.-(C_{g}+{\cal S})$ on the RHS of Eq. (42). This expression has been
derived by using a minimal set of assumptions that quite reasonably describe
the self-consistent evolution of the ZS and, therefore, its generality make it
suitable for various applications since it allows to describe an arbitrary
$\overline{F}_{0}$ while retaining realistic magnetic geometry effects. We
will illustrate some of these applications in A.

In the following, we explore the low-frequency response obtained focusing on
the $l=0$ component. This is clearly the response that is directly connected
with transport. In particular, the hence obtained linear terms can be
expressed in a compact form introducing the plasma polarizability for
$s$-species, $\chi_{zs}$666 $\chi_{zs}$, is connected with the usual
definition of susceptibility, $\chi_{s}$, via the relation
$\chi_{zs}=(1+\chi_{s})k_{r}^{2}\lambda_{D}^{2}$, with
$\lambda_{D}^{2}=T/(4\pi ne^{2})$ the Debye length., defined as:

$\displaystyle\chi_{zs}\left[\delta\phi_{z}\right]_{\psi}$
$\displaystyle\equiv$
$\displaystyle-\frac{4\pi^{2}}{V_{\psi}^{\prime}}\frac{T_{s}}{n_{s}m_{s}}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}\left\\{\overline{\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}\delta\phi_{z}}-\overline{e^{-iQ_{zs}}J_{0}}\overline{e^{iQ_{zs}}J_{0}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}\delta\phi_{z}}\right.$
(43) $\displaystyle\left.+\frac{1/\tau_{bs}}{d\psi/dr}\frac{\partial}{\partial
r}\overline{\left[\frac{2\mu}{\Omega_{s}^{2}}\frac{\partial\bar{F}_{0s}}{\partial\mu}\left(\frac{J_{0}^{2}-1}{\lambda_{s}^{2}}\right)\tau_{bs}\frac{d\psi}{dr}\frac{\partial}{\partial
r}\delta\phi_{z}\right]}\right\\}\;.$

This equation becomes a closed expression for $\chi_{zs}$ once
$\delta\hat{\phi}_{z}\equiv\delta\phi_{z}-\left[\delta\phi_{z}\right]_{\psi}$
is given777Please, note the difference between bounce, $\overline{(...)}$, and
flux averaging, $[...]_{\psi}$, although they are the same at the lowest order
for well circulating particles. The difference between $\widetilde{(...)}$ and
$\widehat{(...)}$ follows consequently.. This can be obtained from the
component of the quasineutrality condition that is varying along the flux
surface, and it can be shown that
$|\delta\hat{\phi}_{z}|\ll|\left[\delta\phi_{z}\right]_{\psi}|$ in the long
wavelength limit, $|Q_{zs}|\ll 1$ (cf., e.g., Ref. [42, 43]). In particular,
$\delta\hat{\phi}_{z}\rightarrow 0$ for $T_{e}/T_{i}\rightarrow 0$. Equation
(43), valid for arbitrary wavelength, generalizes to arbitrary geometry and
distribution functions the plasma polarizability expressions at short
wavelengths studied recently [22, 23, 24, 25]. Introducing the $s$-species
polarization density

$[\delta n_{{\rm
pol}s}]_{\psi}=-\frac{n_{s}e_{s}}{T_{s}}\chi_{zs}\left[\delta\phi_{z}\right]_{\psi}\;,$

the $l=0$ flux surface averaged quasineutrality condition, Eq. (42), can be
rewritten as

$\sum_{s}e_{s}\partial_{t}[\delta n_{{\rm
pol}s}]_{\psi}=\sum_{s}e_{s}\left[\bm{\nabla}\cdot\bm{\Gamma}_{N.L.s}\right]_{\psi}\;,$
(44)

where we have introduced the flux surface averaged divergence of the
$s$-species particle flux due to nonlinear interactions:

$\displaystyle\left[\bm{\nabla}\cdot\bm{\Gamma}_{N.L.s}\right]_{\psi}$
$\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{V_{\psi}^{\prime}}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}\overline{e^{-iQ_{zs}}J_{0}}\,\overline{e^{iQ_{zs}}N.L.}$
(45) $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{V_{\psi}^{\prime}}\frac{\partial}{\partial\psi}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}\left[\overline{\left(e^{-iQ_{zs}}J_{0}\right)}\;\overline{\left(e^{iQ_{zs}}\delta\dot{\bar{\psi}}_{s}\delta
F_{s}\right)}\right.$ $\displaystyle\hskip
50.00008pt\left.-\overline{\frac{\partial\bar{\psi}_{s}}{\partial{\cal
E}}\left(e^{-iQ_{zs}}J_{0}\right)}\;\overline{\left(e^{iQ_{zs}}\delta\dot{\cal
E}\delta F_{s}\right)}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{V_{\psi}^{\prime}}\frac{\partial}{\partial\psi}\left[V_{\psi}^{\prime}\left\langle\left(e^{-iQ_{zs}}J_{0}\right)\;\overline{\left(e^{iQ_{zs}}\delta\dot{\bar{\psi}}_{s}\delta
F_{s}\right)}\right\rangle_{v}\right.$ $\displaystyle\hskip
30.00005pt\left.-V_{\psi}^{\prime}\left\langle\frac{\partial\bar{\psi}_{s}}{\partial{\cal
E}}\left(e^{-iQ_{zs}}J_{0}\right)\;\overline{\left(e^{iQ_{zs}}\delta\dot{\cal
E}\delta F_{s}\right)}\right\rangle_{v}\right]_{\psi}\;.$

The different forms on the RHS are all equivalent and are given here to
illuminate the conservation properties of these expressions. Recalling that,
at the leading order,
$\delta\dot{\psi}=(B_{0}/B_{\parallel}^{*})c\partial_{\zeta}\left\langle\delta
L_{g}\right\rangle$, the last equation demonstrates that only toroidal
symmetry breaking fluctuations drive a finite flux surface averaged particle
transport in tokamaks at the corresponding leading order. Thus, the present
analysis must assume a prescribed spectrum of $n\neq 0$ fluctuations (cf. Sec.
2.3). Physically, Eq. (44) is readily interpreted as the nonlinear charge
density modification compensating the polarization charge to ensure
quasineutrality; that is, $\sum_{s}e_{s}\partial_{t}[\delta n_{{\rm
pol}s}]_{\psi}=-\sum_{s}e_{s}\partial_{t}[\delta n_{N.L.s}]_{\psi}$.
Meanwhile, without summing over all particle species, it is possible to cast
the same equation in the form of flux surface averaged particle continuity
equation:

$\partial_{t}\left[n_{s}\right]_{\psi}=\partial_{t}[\delta n_{{\rm
pol}s}]_{\psi}-\left[\bm{\nabla}\cdot\bm{\Gamma}_{N.L.s}\right]_{\psi}\;,$
(46)

where collisional neoclassical transport in the banana regime as well as
sources/sinks can be readily included in the expression above by letting
$N.L.\rightarrow N.L.-(C_{g}+{\cal S})$, as discussed below Eq. (42). Note
that $\left[n_{s}\right]_{\psi}$ on the LHS of Eq. (46) is the total particle
density of the $s$-species, since the fluctuation induced nonlinear particle
flux includes both micro- as well as meso- and macro-scale spatio-temporal
behaviors, consistent with Eq. (32). Thus, Eq. (46) describes the variety of
spatiotemporal scales involved in particle transport. Consistently with the
analysis of Ref. [7], polarization effects become important only on
sufficiently short scales, $k_{z}L>\delta^{-1/2}$; i.e., the meso-scales, with
$L$ the characteristic plasma macro-scale, $\rho_{L}$ the Larmor radius and
$\delta=\rho_{L}/L$ the gyrokinetic ordering parameter. Meanwhile, the leading
order flux surface averaged particle flux for symmetry breaking fluctuations
becomes

$\left[\bm{\nabla}\cdot\bm{\Gamma}_{N.L.s}\right]_{\psi}=\frac{1}{V_{\psi}^{\prime}}\frac{\partial}{\partial\psi}\left[\left\langle
V_{\psi}^{\prime}\left(e^{-iQ_{zs}}J_{0}\right)\;\overline{\left[ce^{iQ_{zs}}R^{2}\nabla\phi\cdot\nabla\left\langle\delta
L_{gs}\right\rangle\delta g_{s}\right]}\right\rangle_{v}\right]_{\psi}\;.$
(47)

Again, in the $k_{z}L<\delta^{-1/2}$ long wavelength limit, this expression
reduces to the well-known form adopted in classical analyses of fluctuation-
induced evolution of macroscopic plasma profiles [10, 44, 45, 46, 47, 32].
Following Ref. [7], the same argument can be repeated to show that classical
forms of momentum and energy transport equations are reproduced.

This demonstrates that the PSZS transport equations, Eqs. (21) and (2.2)
derived in the previous section, along with the equations for the self-
consistent determination of the ZFs, Eqs. (32) to (34), fully characterize the
ZS and, at the same time, the multi-spatiotemporal-scale nature of phase-space
transport in collisionless burning plasmas and their possible deviation from
local thermodynamic equilibrium. This description reduces to the previous
gyrokinetic theory of phase-space transport [7, 48] within the framework of
Frieman-Chen nonlinear gyrokinetic equation [4] and recovers earlier works in
the proper limit [10, 44, 45, 46, 47, 32]. Based on the nonlinear gyrokinetic
theory with Hamiltonian description of particle motion accurate up ${\cal
O}(\delta^{2})$ [27, 5], the present PSZS transport equations are valid for
even longer than the characteristic transport time scale, ${\cal
O}(\delta^{-3})\Omega^{-1}$, and typically hold on times $<{\cal
O}(\delta^{-4})\Omega^{-1}$.

##  4 Conclusions and discussion

In this article, we have presented a comprehensive study of plasma transport
processes in fusion plasmas using the phase-space zonal structure (PSZS)
transport theory [7, 8]. We addressed the limitations of current numerical
frameworks which are computationally expensive and often limited in their
ability to capture long-time scale dynamics and non-local (global) behaviors.
To overcome these challenges, we developed the PSZS transport theory, which
provides a proper definition of the plasma nonlinear equilibrium distribution
function by considering slowly evolving structures in the phase-space. The
PSZS theory allows for the derivation of the usual plasma transport equations
as a limiting case when the deviation from the local Maxwellian is small, as
demonstrated in previous work [7]. However, in the general case, PSZS is not
associated with a reference Maxwellian, as it results from the competition
between resonantly induced nonlinear transport, sources, and weakly
collisional effects, necessitating a phase-space description.

Applying the PSZS transport theory, we derived the evolution equation for the
zonal state (ZS), representing the renormalized nonlinear equilibrium
consistent with toroidal symmetry breaking fluctuations and transport time
scale ordering. Specifically, we defined the two components of the ZS, namely
the PSZS and the zonal electromagnetic fields (ZFs). Moreover, applying the
Chew Goldberger Low (CGL) description, we derived the self-consistent
modifications of the reference magnetic equilibrium using the push forward
representation of the macro-/meso-scopic component moments of the PSZS.

As an example of the theoretical framework, we discuss the self-consistent
evolution of the ZS with a given spectrum of toroidally symmetry breaking
perturbations and ZFs dominated by the scalar potential response. We derived
expressions for the plasma polarizability that are applicable to arbitrary
geometry and equilibrium distribution functions, and discussed the features of
transport equations on the different spatial scales that are involved in the
problem. Geodesic acoustic mode (GAM) and Energetic particle driven geodesic
acoustic mode (EGAM) do not belong to the ZS due to their fast time variation
and finite collisionless damping/drive. Nonetheless, finite amplitude GAM/EGAM
may nonlinearly impact on the ZS evolution. Thus, we have added a detailed
Appendix on this problem, where interested readers can find a discussion of
the GAM/EGAM peculiar physics. In particular, we give a general expression for
the linear dielectric response of GAMs. Furthermore, we illustrate examples of
GAM/EGAM nonlinear dynamics, which could be readily adopted to investigate
problems of practical interest in general geometry and with arbitrary
energetic particles (EP) distribution functions, such as the EGAM decay into
two GAMs recently observed in low-collisionality LHD plasmas [49] and the
self-consistent EGAM frequency sweeping [50].

In conclusion, the PSZS transport theory provides a promising approach to
understanding EP transport processes in fusion plasmas. The derived equations
for the ZS and the associated modifications to the equilibrium provide a
comprehensive framework for studying plasma nonlinear equilibrium and its
evolution due to transport processes. This theoretical framework opens new
possibilities for developing advanced reduced EP transport models capable of
capturing the long-time scale evolution of burning plasmas and providing
insight into the non-locality of transport processes. Notably, a recent
advancement in this field is the proposed Dyson-Schrödinger transport Model
(DSM) [8]. PSZS fluxes, computed using the DAEPS-FALCON suite of codes [51,
52], have been calculated within the LIGKA EP workflow [53, 54] considering
realistic Tokamak configurations. This is a crucial step towards the practical
implementation of the PSZS transport theory in realistic geometry. Based on a
gyrokinetic description for the underlying perturbations, employing general EP
distributions functions and using saturation rules obtained from gyrokinetic
non-linear codes will allow us to construct a quantitative and predictive
reduced EP transport model for the interpretation of present-day experimental
results and the investigation of future burning plasmas.

Future research directions include the derivation of general orbit-averaged
source and collision terms on the analytical side. On the numerical side, PSZS
diagnostics have been developed for global gyrokinetic and hybrid codes such
as HMGC and ORB5 [19], enabling the study of phase-space transport processes
during nonlinear gyrokinetic simulation. At present, the EP workflow
calculates the PSZS evolution within the kick model [55] approximation.
Further development of reduced transport models for the PSZS involves the
development of a solver for the DSM [8] and the inclusion of nonlinear
corrections into the governing equations of the EP workflow. More generally, a
comprehensive gyrokinetic transport solver on long time scale can be developed
by means of subcycling and restart of nonlinear gyrokinetic simulations in the
updated ZS computed within the present theoretical framework adopting the
numerically computed phase-space fluxes. These advancements will contribute to
a deeper understanding of EP transport in fusion plasmas and facilitate the
development of more accurate predictive models including core turbulent
transport.

This work was carried out within the framework of the EUROfusion Consortium
and received funding from Euratom research and training programme 2014–2018
and 2019–2020 under Grant Agreement No. 633053 (Project No. WP19-ER/ENEA-05).
This work has been carried out within the framework of the EUROfusion
Consortium, funded by the European Union via the Euratom Research and Training
Programme (Grant Agreement No 101052200 – EUROfusion; Projects No. ENR-
MOD.01.MPG and AC-TSVV.10.MPG). Views and opinions expressed are however those
of the author(s) only and do not necessarily reflect those of the European
Union or the European Commission. Neither the European Union nor the European
Commission can be held responsible for them. This work was also supported by
the National Science Foundation of China Project No. 12261131622. This work
was supported in part by the Italian Ministry of Foreign Affairs and
International Cooperation, grant number CN23GR02 and by the MMNLP project CSN4
of INFN, Italy.

##  Appendix A Applications to geodesic acoustic mode physics

In this appendix, as further illustration of the strength and usefulness of
the present theoretical framework, we focus on selected applications of Eq.
(42) derived above. Despite these results are not entirely novel, their
compact derivation and validity for general geometry and distribution
functions demonstrates the practical implications of the present approach.
Firstly, we derive an expression for the linear dielectric response of GAM
(Geodesic Acoustic Mode) oscillations, which yields the well-known results for
Maxwellian distribution functions and circular equilibria [56, 57] and, next,
we extend the results of Ref. [43] describing the generation of zero frequency
ZFs by GAMs oscillations. Additionally, we investigate the modulation of GAMs
by ZFs showing that the presence of energetic particles (EPs) or higher-order
thermal plasma finite orbit width effects are necessary for non-vanishing
nonlinear interactions as well as for second harmonic GAM generation. This
generalizes the results of [43, 50, 58, 59]. Finally, we explore the nonlinear
dynamics of EGAM highlighting the importance of PSZS and the ZS by considering
the nonlinear term induced by PSZS as a nonlinear equilibrium. We describe the
evolution of the zonal state, accounting for the action of sources,
collisions, and the emission and re-absorption of the GAM/EGAM fluctuations.

###  A.1 Linear dielectric response of energetic particle driven geodesic
acoustic mode

In order to calculate the linear dispersion response, we rewrite Eq. (41)
taking the decomposition
$\delta\phi_{z}=\delta\phi_{G}=\left[\delta\phi_{G}\right]_{\psi}+\delta\hat{\phi}_{G}$
explicitly into account, denoting that the scalar potential refers to
GAM/EGAM. More precisely,

$\displaystyle\delta g_{BG}$ $\displaystyle=$
$\displaystyle-e^{iQ_{G}}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\phi_{G}+\sum_{l}\left[\frac{l\omega_{b}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{b}^{2}}i\sin
l\vartheta_{c}+\frac{l^{2}\omega_{b}^{2}}{\omega_{G}^{2}-l^{2}\omega_{b}^{2}}\cos
l\vartheta_{c}\right]$ (48)
$\displaystyle\times\left[-\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}}\right]\left[\delta\phi_{G}\right]_{\psi}+\sum_{l}\left[\frac{l\omega_{b}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{b}^{2}}i\sin
l\vartheta_{c}\right.$
$\displaystyle\left.+\frac{l^{2}\omega_{b}^{2}}{\omega_{G}^{2}-l^{2}\omega_{b}^{2}}\cos
l\vartheta_{c}\right]\left[-\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\hat{\phi}_{G}}\right]$
$\displaystyle+\sum_{l}\left[\frac{l\omega_{b}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{b}^{2}}\cos
l\vartheta_{c}+\frac{l^{2}\omega_{b}^{2}}{\omega_{G}^{2}-l^{2}\omega_{b}^{2}}i\sin
l\vartheta_{c}\right]$ $\displaystyle\times\left[\overline{e^{iQ_{G}}i\sin
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\hat{\phi}_{G}}\right]\;,$

where $\delta g_{BG}$ denotes the $\delta g_{Bz}$ response to GAM, the
subscript $G$ in $Q_{G}$ reminds that the radial shift operator acts on GAM,
and we have assumed an up-down symmetric equilibrium for simplicity but
without loss of generality, since the general case could be readily restored
at the expense of more complicated formal expressions. Substituting Eq. (48)
back into the linearized Eq. (32) for the varying component on the considered
magnetic flux surface, we can write

$\displaystyle\left[-\frac{n_{e}e^{2}}{T_{e}}+\sum_{s}\left\langle\frac{e_{s}^{2}}{m_{s}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}\right\rangle_{v}\right]\delta\hat{\phi}_{G}+\sum_{s}\sum_{l}\left\langle\left[\frac{l\omega_{bs}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}i\sin
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}\right.\right.$ $\displaystyle\hskip
20.00003pt\left.\left.+\frac{\omega_{G}^{2}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\left(\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}-\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\right)\right]\right.$ $\displaystyle\hskip
20.00003pt\left.\times\left[-\overline{e^{iQ_{Gs}}\cos
l\vartheta_{c}\frac{e_{s}^{2}}{m_{s}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}}\right]\right\rangle_{v}\left[\delta\phi_{G}\right]_{\psi}+\sum_{s}\sum_{l}\left\langle\left[\frac{l\omega_{bs}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\right.\right.$
$\displaystyle\hskip 20.00003pt\times\left.i\sin
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}+\frac{\omega_{G}^{2}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\left(\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}-\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\right)\right]$ $\displaystyle\hskip
20.00003pt\left.\times\left[\overline{-e^{iQ_{Gs}}\cos
l\vartheta_{c}\frac{e_{s}^{2}}{m_{s}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\delta\hat{\phi}_{G}}\right]\right\rangle_{v}+\sum_{s}\sum_{l}\left\langle\left[\frac{l\omega_{bs}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\right.\right.$
$\displaystyle\hskip 20.00003pt\times\left.\left(\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}-\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\right)+\frac{\omega_{G}^{2}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}i\sin
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}\right]$ $\displaystyle\hskip
20.00003pt\left.\times\left[\overline{e^{iQ_{Gs}}i\sin
l\vartheta_{c}\frac{e_{s}^{2}}{m_{s}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\delta\hat{\phi}_{G}}\right]\right\rangle_{v}=0\;,$
(49)

where, for simplicity, we have assumed Maxwellian electrons. Equation (49) is
readily solved for $\delta\hat{\phi}_{G}$ as a function of
$\left[\delta\phi_{G}\right]_{\psi}$ reducing to well known results, e.g. [60,
61, 42], for Maxwellian ions in the long wavelength limit. Meanwhile, the flux
surface averaged quasineutrality condition can be written as

$\displaystyle\frac{1}{V_{\psi}^{\prime}}\frac{1}{d\psi/dr}\frac{\partial}{\partial
r}\left[V_{\psi}^{\prime}\sum_{s}\rho_{Ls}^{2}\frac{n_{s}e_{s}^{2}}{T_{s}}D_{Gs}\frac{d\psi}{dr}\frac{\partial}{\partial
r}\left[\delta\phi_{G}\right]_{\psi}\right]=\frac{4\pi^{2}}{V_{\psi}^{\prime}}\sum_{s}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}e_{s}$
$\displaystyle\times\sum_{l}\frac{i}{\omega_{G}-l\omega_{bs}}\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\,\overline{e^{-il\vartheta_{c}+iQ_{Gs}}N.L.}\;,$
(50)

where $\rho_{Ls}^{2}=(T_{s}/m_{s})/\bar{\Omega}_{s}^{2}$, the temperature
$T_{s}$ is defined as $T_{s}\equiv n_{s}^{-1}\left\langle 2m_{s}\mu
B_{0}\bar{F}_{0s}\right\rangle_{v}$ for a generic non-Maxwellian distribution
function, $\bar{\Omega}_{s}$ is the cyclotron frequency computed at the on-
magnetic-axis magnetic field $B_{0}=\bar{B}_{0}$, and $D_{Gs}$ is the
$s$-species contribution to the GAM/EGAM dispersion response, expressed as,
noting Eq. (49):

$\displaystyle
D_{Gs}\left[\delta\phi_{G}\right]_{\psi}=\frac{4\pi^{2}}{V_{\psi}^{\prime}}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}\left\\{\left[\overline{\left[\frac{2\mu\bar{B}_{0}^{2}}{n_{s}B_{0}}\left(\frac{J_{0}^{2}-1}{\lambda_{s}^{2}}\right)\left(\frac{\partial\bar{F}_{0s}}{\partial{\cal
E}}+\frac{1}{B_{0}}\frac{\partial\bar{F}_{0s}}{\partial\mu}\right)\right]}\right.\right.$
$\displaystyle\left.+\sum_{l}\frac{l^{2}\omega_{bs}^{2}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\overline{\left[e^{iQ_{Gs}}\cos
l\vartheta_{c}\frac{\bar{\Omega}_{s}^{2}}{n_{s}k_{r}^{2}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\right]}\right]\left[\delta\phi_{G}\right]_{\psi}$
$\displaystyle\left.+\sum_{l}\frac{\omega_{G}^{2}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\overline{\left[e^{iQ_{Gs}}\cos
l\vartheta_{c}\frac{\bar{\Omega}_{s}^{2}}{n_{s}k_{r}^{2}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\delta\hat{\phi}_{G}\right]}\right]$
$\displaystyle\left.\left.-\sum_{l}\frac{l\omega_{bs}\omega_{G}}{\omega_{G}^{2}-l^{2}\omega_{bs}^{2}}\overline{\cos
l\vartheta_{c}J_{0}e^{-iQ_{Gs}}}\overline{\left[e^{iQ_{Gs}}i\sin
l\vartheta_{c}\frac{\bar{\Omega}_{s}^{2}}{n_{s}k_{r}^{2}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\delta\hat{\phi}_{G}\right]}\right]\right\\}\;.$
(51)

Note that electrons do not give contribution the GAM/EGAM dispersion response
since they cannot respond to $n=0$ perturbations in the GAM/EGAM frequency
range, as it is well known. Equation (51) generalizes previously derived
expressions of the GAM/EGAM dispersion relation [56, 57] (cf. Ref. [50] for a
recent review) to the case of general geometry and distribution functions, and
recovers them in the proper limit; e.g., for circular cross section tokamak
equilibria, where, upon expanding $e^{iQ_{Gs}}\simeq 1+iQ_{Gs}$ in the long
wavelength limit,

$iQ_{Gs}\simeq\left[1+\left(1+\frac{\mu\bar{B}_{0}}{\bar{v}_{\parallel}^{2}}\right)\frac{r}{R_{0}}\cos\theta\right]\frac{qR_{0}}{r}\frac{\bar{v}_{\parallel}}{\bar{\Omega}_{s}}\partial_{r}\;,$
(52)

with $R=R_{0}$ denoting the magnetic axis, $\bar{v}_{\parallel}$ the parallel
velocity at $\bar{B}_{0}$, and we also have
$V_{\psi}^{\prime}=4\pi^{2}qR_{0}/\bar{B}_{0}$ and $\tau_{b}=2\pi
qR_{0}/|\bar{v}_{\parallel}|$ for well circulating particles. In fact,
assuming $|\omega_{G}|\gg\omega_{b}$ and one single ion species, Eq. (49)
yields

$\delta\hat{\phi}_{G}\simeq
2\frac{T_{e}}{T_{i}}\frac{T_{i}/m_{i}}{\bar{\Omega}_{i}\omega_{G}}\frac{i}{R_{0}}\sin\theta\frac{\partial}{\partial
r}\left[\delta\phi_{G}\right]_{\psi}\;,$

having noted that $\vartheta_{c}\simeq\hat{\sigma}\theta$ for well circulating
particles, while Eq. (51) reduces to

$D_{Gi}\simeq
1-\frac{2T_{i}/m_{i}}{R_{0}^{2}\omega_{G}^{2}}\left(\frac{7}{4}+\frac{T_{e}}{T_{i}}\right)\;,$

from which the leading order GAM frequency can be obtained. The GAM
collisionless damping and/or resonant EGAM excitation by phase-space
anisotropic EPs can be obtained from the wave-particle resonances embedded in
Eq. (51). Meanwhile, GAM/EGAM collisional damping is readily restored by
letting $N.L.\rightarrow N.L.-(C_{g}+{\cal S})$ on the RHS of Eq. (50) (cf.
Eq. (42) in Sec. 3). Most importantly, however, the formally nonlinear term on
the RHS of Eq. (50) allows us to discuss the relative role of different
processes contributing to GAM/EGAM nonlinear dynamics [43, 50], which are
addressed in the next two subsections.

###  A.2 Zero frequency zonal flow generation by geodesic acoustic modes

Consider the generation of zero frequency zonal flow by self-modulation of
interacting GAMs. In particular, let us look at the flux surface averaged
quasineutrality condition, Eq. (42), in the form of Eq. (44); i.e.,

$\sum_{s}\frac{n_{s}e^{2}_{s}}{T_{s}}\chi_{zs}\left[\delta\phi_{z}\right]_{\psi}=-\sum_{s}e_{s}\partial_{t}^{-1}\left[\bm{\nabla}\cdot\bm{\Gamma}_{N.L.s}\right]_{\psi}\;,$
(53)

where $\chi_{zs}$ is the general polarizability expression derived above in
Eq. (43). In order to calculate the nonlinear flux due to GAM, we assume that
mode frequency is much larger than bounce/transit frequency,
$|\omega_{G}|\gg\omega_{b}$, and, thus, from Eq. (48),

$\delta
g_{BG}\simeq-e^{iQ_{G}}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\phi_{G}+i\sum_{l}\frac{l\omega_{b}}{\omega_{G}}\sin
l\vartheta_{c}\left[-\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}}\right]\left[\delta\phi_{G}\right]_{\psi}\;,$
(54)

up to first order in the $\omega_{b}/\omega_{G}$ expansion. Note that all
bounce harmonics are retained due to the fact that GAM are characterized by
finite frequency and that, similarly to Eq. (48), we have considered an up-
down symmetric equilibrium for simplicity but without loss of generality. Now,
let’s note that $\delta\hat{\phi}_{G}={\cal
O}(k_{z}\rho_{L})[\delta\phi_{G}]_{\psi}$ for GAM [42, 43] and, thus, that
linear as well as nonlinear dynamics are dominated by finite orbit width
effects. As a consequence, the nonlinear flux due to GAM on the RHS of Eq.
(53) is dominated by the $\propto
e^{iQ_{z}}\dot{\theta}_{z}\partial_{\theta}\delta F_{z}$ term. Noting also
that

$\delta\dot{\theta}_{z}=-\frac{cRB_{\phi}}{{\cal
J}B_{0}B_{\parallel}^{*}(d\psi/dr)}(J_{0}\delta
E_{rz})\simeq-\frac{cRB_{\phi}}{{\cal J}B_{0}^{2}(d\psi/dr)}(J_{0}\delta
E_{rz})$

at the leading order, where $\delta E_{rz}$ is the GAM radial electric field,
we have

$\delta E_{rz}\simeq\frac{1}{2}\left(\delta
E_{rG}(r,t)e^{-i\omega_{G}t}+\delta E_{rG}^{*}(r,t)e^{i\omega_{G}t}\right)\;;$
(55)

and, thus,

$\displaystyle\overline{e^{iQ_{z}}N.L.}$ $\displaystyle=$
$\displaystyle\left[\hat{\sigma}\sum_{l}\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{cRB_{\phi}}{4{\cal J}B_{0}^{2}(d\psi/dr)}(J_{0}\delta
E_{rG})^{*}}\right.$ (56)
$\displaystyle\times\left.\frac{il^{2}\omega_{b}}{(\omega_{G}+i\partial_{t})}\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\phi_{G}}+c.c.\right]$
$\displaystyle\simeq$
$\displaystyle\partial_{t}\sum_{l}\frac{l^{2}\hat{\sigma}\omega_{b}}{\omega_{G}^{2}}\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{cRB_{\phi}}{4{\cal J}B_{0}^{2}(d\psi/dr)}(J_{0}\delta
E_{rG})^{*}}$ $\displaystyle\times\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0}}{\partial\mathcal{E}}J_{0}\delta\phi_{G}}+c.c.\;,$

where we have assumed that $\partial_{\theta}\sin l\vartheta_{c}\simeq
l\hat{\sigma}\cos l\vartheta_{c}$ for well circulating particles [1]888As for
the case of up-down symmetric equilibria, this assumption simplifies notations
but can be generally relaxed when carrying out numerical quadratures, which
allow using the general map $\theta\mapsto\vartheta_{c}$ for given constants
of motion $(P_{\phi},{\cal E},\mu)$., $c.c.$ stands for complex conjugate and
$(\omega_{G}+i\partial_{t})^{-1}$ is the inverse of
$(\omega_{G}+i\partial_{t})$. Equation (53), thus, becomes

$\displaystyle\sum_{s}\frac{n_{s}e^{2}_{s}}{T_{s}}\chi_{zs}\left[\delta\phi_{z}\right]_{\psi}$
$\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{V_{\psi}^{\prime}}\frac{\partial}{\partial\psi}\sum_{s}\sum_{\hat{\sigma}}\int
d\mathcal{E}d\mu\tau_{bs}\overline{e^{-iQ_{zs}}J_{0}}$ (57)
$\displaystyle\times\sum_{l}\frac{l^{2}\hat{\sigma}\omega_{bs}}{\omega_{G}^{2}}\overline{e^{iQ_{Gs}}\cos
l\vartheta_{c}\frac{cRB_{\phi}}{4{\cal J}B_{0}^{2}}(J_{0}\phi_{G})^{*}}$
$\displaystyle\times\overline{e^{iQ_{Gs}}\cos
l\vartheta_{c}\frac{e^{2}_{s}}{m_{s}}\frac{\partial\bar{F}_{0s}}{\partial\mathcal{E}}J_{0}\delta\phi_{G}}\;.$

This expression generalizes that derived in Ref. [43] and reduces to it upon
expanding $e^{iQ_{Gs}}\simeq 1+iQ_{Gs}$ in the long wavelength limit and
noting Eq. (52) for a high aspect-ratio tokamak equilibrium. Consistent with
[43], Eq. (57) suggests that efficient generation of zero frequency zonal flow
by GAM occurs at short wavelength due to finite orbit width effects.
Meanwhile, in the long wavelength limit, the leading order response is finite
only for distribution functions that are not even in $\hat{\sigma}$. For
distribution functions that are symmetric in $\hat{\sigma}$, retaining higher
order contributions in the $e^{iQ_{Gs}}$ and $e^{-iQ_{zs}}$ expansions is
necessary for computing the leading order non-vanishing term on the RHS of Eq.
(57), as it was shown in Ref. [62] for the case of a bi-Maxwellian
$\bar{F}_{0}$, which is readily recovered from Eq. (57) in the proper limit.

###  A.3 Null modification of GAM by ZFs nor GAM second harmonic generation

Let us first consider the GAM modulation by the ZFs generated either the
process discussed in A.2 or by $n\neq 0$ toroidal symmetry breaking
fluctuations. In Eq. (50) with $|\omega_{G}|\gg\omega_{b}$, the relevant
nonlinear term reads

$\displaystyle\frac{i}{\omega_{G}}\int d{\cal
E}\tau_{bs}\overline{J_{0}N.L.}=\frac{i}{\omega_{G}}\int d{\cal E}\oint
d\theta
J_{0}\left\\{\left[\frac{cRB_{\phi}}{B_{0}v_{\parallel}}\frac{\partial}{\partial\psi}\left(J_{0}\left[\delta\phi_{z}\right]_{\psi}\right)\frac{\partial}{\partial\theta}\right.\right.$
$\displaystyle\hskip
20.00003pt\left.-\frac{\partial}{\partial\theta}\left(\frac{cRB_{\phi}v_{\parallel}}{B_{0}}\right)\frac{\partial}{\partial\psi}\left(J_{0}\left[\delta\phi_{z}\right]_{\psi}\right)\frac{\partial}{\partial{\cal
E}}\right]\delta F_{Gs}$ $\displaystyle\hskip
20.00003pt+\left[\frac{cRB_{\phi}}{B_{0}v_{\parallel}}\frac{\partial}{\partial\psi}\left(J_{0}\left[\delta\phi_{G}\right]_{\psi}\right)\frac{\partial}{\partial\theta}-\frac{\partial}{\partial\theta}\left(\frac{cRB_{\phi}v_{\parallel}}{B_{0}}\right)\right.$
$\displaystyle\hskip
20.00003pt\left.\left.\times\frac{\partial}{\partial\psi}\left(J_{0}\left[\delta\phi_{G}\right]_{\psi}\right)\frac{\partial}{\partial{\cal
E}}\right]\delta F_{zs}\right\\}\;,$ (58)

at the leading order in the ${\cal O}(k_{z}q\rho_{L}),{\cal O}(k_{G}\rho_{L})$
expansion, where $\delta F_{G}$ denotes the GAM particle response, while
$\delta F_{z}$ is the low frequency particle response consistent with the ZFs
generated nonlinearly in Eq. (57). Here, the first terms in the square
brackets represent $\delta\dot{\theta}_{z,G}\partial_{\theta}$, respectively,
while the second ones stand for $\delta\dot{\cal E}_{z,G}\partial_{\cal E}$.
Integrating by parts in $\theta$ the first terms and by parts in ${\cal E}$
the second ones, it can be recognized that this expressions vanish at the
leading order, which means that GAM cannot be modulated by ZFs, either
generated by self-modulation or by other by $n\neq 0$ toroidal symmetry
breaking fluctuations on the parallel nonlinearity time scale. This result is
consistent with the findings of Ref. [43].

Let us now reconsider the GAM self-modulation and compute the generation of
GAM second harmonic. Equation (58) with $|\omega_{G}|\gg\omega_{b}$ can be
specialized to this case and becomes, at the leading order in the ${\cal
O}(k_{G}\rho_{L})$ expansion,

$\displaystyle\frac{i}{\omega_{GII}}\int d{\cal
E}\tau_{bs}\overline{J_{0}N.L.}=\frac{i}{\omega_{GII}}\int d{\cal E}\oint
d\theta
J_{0}\left[\frac{cRB_{\phi}}{B_{0}v_{\parallel}}\frac{\partial}{\partial\psi}\left(J_{0}\left[\delta\phi_{G}\right]_{\psi}\right)\frac{\partial}{\partial\theta}\right.$
$\displaystyle\hskip
20.00003pt\left.-\frac{\partial}{\partial\theta}\left(\frac{cRB_{\phi}v_{\parallel}}{B_{0}}\right)\frac{\partial}{\partial\psi}\left(J_{0}\left[\delta\phi_{G}\right]_{\psi}\right)\frac{\partial}{\partial{\cal
E}}\right]\delta F_{Gs}\;.$ (59)

Here, $\omega_{GII}\simeq 2\omega_{G}$ stands for GAM second harmonic possibly
driven by the considered finite amplitude GAM. Again, integrating by parts in
$\theta$ the first term in square brackets and by parts in ${\cal E}$ the
second one, this expression vanishes at the leading order, which means that
GAM self-modulation cannot generate second harmonic GAM on the parallel
nonlinearity time scale. Second harmonic GAM generation becomes possible by
inclusion of EP nonlinear dynamics in the GAM self-modulation or higher order
thermal plasma finite orbit width effects. Equation (59) generalizes to shaped
geometry and arbitrary distribution functions the original result of Ref. [43,
50, 58, 59].

###  A.4 Nonlinear dynamics of energetic particle driven geodesic acoustic
modes

When looking at GAM excited by EPs, the assumption $|\omega_{G}|\gg\omega_{b}$
underlying the derivations in A.3 is not applicable any longer and that bears
consequences for the GAM-ZFs and GAM-PSZS interactions. Let us reconsider Eq.
(50) with the most general nonlinear interaction term on the RHS. We analyze
first the PSZS induced nonlinear term, noting that, in the low-frequency
limit,

$\delta F_{z}=e^{-iQ_{z}}\overline{e^{iQ_{z}}\delta F_{z}}\;.$

This means that the corresponding low-frequency $\delta F_{z}$ is a function
only of $(P_{\phi},{\cal E},\mu)$. Thus, when looking at its nonlinear
interaction with a generic fluctuation structure, including the $n=0$
GAM/EGAM, we have

$\displaystyle\left(\delta\dot{\psi}_{G}\partial_{\psi}+\delta\dot{\theta}_{G}\partial_{\theta}+\delta\dot{\cal
E}_{G}\partial_{\cal E}\right)\delta
F_{z}=-\left(\partial_{t}+\dot{\bm{X}}_{0}\cdot\bm{\nabla}\right)\left(\frac{e}{m}\left\langle\delta
L_{gG}\right\rangle\right)\frac{\partial\delta F_{z}}{\partial{\cal E}}$
$\displaystyle\hskip
20.00003pt+\left(\partial_{t}+\dot{\bm{X}}_{0}\cdot\bm{\nabla}\right)\left(\frac{RB_{\phi}\left\langle\delta
A_{\parallel gG}\right\rangle}{B_{0}}\right)\frac{\partial\delta
F_{z}}{\partial\bar{\psi}}+\frac{e}{m}\partial_{t}\left\langle\delta
L_{gG}\right\rangle\frac{\partial\delta F_{z}}{\partial{\cal E}}\;,$ (60)

where we have noted Eqs. (14) and (15). Now recall Eq. (18) along with Eqs.
(37) and (38). Thus, when computing the contribution of the first term on the
RHS above to the nonlinear interaction term in Eq. (50), we have

$\overline{e^{-il\vartheta_{c}+iQ_{G}}N.L.}=i(\omega_{G}-l\omega_{b})\overline{e^{-il\vartheta_{c}+iQ_{G}}\left(\frac{e}{m}\left\langle\delta
L_{gG}\right\rangle\right)}\frac{\partial}{\partial{\cal
E}}\overline{e^{iQ_{z}}\delta F_{z}}\;.$ (61)

A similar equation can be derived for the nonlinear contribution of the second
term on the RHS in Eq. (60). As a consequence, we can incorporate the low
frequency response into the zonal state, by allowing a fast spatial variation,
consistent with Eq. (24) as well as Eqs. (17) and (23), preserving the
structure of the governing equations. Physically, this means that the low
frequency distortion in the particle distribution function can be treated as a
nonlinear equilibrium and corresponds to the renormalization of particle
response discussed in Sec. 2.2. It also further illuminates the physical
meaning of PSZS and zonal state.

We analyze now the effect of ZFs on GAM/EGAM. At the leading order, we have

$e^{iQ_{G}}\left(\delta\dot{\psi}_{z}\partial_{\psi}+\delta\dot{\theta}_{z}\partial_{\theta}+\delta\dot{\cal
E}_{z}\partial_{\cal E}\right)e^{-iQ_{G}}\delta
F_{BG}=\left[e^{iQ_{z}}\left(\delta\dot{\theta}_{z}\partial_{\theta}+\delta\dot{\cal
E}_{z}\partial_{\cal E}\right)\right]\delta F_{BG}\;.$ (62)

Therefore, the propagator $(\omega_{z}-l\omega_{b})^{-1}$ in Eq. (41) is
renormalized as

$\left(\omega_{G}+i\partial_{t}-l\omega_{b}-\Delta_{1}\right)^{-1}\;,$ (63)

where

$\Delta_{1}=-i\overline{e^{-il\vartheta_{c}}\left[e^{iQ_{z}}\left(\delta\dot{\theta}_{z}\partial_{\theta}+\delta\dot{\cal
E}_{z}\partial_{\cal E}\right)\right]e^{il\vartheta_{c}}}\;.$ (64)

Here, $\\{...\\}^{-1}$ denotes the inverse operator and, for simplicity, we
have assumed isolated resonances, nearby which, given the results of Sec. A.3,
we can expect that the dominant nonlinear effects occur. The overlapping
resonance case can be handled by a similar approach at the price of additional
technical complications. Note that Eq. (64) describes the “shearing” effect of
the ZFs, that is the wave particle decorrelation effect near resonance due to
the poloidal flow as well as the low-frequency axisymmetric energy
redistribution due to ZFs. The latter is typically negligible for toroidal
symmetry breaking fluctuations but, in general, needs to be taken into account
for a proper treatment of the effect of $n=0$ GAM/EGAM. Following a similar
argument and Dupree’s classical approach to resonance broadening theory [63],
we can also evaluate the effect of renormalization of the propagator in Eq.
(63) by generation of second harmonic component in the particle distribution
functions. Noting that, near the $\omega_{G}=l\omega_{b}$ resonance, where
$\delta F_{BG}\simeq e^{il\vartheta_{c}}\delta F_{BG}^{(l)}$,

$\delta
F_{BGII}\simeq-i\sum_{l^{\prime}}\frac{e^{il^{\prime}\vartheta_{c}}}{(\omega_{GII}-l^{\prime}\omega_{b})}\overline{e^{-il^{\prime}\vartheta_{c}}\left[e^{iQ_{G}}\left(\delta\dot{\theta}_{G}\partial_{\theta}+\delta\dot{\cal
E}_{G}\partial_{\cal E}\right)\right]e^{il\vartheta_{c}}}\,\delta
F_{BG}^{(l)}\;,$ (65)

where, noting Eq. (55),

$\displaystyle\delta\dot{\theta}_{G}$ $\displaystyle=$
$\displaystyle-\frac{cRB_{\phi}}{{\cal
J}B_{0}B_{\parallel}^{*}}\frac{J_{0}\delta E_{rG}}{2d\psi/dr}\,,$
$\displaystyle\delta\dot{\cal E}_{G}$ $\displaystyle=$
$\displaystyle\frac{cv_{\parallel}}{{\cal
J}B_{\parallel}^{*}}\frac{\partial}{\partial\theta}\left(\frac{RB_{\phi}v_{\parallel}}{B_{0}}\right)\frac{J_{0}\delta
E_{rG}}{2d\psi/dr}\,.$

Thus, the further renormalization of the propagator in Eq. (63) yields

$\left(\omega_{G}+i\partial_{t}-l\omega_{b}-\Delta_{1}-\Delta_{2}\right)^{-1}\;,$
(66)

with

$\displaystyle\Delta_{2}$ $\displaystyle=$
$\displaystyle-\sum_{l^{\prime}}\overline{e^{-il\vartheta_{c}}\left[e^{iQ_{G}}\left(\delta\dot{\theta}_{G}\partial_{\theta}+\delta\dot{\cal
E}_{G}\partial_{\cal
E}\right)^{*}\right]e^{il^{\prime}\vartheta_{c}}}(\omega_{GII}+i\partial_{t}-l^{\prime}\omega_{b}-\Delta_{1})^{-1}$
(67)
$\displaystyle\times\overline{e^{-il^{\prime}\vartheta_{c}}\left[e^{iQ_{G}}\left(\delta\dot{\theta}_{G}\partial_{\theta}+\delta\dot{\cal
E}_{G}\partial_{\cal E}\right)\right]e^{il\vartheta_{c}}}\;,$

where, for completeness, we have added long time scale dependences in the
propagator together with the effect of ZFs on wave-particle decorrelation.
Equation (67) accounts for nonlinear frequency shift as well as of resonance
broadening [63], acting as spontaneous nonlinear regulation of the minimum
resonance width for a coherent nearly-periodic spectrum [14]. In summary, the
GAM/EGAM nonlinear problem is formally linear and given by Eq. (50) with
vanishing RHS, where, however, in Eq. (51) the PSZS is given by Eq. (24) and
the renormalized propagator by Eq. (66). This is consistent with the findings
of A.3, predicting null nonlinear interactions in the GAM-ZFs and GAM-GAM
system when wave-particle interactions are neglected [43]. The nonlinear
system is closed by the PSZS evolution equation, Eqs. (21) and (2.2), which,
near the $\omega_{G}\simeq l\omega_{b}$ resonance, can be combined as:

$\displaystyle\partial_{t}\overline{e^{iQ_{z}}F_{0*}}=\left.\overline{e^{iQ_{z}}\left[C_{g}+{\cal
S}\right]}\right|_{z}+\frac{1}{2}\left[\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}J_{0}|\delta\phi_{G}|}\right]$ (68)
$\displaystyle\times\frac{\partial}{\partial{\cal
E}}\left\\{l^{2}\omega_{b}^{2}\left[(\omega_{G}-l\omega_{b})^{2}+\partial_{t}^{2}\right]^{-1}\partial_{t}\left[\overline{e^{iQ_{G}}\cos
l\vartheta_{c}\frac{e}{m}\frac{\partial\bar{F}_{0*}}{\partial\mathcal{E}}J_{0}\left|\delta\phi_{G}\right|}\right]\right\\}\;,$

where integration by parts was made to obtain the first $\overline{[...]}$ on
the RHS. Here, for simplicity, we have dropped $\Delta_{1}$ and $\Delta_{2}$
terms in Eq. (66) and noted the result of A.2 to neglect the
$\propto\dot{\theta}_{G}\partial_{\theta}$ contribution to the nonlinear
response for $F_{0*}$ symmetric in $\hat{\sigma}$. We have also assumed
$T_{e}/T_{i}\ll 1$, without loss of generality, in order to drop
$\delta\hat{\phi}_{G}$ with respect to $[\delta\phi_{G}]_{\psi}$ following
Ref. [50]. Equation (68) represents the evolution of the zonal state under the
action of sources and collisions, as well as of emission and re-absorption of
the GAM/EGAM fluctuations. In this respect, neglecting sources and collisions,
Eq. (68) is a Dyson-like equation [64, 65] as noted earlier [1, 2, 3, 8], and
its solution, which can be formally represented as a Dyson series, describes
the evolution of the ZS. Equation (68) is given in time representation and is
the extension to general geometry of the analogous equation for the evolution
of the renormalized fast ion distribution function given in Ref. [50] using
the frequency representation. Equations (50) and (68) are perhaps one of the
simplest possible illustrations of the DSM to the self-consistent evolution of
the ZS. More detailed analyses of Eqs. (50) and (68) are beyond the present
scope of illustration of simple applications of the general theoretical
framework; thus, they will be reported elsewhere.

## References

  * [1] F. Zonca et al.  “Nonlinear dynamics of phase space zonal structures and energetic particle physics in fusion plasmas”  In _New Journal of Physics_ 17.1 IOP Publishing, 2015, pp. 013052
  * [2] L. Chen and F. Zonca  “Physics of Alfvén waves and energetic particles in burning plasmas”  In _Reviews of Modern Physics_ 88.1 APS, 2016, pp. 015008
  * [3] Fulvio Zonca et al.  “Energetic particles and multi-scale dynamics in fusion plasmas”  In _Plasma Physics and Controlled Fusion_ 57.1 IOP Publishing, 2015, pp. 014024
  * [4] EA Frieman and Liu Chen  “Nonlinear gyrokinetic equations for low-frequency electromagnetic waves in general plasma equilibria”  In _Physics of Fluids (1958-1988)_ 25.3 AIP Publishing, 1982, pp. 502–508
  * [5] A J Brizard and T S Hahm  “Foundations of nonlinear gyrokinetic theory”  In _Reviews of Modern Physics_ 79.2 APS, 2007, pp. 421
  * [6] Hideo Sugama  “Gyrokinetic field theory”  In _Physics of Plasmas_ 7.2 AIP, 2000, pp. 466–480
  * [7] M. V. Falessi and F. Zonca  “Transport theory of phase space zonal structures”  In _Physics of Plasmas_ 26.2, 2019, pp. 022305  DOI: 10.1063/1.5063874
  * [8] Fulvio Zonca, Liu Chen, Matteo V. Falessi and Zhiyong Qiu  “Nonlinear radial envelope evolution equations and energetic particle transport in tokamak plasmas”  In _Journal of Physics: Conference Series_ 1785.1 IOP Publishing, 2021, pp. 012005  DOI: 10.1088/1742-6596/1785/1/012005
  * [9] Felix I Parra and Peter J Catto  “Limitations of gyrokinetics on transport time scales”  In _Plasma Physics and Controlled Fusion_ 50.6 IOP Publishing, 2008, pp. 065014  DOI: 10.1088/0741-3335/50/6/065014
  * [10] I G Abel et al.  “Multiscale gyrokinetics for rotating tokamak plasmas: fluctuations, transport and energy flows”  In _Reports on Progress in Physics_ 76.11 IOP Publishing, 2013, pp. 116201
  * [11] H. Sugama  “Modern gyrokinetic formulation of collisional and turbulent transport in toroidally rotating plasmas”  In _Reviews of Modern Plasma Physics_ 1.1 Springer ScienceBusiness Media LLC, 2017  DOI: 10.1007/s41614-017-0010-5
  * [12] Liu Chen  “Theory of magnetohydrodynamic instabilities excited by energetic particles in tokamaks”  In _Physics of Plasmas_ 1.5 AIP Publishing, 1994, pp. 1519–1522  DOI: 10.1063/1.870702
  * [13] X Wang et al.  “Nonlinear dynamics of nonadiabatic chirping-frequency Alfvén modes in tokamak plasmas”  In _Plasma Physics and Controlled Fusion_ 65.7 IOP Publishing, 2023, pp. 074001  DOI: 10.1088/1361-6587/acd71f
  * [14] F. Zonca, X. Tao and L. Chen  “Nonlinear dynamics and phase space transport by chorus emission”  In _Reviews of Modern Plasma Physics_ 5.1, 2021, pp. 8  DOI: https://doi.org/10.1007/s41614-021-00057-x
  * [15] Shaojie Wang, Zihao Wang and Tiannan Wu  “Self-organization of internal transport barrier in turbulent fusion plasmas”  In _arXiv preprint arXiv:2310.01355_ , 2023
  * [16] Ivan Calvo, José Luis Velasco and Félix Parra  “Orbit-averaged approach to fast-ion transport in stellarators”  In _private communications_ , 2023
  * [17] Allan N. Kaufman  “Quasilinear Diffusion of an Axisymmetric Toroidal Plasma”  In _The Physics of Fluids_ 15.6 AIP Publishing, 1972, pp. 1063–1069  DOI: 10.1063/1.1694031
  * [18] Alain J. Brizard and Anthony A. Chan  “Hamiltonian formulations of quasilinear theory for magnetized plasmas”  In _Frontiers in Astronomy and Space Sciences_ 9 Frontiers Media SA, 2022  DOI: 10.3389/fspas.2022.1010133
  * [19] A. Bottino et al.  “Time evolution and finite element representation of Phase Space Zonal Structures in ORB5”  In _Journal of Physics: Conference Series_ 2397.1 IOP Publishing, 2022, pp. 012019  DOI: 10.1088/1742-6596/2397/1/012019
  * [20] GF Chew, ML Goldberger and FE Low  “The Boltzmann equation and the one-fluid hydromagnetic equations in the absence of particle collisions”  In _Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences_ 236.1204, 1956, pp. 112–118  The Royal Society
  * [21] John R Cary and Alain J Brizard  “Hamiltonian theory of guiding-center motion”  In _Reviews of Modern Physics_ 81.2 APS, 2009, pp. 693
  * [22] Lu Wang and T. S. Hahm  “Generalized expression for polarization density”  In _Physics of Plasmas_ 16.6 AIP Publishing, 2009, pp. 062309  DOI: 10.1063/1.3152601
  * [23] Z X Lu et al.  “Theoretical studies and simulations of mode structure symmetry breaking in tokamak plasmas in the presence of energetic particles”  In _Plasma Physics and Controlled Fusion_ 61.4 IOP Publishing, 2019, pp. 044005  DOI: 10.1088/1361-6587/ab00c6
  * [24] Z.X. Lu et al.  “Effects of anisotropic energetic particles on zonal flow residual level”  In _Nuclear Fusion_ 61.8 IOP Publishing, 2021, pp. 086022  DOI: 10.1088/1741-4326/ac0dd2
  * [25] Y. W. Cho and T. S. Hahm  “Effect of temperature anisotropy on residual zonal flow level”  In _Physics of Plasmas_ 28.5, 2021, pp. 052303  DOI: 10.1063/5.0050824
  * [26] L Chen and F Zonca  “Nonlinear equilibria, stability and generation of zonal structures in toroidal plasmas”  In _Nuclear Fusion_ 47.8 IOP Publishing, 2007, pp. 886
  * [27] H Sugama  “Modern gyrokinetic formulation of collisional and turbulent transport in toroidally rotating plasmas”  In _Reviews of Modern Plasma Physics_ 1.1 Springer, 2017, pp. 9
  * [28] John R Cary and Robert G Littlejohn  “Noncanonical Hamiltonian mechanics and its application to magnetic field line flow”  In _Annals of Physics_ 151.1 Elsevier, 1983, pp. 1–34
  * [29] C. S. Chang, Seunghoe Ku and H. Weitzner  “Numerical study of neoclassical plasma pedestal in a tokamak geometry”  In _Physics of Plasmas_ 11.5 AIP Publishing, 2004, pp. 2649–2667  DOI: 10.1063/1.1707024
  * [30] Zhixin Lu et al.  “Full f and deltaf gyrokinetic particle simulations of Alfvén waves and energetic particle physics”  In _Plasma Physics and Controlled Fusion_ 65.3 IOP Publishing, 2023, pp. 034004  DOI: 10.1088/1361-6587/acb010
  * [31] W.W. Lee  “Gyrokinetic approach in particle simulation”  Office of ScientificTechnical Information (OSTI), 1981  DOI: 10.2172/5790272
  * [32] X Garbet, Y Idomura, L Villard and TH Watanabe  “Gyrokinetic simulations of turbulent transport”  In _Nuclear Fusion_ 50.4 IOP Publishing, 2010, pp. 043002
  * [33] A. M. Dimits et al.  “Comparisons and physics basis of tokamak transport models and turbulence simulations”  In _Physics of Plasmas_ 7.3 AIP Publishing, 2000, pp. 969–983  DOI: 10.1063/1.873896
  * [34] R. M. Kulsrud  In _Basic Plasma Physics_ I Amsterdam: North-Holland, 1983, pp. 115–146
  * [35] Liu Chen  “Theory of magnetohydrodynamic instabilities excited by energetic particles in tokamaks*”  In _Physics of Plasmas_ 1.5, 1994, pp. 1519–1522  DOI: 10.1063/1.870702
  * [36] N. Winsor, J. L. Johnson and J. M. Dawson  “Geodesic Acoustic Waves in Hydromagnetic Systems”  In _Physics of Fluids_ 11, 1968, pp. 2448–2450  DOI: 10.1063/1.1691835
  * [37] FL Hinton and MN Rosenbluth  “Dynamics of axisymmetric and poloidal flows in tokamaks”  In _Plasma Physics and Controlled Fusion_ 41.3A IOP Publishing, 1999, pp. A653
  * [38] R M Kulsrud  “MHD description of plasma: handbook of plasma physics”, 1980  DOI: 10.2172/7072860
  * [39] Moahan Murugappan et al.  “Gyrokinetic simulations of turbulence and zonal flows driven by steep profile gradients using a delta-f approach with an evolving background Maxwellian”  In _Physics of Plasmas_ 29.10 AIP Publishing LLC, 2022, pp. 103904
  * [40] H.L. Berk et al.  “Explanation of the JET n = 0 chirping mode”  In _Nuclear Fusion_ 46.10, 2006, pp. S888  DOI: 10.1088/0029-5515/46/10/S04
  * [41] R. Nazikian et al.  “Intense Geodesic Acousticlike Modes Driven by Suprathermal Ions in a Tokamak Plasma”  In _Phys. Rev. Lett._ 101 American Physical Society, 2008, pp. 185001  DOI: 10.1103/PhysRevLett.101.185001
  * [42] F. Zonca and L. Chen  “Radial structures and nonlinear excitation of geodesic acoustic modes”  In _EPL (Europhysics Letters)_ 83.3 IOP Publishing, 2008, pp. 35001  DOI: 10.1209/0295-5075/83/35001
  * [43] Liu Chen, Zhiyong Qiu and Fulvio Zonca  “On nonlinear geodesic acoustic modes in tokamak plasmas”  In _Europhysics Letters_ 107.1 EDP Sciences, IOP PublishingSocietà Italiana di Fisica, 2014, pp. 15003  DOI: 10.1209/0295-5075/107/15003
  * [44] Radu Balescu  “Transport processes in plasmas”  North-Holland, 1988
  * [45] Ker Chung Shaing  “Neoclassical quasilinear transport theory of fluctuations in toroidal plasmas”  In _Physics of Fluids (1958-1988)_ 31.8 AIP Publishing, 1988, pp. 2249–2265
  * [46] H Sugama, M Okamoto, W Horton and M Wakatani  “Transport processes and entropy production in toroidal plasmas with gyrokinetic electromagnetic turbulence”  In _Physics of Plasmas (1994-present)_ 3.6 AIP Publishing, 1996, pp. 2379–2394
  * [47] Hideo Sugama and W Horton  “Nonlinear electromagnetic gyrokinetic equation for plasmas with large mean flows”  In _Physics of plasmas_ 5.7 AIP, 1998, pp. 2560–2573
  * [48] M. V. Falessi and F. Zonca  “Gyrokinetic theory for particle and energy transport in fusion plasmas”  In _Physics of Plasmas_ 25.3, 2018, pp. 032306  DOI: 10.1063/1.5018175
  * [49] Zhiyong Qiu, Liu Chen, Fulvio Zonca and Matteo Valerio Falessi  “Evidence of ‘two plasmon’ decay of energetic particle induced geodesic acoustic mode”  In _New Journal of Physics_ 23.6 IOP Publishing, 2021, pp. 063045  DOI: 10.1088/1367-2630/ac047a
  * [50] Zhiyong QIU, Liu CHEN and Fulvio ZONCA  “Kinetic theory of geodesic acoustic modes in toroidal plasmas: a brief review”  In _Plasma Science and Technology_ 20.9 IOP Publishing, 2018, pp. 094004  DOI: 10.1088/2058-6272/aab4f0
  * [51] Yueyan Li et al.  “Physics of drift Alfvén instabilities and energetic particles in fusion plasmas”  In _Plasma Physics and Controlled Fusion_ , 2023  DOI: 10.1088/1361-6587/acda5e
  * [52] Matteo Valerio Falessi et al.  “Shear Alfvén and acoustic continuum in general axisymmetric toroidal geometry”  In _Physics of Plasmas_ 26.8 AIP Publishing, 2019, pp. 082502  DOI: 10.1063/1.5098982
  * [53] 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 effects on the MHD stability in tokamaks”  In _Journal of Computational Physics_ 226.1 Elsevier BV, 2007, pp. 447–465  DOI: 10.1016/j.jcp.2007.04.019
  * [54] V-A Popa et al.  “EP-Stability-WF: an IMAS-integrated workflow for energetic particle stability”  In _arXiv preprint arXiv:2306.08442_ , 2023
  * [55] Mario Podesta, Marina Gorelenkova and Roscoe White  “Reduced Fast Ion Transport Model For The Tokamak Transport Code TRANSP”  Office of ScientificTechnical Information (OSTI), 2014  DOI: 10.2172/1128924
  * [56] G. Y. Fu  “Energetic-Particle-Induced Geodesic Acoustic Mode”  In _Phys. Rev. Lett._ 101 American Physical Society, 2008, pp. 185002  DOI: 10.1103/PhysRevLett.101.185002
  * [57] Z. Qiu, F Zonca and Liu Chen  “Nonlocal theory of energetic-particle-induced geodesic acoustic mode”  In _Plasma Physics and Controlled Fusion_ 52.9 IOP Publishing, 2010, pp. 095003  DOI: 10.1088/0741-3335/52/9/095003
  * [58] H.S. Zhang, Z. Qiu, L. Chen and Z. Lin  “The importance of parallel nonlinearity in the self-interaction of geodesic acoustic mode”  In _Nuclear Fusion_ 49.12, 2009, pp. 125009  DOI: 10.1088/0029-5515/49/12/125009
  * [59] G. Y. FU  “On nonlinear self-interaction of geodesic acoustic mode driven by energetic particles”  In _Journal of Plasma Physics_ 77.4 Cambridge University Press, 2011, pp. 457–467  DOI: 10.1017/S0022377810000619
  * [60] Zhe Gao, K. Itoh, H. Sanuki and J. Q. Dong  “Multiple eigenmodes of geodesic acoustic mode in collisionless plasmas”  In _Physics of Plasmas_ 13.10, 2006, pp. 100702  DOI: 10.1063/1.2359722
  * [61] Hideo Sugama and T-H Watanabe  “Collisionless damping of zonal flows in helical systems”  In _Physics of Plasmas (1994-present)_ 13.1 AIP Publishing, 2006, pp. 012501
  * [62] Haijun Ren and X. Q. Xu  “Excitation of zonal flow by nonlinear geodesic acoustic mode”  In _Physics of Plasmas_ 27.3, 2020, pp. 034501  DOI: 10.1063/1.5126872
  * [63] Thurman Henry Dupree  “A perturbation theory for strong plasma turbulence”  In _The Physics of fluids_ 9.9 AIP, 1966, pp. 1773–1782
  * [64] F. J. Dyson  “The $S$ Matrix in Quantum Electrodynamics”  In _Phys. Rev._ 75, 1949, pp. 1736–1755
  * [65] Julian Schwinger  “On the Green’s functions of quantized fields. I”  In _Proceedings of the National Academy of Sciences_ 37.7, 1951, pp. 452–455