Hierarchy of second-order gyrokinetic Hamiltonian models for particle-in-cell codes

TL;DR

This paper derives and compares second-order gyrokinetic Hamiltonian particle models for particle-in-cell simulations, including full FLR and long-wavelength approximations, enhancing the theoretical foundation for plasma simulation accuracy.

Contribution

It provides explicit derivations of second-order gyrokinetic Hamiltonian models suitable for particle-in-cell codes, including a comparison of full FLR and long-wavelength models.

Findings

Derived the full second-order FLR Hamiltonian model.

Presented the long-wavelength approximation model.

Showed the equivalence of models in the long-wavelength limit.

Abstract

The reduced-particle model is the central element for the systematic derivation of the gyrokinetic Vlasov-Maxwell equations from first principles. Coupled to the fields inside the gyrokinetic field-particle Lagrangian, the reduced-particle model defines polarization and magnetization effects appearing in the gyrokinetic Maxwell equations. It is also used for the reconstruction of the gyrokinetic Vlasov equation from the particle characteristics. Various representations of reduced-particle models are available according to the choice of the gyrokinetic phase space coordinates. In this paper, the Hamiltonian representation of the reduced particle dynamics at an order suitable for the implementation in particle-in-cell simulations is explicitly derived from the general reduction procedure. The second-order (with respect to the fluctuating electromagnetic fields), full Finite Larmor Radius (FLR) Hamiltonian gyrokinetic particle model as well as the second-order model suitable specifically for the long-wavelength approximation (i.e., containing up to the second-order FLR corrections), are derived and compared to the model recently implemented in the particle-in-cell code ORB5. We show that the same long-wavelength approximate equations can also be derived by taking the proper limit of the full FLR model.

Equations

$$\mu =m{v}_{\perp }^2/2B,$$

$$\mathbf{x}=\mathbf{X}+{\mathbf{\rho }}_0(\mathbf{X},\mu ,\theta )+{\mathbf{\rho }}_1(\mathbf{X},\mu ,\theta ),$$

$${\mathbf{\rho }}_0\equiv \frac{mc}{e}\sqrt{\frac{2\mu }{mB}}\hat{\mathbf{\rho }}\equiv {\rho }_0\hat{\mathbf{\rho }}\sim \mathcal{O}({ϵ}_B^0),$$

$${\mathbf{\rho }}_1=-{ϵ}_{\delta }\frac{m{c}^2}{e{B}^2}{\nabla }_{\perp }\left({\varphi }_1(\mathbf{X})-\frac{p_z}{mc}{A}_{1\Vert }(\mathbf{X})\right)\sim \mathcal{O}({ϵ}_{\delta })\mathcal{O}({ϵ}_{\perp }),$$

$${\hat{\mathbf{b}}}_1=\hat{\mathbf{b}}\cdot \nabla \hat{\mathbf{b}}/|\hat{\mathbf{b}}\cdot \nabla \hat{\mathbf{b}}|,$$

$${\hat{\mathbf{b}}}_2=\hat{\mathbf{b}}\times {\hat{\mathbf{b}}}_1.$$

$$\hat{\mathbf{\rho }}={\hat{\mathbf{b}}}_1\cos \theta -{\hat{\mathbf{b}}}_2\sin \theta ,\hat{\perp }=-{\hat{\mathbf{b}}}_1\sin \theta -{\hat{\mathbf{b}}}_2\cos \theta ,$$

$$\mathbf{v}=v_{\Vert }\hat{\mathbf{b}}+\sqrt{\frac{2\mu B}{m}}\hat{\perp }.$$

$$L_{\mathrm{gc}}\left(\mathbf{X},p_{\Vert },\mu ,\theta \right)=\frac{e}{c}{\mathbf{A}}^{\ast }\cdot \dot{\mathbf{X}}+\frac{mc}{e}\mu \dot{\theta }-H_{\mathrm{gc}},$$

$${\mathbf{A}}^{\ast }=\mathbf{A}+\frac{c}{e}{p}_{\Vert }\hat{\mathbf{b}}.$$

$$H_{\mathrm{gc}}=\frac{p_{\Vert }^2}{2m}+\mu B.$$

$$p_z=m{v}_{\Vert }+{ϵ}_{\delta }\frac{e}{c}{A}_{1\Vert }(\mathbf{X}+{\mathbf{\rho }}_0),$$

$$L_{\mathrm{pert}}\left(\mathbf{X},p_z,\mu ,\theta ;t,w\right)$$

$$H_{\mathrm{pert}}\left(\mathbf{X},p_z,\mu ,\theta ;t,w\right)$$

$$\{F,G{\}}_{\mathrm{gc}}$$

$${\left\{F,G\right\}}_{\mathrm{gc}}=\frac{1}{{ϵ}_{\delta }}{\left\{F,G\right\}}_{-1}+{\left\{F,G\right\}}_0+{ϵ}_{\delta }{\left\{F,G\right\}}_1-\frac{\partial F}{\partial w}\frac{\partial G}{\partial t}+\frac{\partial F}{\partial t}\frac{\partial G}{\partial w},$$

$${\left\{F,G\right\}}_{-1}$$

$${\left\{F,G\right\}}_0$$

$${\left\{F,G\right\}}_1$$

$$H_{\mathrm{pert}}=H_0+{ϵ}_{\delta }e{\psi }_1(\mathbf{X}+{\mathbf{\rho }}_0,p_z)+{ϵ}_{\delta }^2\frac{1}{2m}{\left(\frac{e}{c}\right)}^2{A}_{1\Vert }^2(\mathbf{X}+{\mathbf{\rho }}_0),$$

$$H_0=\frac{p_z^2}{2m}+\mu B+w,$$

$${\psi }_1(\mathbf{X}+{\mathbf{\rho }}_0,p_z)={\varphi }_1(\mathbf{X}+{\mathbf{\rho }}_0)-\frac{p_z}{mc}{A}_{1\Vert }(\mathbf{X}+{\mathbf{\rho }}_0),$$

$$\bar{G}={\mathrm{e}}^{-\epsilon {\pounds }_S}G=G-\epsilon \left\{S,G\right\}+\frac{1}{2}{\epsilon }^2\left\{S,\{S,G\}\right\}+\mathcal{O}({\epsilon }^3),$$

$$\bar{\mathbf{z}}={\mathrm{e}}^{\epsilon {\pounds }_S}\mathbf{z},$$

$$\bar{\mathbf{X}}$$

$${\bar{p}}_z$$

$$\bar{\mu }$$

$$\bar{\theta }$$

$${\mathrm{e}}^{-{ϵ}_{\delta }{\pounds }_S}(\mathbf{X}+{\mathbf{\rho }}_0)=\mathbf{X}+{\mathbf{\rho }}_0+{\mathbf{\rho }}_1,$$

$${\mathbf{\rho }}_1=-{ϵ}_{\delta }\{S,\mathbf{X}+{\mathbf{\rho }}_0{\}}_{\mathrm{gc}},$$

$${\mathrm{e}}^{-{ϵ}_{\delta }{\pounds }_S}$$

$$\Psi (\mathbf{X}+{\mathbf{\rho }}_0,p_z,\mu ,\theta )=⟨\Psi (\mathbf{X}+{\mathbf{\rho }}_0,p_z,\mu ,\theta )⟩+\tilde{\Psi }(\mathbf{X}+{\mathbf{\rho }}_0,p_z,\mu ,\theta ),$$

$$⟨\Psi (\mathbf{X}+{\mathbf{\rho }}_0,p_z,\mu ,\theta )⟩=\frac{1}{2\pi }{\int }_0^{2\pi }\Psi (\mathbf{X}+{\mathbf{\rho }}_0,p_z,\mu ,\theta )d\theta .$$

Full text

1 Introduction

A magnetized plasma represents a complex system with multi-scale dynamics in both time and space, which is a challenge for numerical implementation. For several decades the gyrokinetic dynamical reduction [1], [2], [3] has been of interest as one of the most powerful tools to study this multi-scale problem both numerically and analytically (see for example [4],[5]).

One of the main challenges for a systematic derivation of the gyrokinetic theory comes from the fact that the general derivation [3] includes several groups of small parameters. The first group is connected to the relative amplitude of the background inhomogeneities, while the second group originates from the ratio between the amplitude of fluctuating electromagnetic fields and background quantities. The full gyrokinetic derivation considers both kinds of small parameters of the same order, which makes not only the derivation of the final model but also its numerical implementation rather challenging.

Nowadays, a simplified ordering, which allows one to transfer corrections from the background inhomogeneities at the next higher order in perturbation theory, is widely used for the derivation of models implemented in gyrokinetic codes. In particular, the Vlasov-Maxwell reduced models with linearised polarisation and magnetisation are of great interest.

From the perspective of the approximations performed on the electromagnetic field fluctuations, models for global codes (i.e., ORB5 [6] and GENE [7]) are usually derived in the $\textit{low}-\beta$ approximation, meaning that the parallel component of the perturbed magnetic field is systematically neglected, which is the framework we consider in this article.

The gyrokinetic theory derived from the field theory formulation has a great advantage over the direct asymptotic decomposition of the Vlasov equation first derived in Ref. [1]. In fact, considering the coupling between fields and reduced particle dynamics within the same mathematical framework gives access to the derivation of self-consistently coupled gyrokinetic Vlasov-Maxwell equations.

The reduced-particle model is an essential element of this derivation. First of all, the reduced particle dynamics affects the Maxwell equations via polarization and magnetization terms. At the same time, it defines the gyrokinetic Vlasov equation, which can be directly reconstructed from its characteristics. The coupling between the reduced field and particle equations can be systematically established via a first-principle derivation from the gyrokinetic Lagrangian [8], [9].

A systematic variational framework for gyrokinetic theory has undergone a significant development during the last two decades, while the development of codes started a decade earlier. For this historical reason, some of the major gyrokinetic codes are based on the asymptotic derivation rather than on field theory. A significant effort toward code verification started in 2014 in the framework of the European project VeriGyro. In a recent work, [10] orderings have been identified and a general gyrokinetic field theory has been developed for an explicit derivation of the model implemented in the PIC code ORB5.

This paper focusses on the detailed derivation of Hamiltonian gyrokinetic models for the reduced particle dynamics suitable for the derivation of the gyrokinetic Vlasov-Maxwell system. Such a derivation has been recently presented in a compact and simplified form for the long-wavelength approximation [10]. Here, we extend this derivation by including the full Finite Larmor Radius (FLR) effects in the second-order Hamiltonian model. We show that the very same long-wavelength approximated equations can be retrieved by taking the proper limit of the full FLR model.

This paper is organised as follows. In Sec. 2, we summarize the main idea of the gyrokinetic reduction. In Sec. 3, we set up the general framework for coordinate change and the derivation of the reduced dynamics. In Sec. 4, we present the derivation of the Hamiltonian gyrokinetic models with the full series of FLR corrections and then in Sec. 5, in the long-wavelength approximation. Finally, in Sec. 6, we give expressions for the corresponding characteristics and provide the corresponding gyrokinetic Vlasov equations. The Vlasov-Maxwell gyrokinetic models, derived from the full FLR second-order gyrokinetic particle model and the model truncated up to the second-order FLR corrections suitable for the long-wavelength approximation are currently implemented in the PIC code ORB5 [6].

2 Gyrokinetic dynamical reduction in a nutshell

In magnetised plasmas, the presence of a strong magnetic field induces a separation of the scales of motion. The particle motion is decomposed into a fast rotation around the magnetic field lines and a slow drift motion in the direction perpendicular to the magnetic field. The scale of gyromotion is set by the cyclotron frequency $\Omega=eB/mc$, where $e$ and $m$ are, respectively, the charge and mass of the particles, $B$ is the magnetic field amplitude and $c$ is the speed of light. The gyromotion is described by a fast gyroangle variable $\theta$ to which corresponds a canonically conjugated, slowly varying magnetic moment $\mu$. At the lowest order, it is given by

[TABLE]

where $v_{\perp}$ is the perpendicular velocity of particle with respect to the magnetic field lines. In slab (constant magnetic field) geometry, $\mu$ is an exact dynamical invariant. However, magnetic curvature effects as well as the presence of electromagnetic fluctuations break this exact invariance. The gyrokinetic dynamical reduction uses the fact that, averaged over long times, the magnetic moment is conserved, i.e., $\left\langle\dot{\mu}\right\rangle_{t}=0$.

The goal of the gyrokinetic dynamical reduction consists in building up a new set of phase space variables, such that the $\theta$ dependence is completely uncoupled and $\mu$ has a trivial dynamics, i.e., $\dot{\mu}=0$. Therefore, the reduced particle dynamics on the $4+1$ dimensional phase space with variables $\left(\mathbf{X},p,\mu\right)$, where $\mathbf{X}$ represents the reduced position, $p$ is a scalar momentum coordinate and the new magnetic moment $\mu$ is constant. This change of coordinates is constructed via a perturbative series of near-identity phase space transformations, i.e., these transformations are invertible at each step of the perturbative procedure. The reduced position $\bf X$ has a simple geometrical meaning: It represents an instantaneous center of the fast particle rotation around the magnetic field line. Therefore, from the space coordinate point of view the gyrokinetic transformation is a shift between the initial particle coordinate $\bf x$ and the instantaneous center of its rotation $\bf X$. The difference between both positions is the polarization displacement $\bm{\rho}$, the derivation of which will be discussed later. Performing numerical simulations on the $4+1$ dimensional phase space instead of the $6$ dimensional one and also removal of the fastest time scales enable the drastic reduction of the numerical cost.

The dynamical reduction can be organized in one or two steps. In the framework of the one-step procedure, the contributions from the background geometry non-uniformities and electromagnetic fluctuations to the breaking of the magnetic moment conservation are taken into account simultaneously. The two-step procedure, composed by the guiding-center reduction and the subsequent gyrocenter reduction, allows one to treat those effects separately, which may have some advantages for making a direct link between the coordinate transformation and the polarization effects induced on the reduced particle and field dynamics. Here, we consider the two-step procedure in order to make a clear separation between the polarization contributions associated with each of these transformations at the lowest order.

In the two-step reduction procedure, a small parameter is associated with each transformation: These parameters are for the guiding-center $\epsilon_{B}=\rho_{th}/L_{B}$, where $\rho_{th}$ is the thermal Larmor radius of particle and $L_{B}=\left|\bm{\nabla}B/B\right|^{-1}$ sets the spatial scale for background magnetic field variation, and for the gyrocenter, $\epsilon_{\delta}=\left(k_{\perp}\rho_{th}\right)\ e\phi_{1}/T\equiv\epsilon_{\perp}e\phi_{1}/T$. The dimensionless parameter $\epsilon_{\perp}$ allows one to distinguish between the gyrokinetic theory with $\epsilon_{\perp}\sim\mathcal{O}(1)$ and the drift-kinetic theory with $\epsilon_{\perp}\ll 1$, known also as the * long-wavelength* approximation when only the second-order $\mathcal{O}(\epsilon_{\perp}^{2})$ corrections are included. Following the gyrokinetic ordering relevant for numerical implementation, we consider $\epsilon_{B}\ll\epsilon_{\delta}$, i.e., all the background gradient effects are of higher order with respect to the amplitude of the fluctuations.

In what concerns the FLR or the $\epsilon_{\perp}$- ordering, we consider models derived in the limit with full FLR corrections as well as models truncated up to the second-order in $\epsilon_{\perp}$.

As a perturbative theory, each of the coordinate transformations, the guiding-center and the gyrocenter, represents an infinite series of corrections ordered according to the corresponding small parameter, $\epsilon_{B}$ or $\epsilon_{\delta}$. In particular, for the particle position $\bf x$, this means that the exact gyrokinetic coordinate transformation contains an infinite series of polarization displacements. Roughly speaking, the position of the particle $\bf x$ as a function of the reduced gyrocenter position $\bf X$ is given by

[TABLE]

where we have introduced two polarization displacements: $\bm{\rho}_{0}$ corresponding to the lowest order guiding-center reduction and $\bm{\rho}_{1}$ corresponding to the lowest order of the subsequent gyrocenter reduction.

The lowest order guiding-center displacement (in the guiding-center coordinates) is given by:

[TABLE]

where $\widehat{\bm{\rho}}$ is the unit vector in the plane perpendicular to the background magnetic field; the magnitude of magnetic field $B$ is evaluated at the reduced (guiding-center) position $\bf X$. The general gyrokinetic derivation comes up with a result that all the following guiding-center polarization displacements are at least of order $\mathcal{O}(\epsilon_{B})$ or higher (see, for example, Eq. (36) in [11] or Eqs. (63) and (66) in [12]), these corrections are neglected in the numerical implementation. The lowest order gyrocenter displacement is given by

[TABLE]

where $p_{z}$ is the gyrocenter scalar canonical momentum coordinate related to the parallel guiding-center momentum and will be defined in Eq. (6); $\phi_{1}$ represents the first order perturbative electrostatic potential and $A_{1\|}$ the first order electromagnetic parallel perturbative potential.

In this work we consider the gyrokinetic coordinate transformation together with the derivation of the corresponding reduced Hamiltonian models in two cases: First, in Sec. 4 we present the transformation containing all the FLR corrections, i.e., from the point of view of functional dependencies, containing corrections of all orders related to the guiding-center transformation where the particle is located at ${\bf x}={\bf X}+\bm{\rho}_{0}$. Later, in Sec. 5 we explicit this change of coordinates at the lowest FLR order, which corresponds from the physical point of view to the long-wavelength approximation, expressed in Fourier space with $\epsilon_{\perp}\equiv k_{\perp}\rho_{\mathrm{th}}\ll 1$.

3 Phase-space perturbative procedure

In gyrokinetic theory, the definition of the reduced phase-space coordinates is done within a common perturbative procedure together with the derivation of the reduced dynamics. At the first step, the guiding-center dynamical reduction starts from the local particle coordinates $\left(\mathbf{x},\mathbf{v}\right)$. To access those coordinates, one needs to define two bases of vectors: the static one and the dynamic one. The static basis is related to the background magnetic field line and the dynamic one rotates with the particle. As the static basis we choose the natural Frenet triad : the unit magnetic field vector $\widehat{\bf b}={\bf B}/B$, the normalized curvature vector

[TABLE]

and

[TABLE]

We notice that in the case of a uniform background magnetic field it is possible to choose a Cartesian frame as a static basis: $\widehat{\mathbf{b}}=\widehat{\mathbf{e}}_{z},\widehat{\mathbf{b}}_{1}=\widehat{\mathbf{e}}_{x}$ and $\widehat{\mathbf{b}}_{2}=\widehat{\mathbf{e}}_{y}$. Then, the dynamic basis $(\widehat{\bm{\rho}},\widehat{\mathbf{b}},\widehat{\perp})$ is defined from the static one as follows:

[TABLE]

where $\widehat{\bm{\rho}}$ is used for definition of the guiding-center displacement $\bm{\rho}_{0}$ in Eq. (2) such that the local particle velocity is decomposed in the following way:

[TABLE]

At the lowest order, the guiding-center transformation is defined as follows: the position of the particle is decomposed as ${\bf x}={\bf X}+\bm{\rho}_{0}({\bf X},\mu,\theta)$, with $\bf X$ the reduced particle (i.e., guiding-center) position and $\bm{\rho}_{0}$ the lowest order guiding-center polarization shift; the scalar momentum coordinate is the parallel kinetic momentum $p_{\|}=mv_{\|}$; $\mu$ is the lowest order magnetic moment given by Eq. (1) and $\theta$ is the fast rotation angle.

We consider the guiding-center phase space Lagrangian [2], [13] in the $\left({\bf X},p_{\|},\mu,\theta\right)$ coordinates as a starting point of the derivation:

[TABLE]

where the symplectic part contains the modified magnetic potential:

[TABLE]

The guiding-center Hamiltonian is given by:

[TABLE]

For the second step, i.e., the gyrocenter reduction, we consider the first order fluctuating time-dependent electromagnetic potential fields $\phi_{1}$ and $A_{1\|}$ both of order $\mathcal{O}(\epsilon_{\delta})$. Following the approximations currently performed on global code models, the perpendicular part of perturbed magnetic potential is omitted here. Moreover, as we have already stated in the introduction, the low-$\beta$ approximation is assumed in our derivation, which corresponds to the choice of considering the perpendicular component of the perturbed magnetic field only: $\mathbf{B}_{1}=\widehat{\mathbf{b}}\times\bm{\nabla}A_{1\|}$. These approximations are implemented in the electromagnetic global particle-in-cell code ORB5 as well as in the global version of the Eulerian code GENE.

To account for the time-dependence of the perturbed electromagnetic potentials $A_{1\|}$ and $\phi_{1}$, we extend the phase space from $6$ to $8$ dimensions. Therefore, formally the gyrocenter dynamical reduction is performed on the $8$-dimensional phase space where $(t,w)$ are canonically conjugate: $t$ corresponds to time and $w$ to energy. This extension of the phase space is used to make the dynamical system autonomous so as to perform the coordinate change in a more consistent way (see for example [14]). From the physical point of view, the relevant reduced dynamics is still performed on the $4$ dimensional part with coordinates $\left({\bf X},p_{z}\right)$, where

[TABLE]

is the guiding-center canonical momentum, and we notice that the perturbed part of magnetic potential $A_{1\|}$ enters in the definition of one of the phase space variables. Consequently, the perturbed guiding-center phase-space Lagrangian and Hamiltonian are [3, 10] :

[TABLE]

All the background quantities $\bf A$, $\widehat{\bf b}$ and $B$ are evaluated in the reduced guiding-center position $\bf X$, while the perturbative electromagnetic potentials $A_{1\|}$ and $\phi_{1}$ are evaluated in the reduced particle position ${\mathbf{X}+\bm{\rho}_{0}}$, i.e., the perturbed electromagnetic potentials depend on the gyroangle through the displacement $\bm{\rho}_{0}$. The first two terms of the right hand side of Eq. (7) represent the non-perturbed symplectic part and $H_{\mathrm{pert}}$, given by Eq. (8) is the perturbed Hamiltonian of the system.

This choice of phase space coordinates corresponds to the Hamiltonian [3], or $p_{z}$-, [15] representation of the perturbed guiding-center dynamics. It gives the possibility of keeping the symplectic part of the phase-space Lagrangian free from all electromagnetic perturbations and therefore, gyroangle-independent. The gyroangle-dependent terms are transferred into the expression for the perturbed Hamiltonian $H_{\mathrm{pert}}$. For this reason, the corresponding representation of the reduced dynamics is called ”Hamiltonian representation”. This maneuver has one significant advantage while constructing the phase-space dynamical reduction procedure for the phase-space Lagrangian: the (canonical) transformations have to be applied only on the Hamiltonian part of the phase space Lagrangian, since the corresponding symplectic part is already free of any $\theta$ dependence.

The Hamiltonian representation is the most common choice for the models implemented in particle-in-cell simulations since it avoids the appearance of the inductive electric field (i.e., the explicit time-derivative of the perturbative magnetic potential $A_{1\|}$) in the particle characteristics given in Sec. 6. For example, a control-variate scheme implemented in PIC code ORB5 [16] is using the Hamiltonian representation.

In the next section we show how to eliminate the gyroangle-dependence from the perturbed guiding-center phase-space Lagrangian induced by the second-order $\sim\mathcal{O}(\epsilon_{\delta}^{2})$ perturbed electromagnetic potentials. To that purpose we build a Lie-transform near-identity transformation and move from the guiding-center to the gyrocenter variables by applying it systematically to the phase space variables and at the same time to the perturbed Hamiltonian $H_{\mathrm{pert}}$, given by Eq. (8). Our goal is to clarify the connection between the displacements $\bm{\rho}_{0}$ and $\bm{\rho}_{1}$, and to eliminate the gyroangle-dependence from the reduced dynamics.

4 Full FLR Hamiltonian model

In this section we build a near-identity phase-space transformation aiming to eliminate the gyroangle-dependence from the perturbed phase-space Lagrangian (7). Since within the $p_{z}$- representation the symplectic part of the phase-space Lagrangian does not contain any $\theta$-dependence, the gyrocenter phase-space transformation will only modify its corresponding Hamiltonian part. For the near-identity phase-space transformations we use Lie transforms which have several advantages, among them, they are canonical transformations, meaning that they do not affect the expression of the Poisson backet, only the expression of the Hamiltonian.

To define this technique, we need first of all, a Poisson bracket for the guiding-center dynamics, which can be derived from the symplectic part of the perturbed Lagrangian given by Eq. (7), which coincides with the symplectic part of the unperturbed one defined by Eq. (5) (see [13] for more details):

[TABLE]

where $\mathbf{B}^{*}=\mathbf{B}+\frac{e}{c}p_{z}\bm{\nabla}\times\widehat{\mathbf{b}}$ and $B_{\|}^{*}=\mathbf{B}^{*}\cdot\widehat{\mathbf{b}}$. For physical reasons, we know that there is a time scale separation among the motion described by the different terms of the Poisson bracket. In order to exploit this separation we order the first three terms of the bracket following the formal ordering introduced by Newcomb [17], $\epsilon=e^{-1}$: the first term $\propto e$ is related to the fast rotation around the magnetic field line (fastest scale of motion), the second term $\propto e^{0}$ represents the parallel motion, and the third term $\propto e^{-1}$ is related to the slow drifts in the perpendicular direction. The last term in the bracket corresponds to the extension of the phase space to $8$ dimensions (autonomization of the system). In the following calculations we fix the small parameter $\epsilon=\epsilon_{\delta}$, so that the Poisson bracket (9) is separated in three brackets according to the scale of motion:

[TABLE]

where

[TABLE]

Following the $p_{z}$-representation of the Lagrangian defined in Eq. (7), we rewrite the perturbed Hamiltonian $H_{\mathrm{pert}}$, given by Eq. (8) as follows:

[TABLE]

where the unperturbed guiding-center Hamiltonian now writes as:

[TABLE]

and

[TABLE]

is the linear perturbed gyrocenter potential. Note that the guiding-center displacement $\bm{\rho}_{0}$ given by Eq. (2) is depending on the phase-space coordinates $(\mathbf{X},\mu,\theta)$. To make formulas more compact we omit writing the functional dependencies of the displacement $\bm{\rho}_{0}$ explicitly when this is unambiguous. We notice that the spatial dependence of $\bm{\rho}_{0}$ can be omitted since $\left|\bm{\nabla}\bm{\rho}_{0}\right|\sim\mathcal{O}(\epsilon_{B})$.

4.1 Gyrocenter phase-space coordinate transformation

A Lie transform generated by a scalar differentiable function $S$ is defined by its action on observables $G$ as:

[TABLE]

where $\varepsilon$ is a small parameter of a given problem and $\{\cdot,\cdot\}$ is a Poisson bracket. To this change of observables corresponds an invertible change of coordinates by the scalar invariance: $\bar{G}(\bar{\bf z})=G({\bf z})$, where the new coordinates are given by

[TABLE]

where we have used the fundamental property of a Lie transform $\varphi({\rm e}^{-\varepsilon\pounds_{S}}G)={\rm e}^{-\varepsilon\pounds_{S}}\varphi(G)$ for any scalar function $\varphi$. Another fundamental property is that Lie transforms are canonical changes of coordinates, in the sense that the expression of the Poisson bracket in the new variables is exactly the same as the one in the old variables. This comes from the property $\{{\rm e}^{-\varepsilon\pounds_{S}}F,{\rm e}^{-\varepsilon\pounds_{S}}F\}={\rm e}^{-\varepsilon\pounds_{S}}\{F,G\}$.

The full gyrokinetic coordinate transformation represents an infinite series of near-identity phase-space transformations, aiming to remove the gyroangle dependence from the reduced dynamics at all orders. The key element of the reduction procedure is the identification of a generating function $S$, which plays a double role: first, defining new phase space coordinates and second, the reduced dynamics. The generating function $S$ is constructed from a perturbative series. Keeping the chosen ordering in mind, we are using the Lie transform given by Eq. (15) with the generating function at the first order $S$, and the Poisson bracket given by Eq. (9), we define the gyrocenter phase-space coordinates as:

[TABLE]

Therefore, in principle, there are two sets of coordinates, the guiding-center ones $(\mathbf{X},p_{z},\mu,\theta)$, and the gyrocenter ones $(\bar{\mathbf{X}},\bar{p_{z}},\bar{\mu},\bar{\theta})$. Since, after performing the Lie transforms, it is clear that the coordinates are the gyrocenter ones, we will omit the bars over the new variables where it is not ambiguous, and denote the gyrocenter phase space coordinates as $(\mathbf{X},p_{z},\mu,\theta)$ in what follows.

The gyrocenter displacement is defined from the Lie transform of the position of the particle expressed in the guiding-center coordinates:

[TABLE]

where

[TABLE]

at the first order. An explicit calculation of the polarization displacement $\bm{\rho}_{1}$ requires the knowledge of the generating function $S$, which is also necessary to get the expression for the reduced Hamiltonian and therefore to derive the reduced dynamics. This emphasizes the essential link between the definition of the new phase-space coordinates and derivation of the reduced gyroangle-independent dynamics.

We derive the expression for the reduced Hamiltonian together with the expression for the lowest order generating function $S$ in the following section.

4.2 Full FLR gyrocenter dynamics

In order to accommodate various orderings and small parameters present in the problem, and in particular to eliminate the singularity in $\epsilon_{\delta}$ appearing in the first term of the Poisson bracket (10), we expand the generating function $S=\epsilon_{\delta}S_{1}+\epsilon_{\delta}^{2}S_{2}+\mathcal{O}(\epsilon_{\delta}^{3})$. Taking Eqns. (10,11,12,13) into account, the corresponding Lie transform up to order $\epsilon_{\delta}^{3}$ is given by:

[TABLE]

The purpose of the Lie transform is to eliminate the gyroangle dependence from the reduced dynamics up to the second-order in $\epsilon_{\delta}$. As already mentioned, since the Poisson bracket (9) is gyroangle-independent, we have to eliminate the gyroangle dependence only from the perturbed Hamiltonian (14).

Each function on the gyrocenter phase space evaluated at the position ${\mathbf{X}}+\bm{\rho}_{0}$, and therefore gyroangle dependent, can be decomposed in its gyoaveraged and fluctuating parts:

[TABLE]

where

[TABLE]

First, we apply the Lie transform (22) to the perturbed Hamiltonian (14), and identify a fluctuating and a gyroaveraged part for each term:

[TABLE]

Second, we identify order by order fluctuating contributions, which should be removed by a well-chosen generating function:

[TABLE]

The expression for the lowest order generating function $S_{1}$ is obtained from Eq. (23), which represents a condition that the gyroangle-dependent part of linear electromagnetic perturbation $\widetilde{\psi}_{1}$ is removed from the lowest order gyrocenter Hamiltonian. Since $H_{0}$ does not depend on $\theta$, the equation determining $S_{1}$ reduces to:

[TABLE]

Therefore, the lowest order generating function writes:

[TABLE]

We recover the conventional expression for the generating function, used for the gyrokinetic calculations in [18], [3],[8].

We are now considering Eq. (24). Since $S_{2}$ removes the fluctuating parts from the second-order terms in the Hamiltonian, we only need to evaluate the corresponding gyroaveraged contributions to the reduced Hamiltonian. Taking into account Eq. (30) which indicates that $S_{1}$ is a purely fluctuating function, the second term in Eq. (24), namely ${\left\{S_{1},H_{0}\right\}}_{0}$, is purely fluctuating and will not affect the expression of the second order reduced Hamiltonian. For similar reasons, the third term possesses the following property: $\left\langle\{S_{1},\psi_{1}\}_{-1}\right\rangle=\left\langle\{S_{1},\widetilde{\psi}_{1}\}_{-1}\right\rangle$, i.e., we obtain a partial cancellation of the second-order terms:

[TABLE]

The next order correction to the generating function can be obtained from Eq. (24) using the expression for the generating function $S_{1}$ , which we have now determined in Eq. (31):

[TABLE]

Here we do not explicitly resolve Eq. (32), we just use the fact that $S_{2}$ removes all the fluctuating parts at the order $\mathcal{O}(\epsilon_{\delta}^{2})$. We notice that the perturbative procedure can be expanded to higher orders.

Finally, we get the expression for the second-order Hamiltonian containing the FLR corrections (i.e., with respect to the polarization displacement $\bm{\rho}_{0}$) at all orders:

[TABLE]

where $S_{1}$ is given by Eq. (31). Rewriting the expression for the lowest order Poisson bracket:

[TABLE]

we have

[TABLE]

and taking into account that $\partial_{\theta}S_{1}=\frac{e}{\Omega}\widetilde{\psi}_{1}(\mathbf{X}+\bm{\rho}_{0},p_{z})$, we get:

[TABLE]

Therefore, the reduced Hamiltonian becomes:

[TABLE]

For further convenience, we decompose the expression for the Hamiltonian as $H_{\rm gy}=\overline{H}_{0}+\epsilon_{\delta}\overline{H}_{1}+\epsilon_{\delta}^{2}\overline{H}_{2}$ with:

[TABLE]

where $\overline{H}_{0}$ is the non-perturbed guiding-center Hamiltonian, $\overline{H}_{1}$ is the first order gyrocenter Hamiltonian and $\overline{H}_{2}$ is the second-order gyrocenter Hamiltonian, which contains the FLR corrections at all orders. This result corresponds to the expression obtained in a slab magnetic geometry [19].

We recall that since the expression for the generating function $S$ has been defined within the dynamical reduction procedure for the Hamiltonian, we have obtained the explicit expression for the second-order Hamiltonian without writing explicitly the expression for the new (i.e., gyrocenter) phase space coordinates. We are now completing the dynamical reduction procedure by providing an explicit expression for the gyrocenter phase space coordinates.

4.3 Explicit gyrocenter phase space coordinate transformation

Using the expression for the generating function $S_{1}$ given by Eq. (31) and the guiding-center Poisson bracket given by Eq. (9), we explicitly evaluate the expression for the corresponding lowest order gyrocenter coordinate transformation given by Eqs. (16-19), i.e., we express the leading order modifications of the gyrocenter coordinates $(\bar{\bf X},\bar{p_{z}},\bar{\mu},\bar{\theta})$ as functions of the guiding-center coordinates $({\bf X},p_{z},\mu,\theta)$:

[TABLE]

We notice that the modifications for $\bar{\bf X}$ and $\bar{p_{z}}$ are of order $\epsilon_{\delta}^{2}$, whereas the modifications for $\bar{\mu}$ and $\bar{\theta}$ are of order $\epsilon_{\delta}$. Because of the modification of $\bar{\mu}$ and $\bar{\theta}$ at the order $\epsilon_{\delta}$, the gyrocenter displacement is of order $\epsilon_{\delta}$.

Concerning the coordinate change, we recall that there are two important sets of coordinates in the two-step gyrokinetic reduction, the guiding-center coordinates $({\bf X},p_{z},\mu,\theta)$ and the gyrocenter ones $(\bar{\bf X},\bar{p_{z}},\bar{\mu},\bar{\theta})$ which we have also denoted $({\bf X},p_{z},\mu,\theta)$ after Lie transforming the Hamiltonian. The position of the particle is ${\bf x}={\bf X}+{\bm{\rho}}_{0}$ in the guiding-center coordinates and ${\bf x}={\bf X}+{\bm{\rho}}_{0}+{\bm{\rho}}_{1}$ in the gyrocenter coordinates (which, for clarity, should be denoted ${\bf x}=\bar{\bf X}+\bar{\bm{\rho}}_{0}+\bar{\bm{\rho}}_{1}$). Therefore, in the gyrocenter coordinates ${\bm{\rho}}_{0}+{\bm{\rho}}_{1}$ is the difference between the position of the particle and the position of the gyrocenter, whereas $\bm{\rho}_{0}$ is the difference between the position of the particle and the guiding center in the guiding-center coordinates. It is tempting to conclude that the gyrocenter displacement ${\bm{\rho}}_{1}$ is the distance between the guiding center and the gyrocenter. However this is incorrect. We have seen above that the gyrocenter displacement is of order $\epsilon_{\delta}$ whereas the distance between the guiding center and the gyrocenter is of order $\epsilon_{\delta}^{2}$. The way to solve the apparent contradiction is to carefully evaluate the system in which the positions are expressed. The expressions of ${\bm{\rho}}_{0}$ in the guiding center and in the gyrocenter coordinates are different, even at the leading order: At the leading order, the difference between these two expressions is given by

[TABLE]

where we have to distinguish both sets of coordinates since they are used in the same equation. The shift in ${\bm{\rho}}_{0}$ by the Lie transform generates exactly the gyrocenter displacement ${\bm{\rho}}_{1}$ at the leading order, whereas the reduced positions remain approximately unchanged (at least at the order $\epsilon_{\delta}$).

To account for all FLR corrections, i.e., corrections related to the displacement $\bm{\rho}_{0}$, in the expressions for the reduced Hamiltonian given by Eqs. (38-40) we have evaluated the fields $\phi_{1}$ and $A_{1\|}$ at the position $\bar{\mathbf{X}}+\bm{\rho}_{0}(\bar{\mathbf{X}},\bar{\mu},\bar{\theta})$ instead of the particle position $\mathbf{x}=\bar{\mathbf{X}}+\bm{\rho}_{0}(\bar{\mathbf{X}},\bar{\mu},\bar{\theta})+\bm{\rho}_{1}(\bar{\mathbf{X}},\bar{\mu},\bar{\theta})$. This generates polarization and magnetization corrections related to the displacement $\bm{\rho}_{1}$ into the second term of the second order Hamiltonian $\bar{H}_{2}$ given by Eq. (40). To make it more explicit, we consider the Lie transform of the gyrocenter potential $\psi_{1}$ evaluated at the position $(\mathbf{X}+\bm{\rho}_{0}({\mathbf{X}},\mu,\theta))$:

[TABLE]

which is obtained from Eq. (20). The expansion of the left-hand side leads to

[TABLE]

where $\bar{\bm{\rho}_{0}}=\bm{\rho}_{0}(\bar{\mathbf{X}},\bar{\mu},\bar{\theta})$. The Taylor expansion of the right-hand side leads to , since :

[TABLE]

since $\bm{\rho}_{1}\propto\mathcal{O}(\epsilon_{\delta})$. Therefore we have $\left\{S_{1},\psi_{1}\right\}_{-1}=-\bm{\rho}_{1}\cdot\bm{\nabla}\psi_{1}$ evaluated at $\bar{\mathbf{X}}+\bar{\bm{\rho}_{0}}$. In other words, according to Eqs. (36) and (40), the polarization effects associated with the displacement $\bm{\rho}_{1}$ are contained in the second order gyrocenter Hamiltonian $\bar{H}_{2}$.

5 Hamiltonian model in long-wavelength approximation

In the previous section, we have derived the gyrocenter Hamiltonian model (37) containing up to $\epsilon_{\delta}^{2}$ corrections. In this section we consider an additional ordering $\epsilon_{\perp}=k_{\perp}\rho_{\mathrm{th}}\ll 1$ and derive the Hamiltonian model in the long-wavelength approximation, containing up to the second-order FLR corrections, i.e., terms of order $\mathcal{O}(\epsilon_{\perp}^{2})$. Such a derivation has a direct interest with respect to the numerical implementation of gyrokinetic Vlasov-Maxwell models. Indeed, from the point of view of the variational derivation, the second-order terms in Eq. (37) will enter the gyrokinetic Maxwell equations as the linear (i.e., of order $\mathcal{O}(\epsilon_{\delta})$) polarization corrections (see for example a detailed derivation in [10]).

The gyrokinetic Vlasov-Poisson model with linearized polarization term derived in the long-wavelength limit with $\epsilon_{\perp}=k_{\perp}\rho_{\mathrm{th}}\ll 1$ has been the first implemented model in a gyrokinetic PIC code [20]. It is still widely used for the investigation of MHD modes, and it is also suitable for studies of turbulence generated by the interaction of modes with low toroidal numbers. As it has been shown in the latest linear electromagnetic benchmark [21], the long-wavelength approximation implemented into the ORB5 code allows one to treat modes with $\epsilon_{\perp}=k_{\perp}\rho_{\mathrm{th}}<0.6$.

From the numerical point of view, the long-wavelength approximate model, containing up to $\mathcal{O}(\epsilon_{\perp}^{2})$ corrections is of special interest because of its consistency with respect to the $n$-point gyroaverage operator approximation, implemented in ORB5. The $n$-point operator is an extension of the original four-point gyroaverage proposed in [20].

The expression for the gyroaverage operator, applied to a scalar field $\psi(\mathbf{X}+\bm{\rho}_{0},p_{z},\mu,\theta)$ can be expressed in Fourier space in the following form:

[TABLE]

where $\widehat{\psi}$ is the Fourier-transformed scalar field $\psi$. The last expression shows that the action of the operator $\mathcal{J}_{0}^{\mathrm{gc}}$ in Fourier space is translated into the multiplication of the Fourier coefficients by the Bessel function $\mathrm{J}_{0}(k_{\perp}\rho_{i})$. In order to be implementable, the gyroaverage procedure is approximated by an average over a finite number of points on the gyroring. In the first PIC code [20] the four-point gyroaverage was used, this procedure is equivalent to considering the finite difference approximation of the Taylor expansion $\mathrm{J}_{0}(k_{\perp}\rho_{i})\approx 1-(k_{\perp}\rho_{i})^{2}/4$, which corresponds in real space to:

[TABLE]

A detailed discussion is available in [22].

Recently, a full FLR solver for the Poisson equation has been implemented in ORB5 [23], [24]. Preliminary results show that the new algorithm is twice slower than the long-wavelength solver. This is mostly due to the need for additional integration points for the gyroaverage algorithm as it was already shown in [25].

We present two different ways to obtain the simplified Hamiltonian model, containing up to the $\epsilon_{\perp}^{2}$ terms. First of all, we perform the lowest order FLR series truncation directly on the perturbative electromagnetic potential $\widetilde{\psi}_{1}$ and we keep up to the second-order FLR terms inside the magnetic potential $A_{1||}$, which enter into the expression of the second-order $\mathcal{O}(\epsilon_{\delta}^{2})$ electromagnetic potential given by Eq. (36). Next, we are following the main steps of the gyrocenter coordinate transformation by introducing the FLR truncation at each step: starting with the expression for the generating function $S_{1}$, getting the corresponding gyrocenter change of coordinates and finally the expression for the simplified Hamiltonian.

In this section, we follow the direct $\epsilon_{\perp}$- truncation of Eq. $(\ref{eq:FullFLR})$ and in Sec. 5.2 we detail the main steps of the second derivation, in which we show that the first order $\epsilon_{\perp}$-truncation of the generating function $S_{1}$ is sufficient for recovering the second-order FLR, i.e., $\epsilon_{\perp}^{2}$ corrections in the Hamiltonian.

5.1 Direct full FLR model truncation

Here we evaluate the second-order, i.e., $\mathcal{O}(\epsilon_{\perp}^{2})$, FLR contributions to the second-order Hamiltonian (40). We start by decomposing the first-order fluctuating electromagnetic potential in the FLR series:

[TABLE]

To get the $\epsilon_{\perp}^{2}$- contributions in the Hamiltonian, we only keep the first term, and we calculate:

[TABLE]

where we have used the definition (2), the fact that $\frac{1}{2}\frac{\partial\rho_{0}^{2}}{\partial\mu}=\frac{c^{2}m}{e^{2}B}$ and the dyadic tensors property $\left\langle\widehat{\bm{\rho}}\widehat{\bm{\rho}}\right\rangle=\frac{1}{2}\left(\widehat{\mathbf{b}}_{1}\widehat{\mathbf{b}}_{1}+\widehat{\mathbf{b}}_{2}\widehat{\mathbf{b}}_{2}\right)=\frac{1}{2}\mathbf{1}_{\perp}$.

The magnetic term of the second-order is truncated up to the second FLR correction:

[TABLE]

We note that the second term is missing in the ORB5 model [10], which corresponds to the slab geometry result obtained in [19]. Finally, truncated up to the second-order, the FLR Hamiltonian writes:

[TABLE]

which for further convenience is separated into three parts: the unperturbed part $\bar{H}_{0}$ and the first order contribution $\bar{H}_{1}$ are defined by Eqs. (38) and (39), while the second-order term is defined as:

[TABLE]

5.2 Gyrocenter coordinate transformation in the long-wavelength approximation

In this section we derive the truncated Hamiltonian model (5.1) by performing the gyrocenter coordinate transformation (16-19) at the lowest FLR order. To that purpose, we only take into account the lowest order FLR correction to the generating function $S=\epsilon_{\delta}S_{1}+\mathcal{O}(\epsilon_{\delta}^{2})$. Since the lowest order FLR correction to the fluctuating electromagnetic potential is $\widetilde{\psi}_{1}\left(\mathbf{X}+\bm{\rho}_{0},p_{z}\right)=\rho_{0}\widehat{\bm{\rho}}\cdot\bm{\nabla}\psi_{1}(\mathbf{X},p_{z})$, from Eq. (31) using the property of the rotating basis vectors $\widehat{\perp}=-\int d\theta\ \widehat{\bm{\rho}}$ , we get:

[TABLE]

We calculate the gyrocenter displacement $\bm{\rho}_{1}$ using the definition (21), taking into account the lowest order Poisson bracket (11):

[TABLE]

From the definition of rotating basis vectors (4), and $\partial_{\mu}\rho_{0}^{2}=\frac{2mc^{2}}{e^{2}B}$, we have:

[TABLE]

Using the definition of the dyadic tensor $\bf{1}_{\perp}\equiv\widehat{\bm{\rho}}\widehat{\bm{\rho}}+\widehat{\perp}\widehat{\perp}$, the first order gyrocenter displacement in the long wavelength approximation is:

[TABLE]

At the leading order, only the gyrocenter phase-space coordinates $(\mu,\theta)$ are transformed according to the Eqs. (18,19) using:

[TABLE]

where we have used some properties of the rotating basis vectors: $\partial_{\theta}\widehat{\perp}=-\widehat{\bm{\rho}}$ and $\partial_{\theta}\widehat{\bm{\rho}}=\widehat{\perp}$.

Therefore, truncating the second order gyrocenter Hamiltonian (obtained by Lie transform generated by $S_{1}$ containing all the FLR corrections) up to the second order in FLR corrections is equivalent to the second order gyrocenter Hamiltonian obtained from the Lie transform generated by the first FLR correction of the generating function $S_{1}$.

6 Gyrokinetic Vlasov equation

In this section we derive an expression for the gyrokinetic Vlasov equation corresponding to the Hamiltonian characteristics according to the models obtained in Secs. 4 and 5.

The particle characteristics are calculated as follows:

[TABLE]

where $H_{\rm gy}$ is the Hamiltonian chosen according to the approximation, full FLR model (37) or truncated up to second FLR order model (5.1). The Poisson bracket is given by Eq. (9). Then the corresponding Vlasov equation suitable for the $\delta f$ simulations, i.e., including the decomposition $F=F_{0}+\epsilon_{\delta}F_{1}$ into the background $F_{0}$ and the fluctuating part $F_{1}$ of the particle distribution function, is reconstructed from the characteristics given by Eqs. (43) and (44) by saying that the distribution function is constant along the characteristics:

[TABLE]

and then

[TABLE]

where we have taken into account the functional dependencies of the gyrocenter distribution function $F=F({\mathbf{X}},p_{z},\mu)$ and the fact that $\mu$ is constant on the reduced gyrocenter phase space. We assume that $\left\{F_{0},\bar{H}_{0}\right\}_{\mathrm{gc}}=0$, i.e., the background particle distribution is time-independent; therefore there are two linear and two non-linear terms). Following the Hamiltonian decomposition given by Eqs. (38), (39), the Vlasov equation (45) can further be expanded as:

[TABLE]

where $\bar{H}_{2}$ is the second-order Hamiltonian chosen according to the approximation, full FLR model (40) or FLR truncated model (5.1). The last term refers to the evolution of the background distribution under the second-order Hamiltonian. This term is systematically neglected in most of global codes, but is necessary for full-$f$ nonlinear simulations.

7 Conclusions

We have performed a derivation of the second-order (in $\epsilon_{\delta}$) gyrocenter Hamiltonian models in the case with full FLR corrections and then considering the long-wavelength approximation. The long-wavelength model has been obtained in two different ways: by using a direct truncation of the full FLR model and by constructing the dynamical reduction procedure with the gyrocenter generating function, containing only the lowest order FLR contributions. Whether the FLR truncation is performed a priori on the generating function, or a posteriori on the reduced Hamiltonian does not change the final expressions at the leading orders. In the course of the reduction, we have clarified the origin of the gyrocenter displacement, responsible for polarization and magnetization terms in the reduced gyrokinetic Vlasov-Maxwell equations, using a proper ordering of the generating function of the Lie transform.

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training Programme No. 2014- 2018 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.