Guiding-center kinetic model based on the assumption of homogeneous distribution over gyrophase

TL;DR

This paper introduces a simplified guiding-center kinetic model for magnetized plasmas, assuming homogeneous gyrophase distribution, which simplifies the dynamics by eliminating gyrophase dependence and incorporates Coulomb interactions.

Contribution

It develops a new guiding-center kinetic model based on the homogeneous gyrophase distribution assumption, including Coulomb interactions, and derives a gyrophase-independent fundamental Lagrangian form.

Findings

The model simplifies plasma dynamics by removing gyrophase dependence.

It incorporates Coulomb interactions directly into the guiding-center framework.

The approach enables more efficient plasma simulations with Coulomb effects.

Abstract

The purpose of this paper is to develop a simplified model as the modeling of the magnetized plasmas. The starting point is an assumption that the distribution of the ensemble of charged particles in the same species is homogeneous over gyrophase. The particles in this ensemble are located at the same guiding-cneter position $(\mathbf{X},\mu,U)$. Then, a fundamental Lagrangian differential 1-form is developed. It contains all particles in the magnetized plasma system as well as the Coulomb pair force between particles instead of field-particle interaction used in conventional gyrokinetic models. By transforming the Lagrangian 1-form to the new one on guiding-center coordinate with the ensemble summation over gyrophase, the new fundamental 1-form is naturally independent of gyrophase of each particle based on the homogeneous distribution over gyrophase, and it determines the dynamics of all particles on the new coordinates. By using a coarse-grained scheme, this new 1-form can be modeled by the guiding-center kinetic model.

Equations

$$\varphi \left(\mathbf{x},t\right)=\frac{1}{4\pi {ϵ}_0}\sum _j^{\prime }\left(\frac{q}{|\mathbf{x}-{\mathbf{x}}_{ij}(t)|}-\frac{e}{|\mathbf{x}-{\mathbf{x}}_{ej}(t)|}\right),$$

$$\gamma ={\gamma }_i+{\gamma }_e,$$

$${\gamma_{i}}=\sum\limits_{j}{\left[{\begin{array}[]{*{20}{l}}{\left({q{\bf{A}}\left({{{\bf{x}}_{ij}}}\right)+{m_{i}}{{\bf{v}}_{ij}}}\right)\cdot d{{\bf{x}}_{ij}}}\\ {-\left({\frac{{{m_{i}}{\bf{v}}_{ij}^{2}}}{2}+\frac{{q\phi\left({{{\bf{x}}_{ij}},t}\right)}}{2}}\right)dt}\end{array}}\right]},$$

$${\gamma_{e}}=\sum\limits_{j}{\left[{\begin{array}[]{*{20}{l}}{\left({-e{\bf{A}}\left({{{\bf{x}}_{ej}}}\right)+{m_{e}}{{\bf{v}}_{ej}}}\right)\cdot d{{\bf{x}}_{ej}}}\\ {-\left({\frac{{{m_{e}}{\bf{v}}_{ej}^{2}}}{2}-\frac{{e\phi\left({{{\bf{x}}_{ej}},t}\right)}}{2}}\right)dt}\end{array}}\right]}.$$

$$\gamma=\sum\limits_{o\in\{i,e\}}{\sum\limits_{j}{\left[{\begin{array}[]{*{20}{l}}{\left({{q_{o}}{\bf{A}}\left({{{\bf{x}}_{oj}}}\right)+{m_{o}}{{\bf{v}}_{oj}}}\right)\cdot d{{\bf{x}}_{oj}}}\\ {-\left({\frac{{{m_{o}}{\bf{v}}_{oj}^{2}}}{2}+\frac{{{q_{o}}\phi\left({{{\bf{x}}_{oj}},t}\right)}}{2}}\right)dt}\end{array}}\right]}}.$$

$${\nabla }^2\varphi \left(\mathbf{x},t\right)=-\frac{1}{{ϵ}_0}\sum _j\left[q\delta \left(\mathbf{x}-{\mathbf{x}}_{ij}(t)\right)-e\delta \left(\mathbf{x}-{\mathbf{x}}_{ej}(t)\right)\right].$$

$$\begin{array}[]{l}\int_{V}{{\nabla^{2}}\phi\left({{\bf{x}},t}\right){d^{3}}{\bf{x}}}\\ =-\frac{1}{{{\epsilon_{0}}}}\int_{V}{\sum\limits_{j^{\prime}}{\left[\begin{array}[]{l}q\delta\left({{\bf{x}}-{{\bf{x}}_{ij}}(t)}\right)\\ -e\delta\left({{\bf{x}}-{{\bf{x}}_{ej}}(t)}\right)\end{array}\right]}{d^{3}}{\bf{x}}},\end{array}$$

$$\sum _k{\nabla }^2{\varphi }_k\left({\mathbf{x}}_k,t\right)\Delta V_k$$

$$=-\frac{1}{{ϵ}_0}\sum _k\frac{{\beta }_{ik}(t)-{\beta }_{ek}(t)}{\Delta V_k}\Delta V_k.$$

$$\frac{{\beta }_{ok}(t)}{\Delta V_k}\Delta V_k=q_o{n}_{ok}({\mathbf{x}}_k,t)\Delta V_k.$$

$${\nabla }^2{\varphi }_k({\mathbf{x}}_k,t)=\frac{1}{{ϵ}_0}\left(e{n}_{ek}({\mathbf{x}}_k,t)-q{n}_{ik}({\mathbf{x}}_k,t)\right).$$

$$\varphi \left({\mathbf{x}}_j,t\right)={\varphi }_k({\mathbf{x}}_k,t)+O(l_c).$$

$$M_{ok}\left(\mathbf{z}\right)=\sum _j\delta \left(\mathbf{x}-{\mathbf{x}}_{oj}(t)\right)\delta \left(\mathbf{v}-{\mathbf{v}}_{oj}\left(t\right)\right)$$

$$n_{ok}({\mathbf{x}}_k,t)=\frac{1}{\Delta V_k}{\int }_{\Delta V_k}{M}_{ok}\left(\mathbf{z}\right)d^3\mathbf{x}{d}^3\mathbf{v}.$$

$${\gamma }_{oj}$$

$$\left(\frac{\partial }{\partial t}+\frac{d\mathbf{x}}{dt}\cdot \nabla +\frac{d\mathbf{v}}{dt}\cdot \frac{\partial }{\partial \mathbf{v}}\right)M_{ok}\left(\mathbf{z}\right)=0.$$

$$\left(\frac{\partial }{\partial t}+\frac{d\mathbf{x}}{dt}\cdot \nabla +\frac{d\mathbf{v}}{dt}\cdot \frac{\partial }{\partial \mathbf{v}}\right)f_{ok}\left(\mathbf{z}\right)=0.$$

$$n_e(\mathbf{x},t)=n_{e0}(\mathbf{x})+\frac{e\varphi \left(\mathbf{x},t\right)}{T_e}{n}_{e0}(\mathbf{x}),$$

$${\nabla^{2}}{\phi_{k}}(\mathbf{x}_{k},t)=-\frac{1}{{{\epsilon_{0}}}}\left({\begin{array}[]{*{20}{l}}{q\int{{}_{\Delta{V_{k}}}{{\tilde{f}}_{ik}}\left({\bf{z}}\right)B({\bf{x}})d{\mu_{1i}}d{U_{1i}}}}\\ {-\frac{{{e^{2}}\phi_{k}\left({{\bf{x}}_{k},t}\right)}}{{{T_{e}}}}{n_{e0}}({\bf{x}},t)}\end{array}}\right)$$

$$\frac{d{Z}_{ojf}^m}{d{\epsilon }_o}\left({\mathbf{z}}_{oj},{\epsilon }_o\right)=g_{oj}^m\left({\mathbf{Z}}_{oj}\right),$$

$$\frac{d{\mathbf{z}}_{oj}}{d{\epsilon }_o}=0,$$

$$\frac{\partial z_{oj}^i}{\partial {\epsilon }_o}+\frac{d{Z}_{oj}^k}{d{\epsilon }_o}\frac{d{z}_{oj}^i}{d{Z}_{oj}^k}=0,$$

$${\mathbf{z}}_{oj}\left({\mathbf{Z}}_{oj},{\epsilon }_o\right)=\exp \left(-{\epsilon }_o{g}_{oj}^m\left({\mathbf{Z}}_{oj}\right){\partial }_{Z_{oj}^m}\right){\mathbf{Z}}_{oj}.$$

$$\begin{array}[]{l}{\bar{\bf{z}}_{o}}\equiv\left({{{\bf{z}}_{o1}},{{\bf{z}}_{o2}},{{\bf{z}}_{o3}},\cdots}\right),\\ {\bar{\bf{Z}}_{o}}\equiv\left({{{\bf{Z}}_{o1}},{{\bf{Z}}_{o2}},{{\bf{Z}}_{o3}},\cdots}\right),\\ {\bar{\bf{g}}_{o}}\equiv\left({{{\bf{g}}_{o1}},{{\bf{g}}_{o2}},{{\bf{g}}_{o3}},\cdots}\right).\\ \end{array}$$

$${\bar{\mathbf{z}}}_{oj}\left({\bar{\mathbf{Z}}}_o,{\epsilon }_o\right)=\exp \left(-{\epsilon }_o{\overline{g}}_o^m{\partial }_{{\bar{Z}}_o^m}\right){\bar{\mathbf{Z}}}_o.$$

$$\frac{\partial {\bar{\Gamma }}_m}{\partial {\epsilon }_i}\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{z}}}_e,{\epsilon }_i,{\epsilon }_e\right)$$

$$=-{\bar{g}}_i^n\left({\overline{\mathbf{Z}}}_i\right)\left[\frac{\partial {\bar{\Gamma }}_m}{\partial {\bar{Z}}^n}\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{z}}}_e,{\epsilon }_i,{\epsilon }_e\right)-\frac{\partial {\bar{\Gamma }}_n}{\partial {\bar{Z}}^m}\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{z}}}_e,{\epsilon }_i,{\epsilon }_e\right)\right]$$

$$-\frac{\partial }{\partial {\bar{Z}}^m}\left[{\bar{g}}_i^n\left({\overline{\mathbf{Z}}}_i\right){\bar{\Gamma }}_n\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{z}}}_e,{\epsilon }_i,{\epsilon }_e\right)\right],$$

$$\bar{\Gamma }\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{z}}}_e,{\epsilon }_i,{\epsilon }_e\right)=\exp \left(-{\epsilon }_o{L}_{{\overline{\mathbf{g}}}_i}\right)\gamma \left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{z}}}_e,{\epsilon }_i,{\epsilon }_e\right)+dS.$$

$$\frac{\partial {\Gamma }_m}{\partial {\epsilon }_e}\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{Z}}}_e,{\epsilon }_i,{\epsilon }_e\right)$$

$$-\frac{\partial }{\partial {\bar{Z}}_e^m}\left[{\bar{g}}_e^n\left({\overline{\mathbf{Z}}}_i\right){\Gamma }_n\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{Z}}}_e,{\epsilon }_i,{\epsilon }_e\right)\right],$$

$$\Gamma \left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{Z}}}_e,{\epsilon }_i,{\epsilon }_e\right)=\exp \left(-{\epsilon }_e{L}_{{\overline{\mathbf{g}}}_e}\right)\bar{\Gamma }\left({\overline{\mathbf{Z}}}_i,{\overline{\mathbf{Z}}}_e,{\epsilon }_i,{\epsilon }_e\right)+dS.$$

Full text

1 Introduction

The composited coordinate transform used in conventional gyrokinetic models(CGM)[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] comprises two independent but consecutive coordinate transform. The first one is guiding center transform[2, 3, 4]; the other one is gyrocenter transform[2, 3, 4, 5]. The purpose of this composited transform is to reduce the gyrophase from the orbit equation of the charged particles, so that the equation of the evolution of the distribution function only solves a five-dimensional distribution[12, 13, 14, 15, 16, 17, 18, 19]. The Poisson equation is solved on particle coordinates, so that the distribution needs to be transformed back to the one on particle coordinates to solve the Poisson equation.

In this paper, a simplified model as the modeling of magnetized plasmas is developed. The starting point is the assumption that the distribution of the ensemble of charged particles in the same species is homogeneous over gyrophase. The particles in this ensemble are located at the same guiding-cneter position $(\mathbf{X},\mu,U)$. This assumption could lead to a significant simplification of the gyrokinetic model through the following steps.

Firstly, the ensemble of particle in the same species and of the same guiding-center coordinate $(\mathbf{X},\mu,U)$ but with their $\theta$ homogeneously distributed over $(0,2\pi]$, can be treated as an ensemble of identical particles, which have the same guiding-center orbit.

Secondly, instead of using the Lagrangian differential 1-form of a test particle as down in CGM [3, 4], this paper implements a fundamental Lagrangian 1-form which determines the dynamics of all ions and electrons on particles’ coordinates. The electrostatic potential in this Lagrangian is originated from the mutual interactions between charged particle pairs.

Thirdly, according to Lie transform perturbation theory, by carrying out the pullback transform over this Lagrangian 1-form, a new Lagrangian 1-form on guiding-center coordinate is derived, with all gyrophase cancelled by the summation of the gyrophase for each particle included in the ensemble surrounding the guiding-field magnetic field line. This kind of summation will be called ensemble summation in this paper. Therefore, it doesn’t need to carry out an additional gyrocenter transform to reduce the gyrophase for each particle, which nevertheless are mutually cancelled by the ensemble summation.

Forthly, the new fundamental Lagrangian 1-form determines the dynamics of all ions and electrons on guiding-center coordinates. Then, new guiding-center kinetic models(GCKM) can be derived based on a coarse-grained scheme as the modeling of the new fundamental Lagrangian 1-form.

The rest of this paper is arranged as follows. In Sec.(2), the fundamental Lagrangian 1-form determining the dynamics of all ions and electrons is introduced and is modeled by the Distribution-Poisson models. The modeling procedure will be used to derive new GCKM on guiding center coordinates. In Sec.(3), the new 1-form on the new coordinates approximated up to the second order is derived by pulling the fundamental Lagrangian 1-form on particle’s coordinates back to the one on the new coordinates. In Sec.(4), with the same modelling method given by Sec.(2), the new 1-form on the new coordinates with second order approximation is modeled by new GCKM. In Sec.(5), the equation for the quasi-neutral condition is introduced on guiding-center coordinate. Sec.(6) is the simple introduction of the numerical application of this new model. Sec.(7) is dedicated to summary and discussion.

2 Modelling the fundamental Lagrangian 1-form on particle’s coordinates by Distribution-Poisson models

It’s well-known that the dynamics of a physical system can be determined by the Lagrangian of this system [20, 21, 22]. The force experienced by each entity in the system is given by the various potentials in the Lagrangian. In this paper, we focus on the electrostatic plasma including only one species of ion and electrons. In real physical systems, particles are located at different spatial positions. The electrostatic potential experienced by one charged particle with spatial coordinate $\mathbf{x}$ is

[TABLE]

where the summation is taken for all ions and electrons except the one whose coordinate is $\mathbf{x}$, as the superscript ′ indicates.

The fundamental Lagrangian 1-form is the summation of the Lagrangian 1-form for each electron and ion. Here, we only consider the electrostatic case. The magnetic field is the background field independent of all the particles. The case of electromagnetic perturbations, for which the perturbed magnetic field can not be treated as background field, will be considered in future work. The Lagrangian can be written as the summation of two parts as

[TABLE]

with

[TABLE]

[TABLE]

The reason for the factor $\frac{1}{2}$ in $\frac{{\phi\left({{{\bf{x}}_{oj}}}\right)}}{2}$ is that the potential between each pair has mutual contributions by the two particles. Eq.(4) can also be written as compactly

[TABLE]

$o\in\{i,e\}$ and $q_{i}=q$ and $q_{e}=-e$ are utilized here and will be adopted throughout this paper.

The Lagrangian given by Eq.(5) determines the dynamics of all ions and electrons included in this system. It denotes a completely autonomous system. The equations of motion for each ion and electron can be derived based on Euler-Lagrange equation for each pair of $(\mathbf{x}_{(i,e)j},\mathbf{v}_{(i,e)j})$. The electrostatic potential originates from the mutual interaction of each pair of charged particles. However, it’s nearly impossible to straightforwardly calculate the electrostatic potential based on mutual interaction, since there are so many particle pairs. Noting that the electrostatic force generated by each charged particle is of the inverse-square force form, it’s well-known that based on Gauss’s law, Poisson’s equation can directly associate the potential at a spatial point with the charge at that point as follows[23]

[TABLE]

Therefore, Poisson’s equation plus a boundary condition can be an alternative way to calculate the electrostatic potential particles feel. Eq.(6) is rigorously based on the inverse-square force. If the force generated by the charges changes, this formula also changes, e.g., if mutual interaction force is the magnetostatic force, the equation relating magnetic potential to the current is the Ampere’s Law. This will be used later to model the fundamental Lagrangian on guiding center by GCKM.

The solution of $\phi(\mathbf{x},t)$ given by Eq.(6) can be solved by integrating both sides of Eq.(6) for each charged particle with a given boundary condition. However, straightforward applications of Eq.(6) are almost impossible since there are so many particles in a real plasma system. An effective Distribution-Poisson model is used to replace the fundamental Lagrangian given by Eq.(5).

It’s well-known that to study the phase transition phenomena presented by the Ising model, single spins are replaced by cells containing multiple spins [24, 25, 26], since the correlation length of the spin fluctuation is much longer than the distance between two neighbouring spins at the phase transition point. As for the electrostatic potential generated by the accumulation of charges as described by Eq.(6), if the length scale of electrostatic potential denoted as $l_{p}$ is much longer than the mean distance between charged particles denoted as $l_{d}$, to determine the electrostatic potential, it’s not needed to know the specific position of each particle. With the same treatment of the Ising model, a coarse graining process is possible to divide the spatial space occupied by the plasma into small cells with $l_{p}\gg l_{c}\gg l_{d}$ satisfied, where $l_{c}$ is the size of the cell, as shown in Fig.(1). To determine the electrostatic potential through Poisson’s equation given in Eq.(6), we only need to know the accumulation of charges within each cell. For the theoretical analysis, $l_{c}$ can be much smaller than the length scale of the grid used in the Particle-In-Cell simulation[27].

Then, the spatial integration of Eq.(6) is given by the following formula

[TABLE]

which can be modeled by the following summation of the quantities within each cell

[TABLE]

$\mathbf{x}_{k}$ is the coordinate of the spatial center of the $k$th cell. $\phi_{k}(\mathbf{x}_{k},t)$ is the average electrostatic potential of the $k$th cell and is located at the spatial center of the cell. So the gradient operator $\nabla$ over $\phi_{k}(\mathbf{x}_{k},t)$ can be modeled by the middle-point difference. $\Delta V_{k}$ represents the volume of the $k$th cell. $\beta_{ik}$ and $\beta_{ek}$ are the charge accumulation for ions and electrons within the $k$th cell, respectively. Upon the length scale $l_{c}$, we make the following replacement

[TABLE]

$n_{ok}(\mathbf{x}_{k},t)$ is the average density of the $k$th cell for ions or electrons and located at the center of the $k$th cell. The location of the average electrostatic potential and average density is given by Fig.(1). The location of $\phi_{k}(\mathbf{x}_{k},t)$ and $n_{ok}(\mathbf{x}_{k},t)$ is different from that in Particle-In-Cell simulation, where the physical quantities such as the potential and charge density are all located at grid points[27]. We take the following approximation of Poisson’s equation for the $k$th cell

[TABLE]

This is the wanted edition of Poisson’s equation as the result of coarse graining. The cell’s length scale $l_{c}$ can shrink to be close to $l_{d}$, so that the approximate potential can be infinitely close to the real potential. For the $j$th particle located within the $k$th cell, the relation between the real potential $\phi\left({{{\bf{x}}_{j}}},t\right)$ and the approximate one ${\phi_{k}}(\mathbf{x}_{k},t)$ satisfies the following formula

[TABLE]

The time evolution of $n_{k}(\mathbf{x}_{k},t)$ is needed. $n_{k}(\mathbf{x}_{k},t)$ can be expressed as an integral of distribution function. The distribution of particles can be given by the Klimontovich distribution[28]

[TABLE]

The average density $n_{ok}$ within the $k$th cell is given by an integral

[TABLE]

The evolution of ${M}_{ok}\left(\mathbf{z}\right)$ is based on the equations of motion of charged particle, which is in turn determined by the Lagrangian 1-form. The electrostatic potential all the particles feel in the $k$th cell is approximated to be the same as ${\phi_{k}}(\mathbf{x}_{k},t)$. According to Eq.(11), by shrinking $l_{c}$ to be much smaller than the Larmor radius of ion or electron, we could use ${\phi_{k}}(\mathbf{x}_{k},t)$ to replace the real potential $\phi(\mathbf{x}_{j},t)$ felt by the $j$th particle at the $k$th cell. As discussed in Appendix.(A), the magnetic field is the background field independent of the motion of charged particles and is given in advance. So each particle experiences different magnetic field, since each particle in the $k$th cell has its own position. Then, the test Lagrangian 1-form describing the motion of the $j$th single particle for ions or electrons within the $k$th cell can be extracted out from the fundamental Lagrangian 1-form in Eq.(5) as

[TABLE]

Here, it should be noted that the factor $\frac{1}{2}$ before the electrostatic potential is removed. Eq.(14) gives the trajectory equations of particles at the $k$th cell with the approximate potential $\phi_{k}(\mathbf{x}_{k},t)$. Collisions, the source and the sink are neglected in this paper. So the evolution of ${M}_{ik}\left(\mathbf{z}\right)$ is governed by the Liouville equation

[TABLE]

The obtained $M_{ok}\left(\mathbf{z}\right)$ can be substituted back in Eq.(13) to get the average charge density of the $k$th cell. According to Klimontovich’s theory[28], the ensemble summation of $M_{ok}\left(\mathbf{z}\right)$ leads to a distribution function $f_{ok}(\mathbf{z})$ as a continuous function of argument $\bf{z}$. Then, Eq.(15) changes to be the Vlasov equation

[TABLE]

Usually, to solve electron’s Vlasov equation, the adiabatic approximation of the fluctuation density of electrons can be utilized[29] as

[TABLE]

where $n_{0e}(\mathbf{x})$ is the equilibrium density for electrons on particle’s coordinates. Dividing $f_{o}(\bf{z})$ into equilibrium and perturbation parts as $f_{i}(\mathbf{z})=f_{i0}(\mathbf{z})+\tilde{f}_{i}(\mathbf{z})$, the charge density of ions generated by the equilibrium part is canceled by the equilibrium part of electrons in the Poisson’s equation, the rest of which becomes

[TABLE]

So far, a group of equations comprising Eqs.(10,14,15)forming the Klimontovich-Poisson model or Eqs.(10,14,16) forming the Vlasov-Poisson model are derived to constitute a close system to model the real fundamental Lagrangian 1-form given by Eq.(5). The main difference between the two systems is that the former one utilizes a coarse graining scheme to model the electrostatic potential and charge density, and the electrostatic potential is obtained by solving the Poisson’s equation. When the length scale $l_{c}$ of the cell is small enough, the subscript $k$ of those equations can be deleted, so that these three equations can be treated as defined on continuous spatial space.

As pointed out previously, the modeling is based on the force between charged particles being inverse-square force, so that Poisson’s equation Eq.(10) can be derived to relate the potential to the charge density. If the force changes, Poisson’s equation should be changed accordingly. For example, if the force is of magnetic origin, the equation should be replaced by the Ampere’s law correspondingly. With the same principle, we could develop GCKM as the modeling of the fundamental Lagrangian 1-form on guiding-center coordinates.

3 The fundamental Lagrangian 1-form on guiding-center coordinates

The fundamental Lagrangian differential 1-form for all ions and electrons on particle’s coordinates is given by Eq.(5). In this subsection, a pullback transform is adopted to pulling the 1-form back to a new one on guiding-center coordinates with $\theta$ angle reduced from the whole dynamical system up to order $O(\varepsilon_{i}^{2})$ for ion and $O(\varepsilon_{e}^{2})$ for electrons. Here, the ordering parameters are given as ${\varepsilon_{i}}\equiv\frac{1}{{q{L_{0}}}}\sqrt{\frac{{2{\mu_{it}}{B_{0}}}}{{{m_{i}}}}}$ and ${\varepsilon_{e}}\equiv\frac{1}{{e{L_{0}}}}\sqrt{\frac{{2{\mu_{et}}{B_{0}}}}{{{m_{e}}}}}$, with $\mu_{it}$ and $\mu_{et}$ the magnetic moment of the thermal velocity for ions and electrons, respectively. Before carrying out the calculation, the three assumptions are firstly listed as follows.

(1). The first assumption is that the distribution is homogeneous over gyrophase on particle coordinates.

(2). The second one is that particles’ coordinates only experience guiding center transform, which is resulted directly from the fact (1). As will be shown later, only guiding-center transform could help reduce the gyrophase of all particles from the dynamics of the whole system.

(3). The third assumption is associated with removing the singularity faced by the potential function in the new coordinates. After coordinate transform, it’s inevitable that the spatial part of some particles’ new coordinates could locate at the same spatial position in new coordinate system as shown in Fig.(2). This would introduce singularity into the electrostatic potential, as two identical spatial coordinate appear at the denominate to make it equal zero. To remove the singularity, it’s assumed that the mutual potential between those particles located at the same spatial place are removed from the total potential. This assumption is equivalent to the one that the interactions between particles, whose new spatial coordinates after the coordinate transform are the same, are removed from the total potential in Eq.(1). The removed part occupies only a very small part of the total electrostatic potential.

3.1 Deriving the fundamental Lagrangian 1-form on the new coordinates by Cary-Littlejohn single-parameter Lie transform method

With the second assumptions, we can go on to derive the new fundamental Lagrangian 1-form on the new coordinates. As shown in Ref.([30]), the single-parameter coordinate transform of the coordinates of a single particle is

[TABLE]

[TABLE]

where not as given previously, $\mathbf{g}_{oj}\left({{{\bf{Z}}_{oj}}}\right)$ is normalized here with $\mathbf{g}_{oj}\left({{{\bf{Z}}_{oj}}}\right)\equiv-\bm{\rho}_{o}/\varepsilon_{o}$. The subscript $f$ in Eq.(19) denotes forward transform. Eq.(20) can be further written as

[TABLE]

where repeated indexes denote the Einstein summation. Eqs.(19) and (21) induce a coordinate transform

[TABLE]

Now, we make following definitions

[TABLE]

The following backward coordinate transformation can be derived

[TABLE]

Now, implementing the method in Ref.([30]) with Eq.(21), we first carry out the pullback transform for ions’ coordinates. The following equation can be derived

[TABLE]

which leads to the solution

[TABLE]

Next, the pullback transform for electrons’ coordinates is carried out for the 1-form in Eq.(25), and the following formula is derived

[TABLE]

which leads to the solution

[TABLE]

By substituting Eq.(25) in Eq.(27), the latter becomes

[TABLE]

which can be further written as

[TABLE]

for the simplicity. Since each coordinate pair $(\mathbf{x}_{(i,e)j},\mathbf{v}_{(i,e)j})$ is independent of all others, the operators ${L_{{\bf{g}}_{(i,e)j}^{\bf{X}}}}$ for each $j$ commutes.

3.2 Approximating Eq.(29) to the order of $O(\varepsilon_{o}^{2})$

Expanding Eq.(29) based on the order of $\varepsilon_{o}$ and $\varepsilon_{e}$, the eventual Lagrangian differential 1-form, can be derived. The following rules will be used

[TABLE]

[TABLE]

Here, $\mathbf{f}(\mathbf{Z})$ and $h(\mathbf{Z})$ are any vector function and scalar function on the new coordinates, respectively.

Among the expansions, the following 1-form is denoted as $\Gamma_{0}$

[TABLE]

Here, the factor $\varepsilon_{o}$ for $o\in\{i,e\}$ is used as a symbol to denote the order of the term adjacent to it. This terminology will be used throughout the remaining part of this paper. The electrostatic potential is

[TABLE]

where ′ denotes all particles located at $\mathbf{X}_{oj}$ are removed from the summation based on the third assumption.

Due to the homogeneous property of the distribution of particles over $\theta$, the summation in Eq.(32) cancels terms depending on the gyroangle, leading to the following formula

[TABLE]

with ${\widehat{\bf{v}}_{ol}}={\bf{e}}_{1}\sin{\theta_{ol}}+{{\bf{e}}_{2}}\cos{\theta_{ol}}$. Eq.(34) can be understood in another way

[TABLE]

The reason for the standing of Eq.(34) or Eq.(35) is as follows. Given any group $\{\mathbf{X},\mu,U\}$, we have a ensemble of ions or electrons, in which the spatial part of guiding-center coordinates, the parallel velocity and magnetic moment of each particle equal $\{\mathbf{X},\mu,U$}, respectively. The particles in this ensemble are homogeneously distributed over $\theta$ surrounding the guiding center $\mathbf{X}$. So, the summation of $\sin\theta_{ol}$ or $\cos\theta_{ol}$ for all $l$ equals zero, where subscript $l$ denotes particles of the same $\{\mathbf{X},U,\mu\}$. Then, $\Gamma_{0}$ becomes

[TABLE]

The next one is

[TABLE]

Eq.(43) contains a summation like $\sum\limits_{j}{{\bf{g}}_{oj}^{\bf{X}}\cdot{\nabla_{oj}}\Phi\left({{{\bf{X}}_{oj}},t}\right)}$, which can further be divided as the summation of different category ensemble, which is $\sum\limits_{l}{{\bf{g}}_{ol}^{\bf{X}}\cdot\nabla\Phi\left({{\bf{X}},t}\right)}$. Here, subscript $l$ denotes the ensemble of charged particle all located at the guiding-center point $(\mathbf{X},\mu,U)$. Due to the homogeneous property, for any generator vector ${{\bf{g}}_{oh}^{\bf{X}}}$ for the particle $j$ in this ensemble, there always exists a particle denoted by subscript $k$ in this ensemble, the generator ${{\bf{g}}_{ok}^{\bf{X}}}$ of which equals $-{{\bf{g}}_{oh}^{\bf{X}}}$. Therefore, the following identity can be derived

[TABLE]

To calculate Eq.(43), we also need another identity

[TABLE]

which results from ${\bf{g}}_{oj}^{\bf{X}}\bot{\widehat{\bf{v}}_{oj\bot}}$. Eventually, it’s derived out that only terms like ${L_{{\bf{g}}_{oj}^{\bf{X}}}}\left({\sqrt{\frac{{2B({{\bf{X}}_{oj}}){\mu_{oj}}}}{{{m_{o}}}}}{{\widehat{\bf{v}}}_{oj\bot}}\cdot d{{\bf{X}}_{oj}}}\right)$ in Eq.(43) can generate non-zero results. The summation of these terms leads to

[TABLE]

The next one is

[TABLE]

We only keep the lower order part of $\Gamma_{2}$. The lower order part of $\gamma$ is written as

[TABLE]

[TABLE]

To calculate Eq.(47), it needs to be noted that only the terms with $n=o$ and $j=h$ in $\Gamma_{2}$ can produce nonzero terms based on the homogeneous assumption in the $\theta$ direction. The following formula is needed as well

[TABLE]

To derive the approximate equality in Eq.(50), the formula $\varepsilon_{o}{\bf{g}}_{oj}^{\bf{X}}\times\nabla\times\left({{\bf{g}}_{oj}^{\bf{X}}\times{\bf{B}}\left({{{\bf{X}}_{oj}}}\right)}\right)\cdot d{{\bf{X}}_{oj}}$ is neglected in the $\mathbf{X}_{oj}$ component in the Lagrangian 1-form, since its order is $O(\varepsilon_{o})$, while the terms in the $\mathbf{X}_{oj}$ component in $\Gamma_{0}$ are of order $O(1/\varepsilon_{o})$. In fact, this term can also be cancelled by introducing a generator $\mathbf{g}_{2oj}^{\mathbf{X}}$ of order $\varepsilon_{o}^{2}$ for each $j$. To derive the second equality in Eq.(50), the equation ${\widehat{\bm{\rho}}_{0}}\times{\widehat{\bf{v}}_{\bot}}=-\widehat{\bf{b}}$ is adopted.

Eventually, the rest of $\Gamma_{2}$ is

[TABLE]

In Eq.(51), the factor $\frac{1}{2}$ before $\Phi(\mathbf{X}_{oj},t)$ is removed by combining all potential depending on $\mathbf{X}_{oj}$ together. The third identity in Eq.(51) also comes from the homogeneous assumption. At last, combining Eqs.(36,46,51) together, we could derive the following fundamental Lagrangian 1-form defined on the new coordinates up to the second order approximation

[TABLE]

with

[TABLE]

[TABLE]

All the symbols $\varepsilon_{o}$ are removed from Eq.(52,53). $\bar{\Gamma}$ is the new fundamental Lagrangian 1-form determining the dynamics of all ions and electrons on the new coordinate system up to the second order approximation. $\Pi_{o}\left({{\bf{X}}_{oj},\mu_{oj},t}\right)$ is a FLR term and could introduce the difference to the trajectory equations compared with those equations derived without $\Pi_{o}\left({{\bf{X}}_{oj},\mu_{oj},t}\right)$.

4 Modelling the fundamental Lagrangian 1-form on guiding-center coordinates by GCKM

Eq.(52) is a model based on first-principle force. It becomes untractable by increasing the particle number $N$. A way to simplify Eq.(52) is to adopt the same modeling method used in Sec.(2) by degenerating the pair-wise Coulomb force to the Poisson equation plus a boundary condition. The potential $\Phi\left({{{\bf{X}}_{oj}}},t\right)$ experienced by the particle located at $\mathbf{X}_{oj}$ is given by Eq.(33). The knowledge of $\Pi_{o}\left({{\bf{X}}_{oj},\mu_{oj},t}\right)$ depends on $\Phi\left({{{\bf{X}}_{oj}}},t\right)$. The formalism of $\Phi\left({{{\bf{X}}_{oj}}},t\right)$ makes sure that we can use Poisson’s equation plus a boundary condition to model it. Similar to Eq.(10), the relation between $\Phi\left({{{\bf{X}}_{oj}}},t\right)$ and the local charge density can be written as

[TABLE]

${{N_{ek}}({{\bf{X}}_{k}},t)}$ and ${{N_{ik}}({{\bf{X}}_{k}},t)}$ are the average density of electrons and ions of the $k$th cell on the new coordinates and located at the spatial center of the $k$th cell. $\Phi_{k}({{\bf{X}}_{k}},t)$ is the average electrostatic potential at the center of the $k$th cell. The sketch map is given in Fig.(3)

Just as Eq.(14), with the modeled electrostatic potential $\Phi_{k}(\mathbf{X}_{k},t)$, the dynamics of the $j$th particle in the $k$th cell on the new coordinate can be described by the following test Lagrangian 1-form

[TABLE]

with

[TABLE]

[TABLE]

where the second order derivative ${\nabla^{2}}\Phi_{k}\left({{{\bf{X}}_{k}},t}\right)$ is given by a middle-point discrete derivative with $\mathbf{X}_{k}$ as the center. Here, the factor $\frac{1}{2}$ before the electrostatic potential is removed, since all the mutual interactions depending on $\mathbf{X}_{oj}$ are combined together and is approximated by ${\Phi_{k}}\left({{{\bf{X}}_{k}},t}\right)$. Eq.(56) determines the trajectory equations on the new coordinates. ${\Pi_{ok}}\left({{{\bf{X}}_{k}},{\mu_{oj}},t}\right)$ denotes the FLR term introduced to the trajectory equations.

Now, we need to calculate the evolution of the density of charged particles in the $k$th cell. To do this, the Klimontovich distribution on the coordinate frame $\mathbf{Z}^{\prime}=\{\mathbf{X},\mu,U\}$ is needed

[TABLE]

The independence of $\theta$ is naturally inherited by $\mathbb{M}{{}_{ok}}\left({\mathbf{Z}^{\prime}}\right)$ from the fundamental Lagrangian 1-form Eq.(52) which is independent of $\theta$.

The evolution of $\mathbb{M}{{}_{ok}}\left(\mathbf{Z}^{\prime}\right)$ is given by the Liouville’s equation

[TABLE]

By ensemble summation of $\mathbb{M}{{}_{ok}}$, the Vlasov distribution ${F}_{o}(\mathbf{Z}^{\prime})$ can be derived and Eq.(60) becomes the Vlasov equation

[TABLE]

By shrinking $l_{c}$ to be small enough (much smaller than the Larmor radius of ions), the subscript $k$ can be removed from Poisson’s equation, the Lagrangian 1-form of a test particle, and the Vlasov equation, all of which are rewritten as

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

So far, we derived a guiding-center Klimontovich-Poisson model comprising Eqs.(62,63,4) and a guiding-center Vlasov-Poisson model comprising Eqs.(62-64) as the modeling of the fundamental Lagrangian 1-form in Eq.(52). These two models constitute a close system, respectively, based on which all theoretical analysis and simulations can be done.

5 The equation of quasi-neutral condition in GCKM

The new guiding-center Vlasov-Poisson model is taken as the example to show the equation for the quasi-neutral condition on guiding-center coordiantes. In this paper, the electrostatic potential with $l_{p}\gg\rho_{e}$ while $l_{p}\approx\rho_{i}$ is taken into account. The FLR term for electrons in Eq.(66) denoted by $\Pi\left({{\bf{X}},{\mu_{ej}},t}\right)$ can be ignored, due to the much smaller Larmor radius of electrons. The adiabatic approximation of the fluctuation density of electrons can be written as

[TABLE]

with ${N}_{e0}(\bf{X})$ being the equilibrium density for electrons on the new coordinates. The FLR term denoted by $\Pi\left({{\bf{X}},{\mu_{oj}},t}\right)$ for ions is kept. The distribution function of ions on new coordinates can be decomposed as [29, 31]

[TABLE]

with the perturbed part being

[TABLE]

Here, $H\left(\mathbf{Z}\right)$ is the non-adiabatic distribution. The perturbed density is derived by integrating over the distribution function

[TABLE]

where $\rho_{t}$ is the Larmor radius with the thermal velocity. The equation for the quasi-neutral condition becomes

[TABLE]

The term proportional to $\rho_{t}^{2}{\nabla^{2}}$ is a kind of the polarization density as explained by Eq.(38) of Ref.([7]). But for the practical application, the edition of the equation for the quasi-neutral condition is

[TABLE]

For the numerical application of Eq.(72), $F_{i1}$ can be straightforwardly calculated from the Vlasov equation based on the potential given by the last time step. To get the solution of the new potential on the current step, it’s not needed to solve operator “$\nabla^{2}$” which appears in the equation for the quasi-neutral equation in CGM.

6 Applications of GCKM

For numerical applications of GCKM, the simulations are carried out totally on the new coordinate. Except at the beginning and at the end of the simulations, the calculations of Fourier spectrum of the perturbations and the average of the gyroangle are exempted at each time step. So the numerical time and numerical instabilities can be significantly reduced.

For the simulation based the Vlasov distribution, the initial distribution of ${{F}_{o}}\left({{\bf{X}},\mu,U,t}\right)$ needs to be transformed from the one on particle’s coordinates. If we know an initial distribution function $f_{oin}(\mathbf{x},u_{1},\mu_{1},t)$ for ions on particle’s coordinates, the initial distribution on $\mathbf{Z}^{\prime}$ is given by

[TABLE]

In fact, Eq.(6) is the first order approximation of Eq.(86). The transform given by Eq.(6) naturally makes ${F}_{iin}$ inherit the perturbation wave from $f_{iin}$. Using the initial distribution $F_{iin}\left(\mathbf{Z}\right)$, Poisson’s equation Eq.(62), and trajectory equations derived from Eq.(63), the time evolution of the distribution in Eq.(64) can be calculated to get $F_{iend}\left(\mathbf{Z}\right)$, based on which various quantities, e.g, transport of number density, momentum, energy, and amplitude of potential fluctuation, can be derived on new coordinate $\bf{Z}$. ${F}_{iend}\left(\mathbf{Z}\right)$ can also be transformed back to the one on particle’s coordinates by the following formula

[TABLE]

The perturbation wave included by $\mathbb{F}_{iend}(\mathbf{Z})$ is naturally transmitted back to $f_{iend}(\mathbf{z})$ by Eq.(79).

For particle-in-cell simulations, all equations needed are the trajectory equations derived from Eq.(63), the Poisson’s equation given by Eq.(62) and the evolution equation Eq.(60) for the Klimontovich distribution. The four-point scheme to calculate the density on particle coordinate from the one on the gyrocenter coordinate is exempted at each time step. Eqs.(6) and (79) may be applied to transform the distribution between $\bf{z}$ and $\bf{Z}$ at the beginning and at the end of the simulation.

7 Summary and discussion

This paper presented an assumption of a homogeneous distribution of the ensemble of charged particles over the gyrophase. As a company to this assumption, this paper developed a fundamental Lagrangian 1-form, which contains all particles in the magnetized plasma system as well as the Coulomb pair force between particles instead of field-particle interaction used in conventional gyrokinetic models. This fundamental Lagrangian 1-form perfectly makes use of the property of the homogeneous assumption to reduce the gyrophase of each particle. Therefore, it doesn’t need an additional gyrocenter transform as used in CGM, which induces a polarization density in the quasi-neutral equation in CGM. However, such a polarization density doesn’t appear in GCKM. It’s expected that the numerical application of GCKM could reduce the numerical noise level and guarantee long-term simulation.

In this paper, the fundamental Lagrangian 1-form is only approximated up to the second order $O(\varepsilon_{o}^{2})$. High order approximation can be carried out for specific problems. This paper only considers electrostatic perturbations, whilst the magnetic vector potential is treated as the background field. For magnetic perturbations, such an operation is improper. The non-equilibrium part of the magnetic vector potential needs special treatment similar to the treatment of the electrostatic potential in this paper. This is left for the future work.

8 Acknowledgments

This work was completed at Uji campus of Kyoto University, Japan. The author is indebted to the discussion with Prof. Weixing Wang, Prof. Yasuaki Kishimoto,Prof. Kenji Imadera, Prof. Alain Brizard, Prof. T.S.Hahm, Prof. Guoyong Fu, in particular, is grateful to Prof. Johan Anderson for reading the manuscript.

Appendix A Background field, non-background field and scalar under the coordinate transform

The background field is given by the magnetic vector potential $\mathbf{A}(\mathbf{x})$, which doesn’t explicitly dependent on time. The other field is the non-background field which is the electrostatic potential generated by the separation of ions and electrons. It should be noted that $\Phi(\mathbf{X},t)$ given by Eq.(33) is different from $\phi(\mathbf{x},t)$ given by Eq.(1) in the following way. $\phi(\mathbf{x},t)$ and $\Phi(\mathbf{X},t)$ both depend on the spatial coordinate of all particles in the respective spatial space. The coordinate transform which induces the pullback transform Eq.(29) is approximated as

[TABLE]

${\bf{g}}_{oj}^{{{\bf{X}}}}$ depends on $\mu_{oj}$. If the RHS of Eq.(80) is substituted back to $\phi(\mathbf{x},t)$, it’s found that $\phi(\mathbf{x},t)$ depends on $(\mathbf{X}_{oj},\mu_{oj})$ for all $j$ and $o\in\{i,e\}$ except the new coordinate of this particle. However, $\Phi(\mathbf{X},t)$ only depends on $\mathbf{X}_{oj}$ for all $o$ and $j$, not on any $\mu_{oj}$. This also explains why $\phi(\mathbf{x},t)$ is not a scalar under coordinate transform $\psi$.

The fact that $\phi(\mathbf{x},t)$ is not a scalar can also be explained as follows. As introduced in Sec.(1), the coordinate transform $\psi$ is defined on the phase space as $\psi:\mathbf{z}\equiv(\mathbf{x},\mathbf{v})\to\mathbf{Z}\equiv(\mathbf{X},\mu,U,\theta)$. The Cary-Littlejohn single-parameter Lie transform theory[30] shows that for a coordinate transform described by a group of autonomous equations like

[TABLE]

and

[TABLE]

if there exists a scalar $s(\bf{z})$ with its counterpart in the new coordinate transform being $S(\mathbf{Z})$, then, from the identity $s\left({\bf{z}}\right)=S\left({{\bf{Z}},\varepsilon}\right)$, the following formula can be derived

[TABLE]

Now, if we assume that $\phi(\mathbf{x},t)$ is a scalar under the coordinate transform, and its new form is denoted as $\Phi(\mathbf{X},t,\varepsilon)$ in the new coordinate, similar to Eq.(83), we could derive an equation

[TABLE]

In Eq.(84), ${g^{\mu}}\left({\bf{Z}}\right)$ is a function of variables $(\mathbf{X},\mu,\theta)$. But $\Phi(\bf{X},t)$ is only a function over $\mathbf{X},t$. Therefore, Eq.(84) is not valid. The essential reason is that $\phi(\bf{x},t)$ is not a scalar under coordinate transform.

However, the distribution is a scalar, the transform of which is

[TABLE]

The solution of Eq.(85) is

[TABLE]

Appendix B Energy conservation of the system given by Eq.(52)

The trajectory equations derived from Eq.(52) are

[TABLE]

[TABLE]

[TABLE]

with ${\bf{B}}_{oj}^{*}={\nabla_{oj}}\times\left({{\bf{A}}\left({{{\bf{X}}_{oj}}}\right)+{U_{oj}}{\bf{b}}}\right)$. The energy per particle is formally written as

[TABLE]

Now, we try to prove $\frac{d}{{dt}}\sum\limits_{o,j}{{H_{oj}}}=0$. It’s first to derive the following three identities

[TABLE]

[TABLE]

[TABLE]

In Eq.(93), the property that ${\Psi\left({{{\bf{X}}_{oj}},{\mu_{oj}}}\right)}$ doesn’t explicitly depend on the time is applied. This is basic property of the electrostatic fluctuation. It’s easy to observed that term (1.1) cancels (1.2), term (2.1) cancels (2.2). Term (3.1) cancels (3.2) based on the following two identities

[TABLE]

[TABLE]

Then, the conservation of the total energy of the system is proved.