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,

$$\begin{aligned} \begin{aligned}&\ddot{x}(t) = \dot{x}(t) \times B(x(t),t) + E(x(t),t) \\&\text {with}\quad B(x,t) = \frac{1}{\varepsilon }\, B_0(\varepsilon x ) + B_1(x,t) \quad \text {for }\ 0<\varepsilon \ll 1 . \end{aligned} \end{aligned}$$
(1.1)

This scaling is of interest in particle methods in plasma physics and is called maximal ordering in [2]; see also [12] for a careful discussion of scalings and a rigorous analysis of this model. It is assumed that \({|B_0(0)|\ge 1}\), that \(B_0\), \(B_1\) and E are smooth functions that are bounded independently of \(\varepsilon \) on bounded domains together with all their derivatives, and that the initial data \(x(0)=x^0\), \(\dot{x}(0)=v^0\) are bounded independently of \(\varepsilon \).

In (1.1), \(x(t)\in {\mathbf{R}}^3\) represents the position at time t of a charged particle (of unit mass and charge) that moves in the magnetic field B and the electric field E. The motion is composed of fast rotation around a guiding center (with the Larmor radius proportional to \(\varepsilon \)) and slow motion of the guiding center.

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

$$\begin{aligned} \frac{x^{n+1} - 2x^n + x^{n-1}}{h^2} = \frac{x^{n+1} - x^{n-1}}{2h} \times B(x^n,t^n) + E(x^n,t^n) \end{aligned}$$
(1.2)

with the velocity approximation \(v^n = \frac{1}{2h}\bigl ( x^{n+1} - x^{n-1}\bigr ) \) at time \(t^n=nh\). This algorithm does, however, not behave well for (1.1) with small \(\varepsilon \). 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 Sect. 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.

A different approach to numerical integrators for charged particles in a strong magnetic field is taken in [3], where damping linearly implicit integrators are studied. For step sizes \(h\ge \varepsilon \), those integrators quickly reduce that part of the kinetic energy that comes from the velocity component perpendicular to the magnetic field. As is shown in [3], the numerical solution then approximates a numerical solution of an asymptotic model for guiding-center motion. In contrast, the original Boris method and the filtered Boris method proposed here do not dampen the high oscillations, and the filtered method is shown to give approximations of second-order accuracy to both the actual position and the actual guiding center over an \({\mathcal O}(1)\) time scale for stepsizes \(h\sim \varepsilon \).

In Sect. 3 we state the main theoretical result, Theorem 3.1, which gives an error bound for the filtered Boris method that behaves favourably with respect to \(\varepsilon \). 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 still 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 Sects. 4 and 5, respectively. Combining those results, the proof of Theorem 3.1 is finally completed in Sect. 6.

We remark that the differential equations for the coefficient functions of the modulated Fourier expansions derived explicitly up to \({\mathcal O}(\varepsilon ^2)\) in Sect. 4 also yield the motion of the guiding center up to \({\mathcal O}(\varepsilon ^2)\). They coincide up to \({\mathcal O}(\varepsilon ^2)\) 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 \(h\le C\varepsilon \) with a possibly large constant C. This does not hold true for the standard Boris integrator.

In Sect. 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 Sect. 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,

$$\begin{aligned} v^{n-1/2} = \frac{1}{h} \bigl (x^{n} - x^{n-1} \bigr ) = v^n - \frac{h}{2}\,v^n \times B(x^n,t^n) - \frac{h}{2} E(x^n,t^n), \end{aligned}$$

the Boris algorithm (1.2) is usually written and implemented as a one-step method \((x^n,v^{n-1/2}) \mapsto (x^{n+1},v^{n+1/2})\),

$$\begin{aligned} \begin{array}{rcl} v^{n-1/2}_+ &{}=&{} v^{n-1/2} + \frac{h}{2}\, E(x^n,t^n) \\ v^{n+1/2}_- - v^{n-1/2}_+ &{}=&{} \frac{h}{2}\, \bigl (v^{n+1/2}_- + v^{n-1/2}_+\bigr ) \times B(x^n,t^n)\\ v^{n+1/2} &{}=&{} v^{n+1/2}_- + \frac{h}{2}\, E(x^n,t^n) \\ x^{n+1} &{}=&{} x^n + h\, v^{n+1/2} . \end{array} \end{aligned}$$
(2.1)

To capture the high oscillations in the velocity more accurately, the second line of (2.1) needs to be modified, and one should rather work with averaged velocities \( v^{n+1/2} \approx \tfrac{1}{h}\int _{t^n}^{t^{n+1}} v(t)\, {\mathrm d}t\) and possibly averaged positions. This can be achieved with the help of filter functions like in [4, 9] and [8, Section XIII.2].

For a vector \(B=(b_1,b_2,b_3)^\top \in {\mathbf{R}}^3\) we denote by |B| the Euclidean norm of B and we use the common notation

$$\begin{aligned} v\times B = - \widehat{B} \,v,\qquad \widehat{B} = \begin{pmatrix} 0 &{} - b_3 &{} b_2\\ b_3&{}0&{}-b_1\\ -b_2&{}b_1 &{} 0 \end{pmatrix} . \end{aligned}$$
(2.2)

For real-analytic functions \(\Psi (\zeta )\) [such as \(\exp (\zeta )\)] we will form matrix functions \(\Psi (-h\widehat{B})\), which are efficiently computed by a Rodriguez-type formula as described in the Appendix.

We denote by

$$\begin{aligned} x^n_\odot = x^n + \frac{v^n \times B^n}{|B^n|^2} \end{aligned}$$
(2.3)

with \(B^n=B(x^n,t^n)\) the guiding center approximation at time \(t^n\) (cf. [11]). For the argument of B in the algorithm we choose a point on the straight line connecting \(x^n\) and \(x^n_\odot \):

$$\begin{aligned} \bar{x}^n = \theta ^n x^n + (1-\theta ^n) x^n_\odot \end{aligned}$$
(2.4)

with \(\theta ^n=\theta (h |B^n|)\) for a real function \(\theta \). It turns out that there is a unique choice of \(\theta \) such that a second-order error bound will be obtained:

$$\begin{aligned} \theta (\xi ) = \frac{1}{{{\,\mathrm{sinc}\,}}(\xi /2)^2}, \end{aligned}$$
(2.5)

where \({{\,\mathrm{sinc}\,}}(\xi )=\sin (\xi )/\xi \). We note that with the scaling (1.1), we have \(\bar{x}^n = x^n +{\mathcal O}(\varepsilon )\), provided that \(h|B^n|\) is bounded away from non-zero integral multiples of \(2\pi \).

We consider the following modification of the Boris algorithm.

Algorithm 2.1

(Filtered Boris algorithm) Given \((x^n , v^{n-1/2})\), the algorithm computes \((x^{n+1} , v^{n+1/2})\) as follows, with \(B^n=B(x^n,t^n)\), \(\bar{B}^n=B(\bar{x}^n,t^n)\) with \(\bar{x}^n\) defined by (2.4), and \(E^n=E(x^n,t^n)\):

$$\begin{aligned} \begin{aligned} v^{n-1/2}_+&= v^{n-1/2} + \frac{h}{2}\, \Psi (h\widehat{B^n})\,E^n \\ v^{n+1/2}_-&= \exp \bigl (-h\widehat{\bar{B}^n}\bigr ) v^{n-1/2}_+. \\ v^{n+1/2}&= v^{n+1/2}_- + \frac{h}{2}\, \Psi (h\widehat{B^n})\,E^n\\ x^{n+1}&= x^n + h\, v^{n+1/2} , \end{aligned} \end{aligned}$$
(2.6)

where \(\,\Psi (\zeta )=\mathrm {tanch}(\zeta /2)\,\) with \(\,\mathrm {tanch}(\zeta )=\tanh (\zeta )/\zeta \).

The velocity approximation \(v^n\) is computed as

$$\begin{aligned} v^n = \Phi _1(h\widehat{\bar{B}^n}) \,\frac{x^{n+1}-x^{n-1}}{2h} - h\Upsilon (h\widehat{B^n}) E^n , \end{aligned}$$
(2.7)

where \(~\displaystyle \Phi _1 (\zeta ) = \frac{1}{{{\,\mathrm{sinch}\,}}(\zeta )}~\) with \(~\displaystyle {{\,\mathrm{sinch}\,}}(\zeta ) = \frac{ \sinh (\zeta )}{\zeta }\), and \(~\displaystyle \Upsilon (\zeta ) =\frac{ \Phi _1 (\zeta ) -1 }{\zeta } \). The starting approximation \(v^{1/2}\) is computed from (2.10) below with \(n=0\).

Theorem 3.1 below shows that the choice of filter functions made in Algorithm 2.1 and (2.5) is the unique choice that gives a second-order error bound.

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

For \(\theta ^n = \theta (h|B^n|)\) with \(\theta (\zeta ) \) from (2.5), the algorithm is implicit, because \(\bar{x}^n\) then depends on \(v^n\) and appears in the argument of \(\bar{B}^n\) in the second line. This can be solved by a rapidly convergent fixed-point iteration for \(\bar{x}^n\):

We start the iteration with \(\bar{x}^n = x^n\), then compute \(v_-^{n+1/2}\) from (2.6) and \(v^n\) from (2.7) using

$$\begin{aligned} \tfrac{1}{2h}(x^{n+1} - x^{n-1}) = \tfrac{1}{2}\bigl (v^{n+1/2} + v^{n-1/2}\bigr ) = \tfrac{1}{2}\bigl (v^{n+1/2}_- + v^{n-1/2}_+\bigr ). \end{aligned}$$
(2.8)

This then yields \(x^n_\odot \) from (2.3) and the new \(\bar{x}^n\) from (2.4). We note that all matrix-vector multiplications can be done with a Rodriguez-type formula.

It is readily checked that this iteration map has a Lipschitz constant of size \({\mathcal O}(\varepsilon ^2)\) under the scaling (1.1) and for stepsizes \(h={\mathcal O}(\varepsilon )\). Since the deviation from the solution of the implicit method is therefore reduced by a factor \({\mathcal O}(\varepsilon ^2)\) in each iteration, making just one iteration is sufficient to get the improved second-order accuracy of the implicit method.

We mention that Algorithm 2.1 preserves volume in phase space exactly in the case of constant B (and time-dependent B(t)), but only approximately up to \({\mathcal O}(h\varepsilon )\) in the general case of an inhomogeneous magnetic field (1.1). For the original Boris method it is, however, known from [13] that the map \((x_n,v_{n-1/2})\mapsto (x_{n+1},v_{n+1/2})\) is volume-preserving for general magnetic fields.

6pt

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

$$\begin{aligned} \frac{x^{n+1} - 2x^n + x^{n-1}}{h^2} = \frac{2}{h}\mathrm {tanh}\bigl ( -\tfrac{1}{2} h\widehat{\bar{B}^n}\bigr ) \,\frac{x^{n+1} - x^{n-1}}{2h} + \Psi ( h \widehat{B^n} ) E^n , \end{aligned}$$
(2.9)

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

Starting value The starting value \(v^{1/2}\) is chosen such that formulas (2.6)–(2.7) also hold for \(n=0\). We find, for arbitrary n, that

$$\begin{aligned} v^{n\pm 1/2} = \varphi _1 \bigl (\mp h \widehat{\bar{B}^n} \bigr ) \Bigl ( v^n + h\Upsilon (h\widehat{B^n}) E^n\Bigr ) \pm \frac{h}{2} \Psi (h\widehat{B^n})\,E^n , \end{aligned}$$
(2.10)

where \(\varphi _1(\zeta )=(e^\zeta -1)/\zeta \). Note that, for given \(x^0\) and \(v^0\), the vectors \(x^n_\odot \) and \(\bar{x}^n\) are known, so that (2.7) provides an explicit formula for \(v^{1/2}\). 6pt

One-step map\({(x^n,v^n)\mapsto (x^{n+1},v^{n+1})}\) Using the last formula of (2.6) together with (2.10) for relating \(x^{n+1}\) and \(x^n\), and (2.10) with n and \(+\) and with \(n+1\) and − for relating \(v^{n+1}\) and \(v^n\), the filtered Boris algorithm can be written as

$$\begin{aligned} \begin{array}{rcl} x^{n+1} &{}=&{} \displaystyle x^n + h \Phi ^n_+ v^n + \tfrac{h^2}{2} \,\Psi ^n_+ E^n \\ \Phi ^{n+1}_- v^{n+1} &{}=&{} \Phi ^{n}_+ v^{n} + \tfrac{h}{2} \,\Psi ^{n+1}_- E^{n+1} + \tfrac{h}{2} \,\Psi ^n_+ E^n , \end{array} \end{aligned}$$
(2.11)

where \(\Phi ^n_\pm = \varphi _1(\mp h\widehat{\bar{B}^n})\) and \(\Psi ^n_\pm = \Psi (h\widehat{B^n} )\pm 2\Phi ^n_\pm \Upsilon (h\widehat{B^n})\).

The method is symmetric in the sense that exchanging \(n\leftrightarrow n+1\) and \(h\leftrightarrow -h\) gives the same formulas.

6pt

The integrator in the case of a constant magnetic field For constant B, we note that \((\Phi ^{n+1}_-)^{-1} \Phi ^{n}_+ = \exp (-h\widehat{B})\), and so (2.11) reduces to the exponential integrator (with the notation \(\Psi _\pm (\zeta )=\Psi (\zeta ) \mp 2\varphi _1(\pm \zeta )\Upsilon (\zeta )\))

$$\begin{aligned} \begin{array}{rcl} x^{n+1} &{}=&{} x^n + h \varphi _1 (-h\widehat{B}) v^n + \frac{h^2}{2} \Psi _+ (-h\widehat{B}) E^n \\ v^{n+1} &{}=&{} \exp (-h\widehat{B}) v^n + \frac{h}{2} \bigl (\Psi _0 (-h\widehat{B}) E^n + \Psi _1 (-h\widehat{B}) E^{n+1} \bigr ) \end{array} \end{aligned}$$
(2.12)

with \(\Psi _0(\zeta )=\Psi _+(\zeta )/\varphi _1(-\zeta )\) and \(\Psi _1(\zeta )=\Psi _-(\zeta )/\varphi _1(-\zeta )\). The method is exact for a constant magnetic field B and vanishing electric field E, because

$$\begin{aligned} \exp \left( \begin{array}{c@{\quad }c} 0 &{} tI \\ 0 &{} -t \widehat{B} \\ \end{array} \right) = \left( \begin{array}{cc} I &{} \ \ t\,\varphi _1(-t \widehat{B}) \\ 0 &{} \ \exp (- t \widehat{B}) \\ \end{array} \right) . \end{aligned}$$
(2.13)

Since we have chosen \(\Psi (\zeta ) =\mathrm {tanch}(\zeta /2)\), the method is also exact for constant B and E. This is seen as follows: For constant B, the variation-of-constants formula for the system \(\dot{x}=v, \ \ \dot{v}= x \times B + E(x)\) reads, in view of (2.13),

$$\begin{aligned} x(t_n+h)= & {} x(t_n)+ h\varphi _1(-h \widehat{B}) v(t_n) \\&+\ h^2 \int _{0}^1(1- s ) \varphi _1(-(1- s )h \widehat{B}) E(x(t_n+h s )) {\mathrm d}s ,\\ v(t_n+h)= & {} \exp (-h \widehat{B}) v(t_n)+h \int _{0}^1 \exp (-(1- s )h \widehat{B}) E(x(t_n+h s )) {\mathrm d}s . \end{aligned}$$

For constant E, this becomes (2.12), which yields \(\Psi _\pm (\zeta )=\varphi _2(\pm \zeta )\), where \(\varphi _2(\zeta )=(e^\zeta -1 - \zeta )/(\zeta ^2/2) = \int _0^1 (1- s ) \,\varphi _1((1- s )\zeta ) {\mathrm d}s \).

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)=\dot{x}(t)\),

$$\begin{aligned} v_\parallel (t) = \frac{B(x(t),t)}{|B(x(t),t)|}\, \biggl ( \frac{B(x(t),t)}{|B(x(t),t)|}\cdot v(t) \biggr ), \qquad v_\perp (t)=v(t)-v_\parallel (t), \end{aligned}$$

and similarly for the numerical velocity \(v^n\),

$$\begin{aligned} v_\parallel ^n = \frac{B(x^n,t)}{|B(x^n,t^n)|}\, \biggl ( \frac{B(x^n,t^n)}{|B(x^n,t^n)|}\cdot v^n \biggr ), \qquad v_\perp ^n=v^n-v_\parallel ^n. \end{aligned}$$

We then have the following result.

Theorem 3.1

We assume the following, with arbitrarily chosen positive constants c, C, M and T:

  1. 1.

    The initial velocity satisfies an \(\varepsilon \)-independent bound

    $$\begin{aligned} |v^0| \le M. \end{aligned}$$
    (3.1)
  2. 2.

    The exact solution x(t) of (1.1) stays in a bounded set K (independent of \(\varepsilon \)) for \(0\le t \le T\).

  3. 3.

    The step size satisfies \(h \le C \varepsilon \) and is such that the following non-resonance condition is satisfied:

    $$\begin{aligned} \big |{{\,\mathrm{sinc}\,}}\bigl (\tfrac{1}{2} kh |B(x(t),t)|\bigr )\big | \ge c >0 \quad \text {for } k=1,2,3. \end{aligned}$$
    (3.2)

If in the filtered Boris algorithm,

  • \(\bar{x}^n\) is given by (2.4) with the function \(\theta \) of (2.5), and

  • the filter functions \(\Psi \) and \(\Upsilon \) are defined as in Algorithm 2.1,

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

$$\begin{aligned} x^n - x(t^n)= & {} {\mathcal O}(\varepsilon ^2), \nonumber \\ v_\parallel ^n - v_\parallel (t^n)= & {} {\mathcal O}(\varepsilon ^2), \qquad v_\perp ^n - v_\perp (t^n) ={\mathcal O}(\varepsilon ). \end{aligned}$$
(3.3)

For a different choice of the functions \(\theta \), \(\Psi \) and \(\Upsilon \), the error bounds are not better than \({\mathcal O}(\varepsilon )\) for general problems (1.1). The constants in the \({\mathcal O}\)-notation are independent of \(\varepsilon \) and h and n with \(0\le t^n=nh\le T\), but depend on T, on the velocity bound M and the constants c and C, and on bounds of derivatives of \(B_0\), \(B_1\) and E in a neighbourhood of the set K.

We remark that in view of the error bounds, the non-resonance condition might be required along the numerical solution \(x^n\) instead of the exact solution x(t) as in (3.2).

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

Remark 3.1

The proof also shows that the choice \(\bar{x}^n = x^n\) is sufficient for order 2 if the magnetic field satisfies, for all \(z\in {\mathbf{C}}^3\) and \(x\in K\) and all times t,

$$\begin{aligned} {\mathrm {Im}\,}(z \times \partial _x B(x,t) \bar{z}) \cdot B(x,t) = {\mathcal O}(\varepsilon ). \end{aligned}$$

Remark 3.2

As a consequence of Theorem 3.1, the approximate guiding center \(x^n_\odot \) of the numerical solution, given by (2.3), is an \({\mathcal O}(\varepsilon ^2)\) approximation to

$$\begin{aligned} x_\odot (t) = x(t) + \frac{\dot{x}(t) \times B(x(t),t)}{|B(x(t),t)|^2}, \end{aligned}$$

which is an \({\mathcal O}(\varepsilon ^2)\) approximation to the guiding center of the trajectory x(t) (by Theorem 4.1 below; cf. [11]).

4 Modulated Fourier expansion of the exact solution

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

$$\begin{aligned} x(t) \approx \sum _{k\in {\mathbf{Z}}} z^k (t) \,{\mathrm e}^{\mathrm {i}k\phi (t)/\varepsilon } \end{aligned}$$
(4.1)

with coefficient functions \(z^k(t)\) for which all time derivatives are bounded independently of \(\varepsilon \), where \(\dot{\phi }(t) /\varepsilon = \big | B\bigl (z^0(t),t\bigr )\big |\), and \(z^0(t)\) describes the motion of the guiding center. Such a formal expansion has first been considered in [10] for proving the existence of an adiabatic invariant (essentially the magnetic moment \(\tfrac{1}{2}|\dot{x} \times B(x)|^2/|B(x)|^3\)). It has been used for a rigorous proof of the long-time near-conservation of the magnetic moment in [7], 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 \(\varepsilon \).

Following [7], we diagonalize the linear map \(v\mapsto v\times B(x,t)\), which has eigenvalues \(\lambda _1 = \mathrm {i}|B(x,t)|\), \(\lambda _0 =0\), and \(\lambda _{-1}=-\mathrm {i}|B(x,t)|\). We denote the normalized eigenvectors by \(v_1(x,t), v_0(x,t),v_{-1}(x,t)\), and remark that \(v_0(x,t)\) is collinear to B(x, t). We let \(P_j(x,t)=v_j(x,t)v_j(x,t)^*\) be the orthogonal projections onto the eigenspaces. Furthermore, we write the coefficient functions of (4.1) in the time-dependent basis \(v_j\bigl (z^0(t),t\bigr )\),

$$\begin{aligned} z^k = z_1^k + z_0^k + z_{-1}^k ,\qquad z_j^k(t) =P_j\bigl (z^0(t),t\bigr )z^k(t) . \end{aligned}$$
(4.2)

Since x(t) is real, we assume \(z^{-k} = \overline{z^k}\) for all k. Together with the fact that \(v_{-1}(x,t) = \overline{v_1(x,t)}\) and \(v_0(x,t)\) is real, it follows

$$\begin{aligned} z_{-1}^{-k} = \overline{z_1^k},\quad z_0^{-k} = \overline{z_0^k},\quad z_{1}^{-k} = \overline{z_{-1}^k} . \end{aligned}$$
(4.3)

The following result is a variant of Theorem 4.1 in [7], adapted to the present case of a strong magnetic field of the form (1.1). Note that B in this paper corresponds to \(B/\varepsilon \) in [7].

Theorem 4.1

Let x(t) be a solution of (1.1) with bounded initial velocity (3.1) that stays in a compact set K for \(0\le t \le T\). For an arbitrary truncation index \(N\ge 1\) we then have

$$\begin{aligned} x(t) = \sum _{|k|\le N} z^k (t)\, {\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon } + R_N (t) , \end{aligned}$$
(4.4)

where the phase function satisfies \(\dot{\phi }(t) = \varepsilon | B (z^0(t),t ) |={\mathcal O}(1)\).

  1. (a)

    The coefficient functions \(z^k(t)\) together with their derivatives (up to order N) are bounded as \(z_j^0 = {\mathcal O}(1)\) for \(j\in \{-1,0,1\}\), \(z_1^1 = {\mathcal O}(\varepsilon )\), \(z_{-1}^{-1} = {\mathcal O}(\varepsilon )\), \(z_j^k = {\mathcal O}(\varepsilon ^3)\) for \(|k|=1\), \(j\ne k\), and for the remaining (j, k) with \( |k| \le N\),

    $$\begin{aligned} z_j^k = {\mathcal O}(\varepsilon ^{|k|+1}). \end{aligned}$$
    (4.5)

    They are unique up to \({\mathcal O}(\varepsilon ^{N+2})\). Moreover, we have \(\dot{z}^0 \times B(z^0,t) = {\mathcal O}(1)\).

  2. (b)

    The remainder term and its derivative are bounded by

    $$\begin{aligned} R_N(t) = {\mathcal O}(t^2 \varepsilon ^N), \quad \dot{R}_N(t) = {\mathcal O}(t \varepsilon ^N) \quad \hbox {for}\quad 0\le t\le T. \end{aligned}$$
    (4.6)
  3. (c)

    The functions \(z_0^0, z_{\pm 1}^0, z_1^1,z_{-1}^{-1}\) satisfy the differential equations

    $$\begin{aligned} \ddot{z}_0^0= & {} P_0(z^0,t) E(z^0,t) + 2\, P_0(z^0,t) \, {\mathrm {Re}\,} \Bigl ( \mathrm {i}\frac{\dot{\phi }}{\varepsilon }z_1^1 \times B'(z^0,t)z_{-1}^{-1}\Bigr )\nonumber \\&+ 2\, \dot{P}_0(z^0,t) \dot{z}^0 + \ddot{P}_0(z^0,t) z^0 + {\mathcal O}(\varepsilon ^2 ) , \end{aligned}$$
    (4.7)
    $$\begin{aligned} \dot{z}_{\pm 1}^0= & {} \dot{P}_{\pm 1}(z^0,t) z^0 \pm \mathrm {i}\frac{\varepsilon }{\dot{\phi }} P_{\pm 1}(z^0,t) E(z^0,t) + {\mathcal O}(\varepsilon ^2), \end{aligned}$$
    (4.8)
    $$\begin{aligned} \dot{z}_{\pm 1}^{\pm 1}= & {} - \frac{\ddot{\phi }}{\dot{\phi }}z_{\pm 1}^{\pm 1} + {\mathcal O}(\varepsilon ^2)={\mathcal O}(\varepsilon ^2), \end{aligned}$$
    (4.9)

    where we use the notation \(\dot{P}_j(z^0,t) = \frac{{\mathrm d}}{{\mathrm d}t} P_j\bigl (z^0(t),t\bigr )\) and similar for \(\ddot{P}_j(z^0,t)\). All other coefficient functions \(z_j^k\) are given by algebraic expressions depending on \(z^0, \dot{z}_0^0,z_1^1,z_{-1}^{-1}\).

  4. (d)

    Assuming \(\phi (0)=0\), initial values for the differential equations of item (c) are given by

    $$\begin{aligned} z^0(0)&= x(0) + \frac{\dot{x}(0) \times B(x(0),0)}{|B(x(0),0)|^2} + {\mathcal O}(\varepsilon ^2), \\ \dot{z}_0^0(0)&= P_0(x(0),0) \dot{x}(0)+ \dot{P}_0(x(0),0) x(0)+ {\mathcal O}(\varepsilon ^2 ) , \\ z_{\pm 1}^{\pm 1} (0)&= \mp \,\mathrm {i}\frac{\varepsilon }{\dot{\phi }(0)}P_{\pm 1}(x(0),0)\dot{x}(0) + {\mathcal O}(\varepsilon ^2). \end{aligned}$$

    The constants symbolised by the \({\mathcal O}\)-notation are independent of \(\varepsilon \) and t with \(0\le t \le T\), but they depend on N, on the velocity bound M in (3.1), on bounds of derivatives of B and E, and on T.

Remark 4.1

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

Proof

(a) and (b) Compared to Theorem 4.1 in [7], where a more general strong magnetic field is considered, the time interval of validity of the modulated Fourier expansion is here \({\mathcal O}(1)\) instead of just \({\mathcal O}(\varepsilon )\), and the bound (4.5) is improved by a factor \(\varepsilon \). The improvement of the time scale comes about by observing that a function \( x_*(t)\) that solves (1.1) up to a defect d(t), i.e.,

$$\begin{aligned} \ddot{x}_*(t) = \dot{x}_*(t) \times B(x_*(t),t) + E(x_*(t),t) +d(t), \end{aligned}$$

satisfies an error bound, for \(0\le t \le T\),

$$\begin{aligned} |x_*(t) - x(t) | \le C \Bigl ( |x_*(0)-x(0)| + |\dot{x}_*(0)-\dot{x}(0)| + \int _0^t |d(t)| \, {\mathrm d}t\Bigr ), \end{aligned}$$

where C is independent of \(\varepsilon \) but grows exponentially with T. This is proved by decomposing \(B(x,t)= \varepsilon ^{-1} B_0(0) + \varepsilon ^{-1}(B_0(\varepsilon x)-B_0(0)) + B_1(x,t)\) and using the variation-of-constants formula and the Gronwall inequality. The improvement of the bound (4.5) is a consequence of the fact that the derivatives of B(x, t) are bounded independently of \(\varepsilon \).

(c) For the error bound of Sect. 6 we need precise formulas for the dominant terms of (4.4). Inserting the expansion (4.1) into the differential Eq. (1.1) and comparing the coefficients of \({\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon }\) yields

$$\begin{aligned} \ddot{z}^k + 2\mathrm {i}k \frac{\dot{\phi }}{\varepsilon }\dot{z}^k+ \Bigl ( \mathrm {i}k \frac{\ddot{\phi }}{\varepsilon }- k^2 \frac{\dot{\phi }^2}{\varepsilon ^2}\Bigr ) z^k = F^k, \end{aligned}$$
(4.10)

where, using Taylor series expansion for the nonlinearities,

$$\begin{aligned} F^k = \sum _{k_1+k_2 = k} \Bigl ( \dot{z}^{k_1} + \mathrm {i}k_1 \frac{\dot{\phi }}{\varepsilon } z^{k_1}\Bigr ) \times \sum _{\begin{array}{c} m\ge 0 \\ s(\alpha ) = k_2 \end{array}} \frac{1}{m!} B^{(m)} (z^0,t) \,{\mathbf{z}}^\alpha + \sum _{\begin{array}{c} m\ge 0 \\ s(\alpha ) = k \end{array}} \frac{1}{m!} E^{(m)} (z^0,t) \,{\mathbf{z}}^\alpha . \end{aligned}$$

Here, \(B^{(m)}(x,t)\) and \(E^{(m)}(x,t)\) denote the mth derivative with respect to x, \(\alpha =(\alpha _1,\ldots ,\alpha _m)\) is a multi-index with \(\alpha _j \in {\mathbf{Z}}{\setminus }\{0\}\), \(s(\alpha ) = \alpha _1 + \cdots + \alpha _m\), \(|\alpha | = |\alpha _1| + \cdots + |\alpha _m|\), and \({\mathbf{z}}^\alpha = (z^{\alpha _1}, \ldots , z^{\alpha _m})\).

From (4.10) it follows that the motion of the guiding center \(z^0(t)\) is given by

$$\begin{aligned} \ddot{z}^0 = \dot{z}^0 \times B(z^0,t) + E(z^0,t) + 2 \,{\mathrm {Re}\,}\Bigl ( \mathrm {i}\frac{\dot{\phi }}{\varepsilon }z^1 \times B'(z^0,t)z^{-1}\Bigr ) + {\mathcal O}(\varepsilon ^2 ) . \end{aligned}$$
(4.11)

The solution \(z^0(t)\) is influenced by the functions \(z^{\pm 1}\) which, by (4.10), satisfy

$$\begin{aligned} \pm 2\mathrm {i}\frac{\dot{\phi }}{\varepsilon }\dot{z}^{\pm 1}+ \Bigl ( \pm \mathrm {i}\frac{\ddot{\phi }}{\varepsilon }- \frac{\dot{\phi }^2}{\varepsilon ^2}\Bigr ) z^{\pm 1} = \Bigl (\dot{z}^{\pm 1} \pm \mathrm {i}\frac{\dot{\phi }}{\varepsilon }z^{\pm 1} \Bigr ) \times B(z^0,t) + {\mathcal O}(\varepsilon ) . \end{aligned}$$
(4.12)

Note that, whereas \(B(z^0,t)\) is of size \({\mathcal O}(\varepsilon ^{-1})\), its derivatives are bounded independently of \(\varepsilon \) due to the special form (1.1).

To get solutions with derivatives bounded uniformly in \(\varepsilon \), one has to extract the dominant terms. Multiplying (4.11) with \(P_0(z^0,t)\) eliminates the \(\varepsilon ^{-1}\)-term that is present in \(B(z^0,t)\), and the second derivative \(\ddot{z}_0^0\) becomes dominant. Differentiating the relation \(z_0^0 = P_0(z^0,t) z^0\) with respect to time yields \(\ddot{z}_0^0 = P_0(z^0,t) \ddot{z}^0 + 2\dot{P}_0(z^0,t) \dot{z}^0 + \ddot{P}_0(z^0,t) z^0\). This then gives (4.7). Note that, due to the special form of B(x, t), the time derivatives of \(P_j(z^0,t)\) are of size \({\mathcal O}(\varepsilon )\).

A multiplication of (4.11) with \(P_{\pm 1}(z^0,t)\) gives

$$\begin{aligned} P_{\pm 1}(z^0,t) \ddot{z}^0 = \pm \mathrm {i}\frac{\dot{\phi }}{\varepsilon }P_{\pm 1}(z^0,t) \dot{z}^0 + P_{\pm 1}(z^0,t)E(z^0,t ) + {\mathcal O}(\varepsilon ) . \end{aligned}$$

Substituting \(P_{\pm 1}(z^0,t) \dot{z}^0 = \dot{z}_{\pm 1}^0 - \dot{P}_{\pm 1}(z^0,t) z^0 \), and extracting \(\dot{z}_{\pm 1}^0\) yields (4.8). Note that \(\dot{z}_{\pm 1}^0 ={\mathcal O}(\varepsilon )\), so that also \(\ddot{z}_{\pm 1}^0 ={\mathcal O}(\varepsilon )\), and \( P_{\pm 1}(z^0,t) \ddot{z}^0 = {\mathcal O}(\varepsilon )\).

Since \(\dot{\phi }/ \varepsilon = | B(z^0,t)|\), the \(\varepsilon ^{-2}\)-terms cancel in (4.12) after projection with \(P_{\pm 1}(z^0,t)\). Therefore, the \(\varepsilon ^{-1}\)-terms are dominant and we obtain (4.9).

(d): Assuming \(\phi (0)=0\), initial values are determined from (4.4) by

$$\begin{aligned} \begin{array}{rcl} x(0) &{}=&{} z^0 (0) + \bigl ( z^1 (0) + z^{-1} (0)\bigr ) + {\mathcal O}(\varepsilon ^3) \\[1mm] \dot{x} (0) &{}=&{}\displaystyle \dot{z}^0 (0) + \bigl ( \dot{z}^1 (0) + \dot{z}^{-1} (0)\bigr ) + \mathrm {i}\frac{ \dot{\phi }(0) }{\varepsilon }\bigl ( z^1 (0) - z^{-1} (0)\bigr ) + {\mathcal O}(\varepsilon ^3 ) . \end{array} \end{aligned}$$
(4.13)

This is a nonlinear system for \(z^0(0), \dot{z}_0^0(0), z_1^1(0), z_{-1}^{-1}(0)\). We write the vectors in the basis \(\{v_j(z^0(0),0)\}\), and we select the dominant terms in each equation. They are \(z^0(0)\) in the upper relation of (4.13), and \(\dot{z}_0^0 (0), z_1^1(0), z_{-1}^{-1}(0)\) in the lower relation. Fixed-point iteration then yields the stated equations for the initial values. Note that the relation \(P_j(z^0(0),0) = P_j(x(0),0) + {\mathcal O}(\varepsilon ^2 )\) has been applied. \(\square \)

5 Modulated Fourier expansion of the numerical solution

We consider the two-step formulation (2.9) of the filtered Boris algorithm, and we write the numerical approximation \(x^n\) as

$$\begin{aligned} x^n \approx \sum _{k\in {\mathbf{Z}}} z^k (t) \,{\mathrm e}^{\mathrm {i}k\phi (t)/\varepsilon } ,\qquad t=nh . \end{aligned}$$
(5.1)

We use the same notation for the coefficient functions as in Sect. 4. Note, however, that for the numerical solution these functions are not the same and depend on the additional parameter h. We again consider the basis \(\{ v_j (x,t) \}\) and the corresponding orthogonal projections \(P_j(x,t)\), and we write the coefficient functions \(z^k\) as in (4.2), with the only difference that here the argument \(z^0(t)\) is the non-oscillating part of (5.1) and not that of (4.1).

Theorem 5.1

Let \(\{ x^n \}\) be a numerical solution of the filtered Boris algorithm applied to (1.1) with bounded initial velocity (3.1), and suppose that it stays in a compact set K for \(0\le nh \le T\). We assume the non-resonance condition

$$\begin{aligned} \big |{{\,\mathrm{sinc}\,}}\bigl (\tfrac{1}{2} kh |B(x^n,t^n)|\bigr )\big | \ge c >0 \qquad \text {for } k=1,\dots ,N+1, \end{aligned}$$
(5.2)

for a fixed, but arbitrary truncation index \(N\ge 2\), and (for convenience of presentation) also the bound \(\eta = h/\varepsilon \le C\). Moreover, we assume that the filter function \(\Psi \) in (2.9) is bounded by \(|\Psi (\mathrm {i}\xi )| \le C \,|\mathrm {tanc}(\tfrac{1}{2} \xi )|\) for all real \(\xi \), where \(\mathrm {tanc}(\xi )=\tan (\xi )/\xi \). Then, we have that

$$\begin{aligned} x^n = \sum _{|k|\le N} z^k (t)\, {\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon } + R_N (t) ,\qquad t=nh, \end{aligned}$$
(5.3)

where the phase function is given by \(\dot{\phi }(t) = \varepsilon | B (z^0(t),t ) |\).

  1. (a) and (b)

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

  2. (c)

    The functions \(z_0^0, z_{\pm 1}^0, z_1^1,z_{-1}^{-1}\) satisfy the differential equations (with \(\theta (\xi )\) used in the definition of \(\bar{x}^n\) in (2.4))

    $$\begin{aligned} \ddot{z}_0^0= & {} P_0(z^0,t) E(z^0,t) + 2\, P_0(z^0,t) \, {\mathrm {Re}\,} \Bigl ( \mathrm {i}\frac{\dot{\phi }}{\varepsilon }z_1^1 \times B'(z^0,t)z_{-1}^{-1}\Bigr ) \nonumber \\&\times \theta (\eta \dot{\phi }) {{\,\mathrm{sinc}\,}}(\eta \dot{\phi }/2)^2+ 2 \dot{P}_0(z^0,t) \dot{z}^0 + \ddot{P}_0(z^0,t) z^0 + {\mathcal O}(\varepsilon ^2 ) , \end{aligned}$$
    (5.4)
    $$\begin{aligned} \dot{z}_{\pm 1}^0= & {} \dot{P}_{\pm 1}(z^0,t) z^0 \pm \frac{\displaystyle \Psi \bigl ( \mathrm {i}{\eta \dot{\phi }} \bigr )}{\displaystyle \mathrm {tan c} \Bigl ( \frac{\eta \dot{\phi }}{2} \Bigr )} \,\mathrm {i}\,\frac{\varepsilon }{\dot{\phi }} \, P_{\pm 1}(z^0,t) E(z^0,t) + {\mathcal O}(\varepsilon ^2), \end{aligned}$$
    (5.5)
    $$\begin{aligned} \dot{z}_{\pm 1}^{\pm 1}= & {} -\,\frac{1}{ \displaystyle \mathrm {tan c} \Bigl ( \frac{\eta \dot{\phi }}{2} \Bigr )} \frac{\ddot{\phi }}{\dot{\phi }}\,z_{\pm 1}^{\pm 1} + {\mathcal O}(\varepsilon ^2) = {\mathcal O}(\varepsilon ^2) . \end{aligned}$$
    (5.6)

    All other coefficient functions \(z_j^k\) are given by algebraic expressions depending on \(z^0, \dot{z}_0^0,z_1^1,z_{-1}^{-1}\).

  3. (d)

    Assuming \(\phi (0)=0\), initial values for the differential equations of item (c) are given by the same equations as for the exact solution, up to \({\mathcal O}(\varepsilon ^2)\),

    $$\begin{aligned} \begin{aligned} z^0(0)&= x(0) + \frac{\dot{x}(0) \times B(x(0),0)}{|B(x(0),0)|^2} + {\mathcal O}(\varepsilon ^2),\\ \dot{z}_0^0(0)&= P_0(x(0),0) \dot{x}(0)+ \dot{P}_0(x(0),0) x(0)+ {\mathcal O}(\varepsilon ^2 ) , \\ z_{\pm 1}^{\pm 1} (0)&= \mp \,\mathrm {i}\frac{\varepsilon }{\dot{\phi }(0)}P_{\pm 1}(x(0),0)\dot{x}(0) + {\mathcal O}(\varepsilon ^2). \end{aligned} \end{aligned}$$
    (5.7)

    The constants symbolised by the \({\mathcal O}\)-notation are independent of \(\varepsilon \) and n with \(0\le nh \le T\), but they depend on N, on the velocity bound M in (3.1), on bounds of derivatives of B and E, and on T.

Proof

(a) and (b) We do not present the details of the proof of the existence of the modulated Fourier expansion and the bounds for the coefficient functions and the remainder term, since this uses the same kind of arguments as in previous such proofs, e.g. in [5, 7, 8]. In particular, for \(|k|=1, j\ne k\) and for \(|k| \ge 2\) the construction of the coefficient functions (see part (c) below) shows that \(z_j^k\) is multiplied by

$$\begin{aligned} \frac{4}{\eta ^2} \sin \Bigl ( \frac{k\eta \dot{\phi }}{2} \Bigr ) \sin \Bigl ( \frac{(k-j)\eta \dot{\phi }}{2} \Bigr ) . \end{aligned}$$

Under the non-resonance assumption (5.2) this expression is bounded from below by a positive constant, so that an algebraic relation for \(z_j^k\) can be extracted.

By construction of the coefficient functions the truncated series of (5.3) satisfies the two-step relation (2.9) up to a defect of size \({\mathcal O}(\varepsilon ^N)\). A standard discrete Gronwall argument then gives the bounds on the remainder.

(d) The initial values are obtained from

$$\begin{aligned} x(0) = z^0(0) + \bigl ( z^1(0) + z^{-1}(0)\bigr ) + {\mathcal O}(\varepsilon ^2), \end{aligned}$$
(5.8)

which is a consequence of (5.1), and from

$$\begin{aligned} \dot{x}(0)= & {} \Phi _1\bigl (h\widehat{B}(x(0),0)\bigr ) \dot{z}^0(0) + \frac{\mathrm {i}\dot{\phi }(0)}{\varepsilon }\, z^1_1(0) - \frac{\mathrm {i}\dot{\phi }(0)}{\varepsilon }\, z^{-1}_{-1}(0)\nonumber \\&- h\Upsilon \bigl (h\widehat{B}(x(0),0)\bigr ) E(x(0),0) + {\mathcal O}(\varepsilon ^2), \end{aligned}$$
(5.9)

which follows from (2.7) and Lemma 5.1. As in the proof of Theorem 4.1 this constitutes a nonlinear system for the values \(z^0(0), \dot{z}_0^0(0), z_1^1 (0), z_{-1}^{-1}(0)\). The relation (5.8) yields \(z^0_0(0)\). Multiplication of (5.9) with \(P_j\bigl ( z^0(0),0\bigr ) = P_j \bigr ( x(0),0 \bigr ) + {\mathcal O}(\varepsilon ^2 )\) gives \(\dot{z}_0^0(0)\) for \(j=0\) and \(z_{\pm 1}^{\pm 1}(0) \) for \(j=\pm 1\), where we use in addition that \(\Phi _1(h\widehat{B}(x(0),0))=\Phi _1(h\widehat{B}(z^0(0),0))+{\mathcal O}(\varepsilon ^2)\) and \(P_{\pm 1}\dot{z}^0=\dot{z}^0_{\pm 1} - \dot{P}_{\pm 1} z^0 = {\mathcal O}(\varepsilon )\). Remarkably we get, up to terms of size \({\mathcal O}(\varepsilon ^2 )\), the same formulas for the initial values as for the exact solution.

By the uniqueness of the modulated Fourier expansion [up to \({\mathcal O}(\varepsilon ^N)\)], these relations hold not only at time 0, but for arbitrary times \(t\le T\), except for a phase factor \(e^{\mp \mathrm {i}\phi (t)}\) in the equation for \(z^{\pm 1}_{\pm 1}\). This phase factor did not appear in (5.7) because we chose \(\phi (0)=0\).

(c) To derive the differential equations for the coefficient functions we first expand the perturbed argument of B(x, t) in the filtered Boris algorithm as

$$\begin{aligned} \bar{x}^n \approx \sum _{k\in {\mathbf{Z}}} \zeta ^k (t) \,{\mathrm e}^{\mathrm {i}k\phi (t)/\varepsilon } ,\qquad t=nh . \end{aligned}$$
(5.10)

The coefficient functions \(\zeta ^k(t)\) are obtained as follows: inserting the modulated Fourier expansion (5.1) into (2.7), using Lemma 5.1 below, and replacing \(\Phi _1\bigl ( h\widehat{\bar{B}^n}\bigr )\) by \(\Phi _1\bigl ( h\widehat{B} (z^0(t^n),t^n)\bigr )\) yields with \(t=nh\)

$$\begin{aligned} v^n = \dot{z}_0^0 (t) + \frac{\mathrm {i}\dot{\phi }(t)}{\varepsilon } z_1^1(t)\,{\mathrm e}^{\mathrm {i}\phi (t)/\varepsilon } - \frac{\mathrm {i}\dot{\phi }(t)}{\varepsilon } z_{-1}^{-1}(t)\,{\mathrm e}^{-\mathrm {i}\phi (t)/\varepsilon } + {\mathcal O}(\varepsilon ) , \end{aligned}$$

see also the more detailed computation in Sect. 6. Since we have \(\dot{z}_0^0 (t) = P_0\bigl ( z^0(t),t\bigr ) z^0(t) + {\mathcal O}(\varepsilon )\) and \(B^n = \widehat{B} \bigl (z^0(t^n),t^n\bigr )+ {\mathcal O}(\varepsilon )\), this implies

$$\begin{aligned} \frac{v^n\times B^n}{|B^n|^2} = - z_1^1(t)\,{\mathrm e}^{\mathrm {i}\phi (t)/\varepsilon } - z_{-1}^{-1}(t)\,{\mathrm e}^{-\mathrm {i}\phi (t)/\varepsilon } + {\mathcal O}(\varepsilon ^2 ) , \end{aligned}$$

and consequently \(x^n_\odot = z^0 (t^n) + {\mathcal O}(\varepsilon ^2 )\), which shows that \(x^n_\odot \) is an excellent approximation of the non-oscillating part of the numerical solution \(x^n\). Together with the definition (2.4) of \(\bar{x}^n\) we find the dominating terms of the expansion (5.10) as

$$\begin{aligned} \zeta ^0(t) =z^0(t) + {\mathcal O}(\varepsilon ^2 ), \qquad \zeta ^{\pm 1} (t) = \theta (t)\, z^{\pm 1} (t) + {\mathcal O}(\varepsilon ^2 ), \end{aligned}$$
(5.11)

where \(\theta (t)=\theta (h\dot{\phi }(t)/\varepsilon )=\theta \bigl (h|B(z^0(t),t)|\bigr )\).

After this preparation, we insert the expansions (5.1) for \(x^n\) and (5.10) for \(\bar{x}^n\) into the two-step formulation (2.9) of the filtered Boris algorithm. Using Lemma 5.1 below, expanding the nonlinearities around \(\zeta ^0\) and \(z^0\), and comparing the coefficients of \({\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon }\) yields

$$\begin{aligned} \displaystyle \sum _{l\ge 0} \varepsilon ^{l-2} d_l^k \frac{{\mathrm d}^l}{{\mathrm d}t^l} z^k= & {} \displaystyle \sum _{k_1+k_2 = k} \Bigl ( \sum _{\begin{array}{c} m\ge 0 \\ s(\alpha ) = k_1 \end{array}} \frac{1}{m!} T_{\widehat{B}}^{(m)} (\zeta ^0,t) \,{\varvec{\zeta }}^\alpha \Bigr ) \Bigl ( \,\sum _{l\ge 0} \varepsilon ^{l-1} c_l^{k_2} \frac{{\mathrm d}^l}{{\mathrm d}t^l} z^{k_2} \Bigr ) \\&+ \displaystyle \sum _{k_1+k_2 = k} \Bigl ( \sum _{\begin{array}{c} m\ge 0 \\ s(\alpha ) = k_1 \end{array}} \frac{1}{m!} \Psi _{\widehat{B}}^{(m)} (z^0,t) \,{\mathbf{z}}^\alpha \Bigr ) \Bigl ( \sum _{\begin{array}{c} m\ge 0 \\ s(\alpha ) = k_2 \end{array}} \frac{1}{m!} E^{(m)} (z^0,t) \,{\mathbf{z}}^\alpha \Bigr ) , \end{aligned}$$

where \(T_{\widehat{B}}^{(m)} (x,t)\) denotes the mth derivative of \(T_{\widehat{B}} (x,t) = \frac{2}{h} \tanh \bigl ( -\frac{h}{2} \widehat{B} (x,t)\bigr )\) with respect to x and, similarly, \(\Psi _{\widehat{B}}^{(m)} (x,t)\) is the mth derivative of \(\Psi _{\widehat{B}} (x,t) = \Psi \bigl ( -h \widehat{B} (x,t)\bigr )\) with respect to x. These derivatives are bounded under the assumption that \(\eta = h/\varepsilon \le c\) and the non-resonance condition (3.2).

For \(k=0\) we obtain

$$\begin{aligned} \ddot{z}^0= & {} \displaystyle T_{\widehat{B}}(\zeta ^0,t) \dot{z}^0 + \Psi _{\widehat{B}} (z^0,t) E (z^0,t) \nonumber \\&+ \displaystyle 2\, {\mathrm {Re}\,}\Bigl ( \bigl (T_{\widehat{B}}' (\zeta ^0,t) \zeta ^{-1} \bigr )\frac{\mathrm {i}}{\varepsilon \eta } \sin (\eta \dot{\phi }) z^1 \Bigr ) + {\mathcal O}(\varepsilon ^2 ) , \end{aligned}$$
(5.12)

and for \(k=\pm 1 \) we get

$$\begin{aligned} \varepsilon ^{-2} d_0^{\pm 1} z^{\pm 1} + \varepsilon ^{-1} d_1^{\pm 1} \dot{z}^{\pm 1} = T_{\widehat{B}}(\zeta ^0 ,t) \bigl ( \varepsilon ^{-1} c_0^{\pm 1} z^{\pm 1} + c_1^{\pm 1} \dot{z}^{\pm 1} \bigr ) + {\mathcal O}(\varepsilon ) . \end{aligned}$$
(5.13)

Because of (5.11), the argument \(\zeta ^0\) can be replaced by \(z^0\) in these equations. In the limit \(h\rightarrow 0\), i.e., \(\eta \rightarrow 0\) we have accordance with the Eqs. (4.11) and (4.12) for the exact solution, respectively.

To get the differential equations for the dominant coefficient functions, we shall use the relations

$$\begin{aligned} \begin{array}{l} \displaystyle P_0(z^0,t) T_{\widehat{B}}(z^0,t) = 0, \qquad P_{\pm 1}(z^0,t) T_{\widehat{B}}(z^0,t) = \pm \, \mathrm {i}\frac{2}{h} \tan \Bigl ( \frac{h \dot{\phi }}{2\varepsilon } \Bigr ) P_{\pm 1}(z^0,t) , \\ \displaystyle P_0(z^0,t) \Psi _{\widehat{B}}(z^0,t) = P_0(z^0,t), \qquad P_{\pm 1}(z^0,t) \Psi _{\widehat{B}}(z^0,t) = \Psi \bigl (\pm \mathrm {i}\frac{h \dot{\phi }}{\varepsilon } \Bigr ) P_{\pm 1}(z^0,t) . \end{array} \end{aligned}$$

Multiplying the Eq. (5.12) with \(P_0(z^0,t)\) and applying the differentiation formula of Lemma 5.2 yields the differential equation

$$\begin{aligned} P_0(z^0,t) \,\ddot{z}_0^0= & {} P_0(z^0,t) E(z^0,t) \\&+ 2\, P_0(z^0,t) \, {\mathrm {Re}\,}\biggl ( -\frac{2}{\eta \dot{\phi }}\tan \Bigl (\frac{\eta \dot{\phi }}{2}\Bigr ) \Bigl ( \widehat{B}'(z^0,t)\zeta _{-1}^{-1}\Bigr )\, \mathrm {i}\frac{\dot{\phi }}{\varepsilon }{{\,\mathrm{sinc}\,}}(\eta \dot{\phi }) z_1^1 \biggr ) + {\mathcal O}(\varepsilon ^2 ) . \end{aligned}$$

Using \(\bigl ( \widehat{B} ' (x,t) \Delta x\bigr ) v = - v\times B ' (x,t) \Delta x\), which follows from differentiation of \(\widehat{B}(x,t) v = - v\times B(x,t)\), the trigonometric identity \(\sin (2\alpha ) = 2 \sin (\alpha ) \cos (\alpha )\), and the second relation of (5.11), this equation becomes (5.4).

A multiplication of (5.12) with \(P_{\pm 1} (z^0,t)\) permits to extract the dominant first derivative \(\dot{z}_{\pm 1}^0\) and gives (5.5).

We next consider the Eq. (5.13). The \(\varepsilon ^{-2}\)-terms in the left and right sides are contained in

$$\begin{aligned}&- \frac{4}{\varepsilon ^2 \eta ^2} \sin ^2 \Bigl ( \frac{\eta \dot{\phi }}{2}\Bigr )z^{\pm 1}\quad \hbox {and}\quad \pm \frac{2}{\varepsilon h} \tanh \Bigl ( - \frac{h}{2} \widehat{B}(z^0,t)\Bigr ) \frac{\mathrm {i}}{\eta }\sin (\eta \dot{\phi }) z^{\pm 1} . \end{aligned}$$

After multiplication with \(P_{\pm 1 }(z^0,t)\) these terms cancel because of the above formula for \(P_{\pm 1}(z^0,t) T_{\widehat{B}}(z^0,t)\). The remaining terms lead to

$$\begin{aligned} \dot{z}_{\pm 1}^{\pm 1} = -\frac{\displaystyle \Bigl ( \cos (\eta \dot{\phi }) + \tan \Bigl ( \frac{\eta \dot{\phi }}{2} \Bigr ) \sin (\eta \dot{\phi }) \Bigr ) \ddot{\phi }}{\displaystyle \frac{2}{\eta }\Bigl ( \sin (\eta \dot{\phi }) - \tan \Bigl ( \frac{\eta \dot{\phi }}{2} \Bigr ) \cos (\eta \dot{\phi })\Bigr )} z_{\pm 1}^{\pm 1} + {\mathcal O}(\varepsilon ^2 ) , \end{aligned}$$

which simplifies to (5.6). \(\square \)

In the above proof we referred to the following lemmas.

Lemma 5.1

([7]) For smooth functions \(\phi (t)\) and \(z^k(t)\) let \(y^k(t) = {\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon } z^k(t)\), and denote \(\eta = h/\varepsilon \). The finite differences of \(y^k(t)\) then satisfy

$$\begin{aligned} \delta _{2h}y^k(t)= & {} \displaystyle \frac{y^k(t+h) - y^k(t-h)}{2h} = {\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon } \sum _{l\ge 0} \varepsilon ^{l-1} c_l^k(t) \frac{{\mathrm d}^l}{{\mathrm d}t^l} z^k(t)\\ \delta _h^2 y^k(t)= & {} \displaystyle \frac{y^k(t+h) - 2y^k(t) + y^k(t-h)}{h^2} = {\mathrm e}^{\mathrm {i}k \phi (t)/\varepsilon } \sum _{l\ge 0} \varepsilon ^{l-2} d_l^k(t) \frac{{\mathrm d}^l}{{\mathrm d}t^l} z^k(t) , \end{aligned}$$

where \(c_{2j}^0 =0\), \(c_{2j+1}^0 = \eta ^{2j} /(2j+1)!\), and \(d_0^0=0\), \(d_{2j}^0 = 2\eta ^{2j-2}/(2j)!\), \(d_{2j+1}^0 = 0\). The leading coefficients are

$$\begin{aligned} c_0^k(t)= & {} \displaystyle \frac{\mathrm {i}}{\eta } \sin \bigl ( k\eta \dot{\phi }(t)\bigr ) - \varepsilon \frac{k\eta }{2} \sin \bigl ( k\eta \dot{\phi }(t)\bigr ) \ddot{\phi }(t)+ {\mathcal O}(\varepsilon ^2 ) \nonumber \\ c_1^k(t)= & {} \displaystyle \cos \bigl ( k\eta \dot{\phi }(t)\bigr ) + {\mathcal O}(\varepsilon ) \nonumber \\ d_0^k(t)= & {} \displaystyle -\frac{4}{\eta ^2} \sin ^2 \Bigl ( \frac{k\eta \dot{\phi }(t)}{2}\Bigr ) + \mathrm {i}\, \varepsilon \, k \cos \bigl ( k\eta \dot{\phi }(t)\bigr ) \ddot{\phi }(t) + {\mathcal O}(\varepsilon ^2 ) \nonumber \\ d_1^k(t)= & {} \displaystyle \frac{2\,\mathrm {i}}{\eta }\sin \bigl ( k\eta \dot{\phi }(t)\bigr ) + {\mathcal O}(\varepsilon ) . \end{aligned}$$
(5.14)

Note that these coefficients depend on \(\eta \), \(\varepsilon \), and t via derivatives of \(\phi (t)\).

Proof

Expanding \(\phi (t\pm h)\) and \(z^k(t\pm h)\) into Taylor series around t yields the stated formulas. \(\square \)

Lemma 5.2

Let \(T_{\widehat{B}} (x,t) = \frac{2}{h} \tanh \bigl ( -\frac{h}{2} \widehat{B} (x,t)\bigr )\), and let \(P_j(x,t)\) be the orthogonal projections onto the eigenspace of \(\widehat{B} (x,t)\) corresponding to the eigenvalues \(\lambda _0 =0\) and \(\lambda _1 = - \mathrm {i}|B(x,t)| = - \mathrm {i}\dot{\phi }(x,t)/\varepsilon \), and \(\lambda _{-1} = \mathrm {i}|B(x,t)| = \mathrm {i}\dot{\phi }(x,t)/\varepsilon \), respectively. Omitting the argument (x, t), we then have with \(\eta = h/\varepsilon \),

$$\begin{aligned} P_0 \Bigl ( T_{\widehat{B}}' \Delta x \Bigr ) P_{\pm 1} = \mp \frac{2}{\eta \dot{\phi }} \tan \Bigl ( \frac{\eta \dot{\phi }}{2} \Bigr ) P_0 \Bigl ( \widehat{B}' \Delta x \Bigr ) P_{\pm 1} , \end{aligned}$$

where prime indicates the derivative with respect to x.

Proof

Writing \(\tanh \) as a Taylor series with coefficients \(\gamma _l\) and differentiating term by term, we obtain

$$\begin{aligned} P_0 \Bigl ( T_{\widehat{B}} ' \Delta x \Bigr ) P_{\pm 1}= & {} \displaystyle \frac{2}{h} \sum _{l\ge 1} \gamma _l \Bigl ( - \frac{h}{2} \Bigr )^l P_0 \Bigl ( \widehat{B} ' \Delta x \Bigr ) \widehat{B}^{l-1} P_{\pm 1} \\= & {} \displaystyle \frac{2}{h} \sum _{l\ge 1} \gamma _l \Bigl ( - \frac{h}{2} \Bigr )^l P_0 \Bigl ( \widehat{B} ' \Delta x \Bigr ) \Bigl (\mp \mathrm {i}\, \frac{\dot{\phi }}{\varepsilon }\Bigr )^{l-1}P_{\pm 1}\\= & {} \displaystyle \frac{2\mathrm {i}}{\eta \dot{\phi }} \tanh \Bigl ( \pm \mathrm {i}\frac{\eta \dot{\phi }}{2} \Bigr ) P_0 \Bigl ( \widehat{B} ' \Delta x \Bigr ) P_{\pm 1} . \end{aligned}$$

This proves the statement of the lemma. \(\square \)

6 Proof of Theorem 3.1

Theorems 4.1 and 5.1 show that the coefficient functions \(z^k(t)\) [and also \(\dot{z}^0(t)\)] of the modulated Fourier expansions of the exact and numerical solutions coincide up to \({\mathcal O}(\varepsilon ^2)\) for the choice (2.5) and \(\,\Psi (\zeta )=\mathrm {tanch}(\zeta /2)\). This also shows that the phase functions \(\phi \) (with \(\dot{\phi }(t)= \varepsilon |B(z^0(t),t)|\)) differ only by \({\mathcal O}(\varepsilon ^2)\), respectively. Since all coefficient functions \(z^k\) of the modulated Fourier expansion with the exception of \(z^0\) are of size \({\mathcal O}(\varepsilon )\) or smaller, this yields that all summands \(z^k(t){\mathrm e}^{\mathrm {i}k\phi (t)/\varepsilon }\) still differ only by \({\mathcal O}(\varepsilon ^2)\). So we obtain the \({\mathcal O}(\varepsilon ^2)\) error bound for the positions as stated in Theorem 3.1.

We now turn to the error bound for the velocities. By Theorem 4.1, using that \(\dot{z}^{\pm 1}_{\pm 1}={\mathcal O}(\varepsilon ^2)\) and \(z^k_j={\mathcal O}(\varepsilon ^3)\) for \(|k|=1\) and \(k\ne j\) and for \(|k|\ge 2\) and all \(j=-1,0,1\), together with their derivatives, the velocity of the exact solution satisfies

$$\begin{aligned} v(t)=\dot{x}(t) = \dot{z}^0(t) + \frac{\mathrm {i}\dot{\phi }(t)}{\varepsilon }\, z^1_1(t) \, {\mathrm e}^{\mathrm {i}\phi (t)/\varepsilon } - \frac{\mathrm {i}\dot{\phi }(t)}{\varepsilon }\, z^{-1}_{-1}(t) \, {\mathrm e}^{-\mathrm {i}\phi (t)/\varepsilon } + {\mathcal O}(\varepsilon ^2). \end{aligned}$$
(6.1)

We shall show below that the numerical solution admits the same expansion with functions \(\phi (t), \dot{z}^0(t)\), \(z^1_1(t)\), \(z^{-1}_{-1}(t)\) that correspond to the modulated Fourier expansion (5.1) of the numerical solution and not to (4.1) of the exact solution. By Theorems 4.1 and 5.1, these functions differ only by \({\mathcal O}(\varepsilon ^2)\). Because of the denominator \(\varepsilon \) in the second and third terms on the right-hand side of (6.1), this yields

$$\begin{aligned} v^n - v(t^n)= {\mathcal O}(\varepsilon ),\qquad \text {but}\qquad v^n_\parallel - v_\parallel (t^n)= {\mathcal O}(\varepsilon ^2), \end{aligned}$$
(6.2)

and proves the statement of Theorem 3.1.

Using Lemma 5.1 and \(\ddot{\phi }(t)={\mathcal O}(\varepsilon )\), we have, with \(t=nh\), that

$$\begin{aligned} \frac{x^{n+1} - x^{n-1}}{2h} = \dot{z}^0 (t) + {{\,\mathrm{sinc}\,}}\bigl ( \eta \dot{\phi }(t)\bigr ) \frac{\mathrm {i}\dot{\phi }(t)}{\varepsilon } \Bigl ( z_1^1(t) {\mathrm e}^{\mathrm {i}\phi (t)/\varepsilon } - z_{-1}^{-1}(t) {\mathrm e}^{-\mathrm {i}\phi (t)/\varepsilon } \Bigr ) + {\mathcal O}(\varepsilon ^2 ) . \end{aligned}$$

A consequence of the maximal ordering in (1.1) is that \(\Phi _1\bigl ( h \widehat{B}(\bar{x}^n,t^n)\bigr ) = \Phi _1\bigl ( h \widehat{B}(z^0(t^n),t^n)\bigr ) + {\mathcal O}(\varepsilon ^2)\), and \(\Upsilon \bigl ( h \widehat{B}( x^n,t^n)\bigr ) = \Upsilon \bigl ( h \widehat{B}(z^0(t^n),t^n)\bigr ) + {\mathcal O}(\varepsilon ^2)\). Splitting \(\Phi _1(\cdot ) \dot{z}^0\) into \(\dot{z}^0 + \bigl (\Phi _1(\cdot ) - I \bigr )\dot{z}^0\) and using \(\Upsilon (\zeta ) = \bigl ( \Phi _1 (\zeta ) -1\bigr )/\zeta \), we therefore have

$$\begin{aligned} v^n= & {} \Phi _1 \bigl ( h\widehat{B}(\bar{x}^n,t^n)\bigr ) \frac{x^{n+1} -x^{n-1}}{2h} - h\Upsilon \bigl (h\widehat{B}( x^n,t^n)\bigr ) E(x^n,t^n)\\= & {} \dot{z}^0(t) + \frac{\mathrm {i}\dot{\phi }(t)}{\varepsilon } \Bigl ( z_1^1(t) {\mathrm e}^{\mathrm {i}\phi (t)/\varepsilon } - z_{-1}^{-1}(t) {\mathrm e}^{-\mathrm {i}\phi (t)/\varepsilon } \Bigr )\\&+ h\Upsilon \bigl ( h \widehat{B}(z^0(t),t)\bigr ) \Bigl (\widehat{B}\bigl (z^0(t),t\bigr ) \dot{z}^0(t) -E\bigl (z^0(t),t\bigr )\Bigr ) + {\mathcal O}(\varepsilon ^2 ) . \end{aligned}$$

Since \(\Upsilon (0)=0\) we have \(\Upsilon \bigl ( h \widehat{B}(z^0(t),t)\bigr )P_0(z^0(t),t) = 0\). On the other hand

$$\begin{aligned} P_{\pm 1}(z^0(t),t) \Bigl (\widehat{B}\bigl (z^0(t),t\bigr ) \dot{z}^0(t) -E\bigl (z^0(t),t\bigr )\Bigr ) = {\mathcal O}(\varepsilon ^2 ) \end{aligned}$$

which follows from (5.5) for \(\Psi (\mathrm {i}y) = \mathrm {tanc} (y/2)\). This proves the relation (6.1) also for the numerical solution.

7 A two-point filtered Boris algorithm

Algorithm 2.1 evaluates the magnetic field B at \(\bar{x}^n\) given by (2.4)–(2.5), which can be far from both \(x^n\) and the guiding center approximation \(x^n_\odot \) of (2.3) when \(h|B(x_n)|\) is close to a nonzero integral multiple of \(2\pi \). In the following we propose an alternative filtered Boris algorithm with the same second-order convergence properties as in Theorem 3.1, which evaluates the magnetic field at the two points \(x^n\) and \(x^n_\odot \).

Algorithm 7.1

(Two-point filtered Boris algorithm) Given \((x^n , v^{n-1/2})\), the algorithm computes \((x^{n+1} , v^{n+1/2})\) as follows, with \(B^n=B(x^n,t^n)\), \(B^n_\odot =B(x^n_\odot ,t^n)\) and \(E^n=E(x^n,t^n)\):

$$\begin{aligned} \begin{aligned} v^{n-1/2}_+&= v^{n-1/2} + \frac{h}{2}\, \Psi (h\widehat{B}^n)\,E^n \\ \Phi _2(h\widehat{B}^n_\odot ) (v^{n+1/2}_- - v^{n-1/2}_+)&= \frac{h}{2}\, \Phi _1(h\widehat{B}^n)\bigl (v^{n+1/2}_- + v^{n-1/2}_+\bigr ) \times B^n\\ v^{n+1/2}&= v^{n+1/2}_- + \frac{h}{2}\, \Psi (h\widehat{B}^n)\,E^n \\ x^{n+1}&= x^n + h\, v^{n+1/2},\\ \end{aligned} \end{aligned}$$
(7.1)

where \(\Psi (\zeta ) = {{\,\mathrm{tanch}\,}}(\zeta /2)\) and \(\displaystyle \Phi _1(\zeta ) = \frac{1}{{{\,\mathrm{sinch}\,}}(\zeta )}\) are as in Algorithm 2.1, and \(\displaystyle \Phi _2 (\zeta ) = \frac{1}{{{\,\mathrm{sinch}\,}}(\zeta /2)^2} \).

The velocity approximation \(v^n\) is again computed by (2.7), with \(B^n\) instead of \(\bar{B}^n\).

For constant B, Algorithms 2.1 and 7.1 are identical and explicit. In the general case, both methods are implicit, but this time the fixed-point iteration for \(x^n_\odot \) requires not only the evaluation of matrix functions by the Rodriguez formula, but in addition the solution of a linear system with the \(3\times 3\) matrix \( \Phi _2(h\widehat{B}^n_\odot )+ \frac{1}{2} h \widehat{B}^n \Phi _1(h\widehat{B}^n)\). We further note that in the case of a vanishing electric field, \(E^n=0\), Algorithm 2.1 preserves the velocity norm \(|v^{n+1/2}|=|v^{n-1/2}|\), which is satisfied only approximately up to \({\mathcal O}(h\varepsilon )\) by Algorithm 7.1. While these properties are unfavourable for Algorithm 7.1, our numerical experiments indicate that it yields higher accuracy than Algorithm 2.1 for stepsizes such that h|B| is large, and in particular it is less sensitive to near-resonances where h|B| is close to an integral multiple of \(2\pi \).

The two-step formulation of Algorithm 7.1 is

$$\begin{aligned}&\frac{x^{n+1} - 2x^n + x^{n-1}}{h^2} \nonumber \\&= \,\Phi _2( h\widehat{B}(\bar{x}^n,t^n))^{-1}\Bigl (\Phi _1( h\widehat{B}^n) \frac{x^{n+1} - x^{n-1}}{2h}\times B^n\Bigr )+ \Psi ( h \widehat{B}^n ) E^n . \end{aligned}$$
(7.2)

The starting value\(v^{1/2}\) is chosen such that formulas (7.1) and (2.7) also hold for \(n=0\). With the abbreviations

$$\begin{aligned} \Lambda ^n= & {} \Phi _2( h\widehat{B}^n_\odot )^{-1}\Phi _1( h\widehat{B}^n), \ \Psi ^n = \Psi (h\widehat{B}^n ),\ \Upsilon ^n=\Upsilon (h\widehat{B}^n), \\ \Phi ^n_\pm= & {} ( I \mp \tfrac{1}{2} \, \Lambda ^n h\widehat{B}^n ) {{\,\mathrm{sinch}\,}}(h\widehat{B}^n) , \\ \Psi ^n_\pm= & {} \Psi ^n \pm 2\Phi ^n_\pm \Upsilon ^n, \end{aligned}$$

we find, for arbitrary n, that like in (2.10),

$$\begin{aligned} v^{n\pm 1/2} = \Phi ^n_\pm v^n \pm \tfrac{h}{2} \,\Psi ^n_\pm E^n, \end{aligned}$$

and the one-step map \({(x^n,v^n)\mapsto (x^{n+1},v^{n+1})}\) is then again given by (2.11) with these modified matrices \(\Phi ^n_\pm \) and \(\Psi ^n_\pm \).

For the two-point filtered Boris algorithm, the second-order convergence result of Theorem 3.1 in x and \(v_\parallel \) and the first-order convergence in \(v_\perp \) remain valid, as can be shown by an adaptation of the proof of Theorem 5.1, for which we omit the details.

8 Numerical experiment

As an illustrative numerical experiment, we consider the charged-particle motion in the magnetic field

$$\begin{aligned} B(x,t) =\nabla \times \frac{1}{\varepsilon }\, \left( \begin{array}{c} 0 \\ x_1\\ 0 \\ \end{array} \right) +\nabla \times \, \left( \begin{array}{c} 0 \\ x_1 x_3 \\ 0 \\ \end{array} \right) =\frac{1}{\varepsilon }\, \left( \begin{array}{c} 0\\ 0 \\ 1 \\ \end{array} \right) + \left( \begin{array}{c} - x_1 \\ 0 \\ x_3 \\ \end{array} \right) , \end{aligned}$$

and the electric field \(E(x,t)=-\nabla _x U(x)\) with the potential

$$\begin{aligned} U(x)=\frac{1}{\sqrt{x_1^2+x_2^2}}. \end{aligned}$$

The initial values are chosen as \(x(0)=(\frac{1}{3},\frac{1}{4},\frac{1}{2})^{\intercal }\) and \(v(0)=(\frac{2}{5},\frac{2}{3},1)^{\intercal }\). We solve this problem for \(0\le t \le 1\) with \(h=\epsilon ,4\epsilon ,16\epsilon \) and compare the numerical errors of the following methods:

  • the standard Boris algorithm,

  • Exp-A: the filtered Boris method of Algorithm 2.1 with \(\theta =1\) in (2.4) (where \(\bar{x}^n=x^n\) and the method is explicit),

  • Imp-A: the filtered Boris method of Algorithm 2.1 with \(\theta \) of (2.5),

  • Two P-A: the two-point filtered Boris method of Algorithm 7.1.

Fig. 1
figure 1

The logarithm of the global error against the logarithm of \(\epsilon \)

Full size image
Fig. 2
figure 2

The logarithm of the global error at \(t = 1\) against \(h/\epsilon \) for \(\epsilon =1/2^{10}\) and \(h=1/k\), where \(k=60,61,\ldots ,600\)

Full size image

The errors in x and \( v_\parallel , v_\perp \) against different \(\epsilon =1/2^j\) are displayed in Fig. 1, where \(j=4,\ldots ,13\). Then we fix \(\epsilon =1/2^{10}\) and show the errors at \(t=1\) against \(h/\epsilon \) in Fig. 2. It is observed that all three filtered Boris methods improve considerably over the standard Boris method, and the optimally filtered methods Imp-A and Two P-A show second order, whereas method Exp-A only shows first order. Methods Imp-A and Two P-A behave very similar away from stepsize resonances, but method Two P-A appears more robust near stepsize resonances. For the implicit methods Imp-A and Two P-A, the error behaviour remains essentially unchanged after just one fixed-point iteration.