Nonlinear equilibria and transport processes in burning plasmas

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 description, separating the spatiotemporal microscale structures. The nonlinear physics of ZFs and of geodesic acoustic modes (GAMs)/energetic particle driven GAMs is then analyzed to illustrate the applications of our theory.


Introduction
One fundamental aspect of magnetized fusion plasmas is the study of physics processes that are occurring in burning conditions, where α particles are produced, e.g., by deuterium-tritium fusion reactions.Understanding the behavior of α particles and, more generally, of energetic particles (EPs) ‡ 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 timeconsuming 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-incell 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 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 magneticdrift/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 selfconsistent 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 longlived 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 indepth 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, Appendix 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.

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) ⊥ /2B 0 +. . .and the toroidal angular momentum § For simplicity of the present analysis, we assume that equilibrium radial electric field, if it exists, corresponds to sufficiently slow E × 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 P ϕ [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 v ∥ dl = J(P ϕ , E, µ) [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.

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.Ż = Ż0 +δ Ż, 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: 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 (θ, ζ, P ϕ , E 0 , µ) as phase space coordinates∥ where θ and ζ 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., B 0 = F ∇ϕ + ∇ϕ × ∇ψ where F = RB ϕ and ϕ is the toroidal angle, which is connected to ζ as ζ = ϕ − ν(ψ, θ), with ν(ψ, θ) chosen such that magnetic flux coordinates are characterized by straight magnetic field lines.We now note that, in the zonal component readily dropped.
∥ For 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.
of Eq. ( 1), we can re-write the term describing the equilibrium motion as: where J P ϕ = J (∂P ϕ /∂ψ) −1 , J = (∇ζ • ∇ψ × ∇θ) −1 is the Jacobian in flux coordinates, ψ is the poloidal magnetic flux, and the toroidal symmetry of the reference state has been used along with and the conservation of P ϕ and energy characterizing particle motion in the equilibrium magnetic field, i.e., respectively Ẋ0 • ∇P ϕ = 0 and Ė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 θ on the gyrokinetic equation while using J P ϕ 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: Recalling the governing equation in the absence of fluctuations [28,21] we can recognize that the averaging used to derive Eq. ( 3) is indeed a time averaging along the integrable particle orbit; denoted as with τ b = dθ/ θ.Thus, we obtain the following (equilibrium) orbit averaged kinetic equation: where δ Ṗϕ = δ Ẋ • ∇P ϕ , δ Ė 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 g s , 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 F z (O) , i.e.O) ] S , the relevant evolution equation is obtained by additionally extracting the macro-/meso-scopic component of Eq. (5); i.e.:

PSZS, as
¶ Here, by the ≃ notation, we mean that ψ contains the leading order expression of P ϕ .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].
where [. ..] 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: In the previous expression we have introduced the following decomposition at each instant of time: The aforementioned spatio-temporal averaging over the micro-scales is what allows us to separate F 0 (O) from F z (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 F 0 (O) , and based on Eq. ( 5), include all residual response in δF z (O) ≡ F z (O) − F 0 (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 δ E/∆E with δP ϕ /∆P ϕ ∼ O(1), where ∆E and ∆P ϕ are, respectively, PSZS characteristic scales in energy and toroidal angular momentum space; and δE and δP ϕ 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.: In particular, as already stated, the PSZS, F 0 (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 F 0 (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, δF z (O) , characterizes the transition between neighboring nonlinear equilibria, which are all undamped by collisionless dissipation [7,26] and slightly deviate from the reference F 0 (O) as schematically described in Fig. 1.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 F 0 (O) by means of Eq. ( 6), and "ensemble average" in a statistical sense [7].The distribution function ) , 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 ̸ = 0) are not explicitly mentioned as elements of the ZS, but they are self-consistently accounted for in its definition, as they have determined F z (O) during the dynamic evolution of the system.
We also note that, although PSZS can be obtained separating out the macro-/mesoscopic component of F z (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], 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 microscale component of the equilibrium orbit averaged distribution function: Note that, here, [. ..]F denotes the spatio-temporal micro-scale component of the argument such that [. ..] ≡ [. ..] S + [. ..]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 Appendix A.

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 ϕ 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 (ψ, θ, ζ).We can define the associated change of coordinates between these two representations noting that 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 Ĥ(P ϕ , θ) = H( ψ + δ ψ, θ) reads: where ψ = −(c/e)P ϕ was defined in Sec.2.1 above and δ ψ = RB ϕ v ∥ /Ω at the leading order, with its dependence on θ 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 (. ..) = τ −1 b dθ(...)/ θ, where, now, the closed poloidal orbit integral follows the constant-θ projection of the actual guiding center orbit on the ψ flux surface.This definition of τ b with respect to the orbit averaging approach, is unambiguous since it is uniquely defined being θ a dummy integration variable.Then, (. ..) where, for further clarity, we have explicitly denoted by the additional subscripts the reference value of P ϕ on the LHS, and of ψ on the RHS.Rephrasing this concept, Eq. ( 12) states that orbit average for given P ϕ (implicitly assuming given E and µ) corresponds to a proper shifted bounce averaging on the flux surface labeled by ψ.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 (ψ, θ, ζ, E, µ) 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 −iQz where Q z = RB ϕ v ∥ /Ω k z /(dψ/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 ≡ −i∂ 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 ϕ , consistent with the adopted gyrokinetic description [27,5].In particular, F z can be written as: where the function δ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 F0 is, from now on, a short notation for F (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.: where angular brackets denote gyro-phase averaging, δL = δϕ−v •δA/c, δϕ is the scalar potential, δA is the vector potential with ⊥ 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: Thus, we can re-write the toroidally symmetric component of Eq. (1) as follows: This equation, consistently with Ref. [8,4], suggests to introduce the following definition: where, as radial coordinate, we are using ψ ≡ −(c/e)P ϕ introduced earlier.From the definition above, the role of F0 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. ( 16) and ( 17) between the δF z (O) = e iQz δF z contribution that should be kept distinct from F0 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 −iQz G Bz , substitute this expression in Eq. ( 16) and apply e iQz on both sides.We find: where J ψ = −(e/c)J P ϕ is computed at the actual gyrocenter particle position; and we have noted iQ z = (ψ − ψ)∂ ψ .Considering the effect of the shift operator on D, we obtain the following kinetic equation: In fact, it can be shown that the shift operator commutes with the partial derivatives of nonlinear terms in Eq. ( 16) taking into account formal simplifications among commutators.Note that δ ψ is typically dominated by δ ψ for n ̸ = 0 symmetry breaking perturbations.The nonlinearities caused by ZFs, meanwhile, need a special attention, since δ ψz vanishes at the leading order.We will come back to this technical but important point while describing some applications of this theory in Appendix A.
Integrating over θ on a closed trajectory, the θ derivatives can be annihilated.Recalling the bounce average definition introduced previously, the following expression is finally obtained: ) 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 F0 follows directly from Eq. ( 20) and is consistent with Eq. ( 6): 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 (selfinteractions).The orbit averaged fast spatiotemporal deviation of the plasma response about the PSZS is given by: where δg z = e −iQz δg Bz , consistent with Eq. ( 17), is the nonadiabatic particle response that is connected with δF z by: Similarly, after some lengthy but straightforward algebra, one can obtain the governing equation for δg Bz = δg Bz − δg Bz [8].Equation ( 21), or the equivalent Eq. ( 6), and Eq. ( 22) 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 * ≡ F0 + e −iQz e iQz δF z F .

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: Using this expression, we can write every fluid moment in terms of its so called pushforward 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 −ρ•∇ 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 J z reads: where α is the gyrophase, T −1 gc v represents the guiding-center transformation of the velocity v and the argument of the delta function accounts for the relation between the particle position r and the guiding center position 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 ∝ F0 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: where ⊥ and ∥ denote the components perpendicular and parallel to B 0 and It is well-known that, assuming B 0 = F ∇ϕ + ∇ϕ × ∇ψ, the radial component of this expression reads: where ∆ * is the usual Grad-Shafranov operator and G(ψ) = (σ F ) 2 /2 is a flux function.Meanwhile, pressure components and F (ψ, B 0 ) function are connected by the parallel and bi-normal 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 J z as well as P ⊥ and P ∥ have been computed from.More precisely, P ⊥ and P ∥ can be calculated integrating the PSZS distribution function and, then, F (ψ, B 0 ) is obtained from the expression of the poloidal plasma current, which is also computed from F0 , 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., δϕ z , δA ∥z and δB ∥z are obtained by means of quasineutrality and Ampère equations following the well known theoretical framework described in [2]: Here, s denotes summation on all particle species, ⟨...⟩ v stands for velocity space integration, ∇ −1 ∥ is the inverse operator of ∇ ∥ , κ 0 ≡ b 0 • ∇b 0 is the magnetic field curvature, δJ ⊥z is readily obtained from Eq. ( 26), δB ⊥z and δB ∥z are expressed in terms of the fluctuating vector potential as in Ref. [2], and the Coulomb gauge ∇ • δA = 0 is assumed.Thus, and, therefore, Eqs. ( 32) to (34) can be solved for δϕ z , δA ∥z and δB ∥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 F0 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 ̸ = 0 symmetry breaking fluctuations.Here, the n ̸ = 0 fluctuation spectrum is assumed as given, but can generally be computed by means of standard nonlinear gyrokinetic theory.

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 ̸ = 0 spectrum is assumed.
Following [2], we adopt low-β ordering with good separation of SAW and compressional Alfvén wave frequencies and we calculate δB ∥ from the perpendicular pressure balance equation: where δP ⊥z represents the perpendicular pressure perturbation.Having solved for δB ∥z explicitly, the fluctuation spectrum is entirely described by the scalar potential δϕ z and the parallel vector potential δA ∥z .Furthermore, for the sake of simplicity, we also assume the ZS is predominantly characterized by δϕ 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.Appendix A).To this aim, we rewrite the gyrokinetic equation for the non adiabatic drift/banana center distribution function, i.e.: where, for the sake of brevity, the nonlinear terms δ Ẋ • ∇δF + δ Ė∂ E δF have been indicated as N.L.. Introducing the lifting of a generic scalar field to the particle phasespace [1] and the action angle coordinates ϑ c and ζ c , i.e., such that ω b = θc and ζc = ωd where ω b and ωd are, respectively, the bounce/transit frequency and the precession drift frequency, we obtain: For now, we neglect the nonlinear term and Fourier decompose the RHS with respect to the ϑ c coordinate.Meanwhile, we introduce the δ Ĝl function which is connected to the Fourier series of the scalar potential by the following definition: Here, the coefficients in the Fourier series, δ Ĝl , can be calculated as: It can be readily shown that the spectral representation of the linear solution of Eq. ( 38) reads: where ω z ≡ i∂ t is the ZFs characteristic frequency that must be intended as an operator.Thus, (ω z − lω b ) −1 in Eq. ( 41) must be intended as the inverse operator of (ω z − lω 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 σ = ±, where σ = v ∥ /|v ∥ | for circulating particles, while, for magnetically trapped particles, σ = ± represents the right-/left-handed rotation of the particles on the outer leg of their poloidal orbit, the magnetic flux surface averaged Eq. ( 32) reads + : Note that, here, σ applies to circulating particles only, since it is reabsorbed by the bounce averaging for trapped particles.Furthermore, τ bs = 2π/ω 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. → N.L. − (C g + 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, + By magnetic flux surface average we mean [...] ψ = (2π/V ′ ψ ) 2π 0 J (...)dθ with V ′ ψ = 2π 2π 0 J dθ.The following equation, for simplicity, does not report a factor (4π 2 /V ′ ψ ) that should appear in front of the double sum, s σ , after flux surface averaging of Eq. (32).its generality make it suitable for various applications since it allows to describe an arbitrary F 0 while retaining realistic magnetic geometry effects.We will illustrate some of these applications in Appendix 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, χ zs * , defined as: This equation becomes a closed expression for χ zs once δ φz ≡ δϕ z −[δϕ z ] ψ is given♯.This can be obtained from the component of the quasineutrality condition that is varying along the flux surface, and it can be shown that |δ φz | ≪ | [δϕ z ] ψ | in the long wavelength limit, |Q zs | ≪ 1 (cf., e.g., Ref. [42,43]).In particular, δ φz → 0 for T e /T i → 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 the l = 0 flux surface averaged quasineutrality condition, Eq. ( 42), can be rewritten as where we have introduced the flux surface averaged divergence of the s-species particle flux due to nonlinear interactions: dEdµτ bs e −iQzs J 0 e iQzs N.L.
* χ zs , is connected with the usual definition of susceptibility, χ s , via the relation ) the Debye length.♯ Please, note the difference between bounce, (...), and flux averaging, [...] ψ , although they are the same at the lowest order for well circulating particles.The difference between (...) and (...) follows consequently.
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, δ ψ = (B 0 /B * ∥ )c∂ ζ ⟨δL g ⟩, 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 ̸ = 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, 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: 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. → N.L. − (C g + S), as discussed below Eq. ( 42).Note that [n s ] ψ 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 > δ −1/2 ; i.e., the meso-scales, with L the characteristic plasma macro-scale, ρ L the Larmor radius and δ = ρ L /L the gyrokinetic ordering parameter.Meanwhile, the leading order flux surface averaged particle flux for symmetry breaking fluctuations becomes Again, in the k z L < δ −1/2 long wavelength limit, this expression reduces to the wellknown 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 ( 22) 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 multispatiotemporal-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 O(δ 2 ) [27,5], the present PSZS transport equations are valid for even longer than the characteristic transport time scale, O(δ −3 )Ω −1 , and typically hold on times < O(δ −4 )Ω −1 .

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.
where δg BG denotes the δ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. (A.1) back into the linearized Eq. ( 32) for the varying component on the considered magnetic flux surface, we can write where, for simplicity, we have assumed Maxwellian electrons.Equation (A.2) is readily solved for δ φG as a function of [δϕ G ] ψ 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 where ρ 2 Ls = (T s /m s )/ Ω2 s , the temperature T s is defined as T s ≡ n −1 s 2m s µB 0 F0s v for a generic non-Maxwellian distribution function, Ωs is the cyclotron frequency computed at the on-magnetic-axis magnetic field B 0 = B0 , and D Gs is the s-species contribution to the GAM/EGAM dispersion response, expressed as, noting Eq. (A.2): 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 (A.4) 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 ≃ 1 + iQ Gs in the long wavelength limit, with R = R 0 denoting the magnetic axis, v∥ the parallel velocity at B0 , and we also have V ′ ψ = 4π 2 qR 0 / B0 and τ b = 2πqR 0 /|v ∥ | for well circulating particles.In fact, assuming |ω G | ≫ ω b and one single ion species, Eq. (A.2) yields having noted that ϑ c ≃ σθ for well circulating particles, while Eq.(A.4) reduces to 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. (A.4).Meanwhile, GAM/EGAM collisional damping is readily restored by letting N.L. → N.L. − (C g + S) on the RHS of Eq. (A.3) (cf.Eq. ( 42) in Sec. 3).Most importantly, however, the formally nonlinear term on the RHS of Eq. (A.3) 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.
up to first order in the ω b /ω 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. (A.1), we have considered an up-down symmetric equilibrium for simplicity but without loss of generality.Now, let's note that δ φG = O(k z ρ L )[δϕ G ] ψ 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. (A.6) is dominated by the ∝ e iQz θz ∂ θ δF z term.Noting also that at the leading order, where δE rz is the GAM radial electric field, we have and, thus, where we have assumed that ∂ θ sin lϑ c ≃ lσ cos lϑ c for well circulating particles [1] † †, c.c. stands for complex conjugate and (ω G + i∂ t ) −1 is the inverse of (ω G + i∂ t ).Equation † † As 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 θ → ϑ c for given constants of motion (P ϕ , E, µ).
(A.6), thus, becomes This expression generalizes that derived in Ref. [43] and reduces to it upon expanding e iQ Gs ≃ 1 + iQ Gs in the long wavelength limit and noting Eq. (A.5) for a high aspectratio tokamak equilibrium.Consistent with [43], Eq. (A.10) 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 σ.For distribution functions that are symmetric in σ, retaining higher order contributions in the e iQ Gs and e −iQzs expansions is necessary for computing the leading order non-vanishing term on the RHS of Eq. (A.10), as it was shown in Ref. [62] for the case of a bi-Maxwellian F0 , which is readily recovered from Eq. (A.10) in the proper limit.
at the leading order in the O(k z qρ L ), O(k G ρ L ) expansion, where δF G denotes the GAM particle response, while δF z is the low frequency particle response consistent with the ZFs generated nonlinearly in Eq. (A.10).Here, the first terms in the square brackets represent δ θz,G ∂ θ , respectively, while the second ones stand for δ Ėz,G ∂ E .Integrating by parts in θ the first terms and by parts in 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 ̸ = 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 (A.11) with |ω G | ≫ ω b can be specialized to this case and becomes, at the leading order in the O(k G ρ L ) expansion, Here, ω GII ≃ 2ω G stands for GAM second harmonic possibly driven by the considered finite amplitude GAM.Again, integrating by parts in θ the first term in square brackets and by parts in 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 (A.12) generalizes to shaped geometry and arbitrary distribution functions the original result of Ref. [43,50,58,59].This means that the corresponding low-frequency δF z is a function only of (P ϕ , E, µ).Thus, when looking at its nonlinear interaction with a generic fluctuation structure, including the n = 0 GAM/EGAM, we have 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. (A.A similar equation can be derived for the nonlinear contribution of the second term on the RHS in Eq. (A.13).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 where, for completeness, we have added long time scale dependences in the propagator together with the effect of ZFs on wave-particle decorrelation.Equation (A.20) 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 nearlyperiodic spectrum [14].In summary, the GAM/EGAM nonlinear problem is formally linear and given by Eq. (A.3) with vanishing RHS, where, however, in Eq. (A.4) the PSZS is given by Eq. ( 24) and the renormalized propagator by Eq. (A. 19).This is consistent with the findings of Appendix 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 (22), which, near the ω G ≃ lω b resonance, can be combined as: where integration by parts was made to obtain the first [...] on the RHS.Here, for simplicity, we have dropped ∆ 1 and ∆ 2 terms in Eq. (A. 19) and noted the result of Appendix A.2 to neglect the ∝ θG ∂ θ contribution to the nonlinear response for F 0 * symmetric in σ.We have also assumed T e /T i ≪ 1, without loss of generality, in order to drop δ φG with respect to [δϕ G ] ψ following Ref.[50].Equation (A.21) 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. (A.21) 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 (A.21) 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 (A.3) and (A.21) are perhaps one of the simplest possible illustrations of the DSM to the self-consistent evolution of the ZS.More detailed analyses of Eqs.(A.3) and (A.21) are beyond the present scope of illustration of simple applications of the general theoretical framework; thus, they will be reported elsewhere.

Figure 1 .
Figure 1.Schematic picture describing the equilibrium orbit averaged distribution function F z (O) .The solid line represents the slowly varying component of the orbit averaged distribution function, while the oscillation around it corresponds to δF z (O) .

Appendix 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 Appendix A.2 or by n ̸ = 0 toroidal symmetry breaking fluctuations.In Eq. (A.3) with |ω G | ≫ ω b , the relevant nonlinear term reads

Appendix A. 4 .
Nonlinear dynamics of energetic particle driven geodesic acoustic modesWhen looking at GAM excited by EPs, the assumption |ω G | ≫ ω b underlying the derivations in Appendix A.3 is not applicable any longer and that bears consequences for the GAM-ZFs and GAM-PSZS interactions.Let us reconsider Eq. (A.3) 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, δF z = e −iQz e iQz δF z .