# A conservative implicit-PIC scheme for the hybrid kinetic-ion fluid-electron
plasma model on curvilinear meshes

A. Stanier, L. Chacón  Applied Mathematics and Plasma Physics, Los Alamos
National Laboratory, Los Alamos, NM 87545, United States

###### Abstract

The hybrid kinetic-ion fluid-electron plasma model is widely used to study
challenging multi-scale problems in space and laboratory plasma physics. Here,
a novel conservative scheme for this model employing implicit particle-in-cell
techniques [1] is extended to arbitrary coordinate systems via curvilinear
maps from logical to physical space. The scheme features a fully non-linear
electromagnetic fomulation with a multi-rate time advance - including sub-
cycling and orbit-averaging for the kinetic ions. By careful choice of
compatible particle-based kinetic-ion and mesh-based fluid-electron
discretizations in curvilinear coordinates, as well as particle-mesh
interpolations and implicit midpoint time advance, the scheme is proven to
conserve total energy for arbitrary curvilinear meshes. In the electrostatic
limit, the method is also proven to conserve total momentum for arbitrary
curvilinear meshes. Although momentum is not conserved for arbitrary
curvilinear meshes in the electromagnetic case, it is for an important subset
of Cartesian tensor-packed meshes. The scheme and its novel conservation
properties are demonstrated for several challenging numerical problems using
different curvilinear meshes, including a merging flux-rope simulation for a
space weather application, and a helical $m=1$ mode simulation for magnetic
fusion energy application.

###### keywords:

Hybrid , curvilinear , plasma , particle-in-cell , implicit , conservative ,
space weather , fusion

##  1 Introduction

Magnetized plasmas in nature and the laboratory can be inherently multi-scale,
with extreme scale disparity between the global system sizes and the ion and
electron kinetic orbit scales. From a modeling perspective, the
magnetohydrodynamic (MHD) fluid approach is widely used to model the bulk
plasma flows and the evolution of the magnetic topology at the largest scales
[2], whereas the high-fidelity kinetic Vlasov-Maxwell model is used for
detailed study of plasma micro-scales [3]. However, it is highly desirable to
bridge this scale gap for many systems to enable the study of cross-scale
coupling. This scale-bridging is likely best accomplished through the combined
use of reduced kinetic models (which are of lower fidelity than the full
kinetic model yet capture key physics) along with algorithmic advances to step
over the stiffest space and time scales in a stable manner.

A promising class of reduced models take a hybrid approach, in which some
plasma components are treated with a high-fidelity kinetic description, and
the remainder using a lower-fidelity reduced fluid approximation [4]. Variants
include the kinetic-ion fluid-electron model [5, 6, 7], the fast-particle-
kinetic bulk-fluid model [8, 9, 10], as well as generalized hybrid models that
are flexible to different choices of the kinetic and fluid components [11]. In
this paper, we focus on the kinetic-ion fluid-electron model with the kinetic
ions being solved with the particle-in-cell (PIC) technique [12], and assume a
zero-inertia and scalar pressure electron fluid. This model is widely used to
study the macro- to ion-scale coupling in magnetic reconnection [13, 14, 15],
collisionless shocks [16, 17, 18], collisionless turbulence [19, 20, 21] as
well as global magnetospheric physics [22, 23, 24, 25].

Common numerical schemes for this model employ explicit time-stepping. As
discussed in Ref. [7], a key algorithmic difficulty with explicit schemes
concerns the use of a static Ohm’s law to evaluate the electric field along
with a leap-frog advance for the particles. The time advanced electric field
is needed to advance the particles, but this in turn depends implicitly on the
time advanced particle velocities. To make this explicit, some form of
predictor-corrector [26, 27, 28], moment-based [29, 30], or forward-
extrapolation schemes [31, 32, 13] have been used (see also Refs. [7, 13] for
reviews). These schemes are often highly efficient, but the approximations
made are such that the numerical conservation of momentum and energy is
violated [1]. A further numerical difficulty with explicit schemes concerns
the extremely stiff CFL-condition associated with the quadratic dispersion
relation for whistler waves, namely $\Delta t\propto n(\Delta x)^{2}/B$, where
$\Delta t$ is the timestep, $\Delta x$ the grid spacing, $n$ the density and
$B$ the magnetic field-strength. This can be extremely stiff in strongly
magnetized and near-vacuum regions, and often requires the syb-cycling of the
magnetic field-solve [33, 30, 13] or some other special treatment [34, 26, 35]
to avoid numerical instability.

Implicit methods have been used to solve the electromagnetic field equations
[34, 4] and, more recently, for the kinetic particle advance. In particular,
Ref. [36] explored implicit timestepping for particles featuring sub-stepping
and orbit-averaging techniques in an electrostatic $\delta F$ hybrid model.
More recently, Ref. [1] proposed a fully implicit scheme for the non-linear
electromagnetic kinetic-ion fluid-electron hybrid model, which combines an
implicit-midpoint time advance with a cell-centered finite difference spatial
discretization, and includes particle sub-cycling and orbit averaging. This
latter scheme, which was formulated for uniformly spaced meshes in Cartesian
geometry, features simultaneous discrete conservation of momentum and energy.
It was demonstrated that this scheme also features favorable stability
properties with respect to the finite grid-instability, which can occur in
non-conservative schemes when the ratio of the ion-to-electron temperatures is
small ($T_{i}/T_{e}\ll 1$) [37, 1].

Non-uniform or curvilinear mapped meshes can be used to resolve thin boundary
layer features efficiently, or to model boundary-fitted domains that have
specific geometries. For the hybrid-PIC model, explicit schemes have been
formulated with cylindrical [34, 38], toroidal [39, 36], or spherical meshes
[33, 40] to study space or laboratory plasmas. Adaptive-Mesh Refinement (AMR)
techniques have also been applied to the hybrid-PIC model for global
magnetosphere simulations [41]. The use of non-uniform meshes and/or
curvilinear coordinates within PIC methods can be tricky in general and
careful choices of compatible grid and particle-based discretizations, as well
as particle-grid interpolations, are required to avoid large violations in
momentum and energy conservation. One issue concerns the non-physical particle
self-forces, which can be amplified in regions of strong mesh gradients to
potentially cause artificial particle reflection [42]. A large range of mesh
resolutions across the simulation domain can also cause problems both in terms
of statistical noise – due to a different number of macro-particles in each
cell – as well as by the loss in stability or accuracy in coarse grid regions
for numerical schemes that are not asymptotic preserving. For the latter, it
is well known for standard kinetic Vlasov-Maxwell PIC algorithms [12] that
significant numerical heating can occur if the electron Debye length
$\lambda_{D,e}=v_{\textrm{th},e}/\omega_{pe}$ is not resolved everywhere,
where $v_{\textrm{th},e}$ is the electron thermal speed and $\omega_{pe}$ is
the plasma frequency. For the quasi-neutral hybrid kinetic-ion fluid-electron
PIC model, the effects can be more subtle. It has been recently shown [43]
that significant unphysical wave dispersion errors can occur in hybrid-PIC
when combining the use of higher-order particle-mesh interpolation functions
with under-resolving the ion skin-depth $d_{i}=v_{A}/\Omega_{ci}$, where
$v_{A}$ is the Alfvén speed and $\Omega_{ci}$ is the ion gyro-frequency. These
errors are caused by the particle-mesh interpolations spreading the electric
field that is experienced by the ion particles, and leading to its inexact
cancellation when combining the particle-based ion and grid-based electron
momentum equations. Remarkably, the source of these errors is removed when
using the lowest order nearest grid-point (NGP) interpolation functions.

In this paper, we extend the fully implicit scheme of Ref. [1] to treat
arbitrary coordinate systems via curvilinear maps from logical to physical
space. This novel hybrid-PIC scheme has the following unique conservation
properties: 1) it conserves energy numerically for arbitrary curvilinear maps
(assuming suitable boundary conditions), 2) in the electrostatic limit it also
conserves total linear momentum for general meshes, and 3) for the full
electromagnetic model, it features discrete momentum conservation for
Cartesian meshes with a tensor-product mesh packing function. Although the
latter property only holds for a narrow subset of meshes, it is notable since
the lack of artificial self-forces means that it avoids the issue of spurious
particle reflection mentioned above. Section 2 of the paper reviews the quasi-
neutral hybrid kinetic-ion fluid-electron plasma model. We describe our
implicit-PIC discretization of the model for curvilinear meshes in Section 3,
and prove its discrete conservation properties in Sec. 4. In Section 5, we
detail our particular numerical implementation and solution strategy for the
discrete model equations, and finally, in Section 6, we present numerical
example problems that utilize different mesh geometries and verify the novel
conservation properties of our scheme.

##  2 Electromagnetic kinetic-ion fluid-electron hybrid model

In the electromagnetic and quasi-neutral hybrid kinetic-ion fluid-electron
model, ions with species index $s$ are described by the distribution function
$f_{s}(t,\boldsymbol{v},\boldsymbol{x})$ that is found from the Vlasov
equation as

$\partial_{t}f_{s}+\boldsymbol{\nabla}\cdot\left(f_{s}\boldsymbol{v}\right)+\left(q_{s}/m_{s}\right)\left(\boldsymbol{E}^{*}+\boldsymbol{v}\times\boldsymbol{B}\right)\cdot\boldsymbol{\nabla}_{v}f_{s}=0.$
(1)

Here, $q_{s},m_{s}$ are the ion species charge and mass, respectively,
$\boldsymbol{B}$ is the magnetic field, and $\boldsymbol{E}^{*}$ is the
frictionless electric field (see below). The magnetic field is advanced by
Faraday’s equation as

$\partial_{t}\boldsymbol{B}=-\boldsymbol{\nabla}\times\boldsymbol{E},$ (2)

where $\boldsymbol{E}$ is the total electric field. The latter is evaluated
using an Ohm’s law as

$\boldsymbol{E}=\boldsymbol{E}^{*}+\boldsymbol{F}_{ie}=-\boldsymbol{u}_{i}\times\boldsymbol{B}+\frac{\boldsymbol{j}\times\boldsymbol{B}}{ne}-\frac{\boldsymbol{\nabla}p_{e}}{ne}+\boldsymbol{F}_{ie},$
(3)

where $e$ is the elementary charge, $n=\sum_{s}(q_{s}/e)\int f_{s}d^{3}v$ is
the quasi-neutral density, $\boldsymbol{u}_{i}=\left(\sum_{s}q_{s}\int
f_{s}\boldsymbol{v}d^{3}v\right)/en$ is the ion current carrying drift
velocity, $\boldsymbol{j}=\boldsymbol{\nabla}\times\boldsymbol{B}/\mu_{0}$ is
the current density, $p_{e}$ is a scalar electron pressure,
$\boldsymbol{F}_{ie}=\eta\boldsymbol{j}$ is the resistive friction, and $\eta$
is the plasma resistivity. We note that only the frictionless electric field
$\boldsymbol{E}^{*}=\boldsymbol{E}-\boldsymbol{F}_{ie}$ is used in Eq. (1) as
is required for momentum conservation, see e.g. Ref. [1].

Finally, the scalar electron pressure is calculated from

$\frac{1}{\gamma-1}\left[\partial_{t}p_{e}+\boldsymbol{\nabla}\cdot\left(\boldsymbol{u}_{e}p_{e}\right)\right]+p_{e}\boldsymbol{\nabla}\cdot\boldsymbol{u}_{e}=H_{e}-\boldsymbol{\nabla}\cdot\boldsymbol{q}_{e},$
(4)

where $\gamma=5/3$ is the ratio of specific heats,
$\boldsymbol{u}_{e}=\boldsymbol{u}_{i}-\boldsymbol{j}/en$ is the electron bulk
velocity, $H_{e}=\eta j^{2}$ is the frictional (Joule) heating term,
$\boldsymbol{q}_{e}=-\kappa\boldsymbol{\nabla}(p_{e}/n)$ is the electron heat
flux, and $\kappa$ is the electron heat conductivity.

With suitable choice of boundary conditions, the set of equations (1-4)
conserves the total mass $\left(\int\sum_{s}m_{s}\int
f_{s}d^{3}vd^{3}x\right)$, linear momentum $\left(\int\sum_{s}m_{s}\int
f_{s}\boldsymbol{v}d^{3}vd^{3}x\right)$, and energy

$\int\left[\sum_{s}\int\tfrac{1}{2}m_{s}v^{2}f_{s}d^{3}v+p_{e}/(\gamma-1)+B^{2}/2\mu_{0}\right]d^{3}x,$
(5)

see e.g. Refs. [4, 1] for details.

##  3 Discretized multi-rate scheme on curvilinear mapped meshes

###  3.1 Time discretization

For the time advance, we employ the second-order implicit-midpoint method. For
a variable $\chi$ with known value at timestep $t=n\Delta t$, i.e. $\chi^{n}$,
the new value $\chi^{n+1}$ is found implicitly from
$(\chi^{n+1}-\chi^{n})/\Delta t=f(\chi^{n+1/2})$ where
$\chi^{n+1/2}=\tfrac{1}{2}(\chi^{n+1}+\chi^{n})$. We note that this differs
from the Crank-Nicolson method in the discretization of non-linear terms. This
timestepping scheme is used here both for the grid-based equations (Sec. 3.3),
and for the particle advance (Sec. 3.4) where it gives long term accuracy for
ion orbit integration (it is symplectic). This method was found to give exact
discrete conservation in the Cartesian formulation of this hybrid-PIC kinetic-
ion and fluid-electron scheme [1].

###  3.2 Curvilinear map and basis

Assuming a static and non-singular curvilinear map from logical space
$\xi^{\alpha}=(\xi,\eta,\mu)$ to physical space $\boldsymbol{x}=(x,y,z)$, i.e.
$\boldsymbol{x}(\xi^{\alpha})$, the Jacobian of the mapping is defined as
$J=|\partial\boldsymbol{x}/\partial\xi^{\alpha}|$, the contravariant basis
vectors as $\boldsymbol{\nabla}\xi^{\alpha}$, and the covariant basis vectors
as
$\partial_{\alpha}\boldsymbol{x}\equiv\partial\boldsymbol{x}/\partial\xi^{\alpha}$.
Then, for a generic vector $\boldsymbol{S}$, the contravariant components are
defined as $S^{\alpha}=\boldsymbol{S}\cdot\boldsymbol{\nabla}\xi^{\alpha}$,
and the covariant components as
$S_{\alpha}=\boldsymbol{S}\cdot\partial_{\alpha}\boldsymbol{x}$. Indices can
be lowered as $S_{\alpha}=g_{\alpha\beta}S^{\beta}$ via the metric tensor
$g_{\alpha\beta}\equiv\partial_{\alpha}\boldsymbol{x}\cdot\partial_{\beta}\boldsymbol{x}$,
and raised as $S^{\alpha}=g^{\alpha\beta}S_{\beta}$ where
$g_{\alpha\beta}g^{\beta\gamma}=\delta_{\alpha}^{\gamma}$. A gives some
further definitions and vector identities used herein.

###  3.3 Discretization of the field equations

In hybrid-PIC discretizations of the model described in Section 2, Eqs. (2, 3,
4) for $\boldsymbol{B}$, $\boldsymbol{E}$, and $p_{e}$ are typically solved on
a spatial grid. Here, we use a uniform structured grid in logical space with
all grid variables $\chi_{g}\equiv\chi_{ijk}$ defined at the position of the
cell centers in logical space $(\xi_{i},\eta_{j},\mu_{k})$ for
$i\in[1,N_{\xi}]$, $j\in[1,N_{\eta}]$, $k\in[1,N_{\mu}]$. The subscript $g$
(shorthand for the 3-dimensional grid index $ijk$) is used to denote grid-
based variables herein. The spatial derivatives of these mesh variables are
computed using a second-order centered finite-difference approximation, as
e.g.
$(\partial_{\xi}\chi)_{g}\equiv(\partial_{\xi}\chi)_{ijk}=(\chi_{i+1jk}-\chi_{i-1jk})/(2\Delta\xi)$,
where $\Delta\xi$ is the uniform grid spacing on the logical mesh.

The covariant electric field $E_{\alpha,g}^{n+1/2}$ is found via static
evaluation of Ohm’s law at the $n+1/2$ time level as

$E_{\alpha,g}^{n+1/2}=-J_{g}\epsilon_{\alpha\beta\gamma}(u^{\beta})^{n+1/2}_{g}(B^{\gamma})^{n+1/2}_{g}+J_{g}\epsilon_{\alpha\beta\gamma}\frac{(j^{\beta})_{g}^{n+1/2}(B^{\gamma})^{n+1/2}_{g}}{en^{n+1/2}_{g}}-\frac{1}{en^{n+1/2}_{g}}\left(\frac{\partial
p_{e}}{\partial\xi^{\alpha}}\right)^{n+1/2}_{g}+\eta(j_{\alpha})_{g}^{n+1/2},$
(6)

where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol. The density
and velocity moments of the kinetic ions, $n_{g}$ and $(u^{\beta})_{g}$, are
defined below. $(B^{\gamma})_{g}$ are the contravariant components of the
magnetic field,
$(j^{\beta})_{g}=(1/J_{g})\epsilon^{\beta\sigma\alpha}(\partial_{\sigma}(g_{\alpha\gamma}B^{\gamma}))_{g}$
are the contravariant components of the current density, and $(p_{e})_{g}$ is
the electron scalar pressure.

The contravariant magnetic field components are calculated via Faraday’s
equation as

$\frac{(B^{\alpha})^{n+1}_{g}-(B^{\alpha})^{n}_{g}}{\Delta
t}=-\frac{1}{J_{g}}\epsilon^{\alpha\beta\gamma}(\partial_{\beta}E_{\gamma})^{n+1/2}_{g}.$
(7)

The scalar electron pressure is updated as

$\frac{1}{(\gamma-1)}\left[\frac{(Jp_{e})^{n+1}_{g}-(Jp_{e})^{n}_{g}}{\Delta
t}+\partial_{\alpha}\left(Ju_{e}^{\alpha}p_{e}\right)^{n+1/2}_{g}\right]+(p_{e})^{n+1/2}_{g}(\partial_{\alpha}Ju_{e}^{\alpha})^{n+1/2}_{g}=\eta
J_{g}(j_{\alpha}j^{\alpha})_{g}^{n+1/2}-\left(\partial_{\alpha}J_{g}q_{e}^{\alpha}\right)_{g}^{n+1/2},$
(8)

where $(u_{e}^{\alpha})_{g}=(u^{\alpha})_{g}-(j^{\alpha})_{g}/en_{g}$ is the
electron bulk velocity. Note that we have multiplied Eq. (8) by $J_{g}$, which
we assume to be constant in time. We note that these choices for
discretization of the spatial grid equations are motivated by those used for a
non-staggered resistive-MHD scheme described in Ref. [44].

###  3.4 Discretization of the particle governing equations

####  3.4.1 Choice of coordinates: logical in space, Cartesian in velocity-
space (hybrid push)

As discussed above, the real advantage of using curvilinear mapped meshes is
to either i) selectively resolve sharp boundary layers, or ii) to use a mesh
that aligns with certain features of the simulated problem, such as the
magnetic field direction or the shape of the domain boundaries. In contrast,
the ion kinetic species are discretized as a collection of Lagrangian phase-
space marker particles and are advanced in continuous space. Thus, their
accurate integration does not depend on the choice of coordinate system in the
same manner. There is freedom to formulate the particle equations of motion in
logical space, Cartesian space, or some hybrid method of the two. Here, we
formulate such a hybrid method in which the particle velocities are
represented in Cartesian coordinate directions but the particle positions are
updated in logical space. The advantage of the Cartesian velocity update is to
avoid the numerical difficulty of additional ficticious force terms in the
particle equations of motion, while the logical position update allows the
efficient tracking of the particle positions within curved cells (logical
cells are regular cubes).

####  3.4.2 Multi-rate time integration: sub-cycling of ion orbits

As in Ref. [1], multi-rate time integration is used where the particles can
have smaller timesteps than are used to solve Faraday’s equation (7) and the
electron pressure equation (8). The benefit of this sub-cycling is to give
accurate integration of particle orbits in regions of strong electromagnetic
fields, or within small cells, while avoiding the need to solve the full non-
linear system at each sub-step. This is particularly advantageous for many
problems of interest where the characteristic frequency
$\omega\ll\Omega_{ci}$, where the particle sub-step can be chosen to follow
the fast gyro-frequency $\Omega_{ci}$ and the macro-scale timestep can be
chosen to follow $\omega$. Each particle with index $p$ at substep $\nu$ is
advanced in time by $\Delta\tau_{p}^{\nu}$, which is calculated adaptively
using local error estimation [45, 1]. The number of substeps $N_{\nu,p}$ for
each particle is such that
$\sum_{\nu=0}^{N_{\nu,p}-1}\Delta\tau_{p}^{\nu}=\Delta t$.

####  3.4.3 Particle push: equations of motion

The logical-space position of each particle, $\xi^{\alpha}_{p}$, is updated in
time by the substep $\Delta\tau_{p}^{\nu}$ as

$\frac{\left(\xi^{\alpha}_{p}\right)^{\nu+1}-\left(\xi^{\alpha}_{p}\right)^{\nu}}{\Delta\tau_{p}^{\nu}}=\boldsymbol{v}_{p}^{\nu+1/2}\cdot\left(\boldsymbol{\nabla}\xi^{\alpha}_{p}\right)^{\nu+1/2},$
(9)

where $\boldsymbol{v}_{p}$ is the Cartesian particle velocity vector and the
covariant basis vector $\left(\boldsymbol{\nabla}\xi^{\alpha}_{p}\right)$ is
evaluated at the particle position. To calculate this basis function
numerically, we first construct a locally continuous map
$\boldsymbol{x}(\xi^{\alpha}_{p})=\sum_{g}\boldsymbol{x}_{g}S_{2}(\xi^{\alpha}_{g}-\xi^{\alpha}_{p})\equiv\sum_{ijk}\boldsymbol{x}_{ijk}S_{2}(\xi_{i}-\xi_{p})S_{2}(\eta_{j}-\eta_{p})S_{2}(\mu_{k}-\mu_{p})$,
using the tensor product of second-order B-spline shape functions. Then the
following identity is used, written for the $\xi$-component ($\alpha=1$) as

$\boldsymbol{\nabla}\xi_{p}=\frac{1}{J_{p}}\frac{\partial\boldsymbol{x}}{\partial\eta}\Bigg{|}_{p}\times\frac{\partial\boldsymbol{x}}{\partial\mu}\Bigg{|}_{p},$
(10)

where the partial derivatives are calculated using the analytical derivatives
of the shape functions, which themselves are combinations of lower-order
B-splines [12].

The Cartesian particle velocity vector $\boldsymbol{v}_{p}$ is updated as

$\frac{\boldsymbol{v}_{p}^{\nu+1}-\boldsymbol{v}_{p}^{\nu}}{\Delta\tau_{p}^{\nu}}=\frac{q_{p}}{m_{p}}\left(\boldsymbol{E}^{*,\nu+1/2}_{p}+\boldsymbol{v}_{p}^{\nu+1/2}\times\boldsymbol{B}_{p}^{\nu+1/2}\right),$
(11)

using the Cartesian electromagnetic field vectors at the particle position.

###  3.5 Particle-mesh interpolations

To close the set of equations, electromagnetic fields are scattered from the
spatial grid to the particle positions, and moments of the particle
distributions are gathered for use in Eqs. (6-8). For the momentum
conservation properties of our scheme described in Sections 4.2.1-4.2.2, we
find it necessary to scatter the Cartesian components of the electromagnetic
fields to the particle positions as

$\boldsymbol{E}_{p}^{*,\nu+1/2}=\sum_{g}\boldsymbol{E}^{*,n+1/2}_{g}S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2}),$
(12)
$\boldsymbol{B}_{p}^{\nu+1/2}=\sum_{g}\boldsymbol{B}^{n+1/2}_{g}S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2}),$
(13)

where $S(\xi^{\alpha}_{g}-\xi^{\alpha}_{p})\equiv
S_{m}(\xi_{i}-\xi_{p})S_{m}(\eta_{j}-\eta_{p})S_{m}(\mu_{k}-\mu_{p})$ is the
tensor product of $m$-th order B-spline shape functions. The electromagnetic
fields are assumed to be constant in time over the particle substeps (so are
given superscript $n+1/2$), such that the change in force occurs only due to
the change in particle substep position. The Cartesian vector values of the
fields are found on the spatial mesh as
$\boldsymbol{E}_{g}=E_{\alpha,g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}$,
$\boldsymbol{B}_{g}=(B^{\alpha})_{g}\left(\frac{\partial\boldsymbol{x}}{\partial\xi^{\alpha}}\right)_{g}$
prior to interpolation.

The moments are gathered using the orbit-averaging technique [46], where a
contribution to the moments at each sub-step is weighted by the fraction of
the macro-scale timestep, $\Delta\tau_{p}^{\nu}/\Delta t$. This orbit
averaging technique increases the number of statistical samples in the
calculation of the moments, thus reducing the noise. Also, as shown below, it
is required to give discrete conservation properties in the presence of
subcycling. The quasi-neutral density is gathered as

$n_{g}^{n+1/2}=\frac{1}{J_{g}\Delta
V}\sum_{p}\frac{q_{p}}{e}w_{p}\frac{1}{\Delta t}\sum_{\nu=0}^{N_{\nu
p}-1}S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2})\Delta\tau_{p}^{\nu},$
(14)

where $w_{p}$ is the macro-particle weight, which we assume to be constant in
time for each particle, and $\Delta V=\Delta\xi\Delta\eta\Delta\mu$ is the
constant logical cell volume. The Cartesian bulk momentum is gathered as

$(n\boldsymbol{u})_{g}^{n+1/2}=\frac{1}{J_{g}\Delta
V}\sum_{p}\frac{q_{p}}{e}w_{p}\frac{1}{\Delta t}\sum_{\nu=0}^{N_{\nu
p}-1}S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2})\boldsymbol{v}_{p}^{\nu+1/2}\Delta\tau_{p}^{\nu}.$
(15)

The contravariant component used in the field and moment equations is then
found as
$(nu^{\alpha})_{g}^{n+1/2}=(n\boldsymbol{u})_{g}\cdot(\boldsymbol{\nabla}\xi^{\alpha})_{g}$.

##  4 Conservation properties of the discrete model

###  4.1 Energy conservation

####  4.1.1 Kinetic energy of the ion particles

We define the discrete kinetic energy of the ion particles at time level
$(n+1)\Delta t$ as
$K_{i}^{n+1}=\sum_{p}w_{p}m_{p}(\boldsymbol{v}_{p}^{n+1})^{2}/2$. Noting that
$(n+1)\Delta t=n\Delta t+\sum_{\nu=0}^{N_{\nu,p}-1}\Delta\tau_{p}^{\nu}$, then
the change in kinetic energy over the macro timestep is

$\frac{K_{i}^{n+1}-K_{i}^{n}}{\Delta
t}=\sum_{p}w_{p}m_{p}\sum_{\nu=0}^{N_{\nu}-1}\frac{(\boldsymbol{v}_{p}^{\nu+1})^{2}-(\boldsymbol{v}_{p}^{\nu})^{2}}{2\Delta
t}=\sum_{p}w_{p}m_{p}\sum_{\nu=0}^{N_{\nu}-1}\boldsymbol{v}_{p}^{\nu+1/2}\cdot\left(\frac{\boldsymbol{v}_{p}^{\nu+1}-\boldsymbol{v}_{p}^{\nu}}{\Delta\tau_{p}^{\nu}}\right)\frac{\Delta\tau_{p}^{\nu}}{\Delta
t}.$ (16)

Then, using Eqs. (11,12, and 15),

$\frac{K_{i}^{n+1}-K_{i}^{n}}{\Delta
t}=\sum_{p}w_{p}q_{p}\sum_{\nu=0}^{N_{\nu}-1}\boldsymbol{v}_{p}^{\nu+1/2}\cdot\boldsymbol{E}_{p}^{*,\nu+1/2}\frac{\Delta\tau_{p}^{\nu}}{\Delta
t}=\sum_{g}e(Jn\boldsymbol{u})_{g}^{n+1/2}\cdot\boldsymbol{E}^{*,n+1/2}_{g}\Delta
V.$ (17)

Finally, re-writing this in terms of curvilinear components gives

$\frac{K_{i}^{n+1}-K_{i}^{n}}{\Delta
t}=\sum_{g}e(nu^{\alpha})_{g}^{n+1/2}(E_{\alpha}^{*})^{n+1/2}_{g}J_{g}\Delta
V.$

####  4.1.2 Magnetic energy on the discrete mesh

The total magnetic energy is defined as
$W_{B}^{n+1}=\sum_{g}(B^{\alpha})^{n+1}_{g}(B_{\alpha})^{n+1}_{g}J_{g}\Delta
V/(2\mu_{0})$. The rate of change over a macro-scale timestep can then be
written as

$\frac{W_{B}^{n+1}-W_{B}^{n}}{\Delta
t}=\sum_{g}\frac{(B^{\alpha})^{n+1}_{g}(B_{\alpha})^{n+1}_{g}-(B^{\alpha})_{g}^{n}(B_{\alpha})_{g}^{n}}{2\mu_{0}\Delta
t}J_{g}\Delta
V=\sum_{g}(g_{\alpha\beta})_{g}\left[\frac{(B^{\alpha})_{g}^{n+1}-(B^{\alpha})_{g}^{n}}{\mu_{0}\Delta
t}\right](B^{\beta})_{g}^{n+1/2}J_{g}\Delta V.$ (18)

Using Eq. (7),

$\frac{W_{B}^{n+1}-W_{B}^{n}}{\Delta
t}=\sum_{g}(g_{\alpha\beta})_{g}\left[-\epsilon^{\alpha\gamma\delta}(\partial_{\gamma}E_{\delta})_{g}^{n+1/2}\right](B^{\beta})_{g}^{n+1/2}\Delta
V.$ (19)

For suitable boundary conditions, such as periodic, it can be shown that
$\sum_{g}\chi_{g}(\partial_{\alpha}\phi_{g})=-\sum_{g}(\partial_{\alpha}\chi_{g})\phi_{g}$,
which is a discrete integration by parts. Using this relation gives

$\frac{W_{B}^{n+1}-W_{B}^{n}}{\Delta
t}=\sum_{g}\epsilon^{\alpha\gamma\delta}(\partial_{\gamma}g_{\alpha\beta}B^{\beta})_{g}^{n+1/2}(E_{\delta})_{g}^{n+1/2}\Delta
V.$ (20)

Finally, using the definition of the contravariant current density from Sec.
3.3 gives

$\frac{W_{B}^{n+1}-W_{B}^{n}}{\Delta
t}=\sum_{g}-(j^{\delta})_{g}^{n+1/2}(E_{\delta})_{g}^{n+1/2}J_{g}\Delta V.$
(21)

####  4.1.3 Conservation of total energy

The total electron kinetic energy (in the absence of electron inertia) is
given by $K_{e}^{n}=\sum_{g}(p_{e})_{g}^{n+1}J_{g}\Delta V/(\gamma-1)$.
Assuming suitable boundary conditions so that the flux terms can be neglected,
the total energy change

$\displaystyle\frac{(K_{i}+W_{B}+K_{e})^{n+1}-(K_{i}+W_{B}+K_{e})^{n}}{\Delta
t}=\sum_{g}$
$\displaystyle\left[e(nu^{\alpha})_{g}^{n+1/2}(E_{\alpha}^{*})^{n+1/2}_{g}-(j^{\alpha})_{g}^{n+1/2}(E_{\alpha})_{g}^{n+1/2}+\eta(j_{\alpha}j^{\alpha})_{g}^{n+1/2}\right]J_{g}\Delta
V$ (22) $\displaystyle-$
$\displaystyle\left[(p_{e})^{n+1/2}_{g}(\partial_{\alpha}Ju_{e}^{\alpha})^{n+1/2}_{g}\right]\Delta
V.$ (23)

Then, using Eq. (6) for the first and second terms on the right hand side and
using the discrete integration by parts for the final term, this gives

$\displaystyle\frac{(K_{i}+W_{B}+K_{e})^{n+1}-(K_{i}+W_{B}+K_{e})^{n}}{\Delta
t}$
$\displaystyle=\sum_{g}\left[-(u^{\alpha})_{g}^{n+1/2}+\frac{(j^{\alpha})_{g}^{n+1/2}}{en_{g}^{n+1/2}}+(u_{e}^{\alpha})_{g}^{n+1/2}\right]\left(\frac{\partial
p_{e}}{\partial\xi^{\alpha}}\right)_{g}^{n+1/2}J_{g}\Delta V$ (24)
$\displaystyle=0,$ (25)

using the definition of the electron bulk velocity.

###  4.2 Linear momentum

The discrete linear momentum is defined as
$\boldsymbol{P}^{n+1}=\sum_{p}w_{p}m_{p}\boldsymbol{v}_{p}^{n+1}$. The rate of
change over a macro-scale timestep can then be written as

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta
t}=\sum_{p}w_{p}m_{p}\sum_{\nu=0}^{N_{\nu}-1}\left(\frac{\boldsymbol{v}_{p}^{\nu+1}-\boldsymbol{v}_{p}^{\nu}}{\Delta\tau^{\nu}_{p}}\right)\frac{\Delta\tau^{\nu}_{p}}{\Delta
t}.$ (26)

Using Eqs. (11, 12, 13, 14 and 15), this can be written as

$\displaystyle\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta t}$
$\displaystyle=\sum_{p}w_{p}q_{p}\frac{1}{\Delta
t}\sum_{\nu=0}^{N_{\nu}-1}\left(\boldsymbol{E}_{p}^{\nu+1/2}+\boldsymbol{v}_{p}^{\nu+1/2}\times\boldsymbol{B}_{p}^{\nu+1/2}\right)\Delta\tau^{\nu}_{p}$
(27)
$\displaystyle=\sum_{g}e\left[n_{g}^{n+1/2}\boldsymbol{E}_{g}^{n+1/2}+(n\boldsymbol{u})_{g}^{n+1/2}\times\boldsymbol{B}_{g}^{n+1/2}\right]J_{g}\Delta
V,$

which can be written in curvilinear form as

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta
t}=\sum_{g}e\left[n_{g}^{n+1/2}(E_{\alpha})_{g}^{n+1/2}+J_{g}\epsilon_{\alpha\beta\gamma}(nu^{\beta})_{g}^{n+1/2}(B^{\gamma})_{g}^{n+1/2}\right]\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}J_{g}\Delta
V.$ (28)

Then, using the covariant form of Ohms law from Eq. (6), and cancelling the
convective electric field terms gives

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta
t}=\sum_{g}\left[J_{g}\epsilon_{\alpha\beta\gamma}\left(j^{\beta}\right)^{n+1/2}_{g}\left(B^{\gamma}\right)^{n+1/2}_{g}-\frac{\partial
p_{e}}{\partial\xi^{\alpha}}\right]\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}J_{g}\Delta
V.$ (29)

Finally, using the definition of the current density
$\left(j^{\beta}\right)_{g}^{n+1/2}$ and the contraction of indices for the
product of Levi-Cevita symbols
$\epsilon_{\alpha\beta\gamma}\epsilon^{\beta\sigma\kappa}=-(\delta_{\alpha}^{\sigma}\delta_{\gamma}^{\kappa}-\delta_{\alpha}^{\kappa}\delta_{\gamma}^{\sigma})$
gives

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta
t}=\sum_{g}\left[\left(\partial_{\gamma}B_{\alpha}\right)_{g}^{n+1/2}\left(B^{\gamma}\right)_{g}^{n+1/2}-\left(\partial_{\alpha}B_{\gamma}\right)_{g}^{n+1/2}\left(B^{\gamma}\right)_{g}^{n+1/2}-(\partial_{\alpha}p_{e})_{g}^{n+1/2}\right]\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}J_{g}\Delta
V.$ (30)

####  4.2.1 Linear momentum conservation for general curvilinear meshes in the
electrostatic limit

Firstly, we consider the electrostatic limit for which the first two terms
within the square brackets of Eq. (30) are zero. The total linear momentum is
then given by the term:

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta
t}=\sum_{g}-J_{g}(\boldsymbol{\nabla}\xi^{\alpha})_{g}\left(\frac{\partial
p_{e}}{\partial\xi^{\alpha}}\right)_{g}^{n+1/2}\Delta V.$ (31)

Using the discrete integration by parts with suitable (e.g., periodic)
boundary conditions, this can be written as

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta
t}=\sum_{g}(p_{e})_{g}^{n+1/2}\left(\partial_{\alpha}J_{g}(\boldsymbol{\nabla}\xi^{\alpha})_{g}\right)_{g}\Delta
V.$ (32)

As discussed in detail by Ref. [47], the continuum version of mesh quantity
term inside of the large bracket is identically equal to zero for any type of
curvilinear map. We have calculated this expression numerically, using our
definitions of the covariant derivative and basis functions above, and
confirmed that the discrete form of this expression

$\left(\partial_{\alpha}J_{g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}\right)_{g}=0$
(33)

holds locally at every point in our simulation domain down to numerical round-
off (at double precision). We have confirmed this result holds even for non-
orthogonal meshes. Thus, in the electrostatic limit, we have that

$\frac{\boldsymbol{P}^{n+1}-\boldsymbol{P}^{n}}{\Delta t}=0.$ (34)

Section 6.1 demonstrates a numerical example of an electrostatic problem that
confirms the conservation of total momentum for a non-orthogonal sinusoidal
mesh.

####  4.2.2 Linear momentum conservation for tensor-packed (Cartesian) meshes
in the full electromagnetic model

Next we consider the electromagnetic terms, which are the first two terms in
the square brackets of Eq. (30). We find that both of these terms may be non-
zero for general curvilinear meshes such as, for example, non-orthogonal
meshes or orthogonal meshes in cylindrical or spherical coordinate systems.
However, we will proceed to show that they are equal to zero for the subset of
Cartesian tensor-packed meshes that are defined by the invertable map
$\boldsymbol{x}(\xi^{\alpha})=(x(\xi),y(\eta),z(\mu))$. This is a incremental
generalization of the result from Ref. [1], where momentum conservation was
shown for this electromagnetic model on a uniform Cartesian mesh. However, it
is an important generalization as this form of mesh is frequently used to
study multi-scale problems for which localized regions of the domain require
finer resolution, such as the example discussed in Section 6.3. The
conservation of momentum for this type of mesh eliminates unphysical particle
self-forces that can otherwise lead to artificial particle reflection in
regions of strong mesh gradients, see Ref. [42].

We firstly note that in the case of $\gamma=\alpha$ the two terms trivially
cancel and thus we only need consider the components with $\gamma\neq\alpha$.
Writing the summation notation explicitly and dropping the temporal indices
(which are all equal) for simplicity of notation, the remaining components can
be written as

$\sum_{g}\sum_{\alpha}J_{g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}\Delta
V\sum_{\gamma\neq\alpha}\left[\left(\partial_{\gamma}B_{\alpha}\right)_{g}\left(B^{\gamma}\right)_{g}-\left(\partial_{\alpha}B_{\gamma}\right)_{g}\left(B^{\gamma}\right)_{g}\right].$
(35)

We consider these two terms separately, starting with the second term. We
multiply the local value of the expression in the square brackets by the
identity transformation
$\left(\frac{\partial\boldsymbol{x}}{\partial\xi^{\gamma}}\right)_{g}\cdot\left(\boldsymbol{\nabla}\xi^{\gamma}\right)\equiv
1$, as

$\sum_{\gamma\neq\alpha}\left(B^{\gamma}\right)_{g}\left(\partial_{\alpha}B_{\gamma}\right)_{g}\left(\frac{\partial\boldsymbol{x}}{\partial\xi^{\gamma}}\right)_{g}\cdot\boldsymbol{\nabla}\left(\xi^{\gamma}\right)_{g}.$
(36)

Then, using the map for the Cartesian tensor-packed mesh stated above, it can
be straightforwardly shown that the quantity
$(\boldsymbol{\nabla}\xi^{\gamma})_{g}$ is independent of $\alpha$ for
$\gamma\neq\alpha$ (this step is not true for a general curvilinear mesh), and
so it can be taken inside the convective derivative term to give

$\sum_{\gamma\neq\alpha}\left(B^{\gamma}\right)_{g}\left(\frac{\partial\boldsymbol{x}}{\partial\xi^{\gamma}}\right)_{g}\cdot\left(\partial_{\alpha}B_{\gamma}\boldsymbol{\nabla}\xi^{\gamma}\right)_{g}=\sum_{\gamma\neq\alpha}\boldsymbol{B}_{g}\cdot\left(\partial_{\alpha}\boldsymbol{B}\right)_{g}=\sum_{\gamma\neq\alpha}\left(\partial_{\alpha}\frac{B^{2}}{2}\right)_{g},$
(37)

where $\boldsymbol{B}$ is the Cartesian magnetic field vector. In the final
step, we have re-written the expression as the derivative of a particular flux
that satisfies the discrete chain rule locally and is conservative. It can be
shown that, for example, the definition of the ZIP flux satisfies this
property, see Refs. [48, 44]. Using this expression, the full second term in
Eq. (35) can be written as

$\sum_{g}\sum_{\alpha}J_{g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}\Delta
V\sum_{\gamma\neq\alpha}\left(\partial_{\alpha}\frac{B^{2}}{2}\right)_{g}.$
(38)

Then, if we consider the case of $\alpha=1$ without loss of generality, the
term $J_{g}\left(\boldsymbol{\nabla}\xi\right)_{g}$ can be written using the
definition of the basis functions as

$J_{g}\left(\boldsymbol{\nabla}\xi\right)_{g}=\frac{\partial\boldsymbol{x}}{\partial\eta}\Bigg{|}_{g}\times\frac{\partial\boldsymbol{x}}{\partial\mu}\Bigg{|}_{g}.$
(39)

For a Cartesian tensor-packed mesh, this is independent of the coordinate
$\xi$ and thus the term can be taken inside the convective derivative term to
give

$\sum_{g}\sum_{\alpha}\Delta
V\sum_{\gamma\neq\alpha}\left(\partial_{\alpha}J_{g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}\frac{B^{2}}{2}\right)_{g}=0.$
(40)

For the first term in Eq. (35), we assume a suitably defined (ZIP) flux and
rewrite this using the discrete local chain rule as

$\sum_{g}\sum_{\alpha}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}\Delta
V\sum_{\gamma\neq\alpha}\left[\left(\partial_{\gamma}J_{g}B_{\alpha}B^{\gamma}\right)_{g}-B_{\alpha}\underbrace{\left(\partial_{\gamma}J_{g}B^{\gamma}\right)}_{=0}\right],$
(41)

where the second term in this expression is zero due to the solenoidal
property of the magnetic field, which holds locally for our discretization for
general curvilinear meshes - see Ref. [44]. Then, as before, noting that
$(\boldsymbol{\nabla}\xi^{\alpha})_{g}$ is independent of $\gamma$ for
$\gamma\neq\alpha$, this can be written as a pure derivative of fluxes which
is conservative, i.e.

$\sum_{g}\sum_{\alpha}\Delta
V\sum_{\gamma\neq\alpha}\left[\left(\partial_{\gamma}J_{g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}B_{\alpha}B^{\gamma}\right)_{g}\right]=0.$
(42)

The numerical confirmation of momentum conservation in Cartesian tensor-packed
meshes is demonstrated in the numerical examples of Sections 6.2 and 6.3.

##  5 Numerical implementation

###  5.1 Iterative solution

As in Ref. [1], our implementation solves for the magnetic vector potential
$(A_{\alpha})^{n+1}_{g}$ rather than the magnetic field formulation
$(B^{\alpha})_{g}^{n+1}$ used in Eq. (7). The Weyl gauge is chosen with
electrostatic potential $\phi=0$, such that $(A_{\alpha})^{n+1}_{g}$ is found
from the solution of

$\frac{(A_{\alpha})_{g}^{n+1}-(A_{\alpha})_{g}^{n}}{\Delta
t}=-\left(E_{\alpha}\right)_{g}^{n+1/2}.$ (43)

We note that all of the conservation properties of Sec. 4 also hold for the
$\boldsymbol{A}$-formulation with the definition of
$(B^{\alpha})_{g}=\epsilon^{\alpha\beta\gamma}\left(\partial_{\beta}A_{\gamma}\right)_{g}/J_{g}$.

The non-linear algebraic equations (43, 8) can be written in the form
$\boldsymbol{G}(\boldsymbol{y}^{n+1})=\boldsymbol{0}$ for the solution vector
$\boldsymbol{y}^{n+1}=\left((A_{\alpha})_{g}^{n+1},(p_{e})_{g}^{n+1}\right)$.
The full details of the method of solution closely follow Ref. [1]. The system
is approximately inverted by using a Jacobian-free Newton Krylov (JFNK) method
[49] using flexible-GMRES [50] as a linear solver. The outer Newton iteration
is iterated until the residual is converged to a user specified non-linear
tolerance $\epsilon_{t}$. Note that the particle quantities
$(\xi^{\alpha}_{p},\boldsymbol{v}_{p})$ do not enter the residual directly;
rather, the particles are advanced by Eqs. (9, 11) and moments gathered by
Eqs. (14, 15) in each residual evaluation. To do this, a Picard-iterated
implicit Boris method is used to solve Eqs. (9, 11) - see Ref. [1].

Effective preconditioning plays a primary role in the efficiency of methods
such as these. However, the discussion of a suitable physics-based
preconditioner is outside of the scope of the present paper, and will be
documented elsewhere. For these results, we do not perform preconditioning
and, as such, we typically use a small timestep and do not consider
performance data for the results presented here.

###  5.2 Normalization

The discrete model is implemented in normalized form using natural magnetized
ion units. For a typical magnetic field strength $B_{0}$, density $n_{0}$, and
ion mass $m_{0}$, velocities are normalized by the Alfvén speed
$v_{0}=v_{A0}=B_{0}/\sqrt{m_{0}\mu_{0}n_{0}}$, times by the inverse cyclotron
period $t_{0}=\Omega_{c0}^{-1}=m_{0}/eB_{0}$, and lengths by the ion skin
depth $L_{0}=d_{i0}=v_{A0}/\Omega_{c0}$. The different species of kinetic ions
are distinguished by their normalized charge $Z_{s}=q_{s}/e$, mass
$M_{s}=m_{s}/m_{0}$, and density $N_{s}=n_{s0}/n_{0}$. If they are initialized
with a Maxwellian distribution function with temperature $T_{s0}$, we define
the thermal velocity as $v_{Ts0}=\sqrt{T_{s0}/m_{s}}$ and plasma-beta
$\beta_{s0}=2M_{s}N_{s}(v_{Ts}/v_{A0})^{2}=2\mu_{0}n_{s0}T_{s0}/B_{0}^{2}$.
The electron sound speed is defined as $C_{s}=\sqrt{T_{e0}/m_{0}}$ and the
ratio of ion to electron temperatures is $\tau_{s}=T_{s0}/T_{e0}$.

##  6 Numerical Verification

###  6.1 Electrostatic Ion Acoustic Wave: Momentum conservation

In the first numerical example, we consider the ion Landau damping of an ion
acoustic wave using a 2D non-orthogonal curvilinear mesh. This is an
electrostatic problem and the motivation is to numerically verify the novel
conservation properties of the scheme for this limit (momentum and energy
conservation, see Sec. 4.2.1) for non-orthogonal mapped meshes. The physical
problem set-up considered is similar to that described in Ref. [1]. The
simulation is initialized with a single ion species with $Z_{i}=M_{i}=1$ with
a Maxwellian velocity function with thermal speed $v_{Ti0}=\sqrt{1/3}$. We set
$\gamma=5/3$ and use a temperature ratio of $\tau=0.2$. The ions are given a
perturbation to their velocity such that the initial bulk velocity is given by
$\boldsymbol{\delta
u}_{i}=0.01\cos{\left(k_{x}x+k_{y}y\right)}\left(\boldsymbol{\hat{x}}+\boldsymbol{\hat{y}}\right)$
with $k_{x}=k_{y}=\pi/8$. We use a non-orthogonal sinusoidal mesh with
periodic boundary conditions for this problem, specified by the map

$x(\xi,\eta)=\xi+\sigma\sin{(2\pi\xi/L_{x})}\sin{(2\pi\eta/L_{y})},$ (44)
$y(\xi,\eta)=\eta+\sigma\sin{(2\pi\xi/L_{x})}\sin{(2\pi\eta/L_{y})}.$ (45)

Here, the parameter $\sigma$ describes the distortion of the mesh. Fig. 1
shows the initial conditions for the $\boldsymbol{\hat{x}}$-component of the
ion velocity moment, along with the spatial mesh used for $\sigma=1$.

Additional numerical parameters used are $64\times 64$ cells, with $5\times
10^{4}$ particles/cell initialized with a quasi-quiet start [51] from a low
discrepancy Hammersley sequence [52] to reduce noise. We use second-order
quadratic-spline shape functions and apply two passes of conservative binomial
smoothing, see C. The timestep used is $\Delta t=0.02$, and the non-linear
tolerance used is $\epsilon_{t}=10^{-12}$.

Figure 1: Non-orthogonal sinusoidal mesh used for electrostatic ion acoustic
wave problem. The colored lines show contours of the Cartesian ion velocity
moment $u_{ix}$, and the black $+$ symbols show the centers of the spatial
cells. The mesh map is given by Eqs. (44-45) with distortion factor
$\sigma=1$.

Figure 2 (left) shows the amplitude of the ion velocity moment vs time for the
sinusoidal mesh (red) of Eqs. (44, 45), as well as a uniform Cartesian mesh
with the same resolution (blue). The amplitude of the wave decreases
exponentially with time due to kinetic ion Landau damping. The theoretical
damping rate is found from the dispersion relation $Z^{\prime}(\xi)=2\tau$,
where $Z$ is the plasma dispersion function and
$\xi=(\omega+i\gamma)/k\sqrt{2}v_{Ti0}$ is the normalized complex frequency.
For our initial conditions, the theoretical value is $\gamma=-0.0744636$,
which is shown as the black dashed line on Fig. 2 (left). There is good
agreement for both the Cartesian and curvilinear mesh cases with this
theoretical damping rate until the noise due to a finite number of particles
overcomes the signal.

The conservation errors are shown in Fig. 2 (right) for both the Cartesian
mesh (blue) and curvilinear mesh (red) simulations. The momenta and energy are
conserved to better than $10^{-13}$ and $10^{-12}$, respectively, which is on
the level of the non-linear tolerance ($10^{-12}$) used. This provides
numerical verification of the result of Section 4.2.1, that, for the
electrostatic limit, the total linear momentum and energy is conserved for
general curvilinear mapped meshes.

(a)  (b)

Figure 2:  a) Ion velocity moment perturbation amplitude from the
electrostatic ion acoustic wave test on non-orthogonal sinusoidal (red) and
uniform Cartesian (blue) meshes. The black dashed line shows the theoretical
damping rate with slope $\gamma=-0.0744636$. b) Numerical conservation errors.
Top: $\boldsymbol{\hat{x}}$-component of linear momentum. Middle:
$\boldsymbol{\hat{y}}$-component of linear momentum. Bottom: Total (ion +
electron kinetic) energy. These errors are comparable to the non-linear
tolerance of $10^{-12}$, as expected.

###  6.2 Spatial convergence for 2D whistler wave propagation

In the second numerical example, we simulate the right-hand polarized Alfvén-
whistler wave in a cold and uniform background plasma. This wave has an exact
analytical solution to the cold-plasma hybrid model for arbitrary perturbation
amplitude, and therefore can be used to perform a formal grid convergence
study. All simulations are performed in 2D, with the initial background plasma
specified by $Z_{i}=M_{i}=N_{i}=1$, $T_{i}=T_{e}=0$, and the magnetic field
oriented diagonally in physical space as
$\boldsymbol{B}/B_{0}=(\boldsymbol{\hat{x}}+\boldsymbol{\hat{y}})/\sqrt{2}$.

The eigenmode of the wave is perturbed as

$\boldsymbol{\delta B}=\boldsymbol{\nabla}\times\boldsymbol{\delta
A}=10^{-3}\boldsymbol{\nabla}\times\left(\frac{\cos{(\boldsymbol{k}\cdot\boldsymbol{x})}}{k_{x}}\boldsymbol{\hat{y}}-\frac{\sin{(\boldsymbol{k}\cdot\boldsymbol{x})}}{k}\boldsymbol{\hat{z}}\right),$
(46) $\frac{\boldsymbol{\delta u}}{v_{A0}}=-\frac{\boldsymbol{\delta
B}}{B_{0}}\frac{v_{A0}k}{\omega},$ (47)

where $k_{x}d_{i}=k_{y}d_{i}=2\pi/16$ and
$kd_{i}=\sqrt{k_{x}^{2}+k_{y}^{2}}\approx 0.56$, and
$\omega=v_{A0}k\left(\tfrac{1}{2}d_{i}k+\sqrt{1+\tfrac{1}{4}d_{i}^{2}k^{2}}\right)$.
This intermediate value of $kd_{i}$ is chosen to simultaneously test both the
grid-based electron and particle-based ion discretizations. These simulations
use $400$ particles/cell initialized with the quasi-quiet start, and apply one
pass of conservative binomial smoothing (C). The timestep used is $\Delta
t\Omega_{ci}=0.0125$.

Figure 3: Meshes used for Alfvén-whistler wave propagation problem. The
background magnetic field and direction of propagation are along the diagonal
($45^{\circ}$) of the 2D domain. The colored lines show contours of the
magnetic vector potential component $A_{z}$, and the black $+$ symbols show
the centers of the spatial cells. Left: Cartesian tensor-packed mesh with
packing factor of 4 in the $x$ and $y$ directions. Right: Non-orthogonal
sinusoidal mesh of Eqs. (44-45) with distortion factor $\sigma=0.4$.

Figure 3 shows the two types of mesh that are used for this numerical test.
The left panel shows a Cartesian tensor-packed mesh of the kind discussed in
Section 4.2.2 that is packed by a factor of $4\times$ in the $x$ and $y$
directions, and the right panel shows a non-orthogonal sinusoidal mesh of Eqs.
(44, 45) with distortion factor $\sigma=0.4$.

Figure 4: Spatial grid convergence test for whistler-wave 2D propagation
problem. The uniform-mesh case is given by black squares. Circles are from the
sinusoidal grid runs with $\sigma=0.5$ (red), $\sigma=1$ (magenta), and
$\sigma=1.5$ (green). Triangles are from the stretched planar mesh runs with
stretch factors of $2\times$ (blue), and $4\times$ (cyan). The black and green
dashed lines show ideal $(\Delta x)^{2}$ convergence.

Figure 4 shows a formal spatial convergence study of the numerical scheme for
the Cartesian tensor-packed mesh as triangles (packing factor of $2\times$ in
blue, and $4\times$ in cyan), the sinusoidal mesh as circles ($\sigma=0.5$ in
red, $\sigma=1$ in magenta, and $\sigma=1.5$ in green), as well as a uniform
Cartesian mesh (black squares) for comparison. The black and green dashed
lines are shown for reference, both with $\epsilon\propto(\Delta x)^{2}$. This
order of convergence holds for all of the different spatial meshes.

(a) Change in components of linear momentum.

(b) Change in magnetic energy, ion kinetic energy, and total energy.

Figure 5:  Conservation errors in momentum and total energy for Alfvén-
whistler wave problem. The blue curve shows results for a Cartesian tensor-
packed mesh with packing factor of $2\times$, and the magenta curve is from
the non-orthogonal sinusoidal mesh simulation with distortion factor
$\sigma=1$.

Figure 5 (left) shows the three components of the linear momentum error,
$\boldsymbol{\Delta P}=\boldsymbol{P}^{n}-\boldsymbol{P}^{0}$, for the
Cartesian tensor-packed mesh with a packing factor of $2\times$ (blue) and the
sinusoidal mesh with $\sigma=1$ (magenta). Both simulations have
$n_{\xi}=n_{\eta}=128$, with other parameters the same as for Fig. 4. In
agreement with Section 4.2.2, the Cartesian tensor-packed meshes show very
small momentum errors $(\textrm{max}\left[|\boldsymbol{\Delta
P}|\right]\approx 10^{-13}$) for a non-linear tolerance of $10^{-12}$, whereas
the sinusoidal mesh has a larger momentum error,
$\textrm{max}\left[|\boldsymbol{\Delta P}|\right]\approx 6\times 10^{-10}$.

Figure 5 (right) shows the energies from the same two simulations. The top two
panels show the variation in the magnetic energy ($\Delta
W_{B}=W_{B}^{n}-W_{B}^{0}$) and ion kinetic energy ($\Delta
K_{i}=K_{i}^{n}-K_{i}^{0}$), respectively, and the bottom panel shows the sum
of the two. The total energy errors for both cases are comparable, with values
$\textrm{max}\left[\Delta\epsilon_{T}\right]\approx 10^{-13}$. We have
verified that this error scales with the chosen non-linear tolerance of the
iterative solver, as expected.

###  6.3 Non-linear magnetic island coalescence problem

The third numerical example is the electromagnetic island coalescence problem,
which is a non-linear and multi-scale magnetic reconnection challenge problem
[53, 54, 55, 14, 56, 57, 58, 59]. Recent simulation studies have highlighted
the key role that the ion kinetic physics has in determining the rate of
magnetic reconnection, and the global motions of the interacting magnetic
islands [14].

We assume periodic boundary conditions at $x=0,L_{x}$, and $y=0,L_{y}$, such
that there are four magnetic islands in total within the domain - two islands
cross the periodic boundary at $x=0$ and $x=L_{x}$ as shown in Figure 6. The
initial conditions given below are in an unstable ideal-magnetohydrodynamic
equilibrium state, where magnetic islands with the same sign of the out-of-
plane current density $j_{z}$ attract each-other due to the Lorentz force.

A single ion species is used with $Z_{i}=M_{i}=1$. The ions are initialized
with a Maxwellian distribution function with initial density profile [60]

$N_{i}=0.2+\frac{1-\alpha^{2}}{(\cosh{(x/\lambda)}+\alpha\cos{(4\pi
y/L_{y})})^{2}}+\frac{1-\alpha^{2}}{(\cosh{((x-L_{x}/2)/\lambda)}+\alpha\cos{(4\pi
y/L_{y})})^{2}}+\frac{1-\alpha^{2}}{(\cosh{((x-L_{x})/\lambda)}+\alpha\cos{(4\pi
y/L_{y})})^{2}},$ (48)

where $\alpha=0.4$, $\lambda=2d_{i}$, and $L_{x}=L_{y}=25.6d_{i}$. The
magnetic vector potential is given by

$\displaystyle A_{z}=$
$\displaystyle-\lambda\\{\ln{[\cosh{(x/\lambda)}+\alpha\cos{(4\pi
y/L_{y})}]}+\ln{[\cosh{((x-L_{x})/\lambda)}+\alpha\cos{(4\pi y/L_{y})}]}$ (49)
$\displaystyle-\ln{[\cosh{((x-0.5L_{x})/\lambda)}+\alpha\cos{(4\pi
y/L_{y})}]}\\}.$

The ion plasma beta is $\beta_{i}=5/6$ with an ion-to-electron temperature
ratio $\tau=5$ (total ion $+$ electron $\beta=1$). The plasma resistivity used
is $\eta=5\times 10^{-4}$, and the heat conductivity is $\kappa=0$. Here, we
also include a higher-order hyper-resistive contribution to the resistive
friction term such that
$F_{ie}=\eta\boldsymbol{j}-\eta_{H}\nabla^{2}\boldsymbol{j}$ with
$\eta_{H}=2\times 10^{-4}$. This is used to set a dissipation scale for the
thin current sheets that form in Fig. 6, otherwise these can shrink to the
mesh-scale [61]. We note that the hyper-resistivity is not needed for
numerical stability in this case, and the runs shown in Fig. 7 do not include
the hyper-resistivity. To break the initial symmetry, we apply a sinusoidal
perturbation with amplitude $\delta A_{z}=0.1$, which causes the negative
current islands (red) to move towards the midplane ($y=L_{y}/2$), and the
positive current islands (blue) to move towards the outer boundaries at $y=0$
and $y=L_{y}$. Magnetic reconnection occurs at each of the magnetic
$X$-points, located at $(x,y)=(0,0)$, $(0,L_{y}/2)$, $(L_{x}/2,0)$,
$(L_{x}/2,L_{y}/2)$, and the matching points on periodic boundaries $x=L_{x}$
and $y=L_{y}$.

(a)

(b)

(c)

Figure 6: Pseudocolor plot of out-of-plane current density $j_{z}$ from the
island coalescence problem with $4$ initial islands ($2$ are shown at the
periodic boundary in $x$). Left column: Cartesian tensor-packed mesh (shown as
black lines) with $n_{\xi}=n_{\eta}=128$ and stretch factor of $4\times$ in
$x$ and $y$ directions. Fine resolution is placed along the simulation
boundaries and in the center of the domain. Right column: Uniform Cartesian
mesh with $n_{x}=n_{y}=512$. Both simulations use $n_{ppc}=400$, so the
tensor-packed mesh case has $16\times$ fewer particles.

Figure 6 shows a comparison between a simulation with a Cartesian tensor-
packed mesh with resolution $n_{\xi}=n_{\eta}=128$ (left column) and a
$n_{x}=n_{y}=512$ resolution uniform Cartesian mesh (right column). Both cases
use $400$ particles/cell which are initialized with a quasi-quiet start, and
use quadratic-spline shape functions with two passes of conservative smoothing
(C). The timestep used is $\Delta t\Omega_{ci}=0.0125$. The map used in the
tensor-packed mesh is chosen to give the same resolution at the magnetic
X-points as for the uniform mesh case. However, as the same number of average
particles per cell is used, the tensor-packed mesh case uses $16\times$ fewer
particles with corresponding reductions in both CPU time and memory usage.
Despite this, the rate of coalescence of the magnetic islands is very similar
between the two simulations as shown by the similar size of the islands at
$t\Omega_{ci}=40$ in the bottom row.

(a) Change in components of linear momentum. (b) Change in magnetic energy,
electron kinetic energy, ion kinetic energy, and total energy.

Figure 7:  Conservation errors in momentum and total energy for the tensor-
packed mesh island coalescence simulation at three values of the relative non-
linear tolerance (exact conservation is expected for arbitrary small
tolerance).

Figure 7 shows the momenta and energies from the same tensor-packed mesh
simulation of Fig. 6 (left) but with zero hyper-resistivity ($\eta_{H}=0$) for
three different values of the relative non-linear tolerance: $10^{-4}$
(green), $10^{-6}$ (red), and $10^{-8}$ (blue). For these tensor-packed
Cartesian meshes, both the momentum and total energy errors are reduced for
tighter values of the non-linear tolerance, as expected.

###  6.4 Kink ($m=1$) mode in helical geometry

In the final numerical example, we illustrate the algorithm for a more
challenging geometry with non-periodic boundary conditions. The $m=n=1$ mode
of a cylindrical plasma column is simulated on a 2D non-orthogonal helically
symmetric mesh with an underlying map given by $x=r\cos{\tau}$,
$y=r\sin{\tau}$. Here, $\tau=(\theta-kz)/m$ is the helical angle, $\theta$ the
poloidal angle, $z$ the axial position, $k=n/R$ the axial wavenumber, and $R$
is the major radius of the cylinder which is a topological torus due to
periodic boundary conditions in the poloidal and axial directions. In
addition, we convolve this map with a packing function that is non-uniform in
the radial direction such that resolution is finer at the rational surface
($\Delta r_{\textrm{min}}=0.002$ at $r=0.5$) and decreased in the cells
adjacent to the singular point ($r\approx 0$). The latter is helpful to
alleviate high levels of statistical noise that otherwise can occur in these
cells due to their extremely small cell volumes and thus low numbers of
particles per cell. The mesh used for the simulation is shown in Fig. 8. The
spatial mesh has $[n_{r},n_{\tau}]=[64,64]$ spatial points with an average of
$400$ macro-particles per cell, but these are distributed non-uniformly to
place more particles into regions of low cell volumes that are more
susceptible to noise as the simulation progresses. The boundary conditions are
a singular-point boundary at $r=0$ and a conducting outer boundary at $r=1$,
which are described in B. We use a tensor product of zero and second-order
shape functions (B) and apply two passes of conservative smoothing (C). The
timestep used is $\Delta t=0.01$ and the non-linear tolerance is
$\epsilon_{t}=10^{-3}$.

We use a single ion species with $Z_{i}=M_{i}=1$. The initial profiles are 1D
cylindrically symmetric profiles, which form an exact Vlasov equilibrium [62].
Macroscopically, the columns magnetic pinch force is balanced by a thermal
pressure gradient that points radially outwards. The plasma density and
poloidal magnetic field profiles are given by

$N_{i}=\frac{1}{\left(1+(r/\lambda)^{2}\right)^{2}},$ (50)
$B_{\theta}/B_{0}=\frac{\left(r/\lambda\right)}{\left(1+(r/\lambda)^{2}\right)}$
(51)

with $B_{z}/B_{0}=1$, such that the helical field
$B_{\tau}/B_{0}=m(B_{\theta}/B_{0}r)+k$. Initial equilibrium requires that
$\beta_{T}=2n_{0}(T_{i0}+T_{e0})/B_{0}^{2}=1$. The physical parameters used
are an ion-to-electron temperature ratio $\tau=T_{i}/T_{e}=0$ (chosen to
compare with a cold-ion Hall-MHD model, see below), resistivity $\eta=2\times
10^{-3}$ and heat conductivity $\kappa=2\times 10^{-3}$. We also note that the
normalization for this simulation differs from Sec. 5.2, as we set $L_{0}$ as
the minor radius with $d_{i}/L_{0}=0.1$. Similarly, the time-scale is
normalized by the Alfven time $\tau_{0}=L_{0}/v_{A0}$.

We apply a small perturbation to the radial ion bulk velocity moment

$\delta u_{r}=10^{-3}\sin^{3}{(\pi r)}\cos{(\theta)}$ (52)

to displace the column from equilibrium.

In order to sample the initial density profile from Eq. (50) and velocity
moment perturbation from Eq. (52) to a high degree of accuracy for the non-
uniformly distributed macro-particle profile, we compute the individual
particle weights ($w_{p}$) and velocities ($\boldsymbol{v}_{p}^{0}$) by
inverting a mass-matrix [63]. These particles are initialized with a quasi-
quiet start.

Figure 8: Magnetic vector potential $A_{z}$ at initial ($t=0$) and final time
($t=9.5\tau_{0}$) from the hybrid-PIC simulation of the $m=1$ mode in 2D
helical geometry. Also shown is the 2D helical mesh which is packed in the
radial direction at the location of the initial rational surface.

Figure 8 shows the component of the magnetic vector potential $A_{z}$ (color
scale) along with the spatial mesh (black lines) for the initial and final
time of the simulation. At the final time, the plasma column shows a
characteristic helical kink-like displacement from the origin. This $m=1$ mode
is driven by the gradients in the current density and plasma pressure.

(a) Radial profile of $u_{\theta}$ in linear growth phase (b) Mode amplitudes
of density (top) and magnetic vector potential (bottom) perturbations

Figure 9: Eigenfunction of $m=1$ mode in helical geometry simulation. Green:
Simulation result from non-linear Hall-MHD simulation. Blue: Simulation from
cold ion hybrid simulation. The linear theory result (red) is computed for
resistive-MHD (RMHD), which is in reasonable agreement since the inertial
layer width is larger than the ion kinetic scales for these parameters.

Figure 9 (left) shows a comparison of the radial profile of the poloidal bulk
velocity $u_{\theta}$ from the simulation, along with the linear eigenfunction
of the $m=1$ mode. This eigenfunction and the theoretical linear growth rate
(red lines in the right panel of Fig. 9) are computed numerically by solving
the linearized resistive-MHD equations in cylindrical geometry using the
equilbrium specified above. Although the resistive-MHD theory misses some
important diamagnetic and finite gyro-radius physics – which gives a real
frequency of the mode and rotation observed in Fig. 8 (right) – our choice of
parameters ensures that these differences remain small. In particular, we pick
initial conditions which have a large MHD drive such that the inertial layer
thickness $\Delta_{in}$ is larger than the ion kinetic scales [64], which
allows comparison with resistive-MHD theory. We also verify the results
against a non-linear Hall-MHD simulation with the same initial set-up (shown
in green). In the limit of cold ions, the hybrid-PIC and Hall-MHD models are
physically equivalent but there are significant algorithmic differences,
namely that Hall-MHD model uses a Eulerian momentum equation to compute the
bulk velocity on the spatial mesh rather than using a Lagrangian particle
discretization. As shown in Fig. 9, there is good agreement between these
algorithms. Small differences close to the singular point ($r<0.1$) are due to
statistical noise in the hybrid-PIC simulation, which is larger in this region
due to the small cell volumes close to the singular point. The numerical
profiles are also in agreement with the resistive-MHD theoretical
eigenfunction.

The amplitudes of the density $||\delta n||(t)=||n(t)-n(t=0)||$ and vector
potential perturbations $||\delta A_{z}||(t)=||A_{z}(t)-A_{z}(t=0)||$ are
shown in the right panel of Fig. 9 from the Hall-MHD (green) and hybrid-PIC
(blue) simulations. The growth rate is in good agreement between the two
codes, and these are both in agreement with the resistive-MHD linear theory
(red). We note that the hybrid-PIC simulation is only run up to $t\approx
9.5\tau_{0}$. At this time in the non-linear regime, there are significant
bulk flows across the singular point region causing the number of particles
within these cells to decrease. The simulation then fails to converge due to
the increased noise fluctuation level. We note that achieving the current
simulation results was non-trivial, requiring the use of suitably tailored
radial mesh packing profiles, macro-particle placement, and adjustment of
particle weights to match the initial moment profiles. Future directions will
explore the use of particle remapping and/or $\delta F$ methods in order to
drastically reduce noise and allow the simulation to proceed to saturation.
These are significant modifications to the algorithm and are outside of the
scope of the present paper.

##  7 Summary

In this paper, we have derived a novel conservative scheme for the
electromagnetic hybrid kinetic-ion fluid-electron plasma model by using
implicit particle-in-cell methods. This work extends the work of Ref. [1] –
which was derived for uniform Cartesian meshes – to general curvilinear meshes
by using a smooth and static map from logical to physical space. The equations
for the electromagnetic fields and electron fluid quantities are discretized
on a cell-centered mesh in logical space. For the kinetic-ion particles, the
velocities are calculated in Cartesian space to avoid ficticious force terms
in the equation of motion, while the particle positions are computed in
logical space. The gather and scatter particle-mesh interpolations are
performed using the Cartesian vector components at the logical space
positions. These choices, along with the multi-rate implicit mid-point
timestepping scheme used in Ref. [1], give the unique conservation properties
of this algorithm. Sec. 4.1 demonstrates that the scheme discretely conserves
the total (ion kinetic $+$ magnetic $+$ electron thermal) energy for general
curvilinear meshes, provided suitable boundary conditions are used.
Furthermore, these choices also give linear momentum conservation in the
electrostatic limit as proved in Section 4.2.1, and demonstrated numerically
in Section 6.1. For the full electromagnetic model, total momentum is
conserved for tensor-packed meshes using a Cartesian coordinate basis (Sec.
4.2.2 and verified numerically in Secs. 6.2, 6.3). Although this property only
applies to this restricted subset, these types of meshes are often used to
efficiently study multiscale problems with thin boundary layers such as the
magnetic island coalescence problem of Sec. 6.3. In Sec. 6, we also present a
formal grid convergence study of the algorithm for simulating right-hand
polarized whistler waves, which demonstrates second-order convergence for
different types of meshes, and also present results from a challenging
simulation of the $m=n=1$ mode of a plasma column on a 2D helically symmetric
mesh geometry with non-periodic boundary conditions.

Future work on the kinetic-ion fluid-electron model will include the
development of optimal and robust physics-based preconditioning strategies,
and document the performance of the scheme. We will also explore noise
reduction via $\delta F$ control-variate techniques.

## Acknowledgements

This work is supported by the Applied Scientific Computing Research (ASCR)
program of the U.S. Department of Energy, and used resources provided by the
Los Alamos National Laboratory Institutional Computing Program, which is
supported by the U.S. Department of Energy National Nuclear Security
Administration under Contract No. DE-AC52-06NA25396.

##  Appendix A Discrete vector identities in curvilinear coordinates

The vector product operations used in Sec. 3 are defined in curvilinear form
using index notation as

$\boldsymbol{a}\cdot\boldsymbol{b}=a^{\alpha}b_{\alpha}=a_{\alpha}b^{\alpha},$
(53)
$\left(\boldsymbol{a}\times\boldsymbol{b}\right)^{\alpha}=\frac{1}{J}\epsilon^{\alpha\beta\gamma}a_{\beta}b_{\gamma},$
(54)
$\left(\boldsymbol{a}\times\boldsymbol{b}\right)_{\alpha}=J\epsilon_{\alpha\beta\gamma}a^{\beta}b^{\gamma},$
(55)

where $\epsilon^{\alpha\beta\gamma}$, $\epsilon_{\alpha\beta\gamma}$ are Levi-
Cevita symbols. The differential operators are defined as

$\boldsymbol{\nabla}\chi=\partial_{\alpha}\chi,$ (56)
$\boldsymbol{\nabla}\cdot\boldsymbol{a}=\frac{1}{J}\partial_{\alpha}\left(Ja^{\alpha}\right),$
(57)
$\left(\boldsymbol{\nabla}\times\boldsymbol{a}\right)^{\alpha}=\frac{1}{J}\epsilon^{\alpha\beta\gamma}\partial_{\beta}a_{\gamma}.$
(58)

##  Appendix B Boundary conditions

The helical kink mode example problem discussed in Sec. 6.4 uses non-periodic
boundaries in the radial direction. These are a singular point boundary
condition at $r=0$, and a conducting shell boundary condition at $r=1$. To
simplify the treatment of these boundary conditions, we use Nearest Grid Point
(NGP) interpolation for the radial direction in the scatter and gather
operations of Eqs. (13, 12, 14, 15), while using second order in the periodic
direction of the helical angle. This choice was made to achieve some noise
reduction while ensuring that the particle-mesh interpolations only access
cells interior to the domain boundaries.

###  B.1 Conducting boundary condition

For the particles we use a reflecting boundary condition

$\boldsymbol{v}_{p}\rightarrow\boldsymbol{v}_{p}-2\boldsymbol{v}_{p}\cdot\boldsymbol{\hat{n}}\boldsymbol{\hat{n}},$
(59)

where $\boldsymbol{\hat{n}}$ is the normal vector to the boundary evaluated at
the particle position. This conserves the particle energy, but not the
momentum.

For the conducting wall we set
$\boldsymbol{\hat{n}}\times\boldsymbol{E}=-\partial_{t}\boldsymbol{\hat{n}}\times\boldsymbol{A}=0$
and $\partial_{t}p_{e}=0$. In cylindrical and helical coordinates this gives

$p_{e}(t)=p_{e}(t=0),\quad A_{2}(t)=A_{2}(t=0),\quad A_{3}(t)=A_{3}(t=0).$
(60)

###  B.2 Singular point boundary condition

For mesh quantities, we note that our cell centered description means that no
quantities are defined exactly on the singular point. For a radial mesh
packing function that is symmetric about the singular point, the ghost cell
(index $i=0$) is located at $r_{0}=-r_{1}$ where $r_{1}$ is the cell center of
the first interior cell (index $i=1$). To set the singular point boundary
conditions requires values to be set at this ghost cell. We follow Ref. [65]
and make use of the identity transformation for cylindrical/helical
coordinates about the singular point $(r,\theta)\rightarrow(-r,\theta+\pi)$.
Scalar quantities are thus filled as

$\chi(r_{0},\theta)=\chi(r_{1},\theta+\pi).$ (61)

For our definition of basis functions in Sec. 3, the contravariant components
of a vector $\boldsymbol{\Pi}$ can be expressed in cylindrical co-ordinates as
$\Pi^{1}=\Pi_{r}$, $\Pi^{2}=\Pi_{\theta}/r$ and $\Pi^{3}=\Pi_{z}$. The ghost
cell values are found as follows
$\Pi^{1}(r_{0},\theta)=-\Pi^{1}(r_{1},\theta+\pi)$,
$\Pi^{2}(r_{0},\theta)=\Pi^{2}(r_{1},\theta+\pi)$,
$\Pi^{3}(r_{0},\theta)=\Pi^{3}(r_{1},\theta+\pi)$. The covariant components
are given by $\Pi_{1}=\Pi_{r}$, $\Pi_{2}=r\Pi_{\theta}$, $\Pi_{3}=\Pi_{z}$,
with the ghost cells being set in the same way. We note that this method is
the same as in Ref. [65] after accounting for differences in the basis vector
definitions used in Sec. 3.

##  Appendix C Conservative binomial smoothing

Binomial smoothing is a useful technique for full-kinetic-PIC or hybrid-PIC
schemes to reduce the effects of statistical noise due to a finite number of
macro-particles [12], while being computationally cheaper to apply than using
higher-order shape functions (at least for the typical case where the number
of particles-per-cell is large). The 3D smoothing operator is defined as
$\textrm{SM}_{g}[Q_{g}]\equiv\textrm{SM}_{k}[\textrm{SM}_{j}[\textrm{SM}_{i}[Q_{ijk}]]]$,
which is a composition of 1D smoothing operators defined as

$\textrm{SM}_{i}\left[Q_{ijk}\right]=\frac{Q_{i-1jk}+2Q_{ijk}+Q_{i+1jk}}{4}.$
(62)

It is possible to do this smoothing in a conservative manner for the hybrid-
PIC method, provided that it is done symmetrically to the electromagnetic
field and moment quantities as defined in e.g., Ref. [1]. This follows
directly from the property of Eq. (62), that
$\sum_{i}A_{ijk}\textrm{SM}_{i}\left[B_{ijk}\right]=\sum_{i}\textrm{SM}_{i}\left[A_{ijk}\right]B_{ijk}$
for suitable (e.g., periodic) boundary conditions. In a similar manner, it can
be shown that the discrete conservation properties of Section 4 continue to
hold for curvilinear meshes when the smoothing is done as described below.

For the electromagnetic fields, the Cartesian vector representations of Eqs.
(12, 13) are smoothed on the spatial mesh before scattering to the particles
as

$\boldsymbol{E}_{p}^{*,\nu+1/2}=\sum_{g}\textrm{SM}_{g}\left[\boldsymbol{E}^{*,n+1/2}_{g}\right]S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2}),$
(63)
$\boldsymbol{B}_{p}^{\nu+1/2}=\sum_{g}\textrm{SM}_{g}\left[\boldsymbol{B}^{n+1/2}_{g}\right]S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2}).$
(64)

In addition, the Jacobian-weighted density and Cartesian momentum vectors of
Eqs. (14, 15) are smoothed as

$J_{g}n_{g}^{n+1/2}=\textrm{SM}_{g}\left[\frac{1}{\Delta
V}\sum_{p}\frac{q_{p}}{e}w_{p}\frac{1}{\Delta t}\sum_{\nu=0}^{N_{\nu
p}-1}S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2})\Delta\tau_{p}^{\nu}\right]_{g},$
(65) $J_{g}(n\boldsymbol{u})_{g}^{n+1/2}=\textrm{SM}_{g}\left[\frac{1}{\Delta
V}\sum_{p}\frac{q_{p}}{e}w_{p}\frac{1}{\Delta t}\sum_{\nu=0}^{N_{\nu
p}-1}S(\xi^{\alpha}_{g}-(\xi^{\alpha})_{p}^{\nu+1/2})\boldsymbol{v}_{p}^{\nu+1/2}\Delta\tau_{p}^{\nu}\right]_{g}.$
(66)

We note that, as discussed in Ref. [43], care must be taken to resolve the ion
skin depth ($d_{i}$) when using higher-order particle shape functions or mesh-
based smoothing, to avoid numerical dispersion errors due to the hybrid-PIC
cancellation problem. Also, the compensation filter can be additionally
applied in the same manner as Eqs. (63-66) to reduce the amount of attenuation
caused by smoothing at long wavelengths ($k\Delta x\ll\pi$) where it is
undesirable, see e.g., Refs [12, 1].

## References

  * [1] A. Stanier, L. Chacon, G. Chen, A fully implicit, conservative, non-linear, electromagnetic hybrid particle-ion/fluid-electron algorithm, Journal of Computational Physics 376 (2019) 597–616.
  * [2] E. Priest, Magnetohydrodynamics of the Sun, Cambridge University Press, 2014.
  * [3] K. Schindler, Physics of space plasma activity, Cambridge University Press, 2006\.
  * [4] A. S. Lipatov, The hybrid multiscale simulation technology: an introduction with application to astrophysical and laboratory plasmas, Scientific Computation, Springer-Verlag Berlin Heidelberg, 2002.
  * [5] J. Byers, B. Cohen, W. Condit, J. Hanson, Hybrid simulations of quasineutral phenomena in magnetized plasma, Journal of Computational Physics 27 (3) (1978) 363 – 396.  doi:https://doi.org/10.1016/0021-9991(78)90016-5.
URL http://www.sciencedirect.com/science/article/pii/0021999178900165

  * [6] D. W. Hewett, C. W. Nielson, A Multidimensional Quasineutral Plasma Simulation Model, Journal of Computational Physics 29 (2) (1978) 219–236.  doi:10.1016/0021-9991(78)90153-5.
  * [7] D. Winske, L. Yin, N. Omidi, H. Karimabadi, K. Quest, Hybrid Simulation Codes: Past, Present and Future - A Tutorial, in: J. Büchner, C. Dum, M. Scholer (Eds.), Space Plasma Simulation, Vol. 615 of Lecture Notes in Physics, Berlin Springer Verlag, 2003, pp. 136–165.
  * [8] E. Belova, R. Denton, A. Chan, Hybrid simulations of the effects of energetic particles on low-frequency mhd waves, Journal of Computational Physics 136 (2) (1997) 324–336.  doi:https://doi.org/10.1006/jcph.1997.5738.
URL https://www.sciencedirect.com/science/article/pii/S0021999197957387

  * [9] J. Drake, H. Arnold, M. Swisdak, J. Dahlin, A computational model for exploring particle acceleration during reconnection in macroscale systems, Physics of Plasmas 26 (1) (2019) 012901.
  * [10] F. Holderied, S. Possanner, X. Wang, Mhd-kinetic hybrid code based on structure-preserving finite elements with particles-in-cell, Journal of Computational Physics 433 (2021) 110143.  doi:https://doi.org/10.1016/j.jcp.2021.110143.
URL https://www.sciencedirect.com/science/article/pii/S0021999121000358

  * [11] T. Amano, A generalized quasi-neutral fluid-particle hybrid plasma model and its application to energetic-particle-magnetohydrodynamics hybrid simulation, Journal of Computational Physics 366 (2018) 366–385.
  * [12] C. K. Birdsall, A. B. Langdon, Plasma Physics via Computer Simulation, McGraw-Hill, New York, 1991.
  * [13] H. Karimabadi, D. Krauss-Varban, J. D. Huba, H. X. Vu, On magnetic reconnection regimes and associated three-dimensional asymmetries: Hybrid, Hall-less hybrid, and Hall-MHD simulations, Journal of Geophysical Research (Space Physics) 109 (2004) A09205.  doi:10.1029/2004JA010478.
  * [14] A. Stanier, W. Daughton, L. Chacón, H. Karimabadi, J. Ng, Y.-M. Huang, A. Hakim, A. Bhattacharjee, Role of Ion Kinetic Physics in the Interaction of Magnetic Flux Ropes, Physical Review Letters 115 (17) (2015) 175004\.  doi:10.1103/PhysRevLett.115.175004.
  * [15] A. Le, A. Stanier, W. Daughton, J. Ng, J. Egedal, W. D. Nystrom, R. Bird, Three-dimensional stability of current sheets supported by electron pressure anisotropy, Physics of Plasmas 26 (10) (2019) 102114.
  * [16] D. Winske, Hybrid simulation codes with application to shocks and upstream waves, Space Science Reviews 42 (1-2) (1985) 53–66.
  * [17] M. S. Weidl, D. Winske, F. Jenko, C. Niemann, Hybrid simulations of a parallel collisionless shock in the large plasma device, Physics of Plasmas 23 (12) (2016) 122102.
  * [18] A. Le, D. Winske, A. Stanier, W. Daughton, M. Cowee, B. Wetherton, F. Guo, Astrophysical explosions revisited: collisionless coupling of debris to magnetized plasma, Journal of Geophysical Research: Space Physics (2021) e2021JA029125.
  * [19] L. Franci, S. Landi, L. Matteini, A. Verdini, P. Hellinger, High-resolution hybrid simulations of kinetic plasma turbulence at proton scales, The Astrophysical Journal 812 (1) (2015) 21.
  * [20] S. S. Cerri, F. Califano, F. Jenko, D. Told, F. Rincon, Subproton-scale cascades in solar wind turbulence: driven hybrid-kinetic simulations, The Astrophysical Journal Letters 822 (1) (2016) L12.
  * [21] A. Le, V. Roytershteyn, H. Karimabadi, A. Stanier, L. Chacon, K. Schneider, Wavelet methods for studying the onset of strong plasma turbulence, Physics of Plasmas 25 (12) (2018) 122310.
  * [22] J. Müller, S. Simon, Y.-C. Wang, U. Motschmann, D. Heyner, J. Schüle, W.-H. Ip, G. Kleindienst, G. J. Pringle, Origin of mercury’s double magnetopause: 3d hybrid simulation study with aikef, Icarus 218 (1) (2012) 666–687.
  * [23] H. Karimabadi, V. Roytershteyn, H. X. Vu, Y. A. Omelchenko, J. Scudder, W. Daughton, A. Dimmock, K. Nykyri, M. Wan, D. Sibeck, M. Tatineni, A. Majumdar, B. Loring, B. Geveci, The link between shocks, turbulence, and magnetic reconnection in collisionless plasmas, Physics of Plasmas 21 (6) (2014) 062308.  doi:10.1063/1.4882875.
  * [24] N. Omidi, G. Collinson, D. Sibeck, Structure and properties of the foreshock at venus, Journal of Geophysical Research: Space Physics 122 (10) (2017) 10–275.
  * [25] J. Ng, L.-J. Chen, Y. Omelchenko, Bursty magnetic reconnection at the earth’s magnetopause triggered by high-speed jets, Physics of Plasmas 28 (9) (2021) 092902\.
  * [26] D. S. Harned, Quasineutral hybrid simulation of macroscopic plasma phenomena, Journal of Computational Physics 47 (1982) 452–462.  doi:10.1016/0021-9991(82)90094-8.
  * [27] D. Winske, K. B. Quest, Electromagnetic ion beam instabilities - Comparison of one- and two-dimensional simulations, Journal of Geophysical Research 91 (1986) 8789–8797.  doi:10.1029/JA091iA08p08789.
  * [28] M. W. Kunz, G. Lesur, Magnetic self-organization in Hall-dominated magnetorotational turbulence, Monthly Notices of the Royal Astronomical Society 434 (2013) 2295–2312.  arXiv:1306.5887, doi:10.1093/mnras/stt1171.
  * [29] D. Winske, K. B. Quest, Magnetic field and density fluctuations at perpendicular supercritical collisionless shocks, Journal of Geophysical Research 93 (1988) 9681–9693.  doi:10.1029/JA093iA09p09681.
  * [30] A. P. Matthews, Current advance method and cyclic leapfrog for 2d multispecies hybrid plasma simulations, Journal of Computational Physics 112 (1) (1994) 102 – 116.  doi:http://dx.doi.org/10.1006/jcph.1994.1084.
URL //www.sciencedirect.com/science/article/pii/S0021999184710849

  * [31] M. Fujimoto, Instabilities in the magnetopause velocity shear layer, Ph.D. thesis, Univ. of Tokyo, Tokyo (1990).
  * [32] V. A. Thomas, D. Winske, N. Omidi, Re-forming supercritical quasi-parallel shocks: 1. one- and two-dimensional simulations, Journal of Geophysical Research: Space Physics 95 (A11) (1990) 18809–18819.  doi:10.1029/JA095iA11p18809.
URL http://dx.doi.org/10.1029/JA095iA11p18809

  * [33] D. W. Swift, Use of a hybrid code to model the earth’s magnetosphere, Geophysical Research Letters 22 (3) (1995) 311–314.  doi:10.1029/94GL03082.
URL http://dx.doi.org/10.1029/94GL03082

  * [34] D. W. Hewett, A global method of solving the electron-field equations in a zero-inertia-electron-hybrid plasma simulation code, Journal of Computational Physics 38 (1980) 378–395.  doi:10.1016/0021-9991(80)90155-2.
  * [35] T. Amano, K. Higashimori, K. Shirakawa, A robust method for handling low density regions in hybrid simulations for collisionless plasmas, Journal of Computational Physics 275 (2014) 197–212.  arXiv:1406.6613, doi:10.1016/j.jcp.2014.06.048.
  * [36] B. J. Sturdevant, Y. Chen, S. E. Parker, Low frequency fully kinetic simulation of the toroidal ion temperature gradient instability, Physics of Plasmas 24 (8) (2017) 081207.  doi:10.1063/1.4999945.
  * [37] P. W. Rambo, Finite-Grid Instability in Quasineutral Hybrid Simulations, Journal of Computational Physics 118 (1995) 152–158.  doi:10.1006/jcph.1995.1086.
  * [38] D. Higginson, P. Amendt, N. Meezan, W. Riedel, H. Rinderknecht, S. Wilks, G. Zimmerman, Hybrid particle-in-cell simulations of laser-driven plasma interpenetration, heating, and entrainment, Physics of Plasmas 26 (11) (2019) 112107\.
  * [39] D. Leblond, Simulation des plasmas de tokamak avec xtor: régimes des dents de scie et évolution vers une modélisation cinétique des ions, Ph.D. thesis (2011).
  * [40] S. Dyadechkin, E. Kallio, R. Jarvinen, A new 3-d spherical hybrid model for solar wind interaction studies, Journal of Geophysical Research: Space Physics 118 (8) (2013) 5157–5168.
  * [41] J. Müller, S. Simon, U. Motschmann, J. Schüle, K.-H. Glassmeier, G. J. Pringle, Aikef: Adaptive hybrid model for space plasma simulations, Computer Physics Communications 182 (4) (2011) 946–966.
  * [42] P. Colella, P. C. Norgaard, Controlling self-force errors at refinement boundaries for amr-pic, Journal of Computational Physics 229 (4) (2010) 947–957.
  * [43] A. Stanier, L. Chacon, A. Le, A cancellation problem in hybrid particle-in-cell schemes due to finite particle size, Journal of Computational Physics 420 (2020) 109705.
  * [44] L. Chacón, A non-staggered, conservative, div(B)=0, finite-volume scheme for 3D implicit extended magnetohydrodynamics in curvilinear geometries, Computer Physics Communications 163 (2004) 143–171.  doi:10.1016/j.cpc.2004.08.005.
  * [45] G. Chen, L. Chacón, A multi-dimensional, energy- and charge-conserving, nonlinearly implicit, electromagnetic Vlasov-Darwin particle-in-cell algorithm, Computer Physics Communications 197 (2015) 73–87.  arXiv:1503.01169, doi:10.1016/j.cpc.2015.08.008.
  * [46] B. I. Cohen, R. P. Freis, V. Thomas, Orbit-averaged implicit particle codes, Journal of Computational Physics 45 (3) (1982) 345 – 366.  doi:https://doi.org/10.1016/0021-9991(82)90108-5.
URL http://www.sciencedirect.com/science/article/pii/0021999182901085

  * [47] V. D. Liseikin, Grid generation methods, Springer, 2017.
  * [48] C. W. Hirt, Heuristic Stability Theory for Finite-Difference Equations, Journal of Computational Physics 2 (1968) 339–355.  doi:10.1016/0021-9991(68)90041-7.
  * [49] D. A. Knoll, D. E. Keyes, Jacobian-free Newton-Krylov methods: a survey of approaches and applications, Journal of Computational Physics 193 (2004) 357–397.  doi:10.1016/j.jcp.2003.08.010.
  * [50] Y. Saad, A flexible inner-outer preconditioned gmres algorithm, SIAM Journal on Scientific Computing 14 (2) (1993) 461–469.  arXiv:https://doi.org/10.1137/0914028, doi:10.1137/0914028.
URL https://doi.org/10.1137/0914028

  * [51] R. Sydora, Low-noise electromagnetic and relativistic particle-in-cell plasma simulation models, Journal of computational and applied mathematics 109 (1-2) (1999) 243–259.
  * [52] J. Hammersley, Monte carlo methods, Springer Science & Business Media, 2013.
  * [53] J. M. Finn, P. Kaw, Coalescence instability of magnetic islands, The Physics of Fluids 20 (1) (1977) 72–78.
  * [54] P. Pritchett, The coalescence instability in collisionless plasmas, Physics of Fluids B: Plasma Physics 4 (10) (1992) 3371–3381.
  * [55] H. Karimabadi, J. Dorelli, V. Roytershteyn, W. Daughton, L. Chacón, Phys. Rev. Lett. 107 (2) (2011) 025002.  doi:10.1103/PhysRevLett.107.025002.
  * [56] J. Ng, Y.-M. Huang, A. Hakim, A. Bhattacharjee, A. Stanier, W. Daughton, L. Wang, K. Germaschewski, The island coalescence problem: Scaling of reconnection in extended fluid models including higher-order moments, Physics of Plasmas 22 (11) (2015) 112104.  arXiv:1511.00741, doi:10.1063/1.4935302.
  * [57] A. Stanier, W. Daughton, A. N. Simakov, L. Chacon, A. Le, H. Karimabadi, J. Ng, A. Bhattacharjee, The role of guide field in magnetic reconnection driven by island coalescence, Physics of Plasmas 24 (2) (2017) 022124.
  * [58] K. Makwana, R. Keppens, G. Lapenta, Two-way coupled mhd-pic simulations of magnetic reconnection in magnetic island coalescence, in: Journal of Physics: Conference Series, Vol. 1031, IOP Publishing, 2018, p. 012019.
  * [59] F. Allmann-Rahn, T. Trost, R. Grauer, Temperature gradient driven heat flux closure in fluid simulations of collisionless reconnection, Journal of Plasma Physics 84 (3).
  * [60] V. Fadeev, I. Kvabtskhava, N. Komarov, Self-focusing of local plasma currents, Nuclear fusion 5 (3) (1965) 202.
  * [61] A. Stanier, A. N. Simakov, L. Chacón, W. Daughton, Fluid vs. kinetic magnetic reconnection with strong guide fields, Physics of Plasmas 22 (10) (2015) 101203.  doi:10.1063/1.4932330.
  * [62] W. H. Bennett, Magnetically self-focussing streams, Phys. Rev. 45 (1934) 890–897.  doi:10.1103/PhysRev.45.890.
URL https://link.aps.org/doi/10.1103/PhysRev.45.890

  * [63] G. Chen, L. Chacón, T. B. Nguyen, An unsupervised machine-learning checkpoint-restart algorithm using gaussian mixtures for particle-in-cell simulations, Journal of Computational Physics 436 (2021) 110185.
  * [64] A. Mishchenko, A. Zocco, Global gyrokinetic particle-in-cell simulations of internal kink instabilities, Physics of Plasmas 19 (12) (2012) 122104.
  * [65] G. L. Delzanno, L. Chacón, J. M. Finn, Electrostatic mode associated with the pinch velocity in reversed field pinch simulations, Physics of Plasmas 15 (12) (2008) 122102.


$\displaystyle\geq J_{0}(\pi^{\star})-\mathcal{O}(\epsilon)$

where the second inequality follows by the feasibility of $\pi^{\star}$ for
$\mathcal{P}(\bm{\tau})$; the last inequality follows by
$\eta\|\bar{\bm{w}}\|_{1}\leq\eta D_{w}=\mathcal{O}(\epsilon)$. This proves
near optimality of $\bar{\pi}$. Now we prove that $\bar{\pi}$ is (exactly)
feasible for $\mathcal{P}(\bm{\tau})$.

#### Exact Feasibility

Define $m=\min_{i\in[I]}(J_{i}(\bar{\pi})-\tau_{i})$. If $m\geq 0$ then
$J_{i}(\bar{\pi})-\tau_{i}\geq 0$ for all $i=1,\dots,I$ and exact feasibility
trivially holds. We only consider the case where $m<0$. Define a mixture
policy $\widetilde{\pi}=(1-\zeta)\pi^{\star}+\zeta\pi_{\phi}^{\star}$ where
$\zeta\in(0,1)$ is to be determined later. The mixture policy has the
interpretation of first drawing a policy from
$\\{\pi^{\star},\pi_{\phi}^{\star}\\}$ with probabilities $(1-\zeta)$ and
$\zeta$, then running the drawn policy for the entire trajectory. Since
$L_{\eta}(\cdot,\bar{\bm{w}})$ is linear, a linear combination of (15) and
(16) with coefficients $1-\zeta$ and $\zeta$ respectively, we get

$L_{\eta}(\widetilde{\pi},\bar{\bm{w}})\leq
L_{\eta}(\bar{\pi},\bm{w})+\epsilon.$

Choosing $\bm{w}$ such that $w_{j}=D_{w}$ for
$j=\operatorname*{argmin}_{i\in[I]}(J_{i}(\bar{\pi})-\tau_{i})$ and $w_{j}=0$
for all other indices, we get

$\displaystyle L_{\eta}(\widetilde{\pi},\bar{\bm{w}})$ $\displaystyle\leq
J_{0}(\bar{\pi})+\bm{w}\cdot(\bm{J}(\bar{\pi})-\bm{\tau}-\eta\bm{1})+\epsilon$
$\displaystyle=J_{0}(\bar{\pi})+D_{w}(m-\eta)+\epsilon.$

On the other hand, using the fact that $\widetilde{\pi}$ is feasible for
$\mathcal{P}(\bm{\tau}+\zeta\phi\bm{1})$, we get

$\displaystyle L_{\eta}(\widetilde{\pi},\bar{\bm{w}})$
$\displaystyle=J_{0}(\widetilde{\pi})+\bar{\bm{w}}\cdot(\bm{J}(\widetilde{\pi})-\bm{\tau}-\eta\bm{1})$
$\displaystyle\geq
J_{0}(\widetilde{\pi})+(\zeta\phi-\eta)\|\bar{\bm{w}}\|_{1}.$

Combining the previous two results (upper bound and lower bound of
$L_{\eta}(\widetilde{\pi},\bar{\bm{w}})$) and rearranging, we get

$J_{0}(\widetilde{\pi})-J_{0}(\bar{\pi})\leq
D_{w}(m-\eta)-(\zeta\phi-\eta)\|\bar{\bm{w}}\|_{1}+\epsilon.$ (17)

Now, to get a lower bound of $J_{0}(\widetilde{\pi})-J_{0}(\bar{\pi})$, let
$(\widetilde{\pi}^{\star},\widetilde{\bm{w}}^{\star})$ be a primal-dual
solution of $\mathcal{P}(\bm{\tau}+\zeta\phi\bm{1})$. Note that
$\mathcal{P}(\bm{\tau}+\zeta\phi\bm{1})$ is feasible by the Slater’s condition
assumption A and the fact that $\zeta\phi\in(0,\phi)$. Since
$(\widetilde{\pi}^{\star},\widetilde{\bm{w}}^{\star})$ is a saddle point of
$L_{\zeta\phi}(\pi,\bm{w})=J_{0}(\pi)+\bm{w}\cdot(J_{\bm{C}}(\pi)-\bm{\tau}-\zeta\phi\bm{1})$,
we get

$L_{\zeta\phi}(\bar{\pi},\widetilde{\bm{w}}^{\star})\leq
L_{\zeta\phi}(\widetilde{\pi}^{\star},\widetilde{\bm{w}}^{\star})=J_{0}(\widetilde{\pi}^{\star})\leq
J_{0}(\pi^{\star})\leq\frac{1}{1-\zeta}J_{0}(\widetilde{\pi})$

where the equality follows by the complementary slackness property; the second
inequality follows since the feasibility set of $\mathcal{P}(\bm{\tau})$
contains that of $\mathcal{P}(\bm{\tau}+\zeta\phi\bm{1})$; and the last
inequality follows by
$J_{0}(\widetilde{\pi})=(1-\zeta)J_{0}(\pi^{\star})+\zeta
J_{0}(\pi_{\phi}^{\star})\geq(1-\zeta)J_{0}(\pi^{\star})$. Rearranging, we get

$\displaystyle J_{0}(\widetilde{\pi})-J_{0}(\bar{\pi})$
$\displaystyle\geq-\zeta
J_{0}(\bar{\pi})+(1-\zeta)\widetilde{\bm{w}}^{\star}\cdot(\bm{J}(\bar{\pi})-\bm{\tau}-\zeta\phi\bm{1})$
$\displaystyle\geq-\zeta+(1-\zeta)(m-\zeta\phi)\|\widetilde{\bm{w}}^{\star}\|_{1}$

where the second inequality follows by $J_{0}(\cdot)\leq 1$ and the definition
of $m$. Combining with the upper bound of
$J_{0}(\widetilde{\pi})-J_{0}(\bar{\pi})$ shown in (17) and rearranging, we
get

$(D_{w}-(1-\zeta)\|\widetilde{\bm{w}}^{\star}\|_{1})m\geq
D_{w}\eta+(\zeta\phi-\eta)\|\bar{\bm{w}}\|_{1}-\zeta-(1-\zeta)\zeta\phi\|\widetilde{\bm{w}}^{\star}\|_{1}-\epsilon.$
(18)

Now, we choose our parameters as follows.

$\zeta=\epsilon,\quad D_{w}=\frac{4}{\phi},\quad\eta=\phi\epsilon.$

Since $\widetilde{\bm{w}}^{\star}$ is a dual solution of
$\mathcal{P}(\bm{\tau}+\zeta\phi\bm{1})$, which has a margin of
$\phi-\zeta\phi$, Lemma 9 gives
$\|\widetilde{\bm{w}}^{\star}\|_{1}\leq\frac{1}{\phi-\zeta\phi}$. Hence,

$\zeta\phi\|\widetilde{\bm{w}}^{\star}\|_{1}\leq\frac{\zeta\phi}{\phi-\zeta\phi}\leq
2\zeta=2\epsilon$

where the second inequality uses $\zeta=\epsilon\leq\frac{1}{2}$. Hence,
$\|\widetilde{\bm{w}}^{\star}\|_{1}\leq\frac{2\epsilon}{\zeta\phi}=\frac{2}{\phi}<D_{w}$
so that $D_{w}-(1-\zeta)\|\widetilde{\bm{w}}^{\star}\|_{1}>0$. Hence, the
previous result (18) gives

$\displaystyle(D_{w}-(1-\zeta)\|\widetilde{\bm{w}}^{\star}\|_{1})m$
$\displaystyle\geq
D_{w}\eta+(\zeta\phi-\eta)\|\bar{\bm{\lambda}}\|_{1}-\zeta-(1-\zeta)\zeta\phi\|\widetilde{\bm{w}}^{\star}\|_{1}-\epsilon$
$\displaystyle\geq 4\epsilon+0-\epsilon-2\epsilon-\epsilon$ $\displaystyle=0.$

Since $D_{w}-(1-\zeta)\|\widetilde{\bm{\lambda}}^{\star}\|_{1}>0$, we have
$m\geq 0$ which implies $\tau_{i}-J_{i}(\bar{\pi})\geq 0$ for all
$i=1,\dots,I$. This leads to the following theorem.


Many scientific problems require identifying a small set of covariates that are associated with a target response and estimating their effects. Often, these effects are nonlinear and include interactions, so linear and additive methods can lead to poor estimation and variable selection. Unfortunately, methods that simultaneously express sparsity, nonlinearity, and interactions are computationally intractable — with runtime at least quadratic in the number of covariates, and often worse. In the present work, we solve this computational bottleneck. We show that suitable interaction models have a kernel representation, namely there exists a "kernel trick" to perform variable selection and estimation in $O$(# covariates) time. Our resulting fit corresponds to a sparse orthogonal decomposition of the regression function in a Hilbert space (i.e., a functional ANOVA decomposition), where interaction effects represent all variation that cannot be explained by lower-order effects. On a variety of synthetic and real datasets, our approach outperforms existing methods used for large, high-dimensional datasets while remaining competitive (or being orders of magnitude faster) in runtime.

§ INTRODUCTION

§ A FRAMEWORK FOR NON-LINEAR INTERACTIONS AND SPARSITY

§ USING TWO KERNEL TRICKS TO REDUCE COMPUTATION COST

§ HOW TO GET SPARSITY, INTERACTIONS, AND FAST INFERENCE WHEN COVARIATES ARE DEPENDENT

§ FINAL ALGORITHM AND IMPLEMENTATION DETAILS

§ RELATED WORK

§ EXPERIMENTS

[1000 Genomes Project, 2015]
1000 Genomes Project, C. (2015) A global reference for human genetic variation.
Nature, 526, 68–74.

[Agrawal et al., 2019]
Agrawal, R., Trippe, B., Huggins, J. and Broderick, T. (2019) The kernel
interaction trick: Fast Bayesian discovery of pairwise interactions in high
In ICML.

[Aschard, 2016]
Aschard, H. (2016) A perspective on interaction effects in genetic association
Genetic Epidemiology.

[Bien et al., 2013]
Bien, J., Taylor, J. and Tibshirani, R. (2013) A Lasso for hierarchical
The Ann. Stat., 41, 1111–1141.

[Candes and Tao, 2007]
Candes, E. and Tao, T. (2007) The Dantzig selector: Statistical estimation
when p is much larger than n.
The Ann. Stat., 2313–2351.

[Chen et al., 1998]
Chen, S., Donoho, D. and Saunders, M. (1998) Atomic decomposition by basis
SIAM Journal on Scientific Computing, 33–61.

[Chipman, 1996]
Chipman, H. (1996) Bayesian variable selection with related predictors.
The Canadian Journal of Statistics, 24, 17–36.

[Curtis et al., 2014]
Curtis, S. M., Banerjee, S. and Ghosal, S. (2014) Fast Bayesian model
assessment for nonparametric additive regression.

[Durrande et al., 2013]
Durrande, N., Ginsbourger, D., Roustant, O. and Carraro, L. (2013) ANOVA
kernels and RKHS of zero mean functions for model-based sensitivity
Journal of Multivariate Analysis, 115, 57–67.

[Ferrari and Dunson, 2020]
Ferrari, F. and Dunson, D. B. (2020a) Bayesian factor analysis
for inference on interactions.
JASA, 0.

[Ferrari and Dunson, 2020]
— (2020b) Identifying main effects and interactions among
exposures using Gaussian processes.
Ann. Appl. Stat., 14.

[Friedman, 1991]
Friedman, J. H. (1991) Multivariate adaptive regression splines.
Ann. Stat., 19, 1–67.

[George and McCulloch, 1993]
George, E. and McCulloch, R. (1993) Variable selection via Gibbs sampling.

[Greene et al., 2010]
Greene, C., Sinnott-Armstrong, N., Himmelstein, D. S., Park, P., Moore, J. and
Harris, B. (2010) Multifactor dimensionality reduction for graphics
processing units enables genome-wide testing of epistasis in sporadic ALS.
Bioinformatics, 26, 694–695.

[Griffin and Brown, 2017]
Griffin, J. and Brown, P. (2017) Hierarchical shrinkage priors for regression
Bayesian Analysis, 12, 135–159.

[Gu and Wahba, 1993]
Gu, C. and Wahba, G. (1993) Smoothing spline ANOVA with component-wise
Bayesian “confidence intervals”.
J. Comput. Graph. Stat., 2, 97–117.

[Gunn and Kandola, 2004]
Gunn, S. and Kandola, J. (2004) Structural modelling with sparse kernels.
Machine Learning.

[Hastie et al., 2001]
Hastie, T., Tibshirani, R. and Friedman, J. (2001) The Elements of
Statistical Learning.

[Hooker, 2007]
Hooker, G. (2007) Generalized functional anova diagnostics for high-dimensional
functions of dependent variables.
J. Comput. Graph. Stat., 16, 709–732.

[Hsu et al., 2014]
Hsu, D., Kakade, S. M. and Zhang, T. (2014) Random design analysis of ridge
Foundations of Computational Mathematics, 14,

[Huang, 1998]
Huang, J. Z. (1998) Projection estimation in multiple regression with
application to functional ANOVA models.
Ann. Stat., 26, 242–272.

[Lanckriet et al., 2004]
Lanckriet, G. R. G., Cristianini, N., Bartlett, P., Ghaoui, L. E. and Jordan,
M. I. (2004) Learning the kernel matrix with semidefinite programming.
Journnal of Machine Learning Research, 27–72.

[Lim and Hastie, 2015]
Lim, M. and Hastie, T. (2015) Learning interactions via hierarchical
group-lasso regularization.
J. Comput. Graph. Stat., 24, 627–654.

[Lin and Zhang, 2006]
Lin, Y. and Zhang, H. H. (2006) Component selection and smoothing in
multivariate nonparametric regression.
Ann. Stat., 34, 2272–2297.

[Liu et al., 2008]
Liu, H., Wasserman, L., Lafferty, J. D. and Ravikumar, P. K. (2008) Sparse
additive models.

[Lou et al., 2013]
Lou, Y., Caruana, R., Gehrke, J. and Hooker, G. (2013) Accurate intelligible
models with pairwise interactions.
In KDD.

[Maher, 2008]
Maher, B. (2008) Personal genomes: The case of the missing heritability.

[Micchelli et al., 2006]
Micchelli, C. A., Xu, Y. and Zhang, H. (2006) Universal kernels.
JMLR, 7, 2651–2667.

[Nakagawa et al., 2016]
Nakagawa, K., Suzumura, S., Karasuyama, M., Tsuda, K. and Takeuchi, I. (2016)
Safe pattern pruning: An efficient approach for predictive pattern mining.
In KDD.

[Rasmussen and Williams, 2006]
Rasmussen, C. E. and Williams, C. K. I. (2006) Gaussian Processes for
Machine Learning.
The MIT Press.

[Rendle, 2010]
Rendle, S. (2010) Factorization machines.

[Scheipl et al., 2012]
Scheipl, F., Fahrmeir, L. and Kneib, T. (2012) Spike-and-slab priors for
function selection in structured additive regression models.
JASA, 107, 1518–1532.

[Shah, 2016]
Shah, R. (2016) Modelling interactions in high-dimensional data with
JMLR, 17, 1–31.

[Slim et al., 2018]
Slim, L., Chatelain, C., Azencott, C. and Vert, J. (2018) Novel methods for
epistasis detection in genome-wide association studies.

[Stone, 1994]
Stone, C. J. (1994) The use of polynomial splines and their tensor products in
multivariate function estimation.
Ann. Stat., 22, 118–171.

[Vapnik, 1995]
Vapnik, V. N. (1995) The nature of statistical learning theory.

[Vitezica et al., 2013]
Vitezica, Z., Varona, L. and Legarra, A. (2013) On the additive and dominant
variance and covariance of individuals within the genomic selection scope.
Genetics, 195, 1223–1230.

[Wei et al., 2019]
Wei, R., Reich, B., Hoppin, J. and Ghosal, S. (2019) Sparse Bayesian additive
nonparametric regression with application to health effects of pesticides
Statistica Sinica.

Supporting Materials

§ ALGORITHMS DESCRIBED IN <REF>

§ FIGURES AND TABLES REFERENCED IN <REF>

§ PROOFS

§ LITERATURE REVIEW

§ ZERO MEAN KERNELS AND FINITE-BASIS FUNCTIONS

§ ADDITIONAL ALGORITHMIC DETAILS

§ MARS ANOVA PROCEDURE

§ ADDITIONAL EXPERIMENTAL DETAILS

§ ADDITIONAL EXPERIMENTAL RESULTS