A gyrokinetic model for the plasma periphery of tokamak devices

TL;DR

This paper introduces a comprehensive gyrokinetic model for tokamak plasma periphery, capable of capturing complex electromagnetic fluctuations and flows, advancing the understanding of plasma behavior in fusion devices.

Contribution

It develops a second order gyrokinetic model using Lie perturbation theory and expands the distribution function onto a Hermite-Laguerre basis for detailed plasma analysis.

Findings

Accurately describes plasma behavior with strong flows and electromagnetic fluctuations.

Derives a set of coupled fluid equations from the gyrokinetic Boltzmann equation.

Provides a variational principle-based electromagnetic field evolution framework.

Abstract

A gyrokinetic model is presented that can properly describe strong flows, large and small amplitude electromagnetic fluctuations occurring on scale lengths ranging from the electron Larmor radius to the equilibrium perpendicular pressure gradient scale length, and large deviations from thermal equilibrium. The formulation of the gyrokinetic model is based on a second order description of the single charged particle dynamics, derived from Lie perturbation theory, where the fast particle gyromotion is decoupled from the slow drifts, assuming that the ratio of the ion sound Larmor radius to the perpendicular equilibrium pressure scale length is small. The collective behavior of the plasma is obtained by a gyrokinetic Boltzmann equation that describes the evolution of the gyroaveraged distribution function and includes a non-linear gyrokinetic Dougherty collision operator. The gyrokinetic model is then developed into a set of coupled fluid equations referred to as the gyrokinetic moment hierarchy. To obtain this hierarchy, the gyroaveraged distribution function is expanded onto a velocity-space Hermite-Laguerre polynomial basis and the gyrokinetic equation is projected onto the same basis, obtaining the spatial and temporal evolution of the Hermite-Laguerre expansion coefficients. The Hermite-Laguerre projection is performed accurately at arbitrary perpendicular wavenumber values. Finally, the self-consistent evolution of the electromagnetic fields is described by a set of gyrokinetic Maxwell's equations derived from a variational principle, with the velocity integrals of the gyroaveraged distribution function explicitly evaluated.

Tables

Table 1. Table 1: Comparison between the gyrokinetic model presented herein and previous gyrokinetic theories. The present model is fully electromagnetic (EM), is in both drift-kinetic (DK) and gyrokinetic (GK) regimes, second order accurate in the collisionless part (CLO) , in both ϵ italic-ϵ \epsilon and ϵ δ subscript italic-ϵ 𝛿 \epsilon_{\delta} , and includes the effects of strong flows (SF). The exact structure the phase-space volume element (PSV) is preserved. A closed and analytical form of the gyroaverage operator, ⟨ ⋅ ⟩ delimited-⟨⟩ ⋅ \left<\cdot\right> , in terms of gyro-moments (G-M) is provided, while a FLR model, Γ Γ \Gamma , is often used in other models, or a closed expression is not given. A gyrokinetic Dougherty (D) collision operator (CO) is used, while a Coulomb (C) operator will be developed in future work. The polarization and magnetization corrections in the gyrokinetic Poisson’s equations (GKP) and Ampere’s laws (GKA) are analytically evaluated at arbitrary wavelengths and as functions of G-M. In other theories, these terms appear as velocity integrals ( ∫ \int ), or are O ​ ( ϵ ⟂ 2 ) 𝑂 superscript subscript italic-ϵ perpendicular-to 2 O(\epsilon_{\perp}^{2}) accurate.

	This Work	Qin et al. (2006)	Hahm et al. (2009)	Dimits (2012)	Madsen (2013)	Jorge et al. (2017)	Mandell et al. (2018)
	Full-F	Full-F	Full-F	Full-F	Trunc. Full-F	Full-F	delta-F
DK	EM	EM	ES	ES	ES	ES	ES
GK	EM	EM	EM	EM	EM	-	ES
CLO	$O({ϵ}^2,{ϵ}_{\delta }^2)$	$O(ϵ,{ϵ}_{\delta }^2)$	$O({ϵ}^2,{ϵ}_{\delta }^2)$	$O({ϵ}^2,{ϵ}_{\delta }^2)$	$O(ϵ,{ϵ}_{\delta })$	$O\left(ϵ\right)$	$O(ϵ,{ϵ}_{\delta })$
SF	YES	YES	YES	NO	YES	YES	NO
PSV	$B_{\parallel }^{\ast }$	$B_{\parallel }^{\ast }$	$B_{\parallel }^{\ast }$	$B_{\parallel }^{\ast }$	$B_{\parallel }^{\ast }\simeq B$	$B_{\parallel }^{\ast }$	$B_{\parallel }^{\ast }\simeq B$
	G-M	NO	$J_0$	NO	${\parallel J_0\parallel }_a\simeq \Gamma$	-	G-M
CO	D	-	-	-	-	C	D
GKP	G-M	G-M	$O\left({ϵ}_{⟂}^2\right)$	$\int$	$O\left({ϵ}_{⟂}^2\right)$	G-M	G-M
GKA	G-M	-	$O\left({ϵ}_{⟂}^2\right)$	$\int$	$O\left({ϵ}_{⟂}^2\right)$	-	-

Equations

$$\frac{\omega }{{\Omega }_i}\sim \frac{e{\varphi }_1}{T_e}\sim \frac{k_{\parallel }}{k_{\perp }}\ll 1,$$

$$k_{\perp }{\rho }_i\sim 1,$$

$$ϵ=\frac{{\rho }_s}{L_P}\sim \frac{{\rho }_s}{L_{\varphi }}\ll 1,$$

$${ϵ}_B=\frac{{\rho }_s}{L_B}\sim {ϵ}^3$$

$$k_{\perp }\sim |{\nabla }_{\perp }\ln n_a|\sim |{\nabla }_{\perp }\ln T_a|\sim |{\nabla }_{\perp }\ln P|\sim |{\nabla }_{\perp }\ln \varphi |,$$

$${ϵ}_{\perp }=k_{\perp }{\rho }_s.$$

$${ϵ}_E=\frac{|{\mathbf{u}}_E|}{c_s}\sim {ϵ}_{\perp }\frac{e\varphi }{T_e}\ll 1.$$

$$\varphi ={\varphi }_0+{\varphi }_1,$$

$${ϵ}_{\delta }=\frac{e{\varphi }_1}{T_e}\ll 1.$$

$$\frac{k_{\parallel }}{k_{\perp }}\sim ϵ.$$

$$\frac{\omega }{{\Omega }_i}\sim {ϵ}_{\perp }ϵ,$$

$$\frac{|\delta \mathbf{B}|}{\hat{B}}\sim ϵ,$$

$$\frac{\left|{\partial }_t\hat{\mathbf{A}}\right|}{c_{s}B}\sim ϵ.$$

$$\delta \mathbf{B}=\nabla \times \delta {\mathbf{A}}_{\perp }+\nabla \times \left(\delta A_{\parallel }\hat{\mathbf{b}}\right),$$

$$\delta \mathbf{B}=\delta {\mathbf{B}}_{\perp }\simeq {\nabla }_{\perp }\delta A_{\parallel }\times \hat{\mathbf{b}},$$

$$\delta A_{\parallel }=A_{\parallel 0}+A_{\parallel 1},$$

$$\mathbf{B}=\hat{\mathbf{B}}+{\nabla }_{\perp }{A}_{\parallel 0}\times \hat{\mathbf{b}},$$

$$\frac{E_{\parallel }}{|{\mathbf{E}}_{\perp }|}\sim ϵ.$$

$${ϵ}_{\nu }=\frac{{\nu }_i}{{\Omega }_i}\sim {ϵ}^2.$$

$$\frac{{\nu }_e}{{\Omega }_i}\sim \sqrt{\frac{m_i}{m_e}}{\left(\frac{T_i}{T_e}\right)}^{3/2}{ϵ}_{\nu },$$

$$k_{\parallel }{\lambda }_{mfp}\sim \frac{ϵ{ϵ}_{\perp }}{{ϵ}_{\nu }}.$$

$$\gamma$$

$$\equiv \mathbf{\gamma }\cdot d\mathbf{z}={\gamma }_{\nu }d{\mathbf{z}}^{\nu }.$$

$$\mathbf{a}({\mathbf{x}}_0,{\theta }_0,t)$$

$$\mathbf{c}({\mathbf{x}}_0,{\theta }_0,t)$$

$$\mathbf{v}=\mathbf{U}+{\mathbf{c}}_{\perp }^{\prime },$$

$$\bm{Z}^{\nu}=\phi_{+}^{\nu}\left(\bm{z},\epsilon\right)=\sum_{n=0}^{\infty}\frac{\epsilon^{n}}{n!}\frac{\partial^{n}\phi_{+}^{\nu}\left(\bm{z},\epsilon\right)}{\partial\epsilon^{n}}\bigg{\rvert}_{\epsilon=0},$$

$${\mathbf{z}}^{\nu }={\varphi }_{-}^{\nu }(\mathbf{Z},ϵ).$$

$$\frac{\partial {\varphi }_{+}^{\nu }}{\partial ϵ}(\mathbf{z},ϵ)=g^{\nu }({\varphi }_{+}^{\nu }(\mathbf{z},ϵ)).$$

$$\frac{\partial {\varphi }_{-}^{\nu }}{\partial ϵ}=-g^{\lambda }\frac{\partial {\varphi }_{-}^{\nu }}{\partial {\mathbf{Z}}^{\lambda }},$$

$$F(\mathbf{Z},ϵ)=f(\mathbf{z}).$$

Full text

1 Introduction

The plasma periphery in a tokamak, extends from the external part of the closed flux surface region, typically referred to as the edge, to the scrape-off layer (SOL) region where the magnetic field lines intercept the machine vessel walls. These two regions are separated by the last closed flux surface (LCFS). Understanding the plasma dynamics in the periphery is necessary to address some the most crucial problems fusion is facing today, still undermining our capabilities to make reliable predictions of the performances of future tokamak, such as ITER (Shimada et al., 2007). Indeed, the plasma dynamics in the periphery largely controls the plasma heat exhaust, the refuelling, and helium removal. In addition, it governs the overall confinement performance of a tokamak. In fact, a low-to-high (L-H) confinement mode transition can be triggered if a heat power threshold is exceeded (Wagner, 2007). This results from the formation of a transport barrier, which reduces considerably the turbulent transport, while a pressure pedestal appears in the proximity of the LCFS (Gohil et al., 1994). This barrier is characterized by strong radial electric fields yielding strong $\bm{E}\times\bm{B}$ sheared flows (Burrell et al., 1994; Xia et al., 2006; Schirmer et al., 2006), which are thought to be responsible for the observed reduction of the H-mode turbulence level (Petty et al., 1998; Naulin, 2007). Once formed, the H-mode pedestal can periodically relax because of the presence of edge-localized modes (ELMs) releasing large amplitude bursts of particle and heat into the SOL (Zohm, 1996; Connor & Wilson, 2000; Kirk et al., 2005). ELM events are a major concern on the way to fusion energy.

The phenomena at play in the plasma periphery are mostly a consequence of the presence of electromagnetic fluctuations. These fluctuations determine the anomalous transport level (Stoneking et al., 1994; Scott, 2003; LaBombard et al., 2005; Snyder et al., 2012; Battaglia et al., 2014), the dynamics behind the L-H mode transition (Connor & Wilson, 2000), they might have a role in the physics of the density limit (Rogers & Drake, 1997; Rogers et al., 1998; LaBombard et al., 2001; Scott, 2003) and, furthermore, they are expected to limit the pedestal pressure gradient (Snyder et al., 2011; Dickinson et al., 2012) and to set the SOL width (Scott, 1997; D’Ippolito et al., 2002; Scott, 2003; Halpern et al., 2014; Mosetto et al., 2015; Militello & Omotani, 2016). Turbulent modes can develop into coherent filamentary structures with important perpendicular extension carrying a large amount of particles outwards in the SOL, strongly affecting the heat load on the plasma-facing components (Zweben et al., 2007; Agostini et al., 2011; Zweben et al., 2015; Garcia et al., 2015).These fluctuations occur on perpendicular scale lengths ranging from the particle Larmor radius, $\rho_{a}=v_{tha}/\Omega_{a}$, to a scale length of the order of the typical equilibrium (in the sense of time-averaged) total pressure and potential gradient scale lengths, namely $L_{P}$ and $L_{\phi}$ (Ritz et al., 1987; Endler et al., 1995; LaBombard et al., 2001; D’Ippolito et al., 2002; Garcia et al., 2007a; Nespoli et al., 2017). Here, $a$ is the species subscript (i.e., for electrons $a=e$ and for the ions $a=i$), $v_{tha}^{2}=2T_{a}/m_{a}$ is the particle thermal velocity, $T_{a}$ is the particle temperature, and $\Omega_{a}=q_{a}B/m_{a}$ is the particle gyrofrequency, with $q_{a}$ and $m_{a}$ the particle charge and mass, respectively. The turbulent dynamics is thought to be the result of the complex nonlinear development of electromagnetic ion-temperature-gradient (ITG) modes, drift waves, Kelvin-Helmholtz instabilities, and ballooning modes (Scott, 1997; Zeiler et al., 1998; Rogers & Dorland, 2005; Mosetto et al., 2013). While these modes occur on the $\rho_{i}$ and larger spatial scales (Mazzucato et al., 2008; Guttenfelder & Candy, 2011; Guttenfelder et al., 2013), there are evidences that also turbulent modes, occurring on the electron gyroscale with perpendicular wavenumber, $k_{\perp}$, satisfying $k_{\perp}\rho_{e}\sim 1$, can contribute significantly to the large scale heat flux and interact with modes at larger scales (Dorland et al., 2000; Jenko et al., 2000; Neiser et al., 2018). Moreover, finite Larmor radius (FLR) effects might be important in the description of the pedestal (Snyder et al., 2002). While the turbulent modes yield small amplitude fluctuations in the edge region (Brower et al., 1987; Ritz et al., 1987), order unity fluctuation levels can be found in the SOL region with $k_{\perp}L_{P}\sim 1$ (Ritz et al., 1987; Wootton et al., 1990; Garcia et al., 2007b; Xu et al., 2009; Zweben et al., 2015; Nespoli et al., 2017). In general, these turbulent modes display elongated structures in the direction parallel to the equilibrium magnetic field (Winslow et al., 1997; Zeiler et al., 1998; Thomsen et al., 2002; Halpern et al., 2013; Grulke et al., 2014), and have typical frequencies considerably smaller than $\Omega_{i}$ (Levinson et al., 1984; Ritz et al., 1987; Zweben et al., 2007).

Over the last decades, successful and significant progress has allowed important advances in the simulation of the turbulent plasma dynamics in the periphery. Both fluid (Dudson et al., 2009; Tamain et al., 2009; Ricci et al., 2012; Easy et al., 2014; Ricci, 2015; Halpern et al., 2016; Paruta et al., 2018) and gyrokinetic simulations (Xu et al., 2007; Cohen & Xu, 2008; Ku et al., 2009; Shi et al., 2015; Hakim et al., 2016; Pan et al., 2016; Chang et al., 2017; Pan et al., 2018) have been used to achieve this goal. Fluid models are usually based on the drift-reduced Braginskii equations (Braginskii, 1965; Zeiler et al., 1997; Simakov & Catto, 2005; Scott, 2007), which evolve the three-lowest order fluid moments of the kinetic equation, that are the particle density, velocity, and temperature. These fluid equations rely on high plasma collisionality as the fluid closure assumes that the particle distribution function is locally close to thermal equilibrium, and on the assumption that the perpendicular scale length of the fluctuations is large compared to $\rho_{i}$ (Hazeltine, 1973; Zeiler et al., 1997; Catto & Simakov, 2004). The assumption of local thermal equilibrium is hard to justify in the edge while it appears more reasonable in the SOL, where the particles temperatures are lower than in the core and at the top of the pedestal (Kocan et al., 2008). However, significant deviations from thermal equilibrium can occur in the SOL. These deviations from a Maxwellian distribution function were numerically investigated by kinetic simulations (Lönnroth et al., 2006; Tskhakaya, 2012; Battaglia et al., 2014), and are due to the presence of a tail of suprathermal particles, yielding important modifications to the parallel heat flux, of transient processes (e.g., ELM bursts influencing locally the particle collisionality), and of the sheath (Sigmar et al., 1996; Batishcheva et al., 1996; Batishchev et al., 1997; Kirk et al., 2005). These effects might be enhanced in typical ITER operating H-mode scenarios (Martin & Takizuka, 2008). To include kinetic modifications covering different range of collisional regimes in a fluid description and non-Maxwellian physics, an extension of fluid models based on higher order moments has been recently addressed by Jorge et al. (2017) within a drift-kinetic regime. This model considers potentially large electrostatic fluctuations by truncating an expansion of the electrostatic potential around the guiding-center position, i.e., the center of the particle gyration. Therefore, it is not able to properly describe fluctuations that occur on the particle gyroscale. Such fluctuations are particularly important in the description of anomalous transport in the edge region (Scott, 2003).

To describe fluctuations on the ion Larmor radius scale, gyrokinetic theories - pioneered by Catto (1978); Frieman & Chen (1982) - were successfully used in numerical and theoretical investigations of microturbulences (Furnish et al., 1999; Heikkinen et al., 2006; Brizard & Hahm, 2007; Idomura et al., 2009; Wang & Hahm, 2010; Hatch et al., 2011; Krommes, 2012). Most of the gyrokinetic models are derived within the standard gyrokinetic ordering (Dubin et al., 1983; Hahm, 1988; Parra & Calvo, 2014), i.e.

[TABLE]

and

[TABLE]

where $\phi_{1}$ is small amplitude fluctuating component of the electrostatic potential $\phi$, and $e=-q_{e}$. While allowing for $k_{\perp}\rho_{i}\sim 1$ and being appropriate for core plasma conditions, the ordering given in Eq. 1 breaks down in the presence of large scale and amplitude fluctuations, such as the ones present in the SOL. The analytical treatment of large scale and amplitude electromagnetic fluctuations was addressed by Dimits et al. (1992) in a generalized gyrokinetic ordering, assuming that $k_{\perp}\rho_{i}e(\phi-v_{\parallel}A_{\parallel})/T_{e}\ll 1$ where $A_{\parallel}$ is the parallel component of the fluctuating magnetic vector potential. This ordering is less restrictive than Eq. 1. Indeed, it allows for order unity fluctuations, $e\phi/T_{e}\sim 1$, on long scale lengths $k_{\perp}\rho_{i}\ll 1$, but also for $e\phi/T_{e}\ll 1$ on $k_{\perp}\rho_{i}\sim 1$ scales. Gyrokinetic theories have also been developed to retain the presence of strong $\bm{E}\times\bm{B}$ flows driven by large radial sheared electric fields in the internal and edge transport barrier regions (see, e.g., Bernstein & Catto, 1985; Brizard, 1995; Hahm, 1996; Hahm et al., 2009; Dimits, 2012). However, despite significant progress, the present gyrokinetic formulations do not provide yet an efficient framework for the description of the turbulent plasma periphery dynamics. The limitations of the present models are due to the lack of a proper description of arbitrarily far from equilibrium distribution functions in the presence of both large and small scale fluctuating electromagnetic fields and to the lack of a nonlinear gyrokinetic collision operator in a form that allows for an efficient numerical implementation, overcoming the extreme computational resources needed for the present edge gyrokinetic codes (Chang et al., 2017; Pan et al., 2018). In this regard, we note that gyrofluid models were developed to limit the computational cost of gyrokinetic simulations while overcoming the limitations of fluid models. These models aim to incorporate FLR and kinetic effects in a fluid-based description by taking a finite number of moments of the gyrokinetic equation (see, e.g., Brizard, 1992; Dorland & Hammett, 1993; Waltz et al., 1994; Beer & Hammett, 1996; Snyder & Hammett, 2001; Ribeiro & Scott, 2008; Madsen, 2013; Held et al., 2016). The truncation schemes for closing the fluid models including collisionless Landau damping and perpendicular coupling between moments of the distribution function (associated with parallel streaming, magnetic gradient drifts, and FLR effects) are usually obtained by techniques pioneered by Hammett & Perkins (1990); Hammett et al. (1992); Dorland & Hammett (1993) and by methods developed by Waltz et al. (1994); Beer & Hammett (1996). Unfortunately, the limitations of these closure schemes undermine the use of gyrofluid models to evolve the plasma periphery dynamics due to the wide range of collisionality and the possible large deviations from equilibrium. Recently, an extension of the gyrofluid model to an arbitrary number of moments has been obtained by Mandell et al. (2018) in an Hermite-Laguerre pseudo-spectral velocity formulation of delta-F gyrokinetic theory for electrostatic perturbations in core conditions. This model allows for a dynamical refinement of previous gyrofluid models with a tuneable accuracy, and provide an ideal framework to approach the development of a gyrokinetic model for the plasma periphery.

In the present paper, leveraging previous works (Hammett et al., 1993; Beer & Hammett, 1996; Sugama, 2000; Hahm et al., 2009; Zocco & Schekochihin, 2011; Madsen, 2013; Omotani et al., 2015; Zocco et al., 2015; Schekochihin et al., 2016; Loureiro et al., 2016; Jorge et al., 2017; Mandell et al., 2018), we derive a gyrokinetic model that we develop into a gyrokinetic moment (gyro-moment) hierarchy able to evolve the turbulent dynamics of the plasma in the periphery of tokamak devices. Our model allows to describe far from equilibrium distribution functions, linear and nonlinear FLR effects driven by small scale electromagnetic fluctuations, and effects associated with time-dependent background electromagnetic fields in the presence of a strong sheared radial electric field. To construct our model, we introduce the small expansion parameters $\epsilon=\rho_{s}/L_{P}\ll 1$ and $\epsilon_{\delta}=e\phi_{1}/T_{e}\ll 1$, with $\rho_{s}=c_{s}/\Omega_{i}$ the sound Larmor radius and $c_{s}^{2}=T_{e}/m_{i}$ the sound speed. Collisional effects are treated by assuming the plasma strongly magnetized, and by introducing the small parameter $\epsilon_{\nu}=\nu_{i}/\Omega_{i}\ll 1$ being $\nu_{i}\equiv\nu_{ii}$ the ion-ion collision frequency. The evolution equation of the distribution function, that is the gyrokinetic Boltzmann equation, is then obtained from fully nonlinear electromagnetic gyrokinetic equations of motion that are second order accurate in $\epsilon$ and $\epsilon_{\delta}$. These equations are derived within a perturbation theory known as Lie-transform perturbation theory, pioneered by Littlejohn (1982) and Cary & Littlejohn (1983). Within this formulation, two successive noncanonical phase-space coordinates transformations are performed to decouple and remove the gyrophase dependent part of the single particle motion from the slow drifts. We expand the full gyroaveraged distribution function onto a complete set of velocity-space Hermite-Laguerre polynomials. Projecting the gyrokinetic Boltzmann equation onto the Hermite-Laguerre basis yields an infinite set of coupled fluid equations describing the evolution of the expansion coefficients, termed gyro-moments. Collisional effects are introduced in our model by a nonlinear gyrokinetic Dougherty collision operator (Dougherty, 1964; Abel et al., 2008). Finally, the gyro-moment hierarchy is coupled to a set of gyrokinetic Maxwell’s equations derived self-consistently from a gyrokinetic variational principle. This set includes two coupled gyrokinetic Poisson’s equation and two coupled gyrokinetic Ampere’s laws. Within this framework, polarization and magnetization effects appear from the phase-space coordinate transformations carried out in the single particle dynamics. These polarization and magnetization terms are expressed as functions of gyro-moments, yielding complete and closed analytical expressions of the gyrokinetic Maxwell’s equations.

The rest of the paper is organized as follows. In Section 2, we order the main spatial and temporal scales at play in the plasma periphery and the relative amplitude of the electromagnetic fluctuations. In Section 3, we describe the single particle dynamics by deriving a set of second order accurate nonlinear electromagnetic gyrokinetic equations of motion by performing two phase-space coordinate transformations. In Section 4, we derive the gyrokinetic Boltzmann equation. The expansion of the distribution function into a Hermite-Laguerre polynomial basis and the derivation of the gyro-moment hierarchy is presented in Section 5, while the nonlinear gyrokinetic Dougherty collision operator is obtained and expanded in gyro-moments in Section 6. Then, the gyro-moment hierarchy is self-consistently coupled to a set of gyrokinetic Maxwell’s equations deduced from a gyrokinetic variational principle in Section 7. Finally, our main results are outlined in Section 8. In particular, the improvements over previous gyrokinetic theories are summarized.

2 Ordering Assumptions for the Plasma Periphery

In this section, we order the main spatial and temporal scales at play in the plasma periphery of tokamak devices. We first consider the scale lengths of the time-averaged profiles. We define the time-averaged perpendicular scale lengths of the total pressure $P=\sum_{a}n_{a}T_{a}$, with $n_{a}$ and $T_{a}$ the density and temperature of species $a$, and of the electrostatic potential $\phi$, namely $L_{P}$ and $L_{\phi}$, as $L_{P}\sim\lvert\nabla_{\perp}\ln\left<P\right>_{\tau}\rvert^{-1}$ and $L_{\phi}\sim\lvert\nabla_{\perp}\ln\left<\phi\right>_{\tau}\rvert^{-1}$, with $\nabla_{\perp}$ the gradient in the direction perpendicular to the magnetic field. The time-average operator, $\left<\cdot\right>_{\tau}$, acts over a time $\tau$ such that $\tau\omega\gg 1$, where $\omega$ is the typical frequency of the fluctuations, as defined by $\omega\sim\lvert\partial_{t}\ln n_{a}\rvert\sim\lvert\partial_{t}\ln T_{a}\rvert\sim\lvert\partial_{t}\ln\phi\rvert$.

The definitions of $L_{P}$ and $L_{\phi}$ allow us to introduce the small parameter $\epsilon$ as the ratio of the ion sound Larmor radius $\rho_{s}$ to the equilibrium total pressure scale length $L_{P}$ (Frieman & Chen, 1982),

[TABLE]

being $L_{P}\sim L_{\phi}$. We remark that $\rho_{s}\sim\rho_{i}$, whereas $\rho_{i}/\rho_{e}\sim\sqrt{m_{i}/m_{e}}\simeq 60$ for a deuterium plasma, assuming $T_{i}/T_{e}\sim 1$. Indeed, $1\lesssim T_{i}/T_{e}\lesssim 4$ in the plasma periphery (Kocan et al., 2008; Elmore et al., 2012). While an alternative small parameter, $\rho_{i}/L_{P}\ll 1$, can also be proposed (instead of $\rho_{s}$ in Eq. 3), the presence of turbulent modes in the plasma periphery with turbulent frequency $\omega/\Omega_{i}\sim\rho_{s}/L_{P}$ (Mosetto et al., 2013), that persist in the cold ion limit, motivates the definition in Eq. 3.

The spatial ordering in Eq. 3 is justified in many tokamaks and a wide range of experimental conditions. For example, in the edge region of the TEXT tokamak with plasma parameters $B=2$ T, $I_{p}=200$ kA, $n_{e}\simeq 3\times 10^{19}$ m*-3*, equilibrium scale lengths $L_{n}\sim L_{T_{e}}\sim L_{T_{i}}\sim L_{P}\sim 1.5$ cm and $\rho_{s}\sim\rho_{i}\simeq 0.02$ cm, one has $\epsilon\sim 0.0167$ (Ritz et al., 1987). For a medium size tokamak, such as TCV, with $B=1.5$ T, $T_{e}\sim T_{i}\sim 40$ eV and $n_{e}\simeq 6\times 10^{18}$ m3, it is possible to estimate $\epsilon\sim 0.043$ (Rossel et al., 2012; Halpern et al., 2016). For the edge transport barrier, typical values of $\epsilon\sim 1/30$ can be found in H-mode discharges in the JT60-U tokamak (Gohil et al., 1994). Additionally, Eq. 3 is also in agreement with detailed experimental measurements of the pedestal equilibrium scale lengths and scaling investigations of the H-mode pedestal structure (Burrell et al., 1994; Hubbard, 2000; Wagner, 2007; Zweben et al., 2015). In general, in H-mode experimental conditions of large aspect ratio tokamaks, the equilibrium pedestal pressure scale length is often found to be of the order of the ion poloidal gyroradius, $\rho_{\theta i}=\rho_{i}B_{T}/B_{\theta}$ (where $B_{T}$ and $B_{\theta}$ are the toroidal and poloidal components of the equilibrium magnetic field, respectively) (Snyder et al., 2009a), that is $L_{P}\sim\rho_{\theta i}$, yielding $\epsilon\sim B_{\theta}/B_{T}\sim 1/10$. Although the latter estimate may break down in the edge of low aspect ratio devices such as MAST and NSTX (Ono et al., 2000), the parameter $\epsilon$ remains small also in this kind of devices (Kirk et al., 2004; Militello et al., 2013). Finally, our ordering is valid in the operational conditions of ITER, as expected from predictions of its pedestal height and width (Snyder et al., 2009b; Beurskens et al., 2011).

The equilibrium magnetic scale length $L_{B}\sim\lvert\nabla_{\perp}\ln\hat{B}\rvert^{-1}$ of the equilibrium magnetic field, $\hat{B}$, is ordered by the small parameter

[TABLE]

since $L_{B}\sim R_{0}$, being $R_{0}$ the major radius of the tokamak device. We remark that, while one expects a maximal ordering $\epsilon_{B}\sim\epsilon$ for typical core conditions, in the plasma periphery where $L_{P}$ is steeper compared to the core, the scale length separation obeys $L_{P}/R_{0}<\rho_{s}/L_{P}$ (Burrell et al., 1994; Zweben et al., 2007) such that $\epsilon_{B}<\epsilon^{2}$, an assumption used in previous gyrokinetic edge model (see, e.g., Hahm et al., 2009).

To describe the perpendicular scale length of the fluctuations, we introduce the perpendicular wavenumber,

[TABLE]

and the parameter

[TABLE]

Our model addresses the presence of turbulent fluctuations on scale lengths ranging from $L_{P}$ to the $\rho_{e}$ scale. More precisely, we consider large scale $\epsilon_{\perp}\sim\epsilon$ fluctuations, typically present in the SOL, where they appear with $k_{\perp}L_{P}\sim 1$, and may have large amplitude (Zweben et al., 2015; Nespoli et al., 2017). At the same time, we include $\epsilon_{\perp}\gtrsim 1$ fluctuations (Wootton et al., 1990; Shats et al., 2005). These fluctuations have scale length of the order or smaller than $\rho_{s}$, and include electron gyroscale fluctuations ($k_{\perp}\rho_{e}\sim 1$), which are thought to be important in ETG driven turbulent transport experiments (Dorland et al., 2000; Colyer et al., 2017).

In order to describe the possible coexistence of electromagnetic fluctuations at amplitudes incompatible with the standard gyrokinetic ordering, in Eq. 1, and at scales ranging from the $L_{P}$ scale to $\rho_{s}$ (or smaller), we order the $\bm{E}\times\bm{B}$ drift, $\bm{u}_{E}=\bm{E}\times\bm{B}/B^{2}$, with $\bm{E}=-\nabla\phi$, to be small with respect to the sound speed $c_{s}$, by introducing the parameter,

[TABLE]

The ordering in Eq. 7 allows for fluctuations of the order of the thermal energy on large scale, i.e. $e\phi/T_{e}\sim 1$ for $\epsilon_{\perp}\sim\epsilon_{E}$ scales, and for small amplitude fluctuations, i.e. $e\phi/T_{e}\sim\epsilon_{E}$ for $\epsilon_{\perp}\sim 1$ scales. Also, it allows the presence of fluctuations occurring at intermediate scales between $L_{P}$ and $\rho_{s}$. While a sonic ordering, $\epsilon_{E}\sim 1$ yielding $e\phi/T_{e}\sim 1/\epsilon_{E}$ for $\epsilon_{\perp}\sim\epsilon_{E}$, has been considered in, e.g., Artun & Tang (1994); Brizard (1995); Hahm (1996); Qin et al. (2007), here we use the fact that fluctuations in the SOL satisfy $e\phi/T_{e}\sim 1$ on $\epsilon_{\perp}\sim\epsilon$ scales (Zweben et al., 2015), while in the edge $e\phi/T_{e}\sim\epsilon_{\delta}\sim\epsilon$ on $\epsilon_{\perp}\sim 1$ scales (Brower et al., 1987), to find from Eq. 7 that the $\epsilon_{E}$ and $\epsilon$ parameters can have similar size.

The $\bm{E}\times\bm{B}$ ordering in Eq. 7 is valid in both H- and L-mode, even when strongly equilibrium sheared radial electric fields and steep pressure gradients are present. In the L-H mode experimental conditions observed in the DIII-D tokamak, we observe $\bm{E}\times\bm{B}$ drift to the sound speed of the order of $|\bm{u}_{E}|/c_{s}\sim 1/10$, estimating $T_{i}\simeq 0.3$ KeV and $B\simeq 2$ T (Doyle et al., 1991; Gohil et al., 1994). Similar values can be inferred from measurements of the long wavelength $\bm{E}\times\bm{B}$ flows in the edge of the HT-$7$ tokamak (Xu et al., 2003). The ordering in Eq. 7 is also appropriate to describe the broadband turbulent power spectrum of the plasma periphery (Levinson et al., 1984; Ritz et al., 1987; Shats et al., 2005) and the corresponding fluctuation levels. Indeed, the turbulent dynamics in the SOL is typically characterized by coherent filamentary structures with perpendicular gradient scale lengths $L_{P}\sim L_{\phi}$ ($k_{\perp}\rho_{s}\sim 0.1$) and large amplitude fluctuations in the density and in the electrostatic potential, that are of the order of $T_{e}/e$, (Garcia et al., 2007a; Xu et al., 2009; Zweben et al., 2015; Nespoli et al., 2017). As observed in experimental characterizations of edge plasma turbulence (Ritz et al., 1987; Wootton et al., 1990; Fonck et al., 1993; Wolf, 2002; Wang et al., 2004; Xia et al., 2006), the amplitude of these large fluctuations present in the SOL decreases inside the LCFS to residual levels, while the turbulent power spectrum in the edge peaks at $k_{\perp}\rho_{i}\simeq 0.6$ with $\rho_{i}\sim\rho_{s}$ (Brower et al., 1987; Ritz et al., 1987). Thus, ion FLR effects might become non-negligible in this region (Snyder et al., 2002; Dickinson et al., 2012). On the other hand, we remark that the particle parallel streaming velocity, $v_{\parallel}$, can be comparable to the sound speed. Values of the dimensionless parallel Mach number of order of unity have been reported by experimental investigations of parallel flows (Wang et al., 2004; Pitts et al., 2007).

For our model to allow for fluctuations at scales from $L_{P}$ to $\rho_{s}$, we separate the potential $\phi$ into two components, i.e. a drift-kinetic and a gyrokinetic component, $\phi_{0}$ and $\phi_{1}$, respectively (see Section 3 for their rigorous definitions). Thus, we write (Dimits et al., 1992)

[TABLE]

and order the size of $\phi_{1}$ to be small, at all scales from $L_{P}$ to $\rho_{s}$, compared to the thermal energy by introducing the small parameter

[TABLE]

By imposing that both $\phi_{0}$ and $\phi_{1}$ provide a similar contribution to $\bm{u}_{E}$, we order their gradients such that $\lvert\nabla_{\perp}\phi_{0}\rvert\sim\lvert\nabla_{\perp}\phi_{1}\rvert$. Equations 7 and 9 imply $\lvert\rho_{s}\nabla_{\perp}\ln\phi_{1}\rvert\sim 1$ such that $\phi_{1}$ is allowed to vary on all scales (from $L_{P}$ to $\rho_{s}$), while $\phi_{0}$ is potentially large on large scale if $\epsilon_{\perp}\sim\epsilon$ (comparable to $L_{P}$) since $\epsilon_{\perp}e\phi_{0}/T_{e}\sim\epsilon_{E}\sim\epsilon$.

Following the critical balance conjecture (Goldreich & Sridhar, 1995; Schekochihin et al., 2008, 2009, 2016), we order the parallel scale length $k_{\parallel}\sim\lvert\nabla_{\parallel}\ln P\rvert\sim\lvert\nabla_{\parallel}\ln\phi\rvert$ by assuming that $c_{s}k_{\parallel}\sim k_{\perp}|\bm{u}_{E}|$ where $k_{\perp}|\bm{u}_{E}|$ is the frequency of the $\bm{E}\times\bm{B}$ flow. Following Eq. 7, this implies that $\rho_{s}k_{\parallel}\sim\epsilon_{\perp}\epsilon$, i.e.

[TABLE]

Equation 10 is compatible with experimental evidences of turbulent structures elongated in the parallel direction in the tokamak periphery and theoretical expectations of ballooning modes (Winslow et al., 1997; Bleuel et al., 2002; Garcia et al., 2007a; Naulin, 2007; Halpern et al., 2013; Grulke et al., 2014). We remark that Eqs. 4 and 10 implies that $k_{\parallel}L_{B}\sim 1/\epsilon$ on the $L_{P}$ scales. This encompasses, e.g., ballooning mode structures for which the parallel wavenumber is of the order of the connection length, i.e. $L_{c}\simeq 2\pi qR_{0}$ with $R_{0}\sim L_{B}$ (Halpern et al., 2013), and drift waves turbulent modes which can be characterized by large $k_{\parallel}$ values (Rogers & Ricci, 2010; Mosetto et al., 2013).

The typical turbulent frequency $\omega$ in the plasma periphery is much smaller than the ion gyrofrequency $\Omega_{i}$. More precisely, we assume that the turbulent frequency is comparable to the frequency of the $\bm{E}\times\bm{B}$ flow, i.e. $\omega\sim k_{\perp}|\bm{u}_{E}|$ (Dimits, 2012), yielding

[TABLE]

being $\omega/\Omega_{i}\sim\epsilon^{2}$ for large scale fluctuations (Zeiler et al., 1997; Simakov & Catto, 2005) while $\omega/\Omega_{i}\sim\epsilon$ on the $\epsilon_{\perp}\sim 1$ scale length (Dubin et al., 1983; Hahm, 1988; Hahm et al., 2009). As a result of the orderings in Eqs. 10 and 11, the parallel advection at the sound speed occurs on the same time scale as the turbulent fluctuations, i.e. $k_{\parallel}c_{s}\sim\omega$, in agreement with the critical balance that sets the parallel correlation length $L_{\parallel}\sim 1/k_{\parallel}$ of the perpendicular turbulent structure.

To describe magnetic fluctuations, we introduce a fluctuating magnetic vector potential $\delta\bm{A}$. We assume that the fluctuating magnetic field, $\delta\bm{B}=\nabla\times\delta\bm{A}$, is small compared to the large scale equilibrium magnetic field, $\hat{\bm{B}}=\nabla\times\hat{\bm{A}}$ where $\hat{\bm{A}}$ is the magnetic vector potential associated with the large scale equilibrium magnetic field. More precisely, we impose

[TABLE]

with $\hat{B}=|\hat{\bm{B}}|$, consistently with experimental measurements of edge magnetic fluctuations (Sötckel et al., 1999; Graessle et al., 1991). We remark that both $\hat{\bm{B}}$ and $\hat{\bm{A}}$ vary on the $L_{B}$ scale, and that the temporal scale of $\hat{\bm{A}}$ is ordered as

[TABLE]

This implies that $\hat{\bm{B}}$ evolves on time scales such that $|\partial_{t}\ln\hat{B}|/\Omega_{i}\sim\epsilon\epsilon_{B}$ (Hazeltine & Meiss, 2003; Green et al., 2003; Wagner, 2009). We remark that Eq. 13 allows us to retain the possible effects due to the change of the equilibrium magnetic field (e.g., for plasma shaping) that can occur on a time scale faster that the current diffusion time scale, but smaller that the typical frequency of the fluctuations. As a consequence of the ordering in Eq. 13, the inductive $\bm{E}\times\bm{B}$ velocity, i.e. $\bm{b}\times\partial_{t}\hat{\bm{A}}/B$, remains small compared to $c_{s}$ (see Eq. 7). The total fluctuating magnetic field can be written as

[TABLE]

where $\delta A_{\parallel}\equiv\hat{\bm{b}}\cdot\delta\bm{A}$ and $\delta\bm{A}_{\perp}\equiv\delta\bm{A}-\delta A_{\parallel}\hat{\bm{b}}$ being the parallel and perpendicular components of the fluctuating magnetic vector potential, respectively. Since compressional waves, associated with $\delta B_{\parallel}\equiv\hat{\bm{b}}\cdot\delta\bm{B}$, are much faster than the typical turbulent time scale considered on the scales ranging from $L_{P}$ to $\rho_{s}$ in low-$\beta$ plasmas (with $\beta=8\pi P/B^{2}$ the total thermal beta plasma) - as inferred from experimental measurements (Petty et al., 1998; Doyle et al., 2007; Saibene et al., 2007; Zweben et al., 2007)-, we systematically neglect $\delta B_{\parallel}$. Thus, we describe $\delta\bm{B}_{\perp}$, at the lowest order, by a fluctuating magnetic vector potential along the equilibrium magnetic field, i.e.

[TABLE]

where we neglect the parallel component, $\delta B_{\parallel}=\delta A_{\parallel}\hat{\bm{b}}\cdot\nabla\times\hat{\bm{b}}$ since it is of higher order and $\delta\bm{A}_{\perp}$. This removes completely the compressional Alfvén waves in our model. We remark that this approximation is consistent with low-$\beta$ plasma models currently implemented in global and edge gyrokinetic codes (see, e.g, Goerler et al. (2011); Lanti et al. (2019); Mandell et al. (2019)), and fulfills the plasma periphery conditions. We also show that our low-$\beta$ approximation is in fact compatible with the presence of steep equilibrium pressure gradients in Section 7.4. Finally, we deduce that $\phi\sim c_{s}\delta A_{\parallel}$ by comparing the $\bm{E}\times\bm{B}$ ordering in Eq. 7 with Eq. 12 (Brizard & Hahm, 2007; Dimits, 2012; Tronko et al., 2017a). This allows a general treatment of magnetic fluctuations in a maximal ordering.

Similarly to Eq. 8, the parallel component of the fluctuating magnetic potential, $\delta A_{\parallel}$, can be written as

[TABLE]

with $A_{\parallel 0}$ and $A_{\parallel 1}$ the drift-kinetic and gyrokinetic parts of the fluctuating magnetic potential, respectively. The time variation of $A_{\parallel 0}$ and $A_{\parallel 1}$ obeys Eq. 11. By imposing that both $A_{\parallel 0}$ and $A_{\parallel 1}$ provide a similar contribution to $\delta\bm{B}_{\perp}$, we order $\lvert\nabla_{\perp}A_{\parallel 0}\rvert\sim\lvert\nabla_{\perp}A_{\parallel 1}\rvert$. In addition, the ordering given in Eq. 9 with $\phi_{1}\sim c_{s}A_{\parallel 1}$ yields $c_{s}eA_{\parallel 1}/T_{e}\sim\epsilon_{\delta}$ such that $A_{\parallel 1}$ is small at all scales, while $A_{\parallel 0}$ can be potentially large at large scales since $\epsilon_{\perp}ec_{s}A_{\parallel 0}/T_{e}\sim\epsilon_{E}$. Thus, the magnetic field from $A_{\parallel 0}$ and the large scale component, $\hat{\bm{B}}$, is written as

[TABLE]

while the magnetic field generated by the small amplitude component $A_{\parallel 1}$ is $\bm{B}_{1}=\nabla_{\perp}A_{\parallel 1}\times\hat{\bm{b}}$, such the total magnetic field is given by $\bm{B}+\bm{B}_{1}$. We also remark that Eq. 11 gives the ordering of the inductive part of the parallel electric field, being $\partial_{t}A_{\parallel 0}\sim\epsilon^{2}\Omega_{i}A_{\parallel 0}$ and $\partial_{t}A_{\parallel 1}\sim\epsilon\Omega_{i}A_{\parallel 1}$ yielding

[TABLE]

In order to describe the different collisional regimes in the plasma periphery, we introduce the small parameter $\epsilon_{\nu}=\nu_{i}/\Omega_{i}\ll 1$ where $\nu_{i}\equiv\nu_{ii}$ is the ion-ion collision frequency. By imposing that the plasma remains strongly magnetized, we require that

[TABLE]

From Eq. 19, the electron-electron collision frequency $\nu_{e}\equiv\nu_{ee}$ is ordered as

[TABLE]

being $\nu_{e}\sim\sqrt{m_{i}/m_{e}}(T_{i}/T_{e})^{3/2}\nu_{i}$. In order to estimate the level of collisionality in the plasma periphery, we look at the ratio of the electron mean-free path, $\lambda_{mfp}=v_{the}/\nu_{e}$, to the parallel scale length $L_{\parallel}\sim 1/k_{\parallel}$. Equation 10 and Eq. 20 yield the estimate

[TABLE]

From Eq. 21, one infers that our model is compatible with high ($k_{\parallel}\lambda_{mfp}<1$) and low ($k_{\parallel}\lambda_{mfp}>1$) collisionality regime.

The ordering given in Eq. 19 is compatible with experimental conditions. Indeed, with typical TCV SOL parameters, estimating $B=1.4$ T, $n_{e}\simeq 1\times 10^{18}$ m*-3*, $T_{e}\sim T_{i}\sim 25$ eV, one finds $\epsilon_{\nu}\simeq 0.0003$ (Garcia et al., 2007b; Nespoli et al., 2017), and for JET discharges, estimating the SOL parameters by $B=2.4$ T, $n_{e}\simeq 10^{18}$ m*-3* and $T_{e}\sim T_{i}\sim 40$ eV, one obtains $\epsilon_{\nu}\simeq 8.3\times 10^{-5}$ (Erents et al., 2000; Xu et al., 2009). Lower values of $\epsilon_{\nu}$ are typical of the edge due the presence of ELMs, when temperatures are of the order of $T_{e}\sim T_{i}\sim 100$ eV, and $\epsilon_{\nu}\sim 10^{-6}$ (Pitts et al., 2003; Kirk et al., 2005).

3 Single Particle Dynamics in the Plasma Periphery

Following previous gyrokinetic models (Dimits et al., 1992; Qin et al., 2006; Hahm et al., 2009; Dimits, 2012; Madsen, 2013), we derive a model for the single particle dynamics in the plasma periphery by performing two successive changes of phase-space coordinates. Taking advantage of the low-frequency ordering in Eq. 11, these two transformations allow us to pass from the particle phase-space, described by the coordinates $(\bm{x},\bm{v},t)$ where $\bm{x}$ the particle position, $\bm{v}$ the particle velocity measured in the laboratory reference frame and $t$ the time coordinate, to a set of new coordinates where the fast gyromotion, associated with the gyrophase dependent part of the particle motion, is decoupled and removed from the gyrophase independent low-frequency (compared to $\Omega_{a}$) drifts. In addition, the transformations ensure the adiabatic invariance of one phase-space coordinate. The transformations we perform are obtained by using Lie-transform perturbation theory, a class of continuous change of phase-space coordinates that preserve the Hamiltonian structure of the original system. We remark that, since the coordinate systems are connected by an invertible transformation, the physics is not affected by the use of the new coordinates (Qin & Tang, 2004).

The present section is organized as follows. In Section 3.1, we describe the general methodology we follow to obtain the proper description of the single particle dynamics in the plasma periphery. The mathematical tools associated with Lie-transform perturbation theory are introduced in Section 3.2. Then, the first and second phase-space coordinate transformations are presented in Sections 3.3 and 3.4, respectively.

3.1 General Methodology

Within a coordinate independent formulation, the dynamics of a particle of species $a$ is prescribed by its action $\mathcal{A}_{a}=\int L_{a}dt=\int\gamma_{a}$, where $L_{a}$ is the particle Lagrangian and $\gamma_{a}=L_{a}dt$ is the Lagrangian one-form. In this section, we drop the particle subscript $a$ for simplicity. The Lagrangian one-form $\gamma$ is given, in the $(6+1)$-dimensional particle phase-space described by the coordinates $\bm{z}=(\bm{x},\bm{v},t)$, by (Jackson, 2012)

[TABLE]

where we use the Einstein’s notation and we denote $\bm{A}=\bm{A}(\bm{x},t)$ and $\phi=\phi(\bm{x},t)$ the magnetic vector potential and electrostatic potential evaluated at the particle position $\bm{x}$, and $\bm{z}^{\nu}$ is the $\nu$ component of $\bm{z}$. In general, the Greek index $\nu$ runs from $1$ to $6+1$ (since it includes the time coordinate $t$), whereas the Latin index $i$ runs from $1$ to $6$. Thus, $\bm{z}=\left(\bm{z}^{i},t\right)$ and $\bm{\gamma}=(\Lambda_{i},-\mathcal{H})$ with $\Lambda_{i}$ and $\mathcal{H}$ the symplectic and Hamiltonian components of the one-form $\gamma$, respectively. In particular, in Eq. 22, $\Lambda_{i}=qA_{i}+mv_{i}$ for $i=1,2,3$ with $A_{i}$ the $i$th component of the magnetic vector potential $\bm{A}$, while $\Lambda_{i}=0$ for $i=4,5,6$ and $\mathcal{H}=q\phi+mv^{2}/2$.

To decouple the fast gyromotion from to the low-frequency drifts in the particle dynamics, both contained in $\gamma$, two successive changes of coordinates are performed from the particle phase-space, $\bm{z}$. In the first coordinate transformation, we introduce the lowest order guiding-center coordinates, $\bm{Z}_{0}=(\bm{x}_{0},v_{\parallel 0},\mu_{0},\theta_{0},t)$. In this coordinate system, $\bm{x}_{0}$ is the zeroth-order guiding-center position, i.e. $\bm{x}_{0}\equiv\bm{x}$ corresponds to the particle position, and the coordinates $(v_{\parallel 0},\mu_{0},\theta_{0})$ are introduced as a first step to decouple the gyrophase dependent and independent parts of the particle motion. In particular, $v_{\parallel 0}=\bm{v}\cdot\bm{b}$ is the parallel velocity ($\bm{b}=\bm{B}/B$), $\theta_{0}$ is the lowest order gyrophase angle, and $\mu_{0}$ is the lowest order magnetic moment. More precisely, the coordinates $\theta_{0}$ and $\mu_{0}$ are defined by introducing the right-handed Frenet-Serret orthonormal vector basis $(\bm{e}_{1},\bm{e}_{2},\bm{b})$ with $(\bm{e}_{1},\bm{e}_{2})$ spanning the plane perpendicular to $\bm{b}$ (Frenet, 1852; Littlejohn, 1988). We also introduce the associated cylindrical vector basis $(\bm{a},\bm{b},\bm{c})$, where $\bm{a}\equiv\bm{a}(\bm{x}_{0},\theta_{0},t)$ and $\bm{c}\equiv\bm{c}(\bm{x}_{0},\theta_{0},t)$ are defined with respect the coordinate angle $\theta_{0}$, i.e.

[TABLE]

From Eqs. 23 and 24, we remark that $\bm{c}=\partial_{\theta_{0}}\bm{a}$. The $\theta_{0}$ coordinate is then explicitly introduced by writing the particle velocity as

[TABLE]

with $\bm{U}$ and $\bm{c}^{\prime}_{\perp}\equiv v_{\perp}\bm{c}$ being the gyrophase independent and dependent parts of the particle velocity, respectively (discussed more in the details in Section 3.3). Then, $v_{\perp}$, the perpendicular velocity measured in the frame moving with the velocity $\bm{U}$, is used to define the lowest order adiabatic magnetic moment $\mu_{0}=mv_{\perp}^{2}/(2B)$. Given the velocity decomposition in Eq. 25, the gyrophase dependent part of the velocity, $\bm{c}^{\prime}_{\perp}$, can be isolated in the particle one-form $\gamma$ in Eq. 22. This allows us to perform a perturbative change of coordinates from the $\bm{Z}_{0}$ coordinates to the guiding-center phase-space coordinates denoted by $\bm{Z}=\left(\bm{R},v_{\parallel},\mu,\theta,t\right)$. The guiding-center coordinates, $\bm{Z}$, are constructed such that the dynamics is described by a gyrophase independent guiding-center one-form $\Gamma=\Gamma(\bm{R},v_{\parallel},\mu,t)$, removing the gyrophase terms associated with background gradients. As a consequence, the guiding-center magnetic moment $\mu$ is dynamically conserved. We remark that the time coordinate $t$ does not change under the transformations. We perform this perturbative transformation up to second order in the small parameter $\epsilon$. We note that this is equivalent to carry out an expansion in the $m/q\sim\Omega^{-1}$ parameter, considered in previous litterature (Northrop, 1963; Kruskal, 1965; Littlejohn, 1979; Brizard, 1995; Hazeltine & Meiss, 2003). Then, we order the terms with respect to the parameter $\epsilon_{E}$ [see Eq. 3] up to second order to obtain a suitable guiding-center model for the plasma periphery, considering potentially large amplitude electromagnetic fluctuations present on the $L_{P}$ scale length (Dimits, 2012). We neglect the $\epsilon_{\delta}$ electromagnetic fluctuations in the guiding-center transformation. The second order accuracy is motivated by the importance of retaining the lowest order FLR corrections and the effects associated with a time-dependent $\bm{E}\times\bm{B}$ drift.

The second perturbative change of coordinates is from the guiding-center to the gyrocenter coordinates, denoted by $\overline{\bm{Z}}=\left(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},\overline{\theta},t\right)$. Within this change of coordinates, we consider the presence of small scale electromagnetic fluctuations, $\phi_{1}(\bm{x})$ and $A_{\parallel 1}(\bm{x})$, on the $\epsilon_{\perp}\sim 1$ scales, and perturb the guiding-center one-form. The conservation of the guiding-center adiabatic magnetic moment $\mu$, violated by the presence of the small scale fluctuations, is retrieved in the gyrocenter magnetic moment $\overline{\mu}$, up to second order in the small amplitude parameter $\epsilon_{\delta}$, by the perturbative construction of $\overline{\bm{Z}}$. As a result, we obtain the gyrophase independent gyrocenter one-form $\overline{\Gamma}=\overline{\Gamma}(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t)$ that we use to derive the gyrokinetic second order accurate electromagnetic equations of motion in the plasma periphery of tokamak devices.

3.2 Lie-Transform Perturbation Theory

We present the formalism we use to perform the perturbative coordinate transformations, i.e. a perturbation approach known as Lie-transform perturbation theory (Cary, 1981; Littlejohn, 1982; Cary & Littlejohn, 1983; Brizard & Mishchenko, 2009). This formalism let us pass from the coordinates $\bm{z}$ and the associated one-form $\gamma=\gamma_{\nu}d\bm{z}^{\nu}$ to a new set of coordinates $\bm{Z}$ with the associated new one-form $\Gamma=\Gamma_{\nu}d\bm{Z}^{\nu}$, the two one-forms, $\gamma$ and $\Gamma$, being linked by the $\bm{z}$ to $\bm{Z}$ transformation. We consider a near-identical coordinate transformation around the small parameter $\epsilon\ll 1$ in the form

[TABLE]

where $\phi_{+}^{\nu}=\phi_{+}^{\nu}\left(\bm{z},\epsilon\right)$ is the mapping function that specifies the coordinate transformation, such that $\phi_{+}^{\nu}(\bm{z},0)=\bm{z}^{\nu}$. In Eq. 26, the function $\phi_{+}^{\nu}$ transforms the coordinates $\bm{z}$ to the new coordinates $\bm{Z}$, given $\epsilon$. Indeed, the coordinates $\bm{Z}$ are the values of the function $\phi_{+}^{\nu}$ evaluated at $(\bm{z},\epsilon)$. Symmetrically, we can define the inverse transformation of Eq. 26 by introducing the mapping function $\phi_{-}^{\nu}(\bm{Z},\epsilon)$, such that

[TABLE]

Our perturbation theory is built on the framework given by Eqs. 26 and 27, and the Lie-transform is a special case of Eq. 26 where the function $\phi_{+}^{\nu}$ is specified by introducing a generating function, $g^{\nu}$, and asking that $\phi_{+}^{\nu}$ is solution of

[TABLE]

We remark that Eq. 28 is a functional relation since both sides are evaluated at $(\bm{z},\epsilon)$. An equation for $\phi_{-}^{\nu}$ can be obtained by taking the derivative with respect to $\epsilon$ on both sides of Eq. 27 with $\bm{Z}^{\nu}=\phi_{+}^{\nu}(\bm{z},\epsilon)$ and using Eq. 28 yielding

[TABLE]

being $d\bm{z}^{\nu}/d\epsilon=0$. We now deduce the transformation rule of scalar functions induced by a Lie-transform specified by Eq. 28. Let $f$ be a scalar function of the coordinates $\bm{z}$ and $F$ a scalar function of the new coordinates $\bm{Z}$. The scalar invariance requires that $f(\bm{z})=F(\bm{Z})$. More precisely, since the coordinate transformation in Eq. 26 depends explicitly on $\epsilon$, the function $F$ should also depend explicitly on $\epsilon$. Thus, we write

[TABLE]

Taking the derivative with respect to $\epsilon$ of Eq. 30, while noticing that $df/d\epsilon=0$, and using Eq. 28, we derive

[TABLE]

where we introduce the Lie-derivative $\mathcal{L}_{g}\equiv g^{\nu}\partial_{\nu}$ applied to scalar functions. Since Eq. 31 is a functional relation, it can be evaluated both at $\bm{z}$ and $\bm{Z}$. Indeed, the differential operator $\partial_{\nu}$ acting on $F$ is defined by

[TABLE]

Expanding $F(\bm{Z},\epsilon)$ around $\epsilon$, using Eq. 31 recursively and the fact that $F(\bm{z},0)=f(\bm{z})$, the functional relation between $F$ and $f$ under the transformation in Eq. 26 can be found, i.e.

[TABLE]

The inverse relation follows directly from Eq. 29, i.e.

[TABLE]

We emphasize again that Eqs. 33 and 34 are relations between functions, and, therefore, can be evaluated at both $\bm{z}$ and $\bm{Z}$. We also remark that $e^{\epsilon\mathcal{L}_{g}}$ indicates a linear differential operator. We refer to Eq. 33 as push-forward transformation, and to Eq. 34 as pull-back transformation (Brizard & Hahm, 2007; Brizard & Mishchenko, 2009).

Equations 33 and 34 allow us to derive the functional form of the coordinate transformation in Eq. 26, which is specified by Eq. 28. With the particular choice of scalar functions $F=\phi_{-}^{\nu}$ and $f=I^{\nu}$ [$I^{\nu}$ is the coordinate function, such that $I^{\nu}(\bm{z})=\bm{z}^{\nu}=\phi_{-}^{\nu}(\bm{Z},\epsilon)$], the push-forward transformation in Eq. 33 evaluated at $\bm{Z}$ yields

[TABLE]

The inverse coordinate transformation of Eq. 35 follows directly from the pull-back transformation in Eq. 34 with, in particular, $f=\phi_{+}^{\nu}$ and $F=I^{\nu}$ evaluated at $\bm{z}$, i.e.

[TABLE]

We now derive the transformation rule of a one-form [e.g., $\gamma$ in Eq. 22] under the coordinate transformation in Eq. 26. From the invariance $\Gamma_{\nu}d\bm{Z}^{\nu}=\gamma_{\nu}d\bm{z}^{\nu}$, the components of $\Gamma$ transform as components of a covariant vector,

[TABLE]

with $S$ a gauge function. The gauge function $S$ reflects the invariance of the action $\mathcal{A}=\int\Gamma$ under the addition of a total derivative. Evaluating the derivative with respect to $\epsilon$ on both sides of Eq. 37, using Eq. 29, and finally expanding $\Gamma_{\nu}(\bm{Z},\epsilon)$ around $\epsilon$, we find the following functional relation

[TABLE]

with $\mathcal{L}_{g}$ the Lie-derivative acting on a one-form $\gamma$ being defined by

[TABLE]

Emphasis is made here on the fact that the Lie-derivative in Eq. 39 is not equivalent to the one in Eq. 33, since they act on different mathematical objects.

When treating the single particle dynamics, a series of change of coordinates are performed in order to remove the gyrophase dependent parts of the particle dynamics at each $n$th-order in the expansion in the form of Eq. 34 (Dragt & Finn, 1976; Cary, 1981). In this case, the pull-back transformation in Eq. 34 can be carried out by using a perturbation approach and can be written as a composition of individual Lie-transforms, i.e.

[TABLE]

where we introduce $\mathcal{L}_{n}\equiv\mathcal{L}_{g_{n}}$ as a shorthand notation for the Lie-derivative associated with the generating function $g_{n}^{\nu}$. Similarly to Eq. 33, the inverse transformation of Eq. 40 is

[TABLE]

From Eq. 40 with $T_{\epsilon}=e^{\epsilon\mathcal{L}_{1}}e^{\epsilon^{2}\mathcal{L}_{2}}\dots$, the second order accurate pull-back transformation becomes

[TABLE]

In particular, with $f=\phi_{+}^{\nu}$ and $F=I^{\nu}$ evaluated at $\bm{z}$ in Eq. 42, the second order accurate coordinate transformation is

[TABLE]

At the same time, the hierarchy giving the functional relation between the one-form $\Gamma=\sum_{n}\Gamma_{n}$ and the one-form $\gamma=\sum_{n}\gamma_{n}$ is

[TABLE]

In Eq. 44, the Lie-derivatives act on one-forms and are, therefore, defined by the relation in Eq. 39.

In the following, we use the Lie-transform perturbation theory and solve the hierarchy in Eq. 44 to obtain successively the gyrophase independent guiding-center one-form $\Gamma(\bm{R},v_{\parallel},\mu,t)$ and the gyrocenter one-form $\overline{\Gamma}(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t)$ starting from the fundamental particle one-form $\gamma(\bm{x},\bm{v},t)$ in Eq. 22. The two transformations are systematically obtained up to second order in $\epsilon$ and $\epsilon_{\delta}$, respectively. The inherent degrees of freedom in choosing the generating functions $g^{\nu}_{n}$ allow for different expressions of $\Gamma$ and $\overline{\Gamma}$ depending on the desired properties of the one-forms (see, e.g., Miyato et al., 2009; Madsen, 2010).

3.3 Guiding-Center Transformation

In the first coordinate transformation $\bm{Z}=T_{\epsilon}\bm{Z}_{0}$ with $\bm{\bm{Z}}_{0}=(\bm{x}_{0},v_{\parallel 0},\mu_{0},\theta_{0},t)$ and $\bm{Z}=\left(\bm{R},v_{\parallel},\mu,\theta,t\right)$, we use the Lie-transform method up to second order in $\epsilon$ aiming to obtain the gyrophase independent guiding-center one-form $\Gamma=\Gamma(\bm{R},v_{\parallel},\mu,t)$. Within this transformation, we consider the presence of potentially large amplitude drift-kinetic electromagnetic fluctuations, i.e. $\phi_{0}$ and $A_{\parallel 0}$, while neglecting the gyrokinetic components $\phi_{1}$ and $A_{\parallel 1}$ following a two-step derivation (see, e.g, Brizard, 1995; Hahm, 1996; Hahm et al., 2009). We remark that the consistency of the ordering is maintained up to second order as discussed in Section 3.4. The present derivation follows closely the ones in Brizard (1995); Madsen (2010), but differs in the functional forms of the symplectic components of $\Gamma$ at second order [see Appendix A].

For the equations of motion to be gyrophase independent, we require that

[TABLE]

being $\mathcal{H}_{0}$ the guiding-center Hamiltonian, $\Gamma_{\bm{R}}$ and $\Gamma_{\parallel}$ the $\bm{R}$ and $v_{\parallel}$ symplectic components of the one-form $\Gamma$, respectively. Besides removing the gyrophase dependence from the equations of motion, we require the guiding-center magnetic moment, $\mu$, to be dynamically conserved, i.e. $\dot{\mu}=0$. The Euler-Lagrange equations of motion show that this can be obtained by imposing the following sufficient conditions in addition to Eq. 45,

[TABLE]

while $\partial_{\mu}\Gamma_{\theta}\neq 0$ in general. In the following, we choose the guiding-center generating functions such that the conditions in Eqs. 45 and 46 are satisfied at all the considered orders. We note that the physical time $t$ is not transformed, such that Eq. 43 imposes that $g_{1}^{t}=g_{2}^{t}=0$. We also note that the same conditions in Eqs. 45 and 46 are used in Miyato et al. (2009) and in Madsen (2010).

We first express $\gamma$ in Eq. 22 in the preliminary coordinates $\bm{Z}_{0}=(\bm{x}_{0},v_{\parallel 0},\mu_{0},\theta_{0},t)$ separating the gyrophase dependent and independent parts of the particle velocity, i.e. $\bm{U}$ and $\bm{c}^{\prime}_{\perp}$, respectively. Thus, we write $\bm{v}=\bm{U}+\bm{c}^{\prime}_{\perp}$ [see Eq. 25], imposing (Littlejohn, 1981; Brizard, 1995; Hahm, 1996; Qin et al., 2006; Hahm et al., 2009; Madsen, 2010; Jorge et al., 2017),

[TABLE]

With respect to standard gyrokinetic formalisms valid in the core, the inclusion of the $\bm{E}\times\bm{B}$ drift in $\bm{U}$ is due to the fact that, in the plasma periphery, strong time-dependent sheared radial electric fields can be present with a time-averaged perpendicular scale length $L_{\phi}\sim L_{P}$. The careful analysis of the role of these sheared electric fields is of special interest for stability and transport studies in the plasma periphery, especially when the L-H mode transition occurs. The Lagrangian one-form $\gamma$ in Eq. 22 can, therefore, be written as $\gamma=\gamma_{0}+\gamma_{1}$, where $\gamma_{1}\sim O(\epsilon\gamma_{0})$ (Northrop, 1963; Kruskal, 1965; Littlejohn, 1979; Hazeltine & Meiss, 2003), with

[TABLE]

being $\bm{A}=\hat{\bm{A}}+\hat{\bm{b}}A_{\parallel 0}$. We remark that, in writing Eq. 48, no ordering in the $\bm{E}\times\bm{B}$ drift is considered. This assumption is made at this point of the derivation in order to include the polarization drift along with the magnetic drifts (Wimmel, 1984) and the important FLR corrections in the background fields. We show then that the second order accurate guiding-center model is appropriate to describe situations when, in particular, $|\bm{u}_{E}|/c_{s}\sim\epsilon_{E}$, relevant for the plasma periphery. We also note that the terms in Eqs. 48a and 48b proportional to $m\bm{c}^{\prime}_{\perp}$ and $m\bm{U}\cdot\bm{c}^{\prime}_{\perp}$ are gyrophase dependent.

We now solve the hierarchy in Eq. 44 and Eq. 43 to obtain the guiding-center one-from $\Gamma$ and the guiding-center coordinates $\bm{Z}$, respectively. Using the fact that Eq. 44 are functional relations, we evaluate them at the guiding-center coordinates $\bm{Z}$, and specify later the analytical expressions of $\bm{Z}$ using Eq. 43. Choosing $S_{0}=0$, the zeroth-order guiding-centre one-form $\Gamma_{0}$ is

[TABLE]

Using the vectorial identity $\nabla\bm{A}-(\nabla\bm{A})^{T}\equiv\bm{\epsilon}\cdot\bm{B}$ with $\bm{\epsilon}$ the Levi-Cevita tensor, $(\bm{\epsilon}\cdot\bm{B})_{ij}=\epsilon_{ijk}B_{k}$, $\bm{B}=\nabla\times\bm{A}$ [see Eq. 17], and $\nabla=\nabla_{\perp}+\bm{b}\nabla_{\parallel}\equiv\partial/\partial\bm{R}\equiv\nabla_{\bm{R}}$, we derive the equations of motion $\dot{\bm{R}}=\bm{E}\times\bm{B}/B^{2}$ with $\bm{E}=-\nabla\phi_{0}-\partial_{t}\bm{A}$ from the Euler-Lagrange equations. This corresponds to the zeroth-order $\bm{E}\times\bm{B}$ drift of the guiding-center $\bm{R}$. Here, all the quantities are evaluated at $\bm{R}$, i.e. $\phi_{0}=\phi_{0}(\bm{R},t)$, $\hat{\bm{A}}=\hat{\bm{A}}(\bm{R},t)$, and $A_{\parallel 0}=A_{\parallel 0}(\bm{R},t)$, and notice that $\phi_{0}$ and $A_{\parallel 0}$ are velocity-independent. Using the Lie-derivative in Eq. 39 and Eq. 44b, we obtain the first order guiding-center correction,

[TABLE]

Choosing $S_{1}=0$ to remove one degree of freedom, we cancel the gyrophase dependent terms in Section 3.3, as required by Eq. 45, by choosing the first order generating function $\bm{g}^{\bm{R}}_{1}$ to be

[TABLE]

with $\bm{\rho}\equiv\bm{\rho}(\bm{R},\mu,\theta)$ the rotating $v_{\perp}$-dependent gyroradius vector [see Eq. 23]. By inserting Eq. 51 into Section 3.3, $\Gamma_{1}$ reduces to

[TABLE]

The calculation of the Lie-transform up to $O(\epsilon)$ is detailed in Appendix A. From Eq. 44c, one deduces that

[TABLE]

where $S=\bm{a}\cdot\partial_{t}\bm{c}=\bm{e}_{1}\cdot\partial_{t}\bm{e}_{2}$ and $\bm{T}=(\nabla\bm{c})\cdot\bm{a}=\nabla\bm{e}_{2}\cdot\bm{e}_{1}$ are the gyrogauge field vectors introduced by Littlejohn (1988). While these two quantities are gyrophase independent, they dependent on the relative orientation of the vector basis $(\bm{e}_{1},\bm{e}_{2},\bm{b})$. Indeed, performing a rotation of $(\bm{e}_{1},\bm{e}_{2},\bm{b})$ in a plane perpendicular to $\bm{b}$ by a shift $\alpha(\bm{R},t)$, such that $(\bm{e}_{1},\bm{e}_{2},\bm{b})\mapsto(\bm{e}^{\prime}_{1},\bm{e}^{\prime}_{2},\bm{b})$ linked by $\theta\mapsto\theta^{\prime}=\theta+\alpha(\bm{R},t)$, the $\bm{R}$ and $t$ components of an arbitrary one-from, $\gamma$, transform according to $\gamma_{\bm{R}}^{\prime}=\gamma_{\bm{R}}-\gamma_{\theta}\nabla\alpha$ and $\gamma_{t}^{\prime}=\gamma_{t}-\gamma_{\theta}\partial_{t}\alpha$, with $\nabla\alpha=\bm{T}^{\prime}-\bm{T}$ (and $\partial_{t}\alpha=S^{\prime}-S$). Here, $\gamma^{\prime}$ is the one-from defined with respect to $(\bm{e}^{\prime}_{1},\bm{e}^{\prime}_{2},\bm{b})$. From the invariance of the one-form, we deduce that $\gamma_{\bm{R}}=\overline{\gamma}_{\bm{R}}-\gamma_{\theta}\bm{T}$ (and $\gamma_{t}=\overline{\gamma}_{t}-\gamma_{\theta}S$) to ensure that $\gamma_{\bm{R}}^{\prime}=\overline{\gamma}_{\bm{R}}-\gamma_{\theta}\bm{T}^{\prime}$ (and $\gamma_{t}^{\prime}=\overline{\gamma}_{t}-\gamma_{\theta}S^{\prime}$) where $\overline{\gamma}_{\bm{R}}$ and $\overline{\gamma}_{t}$ are the gyrogauge invariant parts of $\gamma_{\bm{R}}$ and $\gamma_{t}$, respectively. While Littlejohn (1988) shows that one can make $\bm{T}$ vanish along a specific trajectory, it cannot, in general, be chosen to be zero over a finite region. However, Since the $\bm{T}$ and $S$ terms enters in the equations of motion as higher order terms, we neglect them and only keep the gyrogauge invariant parts of $\Gamma_{2}$ at second order that includes the lowest order FLR corrections in the fields, i.e.

[TABLE]

We notice that the last term in Eq. 53 contains the Ba$\tilde{\text{n}}$os drift, proportional to $\mu v_{\parallel}B\bm{b}\cdot\nabla\times\bm{b}/2\Omega$ (Banos, 1967). This term appears in the symplectic component of the guiding-center one-form derived by Madsen (2010), while it is neglected in the guiding-center one form in Hahm et al. (2009). Indeed, the Baños can be transferred from the symplectic to the Hamiltonian components by a proper choice of gauge function [see Appendix A]. Adding the guiding-center corrections $\Gamma_{1}$ and $\Gamma_{2}$, given in Section 3.3 and Eq. 53 respectively, to $\Gamma_{0}$ in Eq. 49, we obtain the $O(\epsilon^{2})$ accurate gyrophase independent guiding-center one-form, $\Gamma=\Gamma_{0}+\Gamma_{1}+\Gamma_{2}$, with $\Gamma_{1}\sim O(\epsilon\Gamma_{0})$ and $\Gamma_{2}\sim O(\epsilon\Gamma_{1})$, which can be written as

[TABLE]

Here, the velocity dependent effective vector potential, $\bm{A}^{*}$, is

[TABLE]

and the guiding-center Hamiltonian $\mathcal{H}_{0}$,

[TABLE]

We notice that $\Gamma$ in Eq. 55 is gyrophase independent thanks to guiding-center transformation. The first order accurate guiding-center coordinates $\bm{Z}=\left(\bm{R},v_{\parallel},\mu,\theta,t\right)$, which we derive from Eq. 43, are

[TABLE]

From the definition of the effective vector potential $\bm{A}^{*}$ in Eq. 56, we introduce an effective magnetic field $\bm{B}^{*}=\nabla\times\bm{A}^{*}$. The Jacobian of the guiding-center coordinate transformation $\mathcal{J}$, such that $d\bm{x}d\bm{v}=\mathcal{J}d\bm{R}dv_{\parallel}d\mu d\theta$ (Cary & Brizard, 2009), is given by $\mathcal{J}=B_{\parallel}^{*}/m$ with

[TABLE]

From Eq. 59, we derive the approximation $B_{\parallel}^{*}\simeq B+O(\epsilon^{2})$. We now order that the second order guiding-center one-form, $\Gamma$, is appropriate the describe situations where $|\bm{u}_{E}|/c_{s}\sim\epsilon_{E}$, that sets the proper amplitude and spatial scale of the fluctuating potential $\phi_{0}$ in the plasma periphery. In Eq. 57, we remark the presence of second order terms in $\epsilon_{E}$, $u_{E}^{2}\sim\epsilon_{E}^{2}c_{s}^{2}$ and $\mu B\bm{b}\cdot\nabla\times\bm{U}/(2\Omega)\sim\epsilon_{\perp}\epsilon_{E}q\phi_{0}$, which is a second order term with $\epsilon_{E}\sim\epsilon_{\perp}\sim\epsilon$. We remark that the $v_{\parallel}$ and $\bm{u}_{E}$ contributions to $\bm{U}$ [see Eq. 47] should be retained. Indeed, the $v_{\parallel}$ contribution, that is the Baños drift, is a second order term. This can be shown by using Eq. 17, i.e. $\bm{B}\simeq\hat{\bm{B}}+\nabla_{\perp}A_{\parallel 0}\times\hat{\bm{b}}$, and developing to find $\mu v_{\parallel}B\bm{b}\cdot\nabla\times\bm{b}/(2\Omega)\simeq\mu v_{\parallel}\hat{\bm{b}}\cdot[(\hat{\bm{b}}\cdot\nabla)\nabla_{\perp}A_{\parallel 0}-\hat{\bm{b}}\nabla_{\perp}^{2}A_{\parallel 0}]/(2\Omega)\sim\epsilon\epsilon_{E}q\phi_{0}$ with $\phi_{0}\sim c_{s}A_{\parallel 0}$. The $\bm{E}\times\bm{B}$ contribution is proportional to $(\mu B/2\Omega)\bm{b}\cdot\nabla\times\bm{u}_{E}\simeq\mu\nabla_{\perp}^{2}\phi_{0}/(2\Omega)$, and is also the order of $\epsilon\epsilon_{E}q\phi_{0}$. The latter term has a simple physical interpretation, being an FLR correction. Indeed, expanding the potential $\phi_{0}(\bm{x})$ at the particle position around the guiding-center position $\bm{R}$ (according the coordinate transformation $\bm{x}=\bm{R}+\bm{\rho}$), the averaged potential $\left<\phi_{0}(\bm{x})\right>\equiv\int_{0}^{2\pi}d\theta\phi_{0}(\bm{R}+\bm{\rho})/(2\pi)$ acting on the particle around its gyro-orbit is $\left<\phi_{0}(\bm{x})\right>=\phi_{0}(\bm{R})+(\mu/2q\Omega)\nabla_{\perp}^{2}\phi_{0}(\bm{R})+O(\epsilon^{3})$. A similar development can be made with the Baños term as being FLR correction in $A_{\parallel 0}(\bm{x})$. We also remark that the guiding-center one-form $\Gamma$ in Eq. 55 simplifies to the one used in Hahm et al. (2009) if the Baños term is neglected, and reduces to the guiding-center model derived in Jorge et al. (2017) if the FLR corrections in the fields are ignored. Finally, we note that the $\Gamma_{\theta}$ component of $\Gamma$, proportional to $\mu B/\Omega$, ensures $\dot{\mu}=0$, as shown by direct application of the Noether’s theorem (Cary & Brizard, 2009).

3.4 Gyrocenter Transformation

In the second coordinate transformation $\overline{\bm{Z}}=T_{\epsilon_{\delta}}\bm{Z}$, which maps the guiding-center coordinates $\bm{Z}=\left(\bm{R},v_{\parallel},\mu,\theta,t\right)$ to the gyrocenter coordinates $\overline{\bm{Z}}=\left(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},\overline{\theta},t\right)$, we perturb the guiding-center one-form $\Gamma$ in Eq. 55 by introducing small amplitude (of the order of $\epsilon_{\delta}$) gyrokinetic electromagnetic fluctuations, $\phi_{1}=\phi_{1}(\bm{x})$ and $A_{\parallel 1}=A_{\parallel 1}(\bm{x})$, that can vary at all scales, i.e. $\epsilon_{\perp}\sim 1$. Since such fluctuations act at the particle position $\bm{x}=\bm{R}+\bm{\rho}$, the conservation of the guiding-center magnetic moment, $\mu$, is broken by the $\theta$-dependence reintroduced by $\bm{\rho}$ contained in the spatial argument of $\phi_{1}$ and $A_{\parallel 1}$. To retrieve the dynamical conservation of $\mu$, we construct the gyrocenter coordinates, $\overline{\bm{Z}}$, that allows us to obtain the gyrophase independent gyrocenter one-form $\overline{\Gamma}=\overline{\Gamma}_{i}d\overline{\bm{Z}}^{i}-\overline{\mathcal{H}}dt$ such that the conditions in Eq. 46 are satisfied. The gyrocenter one-form $\overline{\Gamma}$ is derived by using Lie-transform perturbation theory around the small parameter $\epsilon_{\delta}$ up to second order. This allows us to describe the single particle dynamics in the plasma periphery in the presence of electromagnetic fluctuations at the particle gyroscale. We start by writing the guiding-center one-form $\Gamma$ in Eq. 55 in the gyrocenter coordinates $\overline{\bm{Z}}$ using the fact that Eq. 44 are functional relations. Then, we introduce the small-scale electromagnetic fluctuations in the guiding-center description obtaining the one-form $\Gamma+\delta\Gamma$ with $\delta\Gamma\sim O(\epsilon_{\delta}\Gamma)$ where $\delta\Gamma$ contains the contributions related to $\phi_{1}$ and $A_{\parallel 1}$, i.e.

[TABLE]

where we use that $\bm{x}\simeq\overline{\bm{R}}+\overline{\bm{\rho}}$ with $\overline{\bm{\rho}}=\bm{\rho}(\overline{\bm{R}},\overline{\mu},\overline{\theta})$ and neglect $\overline{\nabla}\overline{\bm{\rho}}\cdot d\overline{\bm{R}}$ since it is proportional to $\lvert\bm{\rho}\cdot\overline{\nabla}\ln\rho\rvert\sim|\overline{\bm{\rho}}\cdot\overline{\nabla}\ln\hat{B}|$ being of higher order. We take advantage of the degrees of freedom in the choice of the gyrocenter generating functions, denoted by $\overline{g}^{\nu}_{1}$ and $\overline{g}_{2}^{\nu}$, to impose that only the functional form of the guiding-center Hamiltonian, $\mathcal{H}_{0}$, is modified by the gyrocenter transformation. In this framework, the gyrocenter symplectic components $\overline{\Gamma}_{i}$ have the same functional form as $\Gamma_{i}$ in Eq. 55, but are evaluated at the gyrocenter coordinates instead, i.e. $\overline{\Gamma}_{i}\left(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t\right)=\Gamma_{i}\left(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t\right)$. This formulation, known as the Hamiltonian representation (Brizard & Hahm, 2007; Miyato & Scott, 2011), imposes that the symplectic components of $\overline{\Gamma}$ vanish at all $\epsilon_{\delta}^{n}$ orders with $n\geq 1$, while the gyrophase independence of $\overline{\Gamma}$ requires that $\partial_{\overline{\theta}}\overline{\mathcal{H}}=0$. The Hamiltonian formulation has the advantage that the guiding-center Jacobian, $B_{\parallel}^{*}/m$, is free from $\epsilon_{\delta}$ electromagnetic fluctuations, such that it preserves its guiding-center functional form, i.e. $B_{\parallel}^{*}d\bm{R}dv_{\parallel}d\mu d\theta/m=B_{\parallel}^{*}d\overline{\bm{R}}d\overline{v}_{\parallel}d\overline{\mu}d\overline{\theta}/m$. We remark that, compared with previous Hamiltonian formulation where the independent velocity variable is the parallel canonical momentum (Hahm, 1988; Brizard & Hahm, 2007), we use $\overline{v}_{\parallel}$ instead.

We now solve the hierarchy in Eq. 44 up to second order in $\epsilon_{\delta}$. From the zeroth-order transformation in Eq. 44a, we find $\overline{\Gamma}_{0}=\Gamma$ with $S_{0}=0$ and retrieve the guiding-center dynamics at the lowest order in $\epsilon_{\delta}$. The first order gyrocenter correction $\overline{\Gamma}_{1}$, given by Eq. 44b that is $\overline{\Gamma}_{1}=\overline{\Gamma}_{0}+\delta\Gamma$, is obtained by computing the Lie-derivative of $\overline{\Gamma}_{0}$, according to Eq. 39. This yields

[TABLE]

where $\overline{\nabla}\equiv\partial/\partial\overline{\bm{R}}$. Here, the overline notations $\overline{\bm{A}^{*}}$ and $\overline{\mathcal{H}}_{0}$ indicate that the guiding-center quantities are evaluated at $\left(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t\right)$, i.e. $\overline{\bm{A}^{*}}=\bm{A}^{*}(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t)$ with $\overline{\mathcal{H}}_{0}=\mathcal{H}_{0}(\overline{\bm{R}},\overline{v}_{\parallel},\overline{\mu},t)$ defined in Eq. 56 and in Eq. 57, respectively. The gyrophase dependent parts of the fluctuations can be isolated by introducing the gyroaverage operator, $\left<\chi\right>_{\overline{\bm{R}}}=\left<\chi\right>_{\overline{\bm{R}}}(\overline{\bm{R}},\overline{\mu},\overline{v}_{\parallel},t)$ acting on a function $\chi=\chi(\bm{x},\overline{\mu},\overline{v}_{\parallel},t)$, defined by

[TABLE]

with $\overline{\nabla}_{\perp}^{2}\equiv\overline{\nabla}\cdot\overline{\nabla}_{\perp}$. In Section 3.4, we performed the Taylor expansion around $\overline{\bm{R}}$ and evaluated the $\overline{\theta}$-integral at constant $\overline{\bm{R}}$. The gyroaverage operator allows us to separate the gyrophase dependent and independent parts of $\chi$, such that $\chi=\left<\chi\right>_{\overline{\bm{R}}}+\widetilde{\chi}$ where $\widetilde{\chi}$ is defined as the gyrophase dependent part of $\chi$, while $\left<\widetilde{\chi}\right>_{\overline{\bm{R}}}=0$ by construction. We can therefore write

[TABLE]

Imposing that the symplectic components in Eq. 61 vanish and that the Hamiltonian component remains gyrophase independent, one finds the first order gyrocenter generating functions with $\bm{A}_{1}=\hat{\bm{b}}A_{\parallel 1}$,

[TABLE]

where the first order gauge function, $S_{1}$, is solved iteratively. Indeed, by expanding $S_{1}=S_{10}+S_{11}+\dots$ with $S_{11}\sim O(\epsilon S_{10})$, we derive $S_{10}=q\overline{\widetilde{\Phi_{1}}}/\Omega$ (with $\overline{\chi}\equiv\int^{\overline{\theta}}d\overline{\theta}\chi)$ and $S_{11}=-d_{tgc}\overline{S_{10}}$. Here, we introduce the guiding-center convective derivative $d_{tgc}\equiv\partial_{t}+(\overline{v}_{\parallel}\bm{b}+\bm{D}_{\perp})\cdot\overline{\nabla}$ with $\bm{D}_{\perp}=\bm{b}\times\overline{\nabla}\>\overline{\mathcal{H}}_{0}/(qB)$, and the first order gyrokinetic potential $\Phi_{1}$, being defined as

[TABLE]

The details of the evaluation of $S_{1}$ are reported in Appendix B. With Eq. 64, the first order gyrocenter correction, $\overline{\Gamma}_{1}$ in Eq. 61, reduces to

[TABLE]

The second order perturbation analysis in $\epsilon_{\delta}$ is carried out by solving Eq. 44c. We report the details in Appendix B. The second order gyrocenter correction, $\overline{\Gamma}=-\overline{\mathcal{H}}_{2}dt$, is given by the second order gyrocenter Hamiltonian,

[TABLE]

where $S_{1}\simeq S_{10}+S_{11}$. In Section 3.4, we notice the presence of mixed terms, i.e., terms that are proportional to products between the large amplitude and small amplitude fluctuations. While it has been recognized that these terms are an important elements in the elaboration of a gyrokinetic models for the plasma periphery (Dimits, 2010), we notice that they are absent in Brizard (1995); Hahm (1988); Qin et al. (2007); Hahm et al. (2009). In these theories, since a large and time-independent potential yielding $\epsilon_{E}\sim 1$ is considered, such terms should be retained. However, within our ordering, these terms can be shown to be effectively smaller by a factor $\epsilon$ than the leading order one, i.e. the last term in Section 3.4 having imposed $S_{1}\simeq S_{10}$. Indeed, ordering the guiding-center convective derivative $d_{tgc}\sim\epsilon\Omega$ for $\epsilon_{\perp}\sim 1$ (since $\omega\sim k_{\perp}\lvert\bm{u}_{E}\rvert\sim c_{s}k_{\parallel}$) and $k_{\parallel}/k_{\perp}\sim\epsilon$ implies that the terms proportional to $d_{tgc}$, and $S_{11}$, are smaller by at least a factor $\epsilon$ compared to the leading order terms in Section 3.4. Thus, the leading order second order gyrocenter correction to $\overline{\Gamma}$ is given by

[TABLE]

where the second order gauge function, $S_{2}$, is solution of Appendix B. The corresponding second order gyrocenter generating functions are given by

[TABLE]

where we have neglected the higher order terms in Eq. 298. Evaluating the functional expressions in Eqs. 66 and B at the gyrocenter coordinates $\overline{\bm{Z}}$, we obtain the gyrocenter one-form $\overline{\Gamma}$, accurate up to $O(\epsilon^{2},\epsilon_{\delta}^{2})$,

[TABLE]

where the gyrokinetic Hamiltonian is

[TABLE]

with the second order gyrokinetic potential, $\left<\Psi_{1}\right>$, given by

[TABLE]

By using Eq. 43 with Eq. 64, the $O(\epsilon_{\delta})$ accurate gyrocenter coordinates $\overline{\bm{Z}}$ are obtained,

[TABLE]

The gyrokinetic potential $\Psi_{1}$ in Section 3.4 is evaluated at the particle position $\bm{x}$ expressed as a function of the gyrocenter coordinates $\overline{\bm{Z}}$, i.e. $\bm{x}(\overline{\bm{Z}})$. To express $\bm{x}$ in terms of the $\overline{\bm{Z}}$ coordinates, we proceed as follows. As a first step, we write the particle position $\bm{x}$ as a function of the $\bm{z}=(\bm{x},\bm{v},t)$ coordinates by introducing the coordinate function $I^{\bm{x}}(\bm{z})$ such that $I^{\bm{x}}(\bm{z})=\bm{x}$. To find the functional form of $I^{\bm{x}}$ in the gyrocenter phase-space, Eq. 41 is used with $f=I^{\bm{x}}$. We derive the guiding-center functional form of $I^{\bm{x}}$ that is $T_{\epsilon}^{-}I^{\bm{x}}(\bm{Z})=\bm{R}+\bm{\rho}+O(\epsilon^{2})$ [see Eq. 58]. As a second step, we consider again Eq. 41 using now the first order gyrocenter generating functions given in Eq. 64. We obtain the function $T_{\epsilon_{\delta}}^{-}T_{\epsilon}^{-}I^{\bm{x}}$, which gives the particle coordinate evaluated at the gyrocenter coordinates $\overline{\bm{Z}}$, i.e.

[TABLE]

with $\overline{\bm{\rho}}=\bm{\rho}(\overline{\bm{R}},\overline{\mu},\overline{\theta})$, being the function $\bm{\rho}$ defined in Eq. 51, and where we use the fact that $|\overline{\bm{g}}^{\bm{R}}_{1}\cdot\overline{\nabla}\overline{\bm{\rho}}/\overline{\rho}|\sim|\overline{\bm{\rho}}\cdot\overline{\nabla}\ln\hat{B}|\sim\epsilon_{B}$. At the leading order, it is sufficient to approximate

[TABLE]

in the argument of $\Psi_{1}$ since the other terms in Eq. 75 are higher order corrections. This leads to evaluate $\left<\Psi_{1}\right>_{\overline{\bm{R}}}\equiv\left<\Psi_{1}(\overline{\bm{R}}+\overline{\bm{\rho}})\right>_{\overline{\bm{R}}}$ consistently with Eq. 60 (Brizard, 1989; Sugama, 2000).

From the variation of the gyrocenter action $\mathcal{A}=\int\overline{\Gamma}$, we obtain the second order accurate electromagnetic gyrokinetic equations of motion,

[TABLE]

where $\overline{v}_{\parallel}=\bm{b}\cdot\dot{\overline{\bm{R}}}+\overline{\mu}_{\parallel}$ (where we introduce $\overline{\mu}_{\parallel}=\overline{\mu}B(\bm{b}\cdot\overline{\nabla}\times\bm{b})/(2m\Omega)$), $\dot{\overline{\mu}}=0$. The $\dot{\overline{\bm{R}}}$ and $\dot{\overline{v}_{\parallel}}$ equations of motion can be obtained by taking the vector and scalar products of Eq. 76 with $\bm{b}$ and $\bm{B}^{*}$, respectively. Using the definition of $B_{\parallel}^{*}$ in Eq. 59 and the fact that $\overline{v}_{\parallel}=\bm{b}\cdot\dot{\overline{\bm{R}}}+\overline{\mu}_{\parallel}$, we derive the following gyrocenter equations of motion,

[TABLE]

where the convective derivative

[TABLE]

is evaluated with the gyrophase independent particle velocity $\overline{\bm{U}}=\bm{u}_{E}+\overline{v}_{\parallel}\bm{b}$. We remark that the parallel electric field, $E_{\parallel}\equiv\bm{E}\cdot\bm{b}=-\overline{\nabla}_{\parallel}\phi_{0}-(\partial_{t}\bm{A})\cdot\bm{b}$, contains both its electrostatic and inductive parts. In particular, the parallel inductive part is $(\partial_{t}\bm{A})\cdot\bm{b}=(\partial_{t}\hat{\bm{A}})\cdot\bm{b}+[\partial_{t}(A_{\parallel 0}\hat{\bm{b}})]\cdot\bm{b}$, where $\partial_{t}\hat{\bm{A}}$ obeys Eq. 13 while $\partial_{t}A_{\parallel 0}$ follows Eq. 11. We notice also the presence of the ratio of the magnetic field strength $B$ and of the effective magnetic field $\bm{B}^{*}$ to the gyrocenter phase-space volume element $B_{\parallel}^{*}$ in the above equations of motion

Equation 77, which describes the motion of a single gyrocenter in the plasma periphery, includes the polarization drift $\bm{b}\times d_{t}\overline{\bm{U}}/\Omega_{a}$, the magnetic gradient drifts, e.g. $\overline{\mu}\bm{b}\times\overline{\nabla}B/\Omega_{a}$, and, finally, a number of transport terms driven by the fluctuations at the particle Larmor radius scale contained in the $\bm{b}\times\overline{\nabla}\left<\Psi_{1}\right>_{\overline{\bm{R}}}/B$ term [see Section 3.4]. In particular, these terms are the perturbed electrostatic $\bm{E}\times\bm{B}$ drift, proportional to $\bm{B}\times\overline{\nabla}_{\perp}\left<\phi_{1}\right>_{\overline{\bm{R}}}/B^{2}$, the shear-Alfvén transport term, proportional to $\overline{v}_{\parallel}\bm{b}\times\overline{\nabla}_{\perp}\left<A_{\parallel 1}\right>_{\overline{\bm{R}}}/B$, also referred to as the magnetic-flutter velocity (Brizard & Hahm, 2007; Hahm et al., 2009), and, finally, a nonlinear electromagnetic term proportional to $\bm{b}\times\overline{\nabla}_{\perp}\partial_{\overline{\mu}}\left<\widetilde{\Phi_{1}}^{2}\right>/B$ that drives ponderomotive effects (Brizard & Hahm, 2007; Hahm et al., 2009; Krommes, 2012). Equation 78 is the parallel momentum equation that includes the parallel forces associated with the parallel electric field $qE_{\parallel}$, mirror force $\overline{\mu}\overline{\nabla}_{\parallel}B$, and a FLR induced parallel force driven by the gyrokinetic potential proportional to $q\overline{\nabla}_{\parallel}\left<\Psi_{1}\right>_{\overline{\bm{R}}}$. The dynamical conservation of $\overline{\mu}$ is given by Eq. 79, whereas Eq. 80 represents the evolution in time of the gyrocenter gyrophase $\overline{\theta}$, which differs from the guiding-center gyroangle $\theta$ due to torsional effects driven by the small-scale electromagnetic perturbations and the FLR field corrections.

We remark that the equations of motion in Eqs. 77 and 78 contain terms proportional to the gradient of the magnetic field strength, being $B=\lvert\bm{B}\rvert$ with $\bm{B}\simeq\hat{\bm{B}}+\overline{\nabla}_{\perp}A_{\parallel 0}\times\bm{b}$. This yields a polarization drift, $\bm{b}\times d_{t}\bm{U}/\Omega\sim\epsilon^{2}\epsilon_{E}c_{s}$, that can be of the same order as the $\overline{\nabla}B$ drift on the $L_{B}$ scale, while the $\overline{\nabla}B$ drift driven by perpendicular magnetic perturbations, $\delta\bm{B}_{\perp}\simeq\overline{\nabla}_{\perp}A_{\parallel 0}\times\bm{b}$, is of the order of $\epsilon^{2}c_{s}$.

Neglecting the field FLR corrections, proportional to $\bm{b}\cdot\overline{\nabla}\times\bm{U}/\Omega$, the gyrokinetic potential $\left<\Psi_{1}\right>$ and the inductive contribution $(\partial_{t}\bm{A})\cdot\bm{b}$, Eqs. 77 and 78 reduce to the guiding-center equations of motion presented in Jorge et al. (2017) and obtained from a direct gyroaveraging of the particle Lagrangian, with the fields formally expanded around the particle position $\bm{x}$. In addition, the equations of motion in Eqs. 77 and 78 constitute an improvement over previous gyrokinetic theories for the edge region. In particular, they are the generalization of the equations developed by Hahm et al. (2009); Dimits (2012); Madsen (2013) to consider time dependent large amplitude and scale electromagnetic fields. They also generalize the models derived by Dimits et al. (1992) and Qin et al. (2006) by describing electromagnetic fluctuations at second order.

4 Gyrokinetic Boltzmann Equation

Having derived the equations of motion of a single gyrocenter in the presence fluctuating electromagnetic fields in the plasma periphery, we can now address their collective dynamics. In the present section, we therefore express the gyrokinetic Boltzmann equation in the gyrocenter phase-space coordinates.

The distribution function $f_{a}(\bm{x},\bm{v},t)$ of particles species $a$ (hereafter, we reintroduce the species subscript $a$) obeys the Boltzmann kinetic equation,

[TABLE]

with $C_{a}(f_{a})$ being the collision operator (see Section 6). To write Eq. 82 in the gyrocenter phase-space coordinates, $\overline{\bm{Z}}$, derived in Section 3.4, we introduce the full gyrocenter distribution function $\overline{F_{a}}(\overline{\bm{Z}})=f_{a}(\bm{x}(\overline{\bm{Z}}),\bm{v}(\overline{\bm{Z}}))$ (the time coordinate $t$ is omitted for simplicity), and use the chain rule to express $(\bm{x},\bm{v})$ in terms of the $\overline{\bm{Z}}$ coordinates, that is

[TABLE]

where the dynamical conservation of the magnetic moment, $\dot{\overline{\mu}}=0$, is used. The gyrocenter equations of motion $\dot{\overline{\bm{R}}}$, $\dot{\overline{v}_{\parallel}}$ and $\dot{\overline{\theta}}$ are given by Eqs. 77, 78 and 80.

In Eq. 83, the full gyrocenter distribution function $\overline{F_{a}}=\overline{F_{a}}(\overline{\bm{Z}})$ can be written as $\overline{F_{a}}=\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}+\widetilde{\overline{F_{a}}}$, where $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}\equiv\left<\overline{F_{a}}(\overline{\bm{Z}})\right>_{\overline{\bm{R}}}=\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\overline{\bm{R}},\overline{\mu},\overline{v}_{\parallel})$ and $\widetilde{\overline{F_{a}}}=\widetilde{\overline{F_{a}}}(\overline{\bm{Z}})$ are the gyrophase independent and dependent parts of $\overline{F_{a}}$, respectively. The operator $\left<\cdot\right>_{\overline{\bm{R}}}$ is the gyroaverage operator defined in Section 3.4. We remark that no assumption is made of the spatial variation of the gyroaveraged distribution function, i.e. $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ is allowed to vary on both $\epsilon_{\perp}\sim 1$ and $\epsilon_{\perp}\sim\epsilon$ scales. The evolution equation of $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ can be obtained by applying the gyroaverage operator to Eq. 83 yielding

[TABLE]

where we use the gyrophase independence of $\dot{\overline{\bm{R}}}$ and $\dot{\overline{v}}_{\parallel}$. Using the fact that the gyrocenter phase-space volume element, $B_{\parallel}^{*}/m_{a}$ in Eq. 59, is conserved along the gyrocenter trajectories (Brizard & Hahm, 2007), i.e.

[TABLE]

Eq. 84 can be written in a conservative form,

[TABLE]

Equation 86 is a formulation more convenient than Eq. 84 to derive the gyro-moment hierarchy equation (see Section 5). Subtracting Eq. 84 to the gyrokinetic Boltzmann equation Eq. 83, one finds the evolution equation of $\widetilde{\overline{F_{a}}}$,

[TABLE]

We notice that $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ and $\widetilde{\overline{F_{a}}}$ are coupled through the collision operator $C_{a}$ in Eqs. 84 and 87, since $C_{a}$ acts on the full gyrocenter distribution function. Aiming to obtain a closed gyrophase independent evolution equation for $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$, we estimate the magnitude of $\widetilde{\overline{F_{a}}}$ with respect to $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$. This can be done by comparing the leading order terms of Eq. 87 at the left- and right-hand sides. At the left-hand side, the leading term is $\dot{\overline{\theta}}\partial\widetilde{\overline{F_{a}}}/\partial\overline{\theta}$, since $\dot{\overline{\theta}}\partial/\partial\overline{\theta}\sim\Omega_{a}$ while $\partial/\partial t\sim\dot{\overline{v}_{\parallel}}\partial/\partial\overline{v}_{\parallel}\sim\dot{\overline{\bm{R}}}\cdot\overline{\nabla}\sim\epsilon\epsilon_{\perp}\Omega_{i}$. At the right-hand side, expanding the full distribution function $\overline{F_{a}}$ in the collision operator $C_{a}$ as $\overline{F_{a}}=\overline{F_{a0}}+\epsilon_{\nu}\overline{F_{a1}}+\epsilon_{\nu}^{2}\overline{F_{a2}}+\dots$ with $\overline{F_{a0}}=\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ and using the ordering of the collision frequencies Eqs. 19 and 20, such that $C_{e}(\overline{F_{e}})\sim\sqrt{m_{i}/m_{e}}(T_{i}/T_{e})^{3/2}\epsilon_{\nu}\Omega_{i}\overline{F_{e}}$ while $C_{i}(\overline{F_{i}})\sim\epsilon_{\nu}\Omega_{i}\overline{F_{i}}$, we obtain the following estimate for the electrons,

[TABLE]

and for the ions,

[TABLE]

As a consequence, up to second order in $\epsilon$ (or $\epsilon_{\delta}$), the gyrophase dependent part $\widetilde{\overline{F_{a}}}$ can be neglected in the gyrokinetic Boltzmann equation in Eq. 86 for both electrons and ions. Further approximations can be done by noticing that the collision operator $C_{a}$, associated with the gyrocenter transformation [see Section 3.4], can be expanded in powers of $\epsilon_{\delta}$, such that $C_{a}(F_{a})=C_{a0}(F_{a})+\epsilon_{\delta}C_{a1}(F_{a})+\dots$ with $\epsilon_{\delta}\sim\epsilon$. Because of the ordering of the collision frequencies in Eqs. 19 and 20 and up to second order in $\epsilon$ (or $\epsilon_{\delta}$), we can neglect the gyrokinetic corrections $C_{a1}(F_{a})$ and higher, since they are $O(\epsilon_{\nu}\epsilon_{\delta})$. We remark that, while neglecting the $\epsilon_{\delta}$ corrections in $C_{a}$, the FLR effects contained in $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ are retained at arbitrary order in $\epsilon_{\perp}$ in the collision operator. To conclude, the gyrokinetic Boltzmann equation in Eq. 86, second order accurate in $\epsilon$ and $\epsilon_{\delta}$, is

[TABLE]

which is a closed equation for $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$.

5 Gyro-Moment Hierarchy

In this section, we address the development of a gyrokinetic moment hierarchy, referred to as gyro-moment hierarchy, that we propose as a technique to evolve the gyrokinetic Boltzmann equation, Eq. 90, that is derived in Section 4. In the present section, we focus on the collisionless part of the gyrokinetic Boltzmann equation, while the collisional part is the subject of Section 6. In Section 5.1, we introduce a velocity-space Hermite-Laguerre decomposition of the gyroaveraged distribution function and relate the coefficients of this expansion to fluid-like quantities. Section 5.2 describes the gyro-moment hierarchy equation that sets their evolution. The hierarchy is obtained by projecting the gyrokinetic Boltzmann equation onto the Hermite-Laguerre basis. In the process, we retain the parallel and perpendicular phase-mixing terms, i.e. the coupling terms between gyro-moments arising from the parallel and perpendicular drifts and FLR corrections. As an example of the evaluation of the gyro-moment hierarchy equation, we derive the lowest order equations that describe the evolution of the fluid moments in Section 5.3. Finally, the gyro-moment hierarchy equation is completed in Section 5.4 by providing a gyro-moment expansion of the terms that contain the gyroaveraged gyrokinetic potential.

5.1 Gyro-Moment Expansion

As a technique to solve the gyrokinetic Boltzmann equation, Eq. 90, we derive a gyro-moment hierarchy based on decomposing the full gyroaveraged gyrocenter distribution function $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ onto a complete velocity-space basis provided by the Hermite-Laguerre polynomials. Therefore, we write

[TABLE]

where $\overline{N}_{a}^{lk}=\overline{N}_{a}^{lk}(\overline{\bm{R}},t)$ are the expansion coefficients of $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$, hereafter referred to as gyro-moments, and $H_{a}^{lk}$ are the Hermite-Laguerre basis elements defined by

[TABLE]

Here, $H_{l}$ denotes the physicits’ Hermite polynomials of order $l$ defined by the Rodrigues’ formula

[TABLE]

while $L_{k}$ are the Laguerre polynomials given by the Rodrigues’ formula (Abramowitz, 1974)

[TABLE]

The Hermite polynomials $H_{l}$ are orthogonal over the interval $]-\infty,\infty[$ weighted by $e^{-x^{2}}$ such that

[TABLE]

with $\delta_{l^{\prime}}^{l}$ the Kronecker delta, whereas the Laguerre polynomials $L_{k}$ are orthogonal over the interval $[0,+\infty[$ weighted by $e^{-x}$ via the relation

[TABLE]

In Eq. 91, we also introduce the gyrocenter shifted Maxwellian distribution function $\overline{F}_{Ma}$ by

[TABLE]

with $\overline{N}_{a}=\overline{N}_{a}(\overline{\bm{R}},t)$ the gyrocenter density defined as $\overline{N}_{a}=\int d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}B\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}/m_{a}$. The overline notation indicates that the fluid quantities are those associated with the gyrocenters, i.e. defined as moments of gyroaveraged gyrocenter distribution $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$. The velocity variables $\overline{s}_{\parallel a}=(\overline{v}_{\parallel}-\overline{u}_{\parallel a})/\overline{v}_{th\parallel a}$ and $\overline{s}_{\perp a}^{2}=\overline{v}_{\perp}^{2}/\overline{v}_{th\perp a}^{2}=\overline{\mu}B/\overline{T}_{\perp a}$ represent the normalized shifted parallel and perpendicular gyrocenter velocity. Indeed, $\overline{u}_{\parallel a}$ is the gyrocenter parallel fluid velocity, while $\overline{v}_{th\parallel a}^{2}=2\overline{T}_{\parallel a}\ /m_{a}$ and $\overline{v}_{th\perp a}^{2}=2\overline{T}_{\perp a}/m_{a}$ are, respectively, the parallel and perpendicular thermal velocities associated with the parallel and perpendicular temperatures $\overline{T}_{\parallel a}$ and $\overline{T}_{\perp a}$. These are defined by $\overline{N}_{a}\overline{T}_{\perp a}=\int d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}B\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}\overline{\mu}B/m_{a}$ and $\overline{N}_{a}\overline{T}_{\parallel a}=\int d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}B\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(v_{\parallel}-\overline{u}_{\parallel a})^{2}$, respectively, with the parallel fluid velocity $\overline{N}_{a}\overline{u}_{\parallel a}=\int d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}B\overline{v}_{\parallel}\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}/m_{a}$. Our definitions of the perpendicular and parallel thermal speeds motivate the choice of the physicists’ Hermite polynomials as basis, being orthogonal with respect to a Maxwellian distribution, instead of the probabilists’ Hermite polynomials, used, e.g., in Mandell et al. (2018), that are orthogonal to a Gaussian function of the form $e^{-x^{2}/2}$.

Since the Hermite-Laguerre polynomials define an orthogonal basis of the space of functions $f$ such that (Wong, 1998)

[TABLE]

then, any distribution function $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ that fulfills the requirement given in Eq. 98 can be decomposed as in Eq. 91. We remark that the number gyro-moments that must be kept in the Hermite-Laguerre decomposition in order to provide a good pseudo-spectral representation is directly related to the deviation of the distribution function from Maxwellian. In general, a large number of coefficients is needed if fine velocity structures are present (Schekochihin et al., 2016). We note that the choice of $\overline{s}_{\parallel a}$ and $\overline{s}_{\perp a}$ as arguments of the basis functions [see Eq. 92], provides an efficient representation in both the strong ($\overline{u}_{\parallel a}/\overline{v}_{th\parallel a}\sim 1$) and the weak ($\overline{u}_{\parallel a}/\overline{v}_{th\parallel a}\ll 1$) flow regimes, ensuring that the Hermite-Laguerre decomposition, in Eq. 91, is also efficient in the case of distribution functions that do not deviate significantly from a Maxwellian with finite $\overline{u}_{\parallel a}$ (Jorge et al., 2017).

Using the orthogonality relations in Eqs. 95 and 96, the gyro-moments $\overline{N}_{a}^{lk}$ are evaluated as

[TABLE]

and correspond, indeed, to generalized moments of the full gyroaveraged gyrocenter distribution function $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$, i.e. to fluid-like quantities. To conveniently derive the gyro-moment hierarchy, which describes the spatial and time evolution of $\overline{N}_{a}^{lk}$, we introduce the Hermite-Laguerre projector operator of order $(l,k)$, $\left\lVert\cdot\right\rVert^{lk}_{a}$ applied to a phase-space function $\chi=\chi(\overline{\bm{Z}})$,

[TABLE]

To simplify the notation, we denote the $(l,k)=(0,0)$ Hermite-Laguerre projector operator simply by $\left\lVert\cdot\right\rVert_{a}$ (i.e. $\left\lVert\cdot\right\rVert_{a}\equiv\left\lVert\cdot\right\rVert^{00}_{a}$). The definition in Eq. 100 allows us to define the gyrocenter density $\overline{N}_{a}=\left\lVert 1\right\rVert_{a}$ and, more in general, $\overline{N}_{a}^{lk}=\left\lVert 1\right\rVert^{lk}_{a}/\overline{N}_{a}$ in terms of the Hermite-Laguerre projector operator. Analogously, the gyrocenter perpendicular and parallel pressures can be defined as $\overline{P}_{\perp a}=\overline{N}_{a}\overline{T}_{\perp a}=\left\lVert\overline{\mu}B\right\rVert_{a}$ and $\overline{P}_{\parallel a}=\overline{N}_{a}\overline{T}_{\parallel a}=m_{a}\left\lVert(\overline{v}_{\parallel}-\overline{u}_{\parallel a})^{2}\right\rVert_{a}$ with $\overline{N}_{a}\overline{u}_{\parallel a}=\left\lVert\overline{v}_{\parallel}\right\rVert_{a}$, and the parallel and perpendicular heat fluxes as $\overline{Q}_{\parallel a}=m_{a}\left\lVert(\overline{v}_{\parallel}-\overline{u}_{\parallel a})^{3}\right\rVert_{a}$ and $\overline{Q}_{\perp a}=\left\lVert\overline{\mu}B(\overline{v}_{\parallel}-\overline{u}_{\parallel a})\right\rVert_{a}$, respectively, yielding the following lowest order expansion coefficients, $\overline{N}_{a}^{01}=\overline{N}_{a}^{10}=\overline{N}_{a}^{20}=0$, $\overline{N}_{a}^{00}=1$, and, finally,

[TABLE]

In the present model, contrary to previous gyro-moment hierarchies (see, e.g., Beer & Hammett, 1996; Snyder & Hammett, 2001; Madsen, 2013) that approximate $B_{\parallel}^{*}$ by neglecting the $O(\epsilon)$ terms, we retain the velocity-dependence of $B_{\parallel}^{*}/m_{a}$ [see Eq. 59] in the gyrokinetic Boltzmann equation Eq. 86. For this purpose, we introduce the star Hermite-Laguerre projector operator of order $(l,k)$, defined as

[TABLE]

To simplify the notation also in this case, the $(l,k)=(0,0)$ star Hermite-Laguerre projector operator is denoted by $\left\lVert\cdot\right\rVert^{*}_{a}\equiv\left\lVert\cdot\right\rVert^{*00}_{a}$. Similarly to $\overline{N}_{a}^{lk}$, we define the star gyro-moments $\overline{N}_{a}^{*lk}=\left\lVert 1\right\rVert^{*lk}_{a}$ with, in particular, $\overline{N}_{a}^{*}\equiv\left\lVert 1\right\rVert^{*}_{a}$. The star gyro-moments $\overline{N}_{a}^{*lk}$ can be expressed in terms of the gyro-moments $\overline{N}_{a}^{lk}$ using the recursive property of the Hermite polynomials $H_{l+1}(x)=2xH_{l}(x)-2lH_{l-1}(x)$ in Eq. 102, obtaining

[TABLE]

where $\overline{B}_{\parallel a}^{*}=\overline{\bm{B}}_{a}^{*}\cdot\bm{b}$ with

[TABLE]

being the normalized effective magnetic field evaluated at $\overline{v}_{\parallel}=\overline{u}_{\parallel a}$ [see Eq. 59]. By using Eq. 103, we derive

[TABLE]

As a final remark, we note that Hermite polynomials are a well-known velocity-space basis in plasma physics (see, e.g., Grant & Feix, 1967; Dorland & Hammett, 1993; Hammett et al., 1993; Scott, 2010; Zocco & Schekochihin, 2011; Zocco et al., 2015; Schekochihin et al., 2016; Adkins & Schekochihin, 2018; Pezzi et al., 2019). On the other hand, the use of Laguerre polynomials is more recent (see, e.g., Sugama & Nishimura, 2008; Belli & Candy, 2012; Omotani et al., 2015; Zocco et al., 2015; Jorge et al., 2017; Mandell et al., 2018; Jorge et al., 2018, 2019). In particular, it is motivated for the perpendicular velocity dynamics by their relation with Bessel functions [see Section 5.4].

The Hermite-Laguerre decomposition derived in this work presents several advantages over velocity-space grid methods for a numerical resolution of the gyrokinetic Boltzmann equation given in Eq. 90. For example, velocity-spectral representations allow the evaluation of the $\overline{\mu}$- and $\overline{v}_{\parallel}$-derivatives present in the gyrokinetic Boltzmann equation given in Eq. 90 with no truncation error, and ensure a spectral convergence property such that the error decreases exponentially with the number of gyro-moments retained in the case of distribution functions that are not discontinuous, while velocity-space grid methods follow a power law convergence rate (Boyd, 2001). A possible equivalence between a velocity-space grid method with resolution $(\Delta v_{\parallel},\Delta v_{\perp})$ and a Hermite-Laguerre decomposition with order $(L,K)$ can be obtained by estimating the smallest velocity-space structures $\delta v_{\parallel}$ and $\delta v_{\perp}$ that both methods need to resolve. By balancing the collisional diffusion term, $C\sim\nu/\delta v_{\parallel}^{2}$, and the time variation, $\partial_{t}\sim\omega$, one finds $\Delta v_{\parallel}\lesssim\delta v_{\parallel}\sim\sqrt{\nu/\omega}v_{th\parallel}$, while the Hermite order $L$, at which the gyro-moments become negligible because of collisional effects, can be estimated by following the procedure described in Zocco & Schekochihin (2011); Loureiro et al. (2016); Jorge et al. (2018). From nonlinear perpendicular phase-mixing (Schekochihin et al., 2008), uncorrelated $\delta v_{\perp}$ structures of the order of $v_{th\perp}/k_{\perp}\rho$ can develop. Thus, one has $\Delta v_{\perp}\lesssim\delta v_{\perp}\sim v_{th\perp}/k_{\perp}\rho$, while a rough estimate for the Laguerre order $K$ is $K\gtrsim(k_{\perp}\rho)^{2}$ since $k_{\perp}\rho\simeq\sqrt{2K}$ is the position of the maximum of the kernel function $\mathcal{K}_{K}$ [see Section 5.4].

5.2 Gyro-Moment Hierarchy Equation

The Hermite-Laguerre decomposition of $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ in Eq. 91 is the key step to provide an efficient technique to approach the solution of the five-dimensional (including time) gyrokinetic Boltzmann equation Eq. 90. In fact, the evolution of $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ can be obtained by solving an infinite set of coupled three-dimensional (and time dependent) equations for the gyro-moment $\overline{N}_{a}^{lk}(\overline{\bm{R}},t)$ and the accuracy of the solution, i.e. the degree of fidelity, is directly related to the number of retained gyro-moments. The gyro-moment equation hierarchy is obtained by multiplying the gyrokinetic Boltzmann equation in Eq. 90 by the Hermite-Laguerre basis element $H_{a}^{lk}$ and performing the integral over the velocity space. Before proceeding, we highlight the velocity dependence of the gyrokinetic equations of motion in Eqs. 77 and 78 in terms of $\overline{s}_{\parallel a}$ and $\overline{s}_{\perp a}^{2}$. Thus, we write

[TABLE]

with the gyrocenter drifts

[TABLE]

and with the parallel forces

[TABLE]

We note that the lowest order convective fluid derivative $d_{0t}=\partial_{t}+\bm{U}_{0a}\cdot\overline{\nabla}$, in Eqs. 108 and 109, is associated with the gyrocenter fluid velocity $\bm{U}_{0a}=\bm{u}_{E}+\overline{u}_{\parallel a}\bm{b}$ and that $\Omega_{a}^{*}=q_{a}B_{\parallel}^{*}/m_{a}$ with $B_{\parallel}^{*}$ defined in Eq. 59.

Performing the integral over the velocity space of Eq. 90, we obtain the gyro-moment equation hierarchy of particle species $a$, describing the spatial and temporal evolution of the gyro-moments $\overline{N}_{a}^{*lk}$,

[TABLE]

with $C_{a}^{lk}$ the Hermite-Laguerre projection of the collision operator (see Section 6), and $\mathcal{F}_{a}^{lk}$ the fluid operator given by

[TABLE]

having introduced the gyro-moment convective fluid derivative

[TABLE]

The fluid operator in Section 5.2, $\mathcal{F}_{a}^{lk}$, defines the evolution of the fluid quantities $\overline{N}_{a}$, $\overline{u}_{\parallel a}$, $\overline{P}_{\perp a}$ and $\overline{P}_{\parallel a}$, similarly as for the drift-kinetic model in Jorge et al. (2017). The Hermite-Laguerre projections of the electromagnetic gyrokinetic equations of motion, appearing in Eq. 110, are

[TABLE]

where the drifts and forces have equivalent definitions to the ones in Eqs. 108 and 109, having replaced all $\Omega_{a}^{*}$ with $\Omega_{a}$. The gyro-moment expansion of $\mathcal{A^{*}}_{a}$, $F_{\omega a}^{*}$ and $F_{Ba}^{*}$ in Eq. 114 are given, respectively, by

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

having used Eq. 104. In Eqs. 113 and 114, the parallel and perpendicular phase-mixing operators, $\mathcal{V}^{lk}_{pj}$, $\mathcal{V}_{apj}^{lk}$, $\mathcal{V}^{2lk}_{pj}$ and $\mathcal{M}^{lk}_{pj}$, are defined, respectively, by

[TABLE]

and $\mathcal{V}_{apj}^{2lk}=\sum_{r,s}\mathcal{V}_{ars}^{lk}\mathcal{V}_{apj}^{rs}$. These phase-mixing operators are derived from the recursive properties of the Hermite and Laguerre polynomials, that are $H_{l+1}(x)=2xH_{l}(x)-2lH_{l-1}(x)$ and $xL_{k}(x)=(2k+1)L_{k}(x)-kL_{k-1}(x)-(k+1)L_{k+1}(x)$. From the definitions in Eqs. 120a and 120d, we see that the phase-mixing operator $\mathcal{V}_{pj}^{lk}$ couples the $(l-1)$ and $(l+1)$ Hermite gyro-moments due to parallel streaming, whereas the $\mathcal{M}_{pj}^{lk}$ couples the $(k-1)$, $k$ and $(k+1)$ Laguerre gyro-moments because of the presence of the curvature and gradient of the magnetic field [see Eq. 108]. We remark that, in previously developed gyrofluid models, closure approximations in $\mathcal{V}_{pj}^{lk}$ are used to model the associated linear response of Landau damping, a technique pioneered by Hammett & Perkins (1990) and Dorland & Hammett (1993). Also, toroidal closures to model perpendicular phase-mixing in $\mathcal{M}_{pj}^{lk}$ are provided for these gyrofluid models (see, e.g., Beer & Hammett, 1996; Snyder & Hammett, 2001; Madsen, 2013). Instead, our model retains the full coupling between gyro-moments in both parallel and perpendicular directions, with the aim of evolving the number of moments necessary to obtain the desired accuracy.

We now turn to the evaluations of the Hermite-Laguerre projections of the FLR induced transport term in Eq. 113, that is

[TABLE]

and of the FLR induced parallel force present in Eq. 78,

[TABLE]

Equations 121 and 5.2 show that the Hermite-Laguerre projections of the FLR induced terms are reduced to the evaluation of $\left\lVert\overline{\nabla}\left<\Psi_{1}\right>_{\overline{\bm{R}}}\right\rVert^{lk}_{a}$. The definition of $\Psi_{1}$ in Section 3.4 and the phase-mixing operators in Eq. 120 lead to

[TABLE]

The first two terms on the right-hand side of Section 5.2 are the $O(\epsilon_{\delta})$ contributions, whereas the following terms represent the $O(\epsilon_{\delta}^{2})$ nonlinear contributions to the gyrocenter dynamics given by $\phi_{1}$ and $A_{\parallel 1}$. We note that the first $O(\epsilon_{\delta})$ term provides the Hermite-Laguerre projection of the nonlinear $\bm{E}\times\bm{B}$ advection, proportional to $\bm{b}\times\overline{\nabla}\left<\phi_{1}\right>_{\overline{\bm{R}}}/B\cdot\overline{\nabla}_{\perp}\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ appearing in Eq. 121. As shown in previous gyrofluids models (see, e.g., Brizard, 1992; Dorland & Hammett, 1993; Beer & Hammett, 1996; Madsen, 2013), the gyroaverage operator present in this term couples the Laguerre gyro-moments, yielding FLR phase-mixing terms (Schekochihin et al., 2008; Tatsuno et al., 2009; Schekochihin et al., 2009).

5.3 Lowest order Moment Equations

The infinite set of coupled fluid equations that the gyro-moment hierarchy leads to is valid for arbitrarily far from equilibrium distribution functions. For practical and numerical application, a finite number of gyro-moments should be retained, still ensuring a proper description of the plasma dynamic and be explicitly evaluated. As an example, we derive here the equations for the six lowest order gyro-moments from the gyro-moment hierarchy equation given in Eq. 110. These equations are obtained by a truncation of the Hermite-Laguerre expansion of $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ in Eq. 91, i.e. retaining only the lowest order Hermite-Laguerre polynomials (i.e. $l\leq 3$ and $k\leq 1$), yielding a closed set of fluid equations. This system states the evolution of the gyrocenter density $\overline{N}_{a}$, of the parallel gyrocenter fluid velocity, $\overline{u}_{\parallel a}$, of the parallel and perpendicular gyrocenter temperature, $\overline{T}_{\parallel a}$ and $\overline{T}_{\perp a}$, and, finally, of the parallel and perpendicular heat fluxes, $\overline{Q}_{\parallel a}$ and $\overline{Q}_{\perp a}$. While the Hermite-Laguerre decomposition is truncated here and, in particular, the gyro-moments $\overline{N}_{a}^{40}$ and $\overline{N}_{a}^{02}$ are neglected, we note that a high-collisional closure to express $\overline{Q}_{\parallel a}$ and $\overline{Q}_{\perp a}$ in terms of $\overline{N}_{a}$, $\overline{u}_{\parallel a}$, $\overline{T}_{\parallel a}$ and $\overline{T}_{\perp a}$ can be developed. This can be done by explicit evaluation of the collisional terms. The procedure can be extended to arbitrary collisionality by including a larger number of gyro-moments (Zocco & Schekochihin, 2011; Loureiro et al., 2016; Jorge et al., 2017). We leave to a future work the development of closures of the gyro-moment hierarchy.

The fluid equation of the gyrocenter density, $\overline{N}_{a}$, is obtained by setting $(l,k)=(0,0)$ in Eq. 110. This yields

[TABLE]

having defined the lowest order convective derivative, $d_{ta}^{*0}=\overline{N}_{a}^{*}\partial_{t}+\bm{u}_{a}^{0}\cdot\overline{\nabla}$, and

[TABLE]

In Section 5.3 we introduce $\tau=\bm{b}\cdot\overline{\nabla}\times\bm{b}$, where $\bm{b}\simeq\hat{\bm{b}}+\overline{\nabla}A_{\parallel 0}\times\hat{\bm{b}}/\hat{\bm{B}}$, as a shorthand notation. Section 5.3 is indeed the $(l,k)=(0,0)$ Hermite-Laguerre projection of $\dot{\overline{\bm{R}}}$ [see Eq. 113]. In Eq. 124, we note the presence of polarization effects due to the time-dependence of $\overline{B}_{\parallel a}^{*}/B$, which can be written, from Eq. 59, at the leading order as

[TABLE]

In addition, a non-vanishing lowest-order collisional Hermite-Laguerre moment $C_{a}^{00}$, typically due to finite FLR effects [see Section 6], may be present, leading to classical diffusion of the gyrocenter density. Finally, we note that Eq. 124 preserves the positivity of $\overline{N}_{a}$.

The parallel momentum equation describing the evolution of $\overline{u}_{\parallel a}$ is derived by setting $(l,k)=(1,0)$ in Eq. 110, i.e.

[TABLE]

with

[TABLE]

The parallel and perpendicular temperature equations for $\overline{T}_{\parallel a}$ and $\overline{T}_{\perp a}$ are obtained by setting $(l,k)=(2,0)$ and $(0,1)$ respectively in Eq. 110. This yields for the parallel temperature $\overline{T}_{\parallel a}$,

[TABLE]

with

[TABLE]

while for the perpendicular temperature, $\overline{T}_{\perp a}$, it yields

[TABLE]

with

[TABLE]

Finally, the parallel heat flux equation for $\overline{Q}_{\parallel a}$ is obtained by setting $(l,k)=(3,0)$ in Eq. 110 yielding

[TABLE]

while the perpendicular heat flux equation for $\overline{Q}_{\perp a}$ follows from $(l,k)=(1,1)$,

[TABLE]

with

[TABLE]

and

[TABLE]

In the above fluid equations, we note the presence of multiple partial time derivatives that emerge from the fluid operator $\mathcal{F}_{a}^{lk}$, given in Section 5.2, and from the gyro-moment convective fluid derivative $d_{at}^{*lk}$ defined in Eq. 112. In a numerical implementation, these partial time derivatives can be handled by introducing generalized fluid quantities that include also field dependencies (see, e.g, Halpern et al., 2016).

This system of Eqs. 124, 127, 5.3, 131, 5.3 and 5.3 constitutes the gyrokinetic electromagnetic extension of the collisional drift-reduced fluid equations obtained in Jorge et al. (2017). Indeed, the present equations contain several new terms, proportional to $\tau$, and terms related to the small scale and amplitude electromagnetic fluctuations through the gradient of the gyroaveraged gyrokinetic potential, $\overline{\nabla}\left<\Psi_{1}\right>$. Also, they retain exact polarization effects associated with $B_{\parallel}^{*}$ [see Eq. 126]. The evaluation of these terms was approximated in previous reduced gyro-moment hierarchy (Madsen, 2013).

5.4 Hermite-Laguerre Representation of Gyroaverage Operator

In this section, we overcome the limitations of previous gyrofluid models (see, e.g., Beer & Hammett, 1996; Snyder & Hammett, 2001; Scott, 2010; Madsen, 2013) deriving closed expressions for the FLR induced transports and forces appearing in Eqs. 121 and 5.2. More precisely, we aim to express the gyroaveraged small-scale electromagnetic fluctuations, evaluated at $\overline{\bm{x}}=\overline{\bm{R}}+\overline{\bm{\rho}}_{a}$ [see Eq. 75], in terms of their values at $\overline{\bm{R}}$, i.e. $\phi_{1}(\overline{\bm{R}})$ and $A_{\parallel 1}(\overline{\bm{R}})$ and of the gyro-moments $\overline{N}_{a}^{lk}$. This requires a suitable analytical form of the gyroaverage operator to perform the Hermite-Laguerre projection, while retaining the linear and nonlinear couplings accurately at all orders in $\epsilon_{\perp}$ and translating the spatial dependence of the fields from $\overline{\bm{\bm{x}}}$ to $\overline{\bm{R}}$. In the following, we first evaluate the Hermite-Laguerre representation of the transport term associated with the small-scale electrostatic potential $\left<\phi_{1}\right>_{\overline{\bm{R}}}$, i.e. $\left\lVert\overline{\nabla}\left<\phi_{1}\right>_{\overline{\bm{R}}}\right\rVert^{pj}_{a}$. We then generalize the approach to the nonlinear terms, such as $\overline{\nabla}\partial_{\overline{\mu}}\left<A_{\parallel 1}^{2}\right>_{\overline{\bm{R}}}$, $\overline{\nabla}\partial_{\overline{\mu}}\left<\phi_{1}^{2}\right>_{\overline{\bm{R}}}$, $\overline{\nabla}\partial_{\overline{\mu}}\left<\phi_{1}A_{\parallel 1}\right>_{\overline{\bm{R}}}$, and $\overline{\nabla}\partial_{\overline{\mu}}(\left<\phi_{1}\right>_{\overline{\bm{R}}}\left<A_{\parallel 1}\right>_{\overline{\bm{R}}})$ appearing in Section 5.2, and report the Hermite-Laguerre projection of the last term, proportional to $\overline{\nabla}\left<(\bm{b}\times\overline{\nabla}\overline{\widetilde{\Phi_{1}}})\cdot\overline{\nabla}\widetilde{\Phi_{1}}\right>_{\overline{\bm{R}}}$ in Section 5.2, in Appendix C.

In order to obtain $\left\lVert\left<\phi_{1}\right>_{\overline{\bm{R}}}\right\rVert^{pj}_{a}$, we first gyroaverage the potential $\phi_{1}$ in Fourier space. This is preferred to the Taylor expansion of $\phi_{1}(\overline{\bm{x}})$ around the gyrocenter position $\overline{\bm{R}}$ since a large number of expansion terms would be needed for $\epsilon_{\perp}\sim 1$. By introducing the Fourier decomposition $\phi_{1}(\overline{\bm{x}})=\sum_{\bm{k}}\phi_{1}(\bm{k})e^{i\bm{k}\cdot\overline{\bm{x}}}$ and using the fact that $\overline{\bm{x}}=\overline{\bm{R}}+\overline{\bm{\rho}}_{a}$ [see Eq. 75], one obtains

[TABLE]

where we have oriented $(\bm{e}_{1},\bm{e}_{2},\bm{b})$ such that $\bm{k}=k_{\perp}\bm{e}_{1}+k_{\parallel}\bm{b}$, and the Jacobi-Anger expansion (Abramowitz, 1974) is used,

[TABLE]

with $b=k_{\perp}v_{\perp}/\Omega_{a}$, and $J_{m}$ the $m$th-order Bessel function of the first kind, yielding $\left<e^{i\bm{k}\cdot\overline{\bm{\rho}}_{a}}\right>_{\overline{\bm{R}}}=J_{0}(b)$. We remark that, from Eq. 137 and the fact $J_{0}(b)=\sum_{l}(-1)^{l}b^{2l}/(2^{2l}l!^{2})$, perpendicular FLR couplings between gyro-moments arise with the gyroaverage operator.

To perform the Hermite-Laguerre projection of Eq. 137 and express $\left<\phi_{1}\right>_{\overline{\bm{R}}}$ in terms of gyro-moments $N_{a}^{lk}$, we expand $J_{m}$ in terms of generalized Laguerre polynomials $L_{n}^{m+1/2}$ defined by (Gradshteyn & Ryzhik, 2007)

[TABLE]

where

[TABLE]

that is

[TABLE]

with $\sigma_{0}=1$ and $\sigma_{m}=\text{sgn}(m)^{m}$ for $m\neq 0$. and where we introduced the thermal velocity based parameter $b_{a}=b/\overline{s}_{\perp a}=k_{\perp}/\Omega_{a}\sqrt{2\overline{T}_{\perp a}/m_{a}}$, and the $n$th-order kernel functions $\mathcal{K}_{n}$,

[TABLE]

From Eq. 142, one can observe that the kernel functions satisfy the recursive formula $b_{a}^{2}\mathcal{K}_{n}(b_{a})=4(n+1)\mathcal{K}_{n+1}(b_{a})$ and the normalization relation $\sum_{n\geq 0}\mathcal{K}_{n}(b_{a})=1$. In addition, the asymptotic behaviour of the kernel function is $\mathcal{K}_{n}(b_{a})\sim b_{a}^{2n}$ for $b_{a}\ll 1$, while $\mathcal{K}_{n}\sim 1/n!$ at large $n$ for all values of $b_{a}$. Also, we notice that, for ions, $b_{i}\sim\epsilon_{\perp}$ for ions, while $b_{e}\sim\sqrt{m_{e}/m_{i}}\epsilon_{\perp}$ for electrons.

Using Eq. 141 to expand the Bessel function $J_{0}$, the gyroaveraged electrostatic potential $\left<\phi_{1}\right>_{\overline{\bm{R}}}$ can be expressed as a Laguerre series,

[TABLE]

which decouples the spatial and velocity dependencies. In fact, the velocity dependence that involves the perpendicular FLR coupling between gyro-moments is embedded into the Laguerre polynomials $L_{n}$, whereas the spatial gyroaveraging is handled by the kernel $\mathcal{K}_{n}$. As a consequence of Eq. 143, the accuracy of FLR effects in our description is directly related to the Laguerre resolution (Zocco et al., 2015; Mandell et al., 2018).

We remark that a spatial (instead of Fourier) representation of gyroaveraged quantities, e.g., $\left<\phi_{1}\right>$, can be obtained by performing the transformation $\bm{k}\to\overline{\nabla}/i$ such that $b_{a}^{2j}\to(-1)^{j}2^{j}(\overline{T}_{\perp a}/m_{a}\Omega_{a}^{2})^{j}\overline{\nabla}_{\perp}^{2j}$ in Eq. 142 . Thus, at the lowest order, we can write

[TABLE]

In Section 5.4, $\overline{T}_{\perp a}$ and $\Omega_{a}$ are spatially dependent, but are held constant when gyroaveraging $\phi_{1}$.

We now apply the gradient operator $\overline{\nabla}$ to Eq. 143 and the Hermite-Laguerre projector. Since $b_{a}$ and $\overline{s}_{\perp a}$ are spatially varying as they depend on $B$ and $\overline{T}_{\perp a}$, we use the property of the Laguerre polynomials, $xL_{n}^{\prime}(x)=n[L_{n}(x)-L_{n-1}(x)]$, and the fact that $2\mathcal{K}_{n}^{\prime}(x)=x[\mathcal{K}_{n-1}(x)-\mathcal{K}_{n}(x)]$ (here, the prime denotes the derivative with respect to the argument). In order to handle the products of Laguerre polynomials arising from the Hermite-Laguerre projector, such as, e.g., $H_{a}^{lk}L_{n}$ [see Eq. 92], we use the following relationship (Gillis & Weiss, 1960),

[TABLE]

The expansion coefficients $\alpha^{kn}_{s}$ are determined by the Laguerre polynomial orthogonality relation [see Eq. 96],

[TABLE]

and their closed formula is given by

[TABLE]

where the summation is performed over all possible values of $m$ such that the factorials are positive. From Eq. 146, we note that the coefficients $\alpha_{s}^{kn}$ are symmetric in all pairs of indices and, in particular, that $\alpha^{kn}_{0}=\delta_{k}^{n}$ [see Eq. 96]. We then derive

[TABLE]

where the FLR gradient operator $\bm{\mathcal{D}}^{lkj}_{an}$ is introduced, that is

[TABLE]

with $\eta=B/\overline{T}_{\perp a}$, $\iota=\sqrt{\overline{T}_{\perp a}}/B$, and

[TABLE]

In absence of temperature and magnetic gradients, the gyroaverage and $\overline{\nabla}$ operators commute, i.e. $\overline{\nabla}\left<\phi_{1}\right>_{\overline{\bm{R}}}=\left<\overline{\nabla}\phi_{1}\right>_{\overline{\bm{R}}}$. In this case, $\bm{\mathcal{D}}^{lkj}_{an}(b_{a},\bm{k})=\mathcal{D}_{an}^{lkj}(b_{a})i\bm{k}$. We note that Eq. 150 corresponds to the expression obtained within a similar Laguerre treatment of the gyroaverage operator by Mandell et al. (2018). In addition, in the case of a Maxwellian equilibrium, such that $\overline{N}_{a}^{lk}=0$ for $l,k>0$ [see Eq. 91], we notice that Eq. 150 reduces to $\mathcal{D}_{a0}^{000}=\mathcal{K}_{0}(b_{a})=e^{-(b_{a}/2)^{2}}\equiv\int d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}BJ_{0}(b)F_{aM}/m_{a}$ (Dorland & Hammett, 1993).The second term of the right-hand side of Section 5.2 can be evaluated similarly. When projected, this shear-Alfvén term has the following expression,

[TABLE]

which is equivalent to Eq. 148.

We now evaluate the Hermite-Laguerre projections of the second order terms in Section 5.2, such as $\overline{\nabla}\partial_{\overline{\mu}}\left<\phi_{1}A_{\parallel 1}\right>_{\overline{\bm{R}}}$, $\overline{\nabla}\left<A_{\parallel 1}^{2}\right>_{\overline{\bm{R}}}$, and $\overline{\nabla}\partial_{\overline{\mu}}(\left<A_{\parallel 1}\right>_{\overline{\bm{R}}}\left<\phi_{1}\right>_{\overline{\bm{R}}})$. We consider, as an example, the term $\overline{\nabla}\partial_{\overline{\mu}}\left<\phi_{1}A_{\parallel 1}\right>_{\overline{\bm{R}}}$. As a first step, we write the product of $\phi_{1}$ and $A_{\parallel 1}$ in Fourier space,

[TABLE]

where $\bm{K}=\bm{k}+\bm{k}^{\prime}$ with $\bm{k}^{\prime}=k^{\prime}_{\parallel}\bm{b}+k_{\perp}^{\prime}(\cos\alpha\bm{e}_{1}+\sin\alpha\bm{e}_{2})$ (the angle $\alpha$ being the phase shift between $\bm{k}$ and $\bm{k}^{\prime}$ in the plane perpendicular to $\bm{b}$, such that $i\bm{k}^{\prime}\cdot\overline{\bm{\rho}}_{a}=ik_{\perp}^{\prime}\overline{v}_{\perp}\cos(\theta-\alpha)/\Omega_{a}$), and $J_{n}$ and $J_{m}$ are the Bessel functions introduced via the Jacobi-Anger identity in Eq. 138, with arguments $b=k_{\perp}v_{\perp}/\Omega_{a}$ and $b^{\prime}=k_{\perp}^{\prime}v_{\perp}/\Omega_{a}$, respectively. Equation 152 is in a suitable form to perform the gyroaverage operator. Indeed, using the fact that $\int_{0}^{2\pi}d\overline{\theta}e^{i(n+m)\overline{\theta}}=2\pi\delta_{-n}^{m}$ (Gradshteyn & Ryzhik, 2007), one obtains

[TABLE]

We now express the Bessel functions $J_{n}(b)$ and $J_{n}(b^{\prime})$ in terms of generalized Laguerre polynomials according to Eq. 141,

[TABLE]

The product between two generalized Laguerre polynomials weighted by $\overline{s}_{\perp a}^{2|n|}$, i.e., $L_{r}^{|n|}(\overline{s}_{\perp a}^{2})L_{s}^{|n|}(\overline{s}_{\perp a}^{2})\overline{s}_{\perp a}^{2|n|}$ appearing in Section 5.4, can be written in terms of single Laguerre polynomials using the two following identities,

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

Using consecutively Eq. 155 and Eq. 156, we derive

[TABLE]

where $K_{rsef}^{n}=(-1)^{n}e^{in\alpha}d_{re}^{|n|}\overline{d}_{sef}^{|n|}r!s!/[(|n|+r)!(|n|+s)!]$. By taking into account the spatial variation of the velocity variable $\overline{s}_{\perp a}^{2}$ and using the fact that $\partial/\partial\overline{\mu}=\eta\partial/\partial\overline{s}_{\perp a}^{2}$ (with $\eta=B/\overline{T}_{\perp a}$), one can write

[TABLE]

We focus on the first term in Eq. 160. Taking the $\overline{s}_{\perp a}^{2}$-derivative of Eq. 153 with $L_{n}^{\prime}(x)=(-1)L^{1}_{n-1}(x)$ for $n>0$ and the property $L_{m}^{1}=\sum_{l=0}^{m}L_{l}$ (Gradshteyn & Ryzhik, 2007), one derives

[TABLE]

Section 5.4 has the proper form to apply the Hermite-Laguerre projector in Eq. 100. Indeed, applying the gradient operator $\overline{\nabla}$, we obtain

[TABLE]

where the FLR gradient operator $\bm{\overline{\mathcal{D}}}_{anrs}^{lkg}$ is introduced, and defined by

[TABLE]

where

[TABLE]

We now consider the second term in Eq. 160. Using the fact that $\partial/\partial\overline{\mu}=\eta\partial/\partial\overline{s}_{\perp a}^{2}$ with Section 5.4 yields

[TABLE]

In conclusion, we obtain, from Sections 5.4 and 5.4, the Hermite-Laguerre representation of the second order nonlinear electromagnetic term $\overline{\nabla}\partial_{\overline{\mu}}\left<\phi_{1}A_{\parallel_{1}}\right>_{\overline{\bm{R}}}$,

[TABLE]

A similar procedure can be used to derive

[TABLE]

We finally consider the terms containing the product of gyroaveraged fluctuating fields, such as $\overline{\nabla}\partial_{\overline{\mu}}(\left<\phi_{1}\right>_{\overline{\bm{R}}}\left<A_{\parallel 1}\right>_{\overline{\bm{R}}})$ appearing in Section 5.2. The product $\left<\phi_{1}\right>_{\overline{\bm{R}}}\left<A_{\parallel 1}\right>_{\overline{\bm{R}}}$ can be written in Fourier space and expanded in Laguerre polynomials as follows

[TABLE]

having made use of Eq. 145. Following the procedure used to evaluate $\overline{\nabla}\partial_{\overline{\mu}}\left<\phi_{1}A_{\parallel 1}\right>_{\overline{\bm{R}}}$, we obtain the Hermite-Laguerre representation of $\overline{\nabla}\partial_{\overline{\mu}}\left(\left<\phi_{1}\right>_{\overline{\bm{R}}}\left<A_{\parallel 1}\right>_{\overline{\bm{R}}}\right)$, i.e.

[TABLE]

where the FLR gradient operator is defined by

[TABLE]

with

[TABLE]

The remaining $O(\epsilon_{\delta}^{2})$ Hermite-Laguerre projections in Section 5.2 are then

[TABLE]

Equations 148 and 151 for the first order and Sections 5.4, 167, 169 and 169 for the second order terms provide the complete Hermite-Laguerre projection of the gradient of the gyroaveraged gyrokinetic potential $\Psi_{1}$ in Section 5.2. By using these results, $\overline{\nabla}\left<\Psi_{1}\right>_{\overline{\bm{R}}}$ can be fully expressed in terms of the gyro-moments $\overline{N}_{a}^{lk}(\overline{\bm{R}})$ and of the Fourier components of the fluctuating fields, $\phi_{1}(\bm{k})$ and $A_{\parallel 1}(\bm{k}^{\prime})$.

We remark that a spatial (instead of Fourier) representation of the different FLR terms can be adopted by performing the inverse Fourier transform in the above expressions. The FLR closure can be obtained by using a Padé-based approximation closure at arbitrary wavelength for a numerical implementation, avoiding cumbersome convolutions in Fourier space (Held et al., 2019). Indeed, a Padé-based approximation can be applied be noticing that the kernel function $\mathcal{K}_{n}$, defined in Eq. 142, can be approximated by

[TABLE]

where $\sqrt{\Gamma_{0}}$ is the FLR operator introduced by Dorland & Hammett (1993), and $\Gamma_{n}(x)=I_{n}(x)e^{-x}$ (with $I_{n}(x)=i^{-n}J_{n}(ix)$ the modified Bessel function). This FLR operator stems from the fact that $\int d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}BJ_{0}(b)\overline{F}_{aM}/m\simeq\sqrt{\Gamma_{0}}$. The approximated higher order kernels follow by using exact recursive properties.

With the Hermite-Laguerre representation of the FLR terms and the phase-mixing operators in Eq. 120, the full perpendicular and parallel FLR couplings between gyro-moments are exactly represented at arbitrary order in $\epsilon_{\perp}$ in the gyro-moment hierarchy, Eq. 110. Our formulation provides an exact benchmark for FLR closures considered by previous models (Dorland & Hammett, 1993; Beer & Hammett, 1996; Snyder & Hammett, 2001; Madsen, 2013), and allows for a possible dynamical refinement of the plasma description based on the Hermite-Laguerre resolution.

6 Gyrokinetic Collision Operator

The development of a proper gyrokinetic collision operator has been subject of large analytical (Catto & Tsang, 1977; Sugama et al., 2009; Li & Ernst, 2011; Madsen, 2013; Hirvijoki et al., 2017; Pezzi et al., 2019; Pan & Ernst, 2019) and numerical (Abel et al., 2008; Barnes et al., 2009; Estève et al., 2015) efforts, since collisions provide a transport mechanism and influence turbulence and its associated transport. We provide herein the gyro-moment expansion of a relatively simple nonlinear inter-species collision operator, the Dougherty collision operator (Dougherty, 1964). The expansion is valid at arbitrary $\epsilon_{\perp}$ values. Despite its functional simplicity, the Dougherty collision operator has the necessary field-particle terms ensuring conservation of particles, momentum, and energy. It contains pitch-angle scattering, and satisfies the H-theorem (see Dougherty, 1964).

Starting from the Fokker-Planck form of the collision operator,

[TABLE]

where $\nu_{a}=4\pi n_{a}q_{a}^{4}\ln\Lambda/(m_{a}^{2}v_{tha}^{3})$ is the velocity-independent collision frequency, $\ln\Lambda$ is the Coulomb logarithm, $\bm{D}(\bm{z})$ is the velocity diffusion tensor, and $\bm{P}(\bm{z})$ is the friction drag from the background particles.The Dougherty collision operator can be expressed in a Fokker-Planck form, with the velocity diffusion tensor taken to be isotropic and velocity-independent, i.e. $\bm{D}(\bm{z})\equiv T_{a}(\bm{x})\bm{1}$, and the friction drag to be proportional to the velocity difference with the mean particle fluid velocity $\bm{u}_{a}$, i.e. $\bm{P}(\bm{z})\equiv\bm{u}_{a}(\bm{x})-\bm{v}$. More precisely, the Dougherty collision operator expressed in the $\bm{z}=(\bm{x},\bm{v})$ coordinates is defined by

[TABLE]

We remark that the temperature $T_{a}$ and the fluid velocity $\bm{u}_{a}$ appearing in Eq. 178 are moments of the particle distribution function $f_{a}$, and are spatially dependent, i.e. $T_{a}=T_{a}(\bm{x})$ and $\bm{u}_{a}=\bm{u}_{a}(\bm{x})$. They are defined in terms of the particle distribution function $f_{a}=f_{a}(\bm{z})$ as $T_{a}=\int d\bm{v}f_{a}m_{a}(\bm{v}-\bm{u}_{a})^{2}/3$ and as $n_{a}\bm{u}_{a}=\int d\bm{v}f_{a}\bm{v}$ with the particle density $n_{a}=\int d\bm{v}f_{a}$. We note that the collision operators expressed in a Fokker-Planck form - and therefore the Dougherty collision operator - conserve the positivity of the distribution function $f_{a}$, and of $\overline{F_{a}}$ from the scalar invariance $f_{a}(\bm{z})=\overline{F_{a}}(\overline{\bm{Z}}(\bm{z}))$. In fact, if $f_{a}\geq 0$ is twice differentiable, and if it exists $\bm{v}^{*}$ such that $f_{a}(\bm{x},\bm{v}^{*})=0$ and $f_{a}(\bm{x},\bm{v})>0$ for all $\bm{v}\neq\bm{v}^{*}$ with $\partial_{\bm{v}}f_{a}(\bm{x},\bm{v}^{*})=0$, then $\partial_{t}f_{a}(\bm{x},\bm{v}^{*})=\nu_{a}\bm{D}:{\partial^{2}_{\bm{v}\bm{v}}}f_{a}(\bm{x},\bm{v}^{*})\geq 0$ provided that $\bm{D}:\partial^{2}_{\bm{v}\bm{v}}f_{a}(\bm{x},\bm{v}^{*})$ is semi positive-definite, such as in the Dougherty collision operator.

The gyrokinetic Dougherty collision operator is obtained as an expansion in the small parameter $\epsilon_{\delta}$ up to $O(\epsilon_{\nu}\epsilon_{\delta})$ (the positivity of the distribution function $\overline{F_{a}}$ is therefore ensured up to $O(\epsilon_{\nu}\epsilon_{\delta})$). We first express the velocity derivatives of $f_{a}(\bm{z})$ appearing in Eq. 178 in terms of the gyrocenter coordinates $\overline{\bm{Z}}$ derivatives by using the chain rule,

[TABLE]

where we use the fact that $f_{a}(\bm{z})=\overline{F_{a}}(\overline{\bm{Z}}(\bm{z}))$. Consistently with the accuracy requirement for the collision operator (see Section 4), we approximate $C_{a}(\overline{F_{a}})\simeq C_{a0}(\overline{F_{a}})$ with $C_{a0}(\overline{F_{a}})$ the lowest order collision operator in $\epsilon_{\delta}$ and, therefore, we neglect the $O(\epsilon_{\delta})$ terms in the derivatives in Eq. 179 [see Eq. 73]. At the lowest-order in $\epsilon_{\delta}$, we have

[TABLE]

Additionally, for the same reason, we approximate $\overline{F_{a}}\simeq\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$ for both electrons and ions [see Eqs. 88 and 89] and, in particular, $\partial_{\overline{\theta}}\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}=0$ in Eq. 179. In order to express the collision operator in Eq. 178 as a function of $\overline{\bm{Z}}$, we first write the fluid quantities $T_{a}$ and $\bm{u}_{a}$ in Fourier harmonics, i.e. $T_{a}(\bm{x})=\sum_{\bm{k}^{\prime}}T_{a}(\bm{k}^{\prime})e^{i\bm{k}^{\prime}\cdot\bm{x}}$, $\bm{u}_{a}(\bm{x})=\sum_{\bm{k}^{\prime}}\bm{u}_{a}(\bm{k}^{\prime})e^{i\bm{k}^{\prime}\cdot\bm{x}}$. We also note that $\bm{x}=\overline{\bm{R}}+\overline{\bm{\rho}}_{a}$. Then, introducing the Fourier expansion $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}=\sum_{\bm{k}}\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\bm{k})e^{i\bm{k}\cdot\overline{\bm{R}}}$ (we use the shorthand notation $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\bm{k})\equiv\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\bm{k},\overline{\mu},\overline{v}_{\parallel})$) and noticing that $-\bm{a}/\Omega_{a}\cdot\overline{\nabla}=-ik_{\perp}\cos\overline{\theta}/\Omega_{a}$ in Eq. 180, the Dougherty collision operator, expressed as a function of the gyrocenter coordinates $\overline{\bm{Z}}$, is

[TABLE]

where $\bm{K}=\bm{k}+\bm{k}^{\prime}$ and $u_{\parallel a}=\bm{b}\cdot\bm{u}_{a}$. The collision operator in Section 6 is in a suitable form to be gyroaveraged. Applying the gyroaverage operator to Section 6 by noticing that the gyrophase dependence is present in $e^{i\bm{k}^{\prime}\cdot\overline{\bm{\rho}}_{a}}$ and in the terms proportional to $\cos\overline{\theta}$ and $\cos^{2}\overline{\theta}$, and by using the Jacobi-Anger expansion, Eq. 138, to evaluate

[TABLE]

where $\Delta_{n}=e^{-in\alpha}+e^{in\alpha}$, we obtain

[TABLE]

To derive the Hermite-Laguerre representation of the gyroaveraged Dougherty collision operator, i.e.

[TABLE]

appearing on the right-hand side of Eq. 90, we relate the Fourier harmonics of particle fluid quantities, $T_{a}(\bm{k}^{\prime})$ and $\bm{u}_{a}(\bm{k}^{\prime})$ in Section 6 to gyrocenter fluid quantities, i.e. to the moments of $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\bm{k}^{\prime})$. First, we note that the velocity integrals can be expressed as phase-space integrals. In fact, by approximating $\overline{\bm{Z}}\simeq\bm{Z}$, we have $f_{a}(\bm{z})=\int d\overline{\bm{R}}\>\overline{F_{a}}(\overline{\bm{Z}})\delta(\overline{\bm{R}}+\overline{\bm{\rho}}_{a}-\bm{x})$ with $\overline{F_{a}}\simeq\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}$. Second, expanding in Fourier harmonics $f_{a}(\bm{z})=\sum_{\bm{k}}f_{a}(\bm{k})e^{i\bm{k}\cdot\bm{x}}$ with $f_{a}(\bm{k})\equiv f_{a}(\bm{k},\bm{v})$, and $\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\overline{\bm{Z}})=\sum_{\bm{k}}\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\bm{k})e^{i\bm{k}\cdot\overline{\bm{R}}}$, we derive that $f_{a}(\bm{k})=\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}(\bm{k})e^{i\bm{k}\cdot\overline{\bm{\rho}}_{a}}$. Then, the particle density $n_{a}(\bm{k}^{\prime})=\int d\bm{v}f_{a}(\bm{k}^{\prime})$ can be written as

[TABLE]

where the Jacobi-Anger expansion in Eq. 138, the Laguerre expansion of Bessel function in Eq. 141, and the definition of $N_{a}^{*lk}$ [see Eq. 102] are used. Analogously, the particle temperature $T_{a}(\bm{k}^{\prime})=[T_{\parallel a}(\bm{k}^{\prime})+2T_{\perp a}(\bm{k}^{\prime})]/3$, with $T_{\parallel a}(\bm{k}^{\prime})=\int d\bm{v}m_{a}f_{a}(\bm{k}^{\prime})(v_{\parallel}-u_{\parallel a})^{2}/n_{a}$ and $T_{\perp a}(\bm{k}^{\prime})=\int d\bm{v}f_{a}(\bm{k}^{\prime})\mu B/n_{a}$, is expressed using

[TABLE]

and

[TABLE]

The parallel fluid velocity $u_{\parallel a}(\bm{k}^{\prime})$ is given by

[TABLE]

Finally, the perpendicular fluid velocity $\bm{u}_{\perp a}(\bm{k}^{\prime})=\bm{u}_{a}(\bm{k}^{\prime})-\bm{b}\cdot\bm{u}_{a}(\bm{k}^{\prime})$, expressed by

[TABLE]

can be written in terms of gyro-moments by projecting the first order Bessel function $J_{1}$ into Laguerre polynomials using Eq. 141. Thus, we derive from Eq. 189 that

[TABLE]

We remark that, even if the direction of $\bm{u}_{\perp a}$ depends on the basis vector $\bm{e}_{2}$ in Eq. 190, the collisional friction term in Section 6 does not. Indeed, the perpendicular component of the collisional friction term appearing in Section 6 is proportional to $(\bm{v}_{\perp}-\bm{u}_{\perp a})\cdot\bm{c}=v_{\perp}-\cos\theta\bm{u}_{\perp a}\cdot\bm{e}_{2}$. Finally, we multiply Section 6 by the Hermite-Laguerre basis element, $H_{a}^{lk}$, and perform the velocity integral by expanding the first and second-order Bessel function, $J_{1}$ and $J_{2}$, in Laguerre polynomials using Eq. 141. Thus, the Hermite-Laguerre representation of the nonlinear gyrokinetic Dougherty collision operator, accurate at arbitrary values of $\epsilon_{\perp}$, is given by

[TABLE]

Here, the test-particle pitch-angle scattering term $C_{apq}^{lk}$ is,

[TABLE]

and the field-particle collisional term $D_{anpq}^{lk}$, associated with velocity diffusion, is

[TABLE]

with

[TABLE]

Finally, in Eq. 192, we introduce the field-particle collision term $P_{anpq}^{lk}$, associated with the fluid friction force, defined by

[TABLE]

where

[TABLE]

respectively. In Eq. 192, the first three terms are associated with pitch-angle and energy scattering due to inter-species collisions, whereas the last term arises from the $\overline{u}_{\parallel a}$ dependence of the velocity coordinate $\overline{s}_{\parallel a}$ [see Eq. 91]. Also, we remark the presence of hyper collisional gyro-diffusion $\sim\nu_{a}b_{a}^{2n}$ in velocity-space in Eqs. 193 and 197. These terms arise from the gyroaverage of the collision operator (Catto & Tsang, 1977; Abel et al., 2008; Li & Ernst, 2011). We notice that this collisional gyro-diffusion yields a classical diffusion in the gyrocenter continuity equation, obtained by setting $(l,k)=(0,0)$ in the gyro-moment hierarchy equation in Eq. 110, associated with FLR effects (Sugama et al., 2015, 2017; Mandell et al., 2018).

Finally, it is instructive to consider the drift-kinetic limit of the gyrokinetic Dougherty collision operator in Eq. 191. We take $b_{a}\sim\epsilon_{\perp}\sim\epsilon$ and neglect the $O(\epsilon_{\perp}^{2})$ terms. Using the fact that $\mathcal{K}_{n}(b_{a})\sim b_{a}^{2n}$ with $b_{a}\ll 1$ for $n>1$, while $\mathcal{K}_{0}(b_{a})\simeq 1-b_{a}^{2}/4$, we derive $n_{a}=\overline{N}_{a}$, $u_{\parallel a}=\overline{u}_{\parallel a}$, $T_{\perp a}=\overline{T}_{\perp a}$, $T_{\parallel a}=\overline{T}_{\parallel a}$ from Eqs. 186, 187, 188 and 190, in agreement with Jorge et al. (2017). Thus, from Eqs. 193 and 197, the field-particle collisional terms reduce to

[TABLE]

and

[TABLE]

As a consequence, the drift-kinetic Doughery collision operator takes the form

[TABLE]

having performed the inverse Fourier transform. We remark that the absence of gyro-diffusion in the gyrocenter density in the drift-kinetic regime, since $C_{a}^{00}=0$ from Eq. 203.

7 Gyrokinetic Maxwell’s Equations

In this section, we derive a set of gyrokinetic Maxwell’s equation that describe the temporal and spatial evolution of $\phi_{0}$, $\phi_{1}$, $A_{\parallel 0}$, $A_{\parallel 1}$. Moreover, we provide an additional field equation to self-consistently obtain the large scale magnetic vector potential $\hat{\bm{A}}$, and therefore $\hat{\bm{B}}$. These equations are derived self-consistently from a variational principal. While previous gyrokinetic field theories (see, e.g., Brizard (2000); Madsen (2010); Tronko et al. (2016)) have been used for collisionless plasmas, Sugama et al. (2015) shows that the conservation laws obtained from collisionless gyrokinetic field theories are still valid in collisional plasmas provided that the collision operator model conserves particle, momentum, and energy. Therefore, similar procedures to the ones applied in, e.g., Sugama (2000); Brizard (2010); Tronko et al. (2016), and based on Noether’s method, can be applied to the present model in order to derive energy, parallel and toroidal angular momentum conservation laws. These invariant properties are preserved by the Hermite-Laguerre projection in the infinite moment limit, and are altered if one considers a truncated gyro-moment hierarchy. However, alternative conservation laws associated with the truncated gyro-moment hierarchy can still be derived within the same formalism, by using the truncated distribution function to evaluate the gyrokinetic functional action (Madsen, 2013). The derivation of conservation laws is outside the scope of the present work.

Contrary to Brizard (2000); Tronko et al. (2016), in our approach, we do not treat the distribution function as a dynamical field. Instead, we derive the equations for the electromagnetic fields, assuming the single particle dynamics model in Section 3. Our field equations include the effects of strong flows and retain full FLR effects at arbitrary $\epsilon_{\perp}$ values to predict accurately the long and short wavelength components of the fluctuating electromagnetic fields (Qin et al., 1998; Parra & Catto, 2008; Miyato et al., 2013), as they might be important in the description of anomalous transport in the periphery. Within this framework, the polarization and magnetization current densities, associated with the particle and gyrocenter difference yields a classical physical interpretation of the gyrokinetic medium (Qin et al., 2000; Krommes, 2012). We compare our results with previous derivations and give a simple physical interpretation of the obtained equations.

The present section is organized as follows. In Section 7.1, we introduce the general formalism of our variational principle. Then, in Sections 7.2 and 7.3, we derive the gyrokinetic Poisson’s equations and gyrokinetic Ampere’s laws from the least action principle and recover known results. Finally, in Section 7.4, we obtain an Ampere’s law that sets the evolution of $\hat{\bm{A}}$, and show that our gyrokinetic formalism encompasses a leading order pressure balance equation.

7.1 Gyrokinetic Field Theory

Following previous gyrokinetic literature (Sugama, 2000; Brizard, 2000, 2010; Squire et al., 2013; Sugama et al., 2014, 2015; Tronko et al., 2016), we introduce the gyrokinetic functional action $\mathcal{A}$,

[TABLE]

where $\mathcal{A}_{p}$ is the gyrocenter functional action associated with the gyrocenter dynamics, $\mathcal{A}_{f}$ is the field functional action containing the fields contributions, and, finally, $\mathcal{A}_{c}$ is a collisional functional action that we assume independent of the dynamical fields. The field functional action $\mathcal{A}_{f}$ is given by

[TABLE]

where the total electric field is $\bm{E}+\bm{E}_{1}$ with $\bm{E}=-\nabla\phi_{0}-\partial_{t}(\hat{\bm{A}}+\hat{\bm{b}}A_{\parallel 0})$ and $\bm{E}_{1}=-\nabla\phi_{1}-\partial_{t}(A_{\parallel 1}\hat{\bm{b}})$, and the total magnetic field is $\bm{B}+\bm{B}_{1}$ where $\bm{B}=\hat{\bm{B}}+\nabla\times(A_{\parallel 0}\hat{\bm{b}})$ with $\hat{\bm{B}}=\nabla\times\hat{\bm{A}}$ and $\bm{B}_{1}=\nabla\times(A_{\parallel 1}\hat{\bm{b}})$ [see Section 2]. In Eq. 205, the electromagnetic fields are evaluated at the particle position $\bm{x}$. We remark that the contribution from the large scale magnetic vector potential $\hat{\bm{A}}$ is included since it is a dynamical field in our approach. We remark that the quasineutrality can be imposed by neglecting the $|\bm{E}+\bm{E}_{1}|^{2}$ term in $\mathcal{A}_{f}$ as usually applied in gyrokinetic models (see, e.g., Krause et al., 2007; Madsen, 2013; Bottino & Sonnendrucker, 2015; Tronko et al., 2016, 2017b). Moreover, the inductive part of the electric field is retained in Eq. 205 since neglecting it in $\mathcal{A}_{f}$ can lead to spurious terms in the local conservation laws (Correa-Restrepo & Pfirsch, 2005).

The gyrocenter functional action, $\mathcal{A}_{p}$, is defined as

[TABLE]

with $d\overline{\bm{Z}}=d\overline{\bm{R}}d\overline{\mu}d\overline{v}_{\parallel}d\overline{\theta}$, and $\overline{L}_{a}$ the single gyrocenter Lagrangian obtained from the gyrocenter one-form in Eq. 70,

[TABLE]

We remark that in Eq. 207, the fields $\phi_{0}$, $\hat{\bm{A}}$ and $A_{\parallel 0}$ are evaluated at $\overline{\bm{R}}$, while $\phi_{1}$ and $A_{\parallel 1}$ at $\overline{\bm{R}}+\overline{\bm{\rho}}_{a}$ being gyroaveraged at constant $\overline{\bm{R}}$ in $\overline{\mathcal{H}}$. This difference in the spatial argument between $\mathcal{A}_{f}$ and $\mathcal{A}_{p}$ yields polarization and magnetization effects. Since the $O(\epsilon^{2},\epsilon_{\delta}^{2})$ single gyrocenter Lagrangian $\overline{L}_{a}$ is gyrophase independent [see Eq. 70], we can expand $\overline{F_{a}}=\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}+\widetilde{\overline{F_{a}}}$ and perform the $\overline{\theta}$-integral in Eq. 206, showing that the contribution of $\widetilde{\overline{F_{a}}}$ in $\mathcal{A}_{p}$ vanishes.

We now derive a set of gyrokinetic field equations, i.e. the gyrokinetic Poisson’s equations and the Ampere’s laws, from the least action principle

[TABLE]

where the total variation of the gyrokinetic action $\mathcal{A}$, in Eq. 204, is given by

[TABLE]

with $(\check{\phi}_{0},\check{\phi}_{1},\check{\hat{\bm{A}}},\check{A}_{\parallel 0},\check{A}_{\parallel 1})$ arbitrary test functions. The $\circ$ notation in Eq. 209 denotes the functional derivative along an arbitrary test function, such as, e.g.,

[TABLE]

with the spatial argument of the test function $\check{\phi}_{0}$ being the same as $\phi_{0}$. Analogously, the variation with respect to a vector function is defined as the sum of the variation with respect to its component, i.e. $\delta\mathcal{A}/\delta\hat{\bm{A}}\circ\check{\hat{\bm{A}}}\equiv\delta\mathcal{A}/\delta\hat{\bm{A}}_{i}\circ\check{\hat{\bm{A}}}_{i}$.

We remark that, from the total variation in Eq. 209, two gyrokinetic Poisson’s equations and three Ampere’s laws are obtained. These coupled equations are necessary to determine the fields $\phi_{0}$, $\phi_{1}$, $A_{\parallel 0}$, and $A_{\parallel 1}$. We note that, by imposing the quasineutrality condition, a plasma vorticity equation is derived. One additional constraint, which corresponds to the pressure balance obtained from the evolution equation of $\hat{\bm{A}}$, can also be derived, as explained in Section 7.4.

7.2 Gyrokinetic Poisson’s Equations

The variation of the action $\mathcal{A}$ with respect to $\phi_{0}$ yields the first gyrokinetic Poisson’s equation, referred to as GKPI, whereas the variation with respect to $\phi_{1}$ produces a second gyrokinetic Poisson’s equation that we denote by GKPII.

We first compute the variation of the action $\mathcal{A}$ with respect to $\phi_{0}$, and obtain GKPI by imposing that

[TABLE]

The functional derivative of the field functional action, $\mathcal{A}_{f}$ in Eq. 205, is given by

[TABLE]

where the test function $\check{\phi}_{0}$ is evaluated at $\bm{x}$, i.e. $\check{\phi}_{0}=\check{\phi}_{0}(\bm{x})$. The first term in Eq. 212 is identified as the dynamical term, i.e. a term entering in the field equation, whereas the second one is a boundary term that vanishes by carrying out the integration and using a proper choice of $\check{\phi}_{0}$.

We now compute the variation of the gyrocenter functional action $\mathcal{A}_{p}$ in Eq. 206. We remark that the single gyrocenter Lagrangian $\overline{L}_{a}$, in Eq. 207 depends on $\phi_{0}$ and on its gradient $\overline{\nabla}\phi_{0}$, with both quantities evaluated at $\overline{\bm{R}}$, through the electrostatic energy, $q_{a}\phi_{0}$, the symplectic components, $q_{a}\overline{\bm{A}^{*}}$, the $\bm{E}\times\bm{B}$ kinetic energy, proportional to $|\overline{\nabla}\phi_{0}|^{2}$, and, finally, the guiding-center FLR correction term, proportional to $(\mu B/2\Omega_{a})\bm{b}\cdot\overline{\nabla}\times(\bm{b}\times\overline{\nabla}\phi_{0}/B)$. The variation of $\mathcal{A}_{p}$ in Eq. 206 along $\check{\phi}_{0}$ is

[TABLE]

where the spatial argument of the test function $\check{\phi}_{0}$ is translated from $\overline{\bm{R}}$ to the particle position $\bm{x}$ by noticing $\check{\phi}_{0}(\overline{\bm{R}})=\int d\bm{x}\delta(\overline{\bm{R}}-\bm{x})\check{\phi}_{0}(\bm{x})$. This is to be consistent with the fact that the variation of field functional action, $\mathcal{A}_{f}$ in Eq. 212, is evaluated at the particle position $\bm{x}$. To evaluate the terms that provide dynamics of $\phi_{0}$ in Section 7.2, we carry out successive integration by parts on terms where derivatives of $\check{\phi}_{0}$ appear. We obtain GKPI from Eq. 211, Eq. 212, and Section 7.2, using the arbitrariness of the test function $\check{\phi}_{0}$ under the assumption that it vanishes on the boundary. Since GKPI is a functional relation, we evaluate GKPI at $\overline{\bm{R}}$ as this is where the gyro-moments are evaluated [see Eq. 110], and in agreement with previous gyrokinetic models (Dannert & Jenko, 2005; Pan et al., 2018). It yields

[TABLE]

In Eq. 214, the gyrocenter charge density $\overline{\varrho}_{a}^{*}$ is

[TABLE]

whereas the polarization charge density, $-\overline{\nabla}\cdot\bm{\mathcal{P}}_{a}^{*}$, is associated with the polarization

[TABLE]

where we have introduced the Pfirsch-Kaufman (PK) polarization (Pfirsch & Morrison, 1985; Kaufman, 1986),

[TABLE]

and the diamagnetic polarization

[TABLE]

We remark that the polarization $\bm{\mathcal{P}}_{a}^{PK*}$ in Eq. 217 is a known term in guiding-center theories (see, e.g., Kaufman, 1986; Pfirsch & Morrison, 1985; Lee & Kolesnikov, 2009; Brizard & Tronko, 2011). This polarization term, proportional to the difference $\dot{\overline{\bm{R}}}-\bm{u}_{E}$ and a consequence of writing the particle velocity $v_{\perp}$ in the frame moving with the background $\bm{E}\times\bm{B}$, is due to the presence of the polarization drift, $\bm{b}\times d_{t}^{0}\bm{U}_{0}/\Omega_{a}$, the magnetic drifts, and the small-scale drifts, $\bm{B}\times\overline{\nabla}_{\perp}\left<\Psi_{1}\right>_{\overline{\bm{R}}}/B^{2}$. The $O(\epsilon^{2})$ guiding-center FLR correction term in Section 7.2 yields the polarization $\bm{\mathcal{P}}_{a}^{D*}$ in Eq. 218, which is associated with the perpendicular fluid pressure.

We now derive the Hermite-Laguerre representation of GKPI in Eq. 214. We perform the velocity integral in Eq. 215 to find

[TABLE]

A similar procedure can be used to obtain the Hermite-Laguerre projections of the polarization terms, in Eqs. 217 and 218, that are

[TABLE]

where we use Eq. 113 with the gyrocenter drifts defined in Eq. 108, and

[TABLE]

respectively.

It is instructive to consider the drift-kinetic limit of GKPI to verify its consistency with previous results. For this purpose, we keep the polarization densities up to $O(\epsilon^{2})$ in Eq. 216. Therefore, only the diamagnetic polarization, $\bm{\mathcal{P}}_{a}^{D*}$, is retained in $\bm{\mathcal{P}}_{a}^{*}$ [the $\bm{\mathcal{P}}_{a}^{PK*}$ polarization is $O(\epsilon^{2},\epsilon_{B})$, neglecting the $O(\epsilon_{\delta})$ and higher corrections]. Approximating $B_{\parallel}^{*}\simeq B+O(\epsilon^{2})$ in Eq. 218 and using $\bm{b}\times(\overline{\nabla}\times\bm{b})\simeq-\overline{\nabla}_{\perp}B/B$, we derive the drift-kinetic Poisson’s equation from Eq. 214,

[TABLE]

Equation 222 corresponds to the drift-kinetic Poisson’s equation used by Jorge et al. (2017) with $\bm{E}=-\overline{\nabla}\phi_{0}$.

We now aim to derive the second gyrokinetic Poisson’s equation, GKPII, from the variation of $\mathcal{A}$ in Eq. 209 with respect to $\phi_{1}$, i.e.

[TABLE]

We first notice that the variation of the field functional action $\mathcal{A}_{f}$ with respect to $\phi_{1}$ has the same functional form as Eq. 212, i.e.

[TABLE]

with the test function $\check{\phi}_{1}$ evaluated at $\bm{x}$. To evaluate the variation of the gyrocenter functional action $\mathcal{A}_{p}$ in Eq. 206, we notice that the $\phi_{1}$ dependent terms in the gyrocenter Lagrangian $L_{a}$ in Eq. 207 are contained only in the gyrokinetic potential $\Psi_{1}$. Due to the complexity of the $O(\epsilon_{\delta}^{2})$ term in $\overline{\mathcal{H}}$ contained in $\left<\Psi_{1}\right>_{\overline{\bm{R}}}$ [see Section 3.4] and for numerical applications, a gyrokinetic long wavelength limit is usually considered, neglecting the $O(\epsilon_{\perp}^{3})$ FLR corrections (see, e.g., Dubin et al., 1983; Lee, 1983; Hahm, 1988; Xu et al., 2007; Cohen & Xu, 2008; Hahm et al., 2009; Madsen, 2013; Tronko et al., 2016, 2017a, 2017b; Shi et al., 2017). However, it has been argued that higher order terms in $\epsilon_{\perp}$ are needed to correctly study also long wavelength modes and to predict the turbulent and neoclassical transport resulting from the nonlinear interactions between scales (Scott, 2003; Parra & Catto, 2008; Lee & Kolesnikov, 2009; Miyato et al., 2013). We therefore retain the full expression of $\left<\Psi_{1}\right>$, and, as a consistency check, we show that the commonly used expression of the linear polarization is recovered in the long wavelength limit.

The variation of the gyrocenter functional action $\mathcal{A}_{p}$ with respect to $\phi_{1}$ is

[TABLE]

In Section 7.2, using that $\check{\phi}_{1}(\overline{\bm{R}}+\overline{\bm{\rho}}_{a})=\int d\bm{x}\delta(\overline{\bm{R}}+\overline{\bm{\rho}}_{a}-\bm{x})\check{\phi}_{1}(\bm{x})$, we introduce the adjoint gyroaverage operator, $\left<\xi\right>_{\bm{x}}^{\dagger}=\left<\xi\right>_{\bm{x}}^{\dagger}(\bm{x},\overline{\mu},\overline{v}_{\parallel},t)$ acting on a function $\xi=\xi(\overline{\bm{R}},\overline{\mu},\overline{v}_{\parallel},t)$, defined by

[TABLE]

and, from Section 3.4, satisfies the property that

[TABLE]

In Section 7.2, the adjoint gyroaverage operator, $\left<\cdot\right>_{\bm{x}}^{\dagger}$, arises from the fact that the gyrocenter distribution function is gyroaveraged in the field equations giving the averaged density (or current density) at the particle position $\bm{x}$. Using Eq. 143, we introduce the adjoint kernel $\mathcal{K}_{n}(b_{a}^{\dagger})$ associated with the Laguerre expansion of the adjoint gyroaverage operator $\left<\cdot\right>_{\bm{x}}^{\dagger}$, where $b_{a}^{\dagger}$ is such that $(b_{a}^{\dagger})^{2j}\to(-1)^{j}2^{j}\overline{\nabla}_{\perp}^{2}(\overline{T}_{\perp a}/m_{a}\Omega_{a}^{2})^{j}$ (Strintzi & Scott, 2004). With the adjoint kernel and Eq. 143, we obtain the spatial representation of the gyroaverage adjoint operator,

[TABLE]

We remark that, consistently with Eq. 227, the spatial dependence of fluid quantities contained in $b_{a}^{\dagger}$, such as $\overline{N}_{a}$ and $\overline{T}_{\perp a}$, and contained in the argument $\overline{s}_{\perp a}^{2}$ of the Laguerre polynomials must be retained when evaluating the adjoint gyroaverage operator. Finally, we notice that it can be shown that the polarization associated with the last term in Section 7.2 is smaller by at least one order compared to the second term in the same equation (Hahm et al., 2009), and can therefore be neglected.

GKPII is obtained from Section 7.2 by first performing the $\overline{\bm{R}}$-integral, that is contained in the $\left<\cdot\right>_{\bm{x}}^{\dagger}$ operator, by expanding the gyroradius, $\overline{\bm{\rho}}_{a}=\bm{\rho}_{a}(\overline{\bm{R}},\overline{\mu},\overline{\theta})$ defined in Eq. 51 that appears in $\delta(\overline{\bm{R}}+\overline{\bm{\rho}}_{a}-\bm{x})$, such as

[TABLE]

where the second term can be neglected being $O(\epsilon_{B})$ [see Eq. 4]. We remark that, this term may become important in the case where $L_{\phi}\sim L_{P}\sim L_{B}$, i.e. $\epsilon\sim\epsilon_{B}$, typical of the tokamak core. Using the arbitrariness of the test function $\check{\phi}_{1}$ in Eqs. 224 and 7.2, Eq. 223 leads to the gyrokinetic Poisson’s equation GKPII, that we evaluate at $\overline{\bm{R}}$,

[TABLE]

In Eq. 230, the gyroaveraged gyrocenter charge density $\left<\overline{\varrho}_{a}^{*}\right>$ are given by

[TABLE]

whereas the polarization charge density $\mathcal{P}_{a}^{\mu*}$ is

[TABLE]

with $\bm{K}=\bm{k}+\bm{k}^{\prime}$. We remark that the expressions of the gyroavereaged gyrocenter and polarization charge densities can be obtained by using the pull-back transformation of the gyrocenter distribution function [see Eq. 40]. Indeed, the charge particle density $\varrho_{a}=q_{a}\int d\bm{v}f_{a}$, can be written as

[TABLE]

where we use that $f_{a}(\bm{z}(\bm{Z}))=T_{\epsilon_{\delta}}\overline{F_{a}}(\bm{Z})\simeq[1+\overline{g}^{\mu}_{1}\partial_{\overline{\mu}}]\left<\overline{F_{a}}\right>_{\overline{\bm{R}}}+O(\epsilon_{\delta}^{2})$ neglecting the $\overline{\bm{g}}^{\bm{R}}_{1}\cdot\nabla\overline{F_{a}}$ and $\overline{g}^{\parallel}_{1}\partial_{v_{\parallel}}\overline{F_{a}}$ terms since they can be shown to be higher order (Hahm et al., 2009). We identify $\left<\overline{\varrho}_{a}^{*}\right>$ and $\mathcal{P}_{a}^{*\mu}$ as the first and second term in the integral of Section 7.2, respectively. This transformation of the distribution function in Section 7.2 is a common technique to introduce the gyrocenter distribution function in the field equations in the gyrokinetic literature (Dubin et al., 1983; Hahm, 1988, 1996; Hahm et al., 2009).

The Hermite-Laguerre projection of GKPII is obtained by performing the velocity integrals in Eqs. 231 and 7.2. Using the Jacobi-Anger identity in Eq. 138 and expanding the Bessel functions in Laguerre polynomials using Eqs. 141, 155 and 156 yields the gyro-moment expansions of the gyroaveraged gyrocenter and polarization charge densities, respectively,

[TABLE]

and

[TABLE]

where we have introduced the Fourier decomposition of $\left<\Phi_{1}\right>_{\overline{\bm{R}}}$ [see Eq. 143]. Here, the generalized gyrocenter parallel current $\overline{J}_{\parallel a}^{*lm}$ is defined by $\overline{J}_{\parallel a}^{*lm}=q_{a}\overline{N}_{a}\left\lVert\overline{v}_{\parallel}^{l}\right\rVert^{*0m}_{a}$.

We now consider the long wavelength limit approximation of $\mathcal{P}_{a}^{\mu*}$, and show that the expression for the linear polarization charge density commonly used in, e.g., gyrokinetic simulations of edge plasma dynamics (see, e.g., Pan et al., 2016; Hakim et al., 2016), is retrieved. Indeed, keeping terms up to $O(\epsilon_{\perp}^{2})$ in Eq. 216, we derive

[TABLE]

Equation 236 corresponds to the linear polarization charge density derived in Hahm et al. (2009) valid for arbitrary distribution function in the long wavelength limit.

We now consider the drift-kinetic limit of GKPII. For this purpose, we neglect the $O(\epsilon_{\delta})$ terms and consider a long wavelength limit approximation of $\left<\overline{\varrho}^{*}\right>$ in Eq. 231. We first notice that $\mathcal{P}_{a}^{\mu*}$, Section 7.2, vanishes since it is $O(\epsilon_{\delta})$. Expanding the lowest-order kernel up to $O(\epsilon^{2})$ with $b_{a}\sim\epsilon$ in Eq. 231, in particular observing that $\mathcal{K}_{0}(b_{a}^{\dagger})\simeq 1-(b_{a}^{\dagger})^{2}/4$, and performing the inverse Fourier transform, GKPII reduces to Eq. 222. Thus, in the drift-kinetic limit, up to $\epsilon^{2}$, GKPI and GKPII are equivalent.

As a final remark, we note that, to solve the field equations, a plasma vorticity equation can be used, which generalises the ones commonly solved in drift and long wavelength models (Zeiler et al., 1997; Ricci et al., 2012; Jorge et al., 2017; Abel & Hallenbert, 2018). The plasma vorticity equation can be obtained from the quasineutrality condition in Eq. 214. Substituting the time derivative of the quasineutrality condition obtained from Eq. 214 into the fluid equation of the gyrocenter density, given in Eq. 124, we obtain the plasma vorticity equation

[TABLE]

where $\bm{U}_{0a}=\bm{E}\times\bm{B}/B^{2}+\overline{u}_{\parallel a}\bm{b}$ (with the electric field $\bm{E}=-\overline{\nabla}\phi_{0}-\partial_{t}\bm{A}$) and $\bm{J}_{a}^{*}=q_{a}\overline{N}_{a}\bm{u}_{a}^{0}$ is the gyrocenter current density. By comparing Section 7.2 with the drift-kinetic vorticity equation found in Jorge et al. (2017), obtained in a small mass ratio approximation, we remark the presence of collisional gyro-diffusion associated with FLR effects in the collision operator that vanish in the drift-kinetic limit. Also, the last two terms in Section 7.2 are not present in Jorge et al. (2017). These are related to the difference between gyro-center and particle densities. Finally, we remark the presence of additional currents driven by higher order drifts [see Section 5.3]. Section 7.2 should be coupled to Eq. 230 that generalises the gyrokinetic quasineutrality condition by including the presence of $\phi_{0}$ polarization effects.

7.3 Gyrokinetic Ampere’s Laws

From the variation of the action $\delta\mathcal{A}$ in Eq. 209, we derive two gyrokinetic Ampere’s laws, referred to as GKAI and GKAII. GKAI is obtained from the variation

[TABLE]

whereas GKAII is deduced from

[TABLE]

We first consider GKAI. Noticing the presence of $A_{\parallel 0}$ in the inductive part of the electric field, the variation of the field functional action $\mathcal{A}_{f}$, Eq. 205, is

[TABLE]

where the test function is evaluated at $\bm{x}$, i.e. $\check{A}_{\parallel 0}=\check{A}_{\parallel 0}(\bm{x})$. In Section 7.3, we identify the first term as the one contributing the field equation of $A_{\parallel 0}$, whereas the two last terms are boundary terms.

The variation of the gyrocenter functional action $\mathcal{A}_{p}$ in Eq. 206 with respect to $A_{\parallel 0}$ is evaluated by noticing that $A_{\parallel 0}$ is contained in the symplectic components $q_{a}\bm{A}^{*}$, in the terms proportional to the magnetic field strength $B$, such as $\overline{\mu}B$, and in the magnetic vector $\bm{B}$, appearing, e.g., in $\bm{u}_{E}$ present in the gyrocenter Lagrangian $\overline{L}_{a}$, Eq. 207. We remark that, while the gyrokinetic potential $\Psi_{1}$ also depends on $A_{\parallel 0}$ through $\bm{b}$ and the factors proportional to $1/\Omega_{a}$ at second order in $\epsilon_{\delta}$ [see Section 3.4], we can neglect their contributions when evaluating the variation since they are of higher order as $\bm{b}\simeq\hat{\bm{b}}+O(\epsilon)$ and $B\simeq\hat{B}+O(\epsilon)$. The variation of $\mathcal{A}_{p}$ is then

[TABLE]

where $\check{\bm{B}}=\nabla\times\check{\bm{A}}$ with $\check{\bm{A}}=\check{A}_{\parallel 0}\hat{\bm{b}}$, and $\check{\bm{B}}_{\perp}=\bm{b}\times\left(\check{\bm{B}}\times\bm{b}\right)$. In Section 7.3, we have introduced $\check{A}_{\parallel 0}(\overline{\bm{R}})=\int d\bm{x}\delta(\overline{\bm{R}}-\bm{x})\check{A}_{\parallel 0}(\bm{x})$, consistently with Section 7.3 where $\check{A}_{\parallel 0}$ is evaluated at $\bm{x}$. To evaluate the terms that provide dynamics of $A_{\parallel 0}$ in Section 7.3, we carry out successive integrations by parts of the terms where derivatives of $\check{A}_{\parallel 0}$ appear. From Eq. 238 with the variations given in Section 7.3 and Section 7.3, we then obtain GKAI, that we evaluate at $\overline{\bm{R}}$, that is

[TABLE]

In Section 7.3, the gyrocenter current density $\overline{\bm{J}}_{a}^{*}$ is given by

[TABLE]

whereas the magnetization current densities, due to the particle and gyrocenter difference, are due to the classical magnetization $\bm{\mathcal{M}}_{a}^{\mu*}$, i.e.

[TABLE]

to the magnetization $\bm{\mathcal{M}}_{a}^{*}$ associated with the polarization charge density $\bm{\mathcal{P}}_{a}^{*}$ [see Eq. 216],

[TABLE]

to the magnetization $\bm{\mathcal{M}}_{a}^{\parallel*}$,

[TABLE]

and, finally, to the Banõs magnetization $\bm{\mathcal{M}}_{a}^{B*}$,

[TABLE]

We remark that the classical magnetization $\bm{\mathcal{M}}_{a}^{\mu*}$ in Eq. 244 [also referred to as the lowest-order intrinsic guiding-center magnetic dipole contribution, see Brizard (2013)] is associated with the perpendicular gyrocenter pressure. Indeed, at the leading order, we have that $\bm{\mathcal{M}}_{a}^{\mu*}\simeq-\overline{N}_{a}^{*}\overline{P}_{\perp a}\bm{B}/B^{2}$. The second term in $\bm{\mathcal{M}}_{a}^{\mu*}$, proportional to $-\overline{\mu}(\overline{\nabla}\times\bm{U})_{\perp}/(2\Omega_{a})$ [see Eq. 244], is a $O(\epsilon^{2})$ correction to the classical magnetization. Similarly to the gyrokinetic Poisson’s equation GKPI in Eq. 214, the magnetization $\bm{\mathcal{M}}_{a}^{*}$, Eq. 245, is due to the fact that the particle velocity $v_{\perp}$ is written in the frame moving with the $\bm{E}\times\bm{B}$ drift. Indeed, this magnetization effect corresponds to the effective magnetization resulting from the $\bm{\mathcal{P}}_{a}^{*}$ polarization (see Jackson, 2012, Eq. (6.100)). The magnetization current $\bm{\mathcal{M}}_{a}^{\parallel*}$, Eq. 246, represents the correction to the magnetic moment $\overline{\mu}$, resulting from the difference between the gyrocenter and $\bm{E}\times\bm{B}$ drifts, i.e. $\dot{\overline{\bm{R}}}-\bm{u}_{E}$. (Pfirsch & Morrison, 1985; Kaufman, 1986). The magnetization $\bm{\mathcal{M}}_{a}^{B*}$ in Eq. 247 is a $O(\epsilon^{2})$ term associated with the presence of the Banõs drift. Finally, we remark the presence of the polarization current, $\partial_{t}\bm{\mathcal{P}}_{a}^{*}$, which originates from the variation of the inductive part of the $\bm{E}\times\bm{B}$ drift (Brizard & Hahm, 2007; Madsen, 2010; Brizard, 2013).

We now obtain the Hermite-Laguerre representation of GKAI by performing the velocity integrals appearing in Eqs. 243, 244, 246 and 247. This yields

[TABLE]

where the Hermite-Laguerre projection of the magnetization $\bm{\mathcal{M}}_{a}^{*}$ is obtained from Sections 7.2 and 221.

The second gyrokinetic Ampere’s law, GKAII, follows from the variation of the action $\mathcal{A}$ with respect to $A_{\parallel 1}$ in Eq. 239. We first notice that the variation of the field functional action, $\mathcal{A}_{f}$, has the same functional form as Section 7.3, i.e.

[TABLE]

with the test function $\check{A}_{\parallel 1}$ evaluated at $\bm{x}$. In the evaluation of the variation of $\mathcal{A}_{p}$, we note that, as a result of our choice to consider the Hamiltonian gyrokinetic formalism, the symplectic components, $q_{a}\overline{\bm{A}^{*}}$ [see Eq. 56], are independent of $A_{\parallel 1}$. Indeed, only the gyrokinetic potential $\Psi_{1}$, Section 3.4, gives a contribution in the variation of $\mathcal{A}_{f}$. Therefore, the variation of the gyrocenter functional action $\mathcal{A}_{p}$ is given by

[TABLE]

where we introduced the adjoint gyroaverage operator, $\left<\cdot\right>_{\bm{x}}^{\dagger}$ defined in Section 7.2, with $\check{A}_{\parallel 1}(\overline{\bm{R}}+\overline{\bm{\rho}}_{a})=\int d\bm{x}\delta(\overline{\bm{R}}+\overline{\bm{\rho}}_{a}-\bm{x})\check{A}_{\parallel 1}(\bm{x})$. Following the similar steps and assumptions to the ones considered for the evaluation of Section 7.2, we obtain GKAII, which we evaluate at $\overline{\bm{R}}$,

[TABLE]

where the parallel gyroaveraged gyrocenter current density is

[TABLE]

while the parallel magnetization current densities are, respectively,

[TABLE]

and

[TABLE]

The Hermite-Laguerre representations of $\left<\overline{J}_{\parallel a}^{*}\right>$, $\overline{J}_{\parallel a}^{\parallel*}$, and $\overline{J}_{\parallel a}^{\mu*}$ are obtained by using the Jacobi-Anger identity in Eq. 138 and by performing the velocity integrals in Eqs. 253, 253 and 7.3. This yields, respectively,

[TABLE]

As a check of our expressions, we now investigate the long wavelength limit of Eq. 255, i.e. the current densities appearing in GKAII. In particular, considering $b_{a}\sim\epsilon_{\perp}\sim\epsilon$, expanding the lowest-order kernel $\mathcal{K}_{0}=1-(b_{a}^{\dagger})^{2}/4$ up to $O(\epsilon^{2})$, and performing the inverse Fourier transform, we derive

[TABLE]

with $\overline{J}_{\parallel a}=q_{a}\overline{N}_{a}\overline{u}_{\parallel a}$ and $\overline{J}_{th\parallel a}=q_{a}\overline{N}_{a}\overline{v}_{th\parallel a}/\sqrt{2}$, while $\overline{J}_{\parallel a}^{\mu*}$ in Eq. 255c reduces to

[TABLE]

Equation 257 corresponds to the Hermite-Laguerre projection of the magnetization current densities obtained by Hahm et al. (2009) in the long wavelength limit and valid for arbitrary distribution functions. Neglecting the $O(\epsilon_{\delta})$ terms in Section 7.3, while keeping term up to $O(\epsilon_{\perp}^{2})$, one obtains the drift-kinetic Ampere’s law, which can be used to provide an electromagnetic extension of the drift-kinetic moment hierarchy derived in Jorge et al. (2017).

7.4 Equilibrium Pressure Balance Equation

Since a MHD-like equilibrium pressure balance is an important element especially in the study of gradient driven modes (Rogers et al., 2018), we show how the spatial and temporal evolution of the large scale magnetic field $\hat{\bm{B}}$ can be determined such that, at leading order, it reduces to an equilibrium pressure balance. The field equation that sets the self-consistent evolution of $\hat{\bm{B}}$ can be deduced from

[TABLE]

that is

[TABLE]

To show that Section 7.4 reduces its pressure balance equation, we consider a leading order approximation by neglecting the $O(\epsilon_{\delta})$ and $O(\epsilon^{2})$ terms. Moreover, the magnetic fluctuations, $\delta\bm{B}$, are ignored. More precisely, we use the fact that $\left|\delta\bm{B}\right|/\hat{B}\sim\epsilon$ [see Eq. 12], which implies that $\bm{B}+\bm{B}_{1}=\hat{\bm{B}}+O(\epsilon)$. Additionally, the plasma is assumed to be quasi-neutral, i.e. $\sum_{a}q_{a}\overline{N}_{a}=0$. Keeping only the leading order classical magnetization, $\bm{\mathcal{M}}_{a}^{*\mu}=-\overline{P}_{\perp a}\hat{\bm{B}}/\hat{B}^{2}+O(\epsilon^{2})$, while neglecting $\partial_{t}\bm{E}$ and $\partial_{t}\bm{\mathcal{P}}_{a}^{*}$ being higher order terms and solving for $\overline{\nabla}\times\hat{\bm{B}}/4\pi$ in Section 7.4, we obtain

[TABLE]

where the leading order gyrocenter current is $\overline{\bm{J}}^{*}_{a}=q_{a}\overline{N}_{a}\left\lVert\dot{\overline{\bm{R}}}\right\rVert^{*}_{a}$, with the leading order gyro-moment expansion

[TABLE]

as deduced from Eq. 113. Then, the leading order pressure balance is obtained by taking the cross product of Eq. 260 with $\hat{\bm{B}}$ and using the quasi-neutrality condition,

[TABLE]

where $\overline{P}_{\parallel}=\sum_{a}\overline{P}_{\parallel a}$ and $\overline{P}_{\perp}=\sum_{a}\overline{P}_{\perp a}$ are the parallel and perpendicular total pressures. We remark that Eq. 262, which contains the anisotropic pressure terms (see, e.g., Chew et al., 1956; Lanthaler et al., 2019), reduces to the MHD equilibrium pressure balance, $(\overline{\nabla}\times\hat{\bm{B}})\times\hat{\bm{B}}=4\pi\overline{\nabla}_{\perp}\overline{P}$, in the static and isotropic limit, i.e. $\overline{P}_{\parallel}=\overline{P}_{\perp}=\overline{P}$.

We now show that the low-$\beta$ approximation is consistent with the leading order equilibrium pressure balance equation, Eq. 262. Indeed, balancing the $\left|\overline{\nabla}\hat{B}^{2}\right|\sim B^{2}/L_{B}$ and $\left|\nabla_{\perp}\overline{P}_{\perp}\right|\sim 1/L_{P}\overline{P}_{\perp}$ terms yields,

[TABLE]

since we assumed $\rho_{s}/L_{P}\sim\epsilon$ [see Eq. 3] and $\epsilon_{B}\sim\epsilon^{3}$ [see Eq. 4]. Thus, steep pressure gradients on scale length $L_{P}$ is consistent with a low-$\beta$ plasma in the presence of an equilibrium magnetic field varying on large scales.

8 Conclusion

In the present work, a gyrokinetic model is derived to evolve the turbulent plasma dynamics in the periphery of tokamak devices. This model takes the form of an infinite set of coupled fluid equations, the gyro-moment hierarchy equation given in Eq. 110 with a nonlinear gyrokinetic Dougherty collision operator in Eq. 191, coupled with a set of gyrokinetic field equations, that are two Poisson’s equations, GKPI in Eq. 214 and GKPII in Eq. 230, and two Ampere’s laws, GKAI in Section 7.3 and GKAII in Eq. 251. The model is based on second order fully electromagnetic gyrokinetic equations of motion of a charged particle. The equations are obtained by using Lie-transform perturbation theory. More precisely, by taking advantage of the scale separation between the equilibrium scales and the particle gyroscales, two changes of phase-space coordinates are performed allowing the description of the single particle motion in the presence of large and small-scale electromagnetic fluctuations and strong flows. Then, the collective behaviour is introduced by deriving a second order accurate gyrokinetic Boltzmann equation. The gyrokinetic equation is further developed into a gyro-moment hierarchy valid for far from equilibrium distribution functions, obtained by projecting the gyrokinetic Boltzmann equation onto a complete set of Hermite-Laguerre polynomials. These polynomials constitute a complete velocity-space basis used to conveniently expand the gyroaveraged distribution function. The result of the projection is an infinite set of coupled fluid equations for the temporal and spatial evolution of the Hermite-Laguerre expansion coefficients, that are referred to as gyro-moments. In the process, the linear and nonlinear coupling between the gyro-moments, associated with parallel streaming along the magnetic field lines, magnetic field gradients, and FLR effects, are analytically treated and retained at arbitrary values of the perpendicular wavenumber. In particular, a closed form of the gyroaveraged operator is derived in terms of gyro-moments. The effects of collisions in the plasma periphery dynamics are introduced by considering a gyrokinetic Dougherty collision operator. This operator is expressed as a function of the gyro-moments, is accurate at arbitrary values of the perpendicular wavenumber, and is nonlinear. We remark that the gyro-moment expansion of a full nonlinear gyrokinetic Coulomb collision operator can be derived within the same formalism as in Section 6 to describe efficiently like and unlike species collisions in more general situations. This will be subject of a future work. Finally, a set of self-consistent gyrokinetic Poisson-Ampere equations are obtained from a variational principle. More precisely, two coupled gyrokinetic Poisson’s equations and two coupled gyrokinetic Ampere’s laws, with an additional Ampere’s law needed to set the evolution of the large scale magnetic field, are derived to provide the necessary closure of the gyro-moment hierarchy by determining the different components of the fluctuating electromagnetic fields. Within our variational formulation, polarization and magnetization corrections, which are related to the change of phase-space coordinates at the particle Lagrangian level, appear self-consistently. The charge and current densities are expressed in terms of the gyro-moments by performing analytically the velocity integrals at arbitrary wavelengths. Thus, a complete generalization of the analytical expressions of the polarization charge and magnetization current densities, present in earlier gyrokinetic derivations, is obtained.

Table 1 illustrates and summarizes the improvements and differences between the present model and previous gyrokinetic theories. In particular, we point out that the present gyro-moment hierarchy is an extension of the drift-kinetic moment hierarchy for the SOL dynamics at arbitrary collisionality developed by Jorge et al. (2017) to fully electromagnetic fluctuations allowing for perpendicular wavenumbers of the order of the ion sound Larmor radius (see Table 1).

Although a numerical implementation of the herein model is outside of the scope of the present work, we remark that closure schemes need to be applied to the gyro-moment hierarchy equation, i.e. a truncation of the series in Eq. 91 is required. A possible closure is based on the truncation of all the gyro-moments $\overline{N}_{a}^{lk}$ of order higher than a given order. Then, the unresolved gyro-moments are comparable to the limitations of velocity-space grid methods to resolve fine velocity structures. However, truncations can lead to poor results especially at low Hermite-Laguerre resolution (Mandell et al., 2018). An exact asymptotic closure, known as semi-collisional closure (Zocco & Schekochihin, 2011; Loureiro et al., 2016), can be rigorously applied even at small (but finite) collisionality. With such a closure, the present hierarchy provides an ideal framework for the simulation of the dynamics of the plasma periphery. While being a rigorous asymptotic limit of the exact gyrokinetic equation [see Eq. 84], the model herein present a tuneable kinetic accuracy which depends on the number of gyro-moments retained.

Acknowledgements

The authors acknowledge helpful discussions with A. Baillod, A. J. Brizard, S. Brunner, W. Dorland, N. F. Loureiro, A. A. Schekochihin, P. B. Snyder, M. Held and N. Tronko. This research was supported in part by the Swiss National Science Foundation, and has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014 - 2018 and 2019 - 2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A Guiding-Center Transformation

This appendix reports on the details of the second order guiding-center transformation whose results are used in Section 3.3. Using Eq. 44c and computing the Lie-derivatives of $\gamma_{1}$ and $\Gamma_{1}$ given in Eqs. 48b and 3.3 with the definition in Eq. 39, we obtain the $O(\epsilon^{2})$ guiding-center correction $\Gamma_{2}$, expressed in the guiding-center coordinates $\bm{Z}$,

[TABLE]

To remove one inherent degree of freedom, we set $S_{2}=0$ and choose $\bm{g}^{\bm{R}}_{1}\cdot\bm{b}=0$. Therefore, $\bm{g}^{\bm{R}}_{1}$ is purely perpendicular to $\bm{b}$ and is given by Eq. 51. This results into $\Gamma_{2\theta}=B\mu/\Omega$, $\Gamma_{2\mu}=0$ and $\Gamma_{2\parallel}=0$. Taking into account that fact that $\bm{c}^{\prime}_{\perp}$ depends on $B$ (indeed, $v_{\perp}=\sqrt{2\mu B/m}$) results the identity

[TABLE]

with $\bm{T}=(\nabla\bm{c})\cdot\bm{a}=\nabla\bm{e}_{2}\cdot\bm{e}_{1}$ being the first Littlejohn’s gyrogauge field vector (Littlejohn, 1988). With $\bm{g}^{\bm{R}}_{1}\times\bm{f}=\left(\bm{g}^{\bm{R}}_{1}\times\bm{b}\right)\left(\bm{b}\cdot\bm{f}\right)+\bm{b}\left(\bm{c}^{\prime}_{\perp}\cdot\bm{f}\right)/\Omega$ where $\bm{f}$ is any arbitrary vector function, and introducing the quantities $h^{\theta}$, $h^{\mu}$ and $h^{\parallel}$ defined by (Brizard, 1995)

[TABLE]

into Appendix A, the perpendicular and parallel components to $\bm{b}$ of the symplectic components $\Gamma_{2\bm{R}}$ can be separated as

[TABLE]

Introducing $h^{\parallel}=\left<h^{\parallel}\right>+\widetilde{h^{\parallel}}$ and $\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}=\left<\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}\right>+\widetilde{\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}}$ where $\left<\chi\right>\equiv\int_{0}^{2\pi}d\theta\chi$ and $\widetilde{\chi}=\chi-\left<\chi\right>$ being the gyrophase dependent part of any gyrophase dependent function $\chi=\chi(\theta)$, we cancel the gyrophase dependent terms parallel to $\bm{b}$ by setting

[TABLE]

and in the perpendicular direction by taking the vector product of Appendix A with $\bm{b}$ to obtain an expression for $\bm{g}_{2\perp}^{\bm{R}}$,

[TABLE]

Thus, $\Gamma_{2\bm{R}}$ reduces to

[TABLE]

At this stage, we have the freedom to choose $\left<h^{\parallel}\right>=-\left<\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}\right>/2$. This particular choice of $\left<h^{\parallel}\right>$ has the effect of transferring the latter term in the Hamiltonian component. An alternative possibility is to set $\left<h^{\parallel}\right>=0$, and so keeping the $m\left<\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}\right>/2$ in the symplectic component (Brizard, 1995; Madsen, 2010). Thus, at second order, we find

[TABLE]

while, $\left<h^{\parallel}\right>=0$ yields $\Gamma_{2\bm{R}}=-B\mu\bm{T}^{*}/\Omega$ where $\bm{T}^{*}=\bm{T}-\bm{b}\bm{b}\cdot\nabla\times\bm{b}/2$.

The Hamiltonian component of $\Gamma_{2}$, namely $\Gamma_{2t}$, can be developed as

[TABLE]

Using the expression for $\bm{g}^{\bm{R}}_{2\perp}$ given in Eq. 269 into Appendix A, terms combine, and $\Gamma_{2t}$ reduces to

[TABLE]

Introducing $h^{\mu}=\left<h^{\mu}\right>+\widetilde{h^{\mu}}$, $h^{\parallel}=\left<h^{\parallel}\right>+\widetilde{h^{\parallel}}$, and $\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}=\left<\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}\right>+\widetilde{\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}}$ in $\Gamma_{2t}$, allows to separate the gyrophase dependent and independent parts by noticing that

[TABLE]

where $S=\bm{a}\cdot\partial_{t}\bm{c}$ being the second gyrogauge field vector introduced by Littlejohn (1988), such that

[TABLE]

where the gyrophase dependent term are cancelled by setting

[TABLE]

Finally, with $\left<h^{\parallel}\right>=-\left<\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{b}\cdot\bm{c}^{\prime}_{\perp}\right>/2$ into Eq. 275, the second order guiding-center correction $\Gamma_{2}$ is given by

[TABLE]

The analytical expression for $\left<h^{\mu}\right>$ in $\Gamma_{2}$ given in Eq. 277 must be obtained by considering a third-order calculation that ensures $\dot{\mu}=O(\epsilon^{4})$. A sufficient condition is to impose that $\Gamma_{3\theta}=0$, $\Gamma_{3\mu}=0$ and $\Gamma_{3\parallel}=0$. As a results, we also obtain the closed expressions for $\bm{b}\cdot\bm{g}^{\bm{R}}_{2}$ and $h^{\theta}$. From Eq. 44d, we remark that only the Lie-derivatives $\mathcal{L}_{1}^{2}\gamma_{1}$, $\mathcal{L}_{1}^{2}\Gamma_{1}$ and $\mathcal{L}_{2}\Gamma_{1}$ contribute in $\Gamma_{3\theta}$, $\Gamma_{3\mu}$ and $\Gamma_{3\parallel}$. In fact, while $\gamma_{3}$ and $\gamma_{2}$ are identically zero [see Eqs. 48a and 48b], $\mathcal{L}_{3}\gamma_{0}$ has only a $\bm{R}$ and a $t$ component. Thus, we derive that

[TABLE]

and

[TABLE]

while

[TABLE]

By explicit evaluation of the different terms, and introducing $h^{\mu}$ given in Eq. 266a, we derive

[TABLE]

Setting $\Gamma_{3\theta}=0$ yields

[TABLE]

and

[TABLE]

where we used the identity $\left<\bm{g}^{\bm{R}}_{1}\cdot\nabla\bm{U}\cdot\bm{c}^{\prime}_{\perp}\right>=\mu\bm{b}\cdot\nabla\times\bm{U}/2$. An closed expression for $h^{\theta}$ can be derived from $\Gamma_{3\mu}=0$, i.e.

[TABLE]

while from $\Gamma_{3\parallel}$ one obtains

[TABLE]

Thus, with Eq. 282, the second order accurate guiding-center one-form, $\Gamma=\Gamma_{0}+\Gamma_{1}+\Gamma_{2}$, is given by

[TABLE]

with $\bm{A}^{*}=\bm{A}+m\bm{U}/q$. We remark that the Baños term, proportional to $\mu v_{\parallel}B\bm{b}\cdot\nabla\times\bm{b}/2\Omega$ contained in the last term of Appendix A, can be removed from the Hamiltonian component with the choice $\left<h^{\parallel}\right>=0$. Indeed, one finds

[TABLE]

where $\bm{T}^{*}=\bm{T}+\bm{b}\bm{b}\cdot\nabla\times\bm{b}/2$ (Brizard, 1995).

Appendix B Gyrocenter Transformation

We describe the second order gyrocenter transformation in $\epsilon_{\delta}$ whose results are presented in Section 3.4. We start to derive $\overline{\Gamma}_{1}$ expressed in the gyrocenter coordinates $\overline{\bm{Z}}$. Within the choice of a Hamiltonian formulation, the symplectic components of $\overline{\Gamma}_{1}$ vanish. Therefore, from Eq. 44b, we obtain that

[TABLE]

The first order gyrocenter correction, $\overline{\Gamma}_{1}=\overline{\Gamma}_{0}-\mathcal{L}_{1}\delta\Gamma$, reduces to

[TABLE]

The closed analytical expressions of the first order generating functions $\overline{g}^{\parallel}_{1}$ and $\overline{\bm{g}}^{\bm{R}}_{1}$ can be obtained by taking the scalar and vector product of Eq. 288a with $\bm{B}^{*}$ and $\bm{b}$, respectively, while the one for $\overline{g}^{\mu}_{1}$ and $\overline{g}^{\theta}_{1}$ follows Eqs. 288c and 288d. Thus, we obtain

[TABLE]

With $\bm{A}_{1}=A_{\parallel 1}\hat{\bm{b}}$, Eq. 290 reduces to Eq. 64. We remark that the contribution of $\bm{A}_{1}$ can be transferred to the symplectic components of $\overline{\Gamma}_{1}$ via the transformation $\bm{A}_{1}\to\widetilde{\bm{A}_{1}}$ in $\overline{\bm{g}}^{\bm{R}}_{1}$ and in $\overline{g}^{\parallel}_{1}$ in Eq. 290, yielding $\overline{\Gamma}_{1\overline{\bm{R}}}=q\left<\bm{A}_{1}\right>_{\overline{\bm{R}}}$ (Brizard & Hahm, 2007; Madsen, 2013; Tronko et al., 2016). Isolating the gyrophase dependent terms in the electromagnetic fluctuations by using Eq. 63, approximating $\partial_{\overline{v}_{\parallel}}\overline{\mathcal{H}}_{0}\simeq m\overline{v}_{\parallel}+O(\epsilon^{2})$, $\partial_{\overline{\mu}}\overline{\mathcal{H}}_{0}\simeq B+O(\epsilon^{2})$ and $B_{\parallel}^{*}\simeq B+O(\epsilon^{2})$ and, finally, neglecting the time derivative of the effective magnetic potential, being $\partial_{t}\bm{A}^{*}\sim\epsilon^{2}\Omega\bm{A}^{*}$, the first order gauge function $S_{1}$ satisfies,

[TABLE]

where we introduce the guiding-center convective derivative $d_{tgc}\equiv\partial_{t}+\left(\overline{v}_{\parallel}\bm{b}+\bm{D}_{\perp}\right)\cdot\overline{\nabla}$ with $\bm{D}_{\perp}=\bm{b}\times\overline{\nabla}\;\overline{\mathcal{H}}_{0}/(qB)$. We remark that $\bm{D}_{\perp}$ reduces to $\bm{u}_{E}$ at the leading order. Ordering the terms on the left-hand side of Eq. 291 according to Section 3, i.e. $\Omega\partial_{\overline{\theta}}S_{1}\sim\Omega S_{1}$, $d_{tgc}S_{1}\sim\epsilon\Omega S_{1}$ while $\partial_{\overline{v}_{\parallel}}S_{1}\overline{\nabla}_{\parallel}\overline{\mathcal{H}}_{0}/m\sim\epsilon^{2}\Omega S_{1}$, shows that a solution for $S_{1}$ can be obtained iteratively by expanding $S_{1}=S_{10}+S_{11}+\dots$ with $S_{11}\sim O(\epsilon S_{10})$. We remark that an analytical solution of the first order gauge function, $S_{1}$, is an important element for the second order analysis. With $\bm{A}_{1}=A_{\parallel 1}\bm{b}$, the leading order solution, $S_{10}$, satisfies $\Omega\partial_{\overline{\theta}}S_{10}=q\widetilde{\Phi_{1}}$ where we introduce the first order gyrokinetic potential,

[TABLE]

being $\phi_{1}\sim\overline{v}_{\parallel}A_{\parallel 1}$, yielding

[TABLE]

Then, the first order correction, $S_{11}$, is given by $\Omega\partial_{\overline{\theta}}S_{11}=-d_{tgc}S_{10}$, such that $S_{11}=-d_{tgc}\overline{S_{10}}/\Omega$. A closed analytical expression for $S_{10}$ can be derived, using Eq. 293, performing the $\overline{\theta}$-integral explicitly in Fourier-space. Indeed, expanding the fluctuation $\widetilde{\Phi_{1}}$ in Fourier harmonics, i.e. $\widetilde{\Phi_{1}}(\bm{k})=\sum_{n\neq 0}i^{n}J_{n}(b)\Phi_{1}(\bm{k})e^{in\overline{\theta}}$ (where $\Phi_{1}(\bm{k})$ is the Fourier harmonics of $\Phi_{1}(\overline{\bm{R}})$). Here, $J_{n}(x)$ is the $n$th-order Bessel function [see Section 5.4]. Then, inserting the latter expression for $\widetilde{\Phi_{1}}$ yields

[TABLE]

With Eq. 291, the first order gyrocenter one-form correction $\overline{\Gamma}_{1}$ in Eq. 289 reduces to

[TABLE]

We now perform the second order analysis in $\epsilon_{\delta}$. From Eq. 44c, the second order gyrocenter one-form correction, $\overline{\Gamma}_{2}$, reads

[TABLE]

with

[TABLE]

Applying the gyrocenter transformation rules of our Hamiltonian formulation, we find the second order generating functions,

[TABLE]

Equation 298 reduces to Eq. 69 at the leading-order with $\bm{A}_{1}=A_{\parallel 1}\hat{\bm{b}}$. Using Eq. 298 and the closed expressions for the first order generating functions, given in Eq. 290 with $\bm{A}_{1}=A_{\parallel 1}\hat{\bm{b}}$ and within the same approximations made in the first order analysis, $-\overline{\mathcal{H}}_{2}$ reduces to

[TABLE]

Noticing that $\left<S_{1}\right>_{\overline{\bm{R}}}=0$ and isolating the gyrophase dependent terms, we obtain the expression given in Section 3.4. In particular, we cancel the gyrophase dependent terms by choosing $S_{2}$ such that it satisfies

[TABLE]

We remark that we do not solve explicitly $S_{2}$ from Appendix B. Approximating with $S_{1}\simeq S_{10}$ and neglecting terms proportional to $d_{tgc}\sim\epsilon\Omega$ and parallel gradients in Section 3.4, we derive the leading order second order gyrocenter correction one-form $\overline{\Gamma}_{2}$,

[TABLE]

Appendix C Hermite-Laguerre Expansion of $\left<(\bm{b}\times\overline{\nabla}S_{1})\cdot\overline{\nabla}\widetilde{\Phi_{1}}\right>_{\overline{\bm{R}}}$

This appendix reports the leading order Hermite-Laguerre projection of the last term in Section 5.2. Using the expression of $\Phi_{1}$ given in Eq. 65 with $S_{1}\simeq S_{10}$ in Eq. 293, the considered term can be written as

[TABLE]

In the evaluation of Appendix C, we pursue a procedure similar to the one applied in Section 5.4 to evaluate the nonlinear terms. In particular, we derive the Hermite-projection of the first term in Appendix C which can then be used to evaluate the remaining terms. Using the closed analytical expression of $S_{10}$ given in Eq. 294, expanding $\widetilde{\phi_{1}}$ in Fourier harmonics, applying the gradient operator and, finally, the gyroaverage operator yields

[TABLE]

where we neglected the terms proportional to $|\overline{\nabla}J_{n}(b)|\sim\lvert\overline{\nabla}_{\perp}B/B\rvert\lvert J_{n}^{\prime}(b)\rvert$ in Eq. 303 since they are smaller by a least a factor $\epsilon$. Expressing the Bessel functions in terms of associated Laguerre polynomials using Eq. 141, and applying consecutively the two following identities

[TABLE]

where

[TABLE]

and

[TABLE]

where

[TABLE]

we derive that the Hermite-Laguerre projection of the last term appearing in Section 5.2, i.e.

[TABLE]

with $\overline{K}^{n}_{rsr_{1}s_{1}}=(-1)^{n}e^{in\alpha}\overline{L}_{rr_{1}}^{|n|}d_{sr_{1}s_{1}}^{|n|}r!s!/[(|n|+r)!(|n|+s)!]$, and where the operator $\bm{\overline{\mathcal{D}}}_{anrs}^{pjs_{1}}(b_{a},b_{a}^{\prime},\bm{K})$ is defined in Section 5.4. Appendix C can be generalized to evaluate the remaining terms in Appendix C.