Drift approximation by the modified Boris algorithm of charged-particle dynamics in toroidal geometry

In this paper, we study the charged-particle dynamics under strong magnetic field in a toroidal axi-symmetric geometry. Using modulated Fourier expansions of the exact and numerical solutions, the long-term drift motion of the exact solution in toroidal geometry is derived and the error analysis of the large-stepsize modified Boris algorithm over long time scales is provided. Numerical experiments illustrate the theoretical results.


Introduction
The time integration of the equations of motion of charged particles is a crucial step in particle methods of plasma physics [1].In the strong magnetic field regime, the charged particles exhibit very fast rotations of small radius around a guiding centre.This often brings stringent restriction on the time stepsize for numerical integrators.There are many works aiming at designing largestepsize integrators with good accuracy for the charged-particle dynamics such as [12,4,3,11,14,7].Among them, a Boris-type integrator with appropriate modifications shows striking numerical results [14] and the rigorous analysis is provided in [10].It is proved that the position and the parallel velocity are approximated with O(h 2 ) accuracy for the modified Boris algorithm with large step sizes h 2 ∼ ε for fixed T = O (1), where ε 1 is a small parameter whose inverse corresponds to the strength of the magnetic field.
In this paper, we are interested in analyzing the long time behavior (over O(ε −1 )) of the modified Boris algorithm in a toroidal axi-symmetric geometry, with a magnetic field everywhere toroidal and an electric field everywhere orthogonal to the magnetic field.This geometry has already been proposed in [5] and a first order description of the slow dynamics for the continuous case is derived.Here we will use a different technique of modulated Fourier expansions [9], which recently has been used for charged-particle dynamics in a strong magnetic field [6,8,7,10,13], to derive the guiding centre drifts of the exact solution in such toroidal geometry.Since this technique can be extended to numerical discretization equally, the analysis of the modified Boris algorithm is also performed.
In Section 2 we formulate the equations of motion in a strongly non-uniform strong magnetic field, describe the concerned toroidal axi-symmetric geometry and the electromagnetic field, and introduce the modified Boris scheme.In Section 3 we state the main results of this paper: Theorem 3.1 states the slow drift motion over O(ε −1 ) in toroidal geometry for the continuous system and Theorem 3.2 states the long-time accuracy of the modified Boris algorithm.Section 4 presents experiments that illustrate the theoretical results.In Section 5 we give the proofs for our main results.

Charged-particle dynamics in toroidal geometry
We consider the differential equation that describes the motion of a charged particle (of unit mass and charge) in a magnetic and electric field, where x(t) ∈ R 3 is the position at time t, v(t) = ẋ(t) is the velocity, B(x) is the magnetic field and E(x) is the electric field.B and E can be expressed via vector potential A(x) ∈ R 3 and scalar potential φ(x) ∈ R as B(x) = ∇ × A(x) and E(x) = −∇φ(x).Here we are interested in the situation of a strong magnetic field where B 1 is smooth and independent of the small parameter ε, with |B 1 (x)| ≥ 1 for all x.The initial values (x(0), ẋ(0)) are bounded independently of ε: for some constants M 0 , M 1 , In this paper, we consider the so called toroidal axi-symmetric geometry which is introduced in [5].To be specific, fix a unitary vector e z , and for any vector x ∈ R 3 , it can be expressed as x = r(x) e r (x) + z(x) e z with z(x) = e z x, r(x) = |e z × x|, and e r (x) = (x − z(x) e z )/r(x).It is assumed that far from the axis e z the magnetic field is stationary, toroidal, axi-symmetric and non vanishing, that is, for some r 0 > 0, when r(x) ≥ r 0 e (x) = e z × x r(x) and for some function b.The electric field satisfies E (x) = 0 and E is axisymmetric when r(x) ≥ r 0 , that is, (2.5) In our proofs, we assume that the functions b, E r , E z and all their derivatives are bounded independently of ε.It is noted that (e r (x), e (x), e z ) forms the orthonormal basis and , 0 e z = (0, 0, 1) The following relations are useful in our proof e r (x) = 1 r(x) e (x)e (x) , e (x) = − 1 r(x) e r (x)e (x) where denotes the Jacobian of the functions considered and ∇ x is the gradient.

Modified Boris method
The modified Boris method proposed in [14] is used to solve the chargedparticle dynamics under a strong magnetic field with large stepsizes.The analysis of accuracy order of such method for general non-uniform strong magnetic field with stepsize h 2 ∼ ε until T = O(1) is recently provided in [10].This algorithm has the following two-step formulation with the initial magnetic moment The velocity is computed as The modified Boris method starts from modified initial values where P (x 0 ) is the orthogonal projection onto the span of B(x 0 ).This means the component of the initial velocity orthogonal to the magnetic field is filtered out, i.e., v 0 ⊥ = P ⊥ (x 0 )v 0 = 0 with P ⊥ (x 0 ) = I − P (x 0 ).We note that the modified Boris method is identical to the standard Boris integrator for the modified electric field E mod (x) = E(x) − µ 0 ∇|B|(x) = −∇(φ + µ 0 |B|)(x) and can be implemented as the common one-step formulation of the Boris algorithm [2].

Main results
Introducing r(t), z(t), ṽ(t) such that they are the solutions of the following initial-value problem for the slow differential equations we then have the following results.
Theorem 3.1 (Drift motion of the exact solution) Let x(t) = r(x(t)) e r (x(t))+ z(x(t)) e z be a solution of (2.1)-(2.3)with (2.4) and (2.5), which stays in a compact set K for 0 ≤ t ≤ cε −1 (with K and c independent of ε) and v (t) = e (x(t)) ẋ(t) be the parallel velocity.Denote r(t) = r(x(t)) and z(t) = z(x(t)), then we have The constant C is independent of ε and t with 0 ≤ t ≤ c/ε, but depends on c and on bounds of derivatives of B 1 and E on the compact set K.
Remark 3.1 A similar result is given by Proposition 5.2 in [5].Here we provide a different proof of the modulated Fourier expansions, which can be extended to the analysis of numerical methods and enables us to obtain the following result.
For the numerical approximation, the nondegeneracy condition is needed as in [10]: For (x, v) along the numerical trajectory, the linear maps have an inverse that is bounded independently of (x, v) and of h and ε with This determines an upper bound C * on the ratio h 2 /ε.Theorem 3.2 (Drift approximation by the numerical solution) Consider applying the modified Boris method to (2.1)-(2.3)with (2.4) and (2.5) and with modified initial values (2.9) using a step size h with h 2 ∼ ε, i.e., for some positive constants c * and C * .Under the nondegeneracy condition (3.2) and provided that the numerical solution The constant C is independent of ε and h and n with 0 ≤ nh ≤ c/ε, but depends on on c, on bounds of derivatives of B 1 and E on the compact set K and on c * and C * .

Numerical experiments
To illustrate the statement of the preceding section we consider the following electromagnetic fields   All the reference solutions are obtained by using standard Boris with small time step size h = 0.05ε.

Proof of main results
The theorems will be proved mainly based on the modulated Fourier expansions for the exact and numerical solutions given in [10].In this section, we will write the guiding centre equations in the toroidal geometry and express all the O(ε) terms explicitly.Following [6], we diagonalize the linear map v → v × B(x) and denote the eigenvalues as λ 1 = i|B(x)|, λ 0 = 0, and λ −1 = −i|B(x)|.The corresponding normalised eigenvectors are denoted by ν 1 (x), ν 0 (x), ν −1 (x) and the orthogonal projections onto the eigenspaces are denoted by P j (x) = ν j (x)ν j (x) * .It is noted that P (x) = P 0 (x) and P ⊥ (x) = I − P (x) = P 1 (x) + P −1 (x).

Proof of Theorem 3.1
The proof is structured into three parts (a)-(c).
(a) The equation of guiding centre motion in cartesian coordinate.
According to Theorem 4.1 of [6], it is known that the solution of (2.1)-(2.2) can be written as where the phase function satisfies φ(t) = |B 1 (y 0 (t))|.The coefficient functions y k (t) together with their derivatives (up to order N ) are bounded as where y k (t) = y k 1 (t) + y k 0 (t) + y k −1 (t) with y k j = P j (y 0 )y k .The remainder term and its derivative are bounded by Similar to Theorem 4.1 of [10], we can divide the interval [0, c/ε] into small intervals of length O(ε) and on each subinterval we consider the above modulated Fourier expansion, which means x(t) can be written as the modulated Fourier expansion for longer time intervals where y k (t) are piecewise continuous with jumps of size O(ε N ) at integral multiples of ε and are smooth elsewhere.The sizes of the coefficients and remainder term are the same as above.Inserting (5.1) into the continuous system and comparing the coefficients of e ikϕ(t)/ε yield the differential equations for y k (t).For k = 0 and k = ±1, we have ÿ0 = ẏ0 × B(y 0 ) + E(y 0 ) + 2Re i|B|y 1 × B (y 0 )y −1 From the first equation of (5.2), it is straightforward to get several slow drifts for P ⊥ ẏ0 (see Remark 4.3 of [10]) and the guiding centre motion of y 0 (t) satisfies with B, ∇|B|, P = e e and E evaluated at the guiding centre y 0 .The initial value of y 0 is (b) The equations in toroidal geometry.
where the functions E r , b, ∂ r b are evaluated at (r 0 , z 0 ).The initial value of z 0 is z 0 (0) = e z y 0 (0) = e z x(0 -By the definition of (5.5) we can derive the equation for v 0 d dt v 0 = d dt (e ẏ0 ) = ė ẏ0 + e ÿ0 . (5.8) The first term on the right hand side is r 0 dr 0 dt (5.9) using (5.6).In the following we will show that the second term e ÿ0 is of size O(ε 2 ).
From the expression of B given in (2.6), it is known that B (y 0 )y k ±1 is parallel to e , thus e II = 0 and e I = e 2 Re(i|B|y The algebraic equation of y ±1 0 can be derived by applying P (y 0 ) to the second equation of (5.2) The dominant term is − φ2 /ε 2 y ±1 0 and the right hand side P ẏ0 × B (y 0 )y ±1 ±1 = 0 using that B (y 0 )y k ±1 is parallel to e .Hence we obtain the following relation for y ±1 0 y ±1 0 = ±2i By differential and substitution the first term on the right hand side of above equation can be removed.Using (5.6), we have The initial value of v 0 is v 0 (0) = e (y 0 (0)) ẏ0 (0) = e (x(0)) ẋ(0) + O(ε).
(c) From short to long time intervals Denoting by y 0,[n] , r 0, [n] , z 0, [n] , v 0, [n] the functions y 0 , r 0 , z 0 , v 0 on time interval nε ≤ t ≤ (n + 1)ε, from (b), it is known that these coefficients satisfy the following equations ) with E r , E z , b, ∂ r b, ∂ z b evaluated at (r 0, [n] , z 0, [n] ).From equation (5.1), on every time interval, we have and thus In view of the factor ε in front of the right hand side of the differential equations (3.1) and (5.10), we have In view of the factor ε in front of the right hand side of the (5.10), we have With the above estimates, we obtain, for nε which is the stated result of Theorem 3.1.

Proof of Theorem 3.2
Similar to the proof of Theorem 3.1, we structure the proof into three parts.
(a) For a general strong magnetic field, the time interval of modulated Fourier expansion for numerical solution is validated over O(h).Using the uniqueness of the modulated Fourier expansion, we can patch together many short-time expansions in the same way as it was done for the exact solution and obtain the expansion for longer time O(1/ε).
From Theorem 4.2 of [10], it is known that the numerical solution x n given by the modified Boris algorithm (2.7)-(2.9)with a step size h satisfying c * ε ≤ h 2 ≤ C * ε can be written as where y 0 = O(1), y 1 = O(h 2 ) are peicewise continuous with jumps of size O(h N ) at integral multiples of h and smooth elsewhere.They are unique up to O(h N ) and P ⊥ (y 0 ) ẏ0 = O(h 2 ), P 0 (y 0 )y 1 = O(h 4 ).
ÿ0 − P±1 ẏ0 .By differentiation and substitution the derivatives of g ±1 can be removed and we have Using that Ṗ1 + Ṗ−1 + Ṗ0 = 0, we obtain ẏ0 = P ẏ0 + P 1 ẏ0 + P −1 ẏ0 (5.13) with B, ∇|B|, P = e e and E evaluated at the guiding centre y 0 .Compared to the guiding centre equation (5.3) of the exact solution, it is noticed that there are additional O(h 2 ) terms in (5.13).
(b) Next we derive the guiding centre equation in toroidal geometry where y 0 (t) can be written as y 0 = r(y 0 )e r (y 0 ) + z(y 0 )e z =: r 0 e r (y 0 ) + z 0 e z .
-Multiplying (5.13) with e z gives e z ẏ0 = − 2h 2 e z Ṗ ÿ0 − h 2 e z P ẏ0 where the first two terms on the right hand side vanish using (5.6) and the orthorgonality of e z , e r , e .Similar to the continuous case, we obtain where the functions E r , b, ∂ r b are evaluated at (r 0 , z 0 ).The initial value of z 0 is z 0 (0) = y 0 (0), e z ) = x 0 , e z ) + O(h 2 ) = z(x 0 ) + O(h 2 ).
-Finally, we need to derive the differential equation for v 0 which can be directly computed as the continuous case dv 0 dt = d dt e ẏ0 = ė ẏ0 + e ÿ0 . (5.17) The first term on the right hand side is the same as (5.9).Multiplying (5.12) with e = e (y 0 ) yields e ÿ0 = −h 2 e ....This gives e ÿ0 = h 2 v 0 r 0 2 dv 0 dt + O(h 4 ).
By patching together the errors together as what was done for the continuous case, we prove that r(x n )−r(t n ) = O(h 2 ), z(x n )−z(t n ) = O(h 2 ), v n −ṽ(t n ) = O(h 2 ), 0 ≤ t ≤ c/ε.

Acknowledgement
The author thanks Professor Christian Lubich for many useful discussions and comments.This work was supported by the Sino-German (CSC-DAAD) Postdoc Scholarship, Program No. 57575640.

Fig. 2 . 1
Fig. 2.1 Toroidal geometry with er(x), e (x), ez the local frame and the magnetic field along e (x).

Fig. 4 . 2
Fig. 4.2 Particle trajectories for t ≤ 1/ε with ε = 10 −3 projected onto the r − z plane as computed by the modified Boris with h = 0.04 (left) and by the Boris method with h = 0.01 (right).

Fig. 4 . 3
Fig. 4.3 Absolute errors |r(x n ) − r ref |, |z(x n ) − z ref | and |v n − v ref | as functions of time, along the numerical solution of the modified Boris algorithm with ε = 10 −3 and three different h.

Fig. 4 . 4
Fig. 4.4 Absolute errors |r(x n ) − r ref |, |z(x n ) − z ref | and |v n − v ref | as functions of time, along the numerical solution of the modified Boris algorithm with ε = 10 −4 and three different h.

Figure 4 .
Figure 4.1 shows the trajectories computed by the standard Boris and modified Boris with final time T = 1/ε.It is observed that the modified Boris method can give the correct trajectories even with very large time step size h = 40ε.For the standard Boris method, the drift motions are totally wrong with large time step size.The projection of the computed particle trajectory onto the (r, z) plane is given in Figure4.2.Figure4.3 shows the absolute errors of r, z, and v along the numerical solution of the modified Boris algorithm with ε = 10 −3 and T = 0.5/ε, which is observed to be size of O(h 2 ) in agreement with our theoretical results.Figure4.4 shows the similar results for ε = 10 −4 .All the reference solutions are obtained by using standard Boris with small time step size h = 0.05ε.

− v 0 r
3) with e r = e r (y 0 E − µ 0 ∇b) × e + O(ε 2 ), (5.4) where e r , e are evaluated at y 0 and v 0 := e ẏ0 .y 0 ) − ė r y 0 = dr 0 dt , and the first term on the right hand side of (5.4) vanishes since e r e × de dt = e r (e × e r ) = 0. Using the fact that e r × e = e z , e z × e = −e r , E = E r e r + E z e z and ∇b = ∂ r b e r + ∂ z b e z , we obtain e r (E − µ 0 ∇b(r 0 , z 0 )) × e = −E z + µ 0 ∂ z b.