# A filtered Boris algorithm for charged-particle dynamics in a strong magnetic field

A modification of the standard Boris algorithm, called filtered Boris algorithm, is proposed for the numerical integration of the equations of motion of charged particles in a strong non-uniform magnetic field in the asymptotic scaling known as maximal ordering. With an appropriate choice of filters, second-order error bounds in the position and in the parallel velocity, and first-order error bounds in the normal velocity are obtained with respect to the scaling parameter. The proof compares the modulated Fourier expansions of the exact and the numerical solutions. Numerical experiments illustrate the error behaviour of the filtered Boris algorithm.

## Authors

• 2 publications
• 12 publications
• 126 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 ...
01/25/2021

### Large-stepsize integrators for charged-particle dynamics over multiple time scales

The Boris algorithm, a closely related variational integrator and a newl...
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...
12/16/2021

### Geometric continuous-stage exponential energy-preserving integrators for charged-particle dynamics in a magnetic field from normal to strong regimes

This paper is concerned with geometric exponential energy-preserving int...
12/19/2018

### An arbitrary order time-stepping algorithm for tracking particles in inhomogeneous magnetic fields

The Lorentz equations describe the motion of electrically charged partic...
11/05/2019

### Numerical methods for antiferromagnetics

Compared with ferromagnetic counterparts, antiferromagnetic materials ar...
10/13/2020

### A simulation framework for particle magnetization dynamics of large ensembles of single domain particles: Numerical treatment of Brown/Néel dynamics and parameter identificati

Magnetic nanoparticles and their magnetization dynamics play an importan...
##### 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

In this paper we propose and analyse a numerical integrator for the equations of motion of a charged particle in a strong inhomogeneous magnetic field,

 ¨x(t)=˙x(t)×B(x(t),t)+E(x(t),t) (1) withB(x,t)=1εB0(εx)+B1(x,t)for  0<ε≪1.

This scaling is of interest in particle methods in plasma physics and is called maximal ordering in brizard07fon ; see also possanner18gfv for a careful discussion of scalings and a rigorous analysis of this model. It is assumed that , that , and are smooth functions that are bounded independently of on bounded domains together with all their derivatives, and that the initial data , are bounded independently of .

In (1), represents the position at time of a charged particle (of unit mass and charge) that moves in the magnetic field and the electric field . The motion is composed of fast rotation around a guiding center (with the Larmor radius proportional to ) and slow motion of the guiding center.

The standard integrator for charged particles in a magnetic field is the Boris algorithm boris70rps (see also, e.g., hairer18ebo ), which in the two-step formulation with step size is given by

 xn+1−2xn+xn−1h2=xn+1−xn−12h×B(xn,tn)+E(xn,tn) (2)

with the velocity approximation at time . This algorithm does, however, not behave well for (1) with small . Here we propose a modification, which we name filtered Boris algorithm. This modified integrator allows us to obtain better accuracy with considerably larger time steps, at minor additional computational cost. It is still a symmetric algorithm. We formulate and discuss this new algorithm in Section 2. It comes in different variants that depend on the choice of a suitable filter function and of the positions where the magnetic field is evaluated, and we identify favourable choices.

In Section 3 we state the main theoretical result, Theorem 3.1, which gives an error bound that behaves favourably with respect to . While most filters only yield a first-order error bound in the positions, for a particular, non-trivial choice of the filter a second-order error bound is obtained. A second-order error bound is also obtained for the component of the velocity that is parallel to the magnetic field. For the normal velocity approximation, there is only a first-order error bound for any filter. The proof of Theorem 3.1 is based on comparing the modulated Fourier expansions of the exact and the numerical solutions, which are derived in Sections 4 and 5, respectively. Combining those results, the proof of Theorem 3.1 is finally completed in Section 6.

We remark that the differential equations for the coefficient functions of the modulated Fourier expansions derived explicitly up to in Section 4 also yield the motion of the guiding center up to . They coincide up to with the guiding center equations of the numerical approximation given by the filtered Boris integrator for an appropriate filter and for non-resonant step sizes with a possibly large constant . This does not hold true for the standard Boris integrator.

In Section 7 we describe a related, but different integrator, called two-point filtered Boris algorithm, which evaluates the magnetic field both in the current position and in the current guiding center approximation in each step, and which has similar convergence properties to the previously considered filtered Boris method.

In Section 8 we present the results of numerical experiments in which we compare the standard and filtered Boris algorithms.

In the Appendix we show how the filters are evaluated efficiently using a Rodriguez-type formula.

## 2 Filtered Boris algorithm

Using the velocity approximation at the mid-point,

 vn−1/2=1h(xn−xn−1)=vn−h2vn×B(xn,tn)−h2E(xn,tn),

the Boris algorithm (2) is usually written and implemented as a one-step method ,

 vn−1/2+=vn−1/2+h2E(xn,tn)vn+1/2−−vn−1/2+=h2(vn+1/2−+vn−1/2+)×B(xn,tn)vn+1/2=vn+1/2−+h2E(xn,tn)xn+1=xn+hvn+1/2. (3)

To capture the high oscillations in the velocity more accurately, the second line of (3) needs to be modified, and one should rather work with averaged velocities and possibly averaged positions. This can be achieved with the help of filter functions like in garcia-archilla99lmf ; hochbruck99agm and (hairer06gni, , Section XIII.2).

For a vector

we denote by the Euclidean norm of and we use the common notation

 v×B=−ˆBv,ˆB=⎛⎜⎝0−b3b2b30−b1−b2b10⎞⎟⎠. (4)

For real-analytic functions (such as ) we will form matrix functions , which are efficiently computed by a Rodriguez-type formula as described in the Appendix.

We denote by

 xn⊙=xn+vn×Bn|Bn|2 (5)

with the guiding center approximation at time (cf. northrop63tam ). For the argument of  in the algorithm we choose a point on the straight line connecting and :

 ¯xn=θnxn+(1−θn)xn⊙ (6)

with for a real function . It turns out that there is a unique choice of such that a second-order error bound will be obtained:

 θ(ξ)=1sinc(ξ/2)2, (7)

where . We note that with the scaling (1), we have , provided that is bounded away from non-zero integral multiples of .

We consider the following modification of the Boris algorithm.

###### Algorithm 2.1 (Filtered Boris algorithm)

Given , the algorithm computes as follows, with , with defined by (6), and :

 vn−1/2+=vn−1/2+h2Ψ(hˆBn)En$$vn+1/2−=exp(−hˆ¯Bn)vn−1/2+.$$vn+1/2=vn+1/2−+h2Ψ(hˆBn)Enxn+1=xn+hvn+1/2, (8)

where with .

The velocity approximation is computed as

 vn=Φ1(hˆ¯Bn)xn+1−xn−12h−hΥ(hˆBn)En, (9)

where with , and . The starting approximation is computed from (12) below with .

For the choice , the algorithm is explicit and requires only matrix-vector multiplications that can be done efficiently with a Rodriguez-type formula (see the Appendix).

For with from (7), the algorithm is implicit, because then depends on and appears in the argument of in the second line. This can be solved by a rapidly convergent fixed-point iteration for , with an error reduction by a factor in each iteration. We start the iteration with , then compute from (8) and from (9) using

 12h(xn+1−xn−1)=12(vn+1/2+vn−1/2)=12(vn+1/2−+vn−1/2+). (10)

This then yields from (5) and the new from (6). In practice, one iteration is sufficient to get the improved accuracy. We note that all matrix-vector multiplications can be done with a Rodriguez-type formula.

We mention that Algorithm 2.1 preserves volume in phase space exactly in the case of constant (and time-dependent ), but only approximately up to in the general case of an inhomogeneous magnetic field (1).

Two-step formulation. The filtered Boris algorithm has the symmetric two-step formulation

 xn+1−2xn+xn−1h2=2htanh(−12hˆ¯Bn)xn+1−xn−12h+Ψ(hˆBn)En, (11)

as is readily obtained by taking two consecutive steps and using (10). This formulation is the basis of our theoretical analysis.

Starting value. The starting value is chosen such that formulas (8)-(9) also hold for . We find, for arbitrary , that

 vn±1/2=φ1(∓hˆ¯Bn)(vn+hΥ(hˆBn)En)±h2Ψ(hˆBn)En, (12)

where . Note that, for given and , the vectors and are known, so that (9) provides an explicit formula for .

One-step map . Using the last formula of (8) together with (12) for relating and , and (12) with and and with and for relating and , the filtered Boris algorithm can be written as

 xn+1=xn+hΦn+vn+h22Ψn+EnΦn+1−vn+1=Φn+vn+h2Ψn+1−En+1+h2Ψn+En, (13)

where and .

The method is symmetric in the sense that exchanging and gives the same formulas.

The integrator in the case of a constant magnetic field. For constant , we note that , and so (13) reduces to the exponential integrator (with the notation )

 xn+1=xn+hφ1(−hˆB)vn+h22Ψ+(−hˆB)Envn+1=exp(−hˆB)vn+h2(Ψ0(−hˆB)En+Ψ1(−hˆB)En+1) (14)

with and . The method is exact for a constant magnetic field and vanishing electric field , because

 exp(0tI0−tˆB)=(I  tφ1(−tˆB)0 exp(−tˆB)). (15)

Since we have chosen , the method is also exact for constant and . This is seen as follows: For constant , the variation-of-constants formula for the system reads, in view of (15),

 x(tn+h)=x(tn)+hφ1(−hˆB)v(tn)+ h2∫10(1−s)φ1(−(1−s)hˆB)E(x(tn+hs))ds,v(tn+h)=exp(−hˆB)v(tn)+h∫10exp(−(1−s)hˆB)E(x(tn+hs))ds.

For constant , this becomes (14), which yields , where .

## 3 Statement of the main result

Our main theoretical result in this paper is the following error bound for the filtered Boris algorithm. Here we denote, for the exact velocity ,

 v∥(t)=B(x(t),t)|B(x(t),t)|(B(x(t),t)|B(x(t),t)|⋅v(t)),v⊥(t)=v(t)−v∥(t),

and similarly for the numerical velocity ,

 vn∥=B(xn,t)|B(xn,tn)|(B(xn,tn)|B(xn,tn)|⋅vn),vn⊥=vn−vn∥.

We then have the following result.

###### Theorem 3.1

We assume the following, with arbitrarily chosen positive constants , , and :

1. The initial velocity satisfies an -independent bound

 |v0|≤M. (16)
2. The exact solution of (1) stays in a bounded set (independent of ) for .

3. The step size satisfies and is such that the following non-resonance condition is satisfied:

 ∣∣sinc(12kh|B(x(t),t)|)∣∣≥c>0for k=1,2,3. (17)

If in the filtered Boris algorithm,

• is given by (6) with the function of (7), and

• the filter functions and are defined as in Algorithm 2.1,

then the errors in the positions and the velocities are bounded by

 xn−x(tn) =O(ε2), (18) vn∥−v∥(tn) =O(ε2),vn⊥−v⊥(tn)=O(ε).

For a different choice of the functions , and , the error bounds are not better than for general problems (1). The constants in the -notation are independent of and and with , but depend on , on the velocity bound and the constants and , and on bounds of derivatives of , and in a neighbourhood of the set .

We remark that in view of the error bounds, the non-resonance condition might be required along the numerical solution instead of the exact solution as in (17).

The proof of this theorem will compare the modulated Fourier expansion of the exact solution (as given in Section 4) with that of the numerical approximation (as given in Section 5). It will be given in Section 6.

###### Remark 1

The proof also shows that the choice is sufficient for order 2 if the magnetic field satisfies, for all and and all times ,

 Im(z×∂xB(x,t)¯z)⋅B(x,t)=O(ε).

## 4 Modulated Fourier expansion of the exact solution

We write the solution of (1) as a modulated Fourier expansion

 x(t)≈∑k∈Zzk(t)eikϕ(t)/ε (19)

with coefficient functions for which all time derivatives are bounded independently of , where , and describes the motion of the guiding center. Such a formal expansion has first been considered in kruskal58tgo

for proving the existence of an adiabatic invariant (essentially the magnetic moment

). It has been used for a rigorous proof of the long-time near-conservation of the magnetic moment in hairer19lta , where this approach was extended to the numerical solution of a variational integrator, for which near-conservation of the magnetic moment and of the energy is rigorously proved over long times that cover arbitrary negative powers of .

Following hairer19lta , we diagonalize the linear map

, which has eigenvalues

, , and

. We denote the normalized eigenvectors by

, and remark that is collinear to . We let

be the orthogonal projections onto the eigenspaces. Furthermore, we write the coefficient functions of (

19) in the time-dependent basis ,

 zk=zk1+zk0+zk−1,zkj(t)=Pj(z0(t),t)zk(t). (20)

Since is real, we assume for all . Together with the fact that and is real, it follows

 z−k−1=¯¯¯¯¯zk1,z−k0=¯¯¯¯¯zk0,z−k1=¯¯¯¯¯¯¯¯zk−1. (21)

The following result is a variant of Theorem 4.1 in hairer19lta , adapted to the present case of a strong magnetic field of the form (1). Note that in this paper corresponds to in hairer19lta .

###### Theorem 4.1

Let be a solution of (1) with bounded initial velocity (16) that stays in a compact set for . For an arbitrary truncation index we then have

 x(t)=∑|k|≤Nzk(t)eikϕ(t)/ε+RN(t), (22)

where the phase function satisfies .

(a) The coefficient functions together with their derivatives (up to order ) are bounded as for , , , for , , and for the remaining with ,

 zkj=O(ε|k|+1). (23)

They are unique up to . Moreover, we have .

(b) The remainder term and its derivative are bounded by

 RN(t)=O(t2εN),˙RN(t)=O(tεN)for0≤t≤T. (24)

(c) The functions satisfy the differential equations

 ¨z00 =P0(z0,t)E(z0,t)+2P0(z0,t)Re(i˙ϕεz11×B′(z0,t)z−1−1) +2˙P0(z0,t)˙z0+¨P0(z0,t)z0+O(ε2), (25) ˙z0±1 =˙P±1(z0,t)z0±iε˙ϕP±1(z0,t)E(z0,t)+O(ε2), (26) ˙z±1±1 =−¨ϕ˙ϕz±1±1+O(ε2)=O(ε2), (27)

where we use the notation and similar for . All other coefficient functions are given by algebraic expressions depending on .

(d) Assuming , initial values for the differential equations of item (c) are given by

 z0(0) =x(0)+˙x(0)×B(x(0),0)|B(x(0),0)|2+O(ε2), ˙z00(0) =P0(x(0),0)˙x(0)+˙P0(x(0),0)x(0)+O(ε2), z±1±1(0) =∓iε˙ϕ(0)P±1(x(0),0)˙x(0)+O(ε2).

The constants symbolised by the -notation are independent of and with , but they depend on , on the velocity bound in (16), on bounds of derivatives of  and , and on .

###### Remark 2

We note that the guiding center motion of the system (1) is given by the non-oscillating term in the modulated Fourier expansion. By the uniqueness of the modulated Fourier expansion up to high powers of , the equations in item (d) hold not only at time 0, but for all .

###### Proof

(a) and (b): Compared to Theorem 4.1 in hairer19lta , where a more general strong magnetic field is considered, the time interval of validity of the modulated Fourier expansion is here instead of just , and the bound (23) is improved by a factor . The improvement of the time scale comes about by observing that a function that solves (1) up to a defect , i.e.,

 ¨x∗(t)=˙x∗(t)×B(x∗(t),t)+E(x∗(t),t)+d(t),

satisfies an error bound, for ,

 |x∗(t)−x(t)|≤C(|x∗(0)−x(0)|+|˙x∗(0)−˙x(0)|+∫t0|d(t)|dt),

where is independent of but grows exponentially with . This is proved by decomposing and using the variation-of-constants formula and the Gronwall inequality. The improvement of the bound (23) is a consequence of the fact that the derivatives of are bounded independently of .

(c): For the error bound of Section 6 we need precise formulas for the dominant terms of (22). Inserting the expansion (19) into the differential equation (1) and comparing the coefficients of yields

 ¨zk+2ik˙ϕε˙zk+(ik¨ϕε−k2˙ϕ2ε2)zk=Fk, (28)

where, using Taylor series expansion for the nonlinearities,

 Fk=∑k1+k2=k(˙zk1+ik1˙ϕεzk1)×∑m≥0s(α)=k21m!B(m)(z0,t)zα+∑m≥0s(α)=k1m!E(m)(z0,t)zα.

Here, and denote the th derivative with respect to , is a multi-index with , , , and .

From (28) it follows that the motion of the guiding center is given by

 (29)

The solution is influenced by the functions which, by (28), satisfy

 ±2i˙ϕε˙z±1+(±i¨ϕε−˙ϕ2ε2)z±1=(˙z±1±i˙ϕεz±1)×B(z0,t)+O(ε). (30)

Note that, whereas is of size , its derivatives are bounded independently of due to the special form (1).

To get solutions with derivatives bounded uniformly in , one has to extract the dominant terms. Multiplying (29) with eliminates the -term that is present in , and the second derivative becomes dominant. Differentiating the relation with respect to time yields . This then gives (25). Note that, due to the special form of , the time derivatives of are of size .

A multiplication of (29) with gives

 P±1(z0,t)¨z0=±i˙ϕεP±1(z0,t)˙z0+P±1(z0,t)E(z0,t)+O(ε).

Substituting , and extracting yields (26). Note that , so that also , and .

Since , the -terms cancel in (30) after projection with . Therefore, the -terms are dominant and we obtain (27).

(d): Assuming , initial values are determined from (22) by

 x(0)=z0(0)+(z1(0)+z−1(0))+O(ε3)˙x(0)=˙z0(0)+(˙z1(0)+˙z−1(0))+i˙ϕ(0)ε(z1(0)−z−1(0))+O(ε3). (31)

This is a nonlinear system for . We write the vectors in the basis , and we select the dominant terms in each equation. They are in the upper relation of (31), and in the lower relation. Fixed-point iteration then yields the stated equations for the initial values. Note that the relation has been applied. ∎

## 5 Modulated Fourier expansion of the numerical solution

We consider the two-step formulation (11) of the filtered Boris algorithm, and we write the numerical approximation as

 xn≈∑k∈Zzk(t)eikϕ(t)/ε,t=nh. (32)

We use the same notation for the coefficient functions as in Section 4. Note, however, that for the numerical solution these functions are not the same and depend on the additional parameter . We again consider the basis and the corresponding orthogonal projections , and we write the coefficient functions as in (20), with the only difference that here the argument is the non-oscillating part of (32) and not that of (19).

###### Theorem 5.1

Let be a numerical solution of the filtered Boris algorithm applied to (1) with bounded initial velocity (16), and suppose that it stays in a compact set  for . We assume the non-resonance condition

 (33)

for a fixed, but arbitrary truncation index , and (for convenience of presentation) also the bound . Moreover, we assume that the filter function in (11) is bounded by for all real , where . Then, we have that

 xn=∑|k|≤Nzk(t)eikϕ(t)/ε+RN(t),t=nh, (34)

where the phase function is given by .

(a) and (b) The coefficient functions together with their derivatives (up to order ) as well as the remainder term and its derivative satisfy the bounds of items (a) and (b) of Theorem 4.1.

(c) The functions