Asymptotically preserving particle-in-cell methods for inhomogenous strongly magnetized plasmas

TL;DR

This paper introduces semi-implicit Particle-In-Cell methods for simulating plasmas with inhomogeneous strong magnetic fields, effectively capturing guiding-center dynamics without restrictive time step constraints.

Contribution

The paper develops and validates a novel class of PIC methods that remain stable and accurate in regimes with large magnetic fields and inhomogeneity, extending previous approaches to more complex plasma conditions.

Findings

The methods accurately reproduce guiding-center motion in inhomogeneous magnetic fields.

Numerical experiments confirm stability and consistency of the proposed schemes.

The approach allows larger time steps without sacrificing accuracy.

Abstract

We propose a class of Particle-In-Cell (PIC) methods for the Vlasov-Poisson system with a strong and inhomogeneous external magnetic field with fixed direction, where we focus on the motion of particles in the plane orthogonal to the magnetic field (so-called poloidal directions). In this regime, the time step can be subject to stability constraints related to the smallness of Larmor radius and plasma frequency. To avoid this limitation, our approach is based on first and higher-order semi-implicit numerical schemes already validated on dissipative systems [3] and for homogeneous magnetic fields [10]. Thus, when the magnitude of the external magnetic field becomes large, this method provides a consistent PIC discretization of the guiding-center system taking into account variations of the magnetic field. We carry out some theoretical proofs and perform several numerical experiments that establish a solid validation of the method and its underlying concepts.

Figures

Equations

$${\mathbf{B}}(t,{\mathbf{x}})\,\,=\,\,\frac{1}{\varepsilon}\,\left(\begin{array}[]{l}0\\ 0\\ b(t,{\mathbf{x}})\end{array}\right),$$

$$\left\{\begin{array}[]{l}\displaystyle{\varepsilon\frac{\partial f^{\varepsilon}}{\partial t}\,+\,\mathbf{v}\cdot\nabla_{\mathbf{x}}f^{\varepsilon}\,+\,\left(\mathbf{E}^{\varepsilon}(t,{\mathbf{x}})\,-\,b(t,{\mathbf{x}})\frac{{\mathbf{v}}^{\perp}}{\varepsilon}\right)\cdot\nabla_{\mathbf{v}}f^{\varepsilon}\,=\,0,}\\ \\ \displaystyle{{\mathbf{E}}^{\varepsilon}=-\nabla_{\mathbf{x}}\phi^{\varepsilon},\quad-\Delta_{{\mathbf{x}}}\phi^{\varepsilon}=\rho^{\varepsilon},\quad\rho^{\varepsilon}=\int_{\mathbb{R}^{2}}f^{\varepsilon}d{\mathbf{v}},}\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle{\varepsilon\frac{\textrm{d}{\mathbf{X}}^{\varepsilon}}{\textrm{d}t}\,=\,{\mathbf{V}}^{\varepsilon},}\\ \\ \displaystyle{\varepsilon\frac{\textrm{d}{\mathbf{V}}^{\varepsilon}}{\textrm{d}t}\,=\,-\frac{1}{\varepsilon}\,{b}(t,{\mathbf{X}}^{\varepsilon})\,({\mathbf{V}}^{\varepsilon})^{\perp}\,+\,{\mathbf{E}}^{\varepsilon}(t,{\mathbf{X}}^{\varepsilon}),}\\ \\ {\mathbf{X}}^{\varepsilon}(t^{0})={\mathbf{x}}^{0},\,{\mathbf{V}}^{\varepsilon}(t^{0})={\mathbf{v}}^{0},\end{array}\right.$$

$$f^{\epsilon }(t,{\mathbf{X}}^{\epsilon }(t),{\mathbf{V}}^{\epsilon }(t))=f^{\epsilon }(t^0,{\mathbf{x}}^0,{\mathbf{v}}^0)$$

$$b(t,\mathbf{x})\ge b_0\forall (t,\mathbf{x})\in {\mathbb{R}}^{+}\times {\mathbb{R}}^2.$$

$$\frac{\text{d}}{\text{d}t}\left({\mathbf{X}}^{\epsilon }-\epsilon \frac{({\mathbf{V}}^{\epsilon }{)}^{\perp }}{b(t,{\mathbf{X}}^{\epsilon })}\right)=\mathbf{F}(t,{\mathbf{X}}^{\epsilon })-({\mathbf{V}}^{\epsilon }{)}^{\perp }\left[{\text{d}}_{\mathbf{x}}\left(\frac{1}{b}\right)(t,{\mathbf{X}}^{\epsilon })({\mathbf{V}}^{\epsilon })+\epsilon \frac{\partial }{\partial t}\left(\frac{1}{b}\right)(t,{\mathbf{X}}^{\epsilon })\right]$$

$$\mathbf{F}(t,\mathbf{x}):=-\frac{1}{b(t,\mathbf{x})}{\mathbf{E}}^{\perp }(t,\mathbf{x}),\forall (t,\mathbf{x})\in {\mathbb{R}}^{+}\times {\mathbb{R}}^2.$$

$$\frac{\text{d}}{\text{d}t}{\mathbf{X}}^{\epsilon }$$

$$\frac{\text{d}\mathbf{Y}}{\text{d}t}=\mathbf{F}(t,\mathbf{Y}).$$

$$\frac{\text{d}}{\text{d}t}\left(\frac{1}{2}\Vert {\mathbf{V}}^{\epsilon }{\Vert }^2-\epsilon \mathbf{F}(t,{\mathbf{X}}^{\epsilon })\cdot {\mathbf{V}}^{\epsilon }\right)=-{\mathbf{V}}^{\epsilon }\cdot {\text{d}}_{\mathbf{x}}\mathbf{F}(t,{\mathbf{X}}^{\epsilon })({\mathbf{V}}^{\epsilon })-\epsilon {\mathbf{V}}^{\epsilon }\cdot \frac{\partial \mathbf{F}}{\partial t}(t,{\mathbf{X}}^{\epsilon }),$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{1}{2}\,\frac{\textrm{\emph{d}}}{\textrm{\emph{d}}t}\left(\,|v^{\varepsilon}_{1}|^{2}-|v^{\varepsilon}_{2}|^{2}\,\right)\,=\,\frac{b(t,{\mathbf{X}}^{\varepsilon})}{\varepsilon^{2}}\,2v^{\varepsilon}_{1}\,v^{\varepsilon}_{2}\,+\,\frac{v^{\varepsilon}_{1}\,E_{1}(t,{\mathbf{X}}^{\varepsilon})-v_{2}\,E_{2}(t,{\mathbf{X}}^{\varepsilon})}{\varepsilon},\\ \\ \displaystyle\frac{1}{2}\frac{\textrm{\emph{d}}}{\textrm{\emph{d}}t}\left(\,2v^{\varepsilon}_{1}\,v^{\varepsilon}_{2}\;\right)\,=\,-\,\frac{b(t,{\mathbf{X}}^{\varepsilon})}{\varepsilon^{2}}\,\left(\,|v^{\varepsilon}_{1}|^{2}\,-\,|v^{\varepsilon}_{2}|^{2}\,\right)\,+\,\frac{v_{2}^{\varepsilon}\,E_{1}(t,{\mathbf{X}}^{\varepsilon})+v_{1}^{\varepsilon}\,E_{2}(t,{\mathbf{X}}^{\varepsilon})}{\varepsilon}\,.\end{array}\right.$$

$$\frac{{\mathbf{V}}^{\epsilon }}{\epsilon }=\mathbf{F}(t,{\mathbf{X}}^{\epsilon })+\frac{\epsilon }{b(t,{\mathbf{X}}^{\epsilon })}{\left(\frac{\text{d}{\mathbf{V}}^{\epsilon }}{\text{d}t}\right)}^{\perp }.$$

$$\frac{1}{2}\frac{\text{d}}{\text{d}t}\Vert {\mathbf{V}}^{\epsilon }{\Vert }^2$$

$$-{\mathbf{V}}^{\epsilon }{\text{d}}_{\mathbf{x}}\left(\frac{1}{b}\right)(t,{\mathbf{X}}^{\epsilon })({\mathbf{V}}^{\epsilon })$$

$$-g{\nabla }_{\mathbf{x}}\left(\frac{1}{b}\right)(t,\mathbf{Y})$$

$$\frac{\text{d}\mathbf{Y}}{\text{d}t}=\mathbf{F}(t,\mathbf{Y})+g\frac{{\nabla }_{\mathbf{x}}^{\perp }b}{b^2}(t,\mathbf{Y}).$$

$$\frac{\text{d}g}{\text{d}t}=-{\text{div}}_{\mathbf{x}}(\mathbf{F})(t,\mathbf{Y})g.$$

$$\left\{\begin{array}[]{l}\displaystyle{\frac{\textrm{d}{\mathbf{Y}}}{\textrm{d}t}\,=\,{\mathbf{F}}(t,{\mathbf{Y}})\,+\,g\,\frac{\nabla_{\mathbf{x}}^{\perp}b}{b^{2}}(t,{\mathbf{Y}})},\\ \\ \displaystyle{\frac{\textrm{d}g}{\textrm{d}t}\,=\,-\textrm{div}_{\mathbf{x}}({\mathbf{F}})(t,{\mathbf{Y}})\,g},\\ \\ {\mathbf{Y}}(t^{0})={\mathbf{x}}^{0},\quad g(t^{0})\,=\,e^{0},\end{array}\right.$$

$$\frac{\text{d}}{\text{d}t}\left(\frac{g}{b(t,\mathbf{Y})}\right)=-\frac{g}{b(t,\mathbf{Y}{)}^2}\left({\partial }_{t}b+\text{curl}(\mathbf{E})\right)(t,\mathbf{Y})$$

$${\partial }_{t}b+\text{curl}(\mathbf{E})=0$$

$$f^{\epsilon }:(t,\mathbf{x},e)\mapsto \frac{1}{2\pi }{\int }_0^{2\pi }{f}^{\epsilon }(t,\mathbf{x},\sqrt{2e}(\cos (\theta ),\sin (\theta )))\text{d}\theta$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial f}{\partial t}+\mathbf{U}\cdot\nabla_{{\mathbf{y}}}f+u_{g}\frac{\partial f}{\partial g}=0,\\ \\ \displaystyle-\Delta_{{\mathbf{y}}}\phi=\rho\,,\quad\rho=2\pi\,\int_{\mathbb{R}^{+}}f\,\textrm{d}g,\end{array}\right.$$

$$\mathbf{U}(t,\mathbf{y},g)=\mathbf{F}(t,\mathbf{y})+g\frac{{\nabla }_{\mathbf{x}}^{\perp }b}{b^2}(t,\mathbf{y}),u_g=-{\text{div}}_{\mathbf{y}}(\mathbf{F})(t,\mathbf{y})g,$$

$${\int }_{{\mathbb{R}}^2\times {\mathbb{R}}^{+}}g{f}^0(\mathbf{y},g)\text{d}\mathbf{y}\text{d}g<\infty .$$

$$\Vert f(t,\cdot ,\cdot ){\Vert }_{L^p}\le \Vert f^0{\Vert }_{L^p}\forall t\in {\mathbb{R}}^{+}$$

$$\mathcal{E}(t):=2\pi {\int }_{{\mathbb{R}}^2\times {\mathbb{R}}^{+}}gf(t,\mathbf{y},g)\text{d}\mathbf{y}\text{d}g+\frac{1}{2}{\int }_{{\mathbb{R}}^2}\Vert \mathbf{E}(t,\mathbf{y}){\Vert }^2\text{d}\mathbf{y}\le \mathcal{E}(0)\forall t\in {\mathbb{R}}^{+}.$$

$$f(t,\mathbf{Y}(t),g(t))=f(t^0,{\mathbf{x}}^0,e^0)\forall (t,t^0,{\mathbf{x}}^0,e^0)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}\times {\mathbb{R}}^2\times {\mathbb{R}}^{+}$$

$${\text{div}}_{\mathbf{y}}\mathbf{U}+\frac{\partial u_g}{\partial g}={\text{div}}_{\mathbf{y}}(\mathbf{F})-{\text{div}}_{\mathbf{y}}(\mathbf{F})=0.$$

$$f_N(t,\mathbf{x},\mathbf{v}):=\sum _{1\le k\le N}{\omega }_k{\phi }_{\alpha }(\mathbf{x}-{\mathbf{X}}_k(t)){\phi }_{\alpha }(\mathbf{v}-{\mathbf{V}}_k(t)),$$

$$\left\{\begin{array}[]{l}\displaystyle\varepsilon\frac{\textrm{d}{\mathbf{X}}_{k}}{\textrm{d}t}\,=\,{\mathbf{V}}_{k},\\ \\ \displaystyle\varepsilon\frac{\textrm{d}{\mathbf{V}}_{k}}{\textrm{d}t}\,=\,-\frac{1}{\varepsilon}\,b(t,{\mathbf{X}}_{k})\,{\mathbf{V}}_{k}^{\perp}\,+\,{\mathbf{E}}(t,{\mathbf{X}}_{k}),\\ \\ {\mathbf{X}}_{k}(0)={\mathbf{x}}^{0}_{k},\quad{\mathbf{V}}_{k}(0)={\mathbf{v}}^{0}_{k},\end{array}\right.$$

Full text

Keywords. High-order time discretization; Vlasov-Poisson system; Inhomogeneous magnetic field; Particle methods.

1. Introduction

We consider a magnetic field acting in the vertical direction and only depending on time $t\geq 0$ and the first two components ${\mathbf{x}}=(x_{1},x_{2})\in\mathbb{R}^{2}$, that is,

[TABLE]

with $b(t,{\mathbf{x}})>0$, where $\varepsilon>0$ is a small parameter related to the ratio between the reciprocal Larmor frequency and the advection time scale. This magnetic field is applied to confine a plasma constituted of a large number of charged particles, modelled by a distribution function $f^{\varepsilon}$ solution to the Vlasov equation coupled with the Poisson equation satisfied by the self-consistent electrical potential $\phi^{\varepsilon}$. Here we focus on the long time behavior of the plasma in a plane orthogonal to the external magnetic field. Therefore we consider the following form of the two dimensional Vlasov-Poisson system with an external strong magnetic field

[TABLE]

where we use notation ${\mathbf{v}}^{\perp}=(-v_{2},v_{1})$ and for simplicity we set all immaterial physical constants to one. The term $\varepsilon$ in front of the time derivative of $f^{\varepsilon}$ stands for the fact that we want to follow the solution on times sufficiently large to observe non trivial averaged evolutions.

We want to construct numerical solutions to the Vlasov-Poisson system (1.1) by particle methods (see [2]), which consist in approximating the distribution function by a finite number of macro-particles. The trajectories of these particles are determined from the characteristic curves corresponding to the Vlasov equation

[TABLE]

where we use the conservation of $f^{\varepsilon}$ along the characteristic curves

[TABLE]

and the electric field is computed from a discretization of the Poisson equation in (1.1) on a mesh of the physical space.

Before describing and analyzing a class of numerical methods for the Vlasov-Poisson system (1.1) in the presence of a strong external inhomogeneous magnetic field, we first briefly expound what may be expected from the continuous model in the limit $\varepsilon\to 0$.

1.1. Asymptotics with external electromagnetic fields

To gain some intuition, we first discuss the case when all electromagnetic fields $({\mathbf{E}}^{\varepsilon},b)$ are known. Conclusions of the present section may actually be completely justified analytically. Yet we do not pursue this line of exposition here to keep the discussion as brief and tight as possible.

Explicitly, in the present section we consider given (independent of $\varepsilon$) smooth electromagnetic fields $({\mathbf{E}}^{\varepsilon},b)$ and we assume that $b$ is bounded below away from zero, that is, we assume that there exists $b_{0}>0$ such that

[TABLE]

In the limit $\varepsilon\to 0$ one expects oscillations occurring on typical time scales $O(1/\varepsilon^{2})$ to coexist with a slow dynamics evolving on a time scale $O(1)$. We sketch now how to identify a closed system describing at main order the slow evolution. To begin with, note that from the second line of system (1.2) it does follow that ${\mathbf{V}}^{\varepsilon}$ oscillates at order $1/\varepsilon^{2}$ thus remains bounded and converges weakly111Though we do not want to be too precise here, let us mention that in the present discussion weakly and strongly refer to the weak-* and strong topologies of $L^{\infty}$ and that the weak convergences that we encounter actually correspond to strong convergence in $W^{-1,\infty}$. to zero. As we detail below, one may also combine both lines of the system to obtain

[TABLE]

where the force field ${\mathbf{F}}$ is the classical electric drift, given by

[TABLE]

Indeed

[TABLE]

This shows that ${\mathbf{X}}^{\varepsilon}$ evolves slowly but, as such, does not provide a closed asymptotic evolution in the limit $\varepsilon\to 0$. Indeed, except in the case when $b$ is constant dealt with in [10] we also need to know what happens to expressions that are quadratic in ${\mathbf{V}}^{\varepsilon}$ and this does not follow readily from the weak convergence of ${\mathbf{V}}^{\varepsilon}$.

Remark 1.1**.**

Incidentally note that, in contrast, when $b$ is constant, as in [10], the formulation of equation (1.4) is already sufficient to identify a guiding center evolution and obtain that ${\mathbf{X}}^{\varepsilon}$ converges to ${\mathbf{Y}}$ solving

[TABLE]

The missing piece of information may be intuited from

Lemma 1.2**.**

Consider ${\mathbf{X}}^{\varepsilon}=(x_{1}^{\varepsilon},x_{2}^{\varepsilon})$ and ${\mathbf{V}}^{\varepsilon}=(v_{1}^{\varepsilon},v_{2}^{\varepsilon})$ solving (1.2). Then we have

[TABLE]

and with ${\mathbf{E}}=(E_{1},E_{2})$

[TABLE]

Proof.

First note that the second equation in (1.2) may also be written as

[TABLE]

Therefore

[TABLE]

Hence (1.6). Likewise we may obtain the first equation in (1.7) — for $(|v^{\varepsilon}_{1}|^{2}-|v^{\varepsilon}_{2}|^{2})$ — by multiplying the second equation of (1.2) by $(v^{\varepsilon}_{1},-v^{\varepsilon}_{2})/\varepsilon$ and the equation for $v^{\varepsilon}_{1}\,v^{\varepsilon}_{2}$ by multiplying it by $(v^{\varepsilon}_{2},v^{\varepsilon}_{1})/\varepsilon$. ∎

This crucial lemma suggests that

• •

the microscopic kinetic energy $\frac{1}{2}\,\|{\mathbf{V}}^{\varepsilon}\|^{2}$ also evolves on a slow scale;

• •

the two variables $2v^{\varepsilon}_{1}\,v^{\varepsilon}_{2}$ and $|v^{\varepsilon}_{1}|^{2}\,-\,|v^{\varepsilon}_{2}|^{2}$ oscillate at scale $1/\varepsilon^{2}$ thus converge weakly to zero.

Indeed, the couple $\left(2v_{1}^{\varepsilon}\,v_{2}^{\varepsilon}\,,\,|v_{1}^{\varepsilon}|^{2}\,-\,|v_{2}^{\varepsilon}|^{2}\right)$ solves a forced spring-mass system with non constant frequency of oscillation of typical size $1/\varepsilon^{2}$. Since the foregoing three variables generates all expressions quadratic in ${\mathbf{V}}^{\varepsilon}$, this is indeed sufficient to conclude.

From previous considerations it follows that when $({\mathbf{X}}^{\varepsilon},{\mathbf{V}}^{\varepsilon})$ solves (1.2) the couple $({\mathbf{X}}^{\varepsilon},\frac{1}{2}\|{\mathbf{V}}^{\varepsilon}\|^{2})$ converges strongly to some $({\mathbf{Y}},g)$ as $\varepsilon\to 0$. To identify an evolution for the limiting $({\mathbf{Y}},g)$, we write relevant quadratic terms of (1.4) and (1.6) in terms of quantities of Lemma 1.2. First

[TABLE]

converges weakly to

[TABLE]

as $\varepsilon\to 0$. Thus, taking the limit $\varepsilon\to 0$ in (1.4) gives

[TABLE]

Likewise, ${\mathbf{V}}^{\varepsilon}\,\cdot\,\textrm{d}_{\mathbf{x}}{\mathbf{F}}(t,{\mathbf{X}}^{\varepsilon})\,({\mathbf{V}}^{\varepsilon})$ converges weakly to $g\,\textrm{div}_{\mathbf{x}}({\mathbf{F}})(t,{\mathbf{Y}})$ so that

[TABLE]

Gathering the latter results, we get that $({\mathbf{X}}^{\varepsilon},\frac{1}{2}\|{\mathbf{V}}^{\varepsilon}\|^{2})$ converges strongly to $({\mathbf{Y}},g)$ solving

[TABLE]

where $e^{0}\,=\,\frac{1}{2}\|{\mathbf{v}}^{0}\|^{2}$.

Remark 1.3**.**

Observe that $\textrm{\emph{div}}_{\mathbf{x}}({\mathbf{F}})=-{\mathbf{E}}^{\perp}\cdot\nabla_{\mathbf{x}}(1/b)+\textrm{\emph{curl}}({\mathbf{E}})/b$. In particular, as expected, the limiting system (1.12) contains the fact that when $b$ depends only on time and ${\mathbf{E}}$ derives from a potential, as in the Vlasov-Poisson case, the microscopic kinetic energy is conserved. Likewise, system (1.12) also implies that

[TABLE]

so that for the limiting system the classical adiabatic invariant $g/b$ becomes an exact invariant when the Maxwell-Faraday equation

[TABLE]

holds. The latter occurs in particular when $b$ depends only on the space variable and ${\mathbf{E}}$ derives from a potential, as in the Vlasov-Poisson case. For consistency’s sake note also that when $b$ is constant the microscopic kinetic energy and the adiabatic invariant essentially coincide and their conservation is equivalent to the irrotationality of ${\mathbf{E}}$.

1.2. Formal asymptotic limit of the Vlasov-Poisson system

We come back to the Vlasov-Poisson system (1.1). Here one cannot anymore remain completely at the characteristic level (1.2). Moreover whereas arguments of the previous subsection could be turned into sound analytic arguments, to the best of our knowledge the present situation does not fall directly into range of the actually available analysis of gyro-kinetic limits. We refer the reader to [13, 14, 20, 4, 19, 15] for a representative sample of such analytic techniques.

Nevertheless the previous subsection strongly suggests for $(f^{\varepsilon},{\mathbf{E}}^{\varepsilon})$ solving the Vlasov-Poisson system (1.1) that in the limit $\varepsilon\rightarrow 0$, the electric field ${\mathbf{E}}^{\varepsilon}$ and the following velocity-averaged version of $f^{\varepsilon}$

[TABLE]

converge to some ${\mathbf{E}}$ and some222We use distinct notation of variables for limiting functions to be consistent with asymptotic analysis at the characteristic level. This is of course completely immaterial. $f:(t,{\mathbf{y}},g)\mapsto f(t,{\mathbf{y}},g)$ solving the following system consisting in a transport equation supplemented with a Poisson equation,

[TABLE]

where the velocity field is given by

[TABLE]

with ${\mathbf{E}}=-\nabla_{\mathbf{y}}\phi$, ${\mathbf{F}}=-{\mathbf{E}}^{\perp}/b$.

Our goal is not to develop a thorough analysis of system (1.13) but let us mention that it generates a reasonable dynamics whose study falls into the scope of classical techniques.

Theorem 1.4**.**

*Consider $b\in W^{1,\infty}(\mathbb{R}^{+}\times\mathbb{R}^{2})$ satisfying (1.3).

Assume $f^{0}\in L^{1}\cap L^{\infty}(\mathbb{R}^{2}\times\mathbb{R}^{+})$, $f^{0}$ is nonnegative and*

[TABLE]

Then system (1.13) possesses a nonnegative weak solution $f$ starting from $f^{0}$ at time [math] such that both Lebesgue norms of $f$ and total energy are non-increasing in time, in particular,

[TABLE]

and

[TABLE]

Furthermore, when $f^{0}\in\mathcal{C}^{1}_{c}(\mathbb{R}^{2}\times\mathbb{R}^{+})$, the solution is unique and preserved along the characteristic curves

[TABLE]

where $({\mathbf{Y}},g)$ solves (1.12).

The foregoing result is not expected to be optimal in any reasonable way and its proof is completely analogous to those for the Vlasov-Poisson system ([1, 9] for weak solutions, [21] for smooth solutions). One key observation that underpins the analysis is that the vector field $({\mathbf{U}},u_{g})$ is divergence-free since

[TABLE]

Note also that for smooth solutions both Lebesgue norms and total energy are exactly preserved by the time evolution.

We now come to our main concern in the present article and seek after a numerical method that is able to capture these expected asymptotic properties, even when numerical discretization parameters are kept independent of $\varepsilon$ hence are not adapted to the stiffness degree of the fast scales. Our objective enters in the general framework of so-called Asymptotic Preserving (AP) schemes, first introduced and widely studied for dissipative systems as in [16, 17]. Yet, in opposition with collisional kinetic equations in hydrodynamic or diffusion asymptotics, collisionless equations like the Vlasov-Poisson system (1.1) involve fast time oscillations instead of fast time relaxation. By many respects this makes the identification of suitable schemes much more challenging.

In this context, among the few families of strategies available, most of them insist in providing also a good description of fast oscillations. To describe fast oscillations without imposing severe restrictions on time steps one needs in some sense to double time variables thus to see $f$ as the trace at $(t,\tau)=(t,t/\varepsilon^{2})$ of a function of variables $(t,\tau,{\mathbf{x}},{\mathbf{v}})$. The corresponding strategies have proved to be very successful when magnetic fields are uniform.

One of the oldest of these strategies, developed by E. Frénod, F. Salvarani and E. Sonnendrücker in [12], is directly inspired by theoretical results on two-scale convergence and relies on the fact that at the limit $\varepsilon\to 0$ the $\tau$-dependence may be explicitly filtered out. Its main drawback is probably that it computes only the leading order term in the limit $\varepsilon\to 0$. In particular it is only available when $\varepsilon$ is very small.

This may be fixed by keeping besides the stiff term to which a two-scale treatment is applied a non-stiff part that is smaller in the limit $\varepsilon\to 0$ but becomes important when $\varepsilon$ is not small. Such a decomposition may be obtained by using a micro-macro approach as in [6] and some references therein. This does allow to switch from one regime to another without any treatment of the transition between those but results in relatively heavy schemes,.

Another approach with similar advantages, developed in [7] and [11], consists in explicitly doubling time variables and seeking higher-dimensional partial differential equations and boundary conditions in variables $(t,\tau,{\mathbf{x}},{\mathbf{v}})$ that contains the original system at the $\varepsilon$-diagonal $(t,\tau)=(t,t/\varepsilon)$. While the corresponding methods are extremely good at capturing oscillations their design require a deep a priori understanding of the detailed structure of oscillations.

Since there are dramatic changes in the structure of oscillations when magnetic fields are not uniform, it seems that the extension of any method capturing oscillations to realistic magnetic fields is either doomed to fail or at least to require a tremendous amount of work and complexity. Observe in particular that in inhomogeneous cases whereas at the limit $\varepsilon\to 0$ the slow evolution still obeys a closed system independent of fast scales (as (1.13) above), fast oscillations do not uncouple anymore since their oscillating frequencies depend now on the slow scales… In any case this extension would constitute a major break-through in the field.

Our goal is somehow more modest as we only require that, in the limit $\varepsilon\to 0$, our schemes capture accurately the non stiff part of the evolution while allowing for coarse discretization parameters. Recently, we have indeed proposed a class of semi-implicit schemes which allow to capture the asymptotic limit of the two dimensional Vlasov-Poisson system with a uniform magnetic field [10]. We stress that by many respects those schemes are remarkably natural and simple.

Here, we show how our approach may be extended to some non uniform cases hence allowing to obtain direct simulations of system (1.1) with time steps large with respect to $\varepsilon$, while still capturing the slow coherent dynamics that emerge from strong oscillations. We develop numerical schemes that are able to deal with a wide range of values for $\varepsilon$ — that is, AP schemes in the terminology mentioned above. Our schemes are consistent with the kinetic model for any positive value of $\varepsilon$, and degenerate into schemes consistent with the asymptotic model (1.13) when $\varepsilon\rightarrow 0$.

We stress that the extension carried out here is by no means a trivial continuation of [10]. To understand inherent difficulties, we now review what are the main ingredients of the asymptotic analysis discussed in Section 1.1, those that we aim at incorporating at the discrete level. As a preliminary warning, it must be stressed however that once the time variable has been discretized the distinction between strong and weak convergences essentially disappear. This is immaterial to the treatment of homogeneous cases since then, as implicitly contained in Remark 1.1, there is only one key-ingredient, the weak convergence of ${\mathbf{V}}^{\varepsilon}/\varepsilon$ to ${\mathbf{F}}(t,{\mathbf{Y}})$, that may be safely replaced with a strong convergence at the discrete level. Here, beyond that requirement we also need three pieces of information essentially contained in Lemma 1.2 : weak convergences of $v_{1}^{\varepsilon}\,v_{2}^{\varepsilon}$ and $|v_{1}^{\varepsilon}|^{2}-|v_{2}^{\varepsilon}|^{2}$ to zero and strong convergence of $\frac{1}{2}(|v_{1}^{\varepsilon}|^{2}+|v_{2}^{\varepsilon}|^{2})$ to $g$ (solving (1.11)). Therefore we need to offer compatible counterparts at the discrete level of the weak convergence of $V^{\varepsilon}$ to zero and the strong convergence of $\|V^{\varepsilon}\|^{2}$ to a non trivial limit, a non trivial task !

2. A particle method for inhomogeneous strongly magnetized

plasmas

The schemes we shall propose here belong to the family of Particle-In-Cell (PIC) methods. Before presenting our specific schemes, we review in a few words the basic ideas underpinning these methods. We refer the reader to [2] for a thorough discussion and other applications to plasma physics.

To keep notation as light as possible, we temporarily omit to denote the dependence of solutions on $\varepsilon$ and discretization parameters. The starting point is the choice of approximating $f$, solving (1.1), by a finite sum of smoothed Dirac masses. More explicitly one would like to compute

[TABLE]

where $\varphi_{\alpha}=\alpha^{-d}\varphi(\cdot/\alpha)$ is a particle shape function with radius proportional to $\alpha$ — usually seen as a smooth approximation of the Dirac measure $\delta_{0}$ — obtained by scaling a fixed compactly supported mollifier $\varphi$, and the set $(({\mathbf{X}}_{k},{\mathbf{V}}_{k}))_{1\leq k\leq N}$ represents the position in phase space of $N$ macro-particles evolving along characteristic curves (1.2) from initial data $({\mathbf{x}}_{k}^{0},{\mathbf{v}}_{k}^{0})$, $1\leq k\leq N$. More explicitly

[TABLE]

where the electric field ${\mathbf{E}}$ solves the Poisson equation (1.1), or, more exactly in the end, it is computed from a solution of a discretization of the Poisson equation on a mesh of the physical space. Parameters are initially chosen to ensure that $f_{N}(0,{\mathbf{x}},{\mathbf{v}})$ is a good approximation of (continuous) initial datum $f^{0}$. To increase the order of approximation of Dirac masses, mollifiers are sometimes chosen with a prescribed number of vanishing moments. Concrete common choices include B-splines. See for instance [18, 5].

The main remaining issue is then to design an accurate and stable approximation of the particle trajectories in phase space. Note moreover that we ultimately want those properties to be uniform with respect to $\varepsilon$ at least for the slow part of the evolution. As already mentioned, to achieve this goal, we need to reproduce at the discrete level both the (weak) convergence of each ${\mathbf{V}}_{k}$ to zero and a non trivial slow dynamics of its microscopic energy $e_{k}=\tfrac{1}{2}\|{\mathbf{V}}_{k}\|^{2}$. We resolve this difficulty by augmenting the dimension of the phase space by one, looking for unknowns $({\mathbf{X}}_{k},{\mathbf{w}}_{k},e_{k})$ starting from $({\mathbf{x}}^{0}_{k},{\mathbf{v}}^{0}_{k},\tfrac{1}{2}\|{\mathbf{v}}_{k}^{0}\|^{2})$ and providing a solution to (2.1) by setting $({\mathbf{X}}_{k},{\mathbf{V}}_{k})=({\mathbf{X}}_{k},\sqrt{2e_{k}}\ {\mathbf{w}}_{k}/\|{\mathbf{w}}_{k}\|)$. The simplest possible choice of such systems would be

[TABLE]

with $({\mathbf{w}}_{k}^{0},e_{k}^{0})$ starting from $({\mathbf{v}}^{0}_{k},\tfrac{1}{2}\|{\mathbf{v}}_{k}^{0}\|^{2})$. However this choice is far from being uniquely determined. For instance one may add to the vector field defining System (2.2) any function of $(t,{\mathbf{X}},e,{\mathbf{w}})$ that vanishes when $e=\tfrac{1}{2}\|{\mathbf{w}}\|^{2}$. In the following we shall rather work with

[TABLE]

with $\chi$ chosen such that

[TABLE]

and

[TABLE]

For concreteness, in the following, we actually choose $\chi$ as

[TABLE]

where $s^{+}=\max(0,s)$. The choice of System 2.3 originates in the fact that it is relatively easy to discretize it so as to ensure that, in the limit $\varepsilon\to 0$, $\varepsilon^{-1}{\mathbf{w}}_{k}$ becomes asymptotically close to

[TABLE]

(as it occurs for $\varepsilon^{-1}{\mathbf{V}}_{k}$, see Section 1.1). Since we start with $e_{k}^{0}=\tfrac{1}{2}\|{\mathbf{w}}_{k}^{0}\|^{2}$ and the constraint is preserved by the evolution, the choice of either System (2.2) or of (2.3) or of any similar system is perfectly immaterial as long as we do not discretize the time variable. However their discretizations will no longer preserve exactly the condition $e_{k}^{0}=\tfrac{1}{2}\|{\mathbf{w}}_{k}^{0}\|^{2}$ thus they will lead to effectively distinct solutions.

Note that the addition of a new phase space variable is much less invasive than doubling time variables (but is insufficient to allow for the description of fast oscillations with coarse grids) and it is also very cheap in the context of PIC methods. In particular, consistently with the former claim, now that once we have chosen an augmented system such as (2.3) we may discretize it by adapting almost readily the strategy developed in [3, 10] and based on semi-implicit solvers for stiff problems.

In the rest of this section, we shall focus on the discretization of System (2.3) and propose several numerical schemes for that purpose. Yet, prior to that, we briefly summarize the above discussion and explain how this discretization will be used to compute solutions to (1.1). Besides already introduced parameters, fix a given time step $\Delta t>0$ and for $n\geq 0$ introduce discrete time $t^{n}=n\,\Delta t$. Then we chose a discretization of System (2.3) expected to provide for any $k\in\{1,\ldots,N\}$, $({\mathbf{x}}^{n}_{k},,e^{n}_{k},{\mathbf{w}}^{n}_{k})$ an approximation of $({\mathbf{X}}_{k}(t^{n}),e_{k}(t^{n}),{\mathbf{w}}_{k}(t^{n}))$, where $({\mathbf{X}}_{k},e_{k},{\mathbf{w}}_{k})$ solves (2.3). Of course equations are not uncoupled thus we compute simultaneously all $({\mathbf{x}}^{n}_{k},,e^{n}_{k},{\mathbf{w}}^{n}_{k})$, $k\in\{1,\ldots,N\}$, and the corresponding electric field computed from a solution on a mesh of the physical space of a discretization of the Poisson equation with macroscopic density corresponding to the distribution function $(f_{N,\alpha}^{n})_{n\geq 0}$, given at time $t^{n}$ by

[TABLE]

where

[TABLE]

2.1. A first-order semi-implicit scheme

Since we only discuss possible discretizations of (2.3), we shall omit from now on the index $k\in\{1,\ldots,N\}$ and consider the electric field as fixed, assuming smoothness of $({\mathbf{E}},b)$ and (1.3).

We start with the simplest semi-implicit scheme for (2.3), which is a combination of the backward and forward Euler schemes. For a fixed time step $\Delta t>0$ it is given by

[TABLE]

Notice that only the third equation on ${\mathbf{w}}^{n+1}$ is really implicit and it only requires the resolution of a two-dimensional linear system. Then, once the value of ${\mathbf{w}}^{n+1}$ has been computed the first and second equations provide explicitly the values of ${\mathbf{x}}^{n+1}$ and $e^{n+1}$.

Proposition 2.1** (Consistency in the limit $\varepsilon\rightarrow 0$ for a fixed $\Delta t$).**

Let us consider a time step $\Delta t>0$, a final time $T>0$ and set $N_{T}=\lfloor T/\Delta t\rfloor$. Assume that $({\mathbf{x}}^{n}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{n},e^{n}_{\varepsilon})_{0\leq n\leq N_{T}}$ is a sequence obtained by (2.6) and such that

• •

for all $1\leq n\leq N_{T}$, $\left({\mathbf{x}}^{n}_{\varepsilon},\varepsilon{\mathbf{w}}^{n}_{\varepsilon},e^{n}_{\varepsilon}\right)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon>0$;

• •

$\left({\mathbf{x}}^{0}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{0},e^{0}_{\varepsilon}\right)_{\varepsilon>0}$* converges in the limit $\varepsilon\rightarrow 0$ to some $({\mathbf{y}}^{0},{\mathbf{w}}^{0},g^{0})$.*

Then, for any $1\leq n\leq N_{T}$, $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{\varepsilon>0}$ converges to some $({\mathbf{y}}^{n},g^{n})$ as $\varepsilon\rightarrow 0$ and the limiting sequence $({\mathbf{y}}^{n},g^{n})_{1\leq n\leq N_{T}}$ solves

[TABLE]

which provides a consistent first-order approximation with respect to $\Delta t$ of the gyro-kinetic system (1.12).

Note that we actually apply this result with $\left({\mathbf{x}}^{0}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{0},e^{0}_{\varepsilon}\right)_{\varepsilon>0}=({\mathbf{x}}^{0},{\mathbf{v}}^{0},\tfrac{1}{2}\|{\mathbf{v}}^{0}\|^{2})$ so that in this case $({\mathbf{y}}^{0},{\mathbf{w}}^{0},g^{0})=({\mathbf{x}}^{0},{\mathbf{v}}^{0},\tfrac{1}{2}\|{\mathbf{v}}^{0}\|^{2})$. Furthermore, in this case, the limiting scheme (2.7) starts from initial data

[TABLE]

More generally, under the assumptions of Proposition 2.1, $({\mathbf{y}}^{1},g^{1})$ differs from $({\mathbf{y}}^{0},g^{0})$ by a first-order term but in general does not satisfy (2.7) with $n=0$. The scheme needs one time-step to reach the expected asymptotic regime in the limit $\varepsilon\to 0$.

Proof.

First, note that, for any $0\leq n\leq N_{T}$, since $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{\varepsilon>0}$ is bounded, ${\mathbf{E}}(t^{n},{\mathbf{x}}^{n}_{\varepsilon})$, $\nabla_{\mathbf{x}}b(t^{n},{\mathbf{x}}^{n}_{\varepsilon})$ and $\chi(e^{n}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{n})$ are also bounded uniformly in $\varepsilon$. Therefore, combined with assumption (1.3) and boundedness of $\left(\varepsilon{\mathbf{w}}^{n}_{\varepsilon}\right)_{\varepsilon>0}$, the third equation of (2.6) written as

[TABLE]

shows that, for any $1\leq n\leq N_{T}$, $\left(\varepsilon^{-1}{\mathbf{w}}^{n}_{\varepsilon}\right)_{\varepsilon>0}$ is bounded and in particular $\left({\mathbf{w}}^{n}_{\varepsilon}\right)_{\varepsilon>0}$ converges to zero. Moreover, for any $1\leq n\leq N_{T}$, since the sequence $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon>0$, we can extract a subsequence still abusively labeled by $\varepsilon$ and find some $({\mathbf{y}}^{n},g^{n})$ such that $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})\rightarrow({\mathbf{y}}^{n},g^{n})$ as $\varepsilon$ goes to zero.

Therefore, by using continuity of $(E,b,\nabla_{\mathbf{x}}b,\chi)$, we first conclude that for any $1\leq n\leq N_{T}-1$, (2.8) yields

[TABLE]

since $\chi(g^{n},0)=g^{n}$. By substituting the foregoing limit in the first and second equations of (2.6) we obtain that the limit $({\mathbf{y}}^{n},g^{n})_{1\leq n\leq N_{T}}$ satisfies (2.7). Moreover the same argument shows that

[TABLE]

Since the limit points $({\mathbf{y}}^{n},g^{n})$ of the extracted subsequence are uniquely determined (by (2.7) and (2.9)), actually all the sequence $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{\varepsilon>0}$ converges (for any $n$). ∎

Remark 2.2**.**

The consistency provided by the latter result is far from being uniform with respect to the time step. However, though we restrain from doing so here in order to keep technicalities to a bare minimum, we expect that an analysis similar to the one carried out in [10, Section 4] could lead to uniform estimates, proving uniform stability and consistency with respect to both $\Delta t$ and $\varepsilon$.

Of course, a first-order scheme may fail to be accurate enough to describe correctly the long time behavior of the solution, but it has the advantage of simplicity. In the following we show how to generalize our approach to second and third-order schemes.

2.2. Second-order semi-implicit Runge-Kutta schemes

Now, we come to second-order schemes with two stages. The scheme we consider is a combination of a Runge-Kutta method for the explicit part and of an $L$-stable second-order SDIRK method for the implicit part.

Remark 2.3**.**

In our original treatment of the homogeneous case in [10] where we design schemes to capture (1.5), the discrete velocity ${\mathbf{v}}^{n}_{\varepsilon}$ was damped to zero in the limit $\varepsilon\to 0$, consistently with weak convergence of ${\mathbf{V}}^{\varepsilon}$ to zero but not with conservation of $\tfrac{1}{2}\|{\mathbf{V}}^{\varepsilon}\|^{2}$. Therefore there was some interest there in tailoring schemes not too dissipative. Here at the discrete level the two distinct behaviors are encoded by two distinct variables ${\mathbf{w}}^{n}_{\varepsilon}$ and $e^{n}_{\varepsilon}$ so that we may rightfully focus solely on implicit $L$-stable choices.

To describe the scheme, we introduce $\gamma>0$ the smallest root of the polynomial $X^{2}-2X+1/2$, i.e. $\gamma=1-1/\sqrt{2}$. Then the scheme is given by the following two stages. First,

[TABLE]

with

[TABLE]

Before the second stage, we first introduce $\hat{t}^{(1)}\,=\,t^{n}+{\Delta t}/{(2\gamma)}$ and explicitly compute $(\hat{{\mathbf{x}}}^{(1)},\hat{{\mathbf{w}}}^{(1)},\hat{e}^{(1)})$ from

[TABLE]

Then the solution of the second stage $({\mathbf{x}}^{(2)},{\mathbf{w}}^{(2)},e^{(2)})$ is given by

[TABLE]

with

[TABLE]

Finally, the numerical solution at the next time step is defined by

[TABLE]

Proposition 2.4** **(Consistency in the

limit $\varepsilon\rightarrow 0$ for a fixed $\Delta t$).

Let us consider a time step $\Delta t>0$, a final time $T>0$ and set $N_{T}=\lfloor T/\Delta t\rfloor$. Assume that $({\mathbf{x}}^{n}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{n},e^{n}_{\varepsilon})_{0\leq n\leq N_{T}}$ is a sequence obtained by (2.10)-(2.15) and such that

• •

for all $1\leq n\leq N_{T}$, $\left({\mathbf{x}}^{n}_{\varepsilon},\varepsilon{\mathbf{w}}^{n}_{\varepsilon},e^{n}_{\varepsilon}\right)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon>0$;

• •

$\left({\mathbf{x}}^{0}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{0},e^{0}_{\varepsilon}\right)_{\varepsilon>0}$* converges in the limit $\varepsilon\rightarrow 0$ to some $({\mathbf{y}}^{0},{\mathbf{w}}^{0},g^{0})$.*

Then, for any $1\leq n\leq N_{T}$, $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{\varepsilon>0}$ converges to some $({\mathbf{y}}^{n},g^{n})$ as $\varepsilon\rightarrow 0$ and the limiting sequence $({\mathbf{y}}^{n},g^{n})_{0\leq n\leq N_{T}}$ solves

[TABLE]

where

[TABLE]

and

[TABLE]

with $\hat{t}^{(1)}=t^{n}+\Delta t/(2\gamma)$ and

[TABLE]

which provides a consistent second-order approximation of the gyro-kinetic system (1.12).

Proof.

We mainly follow the lines of the proof of Proposition 2.1. To make arguments more precise we mark with a suffix n+1,ε intermediate quantities involved in the step from $t^{n}$ to $t^{n+1}$.

To begin with, for $1\leq n\leq N_{T}$, the third equation of (2.10) implies that $(\varepsilon^{-1}{\mathbf{w}}^{(1)}_{n,\varepsilon})_{\varepsilon}$ is bounded, then the first and second equations of (2.12) show that $(\hat{{\mathbf{x}}}^{(1)}_{n,\varepsilon},\hat{e}^{(1)}_{n,\varepsilon})_{\varepsilon}$ is also bounded, from which the third equation of (2.10) shows the boundedness of $(\varepsilon^{-1}{\mathbf{w}}^{(1)}_{n,\varepsilon})_{\varepsilon}$. Then from the third equation of (2.12) we derive that $(\hat{{\mathbf{w}}}^{1}_{1,\varepsilon})_{\varepsilon}$ is bounded and that, for $2\leq n\leq N_{T}$, $(\varepsilon^{-1}\hat{{\mathbf{w}}}^{(2)}_{n,\varepsilon})_{\varepsilon}$ is bounded.

One may then extract converging subsequences from all bounded families. With obvious notation for limits this leads in particular along the chosen subsequence from the third equations of (2.10) and (2.13) to

[TABLE]

and

[TABLE]

when $1\leq n\leq N_{T}-1$, since ${\mathbf{w}}^{n}_{\varepsilon}\to 0$ and ${\mathbf{w}}^{(1)}_{n+1,\varepsilon}\to 0$. With this in hands taking, along the same subsequence, the limit $\varepsilon\to 0$ of first and second equations of (2.12) and (2.13) provides the claimed scheme.

Moreover the same arguments provide explicit expressions of $({\mathbf{y}}^{1},g^{1})$ in terms of $({\mathbf{y}}^{0},{\mathbf{w}}^{0},g^{0},\Delta t)$. This yields uniqueness of limits independently of the chosen subsequence hence convergence of full families. ∎

2.3. Third-order semi-implicit Runge-Kutta schemes

Now we consider a third-order semi-implicit scheme. It is given by a four stages Runge-Kutta method and was introduced in the framework of hyperbolic systems with stiff source terms in [3].

First, we choose $\alpha=0.24169426078821$, $\eta=0.12915286960590$ and set $\beta=\alpha/4$, $\gamma=1/2-\alpha-\beta-\eta$. Then the first stage of the scheme is

[TABLE]

with

[TABLE]

The second stage is

[TABLE]

with

[TABLE]

Then, for the third stage we first compute explicitly intermediate values

[TABLE]

and then construct the new stage $\left({\mathbf{x}}^{(3)},e^{(3)},{\mathbf{w}}^{(3)}\right)$

[TABLE]

with

[TABLE]

Finally, we use explicit intermediate values $\left(\hat{{\mathbf{x}}}^{(3)},\hat{e}^{(3)},\hat{{\mathbf{w}}}^{(3)}\right)$

[TABLE]

to carry out the fourth stage

[TABLE]

with

[TABLE]

where $t^{n+1/2}=t^{n}+\Delta t/2$. The numerical solution at the new time step is finally given by

[TABLE]

We emphasize that in all our schemes at each stage the implicit computation only requires the resolution of a two-by-two linear system. Therefore the computational effort is of the same order as fully explicit schemes.

Proposition 2.5** **(Consistency in the

limit $\varepsilon\rightarrow 0$ for a fixed $\Delta t$).

Let us consider a time step $\Delta t>0$, a final time $T>0$ and set $N_{T}=\lfloor T/\Delta t\rfloor$. Assume that $({\mathbf{x}}^{n}_{\varepsilon},{\mathbf{w}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{0\leq n\leq N_{T}}$ is a sequence obtained by (2.17)-(2.23) and such that

• •

for all $1\leq n\leq N_{T}$, $\left({\mathbf{x}}^{n}_{\varepsilon},\varepsilon{\mathbf{w}}^{n}_{\varepsilon},e^{n}_{\varepsilon}\right)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon>0$;

• •

$\left({\mathbf{x}}^{0}_{\varepsilon},{\mathbf{w}}_{\varepsilon}^{0},e^{0}_{\varepsilon}\right)_{\varepsilon>0}$* converges in the limit $\varepsilon\rightarrow 0$ to some $({\mathbf{y}}^{0},{\mathbf{w}}^{0},g^{0})$.*

Then, for any $1\leq n\leq N_{T}$, $({\mathbf{x}}^{n}_{\varepsilon},e^{n}_{\varepsilon})_{\varepsilon>0}$ converges to some $({\mathbf{y}}^{n},g^{n})$ as $\varepsilon\rightarrow 0$ and the limiting sequence $({\mathbf{y}}^{n},g^{n})_{0\leq n\leq N_{T}}$ solves

[TABLE]

where

[TABLE]

whereas ${\mathbf{U}}^{(1)}$ and $T^{(1)}$ are given by

[TABLE]

and ${\mathbf{U}}^{(2)}$ and $T^{(2)}$

[TABLE]

with $t^{n+1/2}=t^{n}+\Delta t/2$,

[TABLE]

and

[TABLE]

which provides a consistent third-order approximation of the gyro-kinetic system (1.12).

We omit the proof of Proposition 2.5 as almost identical to the one of Proposition 2.4.

Remark 2.6**.**

For both our second and third-order schemes, the limit obtained when $\varepsilon\to 0$ fails to satisfy the first step of the limiting schemes by a first-order error. This is obviously a more serious problem that in the first-order case. However, though we do not pursue this line of investigation here, one should be able to fix this issue by replacing the very first step (from $t^{0}$ to $t^{1}$) of the original ($\varepsilon$-dependent) computation by a step where third equations of intermediate values are also implicited.

3. Numerical simulations

In this section, we provide examples of numerical computations to validate and compare the different time discretization schemes introduced in the previous section.

We first consider the motion of a single particle under the effect of a given electromagnetic field. It allows us to illustrate the ability of the semi-implicit schemes to capture in the limit $\varepsilon\rightarrow 0$ drift velocities due to variations of magnetic and electric fields, even with large time steps $\Delta t$ .

Then we consider the Vlasov-Poisson system with an external non uniform magnetic field. As already mentioned we implement a classical particle-in-cell method but with specific different time discretization techniques to compute the particle trajectories. In particular, in this case, the collection of charged particles moves collectively and gives rise to a self-consistent electric field, obtained by solving numerically the Poisson equation in a space grid.

3.1. One single particle motion without electric field

Before simulating at the statistical level, we investigate on the motion of individual particles in a given magnetic field the accuracy and stability properties with respect to $\varepsilon>0$ of the semi-implicit algorithms presented in Section 2.

Numerical experiments of the present subsection are run with a zero electric field ${\mathbf{E}}=0$, and a time-independent external magnetic field corresponding to

[TABLE]

with $\alpha=0.5$. Moreover we choose for all simulations the initial data as ${\mathbf{x}}^{0}=(5,4)$, ${\mathbf{v}}^{0}=(5,6)$, whereas the final time is $T=2$. In this case, the asymptotic drift velocity predicted by the limiting model (1.12) is explicitly given by

[TABLE]

First, for comparison, we compute a reference solution $({\mathbf{X}}^{\varepsilon},{\mathbf{w}}^{\varepsilon},e^{\varepsilon})_{\varepsilon>0}$ to the initial problem (2.2) thanks to an explicit fourth-order Runge-Kutta scheme used with a very small time step of the order of $\varepsilon^{2}$ and a reference solution $({\mathbf{Y}},g)$ to the (non stiff) asymptotic model (1.12) obtained when $\varepsilon\rightarrow 0$. Recall that the derivation of (1.12) also shows weak convergence of $\varepsilon^{-1}{\mathbf{w}}^{\varepsilon}$ to ${\mathbf{U}}^{0}={\mathbf{U}}({\mathbf{Y}},g)$ in the limit $\varepsilon\to 0$. Then we compute an approximate solution $({\mathbf{X}}^{\varepsilon}_{\Delta t},{\mathbf{w}}^{\varepsilon}_{\Delta t},e^{\varepsilon}_{\Delta t})$ using either (2.10)–(2.15) or (2.17)–(2.23), and compare them to the reference solutions.

Our goal is to evaluate the accuracy of the numerical solution $({\mathbf{X}}^{\varepsilon}_{\Delta t},{\mathbf{w}}^{\varepsilon}_{\Delta t},e^{\varepsilon}_{\Delta t})$ for various regimes when both $\varepsilon$ and $\Delta t$ vary. Computed errors are measured in discrete $L^{1}$ norms. For instance, we set333When necessary and not too confusing, as here, we use the same piece of notation for continuous time-dependent functions and their discrete counterpart, obtained by restriction to discrete times.

[TABLE]

In Figure 1, we present the numerical error on space and velocity444In graphics we use ${\mathbf{v}}$ as a short-hand for $\varepsilon^{-1}{\mathbf{w}}$. This should not be confused with the actual reconstruction of the velocity variable used at the statistical level. variables between the reference solution for the initial problem (1.2) and the one obtained with the third-order scheme (2.17)-(2.23). As expected for a fixed time step $\Delta t$ (taken here between $0.0025$ and $0.01$), the scheme is quite stable even in the limit $\varepsilon\rightarrow 0$ and the error on the space variable, measured by $\|{\mathbf{X}}^{\varepsilon}_{\Delta t}-{\mathbf{X}}^{\varepsilon}\|$, is uniformly small and reach a maximum on intermediate regimes, $\varepsilon\in(0.05,\;0.5)$. In contrast, for a fixed time step, the error on the velocity variable, measured by $\varepsilon^{-1}\,\|{\mathbf{w}}^{\varepsilon}_{\Delta t}-{\mathbf{w}}^{\varepsilon}\|$, is small for large values of $\varepsilon$, but gets very large when $\varepsilon\rightarrow 0$ since the scheme cannot follow high-frequency time oscillations of order $\varepsilon^{-2}$ when $\varepsilon\ll\sqrt{\Delta t}$.

In Figure 2, for the same set of simulations we plot errors with respect to the limiting system (1.12). This shows, see Figure 2 (b), that for a fixed $\Delta t$ we do capture however, in the limit $\varepsilon\to 0$, correct drift velocities. More generally, as expected, errors on both space and velocity variables with respect to limiting asymptotic values gets smaller and smaller as $\varepsilon\to 0$.

Finally, in Figure 3 we show similar quantities for the second-order scheme (2.10)-(2.15) holding the time step fixed to $\Delta t=0.01$. For comparison, we also plot results of the third-order scheme. Errors of the third-order scheme are smaller but behaviors are qualitatively the same : uniform precision on the slow variable, convergence to asymptotic descriptions as $\varepsilon\to 0$ and consistency of the velocity variable with the velocity of the true solution when $\varepsilon\gg\sqrt{\Delta t}$ and with the asymptotic drift when $\varepsilon\ll\sqrt{\Delta t}$.

3.2. One single particle motion with an electromagnetic field

We now add to the otherwise unmodified previous setting a non trivial electric field

[TABLE]

In this case, at the limit $\varepsilon\to 0$ the behavior of the microscopic kinetic energy remains non trivial (in contrast with what happens in the previous subsection) so that we also measure errors on its evaluations. Moreover the electric field also contributes to asymptotic drift velocities, that are now given by

[TABLE]

with $\alpha=1/2$. This simple geometric configuration allows us to clearly distinguish the contribution to the drift velocity due to the interaction of the electric field with the magnetic field from purely magnetic effects due to the gradient of the magnetic field.

Quantitative error measurements are provided in Figure 4. Qualitative features are completely analogous to those of the previous subsection, the microscopic kinetic energy being approximated with essentially the same accuracy as the other slow variable ${\mathbf{x}}$.

To visualize rather than quantify what happens in the limit $\varepsilon\to 0$, we also represent now space trajectories corresponding to different values of $\varepsilon$, but holding the time step fixed to $\Delta t=0.01$. On the reference solution one sees as $\varepsilon\to 0$ faster and faster oscillations — of order $\varepsilon^{-2}$ — but with smaller and smaller amplitude — of order $\varepsilon$ — around guiding center trajectories as predicted by555Note incidentally that this illustrates a classical abuse in terminology, used throughout the present paper. Variables $({\mathbf{X}}^{\varepsilon},e^{\varepsilon})$ are not really slow but they are uniformly close to genuinely slow variables. (1.12). When $\varepsilon\gg 0.1$ the approximate solution follows closely the reference solution, whereas when $\varepsilon\ll 0.1$ they are also almost indistinguishable though the approximate solution does not reproduce the infinitely fast oscillations of the true solution. Main discrepancies are observed for intermediate values, here exemplified by the case $\varepsilon=0.2$, when the true solution is already too fast for the approximate solution but not yet close enough to the guiding center dynamics.

As a conclusion, these elementary numerical simulations confirm the ability of the semi-implicit discretization to capture the evolution of the ”slow” variables $({\mathbf{X}}^{\varepsilon},e^{\varepsilon})$ uniformly with respect to $\varepsilon$ by essentially transitioning automatically to the the guiding center motion encoded by (1.12) when it is not enable anymore to follow the fast oscillations of the true evolution.

3.3. The Vlasov-Poisson system

We now consider the Vlasov-Poisson system (1.1) set in a disk $\Omega=D(0,6)$ centered at the origin and of radius $6$. Our simulations start with an initial data that is Maxwellian in velocity and whose macroscopic density is the sum of two Gaussians, explicitly

[TABLE]

with ${\mathbf{x}}_{0}=(3/2,-3/2)$. The parameter $\varepsilon$ is set either to $\varepsilon=1$, leading to a non-stiff problem or to $\varepsilon=0.05$, where the asymptotic regime is relevant. For both regimes we compute numerical solutions to the Vlasov-Poisson system (1.1) with the third-order scheme (2.17)-(2.23) and time step $\Delta t=0.1$. Note that heuristic arguments based on single-particle simulations suggest that when $\varepsilon=0.05$ this enforces a regime where we are close to the asymptotic limit and our schemes cannot capture fast oscillations of the true solution.

We run two sets of numerical simulations with a time-independent inhomogeneous magnetic field

[TABLE]

that is radial increasing with value one at the origin.

In Figure 6 we represent time evolutions of the total energy

[TABLE]

and of the adiabatic invariant

[TABLE]

As far as smooth solutions are concerned, the total energy is preserved by both the original $\varepsilon$-dependent model and by the asymptotic model (1.13). One goal of our experimental observations is to check that despite the fact that our scheme dissipates some parts of the velocity variable to reach the asymptotic regime corresponding to (1.13) it does respect this conservation. In contrast an essentially exact conservation of the adiabatic invariant is a sign that we have reached the limiting asymptotic regime since it does not hold for the original model but does for the asymptotic (1.13) as $b$ is time-independent and ${\mathbf{E}}$ is curl-free. Observe that, since $b$ is not homogeneous, even in the asymptotic regime the kinetic and potential parts of the total energy are not preserved separately, but the total energy corresponding to the Vlasov-Poisson system is still preserved. Figure 6 shows that all these features are captured satisfactorily by our scheme even on long time evolutions with a large time step and despite its dissipative implicit nature when $\varepsilon=1$ (Figure 6 left column) and also in the asymptotic regime when $\varepsilon=0.05$ (Figure 6 right column).

In Figure 7 and 8 we visualize the corresponding dynamics by presenting some snapshots of the time evolution of the macroscopic charge density. In the first case, with $\varepsilon=1$, the magnetic field is not sufficiently large to provide a good confinement of the macroscopic density, and some macroscopic oscillations occur as already observed on the time evolution of macroscopic quantities. On the other hand, our inhomogeneous case does exhibit some confinement when $\varepsilon\to 0$, since the conservation of $e/b({\mathbf{x}})$ offers (a weak form of) coercivity jointly in $({\mathbf{x}},e)$. Such a confinement is indeed observed when $\varepsilon=0.05$, jointly with the expected eventual merging of two initial vortices in a relatively short time.

4. Conclusion and perspectives

In the present paper we have proposed a class of semi-implicit time discretization techniques for particle-in cell simulations. The main feature of our approach is to guarantee the accuracy and stability on slow scale variables even when the amplitude of the magnetic field becomes large, thus allowing a capture of their correct long-time behavior including cases with non homogeneous magnetic fields and coarse time grids. Even on large time simulations the obtained numerical schemes also provide an acceptable accuracy on physical invariants (total energy for any $\varepsilon$, adiabatic invariant when $\varepsilon\ll 1$) whereas fast scales are automatically filtered when the time step is large compared to $\varepsilon^{2}$.

As a theoretical validation we have proved that under some stability assumptions on numerical approximations, the slow part of the approximation converges when $\varepsilon\rightarrow 0$ to the solution of a limiting scheme for the asymptotic evolution, that preserves the initial order of accuracy. Yet a full proof of uniform accuracy and a classification of admissible schemes — going beyond the exposition of a few possible schemes — remain to be carried out.

From a practical point of view, the next natural step would be to consider the genuinely three-dimensional Vlasov-Poisson system. We point out however that, despite recent progresses [4, 8], the asymptotic limit $\varepsilon\to 0$ is much less understood in this context, even at the continuous level with external electro-magnetic fields.

Acknowledgements

FF was supported by the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

LMR was supported in part by the ANR project BoND (ANR-13-BS01-0009-01). LMR also expresses his appreciation of the hospitality of IMT, Université Toulouse III, during part of the preparation of the present contribution.