# Uniformly accurate methods for three dimensional Vlasov equations under strong magnetic field with varying direction

In this paper, we consider the three dimensional Vlasov equation with an inhomogeneous, varying direction, strong magnetic field. Whenever the magnetic field has constant intensity, the oscillations generated by the stiff term are periodic. The homogenized model is then derived and several state-of-the-art multiscale methods, in combination with the Particle-In-Cell discretisation, are proposed for solving the Vlasov-Poisson equation. Their accuracy as much as their computational cost remain essentially independent of the strength of the magnetic field. The proposed schemes thus allow large computational steps, while the full gyro-motion can be restored by a linear interpolation in time. In the linear case, extensions are introduced for general magnetic field (varying intensity and direction). Eventually, numerical experiments are exposed to illustrate the efficiency of the methods and some long-term simulations are presented.

## Authors

• 6 publications
• 3 publications
• 3 publications
• 3 publications
• 10 publications
05/22/2020

### Error estimates of some splitting schemes for charged-particle dynamics under strong magnetic field

In this work, we consider the error estimates of some splitting schemes ...
02/22/2018

### High order time integrators for the simulation of charged particle motion in magnetic quadrupoles

Magnetic quadrupoles are essential components of particle accelerators l...
09/06/2019

### Magnetically actuated artificial microswimmers as mobile microparticle manipulators

Micro-scale swimming robots have been envisaged for many medical applica...
05/17/2022

### Semi-discretization and full-discretization with optimal accuracy for charged-particle dynamics in a strong nonuniform magnetic field

The aim of this paper is to formulate and analyze numerical discretizati...
10/24/2020

### Towards Real-Time Magnetic Dosimetry Simulations for Inductive Charging Systems

The exposure of a human by magneto-quasistatic fields from wireless char...
04/29/2022

### Long term analysis of splitting methods for charged-particle dynamics

In this paper, we rigorously analyze the energy, momentum and magnetic m...
11/05/2019

### Numerical methods for antiferromagnetics

Compared with ferromagnetic counterparts, antiferromagnetic materials ar...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Vlasov models have been widely considered for modelling the dynamics of plasmas as encountered in magnetic fusion devices known as a tokamaks, where a strong external magnetic field is applied so as to confine the charged particles. In this paper, we consider the three dimensional Vlasov-Poisson equation with a strong non-homogeneous magnetic field whose direction may vary [16, 24, 35]

 ∂tfε(t,x,v)+v⋅∇xfε(t,x,v)+(E(t,x)+1εv×B(x))⋅∇vfε(t,x,v)=0, (1.1a) ∇x⋅E(t,x)=∫R3fε(t,x,v)dv−ni, (1.1b) fε(0,x,v)=f0(x,v), (1.1c)

where, for a given ,

 fε:(t,x,v)∈[0,T]×R3×R3↦fε(t,x,v)∈R

is the unknown,

 f0:(x,v)∈R3×R3↦f0(x,v)∈R

a given initial distribution, where

 B:x∈R3↦B(x)∈R3

denotes the external magnetic field,

 E:(t,x)∈R+×R3↦E(t,x)∈R3

the self-consistent electric-field function, a dimensionless parameter inversely proportional to the strength of the magnetic field and the ion density of the background. The system (1.1) has a lot of invariants and we will be interested in particular in the Hamiltonian defined by

 H(t):=∫R3∫R312|v|2fε(t,x,v)dxdv+12∫R3|E(t,x)|2dx. (1.2)

The above Vlasov-Poisson model (1.1) is derived from the three dimensional Vlasov-Maxwell equations by considering the electrostatic approximation. Unlike some asymptotically reduced models such as the gyrokinetic equations [29, 33] or the drift-kinetic limit equations [3, 22, 18], model (1.1) is set at the kinetic scale and is of paramount importance for studying the plasma dynamics in the tokamak device.

In the strong magnetic field limit regime, the charged particles exhibit very fast rotations with cyclotron period proportional to , while remaining confined along the magnetic line. In such a case, the small parameter renders the solution of (1.1) highly-oscillatory in time. Classical numerical integrators such as splitting or finite-difference schemes thus require time steps smaller than the cyclotron period in order to accurately capture the dynamics, thus implying severe computational burden. Recent efforts have aimed at designing numerical schemes which allow step-sizes much larger than the cyclotron period. Upon assuming that the magnetic field has a fixed direction in space, i.e. that

 B(x)=(0,0,b(x))T,b(x)>0

(a popular choice both for formal and rigorous analyses [3, 16, 24]), several multiscale numerical methods have been proposed [10, 15, 17, 19, 20, 23]. Among them, Filbet et al. constructed Particle-In-Cell schemes in the spirit of asymptotic preserving techniques [27], which, as , are consistent with the drift-limit model [19] or the gyrokinetic model [17, 20]. These schemes are simple and highly accurate, but the gyro-motion is lost in the limit regime. In contrast, the schemes proposed in [10, 15] capture all the information of the kinetic models with an accuracy uniform with respect to . These uniformly accurate (UA) schemes have computational cost as well as accuracy totally independent of (we refer to [14] for a comparison of UA scheme with other multiscale methods). In order to design UA schemes for kinetic models, different numerical approaches may be used: (i) The two-scale formulation technique relies upon an explicit separation of the fast and slow times and allows to smooth out the oscillations [10, 13, 15]. (ii) The multi-revolution composition methods, in the spirit of heterogeneous multiscale method [1], are also UA, as confirmed in the recent paper [11]. Both approaches exploit the periodicity of the solution of the stiff part of the equation. For instance, our recent work [10] isolates the dominant oscillation frequency owing to a confining property in two dimensions. However, the general case of a strong magnetic field with varying direction, i.e.

 B(x)=(b1(x),b2(x),b3(x))T,

has been barely considered so far for the Vlasov-Poisson equation (1.1) due to its complicated highly-oscillatory behaviour in three dimensions. Let us also mention recent developments around symplectic Particle-In-Cell method, which allows for good preservation of invariants for very long time (see [25, 32, 28]).

In this work, we propose efficient numerical schemes for solving the three dimensional Vlasov-Poisson equation (1.1) in the strong magnetic field regime by combining multiscale strategies with the Particle-In-Cell (PIC) discretisation. First, we consider the case of a magnetic field with constant intensity , for which, as already pointed out in [5, 33], the motion induced by the stiff Lorentz term in (1.1) is periodic in time. Taking advantage of this observation, we derive the limit model of (1.1) by using averaging methods [5], and then introduce three UA schemes, namely (i) the multi-revolution composition (MRC) method, (ii) the two-scale formulation (TSF) method and (iii) the micro-macro (MM) method. All three are of uniform second order in time for all , though have specific pros and cons: for instance, MRC methods are phase-space volume preserving, while MM easily allows for the full recovery of the gyro-motion. To the best of our knowledge, this key-feature is new and paves the way for an extension to the case of a magnetic field with varying intensity. In this situation, we indeed introduce, under the PIC discretisation, a reparametrization of time to re-normalise the magnetic field. Within this framework, each particle carries its own fictitious time. Hence, and in order to avoid the occurrence of multiple frequencies, we drop in this situation the Poisson part of (1.1) and consider instead the case of an external electric field (this somehow simplifying assumption is relevant as it marks an important first step towards the solution of the full problem). In order to re-synchronise all particles (a necessary step in order to provide an approximation of ), we then use the interpolation strategy of MM which ensures uniform second order except for the angular variable. Eventually, numerical experiments are presented in order to validate uniform accuracy and to compare the various methods. In particular, we simulate the dynamics of equation (1.1) in a three dimensional screw-pinch setup [29].

The remaining of the paper is now organized as follows. Section 2 considers the limit model of (1.1) with a constant intensity and Section 3 introduces the three aforementioned UA schemes in this situation: Subsection 3.1 is concerned with MRC method, Subsection 3.2 with TSF method and Subsection 3.3 with MM method. Extensions to the case of a varying intensity are presented in Section 4. Finally, numerical results with concluding remarks are exposed in Section 5.

## 2. Averaging

A general assumption throughout this paper is that the magnetic field is bounded from below, i.e. that for all for some independent of . In this section, we further assume that the external magnetic field has constant norm

 |B(x)|≡const>0,x∈R3,

so that the stiff part of equation (1.1) generates periodic motion. This setup has also been considered in [24]. A possible instance of such a is given by where

 B3=√∥B21+B22∥L∞(R2)−B21−B22,

with and . It can be verified that and . In such a case, we are able to apply a recently developed averaging method to quickly obtain the limit model of (1.1) as .

###### Lemma 2.1.

If for some constant , then the solution of

 ∂t~fε(t,x,v)+1εv×B(x)⋅∇v~fε(t,x,v)=0, (2.1)

is -periodic with respect to the fast time-variable .

###### Proof.

The characteristics of (2.1)

 ˙x(t)=0,˙v(t)=1εv(t)×B(x(t)),t>0,

have a periodic solution in which can be obtained, for instance, by Rodrigues’ formula

 x(t) =x(0), v(t) =cos(bt/ε)v(0)+(1−cos(bt/ε))(B(x(0))⋅v(0))B(x(0))−sin(bt/ε)v(0)×B(x(0)). (2.2)

The statement of the lemma is now an immediate consequence. ∎

Using the observation above, we may apply the following theorem from [5]:

###### Theorem 2.2.

Consider a transport equation of form

 ∂tfε(t,y)+[G(y)ε+K(y)]⋅∇yfε(t,y)=0,fε(0,y)=f0(y),

where the flow map of

 ˙y(t)=G(y(t))

is assumed to be -periodic. There exist two formalvector fields and satisfying

 G(y)ε+K(y)=Gε(y)ε+Kε(y) and [Gε,Kε]=0,

such that the system

 ∂τg(t,τ,y)+Gε(y)⋅∇yg(t,τ,y)=0, (2.3a) ∂tg(t,τ,y)+Kε(y)⋅∇yg(t,τ,y)=0, (2.3b) g(0,0,y)=f0(y). (2.3c)

has a unique formal solution independently of the order in which the equations are solved. Moreover, for all positive time we have and the first two terms of may be computed as follows

 K[1]=ΠKτ,K[2]=−12Π∫τ0[Ks,Kτ]ds, with Kτ(y):=(DyΦτ(y))−1(K∘Φτ)(y),

with .

###### Remark 2.3.

If is truncated at order in and , then the order in which the equations in (2.3) are solved does matter. However, the difference between the corresponding two solutions is also of order .

Without loss of generality, we assume in the rest of this section that to derive the averaged model of equation (1.1) as obtained from Theorem 2.2. Details of the derivation may be found in [5] and we thus content ourselves with a sketch of the computations: from (1.1), we have

 G=(0v×B),K=(vE),

and the flow map generated by is given by

 Φτ(y)=(xcos(τ)v+sin(τ)v×B+(1−cos(τ))(B⋅v)B),y=(xv).

Then

 DyΦτ(y)=(I30N1N2),

where

 N1=sin(τ)(v×∇xB)+(1−cos(τ))[(B⋅v)∇xB+B(vT∇xB)], N2=cos(τ)I3−sin(τ)(B×I3)+(1−cos(τ))BBT

where we have denoted and similarly for . A straightforward calculation then leads to

 Π((DyΦτ(y))−1(K⋅Φτ))=((B⋅v)BBv)

where

 Bv= BBT[E−12(v×∇xB)(v×B)+Mv−52M(B⋅v)B]−12M[v−2(B⋅v)B] −12B×I3[M(v×B)+(v×∇xB)(B⋅v)B],

with

 M=(B⋅v)∇xB+B(vT∇xB).

Eventually,

 Kε=((B⋅v)BBv)+O(ε),

and limit model at leading order is

 ∂tg(t,τ,x,v)+(B⋅v)B⋅∇xg(t,τ,x,v)+Bv⋅∇vg(t,τ,x,v)=0,t>0.

By taking and , we get the leading order averaged model of (1.1) for stroboscopic times ,

 ∂tf(t,x,v)+(B⋅v)B⋅∇xf(t,x,v)+Bv⋅∇vf(t,x,v)=0, (2.4a) ∇x⋅E(t,x)=∫R3f(t,x,v)dv−ni, (2.4b) f(0,x,v)=f0(x,v). (2.4c)

As we shall verify later numerically, we have

 fε(t,x,v)−f(t,x,v)=O(ε),0<ε≪1,t∈2πεN.

## 3. Numerical method

In this section, we introduce numerical schemes for equation (1.1) under the assumption . Taking advantage of the periodicity the solution of the stiff part, we apply state-of-art multiscale approaches in combination with PIC discretisation. In this way, we obtain schemes whose accuracy and computational cost are both independent of . Our starting point is the following PIC-representation of as used in, e.g., [19, 23, 26, 35],

 fε(t,x,v)≈Np∑k=1ωkδ(x−xk(t))δ(v−vk(t)),t≥0, x,v∈R2. (3.1)

The characteristic equations of model (1.1) for are then of the form

 ˙xk(t)=vk(t), (3.2a) ˙vk(t)=E(t,xk(t))+1εvk(t)×B(xk(t)),t>0, (3.2b) xk(0)=xk,0,vk(0)=vk,0. (3.2c)

Noticing that

 ∇x⋅E(t,x)=Np∑k=1wkδ(x−xk(t)),

we observe that the electric field in (3.2) has in fact no explicit dependence on time, i.e. where . We are in a position to briefly present three different UA methods.

### 3.1. Multi-revolution composition method

For a general exposition of multi-revolution composition (MRC), we refer to [8]. Here, we focus on a uniformly accurate second-order method.

MRC framework. Suppose that we wish to solve equation (1.1) on for some . Rescaling time in (3.2) leads to (we omit the particle index for brevity)

 ˙x(t)=εv(t), (3.3a) ˙v(t)=εE[X(t)](x(t))+v(t)×B(x(t)),0

Since the stiff part of (3.3) generates a -periodic motion, equation (3.3) is amenable to MRC [8, 9]. To do so, we write

 Tfε=2πMf+Tr,Mf=⌊Tf2πε⌋∈N,0≤Tr<2π. (3.4)

The 2nd order MRC method begins by choosing an integer and defining

 α=12(1+1M0),β=12(1−1M0),M=MfM0,H=εM0. (3.5)

Denoting , the MRC scheme proceeds as follows

 (xn+1vn+1)=Eβ(−2π)Eα(2π)(xnvn),0≤n≤M−1, (3.6)

where denotes the value at time of the flow of

 ˙x(t)=αHv(t), ˙v(t)=αHE[X(t)](x(t))+v(t)×B(x(t)), (3.7)

and the value at time of the flow of

 ˙x(t)=−βHv(t), ˙v(t)=−βHE[X(t)](x(t))+v(t)×B(x(t)). (3.8)

The solution at final time is then obtained by applying to the flow at time of

 ˙x(t)=εv(t), ˙v(t)=εE[X(t)](x(t))+v(t)×B(x(t)). (3.9)

Splitting scheme. The full MRC scheme calls for the numerical evaluation of the sub-flows , and . This is done here through a splitting, for instance of , in

 Exα(t):{˙x(s)=αHv(s),˙v(s)=0,0

Note that both and can be exactly integrated. The exact flow of is clearly

 x(t)=x(0)+tαHv(0),v(t)=v(0),t≥0,

while the exact flow of , by using the Rodrigues’ rotation formula, can also be written explicitly

 x(t)= x(0), v(t)= cos(t)v(0)+sin(t)v(0)×B+αHsin(t)E+αH(t−sin(t))(B⋅E)B+αH(1−cos(t))E×B +(1−cos(t))(B⋅v(0))B,t≥0,

where and . In our experiments, we shall take the value of the (micro) time step , so that

 Eα(2π)≈(Exα(h/2)Evα(h)Exα(h/2))M.

Approximations for and are obtained in a similar way. It may then be proved (see (3.6)) that the error of MRC is of size for a computational cost of size , making the overall scheme of order one111Under the assumption that and in (3.2)..

It remains to comment on what happens when the user-controlled increases to the limit where reaches the critical value , for which (3.5) and . In this case, the full MRC scheme may be regarded as just the discretisation of equation (3.3) by Strang’s method with time step . Therefore, as soon as implies , we replace MRC method by Strang splitting with time step . Finally, note that all vector fields involved in MRC are divergence free so that their exact flows are phase-space volume preserving as is the MRC method itself.

### 3.2. Two-scale formulation method

Two-scale formulation (TSF) methods have been developed in [6, 12]. Their underlying rationale is to consider the fast time as an additional variable. In order to isolate the fast time, we apply the change of unknowns where

 y(t)=cos(t/ε)v(t)+(1−cos(t/ε))(B(x(t))⋅v(t))B(x(t))−sin(t/ε)v(t)×B(x(t)). (3.11)

 ˙x(t)=Fx(t/ε,x(t),y(t)), ˙y(t)=Fy(t/ε,x(t),y(t)),t>0, (3.12) x(0)=x0,y(0)=v0,

where

 Fx(τ,x,y):=cos(τ)y+(1−cos(τ))(B(x)⋅y)B(x)+sin(τ)y×B(x),

and

 Fy(τ,x,y)= cos(τ)E[X](x)+(1−cos(τ))(B(x)⋅E[X](x))B(x)−sin(τ)E[X](x)×B(x) −12sin(2τ)qτ(y)−12(2sin(τ)−sin(2τ))qτ((B(x)⋅y)B(x)) −12(1−cos(2τ))qτ(y×B(x))+12(2cos(τ)−cos(2τ)−1)pτ(y)+12(3−4cos(τ) +cos(2τ))pτ((B(x)⋅y)B(x))+12(2sin(τ)−sin(2τ))pτ(y×B(x)),

with the vector fields

 pτ(z):=((∇xB(x)Fx(τ,x,y))⋅z)B(x)+(B(x)⋅z)(∇xB(x)Fx(τ,x,y)), qτ(z):=z×(∇xB(x)Fx(τ,x,y)),z∈R3.

Denoting and , the two-scale formulation of system (3.12) now reads

 ∂tU(t,τ)+1ε∂τU(t,τ)=F(τ,U(t,τ)),t>0, τ∈T, (3.13) U(0,0)=u(0),

where , and one recovers the solution of (3.12) by taking the diagonal, i.e.

 U(t,t/ε)=u(t),t≥0.

It remains to prescribe an appropriate initial data to (3.13) so that the solution has its derivatives uniformly bounded up to some order.

Initial data. In order to derive , we follow the Chapman-Enskog procedure. From the decomposition

 U––(t)=ΠU(t,⋅),h(t,τ)=U(t,τ)−–U(t), with ΠU(t,⋅)=12π∫2π0U(t,τ)dτ,

we split (3.13) into

 ˙U––(t)=ΠF(⋅,U––(t)+h(t,⋅)),t>0, ∂th(t,τ)+1ε∂τh(t,τ)=(I−Π)F(τ,U––(t)+h(t,τ)),t>0, τ∈T.

Denote and we have

 h(t,τ)=εAF(τ,U––(t)+h(t,τ))−εL−1∂th(t,τ).

Differentiate the above with respect to on both sides:

 ∂th(t,τ)=εA∇F(τ,U––+h)(˙U––+∂th)−εL−1∂2th(t,τ).

By assuming that for , one gets and has the first order asymptotic expansion:

 h(t,τ)=εAF(τ,U––(t))+O(ε2).

Using the fact one gets at initial time

 h(0,τ)=h1st(τ)+O(ε2),% withh1st(τ):=εAF(τ,u(0)),

and we denotes the first order initial data as:

 U1st(τ):=u(0)+h1st(τ)−h1st(0). (3.14)

In fact, one can show rigorously that the equation (3.13) with the well-prepared initial data offers

 ∂tU(t,τ),∂2tU(t,τ)=O(1),ε∈]0,1]. (3.15)

We refer the readers to [6] for the mathematical justification. The boundedness of the time derivatives (3.15) is the key to design UA schemes.

Exponential integrator. Thanks to the two-scale formulation (3.13) with the well-prepared initial data from (3.14), we can now directly apply the second order exponential integrator scheme proposed in [15] for integrating (3.13): Choose an even integer to uniformly discretize on and take a to define . Denote for and let . We update the for as

 (3.16a) ˆ(U)n+1l=e−ilΔtεˆ(U)nl+plˆ(F)nl+ql1Δt(ˆ(F)nl−ˆ(F)n−1l),n≥1, (3.16b)

where for ,

 Un(τ)=Nτ/2−1∑l=−Nτ/2ˆ(U)nleilτ,Fn(τ)=Nτ/2−1∑l=−Nτ/2ˆ(F)nleilτ,F∗,1(τ)=Nτ/2−1∑l=−Nτ/2ˆ(F)∗,1leilτ,

and with

 ˆ(U)∗,1l=e−ilΔtεˆ(U)0l+plˆ(F)0l,U∗,1(τ)=Nτ/2−1∑l=−Nτ/2ˆ(U)∗,1leilτ,

and

Suppose we have the numerical solution from the above scheme, then the numerical solution of the original characteristics (3.2) reads:

 xn=Xn(tn/ε),n≥1, vn=cos(tn/ε)Yn(tn/ε)+(1−cos(tn/ε))(B(xn)⋅Yn(tn/ε))B(xn)+sin(tn/ε)Yn(tn/ε)×B(xn).

The derivation and convergence analysis of the above scheme can be found in [15]. Since the filter (3.11) involves the magnetic field , the filtered system (3.12) which is less smooth than the original form (3.2), needs more regularity for optimal convergence of the algorithm. Assuming that and in (3.2), the two-scale formulation (TSF) exponential integrator (3.16) gives uniform second order accuracy in terms of for all and uniform spectral accuracy in terms of (due to periodicity):

 O(Δt2+N−m0τ).

The total cost of the TSF method is .

### 3.3. Micro-macro method

Now, we present the main new method of this work. It is based on the micro-macro (MM) decomposition that has been proposed very recently in [7]. We shall for the first time consider this approach for the Vlasov-Poisson equation and propose a second order UA scheme. The same notations will be adopted from the previous subsection.

MM decomposition. By the averaging theory [34], it is known that for general oscillatory problem

 ˙u(t)=F(t/ε,u(t)),t>0, (3.17)

with -periodicity in of , the solution can be written as a composition

 u(t)=Φt/ε∘Ψt∘Φ−10(u(0)), (3.18)

where is a change of variable with -periodicity in for some , and is the flow map of the autonomous equation with initial value :

 ˙Ψ