# Asymptotic-preserving gyrokinetic implicit particle-orbit integrator for
arbitrary electromagnetic fields111This 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 (LFR) (17-SC-20-SC), a collaborative effort of the U.S. Department of
Energy Office of Science and the National Nuclear Security Administration, and
the U.S. DOE Office for Advanced Scientific Computing Research (LC).

L.F. Ricketson<EMAIL_ADDRESS>L. Chacón<EMAIL_ADDRESS>Lawrence
Livermore National Laboratory, Livermore, CA 94550, United States  Los Alamos
National Laboratory, Los Alamos, NM 87545, United States

###### Abstract

We extend the asymptotic preserving and energy conserving time integrator for
charged-particle motion developed in [Ricketson & Chacón, JCP, 2020] to
include finite Larmor-radius (FLR) effects in the presence of electric-field
length-scales comparable to the particle gyro-radius (the gyro-kinetic limit).
We introduce two modifications to the earlier scheme. The first is the
explicit gyro-averaging of the electric field at the half time-step, along
with an analogous modification to the current deposition, which we show
preserves total energy conservation in implicit PIC schemes. The number of
gyrophase samples is chosen adaptively, ensuring proper averaging for large
timesteps, and the recovery of full-orbit dynamics in the small time-step
limit. The second modification is an alternating large and small time-step
strategy that ensures the particle trajectory samples gyrophases evenly. We
show that this strategy relaxes the time-step restrictions on the scheme,
allowing even larger speed-ups than previously achievable. We demonstrate the
new method with several single-particle motion tests in a variety of
electromagnetic field configurations featuring gyro-scale variation in the
electric field. The results demonstrate the advertised ability to capture FLR
effects accurately even when significantly stepping over the gyration time-
scale.

##  1 Introduction

The gyrokinetic approximation of the Vlasov equation (GK) [1, 4, 18, 24, 29]
has enabled tremendous advances in both theoretical and computational
understanding of strongly magnetized plasmas. A key advantage of the
gyrokinetic equations is that they allow stepping over the gyration time-
scale, $\Omega_{c}^{-1}=(qB/m)^{-1}$, even when length scales of the order of
the gyroradius are present in the electromagnetic fields (leading to so-called
finite-Larmor-radius effects, FLR). This, in turn, affords enormous
computational speed-ups of the GK model vs. full-orbit Vlasov models. One key
application is magnetic confinement fusion, where gyrokinetic simulations
remain the state of the art.

However, the gyrokinetic approach is not without drawbacks. Realistic boundary
conditions for the gyrokinetic equations are difficult to formulate and
implement due to the non-zero spatial extent of the simulated particles, which
are effectively rotating rings of charge. There is also increasing concern
that the gyrokinetic approximation is insufficiently accurate in portions of
the domain for problems of interest. Examples include field-reversed
configurations (FRCs) and the scrape-off layer (SOL) in certain tokamaks [15,
17, 23, 26, 30, 36, 37, 38, 43].

As a result of these limitations, there is a growing interest in enabling the
simulation of full-orbit kinetic systems when required, at least for ions
[36], without losing the computational advantage of GK models when
appropriate. This can be accomplished by developing so-called asymptotic-
preserving (AP) algorithms that can track the full particle orbit when the
time-step $\Delta t$ is small compared to the gyroperiod, $\Omega_{c}\Delta
t\lesssim 1$, and revert to the gyrokinetic treatment when $\Omega_{c}\Delta
t\gg 1$. The desire to step over a physical time-scale clearly motivates
looking at implicit time integration. We are also motivated to produce a
scheme that can eventually be coupled with strictly charge- and energy-
conserving implicit particle-in-cell (PIC) schemes [6, 9, 10, 11, 12, 13, 39]
without loss of conservation properties.

Recently, the authors introduced [33] an implicit, strictly energy conserving,
AP particle-push algorithm able to recover guiding-center drift-kinetic motion
(i.e., correct gyroradius $\rho$ and all first-order drifts) even when
$\Omega_{c}\Delta t\gg 1$. The scheme overcomes many of the obstacles
mentioned above; namely, the simulated particles can deal with realistic
boundary conditions as in the full-orbit kinetic model, and it enables a
unified treatment of magnetized and unmagnetized regions in a single scheme.
The scheme has been recently reformulated for optimal nonlinear convergence
performance [28] when inverting the associated nonlinear algebraic system
stemming from the implicit discretization. However, the scheme tacitly assumes
that the gyroradius is much smaller than any other scales of interest (because
fields are not gyroaveraged), and therefore it is not expected to capture the
finite Larmor radius (FLR) effects of the gyro-kinetic limit. In fact, our
numerical results here will confirm that it does not.

There have been several recent studies [20, 21, 7, 8, 25, 22] attempting to
develop AP particle-orbit integration schemes similar in spirit to that in
Ref. [33], with varying degrees of generality and complexity. However, none of
these schemes conserve energy exactly (and therefore are not suitable for
coupling with modern strictly energy-conserving implicit PIC codes), and to
our knowledge none of them are able to deal with FLR effects to reproduce the
gyrokinetic limit.

The purpose of the present work is to generalize the AP energy-conserving
scheme in Ref. [33] to capture the necessary FLR effects in the gyrokinetic
limit. In particular, we will show that our generalization can handle an
electric field $\mathbf{E}$ with spatial scales comparable to the gyroradius
$\rho$, i.e., $k_{\perp}\rho\sim 1$, with $k_{\perp}$ a measure of the
characteristic inverse electric-field length scale. Our basic strategy for
capturing FLR effects is relatively standard: rather than evaluating the
electric field at the instantaneous particle location, we compute a gyro-
averaged field using some number $n_{g}$ of equispaced samples along a gyro-
orbit estimated from the instantaneous particle state. The subtlety arises
from two factors: our desire to maintain exact total energy conservation, and
the need to preserve second-order convergence to the full-orbit motion in the
limit $\Omega_{c}\Delta t\rightarrow 0$. To the first point, we show that an
analogous and physically reasonable modification to current deposition
preserves exact energy conservation in implicit PIC schemes. To the second, we
describe a method of adaptively choosing $n_{g}$ based on local field
conditions such that it tends to one in the small time-step limit, and we show
that this single gyrophase sample can be chosen in such a way that second-
order accuracy is preserved.

In the development of this FLR capturing scheme, we have also discovered a
technique for relaxing previously derived time-step constraints on the scheme
[33]. Recall from the reference that, in certain regimes, the scheme’s time-
step is constrained by

$\Omega_{c}\Delta t\lesssim\min\left\\{\left(\rho\frac{\|(\nabla
B)_{\perp}\|}{B}\right)^{-1/2},\sqrt{\Omega_{c}\tau_{res}}\right\\},$ (1)

where $\tau_{res}$ is the smallest time-scale one wishes to resolve in the
simulation. While this does permit stepping over the gyration time-scale, the
presence of square roots means the scheme still requires smaller time-steps
than a truly gyrokinetic particle push. We show here that, by alternating a
large ($\Omega_{c}\Delta t\gg 1$) time-step followed by a smaller one
($\Omega_{c}\delta t\sim 1$), we can relax these time-step constraints to:

$\Omega_{c}\Delta t\lesssim\min\left\\{\left(\rho\frac{\|(\nabla
B)_{\perp}\|}{B}\right)^{-1},\Omega_{c}\tau_{res}\right\\}.$ (2)

In strongly magnetized plasmas, both arguments in the minimum function are
much larger than one, so removing the square roots represents a considerable
potential speedup that is well worth the approximate factor of two penalty
paid for taking a smaller time-step at every other step.

We selectively compare the new scheme (as a test-particle scheme, not coupled
with PIC at this point) against the Boris push, the present scheme’s precursor
[33], the scheme of Brackbill, Forslund and Vu (BFV) [3, 42], and the standard
Crank-Nicolson scheme in a variety of field configurations. We find that the
modifications do enable accurate prediction of single-particle orbits also
when $k_{\perp}\rho\sim 1$, making it the first full-orbit AP scheme to do so
while stepping over the gyration time-scale. We also show that the alternating
time-step approach enables accurate results at larger time-steps than
previously achievable.

The remainder of the article is structured as follows. In Section 2, we review
the previous AP scheme including its motivation, derivation, time-step
restrictions, and relationship to implicit PIC. In Section 3, we motivate and
present the new scheme. We show how to modify implicit PIC methods to be
compatible with both the new scheme and energy conservation, and we derive
new, relaxed time-step restrictions on the scheme and an adaptive time-
stepping strategy that respects those restrictions. We present numerical
results in Section 4 and conclude in Section 5.

##  2 Background

###  2.1 Prior energy-conserving and asymptotic-preserving time integrator

The authors’ previous energy-conserving and asymptotic-preserving (AP) scheme
[33] is built upon Crank-Nicolson while borrowing ideas from [3, 23, 42]. The
scheme is defined by (using the reformulation of the effective force proposed
in Ref. [28]):

$\begin{split}\mathbf{x}^{n+1}&=\mathbf{x}^{n}+\Delta t\mathbf{v}^{n+1/2},\\\
\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}+\mathbf{F}_{eff}\right],\\\
\mathbf{F}_{eff}&=\mathbf{v}\times\left(\frac{\mathbf{G}^{n+1/2}\times\mathbf{v}_{\perp}}{v_{\perp}^{2}}\right)\\\
\mathbf{G}^{n+1/2}&=F_{B,\parallel}\mathbf{b}+\frac{(\mathbf{I}-\widehat{\mathbf{v}}_{E}\widehat{\mathbf{v}}_{E})\cdot\mathbf{F}_{B,\perp}}{1-\eta^{2}/2}+\frac{2}{\eta^{2}}\widehat{\mathbf{v}}_{E}\left(\widehat{\mathbf{v}}_{E}\cdot\mathbf{F}_{B,\perp}+\mathbbm{1}_{v_{E}>u}F_{B,\parallel}\frac{v_{\parallel}}{v_{\perp}}\right),\\\
\mathbf{F}_{B}&=-\tilde{\mu}\nabla
B^{n+1/2},\qquad\tilde{\mu}=m\frac{\left\|\mathbf{v}_{\perp}^{n+1}-\mathbf{v}_{\perp}^{n}\right\|^{2}}{8B^{n+1/2}}.\end{split}$
(3)

In the above, quantities at the half-step are defined by

$\begin{split}\mathbf{v}^{n+1/2}&=\left(\mathbf{v}^{n}+\mathbf{v}^{n+1}\right)/2,\quad\mathbf{x}^{n+1/2}=\left(\mathbf{x}^{n}+\mathbf{x}^{n+1}\right)/2,\\\
\mathbf{E}^{n+1/2}&=\mathbf{E}\left(\mathbf{x}^{n+1/2},t^{n+1/2}\right),\quad\mathbf{B}_{p}^{n+1/2}=\mathbf{B}\left(\mathbf{x}^{n+1/2},t^{n+1/2}\right).\end{split}$
(4)

Moreover,
$\mathbf{v}_{E}=\mathbf{v}_{E}^{n+1/2}=\mathbf{E}^{n+1/2}\times\mathbf{B}^{n+1/2}/(B^{n+1/2})^{2}$,
$\widehat{\mathbf{v}}_{E}=\mathbf{v}_{E}/v_{E}$,
$\mathbf{b}=\mathbf{b}^{n+1/2}=\mathbf{B}^{n+1/2}/B^{n+1/2}$,
$\mathbf{u}=\mathbf{u}^{n+1/2}=\mathbf{v}_{\perp}^{n+1/2}-\mathbf{v}_{E}$, and
$\eta=\min\left\\{1,u/v_{E}\right\\}$. All velocities lacking a superscript
are evaluated at time $t^{n+1/2}$. The subscripts $\perp$ and $\parallel$
refer to vector components orthogonal and parallel to $\mathbf{b}$,
respectively. Finally, $\mathbbm{1}$ denotes the indicator function - that is,
$\mathbbm{1}_{A}$ is one if $A$ is true and zero otherwise.

The full derivation of the scheme is found in [33], but we briefly summarize
the motivation and overarching approach here. If one were to set
$\mathbf{F}_{eff}=0$, the first two lines above are just Crank-Nicolson (CN).
Without modification, CN achieves the correct gyroradius for arbitrary $\Delta
t$ (the reasoning is the same as that used to conclude that the Boris
integrator features an artifically large gyroradius for $\Omega_{c}\Delta t\gg
1$, see e.g. [2, 31]). Moreover, the analysis of Boris in [31] can be
straightforwardly generalized to Crank-Nicolson to show that it also captures
all first-order guiding center drifts except for the $\nabla B$ drift and
mirror force.

To extend Crank-Nicolson to capture the missing $\nabla B$ force, Brackbill,
Forslund and Vu (henceforth BFV) [3, 42] proposed introducing the effective
force $\mathbf{F}_{B}$. This force has the properties:

$\begin{split}\mathbf{F}_{B}\approx-\mu\nabla
B&\quad\text{for}\quad\Omega_{c}\Delta t\gg 1,\\\ \mathbf{F}_{B}=O\left(\Delta
t^{2}\right)&\quad\text{for}\quad\Omega_{c}\Delta t\ll 1.\end{split}$ (5)

Here, $\mu=mu^{2}/2B$ is the particle’s magnetic moment. As a result, this
force is effective in capturing the missing drift motions when
$\Omega_{c}\Delta t\gg 1$ and preserves the second-order convergence of the
scheme to full orbit dynamics as $\Delta t\rightarrow 0$. The drawback is that
$\mathbf{F}_{B}$ can do work on the particle, since there is no guarantee that
$\mathbf{v}^{n+1/2}\cdot\mathbf{F}_{B}=0$. This can lead to $O(1)$ changes in
energy that have significant deleterious effects on the accuracy of the
computed orbits; see the numerical examples in [33].

The scheme (3) improves upon BFV by enforcing energy conservation,
$\mathbf{v}^{n+1/2}\cdot\mathbf{F}_{eff}=0$, in addition to capturing all
first-order drifts for $\Omega_{c}\Delta t\gg 1$. Energy conservation is
achieved for arbitrary $\mathbf{G}$ by the third line of (3). It can then be
shown that the form of $\mathbf{G}$ given is the unique expression satisfying
both

$\begin{split}\left\langle\mathbf{F}_{eff}\right\rangle&=\mathbf{F}_{B},\\\
\left\langle\mathbf{G}\right\rangle&=\mathbf{G},\end{split}$ (6)

where $\langle\cdot\rangle$ denotes the gyro-average operator. The first line
above ensures the AP properties we desire, while the second merely makes the
algebraic derivation tractable. The expression for $\mathbf{G}$ is valid under
the assumption that the particle motion is dominated jointly by gyration and
the $\mathbf{E}\times\mathbf{B}$ drift. It is assumed that if this is not the
case, the gyro-period must be resolved, meaning that
$\mathbf{F}_{eff}=O\left(\Delta t^{2}\right)$ is negligible in any case.

A key observation that will be leveraged in our subsequent analysis is that
$\mathbf{F}_{eff}$ cannot be made to approximate $-\mu\nabla B$ pointwise in
time while simultaneously conserving energy. This is clear in the case when
gyration dominates the particle motion: $\mathbf{v}^{n+1/2}$ rotates rapidly
at the gyrofrequency, while $\mu\nabla B$ varies over much slower time-scales,
so one cannot enforce $\mathbf{v}^{n+1/2}\cdot\mu\nabla B=0$ at all times.
This is what forces us to settle in (6) for imposing constraints on the gyro-
average of $\mathbf{F}_{eff}$ rather than $\mathbf{F}_{eff}$ itself. The
scheme relies on a sequence of consecutive time-steps to sample the particle
at various gyrophases to accumulate an approximate gyro-average of
$\mathbf{F}_{eff}$. This idea motivates the use of alternating large and small
time-steps presented below, as this can force the particle to traverse a
prescribed gyrophase at each time-step.

Clearly, the reliance on time-stepping to compute gyro-averages introduces a
new numerical time-scale that must be resolved in addition to physical scales.
Namely, the time it takes to accumulate an accurate gyro-average should be
small compared to time-scales of interest. Additionally, on shorter time-
scales the particle may drift due to the fact that $\mathbf{F}_{eff}$ deviates
from its desired gyro-average. One should ensure that the distance the
particle drifts due to this unphysical effect is small compared to length-
scales of interest. These constraints were analyzed quantitatively in [33] to
arrive at upper bounds on $\Omega_{c}\Delta t$ based on physical parameters in
two asymptotic limits: $u^{n+1/2}\gg v_{E}$ and $u^{n+1/2}\ll v_{E}$. These
results are summarized in Table 1. Therein, $\tau_{\textrm{res}}$ is a
shorthand for the fastest time-scale one wishes to resolve, and

$\delta_{\perp}=\rho\frac{\left\|(\nabla
B)_{\perp}\right\|}{B},\qquad\delta_{\parallel}=\frac{v_{\parallel}}{\Omega_{c}}\frac{\left\|(\nabla
B)_{\parallel}\right\|}{B},\qquad\delta_{E}=\frac{v_{E}}{\Omega_{c}}\frac{\left\|(\nabla
B)_{\perp}\right\|}{B}.$ (7)

These three non-dimensional parameters have simple interpretations:
$\delta_{\perp}$ is the fractional change in $B$ over a gyro-orbit;
$\delta_{\parallel}$ is the fractional change in $B$ over a gyration time due
to parallel motion; $\delta_{E}$ is the fractional change in $B$ over a
gyration time due to $\mathbf{E}\times\mathbf{B}$ motion. As before, the
expressions in Table 1 are valid under the assumption that the particle’s
velocity is jointly dominated by gyration and $\mathbf{E}\times\mathbf{B}$
motion.

Table 1: The maximum allowable value of $\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 $\tau_{\textrm{res}}$. Expressions are valid in the asymptotic limits (i) $u^{n+1/2}\gg v_{E}$ and (ii) $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\left\\{\sqrt{2}(\delta_{\perp})^{-1/2},(\delta_{\parallel})^{-1/2}\right\\}$ | $2\left(\frac{\Omega_{c}\tau_{\textrm{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_{\textrm{res}}$

###  2.2 Implicit PIC: Energy and local charge conservation

The scheme described above was designed to function as a piece of the implicit
PIC schemes developed over the last decade [6, 9, 10, 11, 12, 13, 34, 39]. The
essence of these schemes can be stated compactly in the simple collisionless,
electrostatic case:

$\begin{split}\mathbf{x}^{n+1}_{p}&=\mathbf{x}^{n}_{p}+\mathbf{v}_{p}^{n+1/2},\\\
\mathbf{v}^{n+1}_{p}&=\mathbf{v}^{n}_{p}+\frac{q_{p}}{m_{p}}\left[\mathbf{E}^{n+1/2}_{p}+\mathbf{v}^{n+1/2}_{p}\times\mathbf{B}^{n+1/2}_{p}\right],\\\
\mathbf{E}_{p}^{n+1/2}&=\sum_{i}\mathbf{E}^{n+1/2}_{i}S\left(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{i}\right),\\\
\mathbf{E}_{i}^{n}&=-\nabla_{h}\phi_{i}^{n},\\\
\nabla_{h}^{2}\phi_{i}^{n+1}&=\nabla_{h}^{2}\phi_{i}^{n}+\frac{\Delta
t}{\epsilon_{0}}\nabla_{h}\cdot\mathbf{j}^{n+1/2}_{i},\\\
\mathbf{j}_{i}^{n+1/2}&=\frac{1}{\lvert\mathbf{h}\rvert}\sum_{p}\mathbf{v}_{p}^{n+1/2}S\left(\mathbf{x}_{i}-\mathbf{x}_{p}^{n+1/2}\right).\end{split}$
(8)

Here, the subscript $p$ indicates a particle quantity while $i$ indicates a
grid quantity. $S$ is the so-called “shape function” that interpolates between
the grid and particles, and $\nabla_{h}$, $\nabla_{h}\cdot$, $\nabla^{2}_{h}$
are discrete gradient, divergence, and Laplacian operators, respectively.
$\mathbf{B}$ is assumed to be a known function that may be evaluated at
particle locations at will, $\lvert\mathbf{h}\rvert$ denotes the cell volume,
and $\epsilon_{0}$ is the permittivity of free space. Note that the first two
lines are simply the Crank-Nicolson time integrator for particles, while the
fifth line is a discretization of the divergence of Ampere’s equation in the
electrostatic limit.

This implicit discretization of the Vlasov-Ampere system has been carefully
constructed to ensure exact energy conservation for any spatial resolution
$\mathbf{h}$ and time-step $\Delta t$. Stated precisely, one has:

$\left[\sum_{p}\frac{1}{2}m_{p}v_{p}^{2}+\frac{\epsilon_{0}}{2}\lvert\mathbf{h}\rvert\sum_{i}E_{i}^{2}\right]^{n+1}=\left[\sum_{p}\frac{1}{2}m_{p}v_{p}^{2}+\frac{\epsilon_{0}}{2}\lvert\mathbf{h}\rvert\sum_{i}E_{i}^{2}\right]^{n}.$
(9)

The proof of this fact in several different contexts is well-described in the
citations above, but we re-print it here because following its logic instructs
later modifications arising from our new time integrator. The proof reads:

$\begin{split}\frac{1}{2}\sum_{p}m_{p}\left(\left\|\mathbf{v}_{p}^{n+1}\right\|^{2}-\left\|\mathbf{v}_{p}^{n}\right\|^{2}\right)&=\sum_{p}m_{p}\mathbf{v}_{p}^{n+1/2}\cdot\left(\mathbf{v}_{p}^{n+1}-\mathbf{v}_{p}^{n}\right)\\\
\textrm{CN definition}\rightarrow\quad&=\Delta
t\sum_{p}q_{p}\mathbf{v}_{p}^{n+1/2}\cdot\left(\mathbf{E}^{n+1/2}_{p}+\mathbf{v}_{p}^{n+1/2}\times\mathbf{B}_{p}^{n+1/2}\right)\\\
\textrm{cross prod. + def. of }\mathbf{E}_{p}\rightarrow\quad&=\Delta
t\sum_{p}\sum_{i}q_{p}\mathbf{v}_{p}^{n+1/2}\cdot\mathbf{E}^{n+1/2}_{i}S\left(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{i}\right)\\\
\textrm{swap sums + }\mathbf{j}\textrm{ def.}\rightarrow\quad&=\Delta
t|\mathbf{h}|\sum_{i}\mathbf{E}^{n+1/2}_{i}\cdot\mathbf{j}_{i}^{n+1/2}\\\
\mathbf{E}=-\nabla_{h}\phi\rightarrow\quad&=-\Delta
t|\mathbf{h}|\sum_{i}\mathbf{j}^{n+1/2}_{i}\cdot\nabla_{h}\phi_{i}^{n+1/2}\\\
\textrm{integrate by parts}\rightarrow\quad&=\Delta
t|\mathbf{h}|\sum_{i}\left(\nabla_{h}\cdot\mathbf{j}^{n+1/2}_{i}\right)\phi_{i}^{n+1/2}\\\
\textrm{def. of Ampere
update}\rightarrow\quad&=\epsilon_{0}|\mathbf{h}|\sum_{i}\left(\nabla^{2}_{h}\phi_{i}^{n+1}-\nabla^{2}_{h}\phi_{i}^{n}\right)\phi_{i}^{n+1/2}\\\
\textrm{integrate by
parts}\rightarrow\quad&=-\epsilon_{0}|\mathbf{h}|\sum_{i}\left(\nabla_{h}\phi_{i}^{n+1}-\nabla_{h}\phi_{i}^{n}\right)\cdot\nabla_{h}\phi_{i}^{n+1/2}\\\
&=-\frac{\epsilon_{0}}{2}|\mathbf{h}|\sum_{i}\left\\{\left\|\mathbf{E}_{i}^{n+1}\right\|^{2}-\left\|\mathbf{E}_{i}^{n}\right\|^{2}\right\\}.\end{split}$
(10)

Simple rearrangement of the first and last expressions yields (9). Note that
we have assumed the discrete differentiation operators preserve the
integration-by-parts properties enjoyed by their continuum analogues. In
practice, some care is needed in choosing $\nabla_{h}$, $\nabla_{h}\cdot$, and
$\nabla^{2}_{h}$ so that this is true. However, as the subject of this article
is the particle advance, we leave those details to the references above.

Note that in moving from the second line to the third in (10), any effective
force $\mathbf{F}_{eff}$ that is added to Crank-Nicolson satisfying
$\mathbf{v}_{p}^{n+1/2}\cdot\mathbf{F}_{eff}=0$ vanishes. Thus, the previously
developed AP time-integrator (3) may be substituted into implicit PIC for
Crank-Nicolson while preserving discrete energy conservation. Similarly, we
will show that the modifications introduced here will preserve this desirable
property when benign changes to relevant definitions are made.

A second important conservation property that should be respected is local
charge conservation. That is, a discrete continuity equation should be
satisfied. With a second-order, centered finite difference spatial
discretization, this takes the form

$\rho_{i+1/2}^{n+1}=\rho_{i+1/2}^{n}-\Delta
t\nabla_{h}\cdot\mathbf{j}_{i}^{n+1/2},$ (11)

where $\rho$ is a charge density. The combination of Ampere’s law and this
continuity equation enforces Gauss’ law. Thus, enforcement of continuity
allows the discrete Vlasov-Ampere system to avoid the need for so-called
“cleaning” of the field at each time-step to enforce Gauss’ law. Instead,
Gauss’s law is automatically satisfied at every time-step.

There are multiple methods in the literature for enforcing this discrete
continuity equation [11, 14, 19], but the vast majority are based on the
observation that within a single cell it is often possible to find a shape
function $\tilde{S}$ (typically different from the $S$ appearing above, and in
fact it may be necessary for $S$ to be different for different components of
$\mathbf{j}$) such that:

$\tilde{S}(\mathbf{x}_{i+1/2}-\mathbf{x}_{p}^{n+1})-\tilde{S}(\mathbf{x}_{i+1/2}-\mathbf{x}_{p}^{n})+\nabla_{h}\cdot\left[(\mathbf{x}_{p}^{n+1}-\mathbf{x}_{p}^{n})S(\mathbf{x}_{i}-\mathbf{x}_{p}^{n+1/2})\right]=0.$
(12)

This identity often requires leveraging the fact that a centered finite
difference is exact for quadratic functions, which places constraints on the
choices of $S$ and $\tilde{S}$ that can be made. Recalling that
$\mathbf{x}_{p}^{n+1}-\mathbf{x}_{p}^{n}=\Delta t\mathbf{v}_{p}^{n+1/2}$,
dividing through by $|\mathbf{h}|$ and summing over all particles, one
observes that this is sufficient to guarantee charge conservation with the
appropriate definition of $\rho$:

$\rho^{n}_{i+1/2}=\frac{1}{|\mathbf{h}|}\sum_{p}\tilde{S}(\mathbf{x}_{i+1/2}-\mathbf{x}_{p}^{n}).$
(13)

Of course, one still must reckon with cell crossings, across which shape
functions are non-smooth and (12) is generally violated. One typically deals
with this by subdividing particle trajectories within a time-step into sub-
steps that terminate on cell boundaries, and performing current deposition and
field interpolation directly on these sub-steps. See the references above for
more detail, including demonstrations that the energy conservation derivation
above carries through without meaningful changes due to the sub-steps.

##  3 Finite Larmor Radius Effects and Time-Step Optimization

As discussed in the preceding section, the time-integration scheme (3)
accurately reproduces the particle gyroradius and all guiding center drifts up
to first order for $\Omega_{c}\Delta t\gg 1$. It does this while still
converging to the full-orbit dynamics at second order as $\Delta t\rightarrow
0$ and conserving energy. While this is a promising step, the guiding center
drifts in question are derived under the assumption of infinitesimal
gyroradius. That is, it is assumed that $\rho/L\ll 1$, where $L$ is the
smallest length scale of interest. This is an assumption of the drift-kinetic
limit, but not of the gyro-kinetic limit.

To extend to the gyrokinetic limit, we must allow $\mathbf{E}$ to vary on
spatial scales comparable to $\rho$ in the perpendicular direction. This
spatial scale length is often written in terms of a perpendicular wave number
$k_{\perp}$. In this notation, the case we are concerned with is
$k_{\perp}\rho\sim 1$. While derivations of gyrokinetic Vlasov systems can be
quite complex, the modification to single particle drift motion to handle this
case is relatively simple: wherever $\mathbf{E}$ appears, it should be
replaced by its gyro-average $\langle\mathbf{E}\rangle$. That is, it is no
longer sufficient to evaluate $\mathbf{E}$ at the gyrocenter; it must be
averaged along the particle’s gyro-orbit.

In previous work, we relied on a sequence of consecutive time-steps to compute
a gyro-average. In particular, we used the fact that

$\frac{1}{N}\sum_{n=1}^{N}\mathbf{F}_{eff}^{n+1/2}\approx\left\langle\mathbf{F}_{eff}\right\rangle.$
(14)

A potential strategy, then, to evaluate $\langle\mathbf{E}\rangle$ is to again
rely on a sequence of consecutive time-steps. The challenge would then be to
arrange the scheme so that

$\frac{1}{N}\sum_{n=1}^{N}\mathbf{E}^{n+1/2}\approx\left\langle\mathbf{E}\right\rangle.$
(15)

This is considerably more difficult than what was previously done with
$\mathbf{F}_{eff}$. There, no field contained structure on the scale of the
gyro-radius, so it sufficed to allow the time-stepping to sample an
uncontrolled sequence of gyrophases. When $\mathbf{E}$ varies on scales
comparable to $\rho$, the gyro-average must be performed with considerably
more care if accuracy is to be expected. Moreover, the electric field may have
structure in the parallel direction, so sampling distinct gyrophases at
distinct times in scenarios that feature non-trivial parallel particle motion
can lead to further difficulties.

As a result of these considerations, we instead propose to explicitly compute
a gyro-average of $\mathbf{E}$ at each time-step and to use this gyro-average
to push particles rather than the pointwise value of $\mathbf{E}$. However,
care is needed to ensure that exact energy conservation is preserved, as well
as convergence to full-orbit motion. We show how this is accomplished in the
following two subsections. Afterward, we describe our strategy for alternating
large and small time-steps, which results in a more accurate approximate
gyroaverage of $\mathbf{F}_{eff}$ after a smaller number of time-steps $N$.
This results in relaxed time-step constraints, which we derive, in addition to
showing that the numerical gyroradius is unaffected by radical changes in the
sizes of consecutive time-steps.

Armed with new, relaxed time-step constraints, we present an adaptive time-
stepping procedure. This also requires reckoning with FLR effects. Next, we
summarize the scheme and its justification in complex-field configurations.
Potential strategies for maintaining exact local charge conservation in a PIC
context are discussed along with energy conservation, but their full
development and implementation is left to future work.

###  3.1 Electric Field Evaluation

Traditionally, in implicit PIC schemes, the electric field used to push
particles $\mathbf{E}_{p}^{n+1/2}$ is evaluated at the midpoint of the time-
step in both space in time. For single-particle motion with an analytically
prescribed electric field, this is written mathematically as:

$\mathbf{E}^{n+1/2}=\mathbf{E}\left(\mathbf{x}^{n+1/2},t^{n+1/2}\right).$ (16)

In the context of a PIC scheme with $\mathbf{E}$ defined on a spatial grid,
this is written as:

$\mathbf{E}_{p}^{n+1/2}=\sum_{i}\mathbf{E}_{i}^{n+1/2}S\left(\mathbf{x}_{i}-\mathbf{x}_{p}^{n+1/2}\right),$
(17)

where $S$ is the particle shape function, $\mathbf{x}_{i}$ the grid points,
and $\mathbf{E}_{i}^{n+1/2}=(\mathbf{E}_{i}^{n}+\mathbf{E}_{i}^{n+1})/2$ is
the half-step field at the grid points.

The desire to capture FLR effects necessitates a change to this midpoint
evaluation. The simplest way to see this is to note that for $\Omega_{c}\Delta
t\gg 1$ the particle roughly traverses a diameter of the gyro-orbit during a
time-step. Thus, $\mathbf{x}_{p}^{n+1/2}$ is approximately the gyro-center
position, rather than residing anywhere on the gyro-orbit. An average over
consecutive time-steps will thus not approximate the gyro-average, but rather
be biased in favor of the field at the gyro-center.

We wish instead to approximate the gyroaverage of $\mathbf{E}$ at the half-
step. To do this, one first needs to use the current particle state to
estimate the gyro-center and gyro-radius, $\mathbf{x}_{c}$ and $\rho$,
respectively. We again assume that the particle’s motion is jointly dominated
by gyration and the $\mathbf{E}\times\mathbf{B}$ drift, so to a good
approximation:

$\mathbf{u}=\mathbf{v}_{\perp}-\mathbf{v}_{E},\qquad\mathbf{x}_{c}=\mathbf{x}+\Omega_{c}^{-1}\mathbf{u}\times\mathbf{b},\qquad\rho=\frac{u}{\Omega_{c}}$
(18)

at any given instant.

As the rest of our scheme is time-centered, we choose time-centered
approximations of $\mathbf{x}_{c}$ and $\rho$. Define an estimate of the
gyrocenter location and gyroradius at integer times $t^{n}$ by

$\mathbf{x}_{c}^{n}=\mathbf{x}^{n}+\frac{1}{\Omega_{c}^{n}}\mathbf{u}^{n}\times\mathbf{b}^{n},\qquad\rho^{n}=\frac{u^{n}}{\Omega_{c}^{n}}$
(19)

and at half-steps by

$\mathbf{x}_{c}^{n+1/2}=\frac{1}{2}\left(\mathbf{x}_{c}^{n}+\mathbf{x}_{c}^{n+1}\right),\qquad\rho^{n+1/2}=\frac{1}{2}\left(\rho^{n}+\rho^{n+1}\right).$
(20)

We then define the estimated gyro-averaged $\mathbf{E}$-field as:

$\left\langle\mathbf{E}\right\rangle^{n+1/2}=\frac{1}{n_{g}}\sum_{k=1}^{n_{g}}\mathbf{E}\left(\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{k},t^{n+1/2}\right),$
(21)

where the $\mathbf{n}_{k}$ are unit vectors such that
$\mathbf{n}_{k}\cdot\mathbf{b}^{n+1/2}=0$ and
$\mathbf{n}_{k+1}=\mathcal{R}_{2\pi/n_{g}}\mathbf{n}_{k}$. Here,
$\mathcal{R}_{\theta}$ is the rotation operator about $\mathbf{b}$ by angle
$\theta$. Clearly, $n_{g}$ is the number of gyrophase samples used.
Additionally, it is clear from Eq. (17) how to generalize this statement to an
electric field defined on a mesh in a self-consistent PIC scheme.

It is immediately evident that this evaluation approximates the gyro-average
of the electric field. However, it is less clear how to guarantee that this
modification preserves second-order convergence to the full-orbit motion in
the small time-step limit. We achieve this in two stages. First, we adaptively
choose $n_{g}$ such that it tends to one as $\Delta t\rightarrow 0$. Second,
we choose the first unit vector $\mathbf{n}_{1}$ intelligently so that
$\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{1}\approx\mathbf{x}^{n+1/2}$.

The first step is achieved by choosing $n_{g}$ adaptively according to:

$n_{g}=\left\lceil\min\left\\{a(k_{\perp}\rho)^{c},b\Omega_{c}\Delta
t\right\\}\right\rceil,$ (22)

with $\lceil~{}\rceil$ the ceiling function, and $a$, $b$ and $c$ free
parameters to be determined later. In this way, if we are either taking small
time-steps or the spatial structure is on scales much larger than the
gyroradius, we choose $n_{g}=1$. We use increasing values of $n_{g}$ as the
time-step gets larger and the spatial scales become more comparable to $\rho$.
We discuss the specification of $a$, $b$, and $c$ in more detail in the
adaptive time-stepping and numerical results sections, as well as the
estimation of $k_{\perp}$ from local field data.

It just remains to choose the proper initial unit normal vector
$\mathbf{n}_{1}$ such that in the limit of small $\Delta t$ and $n_{g}=1$, we
have:

$\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{1}=\mathbf{x}^{n+1/2}+O(\Delta
t^{2}).$ (23)

This will ensure that the scheme is still convergent to the full-orbit
solution at second order as $\Delta t\rightarrow 0$. The last equation
suggests that
$\mathbf{n}_{1}\propto\mathbf{x}_{c}^{n+1/2}-\mathbf{x}^{n+1/2}$, and
therefore we choose:

$\mathbf{n}_{1}=\frac{\mathbf{x}^{n+1/2}-\mathbf{x}_{c}^{n+1/2}}{\left\|\mathbf{x}^{n+1/2}-\mathbf{x}_{c}^{n+1/2}\right\|}=-\frac{\mathbf{u}^{n+1/2}\times\mathbf{b}^{n+1/2}}{u^{n+1/2}},$
(24)

which manifestly has unit magnitude and is orthogonal to $\mathbf{b}$, as
required. Moreover, elementary algebra and the definitions in (20) reveal that

$\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{1}=\mathbf{x}^{n+1/2}+\left(\mathbf{u}^{n+1/2}\times\mathbf{b}^{n+1/2}\right)\Omega_{c}^{-1}\left(1-\frac{u^{n+1}+u^{n}}{2u^{n+1/2}}\right).$
(25)

Note that, when considering the small time-step limit, we may ignore the
effect of $\mathbf{F}_{eff}$ when seeking second-order convergence, as it is
second-order in this limit by construction – see (5)-(6), and [33] for more
details. Thus, up to second order, the update for $\mathbf{u}$ is simply
Crank-Nicolson, which satisfies
$\mathbf{u}^{n+1}=\mathcal{R}_{\theta}\mathbf{u}^{n}$, with the rotation angle
$\theta$ given by [33]

$\cos\theta=\frac{1-\Omega_{c}^{2}\Delta t^{2}/4}{1+\Omega_{c}^{2}\Delta
t^{2}/4}.$ (26)

Thus,

$\left(\frac{u^{n+1}+u^{n}}{2u^{n+1/2}}\right)^{2}=\left(u^{n}\right)^{2}\frac{4}{\left\|R_{\theta}\mathbf{u}^{n}+\mathbf{u}^{n}\right\|^{2}}=\frac{2}{1+\cos\theta}=1+\Omega_{c}^{2}\Delta
t^{2}/4.$ (27)

As a result,

$1-\frac{u^{n+1}+u^{n}}{2u^{n+1/2}}=1-\sqrt{1+\Omega_{c}^{2}\Delta
t^{2}/4}=-\Omega_{c}^{2}\Delta t^{2}/8+O(\Delta t^{4}).$ (28)

This, along with (25), confirms that indeed $\mathbf{x}^{n+1/2}$ and
$\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{1}$ differ by $O(\Delta
t^{2})$ in the small step limit, so second-order accuracy of the overall
scheme is preserved as a consequence of the second-order accuracy of Crank-
Nicolson itself.

###  3.2 Energy and Charge Conservation

Considerable care has been taken in developing implicit PIC methods of the
type (8) to ensure that they feature exact energy and local charge
conservation. It may initially appear that the redefinition of
$\mathbf{E}_{p}^{n+1/2}$ proposed here [see (21)] breaks these properties.
However, we show next that, with analogous changes to the current and charge
deposition processes, energy and charge conservation can be recovered. For the
moment, our proof of charge conservation is limited to the case in which
$n_{g}$ is independent of the time-step, but future work will seek to
generalize this.

To recover energy conservation, we make the following modification to the
current deposition definition in (8):

$\mathbf{j}_{i}^{n+1/2}=\frac{1}{|\mathbf{h}|}\sum_{p}\mathbf{v}_{p}^{n+1/2}\frac{1}{n_{g}}\sum_{k=1}^{n_{g}}S\left(\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{k}-\mathbf{x}_{i}\right).$
(29)

Along with the definition of $\mathbf{E}_{p}^{n+1/2}$ in (51), the logic in
(10) carries through unchanged to recover exact energy conservation. Note that
this definition of current density makes particular physical sense in the
gyrokinetic context, as it distributes current along the gyro-orbit. Moreover,
in the large time-step limit, $\mathbf{v}_{p}^{n+1/2}$ is dominated not by the
particle gyration velocity [at the half-step the gyration velocity vanishes in
the large step limit; see [33] and (62)], but by drift velocities. Thus, to
leading order for large time-steps, Eq. (29) approximates the current density
due to drift motion distributed over the gyro-orbit, as physically desired.

Turning to charge conservation, we recall that what is required is adherence
to some discrete continuity equation along the lines of (11), and that this is
typically achieved by finding a shape function satisfying (12) within a cell
and decomposing particle displacements into substeps that lie within a cell.
We go into more detail here on strategies for division into substeps. In
previous implicit PIC schemes, one either directly chooses time-step sizes so
that particles stop on cell boundaries [11] or one chooses time-steps freely
and interprets particle trajectories (for the purposes of current deposition
and $\mathbf{E}$ evaluation) within a time-step as composed of segments, each
of which lies within a cell [14, 19, 40, 41].

The former strategy will prove clearly incompatible with the method of
alternating large and small time-steps presented below. Indeed, we will
require detailed control over time-step size in order to maximize efficiency.
We will thus focus the remainder of our discussion on the latter strategy.
Note that if we have found an $\tilde{S}$ such that (12) holds for any pair of
points $\mathbf{x}_{p}^{n+1}$, $\mathbf{x}_{p}^{n}$ which lie in a single
cell, then the same relation holds equally well for the pair of points
$\mathbf{x}_{c}^{n+1}+\rho^{n+1}\mathbf{n}_{k}$,
$\mathbf{x}_{c}^{n}+\rho^{n}\mathbf{n}_{k}$. Also note that, critically, the
point at which we evaluate the current’s shape function is the mean of these
two points. Summing the resulting version of (12) from $k=1$ to $n_{g}$ and
dividing by $n_{g}$ as well as $|\mathbf{h}|$, one finds that the discrete
continuity equation is indeed satisfied with the definition

$\rho_{i+1/2}^{n}=\frac{1}{|\mathbf{h}|}\sum_{p}\frac{1}{n_{g}}\sum_{k=1}^{n_{g}}\tilde{S}(\mathbf{x}_{c}^{n}+\rho^{n}\mathbf{n}_{k}-\mathbf{x}_{i+1/2}).$
(30)

However, the logic above implicitly assumes common definitions of $n_{g}$ and
$\mathbf{n}_{k}$ at time-step $n$ and $n+1$. Our proposal to adaptively choose
$n_{g}$ based on local field conditions is inconsistent with this assumption.
The topic of exact charge conservation with adaptive gyro-averaging thus
represents a challenge for future work in applying this scheme to self-
consistent PIC simulations.

###  3.3 Alternating Time-steps

In developing the FLR-capturing scheme above, the authors attempted several
other solutions that, while ultimately unsatisfactory for realistic field
configurations, led to a useful technique for relaxing time-step restrictions
on the scheme that we describe here. We begin in a simplified scenario with
fixed, spatially uniform $\mathbf{E}$ and $\mathbf{B}$. We will argue later
that our conclusions extend to more general field configurations. We recall
from the analysis in [33] – which borrows from ideas in [31] – that the
perpedicular velocity update for Crank-Nicolson (and thus for (3) in uniform
$\mathbf{B}$) may be written

$\mathbf{v}_{\perp}^{n+1}=\mathbf{v}_{E}+\mathcal{R}_{\theta}\left[\mathbf{v}_{\perp}^{n}-\mathbf{v}_{E}\right],$
(31)

where $\mathbf{v}_{E}$ and $\mathcal{R}_{\theta}$ are as defined above, and
the rotation angle $\theta$ is again defined by (26).

In this simple circumstance, the gyration velocity $\mathbf{u}$ is exactly
$\mathbf{v}_{\perp}-\mathbf{v}_{E}$, so this expression may be shortened to

$\mathbf{u}^{n+1}=\mathcal{R}_{\theta}\mathbf{u}^{n}.$ (32)

We now propose to alternate between a large time-step $\Delta t$ and a smaller
one $\delta t$. Let these time-steps give rise to gyration velocity rotation
angles $\theta$ and $\delta\theta$, respectively. Without loss of generality,
let $t^{n+1}-t_{n}=\Delta t$ and $t^{n+2}-t^{n+1}=\delta t$. It is
straightforward to see, geometrically, that the angle between
$\mathbf{u}^{n+1/2}$ and $\mathbf{u}^{n+3/2}$ is $(\theta+\delta\theta)/2$.
This observation is independent of the ordering of the large and small steps,
so we see that the half-step gyration velocity rotates by
$(\theta+\delta\theta)/2$ at every time-step.

It is this simple observation that motivates our alternating-step approach.
Note that as a consequence of (26), for $\Omega_{c}\Delta t\gg 1$ we have
$\theta\approx\pi$. When stepping over the gyration time-scale by a large
factor, a particle thus oscillates between two nearly diametrically opposed
gyrophases at each time-step, meaning it only very slowly accumulates a
representative sample of phases from the entire orbit. Noting, however, that
the half-step velocity is that which actually updates the particle’s position
[see (3)] we see that alternating time-steps allows us some control over the
displacement rotation angle, and thus over which gyrophases are sampled. While
the scheme now explicitly computes a gyro-average of $\mathbf{E}$, it still
relies on a sequence of consecutive time-steps to approximate a gyro-average
of $\mathbf{F}_{eff}$. If we can accumulate a representative gyrophase
sampling more quickly by controlling gyrophase traversal, one expects this to
allow us to take larger time-steps.

We describe how to exercise this control to achieve arbitrary numbers of
equispaced gyrosamples in what follows. First, though, it is important to
confirm that this new alternating-step strategy does not break the desirable
properties enjoyed by (3). To that end, recall that the geometric derivations
in [2, 31, 33] of the numerical gyroradius of the Boris and Crank-Nicolson
integrators assume a fixed time-step. Indeed, they assume the particle
displacement vector has fixed magnitude at each step. This assumption is
clearly violated by the proposed alternating time-step scheme, so it is not
obvious that Crank-Nicolson (and thus the AP integrator (3) above) still
recovers the correct gyroradius when alternating between large and small time-
steps.

Admittedly, it would be surprising if the numerical gyroradius deviated
significantly from the physically correct gyroradius as a result of varying
time-step size. Nevertheless, it is a valuable exercise to show that, in fact,
the exact correct gyroradius is still recovered for an arbitrary sequence of
time-steps. This is shown in Lemma 1 below, and the tools used in this proof
will be instrumental in deriving new time-step restrictions on the scheme.

######  Lemma 1.

Let the particle position and velocity be updated according to

$\begin{split}\mathbf{x}^{n+1}&=\mathbf{x}^{n}+\Delta
t_{n}\mathbf{u}^{n+1/2},\\\ \mathbf{u}^{n+1}&=\mathbf{u}^{n}+\frac{q}{m}\Delta
t_{n}\mathbf{u}^{n+1/2}\times\mathbf{B},\end{split}$ (33)

with $\mathbf{B}$ constant in space and time, and for an arbitrary sequence of
time-steps $\Delta t_{n}$. Further, let $\mathbf{u}^{0}\cdot\mathbf{B}=0$.
Then all particle locations $\mathbf{x}^{n}$ lie on a single circle of radius
$\rho=u/\Omega_{c}$, where $u$ is the common magnitude of all the vectors
$\mathbf{u}^{n}$.

###### Proof.

See Appendix A. ∎

It is also important to understand the relationship of the numerical
gyrocenter to the analytic gyrocenter. Indeed, simply having the correct
gyroradius is insufficient to guarantee the correct gyro-orbit. To that end,
we note a surprising consequence of Lemma 1: Under the evolution (33), the
first three particle locations
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$ define a unique circle.
According to Lemma 1, this circle has radius $\rho$ and all other future
particle locations lie on it. Thus, in this simple scenario, the gyro-orbit –
including the gyro-center – is completely specified by the first two time-
steps, independent of the size of all future time-steps!

This simple observation is the basis for the following corollary.

######  Corollary 1.

Let particle position and velocity again be updated according to (33) with
$\mathbf{u}^{0}\cdot\mathbf{B}=0$. As before, assume that the sequence $\Delta
t_{n}$ is arbitrary and $\mathbf{B}$ is independent of $\mathbf{x}$ and $t$.
Then the particle’s gyrocenter – and thus the entire gyro-orbit – is recovered
exactly in the two limits

  1. (i)

$\displaystyle\Omega_{c}\max\\{\Delta t_{0},\Delta t_{1}\\}\rightarrow 0$,

  2. (ii)

$\displaystyle\Omega_{c}\Delta t_{0}\rightarrow\infty$.

###### Proof.

See Appendix B. ∎

We thus see that if one wishes to recover the gyro-orbit accurately, one only
needs to place a constraint on the initial time-step(s). One may either
initialize a simulation with two very small time-steps satisfying
$\Omega_{c}\max\\{\Delta t_{0},\Delta t_{1}\\}\ll 1$, or one may insist on a
very large initial step satisfying $\Omega_{c}\Delta t_{0}\gg 1$.

Having confirmed that we do not sacrifice our ability to capture gyro-motion
by varying $\Delta t$ from time-step to time-step, we return to our strategy
of alternating large steps $\Delta t$ with smaller steps $\delta t$. As
already noted, this results in the half-step gyration velocity rotating by
angle $\psi=(\theta+\delta\theta)/2$ at every time-step. Under pure gyration,
the angle connecting
$\mathbf{x}^{n}\rightarrow\mathbf{x}^{n+1}\rightarrow\mathbf{x}^{n+2}$ is thus
$\pi-\psi$ and is inscribed in the circle of radius $\rho$ on which these
three particle positions lie. By the inscribed-angle theorem, the arc
connecting $\mathbf{x}^{n}\rightarrow\mathbf{x}^{n+2}$ subtends a central
angle $2(\pi-\psi)$. This state of affairs is illustrated in Figure 1.

Figure 1: Illustration of the relationship between half-step velocity rotation
angle between two time-steps - i.e. $\psi$ \- and the gyrophase angle
subtended by the resulting particle positions.

We find it intuitive to think in terms of asking the particle to sample a
number $m$ of equispaced gyrophases over a sequence of time-steps. This $m$
need not be the same as $n_{g}$. Indeed, it is not even necessary that $m$ be
an integer. Ultimately, we will choose $m$ to make the time-step constraints
that result from this alternating-step strategy as loose as possible. Using
the above results, we see that if we wish for the time-step particle locations
to sample $m$ evenly spaced points around the orbit, we must require
$2(\pi-\psi)m=2\pi$. In this way, we traverse the entire circle over $m$ time-
steps. Solving this expression for $\delta\theta$ gives:

$\delta\theta=\left(1-\frac{1}{m}\right)2\pi-\theta.$ (34)

Recall that the large time-step $\Delta t$ sets the large rotation angle
$\theta$. Given that, the formula above tells how to pick the small time-step.
Indeed, (26) may be inverted to find:

$\Omega_{c}\delta t=2\sqrt{\frac{1-\cos\delta\theta}{1+\cos\delta\theta}},$
(35)

so the angle $\delta\theta$ uniquely specifies the small time-step value. It
is instructive to substitute (34) into (35), yielding:

$\Omega_{c}\delta
t=2\sqrt{\frac{1-\cos\left(\theta+2\pi/m\right)}{1+\cos\left(\theta+2\pi/m\right)}}.$
(36)

One can in theory use trigonometric identities and (26) to directly write
$\Omega_{c}\delta t$ in terms of $\Omega_{c}\Delta t$, but the expression is
cumbersome and unenlightening. However, we note that in the limit
$\Omega_{c}\Delta t\gg 1$, we have $\theta\approx\pi$, and therefore:

$\Omega_{c}\delta t\approx
2\sqrt{\frac{1+\cos\left(2\pi/m\right)}{1-\cos\left(2\pi/m\right)}}.$ (37)

Taylor-expanding the cosine, we find:

$\Omega_{c}\delta t\approx\frac{2}{\pi}m,$ (38)

which is reasonably accurate for $m\geq 8$. While these limits are useful in
building intuition, we use the exact expression (Eq. 36) when implementing the
scheme.

###  3.4 Time-step Constraints

As ever, the ability to step over the gyration time-scale does not permit
arbitrarily large time-steps. An understanding of the restrictions on time-
step size is crucial in practice to avoid accepting erroneous results as
legitimate. In this section, we extend the results on time-step restrictions
developed for the previous scheme (3) in [33]. These prior results were
summarized in Table 1.

Here, we will restrict ourselves to the case $u^{n+1/2}\gg v_{E}$ when
$k_{\perp}\rho\sim 1$. This is so because we require $v_{E}\Delta t\ll
k_{\perp}^{-1}\approx\rho$ to resolve the $\mathbf{E}\times\mathbf{B}$ drift
motion, and $\Omega_{c}\Delta t\gg 1$ implies $u^{n+1/2}\Delta t\approx
2\rho$, and therefore $u^{n+1/2}\gg v_{E}$. Indeed, the very assumption that
the gyro-orbit is circular (rather than elliptical) implies that the gyration
velocity is large compared to the $\mathbf{E}\times\mathbf{B}$ drift velocity.
Promising steps to relax this assumption have been taken in [27], but we only
consider here the standard gyrokinetic approximation. In cases with
$k_{\perp}\rho\ll 1$, the time-step constraints previously derived in [33] in
the limit $u^{n+1/2}\ll v_{E}$ still apply.

As in [33], there are two constraints one must obey if there is to be hope of
accurate integration. First, the time required to accumulate an accurate gyro-
average of $\mathbf{F}_{eff}$ must be small compared to time-scales of
interest, denoted $\tau_{res}$. Second, the anomalous displacement due to
imperfect gyro-averaging on shorter time-scales should be small compared to
length scales of interest, denoted here by $L_{res}$.

The mathematical formulation of the first constraint is considerably
simplified for the new scheme due to the alternating large and small time-
steps. The time required to accumulate an accurate gyro-average is simply
$m(\Delta t+\delta t)/2$. We thus simply have the constraint:

$m\frac{\Delta t+\delta t}{2}<\tau_{res}.$ (39)

Of course, if $\Omega_{c}\tau_{res}\gg 1$, this still permits
$\Omega_{c}\Delta t\gg 1$ for modest $m$ values. Indeed, isolating
$\Omega_{c}\Delta t$ and using the approximate value of $\delta t$ from (38),
this constraint may be simplified to

$\Omega_{c}\Delta t<\frac{2}{m}\Omega_{c}\tau_{res}-\frac{2}{\pi}m.$ (40)

More complex is the derivation of time-step restrictions due to anomalous
displacements arising from short-time-scale errors in gyro-averaging. As
before, such anomalous displacements arise from the fact that the $\nabla B$
force only approximates the physically correct value in a gyro-averaged sense,
and thus deviates from it on short time-scales, leading to anomalous
displacement. Again following [33], we recall that the $\nabla B$ drift has
the physical velocity

$\mathbf{v}_{\nabla B}=\frac{\mathbf{b}}{\Omega_{c}}\times\frac{\mu}{m}\nabla
B,$ (41)

while the AP scheme (3) features the drift velocity

$\mathbf{v}^{\text{AP}}_{\nabla
B}=\frac{\mathbf{b}}{\Omega_{c}}\times\frac{1}{m}\left(\mu_{\text{eff}}^{\text{CN}}\nabla
B-\mathbf{F}_{eff}\right),$ (42)

where $\mu_{\text{eff}}^{\text{CN}}=\mu/(1+\Omega_{c}^{2}\Delta t^{2}/4)$. The
scheme thus features an anomalous displacement with velocity

$\Delta\mathbf{v}_{\nabla B}=\mathbf{v}_{\nabla
B}-\mathbf{v}^{\text{AP}}_{\nabla B}.$ (43)

By virtue of the scheme’s construction, $\left\langle\Delta\mathbf{v}_{\nabla
B}\right\rangle=0$ to a good approximation, but it may deviate from zero on
shorter time-scales. In Appendix C, we show that this displacement may be
analyzed in the presence of alternating large and small time-steps using the
tools developed in Appendix A. The end result of the analysis in Appendix C is
another constraint on time-step: if the anomalous displacement’s magnitude is
to be bounded above by $L_{res}$, then one requires

$\Omega_{c}\left(\frac{\Delta t+\delta
t}{2}\right)<\min\left\\{\frac{L_{res}}{\rho}\frac{\sin\left(2\pi/m\right)}{\delta_{\perp}},\left[\frac{L_{res}}{\rho}\frac{\sin\left(\pi/m\right)}{\delta_{\parallel}}\right]^{1/2}\right\\}.$
(44)

Typically, it will be the case that $L_{res}\approx k_{\perp}^{-1}$.
Additionally, $\delta_{\perp}$ and $\delta_{\parallel}$ are as defined in (7).
Thus, we note that we are still permitted to step over the gyration time-scale
even when $k_{\perp}\rho\sim 1$ so long as $\delta_{\perp}$ and
$\delta_{\parallel}$ are small. Again, if we assume $\Omega_{c}\Delta t\gg 1$,
we can use our approximation of $\Omega_{c}\delta t$ to write a restriction on
the large time-step $\Delta t$ alone:

$\Omega_{c}\Delta
t<2\min\left\\{\frac{L_{res}}{\rho}\frac{\sin\left(2\pi/m\right)}{\delta_{\perp}},\left[\frac{L_{res}}{\rho}\frac{\sin\left(\pi/m\right)}{\delta_{\parallel}}\right]^{1/2}\right\\}-\frac{2}{\pi}m.$
(45)

Note that the restrictions (39) and (44) each represent an improvement over
the previously derived restrictions in Table 1. Previously, time-steps were
bounded above by multiples of $\sqrt{\Omega_{c}\tau_{res}}$ and
$1/\sqrt{\delta_{\perp}}$, while now they are bounded by multiples of
$\Omega_{c}\tau_{res}$ and $1/\delta_{\perp}$. This arises directly from the
improved gyro-averaging, which for very large time-steps shortens the number
of steps required to accumulate a representative sample of gyrophases. In
typical, strongly magnetized applications, we expect
$\Omega_{c}\tau_{res},\delta_{\perp}^{-1}\gg 1$, so this improvement can lead
to significant computational gains.

Of course, we still find that the time-step is constrained to be
$O(1/\sqrt{\delta_{\parallel}})$. Thus, in situations in which parallel motion
represents the active constraint on time-step size, no improvement has been
achieved. However, it is typically the case in fusion plasmas that
characteristic length scales in the parallel direction are much longer than
those in the perpendicular direction. One thus expects it to be somewhat rare
that resolution of parallel motion represents a more stringent constraint than
perpendicular motion, and there is thus reason to hope that the relaxed
constraints in the perpendicular direction will be applicable to the majority
of particles in typical simulations.

###  3.5 Adaptive Time-Stepping and Gyro-averaging

It is clearly desirable to select the largest possible time-step allowed by
the constraints above, so that the orbit may be captured accurately with
minimal computational cost. In this section, we present a slight modification
of the adaptive time-stepping scheme of [33]. The modifications are motivated
by two considerations: (a) the relaxed time-step constraints of the new scheme
derived in the section above, (b) the simplifying assumption that
$v_{E}^{n+1/2}\ll u^{n+1/2}$ is, as argued above, necessary for stepping over
the gyroperiod when $k_{\perp}\rho\sim 1$. Future work will seek to derive a
more general adaptive time-stepping strategy that smoothly transitions between
gyrokinetic, drift-kinetic (in which it may be true that $u^{n+1/2}\leq
v_{E}$), and unmagnetized regimes. Additionally, the adaptive gyro-averaging
in the new scheme requires an estimate of a local value of $k_{\perp}$ in
order to choose the number of gyrophase samples $n_{g}$. There are many
possible methods to achieve this, but we choose a particularly simple one
below.

In principle, adaptive time-step selection is a simple process. One estimates
local length scales, both in the parallel and perpendicular directions, as
well as particle velocity in both directions. Dividing a characteristic length
by a velocity gives a characteristic time. One then chooses a time-step that
is small compared to both parallel and perpendicular characteristic times, and
also respects the time-step constraints derived in the previous section. The
challenge, of course, is in the details.

First, we require an estimate of $k_{\perp}$ to use (22) to specify a a number
of gyro-average samples. The simplest method would be to estimate $k_{\perp}$
via $\|\nabla_{\perp}\mathbf{E}\|/E$, but this of course fails at local
extrema of the electric field. We thus use the somewhat more robust estimate

$k_{\perp}=\max\left\\{\frac{\|\nabla_{\perp}\mathbf{E}\|}{E},\left(\frac{\|D^{2}\mathbf{E}\|}{E}\right)^{1/2}\right\\},$
(46)

where $\|D^{2}\mathbf{E}\|$ is defined by

$\|D^{2}\mathbf{E}\|^{2}=\|\partial_{x}^{2}\mathbf{E}\|^{2}+\|\partial_{y}\mathbf{E}\|^{2},$
(47)

where $x$ and $y$ here are Cartesion coordinates orthogonal to $\mathbf{B}$ at
the particle location. In our code, $x$ is oriented in the
$\mathbf{E}\times\mathbf{B}$ direction, and $y$ is in the
$(\mathbf{E}\times\mathbf{B})\times\mathbf{B}$ direction. Intuitively, this
procedure uses both first and second derivatives to estimate $k_{\perp}$ and
trusts whichever estimate is more pessimistic. Of course, there are numerous
methods to estimate characteristic perpendicular length scales in
$\mathbf{E}$, each with benefits and drawbacks. Implementation of this method
in self-consistent PIC codes in the future will shed additional light on
appropriate choices.

We now consider the three free parameters in Eq. (22). Note that when the
nondimensional step size $\Omega_{c}\Delta t\leq b^{-1}$, we are guaranteed
$n_{g}=1$. One expects to need no gyro-averaging when the gyroperiod is
reasonably well resolved, so we choose $b^{-1}=1/2$. For large time-steps, $a$
sets the number of gyro-samples we use when $k_{\perp}\rho=1$. We find $a=16$
samples to be reasonably accurate for these parameters. The parameter $c$ sets
the sensitivity of $n_{g}$ to variations in $k_{\perp}$. Because our estimate
of $k_{\perp}$ may vary somewhat from time-step to time-step, we find it
desirable to keep $n_{g}$ more consistent in size by choosing $c=1/2$.
Finally, we impose the absolute maximum value of $n_{g}$ to be 64.

We next discuss how we estimate characteristic length scales. We choose

$\begin{split}L_{res,\parallel}&=\Gamma\min\left\\{\frac{B}{|\nabla_{\parallel}B|},\kappa^{-1},\frac{v_{E}}{|\nabla_{\parallel}v_{E}|}\right\\},\\\
L_{res,\perp}&=\Gamma\min\left\\{\frac{B}{\|\nabla_{\perp}B\|},\frac{v_{E}}{\|\nabla_{\perp}v_{E}\|}\right\\}.\end{split}$
(48)

Here, $\kappa=\|(\mathbf{b}\cdot\nabla)\mathbf{b}\|$ is the local magnetic
field curvature. $v_{E}$ is computed using a gyro-averaged electric field, so
this estimate does not explicitly include any gyro-scale structure of
$\mathbf{E}$. It is expected that such structure is captured by the adaptive
gyro-averaging procedure. $\Gamma$ is a free parameter that governs what
fractional change in each quantity we are willing to tolerate within a time-
step. We follow [33] by choosing $\Gamma=0.1$.

Pseudo-code for the adaptive time-stepping procedure is as follows:

  1. 1.

Compute $k_{\perp}$, $n_{g}$, $L_{res,\parallel}$, and $L_{res,\perp}$ as
described above.

  2. 2.

Let

$\tau_{res}=\min\left\\{\frac{L_{res,\perp}}{v_{E}},\frac{L_{res,\parallel}}{v_{\parallel}}\right\\}.$
(49)

  3. 3.

Set $F=(L_{res,\perp}/\rho)\sin\left(2\pi/m\right)$, and use time-step

$\Delta
t=2\Omega_{c}^{-1}\alpha\min\left\\{\frac{F}{\delta_{\perp}},\sqrt{\frac{F}{\delta_{\parallel}}},\frac{\Omega_{c}\tau_{res}}{m}\right\\}$
(50)

As in [33], $\alpha$ is a parameter that controls how close to the time-step
restrictions we are willing to approach. As before, we use $\alpha=0.9$.

Clearly, this procedure is only used to pick the larger time-steps $\Delta t$
in the alternating time-step procedure. Note that if $\Delta t$ is
sufficiently small, it is actually possible that the “smaller” time-step
$\delta t$ defined by (36) and (26) can be larger than $\Delta t$, and thus
inadmissible. For this reason, at the subsequent time-step, we will choose
$\min\\{\Delta t,\delta t\\}$, where $\Delta t$ is specified using the
procedure above.

In this initial proof-of-principle study, we additionally specify an absolute
cap on the size of $\Omega_{c}\Delta t$. In our tokamak test case below, this
cap is set to 70, but results are not substantially impacted by small changes
to this parameter.

###  3.6 Algorithm Summary and Complex Field Configurations

In this section, we outline the proposed algorithm and provide an intuitive
justification for its use in complex field configurations. The overall
structure of the procedure is as follows:

  1. 1.

Given $\mathbf{x}^{n}$, $\mathbf{v}^{n}$, choose an appropriate time-step
$\Delta t$ and number of gyrophase samples $n_{g}$ using the procedure in the
previous subsection.

  2. 2.

Advance in time by one step of size $\Delta t$ using the scheme (3), but with
the modification that

$\mathbf{E}^{n+1/2}=\frac{1}{n_{g}}\sum_{k=1}^{n_{g}}\mathbf{E}\left(\mathbf{x}_{c}^{n+1/2}+\rho^{n+1/2}\mathbf{n}_{k},t^{n+1/2}\right).$
(51)

  3. 3.

Now, having $\mathbf{x}^{n+1}$, $\mathbf{v}^{n+1}$, choose the next time-step
$\delta t$ satisfying

$\Omega_{c}^{n+1}\delta
t=2\sqrt{\frac{1-\cos\left(\theta+2\pi/m\right)}{1+\cos\left(\theta+2\pi/m\right)}}$
(52)

where $\theta$ is defined by $\cos\theta=(1-(\Omega_{c}^{n})^{2}\Delta
t^{2}/4)/(1+(\Omega_{c}^{n})^{2}\Delta t^{2}/4)$, and
$\Omega_{c}^{n}=qB(\mathbf{x}^{n})/m$.

  4. 4.

Advance by one step of size $\delta t$, again using (3) with the modification
(51).

  5. 5.

$n\leftarrow n+2$ and return to step 1.

We note that in our derivations we took advantage of the fact that the
gyration velocity $\mathbf{u}$ is exactly $\mathbf{v}-\mathbf{v}_{E}$. This is
still true to good approximation in varying fields, for all other drift
motions are small in the magnetized regime. When these drift motions grow to a
substantial size (as can easily be tested at each step based on local field
quantities), the gyro-motion must be resolved in any case.

One should choose $m$ so that time-steps can be taken as large as possible on
average. Even casual investigation of (39) and (44) reveals that the optimal
$m$ that allows the largest possible average time-step depends on which
constraint is active. We take a simplified approach and uniformly choose
$m=5$.

##  4 Numerical Results

The new scheme is tested in four configurations. The first features constant
$\mathbf{B}$, with $\mathbf{E}$ varying sinusoidally in one coordinate. This
case is chosen for the ability to compare with asymptotic theory predicting
the $\mathbf{E}\times\mathbf{B}$ drift velocity as a function of
$k_{\perp}\rho$ for $k_{\perp}\rho\lesssim 1$. In the second case, we let
$\mathbf{E}$ vary sinusoidally in both coordinates and confirm that the
alternating time-step methodology permits enlarged time-steps. The third case
features linearly varying $\mathbf{B}$ with $\mathbf{E}$ again varying
sinusoidally in both perpendicular coordinates, although this time with
$k_{\perp}\rho>1$. The fourth case uses a tokamak geometry with Solov’ev
equilibrium $\mathbf{B}$ field with a rapidly varying perpendicular
$\mathbf{E}$, chosen to test the scheme’s ability to capture parallel dynamics
in conjunction with complex geometry and FLR effects.

For clarity, we work in a dimensionless formulation in which $q_{p}/m_{p}=1$
(implying $B=\Omega_{c}$). To measure the accuracy of the present scheme
(which we term CBFV+FLR), we compare with asymptotic theory when available and
with the results of the widely-used Boris integrator using $\Omega_{c}\Delta
t=0.1$. We also compare against other AP integrators such as the scheme of
Brackbill, Forsluand and Vu (BFV) [3, 42], the present scheme’s precursor
(CBFV) [33], and the standard Crank-Nicolson (CN) scheme.

The systems of equations resulting from our implicit integrator are solved
using the preconditioned Picard iteration introduced in [28] with nonlinear
tolerance set to $10^{-12}$. In the first three cases, the particle is
initialized in the $x$-$y$ plane at

$\mathbf{x}(t=0)=\frac{1}{B}\left(\begin{array}[]{c}\cos\omega\\\
-\sin\omega\end{array}\right),\qquad\mathbf{v}(t=0)=\left(\begin{array}[]{c}-\sin\omega\\\
-\cos\omega\end{array}\right)$ (53)

for some $\omega$, and the $B$ above is evaluated at the origin. This is done
so that the particle’s initial gyrocenter is located at the origin. While we
do not report them here in the interest of brevity, all tests are performed
for multiple values of $\omega$ to confirm the gyrophase-independence of the
results.

The code is implemented in Python, and beyond the use of the efficient solver
of [28], no effort has been made at performance optimization in this proof-of-
principle study. This has two noteworthy consequences. Firstly, direct
computation-time comparisons between the new method and Boris (or any other
method) are not to be trusted if one wishes to extrapolate to optimized,
compiled code running on HPC systems. We thus do not report specific numbers,
instead only noting that our implicit implementation is typically between
10$\times$ and 120$\times$ faster than well-resolved Boris in the experiments
performed, lending confidence that considerable speed-ups are possible. Note
that in practice, speed-ups will be problem-dependent, and in particular
depend on the non-dimensional parameters appearing in (7).

Secondly, the exact energy conservation properties the scheme provably enjoys
in the implicit-PIC context do not carry over to single particle motion in
prescribed analytic fields. Indeed, the Crank-Nicolson scheme on which it is
based only conserves energy to $O(\Delta t^{2})$ in analytically prescribed
fields, and our modifications do not change that fact. In [33], this was
corrected using the scheme of [35], in which the electric field is modified
according to

$\mathbf{E}^{n+1/2}\rightarrow\frac{\phi^{n}-\phi^{n+1}}{\left(\mathbf{x}^{n+1}-\mathbf{x}^{n}\right)\cdot\mathbf{E}^{n+1/2}}\mathbf{E}^{n+1/2}$
(54)

to recover exact energy conservation, where $\phi$ is the electrostatic
potential associated with $\mathbf{E}$. However, in the highly oscillatory
electric fields used here, the exactly-conservative scheme has tremendous
difficulty converging without sophisticated preconditioning. As such, when
measuring energy conservation, we will be satisfied with observing that energy
errors are small and do not feature secular growth in time. Future work will
empirically confirm whether exact conservation is realized when coupling to an
implicit PIC solver.

###  4.1 1-D Sinusoidal Variation featuring asymptotic theory

We set $\mathbf{B}=100\widehat{\mathbf{z}}$, with

$\mathbf{E}=\cos(k_{\perp}y)\widehat{\mathbf{y}}.$ (55)

The particle is initialized with $\omega=\pi/2$. Asymptotic theory [16]
predicts that the $\mathbf{E}\times\mathbf{B}$ drift velocity in this scenario
is given by

$\mathbf{v}_{E}\approx\mathbf{v}_{E,gc}\left(1-\frac{k_{\perp}^{2}\rho^{2}}{4}\right),$
(56)

where $\mathbf{v}_{E,gc}$ is the drift velocity computed with $\mathbf{E}$ at
the gyrocenter, i.e., the drift-kinetic limit. This expression, being derived
via Taylor expansion in $k_{\perp}\rho$, is valid for $k_{\perp}\rho\lesssim
1$.

As the particle is initialized with its gyrocenter at the origin (where
$\mathbf{E}=\widehat{\mathbf{y}}$) we have
$\mathbf{v}_{E,gc}=\widehat{\mathbf{x}}/100$. We simulate to a final time
$T=100$, so that the predicted final $x$-position of the particle is simply
$x(T)=1-k_{\perp}^{2}\rho^{2}/4$. This asymptotic theory result is compared in
Figure 2 to Boris, CBFV, and the new scheme (CBFV+FLR). For this simple test,
we do not use adaptive time-stepping or adaptive gyro-averaging. The larger
time-step is fixed at $\Omega_{c}\Delta t=100$, with alternating steps chosen
in the manner described in Section 3.4. We show results using fixed values of
$n_{g}$ of 8 and 4.

Figure 2: Final $x$-position of the particle resulting from Boris (blue)
integration, the new scheme with $n_{g}=8$ gyrophase samples (CBFV+FLR, red),
new scheme with $n_{g}=4$ (CBFV+FLR, teal), asymptotic theory (black), and a
prior AP scheme with no FLR corrections (CBFV, orange). We see the ability of
the new scheme to match asymptotic theory when $k_{\perp}\rho\lesssim 1$, as
well as good agreement with fully-resolved Boris out to even large
$k_{\perp}\rho$ when 8-point gyroaveraging is used.

One immediately sees agreement with the asymptotic theory in its range of
validity, $k_{\perp}\rho<1$. Outside that range, the well-resolved Boris
simulation may be regarded as the ground truth. We thus also see excellent
agreement of new scheme well beyond $k_{\perp}\rho=1$ when $n_{g}=8$. As is to
be expected, more accurate gyro-averaging is necessary to improve agreement
for large $k_{\perp}\rho$. Also as expected, the prior scheme without FLR
corrections fails to predict any change in the drift velocity as a function of
$k_{\perp}$.

###  4.2 Alternate Time-stepping test

Next, we perform a similar test to the last one, but focusing on the benefits
of the proposed strategy of alternating large and small time-steps. We use
fields with locally similar qualitative properties but different parameters:

$\mathbf{B}=\left(100+3\sin
x\right)\widehat{\mathbf{z}},\qquad\mathbf{E}=\frac{1}{10}\left(\cos
k_{x}x+\sin k_{y}y\right).$ (57)

We choose $k_{x}=5\pi$ and $k_{y}=1$. We initialize the particle at the origin
with $\mathbf{v}^{0}=\widehat{\mathbf{x}}$. These parameters result in
$k_{\perp}\rho\approx 0.15$, which is an optimal regime to demonstrate the
benefits of the alternate-stepping strategy.

Because we are specifically interested in pushing the limits of the time-step
size in the alternating-step method, for this test we do not use the adaptive
time-stepping scheme described in Section 3.5. Instead, we manually choose the
large time-step as $\Omega_{c}\Delta t=800$, while the small step $\delta t$
is chosen as described in Section 3.5. When not using the alternating-step
technique, we fix the time-step at $\Omega_{c}\Delta t=400$ so that the the
two methods have a roughly equal average time-step.

We plot the resulting orbits in the $x$-$y$ plane as well as their time
histories in Figure 3. It is immediately evident that the alternating-step
approach dramatically improves the accuracy of the orbit, as predicted.
Indeed, the correct orbit is periodic, while without alternating steps the
orbit is non-periodic on the time-scales observed.

Figure 3: Depiction of orbit and time-histories for alternating-time-step test
problem. For comparable average time-steps, the alternating-step method
reproduces the correct closed, elliptic orbit, while the non-alternating
approach does not.

###  4.3 Sinusoidal Variation in Perpendicular Plane with $\nabla B$ drift

For this test, we set

$\mathbf{B}=100(1-x/20)\widehat{\mathbf{z}},\qquad\mathbf{E}=\frac{1}{2}\cos(k_{x}x)\widehat{\mathbf{x}}+\cos(k_{y}y)\widehat{\mathbf{y}}.$
(58)

While several values of $k_{x},k_{y}$ were tested, we report here results for
$k_{x}\rho=\pi/4$, $k_{y}\rho=\pi/6$, where $\rho$ here is the initial
gyroradius of $1/100$. These wavenumbers are chosen because they lead to a
complex orbit arising from competition between the $\nabla B$ drift (in the
positive $y$-direction) and $\mathbf{E}\times\mathbf{B}$ drift (in the
negative $y$-direction), leading to a more challenging test case.

For this test, we also compare against BFV, CBFV and CN. Results appear in
Figure 4. As expected, CN completely fails to see the $\nabla B$-drift. BFV
and CBFV see the $\nabla B$-drift, but incorrectly predict the confinement in
the $x$-direction. Only CBFV+FLR gives orbits that match Boris. Lastly, in
Figure 5, we display energy conservation properties of each scheme. As
discussed above, only $O(\Delta t^{2})$ energy conservation is expected in the
single-particle context, so it is unsurprising that each scheme has larger
energy errors than Boris purely as a result of the much larger time-steps. It
is encouraging, however, to observe that the new modifications presented here
do not degrade energy conservation in the single-particle context relative to
previous schemes.

Figure 4: Depiction of orbit for 2-dimensional problem with sinusoidal, gyro-
scale variation in electric field. The FLR corrections presented here are
necessary to reproduce the fully-resolved orbit. Figure 5: Fractional energy
change as a function of time for two-dimensional problem with sinusoidal,
gyro-scale variation in electric field. Exact conservation is not expected
here, but the new scheme does not feature secular energy growth and conserves
energy no worse than previous schemes.

###  4.4 Tokamak Geometry

As a final test, we consider a tokamak-like geometry with fixed electric and
magnetic fields. The setup is identical to that in [33] except for the
addition of an electric field with gyration-scale fluctuations, but we
reproduce the description here for completeness. The poloidal magnetic field
is based on a simple analytic solution $\psi$ of the Grad-Shafranov equation
[5]:

$\psi(r,z)=\frac{C}{8}r^{4}+d_{1}+d_{2}r^{2}+d_{3}\left(r^{4}-4r^{2}z^{2}\right),$
(59)

where $C,d_{i}$ are constants and $(r,\theta,z)$ is the standard polar
coordinate system.

The flux function $\psi$ specifies the poloidal field by
$\mathbf{B}_{p}=\nabla\psi\times\mathbf{e}_{\theta}/r$. We let $C=300$, and we
take advantage of the fact, shown in [32], that there is a one-to-one linear
map between the $d_{i}$ and the tokamak shape parameters $\varepsilon$
(inverse aspect ratio), $\kappa_{T}$ (elongation), and $\delta$
(triangularity). Thus, the poloidal field is uniquely specified by choosing
the shape parameters of the International Thermonuclear Energy Reactor (ITER):
$\varepsilon=0.32$, $\kappa_{T}=1.7$, $\delta=0.33$. The toroidal field is
chosen to be $\mathbf{B}_{tor}=800\mathbf{e}_{\theta}/r$.

We additionally impose an electrostatic potential that is a flux function
(i.e. expressible as a function of $\psi$) such that
$\mathbf{E}\cdot\mathbf{B}=0$ according to the ideal MHD Ohm’s law,
$\mathbf{E}+\mathbf{v}\times\mathbf{B}=0$. That flux function is chosen to be
$\phi=\sin(K\psi)/2K$. $K$ is chosen to target a specific characteristic value
of $k_{\perp}$ in the electric field by computing
$K=k_{\perp}/\|\nabla\psi(\mathbf{x}^{0})\|$.

The particle is initialized at the Cartesian location
$\mathbf{x}^{0}=1.2\widehat{\mathbf{x}}$, at which point the magnitude of the
magnetic field is $B\approx 667.7$. This knowledge allows us to compute the
gyroradius $\rho$ at the initial particle location, and to choose $\kappa$ so
that $k_{\perp}\rho=1.5$ at the initial particle location. The value of
$k_{\perp}\rho$ varies in space, but remains $O(1)$ along the entire particle
trajectory. The initial velocity is given in Cartesian coordinates by
$\mathbf{v}^{0}=(1,v_{y},0)^{T}$. Note that as the field is predominantly
toroidal (i.e., along the $y$-axis) at the initial location, $v_{y}$ may be
taken as a reasonable proxy for the particle’s initial parallel velocity. We
choose $v_{y}=0.6$ because it results in an orbit near the trapped-passing
boundary. Simulations are run to a final time of $T=400$, which corresponds to
several traversals of the banana orbit.

Figure 6: Projection into the poloidal plane of particle orbits for the
tokamak equilibrium test case. Different curves are trajectories resulting
from well-resolved Boris integration (blue), the new scheme with FLR
corrections (red), and the previous scheme without FLR corrections (orange).
Single lines represent flux surfaces of the Solov’ev equilibrium magnetic
field. Full orbit appears at left, with a zoom-in of the upper turning point
on the right. The modified structure of the orbit due to fine-scale variation
in the electric field is clearly visible. In the poloidal projection of the
orbit, the FLR corrections seem to afford only mild improvements in accuracy
in capturing the turning points, although this betrays catastrophic errors in
their timing, as demonstrated in the next figure.

To make plots more readable, we only show results from Boris, CBFV, and
CBFV+FLR for this test (although we have confirmed that the other integrators
fail as catastrophically as in the previous test). Only slight improvements
due to the FLR corrections are visible from a poloidal projection of the
particle orbit, shown in Figure 6. However, without FLR corrections, CBFV
results in the particle traversing this orbit at the wrong rate. This is shown
in Figure 7, in which time-histories of $r$ and $z$ are reported. The rate of
traversal is dramatically improved by the FLR corrections, resulting in only a
slight deviation from fully-resolved Boris over several banana orbits.

Figure 7: Time histories of $r$ and $z$ for the tokamak equilibrium test case.
Note the dramatically improved temporal accuracy of the new integrator with
FLR corrections.

##  5 Conclusions

In this work, we have built on the asymptotic-preserving and energy-conserving
integrator of [33] with two modifications that allow it to capture FLR effects
even when stepping over the gyration time-scale. Those modifications are (a)
introducing an explicit and adaptive gyro-average computation of $\mathbf{E}$,
and (b) a strategy for alternating large and small time-steps so that the
particle samples evenly spaced gyrophases over successive time-steps. We have
shown that with the first modification, the scheme still reproduces the proper
gyro-orbit while the time-step constraints on the scheme are considerably
relaxed due to the second. We have also shown that with an analogous
modification to current deposition, the scheme remains compatible with the
exact energy conservation enjoyed by existing implicit PIC schemes. The scheme
was tested on single-particle motion with several specified electromagnetic
fields, demonstrating the ability of the scheme to accurately predict orbits
even when the electric field varies on scales comparable to the gyroradius.

There are, of course, opportunities for further development. Of obvious
interest is the implementation of this scheme in implicit PIC codes, and its
subsequent use in simulating physically relevant phenomena. Such an
implementation will certainly benefit from a more robust and general adaptive
time-stepping scheme that can smoothly and accurately handle transitions
between magnetized and unmagnetized regimes in the presence of the alternating
time-step approach. Additionally, local charge conservation in the presence of
variable $n_{g}$ requires further investigation. Each of these is an active
area of current research.

## Acknowledgements

The authors wish to acknowledge private communication with Matthew
Miecnikowski, whose input directly led to the idea of alternating large and
small time-steps. Additional discussions with C-S Chang, Guangye Chen, Joshua
Burby, and Davide Curreli were valuable in the development of this work.

##  Appendix A Proof of Lemma 1

It suffices to show that any four consecutive particle locations lie on a
single circle of radius $\rho=u/\Omega_{c}$. This extends to all points by
induction as follows. Suppose
$\\{\mathbf{x}^{k},\mathbf{x}^{k+1},\mathbf{x}^{k+2},\mathbf{x}^{k+3}\\}$ lie
on a circle of radius $\rho$, as do
$\\{\mathbf{x}^{k+1},\mathbf{x}^{k+2},\mathbf{x}^{k+3},\mathbf{x}^{k+4}\\}$.
Then, $\\{\mathbf{x}^{k+1},\mathbf{x}^{k+2},\mathbf{x}^{k+3}\\}$ lie on a
shared circle of radius $\rho$. Three non-colinear points define a unique
circle. Thus, all five points
$\\{\mathbf{x}^{k},\mathbf{x}^{k+1},\mathbf{x}^{k+2},\mathbf{x}^{k+3},\mathbf{x}^{k+4}\\}$
lie on that same unique circle of radius $\rho$, and the induction carries.

Without loss of generality, we proceed by showing that the first four points
$\\{\mathbf{x}^{k}:\,k=0,1,2,3\\}$ lie on a single circle of radius $\rho$. As
already noted, the velocity update in the perpendicular direction may be
written

$\mathbf{u}^{n+1}=R_{\theta_{n}}\mathbf{u}^{n},$ (60)

where $\theta_{n}$ is defined by

$\cos\theta_{n}=\frac{1-\Omega_{c}^{2}\Delta
t_{n}^{2}/4}{1+\Omega_{c}^{2}\Delta t_{n}^{2}/4}.$ (61)

This already demonstrates that the magnitude of $\mathbf{u}^{n}$ is
independent of $n$. We denote this common magnitude by $u$ for the remainder
of the proof.

Moreover, it is a simple algebraic exercise – see, e.g. [33] – to show that

$\mathbf{u}^{n+1/2}=\frac{\mathbf{u}^{n}+\frac{1}{2}\Omega_{c}\Delta
t_{n}\mathbf{u}^{n}\times\mathbf{b}}{1+\Omega_{c}^{2}\Delta t_{n}^{2}/4}.$
(62)

Taking the magnitude of this expression gives

$u^{n+1/2}=\frac{u}{\sqrt{1+\Omega_{c}^{2}\Delta t_{n}^{2}/4}}.$ (63)

One may thus write an expression for the particle displacement magnitude
$\Delta
x^{n}=\left\|\Delta\mathbf{x}^{n}\right\|=\left\|\mathbf{x}^{n+1}-\mathbf{x}^{n}\right\|$
at each time-step:

$\Delta x^{n}=u^{n+1/2}\Delta t_{n}=\frac{u\Delta
t_{n}}{\sqrt{1+\Omega_{c}^{2}\Delta t_{n}^{2}/4}}.$ (64)

Using the trigonometric identity
$\sin\left(\theta_{n}/2\right)=\sqrt{(1-\cos\theta_{n})/2}$ and the definition
of $\theta_{n}$ in (61), we find

$\Delta x^{n}=2\rho\sin\left(\frac{\theta_{n}}{2}\right).$ (65)

Further, as discussed in the main text, the angle between the displacements
$\Delta\mathbf{x}^{n}$ and $\Delta\mathbf{x}^{n+1}$ is necessarily
$(\theta_{n}+\theta_{n+1})/2$. For brevity, we define
$\theta_{n+1/2}=(\theta_{n}+\theta_{n+1})/2$. The law of sines thus gives the
following expression for the radius $r$ of the unique circle containing
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$:

$r=\frac{1}{2}\frac{\left\|\mathbf{x}^{2}-\mathbf{x}^{0}\right\|}{\sin\left(\pi-\frac{\theta_{0}+\theta_{1}}{2}\right)}=\frac{1}{2}\frac{\left\|\mathbf{x}^{2}-\mathbf{x}^{0}\right\|}{\sin\theta_{1/2}}.$
(66)

The numerator may be simplified by applying the law of cosines:

$\begin{split}\left\|\mathbf{x}^{2}-\mathbf{x}^{0}\right\|&=\left[\left(\Delta
x^{0}\right)^{2}+\left(\Delta x^{1}\right)^{2}+2\Delta x^{0}\Delta
x^{1}\cos\theta_{1/2}\right]^{1/2}\\\
&=2\rho\left[\sin^{2}\frac{\theta_{0}}{2}+\sin^{2}\frac{\theta_{1}}{2}+2\sin\frac{\theta_{0}}{2}\sin\frac{\theta_{1}}{2}\cos\theta_{1/2}\right]^{1/2}.\end{split}$
(67)

This may be further simplified by applying an angle sum identity to
$\cos\theta_{1/2}$:

$\begin{split}\left\|\mathbf{x}^{2}-\mathbf{x}^{0}\right\|&=2\rho\left[\sin^{2}\frac{\theta_{0}}{2}+\sin^{2}\frac{\theta_{1}}{2}+2\sin\frac{\theta_{0}}{2}\sin\frac{\theta_{1}}{2}\left(\cos\frac{\theta_{0}}{2}\cos\frac{\theta_{1}}{2}-\sin\frac{\theta_{0}}{2}\sin\frac{\theta_{1}}{2}\right)\right]^{1/2}\\\
&=2\rho\left[\sin^{2}\frac{\theta_{0}}{2}\cos^{2}\frac{\theta_{1}}{2}+\cos^{2}\frac{\theta_{0}}{2}\sin^{2}\frac{\theta_{1}}{2}+2\sin\frac{\theta_{0}}{2}\sin\frac{\theta_{1}}{2}\cos\frac{\theta_{0}}{2}\cos\frac{\theta_{1}}{2}\right]^{1/2}\\\
&=2\rho\left[\left(\sin\frac{\theta_{0}}{2}\cos\frac{\theta_{1}}{2}+\cos\frac{\theta_{0}}{2}\sin\frac{\theta_{1}}{2}\right)^{2}\right]^{1/2}\\\
&=2\rho\left[\sin^{2}\theta_{1/2}\right]^{1/2}\\\
&=2\rho\sin\theta_{1/2}.\end{split}$ (68)

Substituting this into (66), we have $r=\rho$. Thus,
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$ lie on a circle of radius
$\rho$. As we have not used any special properties of the first three points,
it also follows that $\\{\mathbf{x}^{1},\mathbf{x}^{2},\mathbf{x}^{3}\\}$ lie
on a circle of radius $\rho$.

Now, $\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$ and
$\\{\mathbf{x}^{1},\mathbf{x}^{2},\mathbf{x}^{3}\\}$ each lie on a circle of
radius $\rho$ and they share two points. There are thus at most two circles of
radius $\rho$ on which these four points can lie. This is depicted graphically
below in Figure 8.

Figure 8: Illustration of the two possible (solid) circles on which
$\mathbf{x}^{3}$ may lie. In the limit $\Omega_{c}\Delta t\rightarrow\infty$,
$\mathbf{x}^{3}$ must also lie on the dashed circle, so it must lie on one of
the two points of tangency. In the text, it is shown that only the “correct”
tangency point (the leftward one) is possible.

Our aim is to show that $\mathbf{x}^{3}$ lies on the same circle as
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$. To this end, consider
first the limit $\Omega_{c}\Delta t_{2}\rightarrow\infty$. Clearly, according
to (61) and (65), we have $\Delta x^{2}\rightarrow 2\rho$ in this limit. Thus,
in this limit, $\mathbf{x}^{3}$ must lie on one of the two circles of radius
$\rho$ appearing in Figure 8, and on a circle of radius $2\rho$ centered at
$\mathbf{x}^{2}$ – the dashed circle in Figure 8. $\mathbf{x}^{3}$ thus lies
on one of the points of tangency between the circles of radius $\rho$ and that
of $2\rho$. Note now that, in order to lie on the point of tangency located on
the circle of radius $\rho$ that does not contain $\mathbf{x}^{0}$, it must be
the case that the particle velocity rotates in the opposite direction in
moving from $t_{2}$ to $t_{3}$ as it has previously (or to rotate by an angle
larger than $\pi$, which is also impossible). This would require $\mathbf{B}$
to change direction, which is ruled out by assumption.

Thus, at least in the $\Omega_{c}\Delta t_{2}\rightarrow\infty$ limit,
$\mathbf{x}^{3}$ lies on the circle of radius $\rho$ defined by
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$. We again appeal to (61)
and (65) to see that the particle displacement is a smooth function of $\Delta
t_{2}$. For $\mathbf{x}^{3}$ to jump from one of the two radius-$\rho$ circles
in Figure 8 to the other as $\Delta t$ varies would require a discontinuity at
least in the derivative $d\Delta\mathbf{x}^{2}/d\Delta t_{2}$. As no such
discontinuity exists, we conclude that
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2},\mathbf{x}^{3}\\}$ lie on a
single circle of radius $\rho$ independent of $\Delta t_{n}$.

As we have nowhere exploited the fact that $\mathbf{x}^{0}$ is the initial
particle location, the same conclusion holds for
$\\{\mathbf{x}^{n},\mathbf{x}^{n+1},\mathbf{x}^{n+2},\mathbf{x}^{n+3}\\}$.
Thus, the induction carries and the proof is complete. ∎

##  Appendix B Proof of Corollary 1

The first limit (i) is clear. In the small step limit, the orbit converges to
the analytic orbit. The first two time-steps specify
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$, which in turn specify a
unique circle of radius $\rho$ on which all future points lie. Thus, if
$\\{\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2}\\}$ converge to their true
values, then the numerical orbit they specify coincides with the true orbit.

In considering the second limit (ii), we first note that the true gyro-center
is located at

$\mathbf{x}_{gc}=\mathbf{x}_{0}+\Omega_{c}^{-1}\mathbf{u}^{0}\times\mathbf{b}=\mathbf{x}^{0}+\rho\widehat{\mathbf{u}}^{0}\times\mathbf{b},$
(69)

where $\widehat{\mathbf{u}}^{0}=\mathbf{u}^{0}/u^{0}$. Next, we return to (62)
and note that

$\lim_{\Omega_{c}\Delta
t_{0}\rightarrow\infty}\mathbf{u}^{1/2}=2\frac{\mathbf{u}^{0}\times\mathbf{b}}{\Omega_{c}\Delta
t_{0}}=2\frac{\rho}{\Delta t_{0}}\widehat{\mathbf{u}}^{0}\times\mathbf{b}.$
(70)

Thus, the next particle location is

$\lim_{\Omega_{c}\Delta
t_{0}\rightarrow\infty}\mathbf{x}^{1}=\mathbf{x}^{0}+\Delta
t_{0}\mathbf{u}^{1/2}=\mathbf{x}^{0}+2\rho\widehat{\mathbf{u}}^{0}\times\mathbf{b}.$
(71)

We see that, in this limit, $\mathbf{x}^{1}$ and $\mathbf{x}^{0}$ are
connected by a diameter of the true gyro-orbit. Moreover, there is a unique
circle of radius $\rho$ that contains this diameter. Thus, all future particle
locations share this circle, and thus lie on the true gyro-orbit. ∎

##  Appendix C Derivation of Time-Step Restrictions

As noted in the main text, we consider only the case $u^{n+1/2}>v_{E}$. The
expression for $\mathbf{F}_{eff}$ simplifies considerably in this case,
allowing us to write

$\Delta\mathbf{v}_{\nabla
B}=\frac{\mathbf{b}}{\Omega_{c}}\times\frac{\tilde{\mu}}{m}\left\\{\frac{v_{\parallel}}{v_{\perp}}\widehat{\mathbf{v}}_{\perp}(\nabla
B)_{\parallel}+(\mathbf{I}-2\widehat{\mathbf{v}}_{\perp}\widehat{\mathbf{v}}_{\perp})\cdot(\nabla
B)_{\perp}\right\\},$ (72)

where

$\tilde{\mu}=\frac{\Omega_{c}^{2}\Delta t^{2}/4}{1+\Omega_{c}^{2}\Delta
t^{4}/4}\mu$ (73)

for the large time-step, and analogous for the small time-step. We break this
anomalous velocity into two components to be analyzed separately:

$\Delta\mathbf{v}_{\nabla
B}^{1}=\frac{\mathbf{b}}{\Omega_{c}}\times\frac{\tilde{\mu}}{m}\frac{v_{\parallel}}{v_{\perp}}(\nabla
B)_{\parallel}\widehat{\mathbf{v}}_{\perp},\qquad\Delta\mathbf{v}_{\nabla
B}^{2}=\frac{\mathbf{b}}{\Omega_{c}}\times\frac{\tilde{\mu}}{m}\left(\mathbf{I}-2\widehat{\mathbf{v}}_{\perp}\widehat{\mathbf{v}}_{\perp}\right)\cdot(\nabla
B)_{\perp}.$ (74)

As before, quantities without superscripts are assumed to be evaluated at the
half-step.

We will denote the anomalous displacements arising from these anomalous
velocities by $\mathbf{a}^{1}$ and $\mathbf{a}^{2}$, defined respectively by

$\mathbf{a}^{i,n+1}=\mathbf{a}^{i,n}+\Delta t_{n}\Delta\mathbf{v}_{\nabla
B}^{i,n+1/2}$ (75)

for $i=1,2$. In the uniform time-step case of [33], it was noted that these
respective anomalous displacements have constant magnitude and rotate by a
fixed angle at every time-step - $\theta$ for $\mathbf{a}^{1}$, and $2\theta$
for $\mathbf{a}^{2}$. A simple geometric argument then gave formulae for the
radii of the circles traversed by $\mathbf{a}^{1}$ and $\mathbf{a}^{2}$.
Bounding these radii in terms of physical parameters then led immediately to
time-step restrictions.

We follow the same general approach here, but the geometric arguments differ
because the displacement magnitude is no longer constant due to the variation
in time-step size. As such, we will use similar arguments to those used in
Appendix A. To that end, note that the rotation angle of each displacement is
still constant - $\psi$ for $\mathbf{a}^{1}$, and $2\psi$ of $\mathbf{a}^{2}$.
Here $\psi=\pi(1-1/m)$ is as defined in the main text. We will denote
displacements in the large and small time-steps by $\Delta\mathbf{a}^{i}$ and
$\delta\mathbf{a}^{i}$, respectively.

We begin with the second anomalous motion $\mathbf{a}^{2}$. By taking
magnitudes, we have

$\Delta a^{2}=\frac{\Omega_{c}^{2}\Delta t^{2}/4}{1+\Omega_{c}^{2}\Delta
t^{2}/4}\frac{\mu}{m\Omega_{c}}\left\|\left(\nabla
B\right)_{\perp}\right\|\Delta t,\qquad\delta a^{2}=\frac{\Omega_{c}^{2}\delta
t^{2}/4}{1+\Omega_{c}^{2}\delta
t^{2}/4}\frac{\mu}{m\Omega_{c}}\left\|\left(\nabla
B\right)_{\perp}\right\|\delta t.$ (76)

Some algebra yields

$\Delta
a^{2}=\left(1-\cos\theta\right)\frac{\delta_{\perp}}{4}\Omega_{c}\Delta
t\rho,\qquad\delta
a^{2}=\left(1-\cos\delta\theta\right)\frac{\delta_{\perp}}{4}\Omega_{c}\delta
t\rho,$ (77)

where we have used

$\frac{\Omega_{c}^{2}\Delta t^{2}/4}{1+\Omega_{c}^{2}\Delta
t^{2}/4}=\frac{1}{2}\left[1-\frac{1-\Omega_{c}^{2}\Delta
t^{2}/4}{1+\Omega_{c}^{2}\Delta t^{2}/4}\right]=\frac{1-\cos\theta}{2}$ (78)

along with the analogous reasoning for the small time-step. As in Appendix A,
we now use the law of sines to write an expression for the radius $r$ of a
circle containing three consecutive $\mathbf{a}^{2}$ points:

$2r=\frac{\left[\left(\Delta a^{2}\right)^{2}+\left(\delta
a^{2}\right)^{2}+2\Delta a^{2}\delta
a^{2}\cos\left(2\psi\right)\right]^{1/2}}{\lvert\sin\left(2\psi\right)\rvert}.$
(79)

Evidently, $r$ does not depend on time-step number, so we can use identical
logic to Appendix A to conclude that all $\mathbf{a}^{2}$ points lie on a
single circle with this radius. By using trivial bounds on geometric
functions, we can find a simple upper bound on the radius of this circle:

$r\leq\frac{1}{2}\frac{\rho\delta_{\perp}}{\sin\left(\frac{2\pi}{m}\right)}\Omega_{c}\left(\frac{\Delta
t+\delta t}{2}\right).$ (80)

where we have substituted in the definition of $\psi$ in terms of the number
of desired gyrophase samples $m$. If the length scale we wish to resolve is
$L_{res}$, a sufficient condition to make sure the anomalous displacement is
smaller than that scale is thus

$\Omega_{c}\left(\frac{\Delta t+\delta
t}{2}\right)<\left(\frac{L_{res}}{\rho}\right)\frac{\sin\left(\frac{2\pi}{m}\right)}{\delta_{\perp}}.$
(81)

We now return to $\mathbf{a}^{1}$, which is analyzed in the same way. Taking
magnitudes, we have

$\Delta a^{1}=\frac{1}{2}\frac{\Omega_{c}^{3}\Delta
t^{3}/4}{\sqrt{1+\Omega_{c}^{2}\Delta
t^{2}/4}}\delta_{\parallel}\rho,\qquad\delta
a^{1}=\frac{1}{2}\frac{\Omega_{c}^{3}\delta
t^{3}/4}{\sqrt{1+\Omega_{c}^{2}\delta t^{2}/4}}\delta_{\parallel}\rho.$ (82)

The root in the denominator arises from the evaluation of $v_{\perp}$ at the
half-step, and the assumption that perpendicular motion is dominated by
gyration, giving $v_{\perp}\approx u^{n}/\sqrt{1+\Omega_{c}^{2}\Delta
t^{2}/4}$ and similarly for the small step.

Anticipating that we’ll resort to inequalities as before, the expressions for
the displacements can be reduced to upper bounds:

$\Delta x^{1}\leq\frac{1}{4}\Omega_{c}^{2}\Delta
t^{2}\delta_{\parallel}\rho,\qquad\delta
x^{1}\leq\frac{1}{4}\Omega_{c}^{2}\delta t^{2}\delta_{\parallel}\rho.$ (83)

Again substituting into the law of sines with the fact that the rotation angle
is now $\psi$, we have

$r\leq\frac{1}{2}\frac{\rho\delta_{\parallel}}{\sin\left(\frac{\pi}{m}\right)}\left[\Omega_{c}\left(\frac{\Delta
t+\delta t}{2}\right)\right]^{2}.$ (84)

We again requre that $2r<L_{res}$. A sufficient condition for satisfying this
is

$\Omega_{c}\left(\frac{\Delta t+\delta
t}{2}\right)<\left[\left(\frac{L_{res}}{\rho}\right)\frac{\sin\left(\frac{\pi}{m}\right)}{\delta_{\parallel}}\right]^{1/2}.$
(85)

The time-step restriction (44) in the main text simply enforces both (81) and
(85).

## References

  * [1] IG Abel, GG Plunk, E Wang, M Barnes, SC Cowley, W Dorland, and AA Schekochihin.  Multiscale gyrokinetics for rotating tokamak plasmas: fluctuations, transport and energy flows.  Reports on Progress in Physics, 76(11):116201, 2013.
  * [2] C.K. Birdsall and A.B. Langdon.  Plasma physics via computer simulation.  CRC press, 2004.
  * [3] J.U. Brackbill and D.W. Forslund.  Simulation of low-frequency, electromagnetic phenomena in plasmas.  In Multiple time scales, pages 271–310. Elsevier, 1985.
  * [4] Alain J Brizard and Taik Soo Hahm.  Foundations of nonlinear gyrokinetic theory.  Reviews of modern physics, 79(2):421, 2007.
  * [5] A.J. Cerfon and J.P. Freidberg.  “one size fits all” analytic solutions to the Grad–Shafranov equation.  Physics of Plasmas, 17(3):032502, 2010.
  * [6] L. Chacón, G. Chen, and D.C. Barnes.  A charge-and energy-conserving implicit, electrostatic particle-in-cell algorithm on mapped computational meshes.  Journal of Computational Physics, 233:1–9, 2013.
  * [7] Philippe Chartier, Nicolas Crouseilles, Mohammed Lemou, Florian Méhats, and Xiaofei Zhao.  Uniformly accurate methods for vlasov equations with non-homogeneous strong magnetic field.  Mathematics of Computation, 88(320):2697–2736, 2019.
  * [8] Philippe Chartier, Nicolas Crouseilles, Mohammed Lemou, Florian Méhats, and Xiaofei Zhao.  Uniformly accurate methods for three dimensional vlasov equations under strong magnetic field with varying direction.  SIAM Journal on Scientific Computing, 42(2):B520–B547, 2020.
  * [9] G. Chen and L. Chacón.  An energy-and charge-conserving, nonlinearly implicit, electromagnetic 1D-3V Vlasov–Darwin particle-in-cell algorithm.  Computer Physics Communications, 185(10):2391–2402, 2014.
  * [10] G. Chen and L. Chacón.  A multi-dimensional, energy-and charge-conserving, nonlinearly implicit, electromagnetic Vlasov–Darwin particle-in-cell algorithm.  Computer Physics Communications, 197:73–87, 2015.
  * [11] G. Chen, L. Chacón, and D.C. Barnes.  An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm.  Journal of Computational Physics, 230(18):7018–7036, 2011.
  * [12] G. Chen, L. Chacón, and D.C. Barnes.  An efficient mixed-precision, hybrid CPU–GPU implementation of a nonlinearly implicit one-dimensional particle-in-cell algorithm.  Journal of Computational Physics, 231(16):5374–5388, 2012.
  * [13] G. Chen, L. Chacón, C.A. Leibs, D.A. Knoll, and W. Taitano.  Fluid preconditioning for Newton–Krylov-based, fully implicit, electrostatic particle-in-cell simulations.  Journal of computational physics, 258:555–567, 2014.
  * [14] G. Chen, L. Chacón, L. Yin, B.J. Albright, D.J. Stark, and R.F. Bird.  A semi-implicit, energy-and charge-conserving particle-in-cell algorithm for the relativistic vlasov-maxwell equations.  arXiv preprint arXiv:1903.01565, 2019.
  * [15] Y. Chen and S. E. Parker.  Particle-in-cell simulation with Vlasov ions and drift kinetic electrons.  Physics of Plasmas, 16(5):052305, May 2009.
  * [16] D. Curreli.  private communication.
  * [17] S. De, T. Singh, A. Kuley, J. Bao, Z. Lin, G. Y. Sun, S. Sharma, and A. Sen.  Kinetic particle simulations in a global toroidal geometry.  Physics of Plasmas, 26(8):082507, Aug 2019.
  * [18] Daniel HE Dubin, John A Krommes, C Oberman, and WW Lee.  Nonlinear gyrokinetic equations.  The Physics of fluids, 26(12):3524–3535, 1983.
  * [19] T Zh Esirkepov.  Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor.  Computer Physics Communications, 135(2):144–153, 2001.
  * [20] F. Filbet and L.M. Rodrigues.  Asymptotically stable particle-in-cell methods for the vlasov–poisson system with a strong external magnetic field.  SIAM Journal on Numerical Analysis, 54(2):1120–1146, 2016.
  * [21] F. Filbet and L.M. Rodrigues.  Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas.  SIAM Journal on Numerical Analysis, 55(5):2416–2443, 2017.
  * [22] Francis Filbet and Luis Miguel Rodrigues.  Asymptotically preserving particle methods for strongly magnetized plasmas in a torus.  Journal of Computational Physics, 480:112015, 2023.
  * [23] T.C. Genoni, R.E. Clark, and D.R. Welch.  A fast implicit algorithm for highly magnetized charged particle motion.  Open Plasma Physics Journal, 3:36–41, 2010.
  * [24] TS Hahm, WW Lee, and A Brizard.  Nonlinear gyrokinetic theory for finite-beta plasmas.  The Physics of fluids, 31(7):1940–1948, 1988.
  * [25] Ernst Hairer, Christian Lubich, and Yanyan Shi.  Large-stepsize integrators for charged-particle dynamics over multiple time scales.  Numerische Mathematik, 151(3):659–691, 2022.
  * [26] Y. Hu, M. Miecnikowski, Y. Chen, and S. E. Parker.  Fully kinetic simulation of ion-temperature-gradient instabilities in tokamaks.  Plasma, 1(1):105–118, May 2018.
  * [27] Ilon Joseph.  Guiding center and gyrokinetic orbit theory for large electric field gradients and strong shear flows.  Physics of Plasmas, 28(4), 2021.
  * [28] Oleksandr Koshkarov, Luis Chacon, Guangye Chen, and Lee Forrest Ricketson.  Fast nonlinear iterative solver for an implicit, energy-conserving, asymptotic-preserving charged-particle orbit integrator.  Journal of Computational Physics, 459:111146, 2022.
  * [29] WW Lee.  Gyrokinetic approach in particle simulation.  The Physics of Fluids, 26(2):556–562, 1983.
  * [30] R. Pandit, D. Curreli, J. Choi, R. Hager, M. Churchill, and C-S Chang.  Conversion from gyrokinetic to full-orbit particles for plasma sheath boundary conditions.  In APS Meeting Abstracts, 2018.
  * [31] S..E. Parker and C.K. Birdsall.  Numerical error in electron orbits with large $\omega_{ce}\delta t$.  Journal of Computational Physics, 97(1):91–102, 1991.
  * [32] A. Pataki, A..J Cerfon, J.P. Freidberg, L. Greengard, and M. O’Neil.  A fast, high-order solver for the Grad–Shafranov equation.  Journal of Computational Physics, 243:28–45, 2013.
  * [33] L.F. Ricketson and L. Chacón.  An energy-conserving and asymptotic-preserving charged-particle orbit implicit time integrator for arbitrary electromagnetic fields.  arXiv preprint arXiv:1904.09478, 2019.
  * [34] L.F. Ricketson and G. Chen.  A pseudospectral implicit particle-in-cell method with exact energy and charge conservation.  Computer Physics Communications, 291:108811, 2023.
  * [35] J.C. Simo, N. Tarnow, and K.K. Wong.  Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics.  Computer methods in applied mechanics and engineering, 100(1):63–116, 1992.
  * [36] B. J. Sturdevant.  Fully kinetic ion models for magnetized plasma simulations.  Doctoral Dissertation at CU Boulder, 2016.
  * [37] B. J. Sturdevant, Y. Chen, and S. E. Parker.  Low frequency fully kinetic simulation of the toroidal ion temperature gradient instability.  Physics of Plasmas, 24(8):081207, Aug 2017.
  * [38] B. J. Sturdevant, S. E. Parker, Y. Chen, and B. B. Hause.  An implicit $\delta f$ particle-in-cell method with sub-cycling and orbit averaging for lorentz ions.  Journal of Computational Physics, 316:519 – 533, 2016.
  * [39] W.T. Taitano, D.A. Knoll, L. Chacón, and G. Chen.  Development of a consistent and stable fully implicit moment method for Vlasov–Ampère particle in cell (PIC) system.  SIAM Journal on Scientific Computing, 35(5):S126–S149, 2013.
  * [40] T. Umeda, Y. Omura, T. Tominaga, and H. Matsumoto.  A new charge conservation method in electromagnetic particle-in-cell simulations.  Computer Physics Communications, 156(1):73–85, 2003.
  * [41] J. Villasenor and O. Buneman.  Rigorous charge conservation for local electromagnetic field solvers.  Computer Physics Communications, 69(2-3):306–316, 1992.
  * [42] H.X. Vu and J.U. Brackbill.  Accurate numerical solution of charged particle motion in a magnetic field.  Journal of Computational Physics, 116(2):384–387, 1995.
  * [43] Z. Wang, H. Qin, B. Sturdevant, and C-S Chang.  Geometric electrostatic particle-in-cell algorithm on unstructured meshes.  12 2020.