An energy-conserving and asymptotic-preserving charged-particle orbit implicit time integrator for arbitrary electromagnetic fields

TL;DR

This paper introduces an energy-conserving, asymptotic-preserving implicit integrator for charged-particle orbits in electromagnetic fields, capable of accurately capturing full or guiding center motion with adaptive time-stepping.

Contribution

The paper presents a novel implicit scheme based on Crank-Nicolson that conserves energy and handles multiple physical regimes with adaptive time-stepping.

Findings

Significantly improved accuracy due to energy conservation.

Smooth transition between magnetized and unmagnetized regimes.

Effective over gyration time-scale with adaptive stepping.

Abstract

We present a new implicit asymptotic preserving time integration scheme for charged-particle orbit computation in arbitrary electromagnetic fields. The scheme is built on the Crank-Nicolson integrator and continues to recover full-orbit motion in the small time-step limit, but also recovers all the first-order guiding center drifts as well as the correct gyroradius when stepping over the gyration time-scale. In contrast to previous efforts in this direction, the new scheme also features exact energy conservation. In the derivation of the scheme, we find that a new numerical time-scale is introduced. This scale is analyzed and the resulting restrictions on time-step are derived. Based on this analysis, we develop an adaptive time-stepping strategy the respects these constraints while stepping over the gyration scale when physically justified. It is shown through numerical tests on single-particle motion that the scheme's energy conservation property results in tremendous improvements in accuracy, and that the scheme is able to transition smoothly between magnetized and unmagnetized regimes as a result of the adaptive time-stepping.

Tables

Table 1. Table 1. The maximum allowable value of Ω c ​ Δ ​ t subscript Ω 𝑐 Δ 𝑡 \Omega_{c}\Delta t to ensure (a) the anomalous displacement does not exceed gyroradius and (b) the time-scale over which accurate implicit gyroaveraging is achieved is shorter than τ res subscript 𝜏 res \tau_{\textrm{res}} . Expressions are valid in the asymptotic limits (i) u n + 1 / 2 ≫ v E much-greater-than superscript 𝑢 𝑛 1 2 subscript 𝑣 𝐸 u^{n+1/2}\gg v_{E} and (ii) u n + 1 / 2 ≪ v E much-less-than superscript 𝑢 𝑛 1 2 subscript 𝑣 𝐸 u^{n+1/2}\ll v_{E} , respectively.

$\max \left({\Omega }_c\Delta t\right)$	(a) Displacement restriction	(b) Avg. restriction
(i) $u^{n+1/2}\gg v_E$	$2\min \{\sqrt{2}{\left({\delta }_{⟂}\right)}^{-1/2},{\left({\delta }_{\parallel }\right)}^{-1}\}$	$2{\left(\frac{{\Omega }_c{\tau }_{\text{res}}}{\pi }\right)}^{1/2}$
(ii) $u^{n+1/2}\ll v_E$	$\sqrt{2}{\left({\delta }_E+{\delta }_{\parallel }\right)}^{-1/2}$	${\Omega }_c{\tau }_{\text{res}}$

Table 2. Table 2. Mean and standard deviation in final μ 𝜇 \mu value for the third test case.

	Boris (${\Omega }_c\Delta t\ll 1$)	CBFV	BFV	MI
Mean $\mu$	0.099487	0.100280	0.026340	0.0172801
Std. dev. in $\mu$	0.051293	0.058710	0.025910	0.014327

Table 3. Table 3. Critical values of v y 0 superscript subscript 𝑣 𝑦 0 v_{y}^{0} that separate trapped and passing regimes found by each scheme. Boris is taken as the reference case. Values for MI are not reported due to the unphysical character of the orbits.

	Boris (${\Omega }_c\Delta t\ll 1$)	CBFV	BFV
Critical $v_y^0$	0.73107	0.73110	0.62725
% Error	$-$	0.004	14.2

Equations

$$\begin{array}{rl}{\mathbf{v}}_p^{n+1}& ={\mathbf{v}}_p^n+\Delta t\frac{q_p}{m_p}\left({\mathbf{E}}_p^{n+1/2}+{\mathbf{v}}_p^{n+1/2}\times {\mathbf{B}}^{n+1/2}\right),\\ {\mathbf{x}}_p^{n+1}& ={\mathbf{x}}_p^n+\Delta t{\mathbf{v}}_p^{n+1/2}.\end{array}$$

$$\begin{array}{rl}{\mathbf{v}}_p^{n+1/2}& =\left({\mathbf{v}}_p^n+{\mathbf{v}}_p^{n+1}\right)/2,{\mathbf{x}}_p^{n+1/2}=\left({\mathbf{x}}_p^n+{\mathbf{x}}_p^{n+1}\right),\\ {\mathbf{E}}_p^{n+1/2}& =\mathbf{E}\left({\mathbf{x}}_p^{n+1/2},t^{n+1/2}\right),{\mathbf{B}}_p^{n+1/2}=\mathbf{B}\left({\mathbf{x}}_p^{n+1/2},t^{n+1/2}\right).\end{array}$$

$$\begin{array}{rl}\sum _p\frac{m_p}{2}{\Vert {\mathbf{v}}_p^{n+1}\Vert }^2-\sum _p\frac{m_p}{2}{\Vert {\mathbf{v}}_p^n\Vert }^2& =\frac{1}{2}\sum _p{m}_p{\mathbf{v}}_p^{n+1/2}\cdot \left({\mathbf{v}}_p^{n+1}-{\mathbf{v}}_p^n\right)\\ & =\frac{\Delta t}{2}\sum _p{q}_p{\mathbf{v}}_p^{n+1/2}\cdot {\mathbf{E}}^{n+1/2}.\end{array}$$

$$\begin{array}{rl}\mathbf{v}& \approx \mathbf{u}+v_{\parallel }\mathbf{b}+{\mathbf{v}}_E+{\mathbf{v}}_I+{\mathbf{v}}_{\nabla B},\\ \mathbf{u}& =u_{\perp }(\cos {\Omega }_{c}t{\mathbf{e}}_1+\sin {\Omega }_{c}t{\mathbf{e}}_2),\\ {\mathbf{v}}_E& =\frac{\mathbf{E}\times \mathbf{B}}{B^2},\\ {\mathbf{v}}_I& ={\mathbf{v}}_P+{\mathbf{v}}_C,\\ {\mathbf{v}}_P& =\frac{\mathbf{b}}{{\Omega }_c}\times \frac{d{\mathbf{v}}_E}{dt},\\ {\mathbf{v}}_C& =\frac{\mathbf{b}}{{\Omega }_c}\times v_{\parallel }\frac{d\mathbf{b}}{dt},\\ {\mathbf{v}}_{\nabla B}& =\frac{\mathbf{b}}{{\Omega }_c}\times \frac{\mu }{m}\nabla B.\end{array}$$

$$m\frac{d{v}_{\parallel }}{dt}=q{E}_{\parallel }-\mu \nabla B\cdot \mathbf{b}-m\frac{d{\mathbf{v}}_E}{dt}\cdot \mathbf{b}$$

$${\tilde{\mathbf{v}}}^n={\mathbf{v}}^n-{\mathbf{v}}_E.$$

$${\tilde{\mathbf{v}}}^{n+1/2}={\tilde{\mathbf{v}}}^n+\frac{\Delta t}{2}\frac{q}{m}\left({\mathbf{E}}_{\parallel }+{\tilde{\mathbf{v}}}^{n+1/2}\times \mathbf{B}\right).$$

$${\tilde{\mathbf{v}}}^{n+1/2}\times \mathbf{b}={\tilde{\mathbf{v}}}^n\times \mathbf{b}-\frac{\Delta t}{2}\frac{q}{m}{\tilde{\mathbf{v}}}_{\perp }^{n+1/2}B.$$

$${\tilde{\mathbf{v}}}_{\perp }^{n+1/2}=\frac{{\tilde{\mathbf{v}}}_{\perp }^n+\frac{1}{2}{\Omega }_c\Delta t{\tilde{\mathbf{v}}}^n\times \mathbf{b}}{1+{\Omega }_c^2\Delta t^2/4}.$$

$${\tilde{\mathbf{v}}}_{\perp }^{n+1}=\frac{\left(1-{\Omega }_c^2\Delta t^2/4\right){\tilde{\mathbf{v}}}_{\perp }^n+{\Omega }_c\Delta t\left({\tilde{\mathbf{v}}}^n\times \mathbf{b}\right)}{1+{\Omega }_c^2\Delta t^2/4}.$$

$$\cos \theta =\frac{1-{\Omega }_c^2\Delta t^2/4}{1+{\Omega }_c^2\Delta t^2/4}.$$

$${\mathbf{v}}_{\perp }^{n+1}={\mathbf{v}}_E+R_{\theta }\left[{\mathbf{v}}_{\perp }^n-{\mathbf{v}}_E\right],$$

$$\left(\mathbf{I}-{\mathbf{b}}^{n+1/2}{\mathbf{b}}^{n+1/2}\right)\cdot {\mathbf{v}}^{n+1}={\mathbf{v}}_E^{n+1/2}+R_{\theta }^{n+1/2}\left[\left(\mathbf{I}-{\mathbf{b}}^{n+1/2}{\mathbf{b}}^{n+1/2}\right)\cdot {\mathbf{v}}^n-{\mathbf{v}}_E^{n+1/2}\right],$$

$$\begin{array}{rl}{\mathbf{x}}^{n+1/2}& ={\mathbf{x}}^{n-1/2}+\Delta t{\mathbf{v}}^n,\\ {\mathbf{v}}^{n+1}& ={\mathbf{v}}^n+\Delta t\frac{q}{m}\left({\mathbf{E}}^{n+1/2}+{\mathbf{v}}^{n+1/2}\times {\mathbf{B}}^{n+1/2}\right)\end{array}$$

$${\rho }_{\text{eff}}=\frac{v_{\perp }\Delta t}{2\left|\sin \left(\varphi /2\right)\right|}.$$

$${\rho }_{\text{eff}}^{\text{Boris}}={\rho }_{\text{true}}\sqrt{1+\frac{{\Omega }_c^2\Delta t^2}{4}},$$

$$v_{\perp }^{n+1/2}=\frac{v_{\perp }^n}{\sqrt{1+{\Omega }_c^2\Delta t^2/4}}.$$

$${\rho }_{\text{eff}}^{\text{CN}}={\rho }_{\text{true}}.$$

$${\mu }_{\text{eff}}=\frac{q{\rho }_{\text{eff}}^2}{2\Delta t}\sin \theta .$$

$${\mu }_{\text{eff}}^{\text{CN}}=\frac{{\mu }_{\text{true}}}{1+{\Omega }_c^2\Delta t^2/4}.$$

$${\mathbf{F}}_{BFV}^{n+1/2}=-\tilde{\mu }\nabla B,\tilde{\mu }=m\frac{{\Vert {\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_{\perp }^n\Vert }^2}{8B},$$

$$\begin{array}{rl}{\Vert {\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_{\perp }^n\Vert }^2& ={\Vert \left({\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_E^{n+1/2}\right)-\left({\mathbf{v}}_{\perp }^n-{\mathbf{v}}_E^{n+1/2}\right)\Vert }^2,\\ & ={\Vert \left[R_{\theta }^{n+1/2}-\mathbf{I}\right]\left({\mathbf{v}}_{\perp }^n-{\mathbf{v}}_E^{n+1/2}\right)\Vert }^2,\\ & ={\Vert {\mathbf{v}}_{\perp }^n-{\mathbf{v}}_E^{n+1/2}\Vert }^2\frac{{\Omega }_c^2\Delta t^2}{1+{\Omega }_c^2\Delta t^2/4}.\end{array}$$

$$\begin{array}{rl}{\Omega }_c\Delta t\ll 1& ⟹\tilde{\mu }=O(\Delta t^2),\\ {\Omega }_c\Delta t\gg 1& ⟹\tilde{\mu }\approx \frac{m{\Vert {\mathbf{v}}_{\perp }^n-{\mathbf{v}}_E^{n+1/2}\Vert }^2}{2B}.\end{array}$$

$$\begin{array}{rl}-\left({\mu }_{\text{eff}}^{\text{CN}}+\tilde{\mu }\right)\nabla B& =-{\mu }_{\text{true}}\nabla B\left\{\frac{1}{1+{\Omega }_c^2\Delta t^2/4}+\frac{{\Omega }_c^2\Delta t^2/4}{1+{\Omega }_c^2\Delta t^2/4}\right\}\\ & =-{\mu }_{\text{true}}\nabla B.\end{array}$$

$${\mathbf{F}}_{MI}^{n+1/2}=\frac{q}{4}\left({\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_{\perp }^n\right)\times \left({\mathbf{B}}^{n+1}-{\mathbf{B}}^n\right).$$

$${\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_{\perp }^n={\Omega }_c\Delta t\left({\mathbf{v}}_{\perp }^{n+1/2}-{\mathbf{v}}_E^{n+1/2}\right)\times {\mathbf{b}}^{n+1/2}.$$

$$\begin{array}{rl}{\mathbf{B}}^{n+1}-{\mathbf{B}}^n& =\left({\mathbf{x}}^{n+1}-{\mathbf{x}}^n\right)\cdot \nabla {\mathbf{B}}^{n+1/2}+O\left({\Vert {\mathbf{x}}^{n+1}-{\mathbf{x}}^n\Vert }^3\right)\\ & =\Delta t\left({\mathbf{v}}^{n+1/2}\cdot \nabla \right){\mathbf{B}}^{n+1/2}+O\left(\Delta t^3\right).\end{array}$$

$$\frac{1}{4}\left({\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_{\perp }^n\right)\times \left({\mathbf{B}}^{n+1}-{\mathbf{B}}^n\right)\approx \frac{1}{4}{\Omega }_c\Delta t^2\left[\left({\mathbf{v}}_{\perp }-{\mathbf{v}}_E\right)\times \mathbf{b}\right]\times \left(\mathbf{v}\cdot \nabla \right)\mathbf{B},$$

$$⟨\left[\left({\mathbf{v}}_{\perp }-{\mathbf{v}}_E\right)\times \mathbf{b}\right]\times \left(\mathbf{v}\cdot \nabla \right)\mathbf{B}⟩\approx ⟨\left({\tilde{\mathbf{v}}}_{\perp }\times {\mathbf{b}}_{gc}\right)\times \left({\tilde{\mathbf{v}}}_{\perp }\cdot \nabla \right){\mathbf{B}}_{gc}⟩,$$

$$\frac{q}{4}⟨\left({\mathbf{v}}_{\perp }^{n+1}-{\mathbf{v}}_{\perp }^n\right)\times \left({\mathbf{B}}^{n+1}-{\mathbf{B}}^n\right)⟩\approx -\frac{{\Omega }_c^2\Delta t^2}{4}\frac{q{\Vert {\tilde{\mathbf{v}}}_{\perp }^{n+1/2}\Vert }^2}{2{\Omega }_c}\nabla B^{n+1/2}.$$

Full text

1. Introduction

In recent years, considerable progress has been made on implicit particle-in-cell (PIC) methods - e.g. [5, 7, 8, 9, 10, 11, 25]. Prominent among the advantages of these methods is their exact total energy conservation for any time-step and cell size. Due to this conservation, implicit PIC schemes have been shown to be more robust against the finite-grid instability than their explicit counterparts. In addition to the obvious benefits of exact energy conservation for long-time accuracy, resistance to this instability facilitates enormous speed-ups when solution structures of interest are much larger than the Debye length.

However, in strongly magnetized plasmas, there is another onerous time-step constraint not overcome by these implicit schemes - namely, the time-step must resolve the gyroperiod. Traditionally, this difficulty has been overcome using gyrokinetic plasma models [6, 12, 14, 17, 21], in which the gyroperiod time-scale has been analytically eliminated by an asymptotic expansion.

This gyrokinetic approach presents two difficulties. Firstly, the exact energy conservation property enjoyed by implicit schemes is not easily carried over to the gyrokinetic context. Indeed, some gyrokinetic models do not conserve energy even at the continuum level, and those that do conserve only an approximation of the total full-orbit energy [3]. Secondly, in many problems of interest, the gyrokinetic approximation may be valid only in a subset of the spatial domain. If the transition region - between strongly and weakly magnetized portions of the domain - is narrow and its location is known a priori, then a domain decomposition method can be effective. However, in challenging problems of scientific interest, the transition region may be wide and its location may be unknown and/or time varying.

An alternate approach to overcoming the gyroperiod time-step restriction that side- steps these difficulties is to derive an asymptotic preserving (AP) time-stepping scheme. That is, one seeks a time-stepping scheme for particle evolution that (a) recovers the exact dynamics as $\Delta t\rightarrow 0$, (b) recovers the guiding-center particle motion when $\Omega_{c}\Delta t\gg 1$, where $\Omega_{c}=qB/m$ is the gyro-frequency, and (c) transitions seamlessly between these two regimes. Several such schemes exist in the literature [2, 13, 15, 16, 18, 26], but none feature exact total energy conservation in the context of a PIC scheme.

The purpose of the present article is to propose a new implicit, AP time-stepping method that exactly conserves total energy. The starting point is the Crank-Nicolson integrator used in prior implicit PIC work - see citations above. Building on asymptotic and empirical results in [2], we show that, in contrast to the classical Boris scheme (which is known to feature an artificially enlarged gyroradius for $\Omega_{c}\Delta t\gg 1$ [1, 22]), Crank-Nicolson recovers the correct gyroradius for arbitrary time-step. We then add modifications that borrow ideas from both [2, 26] and [18] to ensure that the scheme captures the $\nabla B$ drift while preserving energy conservation. We refer to prior work [22] to find that the scheme captures the inertial drifts for arbitrary time-step with no additional modifications.

Additionally, we derive expressions for the time-step restrictions on our new scheme. It is shown that when the assumptions necessary for drift motion are satisfied, these time- step limits do in fact allow for $\Omega_{c}\Delta t\gg 1$. We present an adaptive time-stepping scheme that ensures that these new restrictions are obeyed at each time step in spatiotemporally varying fields.

Finally, our scheme is compared against the standard Boris push and existing AP integrators. This is done for single particle motion in several field configurations of varying complexity. We find that the enforcement of energy conservation gives the new scheme dramatically improved long-time accuracy compared to previous efforts. For instance, while adiabatic invariance of magnetic moment is not explicitly enforced, we find that the new scheme conserves it to a good approximation in magnetized regions while previous schemes fail to do so.

The remainder of the article is structured as follows. In Section 2, we review background material, including the requirements for energy conservation in implicit PIC, guiding center motion, properties of Crank-Nicolson applied to charged particle motion, and existing AP integrators. In Section 3, we develop our energy-conserving modification to capture the $\nabla B$ drift; we also derive time-step restrictions for our scheme and present an adaptive time- stepping strategy. We present numerical examples in Section 4 and conclude in Section 5.

2. Review

2.1. Implicit PIC

At the heart of recent energy-conserving implicit PIC schemes is the Crank-Nicolson integrator for updating particle position and velocity at time $t_{n}=n\Delta t$, denoted by $(\mathbf{x}^{n}_{p},\mathbf{v}_{n}^{p})$, to time $t^{n+1}$:

[TABLE]

where the quantities at the half-time-step are defined by

[TABLE]

Here, $\mathbf{E}(\mathbf{x},t)$ and $\mathbf{B}(\mathbf{x},t)$ are computed self-consistently using the particle data as inputs.

Exact energy conservation is shown under various conditions by noting that

[TABLE]

From here, the details vary slightly depending on the context (e.g. electrostatic [9] or Vlasov-Darwin [8]), but the common theme is that if (a) $\mathbf{E}$ is updated using Ampere’s law rather than Poisson’s equation, and (b) the current density $\mathbf{j}$ is interpolated to the grid using the same shape function that is used to interpolate $\mathbf{E}$ from the grid to the particles, then one can choose the discretization such that the second line of (3) is equivalent to the negation of the change in discrete potential energy from step $n$ to $n+1$. Thus, the total energy is unchanged by a time-step. See [8, 9] for considerably more detail.

For our purposes, the key step in this analysis lies in going from the first to second line in (3). There, we have used the fact that the discrete magnetic-field force is necessarily orthogonal to $\mathbf{v}_{p}^{n+1/2}$. Any additional forces added to the first line of (1) that are also orthogonal to $\mathbf{v}_{p}^{n+1/2}$ will thus not break the energy conservation property. This observation is crucial in the derivation of our energy-conserving AP integrator, as it allows us to work without considering the details of the field solve and its relation to particle data. As such, for the remainder of the article we are safe in dropping the $p$ subscript everywhere and considering a single particle in prescribed electromagnetic fields.

2.2. Guiding Center Motion and Crank-Nicolson

Single particle motion in prescribed electromagnetic fields is well-understood in the asymptotic limit in which the gyroperiod $\Omega_{c}^{-1}$ is small compared to all other time scales. The standard derivation - see e.g. [20, 19] - uses the method of averaging to show that, to leading order, the particle velocity is composed of a rapid gyration $\mathbf{u}$ superimposed on several slower “drifts”:

[TABLE]

Above, $\mathbf{b}=\mathbf{B}/B$, $\{\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{b}\}$ form an orthonormal basis for $\mathbb{R}^{3}$, $\mu=mu_{\perp}^{2}/2B$ is the magnetic moment, and the evolution of $u_{\perp}$ is specified by conservation of $\mu$ - i.e. $\dot{\mu}=0$. The parallel velocity solves the ODE

[TABLE]

to leading order.

These drifts are often given names, which we will find convenient to use: $\mathbf{v}_{E}$ is called the “$\mathbf{E}\times\mathbf{B}$ drift”, $\mathbf{v}_{I}$ the “inertial drift”, and $\mathbf{v}_{\nabla B}$ the “magnetic drift” or “$\nabla B$ drift”. The inertial drift is often decomposed into $\mathbf{v}_{P}$ - the “polarization drift” - and $\mathbf{v}_{C}$, which we will refer to as the “curvature drift”. Traditionally, the curvature drift is the name for a particular component of $\mathbf{v}_{C}$, but it is frequently the most important one, so we will use this name for all of $\mathbf{v}_{C}$.

Additionally, the effective force $-\mu\nabla B\cdot\mathbf{b}$ appearing in (5) is often called the “mirror force”, so named because it gives rise to particle confinement in a magnetic mirror. In the course of the derivation, it has the same origin as the $\nabla B$ drift - namely, an effective force $-\mu\nabla B$ acts on the guiding center. The perpendicular component of this force results in the $\nabla B$ drift, while the parallel component manifests as the mirror force. For this reason, we shall use the term “$\nabla B$ drift” to refer both to the explicit drift $\mathbf{v}_{\nabla B}$ and motion induced by the mirror force when no confusion results.

It is our aim to find a time integration scheme that accurately reproduces all of these drifts when $\Omega_{c}\Delta t\gg 1$. It is straightforward to show that the Crank-Nicolson scheme (1) already captures the $\mathbf{E}\times\mathbf{B}$ drift for arbitrary time-step. However, we present the argument in detail here because it establishes notation that we will find useful in our later analysis. It will also lead us to the conclusion that Crank-Nicolson recovers the correct gyroradius for arbitrary $\Delta t$.

To that end, let us temporarily assume that $\mathbf{E}$ and $\mathbf{B}$ are uniform in time and space so that we may drop their $n+1/2$ superscripts. We define the new discrete velocity

[TABLE]

Elementary manipulation shows that

[TABLE]

We wish to focus on perpendicular motion, so we cross the above with $\mathbf{b}$ to find

[TABLE]

Substituting this into the perpendicular component of (7) gives

[TABLE]

From here, simple manipulation gives

[TABLE]

Dotting both sides with themselves and doing some algebra, we can see that $\left\|\tilde{\mathbf{v}}_{\perp}^{n+1}\right\|=\left\|\tilde{\mathbf{v}}^{n}_{\perp}\right\|$. This reveals that the transformation $\tilde{\mathbf{v}}^{n}\rightarrow\tilde{\mathbf{v}}^{n+1}$ is a rotation. Dotting with $\tilde{\mathbf{v}}^{n}$ shows that the rotation angle $\theta$ satisfies

[TABLE]

In constant, uniform fields, we can thus write the perpendicular component of the velocity update as

[TABLE]

where $R_{\theta}$ denotes the rotation matrix about $\mathbf{b}$ by angle $\theta$ given above - i.e. it can be defined by the relation $\tilde{\mathbf{v}}_{\perp}^{n+1}=R_{\theta}\tilde{\mathbf{v}}_{\perp}^{n}$ along with (10). In uniform fields, this would be the exact solution if it were the case that $\theta=\Omega_{c}\Delta t$. However, $\theta$ only approximates $\Omega_{c}\Delta t$ for $\Omega_{c}\Delta t\ll 1$, as can be easily seen from Taylor expansion of the right side of (11). It bears noting that for $\Omega_{c}\Delta t\gg 1$, the rotation angle $\theta$ tends to $\pi$. This observation lends intuition to later analysis.

Thus, we see that at least in constant, homogeneous fields, Crank-Nicolson exactly recovers the correct $\mathbf{E}\times\mathbf{B}$ drift, but makes errors in the gyro-phase of the particle that grow with time-step. However, in the strongly magnetized case, errors in gyro-phase are of little concern, since we assume that the gyroradius is small compared to length scales of interest.

Finally, we note that nothing in this derivation actually requires $\mathbf{E}$ and $\mathbf{B}$ to be spatially or temporally uniform - it was simply convenient to drop the superscripts. We may thus reintroduce the $n+1/2$ superscripts on all fields now. The only additional subtlety is that the $\perp$ component is always with respect to $\mathbf{b}^{n+1/2}$, regardless of the time-level at which we are evaluating the velocity - this is a consequence of (8), in which the $\mathbf{b}$ in the cross product is evaluated at $t^{n+1/2}$. The general form of (12) is thus

[TABLE]

where $\mathbf{I}$ is the identity tensor. Thus, even in varying fields, the $\mathbf{E}\times\mathbf{B}$ drift is approximately recovered so long as $\mathbf{b}$ and $\mathbf{v}_{E}$ change little within a time-step. Since each is assumed to vary on a scale much longer than the gyro-period, this is a much less restrictive constraint than $\Omega_{c}\Delta t\ll 1$.

2.3. Numerical Gyroradius

It is well known that the classical, explicit Boris time integrator, which may be defined by

[TABLE]

features a numerical gyroradius much larger than the true gyroradius for $\Omega_{c}\Delta t\gg 1$ [1, 22]. In particular, the velocity update is identical to the Crank-Nicolson scheme defined above, so in a homogeneous magnetic field with $\mathbf{E}=0$, the velocity update still represents rotation by the angle $\theta$ defined in (11). Geometric arguments show that when the position is updated using a velocity of constant magnitude $v_{\perp}$ rotating at each time step by an angle $\phi$, the radius of the circular motion is given by [1, 22] (see particularly Figure 4-3b in [1])

[TABLE]

In the Boris push, the perpendicular velocity used to update position is simply $\mathbf{v}_{\perp}^{n}$. Making use of the trigonometric identity $\left\lvert\sin\left(\phi/2\right)\right\rvert=\sqrt{(1-\cos\phi)/2}$ and (11), one finds that

[TABLE]

where $\rho_{\textrm{true}}=v_{\perp}^{n}/\Omega_{c}$ is the physically correct gyroradius. This clearly has negative consequences for accuracy in spatially varying fields when $\Omega_{c}\Delta t\gg 1$, as discussed in [26].

In contrast, when using the Crank-Nicolson scheme defined in (1), the velocity used in updating position is not $\mathbf{v}^{n}$ but $\mathbf{v}^{n+1/2}$. By returning to (9) and dotting $\mathbf{v}_{\perp}^{n+1/2}$ with itself (we may drop the tilde since we are assuming $\mathbf{E}=0$), we find

[TABLE]

This statement may initially be counter-intuitive, but may be explained in the following way: We have already established that for $\Omega_{c}\Delta t\gg 1$, the angle $\theta$ through with $\mathbf{v}_{\perp}$ rotates in a time step is nearly $\pi$, so that $\mathbf{v}_{\perp}^{n+1}\approx-\mathbf{v}_{\perp}^{n}$. Thus, $\mathbf{v}_{\perp}^{n+1/2}$, being an average of two vectors which are approximately the negations of each other, should be much smaller than either. Indeed, $v_{\perp}^{n+1/2}$ should tend to zero as $\Omega_{c}\Delta t$ diverges, which is precisely what (17) states.

Substituting this expression for $\mathbf{v}_{\perp}$ into (15), we find that when the Crank-Nicolson scheme is used, we have

[TABLE]

Thus, remarkably, the Crank-Nicolson scheme recovers the exact correct gyroradius for any value of $\Omega_{c}\Delta t$. This was shown in the limit $\Omega_{c}\Delta t\gg 1$ in [26], and observed empirically in the general case in the numerical experiments of [2, 26]. Thus, the implicit Crank-Nicolson scheme begins at an advantage over the Boris push when investigating particle dynamics with $\Omega_{c}\Delta t\gg 1$.

2.4. Preservation of Inertial Drifts and Non-Preservation of $\nabla B$ Drift

In [22], Parker and Birdsall show that the Boris push recovers all first order drift motions for arbitrary $\Omega_{c}\Delta t$, albeit with an artificially enlarged gyroradius. The discretization is shown to only introduce an additional force of $O(q\rho_{\textrm{eff}}^{2}\|D^{2}\mathbf{B}\|/m)$, where $\|D^{2}\mathbf{B}\|$ is a norm of the Hessian of $\mathbf{B}$. That is, if the effective gyroradius is small compared to the scale over which $\mathbf{B}$ varies, this force is small.

In Parker and Birdsall’s numerical guiding center equation derivation, the numerical $\nabla B$ drift arises from an effective force $-\mu_{\textrm{eff}}\nabla B$, where

[TABLE]

One may use the expressions above for $\theta$ and $\rho_{\textrm{eff}}^{\textrm{Boris}}$ to show that, for the Boris push, $\mu_{\textrm{eff}}^{\textrm{Boris}}=\mu_{\textrm{true}}=mu_{\perp}^{2}/2B$.

However, we have just shown that, for Crank-Nicolson, $\rho_{\textrm{eff}}$ is much smaller than for Boris when $\Omega_{c}\Delta t\gg 1$. Thus, Crank-Nicolson dramatically underestimates the $\nabla B$ drift for $\Omega_{c}\Delta t\gg 1$. To be precise, since $\rho_{\textrm{eff}}^{\textrm{CN}}=\rho_{\textrm{eff}}^{\textrm{Boris}}/\sqrt{1+\Omega_{c}^{2}\Delta t^{2}/4}$ - see (16) and (18) - and the rotation angles $\theta$ are identical, we have

[TABLE]

This, though, is the only modification of the Parker/Birdsall analysis. Thus, Crank-Nicolson still captures the inertial drifts - this is verified empirically in our numerical experiments in Section 4. For this reason, we follow the lead of previous efforts toward AP schemes by focusing on modifying Crank-Nicolson to capture the $\nabla B$ drift.

2.5. Prior AP Schemes

Previous works have succeeded in capturing the $\nabla B$ drift, although without the energy conservation property we seek. These schemes can be sorted into three categories. In the first category are schemes [2, 18, 26] that modify the velocity update equation by adding an “effective” force that approximates $-\mu\nabla B$ for $\Omega_{c}\Delta t\gg 1$ and is negligible for $\Omega_{c}\Delta t\ll 1$. In the second category are schemes [13, 18] that modify the position update equation by adding an “effective” velocity that approximates $\mathbf{v}_{\nabla B}$ for $\Omega_{c}\Delta t\gg 1$ and is negligible for $\Omega_{c}\Delta t\ll 1$. In the third category are the semi-implicit Runge-Kutta schemes of [16] - these schemes are restricted to magnetic fields with fixed direction but variable magnitude.

Note that the scheme of Genoni et al. [18], called the Magnetized Implicit (MI) scheme, falls into both of the first two categories in the sense that it uses an effective force to capture the mirror force ($-\mu\nabla\mathbf{B}\cdot\mathbf{b}$) and an effective velocity to capture the perpendicular drift $\mathbf{v}_{\nabla B}$.

In summarizing previous efforts, we will focus on the first category because it is more compatible with energy conservation than the second category - we make the reasons for this clear in Section 3 - and applicable to more general field geometries than the third category.

2.5.1. Brackbill-Forslund-Vu

The scheme developed by Brackbill, Forslund, and Vu in [2, 26] (denoted BFV henceforth) introduces an effective force $\mathbf{F}_{BFV}$ defined by

[TABLE]

where all quantities without superscripts are evaluated at $t^{n+1/2}$. Using the notation in Section 2.2, straightforward manipulation reveals that

[TABLE]

From here, one can read off the following limits:

[TABLE]

If $\mathbf{v}_{E}$ changes little in a time step, then $\left\|\mathbf{v}_{\perp}^{n}-\mathbf{v}_{E}^{n+1/2}\right\|\approx u_{\perp}^{n}$, since $\mathbf{u}$ and $\mathbf{v}_{E}$ dominate other drifts in the guiding center limit. Thus, $\tilde{\mu}\approx\mu$ for $\Omega_{c}\Delta t\gg 1$, and $\mathbf{F}_{BFV}$ has the features necessary to capture the $\nabla B$ drift for large $\Omega_{c}\Delta t$.

We can actually go a step further than this. Recall that applying the analysis in [22] to Crank-Nicolson, we were able to show that unmodified Crank-Nicolson features an effective $\nabla B$ force of $-\mu_{\textrm{eff}}^{\textrm{CN}}\nabla B$, with $\mu_{\textrm{eff}}^{\textrm{CN}}$ given in (20). The BFV scheme proposes to modify Crank-Nicolson so that the total force due to $\nabla B$ is

[TABLE]

In the first equality, we have just used (22) and retained our assumption that the gyration velocity is approximated $\mathbf{v}_{\perp}-\mathbf{v}_{E}$. Thus, remarkably, the BFV scheme recovers - to leading order - the correct $\nabla B$ force not merely in the limits $\Omega_{c}\Delta t\ll 1$ and $\Omega_{c}\Delta t\gg 1$, but also for every value of $\Omega_{c}\Delta t$!

This scheme has two noteworthy drawbacks. First, and most important for us, is that it does not satisfy the constraint $\mathbf{v}^{n+1/2}\cdot\mathbf{F}_{BFV}=0$, so it does not conserve energy. Second, it requires the direct computation of $\nabla B$ at every time-step - in fact, since an iterative solver is typically used to compute the implicit particle update, it must be computed many times per step.

2.5.2. Magnetized Implicit Scheme

The MI scheme of Genoni et al. [18] introduces the effective force

[TABLE]

Strictly speaking, MI uses only the parallel component, $\left(\mathbf{F}_{MI}^{n+1/2}\cdot\mathbf{b}^{n+1/2}\right)\mathbf{b}^{n+1/2}$ - the perpendicular drift is captured by an effective velocity instead. However, we show that the full $\mathbf{F}_{MI}$ approximates all components of $-\mu\nabla B$ - note that this is not shown in [18].

Note that, using results from Section 2.2, we may rewrite the first term in the cross product defining $\mathbf{F}_{MI}$ as

[TABLE]

By Taylor expansion, the second term in the cross product may be approximated by

[TABLE]

Thus, we can write

[TABLE]

where all quantities without superscripts are evaluated at $t^{n+1/2}$.

Let us now consider the gyroaverage, denoted by $\langle\cdot\rangle$, of this expression. To leading order, we have

[TABLE]

where $\mathbf{B}_{gc}$ denotes the field at the guiding center position, we assumed that $\mathbf{B}$, the drift velocities, and parallel motion change little in a gyroperiod, and the notation $\tilde{\mathbf{v}}=\mathbf{v}-\mathbf{v}_{E}$ has been reintroduced. We may then follow the standard $\nabla B$ drift derivation in [20, 19] - namely, expand the triple product and note that $\left\langle\tilde{\mathbf{v}}_{\perp}\tilde{\mathbf{v}}_{\perp}\right\rangle\approx\left\langle\mathbf{u}\mathbf{u}\right\rangle=u_{\perp}^{2}(\mathbf{I}-\mathbf{b}\mathbf{b})/2$ - to arrive at

[TABLE]

By once again using (17), we find

[TABLE]

This clearly tends to $-\mu\nabla B$ for large $\Omega_{c}\Delta t$, and is $O(\Delta t^{2})$ for $\Omega_{c}\Delta t\ll 1$, as desired. Additionally, note that similarly to the BFV scheme, this effective force has the property that

[TABLE]

to leading order. Thus, at least in a gyroaveraged sense, the MI scheme recovers the correct leading order $\nabla B$ force for every value of $\Omega_{c}\Delta t$.

A key difference, of course, between $\mathbf{F}_{MI}$ and $\mathbf{F}_{BFV}$ is that the former only approximates the correct $\nabla B$ force in a gyroaveraged sense. The MI scheme thus must rely on the time-stepping process to, in some sense, perform the gyro-averaging implicitly by sampling the particle at different gyro-phases at different time-steps. This works well in practice for moderate $\Omega_{c}\Delta t$ - see [18] and our numerical results below - but introduces anomalous drifts for very large $\Omega_{c}\Delta t$.

This limitation will be analyzed in detail in the context of our scheme in the following section, but may be intuitively understood as follows: For very large $\Omega_{c}\Delta t$, the gyration velocity rotates by approximately $\pi$ at every step. It thus takes many time steps to accumulate a representative sample of gyrophases. Before that number of steps is reached, a biased estimate of the gyroaverage is implicitly computed, which gives rise to anomalous displacements. At smaller time steps, an accurate gyroaverage is computed in just a few time steps, and the size of the anomalous displacements is reduced so as to be smaller than the gyroradius.

Of course, this effective force also does not conserve energy. In the full MI scheme, a modification to the perpendicular motion is made to compensate for this, so that it remains true that the magnetic field does no work. Even with that correction, though, the inclusion of an effective velocity to capture $\mathbf{v}_{\nabla B}$ prevents this scheme from conserving even individual particle energy, much less total energy in the context of an implicit PIC scheme.

This scheme has two benefits. First, it does not require an explicit computation of $\nabla B$. Second, as it is presented in [18], it is actually an explicit scheme, using a two-stage predictor-corrector procedure with the quantities evaluated at the $t^{n+1}$ using the values from the predictor stage. This has obvious benefits for computation speed and implementation simplicity, but further reduces accuracy and conservation as seen in our numerical tests.

3. The $\nabla B$ Drift and Energy Conservation

In seeking to preserve total energy conservation in the context of implicit PIC simulation, we immediately rule out the use of an effective velocity, favoring instead an effective force. This is because it is clear that an effective velocity is likely to break energy conservation. We see this from the following logic: Denote a particle’s change in kinetic energy in a single time step by $\Delta K_{p}$, and its change in potential energy by $\Delta P_{p}$. The energy conservation property we seek is written $\sum_{p}(\Delta K_{p}+\Delta P_{p})=0$. The velocity update is ignorant of modifications to the position update to leading order, so the addition of an effective velocity will not alter $\Delta K_{p}$, but will alter $\Delta P_{p}$ (the particle’s position, and thus electrostatic potential energy, is modified).

One has very little hope of maintaining conservation when $\Delta P$ is modified while $\Delta K$ is fixed. Indeed, each particle’s new effective velocity would have to conspire to modify its potential energy in such a way that the sum over all modifications is zero - and this must be true for arbitrary field configurations! This circumstance is made even more unlikely by the fact that the effective velocity is meant to capture the $\nabla B$ drift, which is agnostic with respect to the electric field. Altering $\mathbf{E}$ thus should not alter the effective velocity, meaning it must surely alter $\Delta P$.

In contrast, the addition of an effective force modifies both the position and velocity updates in a self-consistent manner - the position is modified indirectly by the modification of the velocity. Aside from the constraint mentioned above that the new force be orthogonal to $\mathbf{v}^{n+1/2}$, the only other place in the implicit PIC energy conservation proofs of [8, 9] that the particle velocity arises is in the computation of the current density. All that is relied upon there is that the position and velocity are related by the Crank-Nicolson position update, which is unmodified by an effective force. Thus, we see that the introduction of an effective force holds the promise of energy conservation (it conserves energy if and only if it is orthogonal to $\mathbf{v}^{n+1/2}$), while the introduction of an effective velocity is extremely unlikely to do so.

3.1. Conservative Effective Force

Our guiding principle for finding an effective force that both approximates $-\mu\nabla B$ and is orthogonal to $\mathbf{v}^{n+1/2}$ is as follows: the effective force is postulated to be the projection of some vector $\mathbf{G}$ onto the orthogonal complement of $\mathbf{v}^{n+1/2}$. That is,

[TABLE]

for some suitably chosen $\mathbf{G}$. An intuitively reasonable choice for $\mathbf{G}$ would seem to be a scalar multiple of either $\mathbf{F}_{BFV}$ or $\mathbf{F}_{MI}$. This will turn out to be true in certain cases, which we show below.

This simplistic view, unfortunately, has the problem that it artificially mixes parallel and perpendicular motion. Consider particle motion in the fields $\mathbf{E}=0$, $\mathbf{B}=B(x,y)\widehat{\mathbf{z}}$. In reality, a particle’s velocity in the $z$-direction is constant in these fields, since $\nabla B\cdot\mathbf{b}=0$. However, the effective force in (33) will in general have a non-zero component in the $z$-direction. In particular, the $z$-component is non-zero if the particle initially has non-zero $v_{z}$ and $\mathbf{v}^{n+1/2}$ is not orthogonal to $\mathbf{G}$. Thus, introduction of this force will modify the parallel velocity even when parallel velocity is physically constant.

This issue is fixed by breaking our projection into parallel and perpendicular components. Instead of (33), we use the modified projection

[TABLE]

where all instances of $\mathbf{b}$ and $\mathbf{v}$ are evaluated at $t^{n+1/2}$, and the $\perp$ and $\parallel$ symbols refer to components orthogonal and parallel, respectively, to $\mathbf{b}^{n+1/2}$. Note that each of the two terms appearing in (34) is independently orthogonal to $\mathbf{v}^{n+1/2}$. The first term is responsible for the mirror force - i.e. $-\mu\nabla B\cdot\mathbf{b}$ - with a modification to the perpendicular motion to conserve energy, and the second is intended to give rise to $\mathbf{v}_{\nabla B}$. In the example above in which parallel velocity should be constant, this effective force behaves correctly so long as $G_{\parallel}=0$, since this implies $F_{cons,\parallel}=0$.

We now return to our goal, which is that $\mathbf{F}_{cons}$ should, in some sense, approximate $-\mu\nabla B$ when $\Omega_{c}\Delta t\gg 1$. One consequence of this goal is that we require $\mathbf{F}_{cons}$ to be independent of gyrophase. However, we are also requiring $\mathbf{F}_{cons,\perp}$ to be orthogonal to $\mathbf{v}_{\perp}$, which manifestly does depend on gyrophase. It is thus unrealistic to expect $\mathbf{F}_{cons}$ to approximate the $\nabla B$ force at every time-step.

We can, however, enforce that $\mathbf{F}_{cons}$ approximate the correct $\nabla B$ force in a gyro-averaged sense, as in the MI scheme. That is, we will specify the vector $\mathbf{G}$ by insisting that

[TABLE]

and requiring additionally that $\mathbf{G}$ be independent of gyrophase. Note that we have chosen the gyroaverage to equal $\mathbf{F}_{BFV}$ so that we retain the desirable property enjoyed by both the BFV and MI schemes that the $\nabla B$ force is recovered correctly for every value of $\Omega_{c}\Delta t$. We have chosen $\mathbf{F}_{BFV}$ instead of $\mathbf{F}_{MI}$ because the latter is dependent on gyrophase, and the gyroaverages of the two are identical to leading order as shown above.

It remains only to compute $\mathbf{G}$ given the constraints above. We will assume that, as in the particle drift derivation, $\mathbf{v}$ is dominated jointly by gyration $\mathbf{u}$, the $\mathbf{E}\times\mathbf{B}$ drift $\mathbf{v}_{E}$, and parallel velocity $v_{\parallel}$. We write $\mathbf{v}_{\perp}=\mathbf{u}+\mathbf{v}_{E}$ to leading order, and assume that $v_{\parallel}$, $\mathbf{v}_{E}$, and $\mathbf{b}$ all change on time-scales much longer that $\Omega_{c}^{-1}$. We can then write the gyroaverage condition above as

[TABLE]

Only one of the terms on the right is parallel to $\mathbf{b}$, so we immediately find that $G_{\parallel}=F_{BFV,\parallel}$. Note that this does recover constant parallel velocity in the simple case described above.

By evaluating the remaining three distinct gyroaveraged expressions - see Appendix A for more detail - and substituting in the known value of $G_{\parallel}$, one can simplify the perpendicular components of the expression above to find

[TABLE]

where $\eta=\min\{1,u/v_{E}\}$ and we have introduced $\widehat{\mathbf{v}}_{E}=\mathbf{v}_{E}/v_{E}$ for brevity. Note that it is not surprising that the cases $u>v_{E}$ and $u<v_{E}$ are qualitatively different. Indeed, for $u>v_{E}$, $\mathbf{v}_{\perp}/v_{\perp}$ traverses the full unit circle in the course of a gyration, while in the oposite case it has a positive component in the $\mathbf{v}_{E}$ direction for every gyrophase.

It is thus sensible to break further analysis into cases: $u>v_{E}$ (i.e. $\eta=1$) and $u<v_{E}$ (i.e. $\eta=u/v_{E}<1$). In the former case, the entire first term vanishes and we find

[TABLE]

Again, this is unsurprising because $\mathbf{v}_{\perp}/v_{\perp}$ traverses the full unit circle in this case, so the average of the projection operator onto its orthogonal complement is $\mathbf{I}/2$.

The case $u<v_{E}$ is more algebraically cumbersome. After some simplification - in particular, noting that $v_{E}/u=\eta^{-1}$ in this case - we can simplify (37) to

[TABLE]

This vector equation can be solved for $\mathbf{G}_{\perp}$ by decomposing it into components parallel and orthogonal to $\widehat{\mathbf{v}}_{E}$. By dotting with $(\mathbf{I}-\widehat{\mathbf{v}}_{E}\widehat{\mathbf{v}}_{E})$ and $\widehat{\mathbf{v}}_{E}$ respectively then solving for the relevant component of $\mathbf{G}_{\perp}$, we find

[TABLE]

From here, it is straightforward to write an expression for $\mathbf{G}_{\perp}$ in the case $u<v_{E}$:

[TABLE]

We can now see by comparing (41) to (38) - along with the knowledge that $G_{\parallel}=F_{BFV,\parallel}$ - that the following is a uniformly valid expression for $\mathbf{G}$:

[TABLE]

where $\mathbbm{1}_{A}$ denotes the indicator function, i.e., $1$ if $A$ holds and zero otherwise.

For completeness, we summarize here the full form of our proposed energy-conserving, asymptotic preserving particle update:

[TABLE]

where $\mathbf{G}$ is given by (42) with all velocities evaluated at $t^{n+1/2}$, and $\mathbf{F}_{BFV}$ is defined by (21).

There are three interesting points that warrant elaboration here. Firstly, the relative size of $u$ and $v_{E}$ is dependent on $\Delta t$. This is because the relevant values of $u$ and $v_{E}$ are those at $t^{n+1/2}$; while $\mathbf{v}_{E}$ changes little in a time-step, the rapid variation of $\mathbf{u}$ has already been shown to lead it to shrink for increading $\Omega_{c}\Delta t$ - see (17), where we have assumed $\mathbf{E}=0$ so $\mathbf{v}_{\perp}\approx\mathbf{u}$. The consequences of this fact will be expanded upon in the following sections.

Secondly, note that $\mathbf{G}_{\perp}$ diverges as $u/v_{E}\rightarrow 0$ (i.e. $\Omega_{c}\Delta t\rightarrow\infty$ if $v_{E}\neq 0$). This may initially appear alarming, but is accounted for by the fact that $\mathbf{G}_{\perp}$ itself is never actually used in the velocity update - only its projection onto the orthogonal complement of $\mathbf{v}_{\perp}$. In the same limit $u/v_{E}\rightarrow 0$, $\mathbf{v}_{E}$ dominates $\mathbf{v}_{\perp}$, so that this projection is very nearly onto the orthogonal complement of $\mathbf{v}_{E}$. Since the divergence of $\mathbf{G}_{\perp}$ is in the $\mathbf{v}_{E}$ direction, this projection results in the actual effective force remaining well-behaved as $\eta\rightarrow 0$. The reliance of the scheme on this precise cancellation, though, does indicate that it is crucial that the implicit solve be performed accurately. In particular, $\eta$ should be computed self-consistently (as opposed to using some explicit estimation of $\eta$) since the scheme is extremely sensitive to its precise value when $\eta\ll 1$.

Thirdly, one can see that $\mathbf{G}$ is in general discontinuous on the manifold defined by $u=v_{E}$. Indeed,

[TABLE]

due to the presence of the indicator function. This, again, is a consequence of the qualitative difference in the behavior of $\mathbf{v}_{\perp}/v_{\perp}$ over a gyro-orbit depending on which of $u$, $v_{E}$ is larger. Of course, by construction, the gyroaverage of the effective force arising from $\mathbf{G}$ is continuous, so this discontinuity is only expected to be an issue for the small subset of particles for which $u$ and $v_{E}$ are very close to equal, so that there is a risk of their relative size changing from time-step to time-step. Our adaptive time-stepping strategy outlined below will endeavor to choose time-steps so that this is not the case.

3.2. Time-step Restrictions

While an asymptotic preserving scheme admits the use of time-steps with $\Omega_{c}\Delta t\gg 1$, one of course cannot expect accurate results for arbitrarily large time-steps. Understanding the restrictions on $\Delta t$ is crucial when applying the scheme to problems of scientific interest, as failure to do so risks viewing erroneous results as legitimate.

An obvious restriction is that the time-step should be sufficiently small that variations in $\mathbf{E}$ and $\mathbf{B}$ are well-resolved. The scheme proposed here features two additional time-step restrictions that we now describe. After deriving the time-step restrictions on the scheme, we present an adaptive time-stepping scheme designed to ensure these constraints are respected at every step.

Like the MI scheme but unlike BFV, our proposed scheme (43) has an effective force that only approximates $-\mu\nabla B$ in a gyroaveraged sense. That is, we must rely on the time-stepping process to perform a sort of “implicit” gyroaverage for us - here we use the term “implicit” not in the sense that a system of equations must be solved, but in the sense that no gyroaverage operator explicitly appears in our time integration scheme. This fact was alluded to when discussing MI; we discuss it in considerably more detail here.

Recall that the $\nabla B$ drift in the perpendicular direction that we seek to capture has the form

[TABLE]

The scheme proposed here, on the other hand, gives rise to a drift of the form

[TABLE]

The scheme thus features an anomalous drift with velocity

[TABLE]

By construction, the gyroaverage of this anomalous displacement vanishes. However, on shorter time-scales, it has the potential to impact the particle’s long-term trajectory. This can be avoided by enforcing two conditions: (a) the anomalous displacement due to this anomalous drift should be no larger than the gyroradius, and (b) the time-scale over which this anomalous displacement averages to zero should be small compared to time-scales of interest, denoted henceforth by $\tau_{\textrm{res}}$.

Deriving the restrictions on time-step arising from these conditions is quite challenging in general, but can be done analytically in two important limits: (i) $u^{n+1/2}\gg v_{E}$, and (ii) $v_{E}\gg u^{n+1/2}$. The derivations of these restrictions are lengthy, and therefore confined to Appendix B. We summarize the results in Table 1, in which we record the maximum value of $\Omega_{c}\Delta t$ permitted by each restriction in each of the two asymptotic limits above.

The relevant non-dimensional quantities in determining these restrictions on $\Omega_{c}\Delta t$ are defined by

[TABLE]

The quantity $\tau_{\textrm{res}}$ denotes the shortest time-scale in our problem that we wish to resolve. For brevity, we quote here only simplified restrictions under the assumptions $\delta_{\parallel},(\Omega_{c}\tau_{\textrm{res}})^{-1}\ll 1$. The more general expressions may be found in Appendix B.

3.3. Adaptive Time-stepping

Having understood the limitations on our scheme’s time-step, we propose an adaptive time-stepping strategy that respects these limits. In addition, the time-stepping strategy should endeavor to avoid the $u^{n+1/2}=v_{E}$ manifold on which our effective force is discontinuous. We specify parameters $\alpha$, $\beta$, $\epsilon\in(0,1)$ and $\Gamma>0$. We then select a time-step as follows:

1. (1)

Compute all relevant quantities - $\mathbf{E}$ and $\mathbf{B}$, their gradients, $\mathbf{v}_{E}$, etc. - at the current particle position. 2. (2)

Set $(\Omega_{c}\Delta t)_{E}=\alpha\min\left\{\sqrt{2}\left(\delta_{E}+\delta_{\parallel}\right)^{-1/2},\Omega_{c}\tau_{res},\Gamma\delta_{\perp}^{-1}\right\}$. 3. (3)

Estimate $\mathbf{u}$ at the half time-step with $\tilde{u}^{n+1/2}=\left\|\mathbf{v}_{\perp}-\mathbf{v}_{E}\right\|/\sqrt{1+(\Omega_{c}\Delta t)_{E}^{2}/4}$. 4. (4)

If $v_{E}>(1+\beta)\tilde{u}^{n+1/2}$, set $\Delta t=(\Omega_{c}\Delta t)_{E}/\Omega_{c}$ and proceed with the time-step. 5. (5)

Else, if $v_{E}>(1-\beta)u^{n}$, set $\Delta t=\min\{\tau_{\textrm{res}},\Gamma\Omega_{c}^{-1}\}$ and proceed with the time-step. 6. (6)

Else, set $(\Omega_{c}\Delta t)_{u}=\alpha\min\left\{2\min\{\sqrt{2}/\sqrt{\delta_{\perp}},1/\delta_{\parallel}\},2\sqrt{\Omega_{c}\tau_{\textrm{res}}/\pi}\right\}$. Then set $\Delta t=(\Omega_{c}\Delta t)_{u}/\Omega_{c}$ and proceed with the time-step. 7. (7)

If, upon completion of the time-step, it is found that the fractional change in magnetic moment $\Delta\mu/\mu$ during the time-step exceeds $\epsilon$, shrink the step size by a factor of $\alpha\epsilon\mu/\Delta\mu$ and recompute the step.

Several details of this time-step selection warrant elaboration. First, the parameters $\alpha$, $\beta$ and $\Gamma$ have simple interpretations. $\alpha$ controls how close to the time-step restrictions of Table 1 one is willing to get. $\alpha=1$ corresponds to steps at the limit of what is valid, while smaller values of $\alpha$ correspond to more conservative strategies. In our numerical experiments, we use $\alpha=0.9$.

The parameter $\beta$ controls how close to the $u^{n+1/2}=v_{E}$ manifold - on which our effective force is discontinuous - we allow ourselves to approach. If there is no way to respect the time-step constraints while keeping either $v_{E}/u^{n+1/2}>1+\beta$ or $v_{E}/u^{n+1/2}<1-\beta$, then the gyroperiod is resolved (see step 5) so that the discontinuity is no longer a concern. In our numerical experiments, we use $\beta=0.2$. Note that the circumstances in which this restriction is active are extremely rare. Indeed, being required to resolve $\Omega_{c}^{-1}$ requires both that $v_{E}/u^{n}\in(1-\beta,1+\beta)$ and the maximum allowable time-step to be not much larger than $\Omega_{c}^{-1}$. In quasi-neutral plasmas, one typically expects $v_{E}\ll u^{n}$, so this constraint only applies to a select few particles for which a large time-step was impossible in the first place.

The parameter $\Gamma$ measures the accuracy with which we wish to resolve spatial variations in the magnetic field. Indeed, note that one component of the constraint in step 2 reads

[TABLE]

The right-side measures fractional variation in the magnetic field due to perpendicular velocity. $\Gamma=1$ corresponds to resolving the spatial scale over which $B$ doubles. In our numerical experiments, we choose $\Gamma=0.1$, corresponding to resolving $10\%$ variations in $B$. Note that $\Gamma$ also appears in the time-step chosen when $\Omega_{c}^{-1}$ must be resolved (step 5). It is not necessary that these two constants be equal, but $0.1$ turns out to be a reasonable value for each.

The parameter $\epsilon$ controls the maximum permissible fractional change in $\mu$ within a time-step. It appears in step 7 as a safeguard for situations in which the particle is weakly magnetized. When $\mu$ varies non-adiabatically, it is typically important to resolve these changes in $\mu$. This is ensured by choosing a relatively small value of $\epsilon$. Our numerical tests use $\epsilon=0.15$.

A final note regarding $\tau_{\textrm{res}}$: In general, this is a problem-dependent quantity to be determined by the user. A reasonable starting point, however, that is used in the majority of our numerical tests is as follows. Begin by defining lengths scales

[TABLE]

arising from perpendicular and parallel variations in $B$ and $v_{E}$. A time-scale to be resolved can then be defined by

[TABLE]

4. Numerical Examples

We test our scheme by investigating single-particle motion in four electromagnetic field configurations. We use a simple magnetic mirror configuration, a case with transverse magnetic field gradient in a uniform electric field, a case in which the particle passes through a weakly magnetized region, and a Solov’ev equilibrium in tokamak geometry with a flux-function electric field.

In all cases, we work in a dimensionless formulation in which $q/m=1$, and measure both accuracy and conservation properties of the scheme. To measure accuracy, we compare to a simulation using the standard Boris integrator - used by the majority of PIC schemes - with $\Omega_{c}\Delta t\leq 0.1$. We solve the systems of equations arising from our implicit discretizations using a Jacobian-Free Newton Krylov (JFNK) method with tolerance set to $10^{-12}$.

The adaptive time-stepping scheme and definition of $\tau_{\textrm{res}}$ found above are used except where otherwise noted. When individual values of $\Omega_{c}\Delta t$ are reported, these are averages over the entirety of the simulation. When noteworthy, we plot the value of $\Omega_{c}\Delta t$ as a function of time to illuminate the behavior of the adaptive time-stepping scheme.

In each test case, we compare our scheme against not only a well-resolved Boris run, but also the BFV scheme, MI, and unmodified Crank-Nicolson. Each scheme uses the adaptive time-stepping scheme described above for a fair comparison.

4.1. Magnetic Mirror

As is standard, our magnetic mirror is formed by two circular current loops. Let the loops be oriented normal to the $z$-axis, centered at $(0,0,\pm L/2)$, with radius $r$ and carrying current $I$. The resulting magnetic field (up to constants) near the $z$-axis is given by

[TABLE]

We work with $r=3$, $L=15$, $I=3\times 10^{4}$, leading to a mirror ratio of $R_{mirr}=9.8342$. We initialize the particle at $\mathbf{x}^{0}=0.02\widehat{\mathbf{y}}$ with $\mathbf{v}^{0}=\widehat{\mathbf{x}}+v^{0}_{\parallel}\widehat{\mathbf{z}}$. To a good approximation, this results in gyration about the $z$-axis. Having fixed the perpendicular velocity at $1$, the trapped-passing boundary may be expressed purely in terms of a critical initial parallel velocity $v_{crit}$ given by

[TABLE]

We perform a parameter scan in the initial parallel velocity $v_{\parallel}^{0}$ to test how accurately each scheme captures the trapped-passing boundary. The particle trajectories in the $x$-$z$ plane for several particular values of $v_{\parallel}^{0}$ are found in Figure 1.

As expected, Crank-Nicolson dramatically underestimates the mirror force for $\Omega_{c}\Delta t\gg 1$, and confinement is not achieved for any of the parallel velocities tested. BFV correctly predicts confinement for parallel velocities relatively small compared to $v_{crit}$, but underestimates the trapped-passing boundary by more that $40\%$.

In this case, both the present scheme and MI perform quite well, predicting the trapped-passing boundary quite accurately. An understanding of the reason for this can be gleaned from Figure 2, in which we measure energy conservation during a sample simulation.

As expected, our scheme and Crank-Nicolson conserve energy up to the tolerance of the nonlinear solver. In the absence of an electric field, MI and Boris conserve energy up to numerical round-off. BFV, on the other hand, features $O(1)$ errors in energy due to the form of the effective force introduced, leading to its poor performance.

In addition to the $x$-$z$ plane, we plot $z$ as a function of time for each method in Figure 3, this time focusing on values of $v_{\parallel}^{0}$ near $v_{crit}$.

When focusing on values of $v_{\parallel}^{0}$ near the trapped-passing boundary, the improved performance of the new scheme, even relative to MI, is readily visible.

4.2. Transverse Magnetic Field Gradient

The previous test exercises only the mirror force, in the absence of an electric field. Here, we test each scheme’s ability to capture the transverse $\nabla B$ drift in a homogeneous electric field. We set $\mathbf{B}=(10+x)\widehat{\mathbf{z}}$, $\mathbf{E}=\widehat{\mathbf{y}}$, and initialize the particle at the origin with $\mathbf{v}^{0}=\widehat{\mathbf{x}}$. Simulations are run to the final time $T=300$.

These fields induce an $\mathbf{E}\times\mathbf{B}$ drift in the positive $x$-direction and a $\nabla B$ drift in the positive $y$-direction. This test case affords the opportunity to measure magnetic moment conservation and ability to capture the polarization drift, the latter arising because the $\mathbf{E}\times\mathbf{B}$ speed is time-varying. The trajectory in the $x$-$y$ plane as well as $\mu$ as a function of time appear in Figure 4.

As is readily apparent, the new scheme features improved accuracy in capturing the trajectory. This can be directly linked to improved magnetic moment conservation: other schemes underestimate $\mu$, so they also underestimate the $\nabla B$ drift, for it is proportional to $\mu$. The reason for the failure of earlier schemes - especially BFV - to conserve $\mu$ is simple to understand. Fundamentally, adiabatic invariance of $\mu$ is a result of the energy equation

[TABLE]

followed by substituting in the drift motion for $\mathbf{v}$, gyroaveraging, and neglecting small terms. When the “effective” force $-\tilde{\mu}\nabla B$ is introduced as in BFV, a new term appears in the energy equation above:

[TABLE]

In addition to leading to $O(1)$ reduction in total energy - see Figure 5 - this term also means that when retracing the steps in the $\mu$-conservation derivation, one must find $\dot{\mu}\approx-\tilde{\mu}\mathbf{v}_{E}\cdot\nabla B$, where we have again assumed that $\mathbf{v}_{E}$ dominates the drift motion. In our case, $\mathbf{v}_{E}$ and $\nabla B$ are parallel, so for $\Omega_{c}\Delta t\gg 1$ (i.e. $\mu\approx\tilde{\mu}$), one expects exponential decay of $\mu$ at rate $v_{E}\left\|\nabla B\right\|$, just as observed.

Crank-Nicolson, on the other hand, does conserve energy. However, its underestimation of the $\nabla B$ drift keeps the particle in a region of artificially high electrostatic potential. This must be compensated for by robbing the particle of kinetic energy. Since we’ve already shown that the $\mathbf{E}\times\mathbf{B}$ drift is accurately captured, that kinetic energy must come from the particle’s gyration, thus leading to decay of $\mu$.

Note also that while MI conserves energy to numerical roundoff in the absence of an electric field - see magnetic mirror test case above - it features $O(1)$ errors in the presence of even a uniform electric field due to its usage of an effective velocity. Finally, we note that our adaptive time-stepping strategy chooses larger time-steps for BFV, MI, and Crank-Nicolson than for our scheme near the end of the simulation. This is a result of the shrinking $\mu$ value those schemes experience, which shrinks the gyroradius and leads the scheme to believe a larger time-step is possible.

4.3. Weakly Magnetized Region

A strong motivation for an AP scheme capable of stepping over the gyrofrequency is that, unlike drift- or gyro-kinetics, it can handle weakly magnetized regions in which the guiding-center approximation is invalid. We test this capability in a simple geometry by setting

[TABLE]

Having already convinced ourselves that C-N fails to capture the $\nabla B$ drift for $\Omega_{c}\Delta t\gg 1$, we omit it in plots here to reduce clutter.

We initialize a particle at the origin, where it experiences an $\mathbf{E}\times\mathbf{B}$ drift in the positive $x$-direction. This pushes it through the narrow region surrounding $x=12$ in which the particle is quite weakly magnetized and the guiding center approximation is violated. Therein, it experiences a large, non-adabatic jump in magnetic moment $\mu$ before continuing to drift to the right into another strongly magnetized region. A sample trajectory and time history of $\mu$ for each scheme appear in Figure 6.

Note that, while our proposed scheme is the only one able to reproduce the adiabatic invariance of $\mu$ in the magnetized regions far from $x=12$, it still does not recover the particle trajectory accurately. This is because the rapid, non-adiabatic jump in $\mu$ that occurs as the particle passes through the weakly magnetized region is very sensitive to the particle’s gyrophase as it enters the weakly magnetized region. Since stepping over the gyroperiod inherently throws away gyrophase information, no scheme using $\Omega_{c}\Delta t\gg 1$ can hope to recover the details of any given particle trajectory in this case.

What one can hope for, however, is that the average jump in magnetic moment is reproduced accurately when averaged over initial particle gyrophase. Indeed, a fundamental assumption in the drift- and gyro-kinetic approximations is that the distribution function is independent of gyrophase. In this case, recovering the average $\mu$ jump is sufficient to predict the ensemble dynamics. To this end, we initialize particles at the origin with velocities $\mathbf{v}^{0}=\left(\cos\gamma,\sin\gamma\right)$, $\gamma=2\pi n/64$, $n=0,...,63$ and measure the final magnetic moment of the particle averaged over the interval $t\in[100,120]$. The final magnetic moment as a function of initial gyrophase for each scheme is plotted in Figure 7.

As expected, the magnetic moment is not accurately reproduced for most particular gyrophases. However, the present scheme reproduces the mean $\mu$ of the fully-resolved Boris run to within $0.8\%$, and the standard deviation to within $14\%$. For detailed values, see Table 2.

Recall that a drift- or gyro-kinetic approximation would predict no change in $\mu$ in this (or any other) case.

As a final note, we plot $\Omega_{c}\Delta t$ as a function of time and the energy errors for this test case in Figure 8. We use initial gyrophase $\gamma=2\pi\times 63/64$ for these plots, as in Figure 6.

As expected, our adaptive time-step strategy permits large time-steps except when the particle is in the weakly magnetized region, in which all scales are resolved. Observe that, as in the previous test case, the underestimation of $\mu$ by BFV and MI leads them to take larger time-steps toward the end of the simulation.

4.4. Tokamak Equilibrium

As a final test, we consider an equilibrium magnetic field in a tokamak-like geometry. Simple analytic solutions of the Grad-Shafranov solutions described in [4] specify the poloidal field, and the toroidal field is assumed to arise from a line current along the $z$-axis. In particular, a Solov’ev equilibrium solution of the Grad-Shafranov equation for the flux function $\psi$ can be written as

[TABLE]

where $C$, $d_{i}$ are constants and $(r,\phi,z)$ is the standard polar coordinate system. The poloidal magnetic field is then given by $\mathbf{B}_{p}=\nabla\psi\times\mathbf{e}_{\phi}/r$. In [23], it is shown that the coefficients $d_{i}$ may be found uniquely by specifying the tokamak parameters $\varepsilon$ (inverse aspect ratio), $\kappa$ (elongation), and $\delta$ (triangularity). Our poloidal field is thus uniquely specified by setting $C=300$ ($C$ simply controls the overall strength of the poloidal field) and using the shape parameters of the International Thermonuclear Energy Reactor (ITER): $\varepsilon=0.32$, $\kappa=1.7$, and $\delta=0.33$. The toroidal field is chosen to be $\mathbf{B}_{tor}=800\mathbf{e}_{\phi}/r$.

Additionally, we introduce an electrostatic potential, which is itself a flux function. This choice is motivated by the fact that, in steady state, the ideal MHD Ohm’s law ($\mathbf{E}+\mathbf{v}\times\mathbf{B}=0$) requires that $\mathbf{B}\cdot\nabla\phi=0$. We define the electrostatic potential by $\phi=-\psi/5$. This non-trivial spatial dependence of the electric field necessitates an additional subtlety in order to retain energy conservation in the single-particle case. Indeed, Crank-Nicolson only features exact energy conservation when the electric field varies linearly in space, so that $\left(\mathbf{x}^{n+1}-\mathbf{x}^{n}\right)\cdot\mathbf{E}^{n+1/2}=\phi^{n}-\phi^{n+1}$ identically. In [24], a trick is introduced to recover exact energy conservation in the general case; in the present context, the trick consists of the replacement

[TABLE]

in the velocity update equation. Note that this replacement results in a scheme that is still second-order accurate, as the new pre-factor tends to one at second order. Note further that this trick is not necessary for the total energy conservation theorems of implicit PIC - again, see [8, 9]. We simply use it here demonstrate that exact conservation can be achieved in the single particle case as well.

We initialize a particle in these fields at $\mathbf{x}=1.2\widehat{\mathbf{x}}$ with velocity $\mathbf{v}=\widehat{\mathbf{x}}+v_{y}^{0}\widehat{\mathbf{y}}+\widehat{\mathbf{z}}$. The $y$-velocity is left variable so as to again study the trapped-passing boundary - at the initial particle location, the magnetic field is predominantly toroidal, making $v_{y}$ a reasonable proxy for $v_{\parallel}$. A sample trajectory projected onto the $r$-$z$ plane with $v_{y}^{0}=0.7$ can be found in Figure 9.

Note that BFV incorrectly predicts a passing particle, while MI predicts an unphysical trajectory even after reducing the time-step by a factor of four relative to CBFV and BFV.

We use a binary search to approximate the critical value of $v_{y}^{0}$ that separates the trapped and passing regimes. We report the resulting values and percentage errors - the well-resolved Boris run is taken as the baseline in the absence of analytic theory - in Table 3.

As in previous cases, the improved accuracy of the new method is attributable to improved energy and magnetic moment conservation, each of which is plotted in Figure 10.

5. Conclusions

In this work, we have introduced a first-of-kind numerical method for particle orbit integration that both (a) reproduces guiding center motion, including all first-order drifts, even when stepping over the gyration time-scale, and (b) conserves energy. The time integrator is implicit, built upon the classical Crank-Nicolson algorithm, and is thus appropriate for incorporation into implicit PIC methods. We have presented numerical results for single particle motion in several field configurations that demonstrate the dramatically improved conservation properties and accuracy of the scheme relative to earlier efforts.

Opportunities for further development are myriad. Practical implementation of the method in PIC codes stands to benefit tremendously from effective preconditioning of the nonlinear solve of the orbit integral. Work in this direction, building on results for BFV in [26] that analytically eliminate portions of the nonlinearity, is underway. Additionally, the present scheme can only be said to capture the drift-kinetic limit - truly capturing the gyrokinetic limit requires removing the assumption that the gyroradius is small compared to length scales of interest. Progress in this direction is of obvious interest. Further refinement and generalization of the adaptive time-stepping strategy is also of importance for the scheme’s future practical utility.

Acknowledgements

This work was performed under the auspices of the U.S. Department of Energy by LLNL and LANL under contracts DE-AC52-07NA27344, DE-AC52-06NA25396, and supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. The authors wish to acknowledge valuable conversations regarding this work with Guangye Chen, William Taitano, C-S Chang, Joshua Burby, and Matthew Miecnikowski.

Appendix A Gyroaverage computations

When computing the gyroaverages on the right of (36), we may without loss of generality align $\mathbf{v}_{E}$ with the $x$-axis. Denoting gyrophase by $\gamma$ and writing $\mathbf{u}=u(\cos\gamma,\sin\gamma)^{T}$, we have

[TABLE]

A computer algebra system can evaluate all of these integrals analytically. We find

[TABLE]

where $\eta=\min\{1,u/v_{E}\}$ as in the main text. Recalling that we assumed $\mathbf{v}_{E}$ was in the $x$-direction, we can make these expressions general again with

[TABLE]

Here, again as in the main text, $\widehat{\mathbf{v}}_{E}=\mathbf{v}_{E}/v_{E}$.

Appendix B Derivation of Time-step restrictions

As noted in the main text, the fact that our scheme features a $\nabla B$ drift that differs from the true $\nabla B$ drift velocity induces an anomalous drift that (a) should be kept small and (b) should average to zeros on time-scales of interest. Here, we derive the resulting restrictions on time-step size in the limits (i) $u^{n+1/2}\gg v_{E}$ and (ii) $u^{n+1/2}\ll v_{E}$.

B.1. The case $u^{n+1/2}\gg v_{E}$

As noted in the main text - see (38) - in the case $u^{n+1/2}>v_{E}$, the expression for $\mathbf{G}$ simplified considerably:

[TABLE]

Substituting this into (34) , we find an expression for $\mathbf{v}_{\nabla B}^{\textrm{AP}}$:

[TABLE]

Recalling the relationship between $\mu$ and $\tilde{\mu}$ implied by (22) and subtracting the above from $\mathbf{v}_{\nabla B}$, we find the following expression for the anomalous velocity:

[TABLE]

At this point, it is useful to break $\Delta\mathbf{v}_{\nabla B}$ into two pieces. We write $\Delta\mathbf{v}_{\nabla B}=\Delta\mathbf{v}_{\nabla B}^{1}+\Delta\mathbf{v}_{\nabla B}^{2}$, with

[TABLE]

Note that in the present limit in which gyration dominates the perpendicular velocity, each piece of $\Delta\mathbf{v}_{\nabla B}$ has constant magnitude over time-scales short compared to variations in $\mathbf{B}$ and $v_{\parallel}$. The first piece simply rotates by angle $\theta$ given in (11). In the second piece, the operator $(\mathbf{I}-2\widehat{\mathbf{v}}_{\perp}\widehat{\mathbf{v}}_{\perp})\cdot$ corresponds to a reflection about the $\widehat{\mathbf{v}}_{\perp}$ axis. Elementary geometry can thus be used to show that the second piece rotates by angle $2\theta$ at each time-step.

As already noted, motion whose velocity has constant magnitude with velocity rotating by fixed angle transcribes a circle of radius given in (15). Thus, the displacement induced by $\Delta v_{\nabla B}$ in this limit is a superposition of two circles of radius $R_{1}$ and $R_{2}$, respectively. Elementary manipulation leads to formulae for these radii:

[TABLE]

where $\rho=u/\Omega_{c}$ is the gyroradius and $\delta_{\parallel}$ and $\delta_{\perp}$ are as defined in (48).

To guarantee that the total anomalous displacement never exceeds the gyroradius, it thus suffices to impose $\max\{R_{1},R_{2}\}\leq\rho/2$. This leads to the simultaneous time-step constraints

[TABLE]

Taking the minimum of these two restrictions in the limit $\delta_{\parallel}\ll 1$ gives the time-step restriction (ia) in Table 1.

We get the other restriction by analyzing the time required to traverse each circle and asking that it be less than the smallest time-scale of interest, denoted $\tau_{\textrm{res}}$. The traversal time is simply the circumference of the circle - approximately $2\pi R$ - over the anomalous drift speed. For each of the two circles in question, we thus find traversal times

[TABLE]

Insisting that both these quantities be smaller than $\tau_{\textrm{res}}$ leads to

[TABLE]

In the limit $\Omega_{c}\tau_{\textrm{res}}\gg 1$, this reduces to the restriction (ib) in Table 1. We note that no time-step satisfies this constraint if $\Omega_{c}\tau_{\textrm{res}}<\pi$. However, in this case there is no separation in time-scales between gyration and $\tau_{\textrm{res}}$ and we are forced to resolve the gyration scale at any rate.

B.2. The case $u^{n+1/2}\ll v_{E}$

We note that this limit corresponds to $\eta\ll 1$. We thus keep only the terms proportional to $\eta^{-2}$ in evaluating the anomalous drift velocity in this case. Using the general form of $\mathbf{G}$ in (42), we find that to leading order

[TABLE]

Next, we note that

[TABLE]

where in the second approximation we’ve again assumed $u\ll v_{E}$, and $\cos\psi=\widehat{\mathbf{u}}\cdot\widehat{\mathbf{v}}_{E}$ defines $\psi$ as the angle between $\mathbf{u}$ and $\mathbf{v}_{E}$. This leads us to the conclusion

[TABLE]

to leading order in $\eta$, where we’ve used $\widehat{\mathbf{v}}_{\perp}-\widehat{\mathbf{v}}_{E}\approx\eta\widehat{\mathbf{u}}$ in the second approximation.

Substituting this into (71), we have

[TABLE]

Note that $\widehat{\mathbf{u}}$ rotates by $\theta$ at every time-step, and the angle $\psi$ changes by $\theta$ as well so long as $\mathbf{v}_{E}$ changes little in a time-step. Thus, for $\Omega_{c}\Delta t\gg 1$ - i.e. $\theta\approx\pi$ - $\Delta\mathbf{v}_{\nabla B}$ approximately flips sign at every time-step. More generally, for $\theta$ near but not equal to $\pi$, we my reasonably estimate that the anomalous displacement averages to zero over $2\pi/\theta$ time-steps. As such, we may bound the maximum anomalous displacement by its value after $\pi/\theta$ time-steps, and insist that this value be less than $\rho$. That is, we insist that $\|\Delta\mathbf{v}_{\nabla B}\|\Delta t\pi/\theta\leq\rho$. For large $\Omega_{c}\Delta t$, Taylor expansion shows that $\theta\approx\pi-4/\Omega_{c}\Delta t$. Taking the worst case scenario in which $\|\widehat{\mathbf{u}}-\widehat{\mathbf{v}}_{E}\cos\psi\|=2$ and $\nabla B$ is parallel to $\mathbf{v}_{E}$, and recalling the subtlety that the relevant $u$ in the definition of $\eta$ is $u^{n+1/2}$ while in $\mu$ it is $u^{n}$, this constraint reduces to

[TABLE]

Again assuming $\Omega_{c}\Delta t\gg 1$ gives the constraint (iia) in Table 1.

The constraint on the time-scale over which the anomalous displacement averages to zero is in this case simply $\pi\Delta t/\theta\leq\tau_{\textrm{res}}$. Again in the limit $\Omega_{c}\Delta t\gg 1$, this reduces to constraint (iib) in Table 1.