turbulence well describes the turbulent-cascade process at play. For the
small- and medium-box cases, the physics involved in the nonlinear stage is
substantially different, and our estimate fails to accurately describe the
energy-injection process. In 3D, we expect the physics of this injection to
include additional dynamics along $y$, which allows $B_{y}$ to participate in
the cascade process. The absence of this mechanism in 2D is a potential source
of discrepancy between 2D and 3D results.

Figure 9: Estimate for injection of magnetic energy and conversion into
kinetic energy assuming Alfvénic turbulence, for the three simulations
presented in section 4.5. In each case, the end of the nonlinear phase is
marked with a black cross.

###  4.6 Stresses and angular-momentum transport in 2D

Figure 10: Left column: evolution of the $\alpha$-parameter and of the volume-
averaged Maxwell, Reynolds, and anisotropic stresses for a small box/large
$\omega_{\rm C}/\Omega_{0}$ (top), medium box/medium $\omega_{\rm
C}/\Omega_{0}$ (center), and large box/small $\omega_{\rm C}/\Omega_{0}$
(bottom) simulation of the 2D pair-plasma MRI. Right column: evolution of the
volume-averaged pressure anisotropy and $\beta_{\|}$, $\beta_{\perp}$ for the
same simulations.

In this section, we analyze the effective viscosity that originates from the
MRI nonlinear dynamics in 2D pair-plasma simulations. We define an effective-
viscosity parameter (Shakura & Sunyaev 1973),

$\alpha\equiv\frac{\langle R_{xy}+M_{xy}+A_{xy}\rangle}{\langle p\rangle},$
(63)

where $\langle...\rangle$ denotes a volume average. Here,
$R_{xy}=mnU_{x}U_{y}/\Gamma^{2}$ is the Reynolds stress, with
$\Gamma=\sqrt{1+U^{2}/c^{2}}$ and $\mbox{\boldmath{$U$}}=(U_{x},U_{y},U_{z})$
the spatial part of the bulk 4-velocity of the plasma in the frame moving with
velocity $\mbox{\boldmath{$v$}}_{\mathrm{s}}$; $M_{xy}=-B_{x}B_{y}/(4\pi)$ is
the Maxwell stress; and $A_{xy}=-(p_{\perp}-p_{\|})B_{x}B_{y}/B^{2}$ is the
anisotropic stress (e.g. Sharma et al. 2006). Our PIC simulations, where the
plasma can be relativistically hot, require special care when measuring
$\alpha$ and its components. Given a distribution function
$f(\mbox{\boldmath{$x$}},\mbox{\boldmath{$u$}},t)$, we first compute the bulk
velocity
$\mbox{\boldmath{$U$}}=(1/n)\int\mbox{\boldmath{$u$}}f(\mbox{\boldmath{$x$}},\mbox{\boldmath{$u$}},t)\textrm{d}^{3}\mbox{\boldmath{$u$}}$
for each particle species separately. Then, we apply a Lorentz boost to obtain
$\widetilde{\mbox{\boldmath{$u$}}}$, the spatial part of the 4-velocity of
each particle in the frame moving with $U$. We subsequently define a pressure
tensor for each particle species as

$\textbf{P}=m\int\frac{\widetilde{\mbox{\boldmath{$u$}}}\widetilde{\mbox{\boldmath{$u$}}}}{\widetilde{\gamma}}f(\mbox{\boldmath{$x$}},\widetilde{\mbox{\boldmath{$u$}}},t)\textrm{d}^{3}\widetilde{\mbox{\boldmath{$u$}}},$
(64)

where $\widetilde{\gamma}=\sqrt{1+\widetilde{u}^{2}/c^{2}}$. The isotropic,
parallel, and perpendicular pressure are then $p=\mathrm{tr}(\textbf{P})/3$,
$p_{\|}=\textbf{P}\mbox{\boldmath{$\,\mbox{\boldmath{$:$}}\,$}}\hat{\mbox{\boldmath{$b$}}}\hat{\mbox{\boldmath{$b$}}}$,
and
$p_{\perp}=(1/2)\textbf{P}\,\mbox{\boldmath{$:$}}\,(\mathbb{I}-\hat{\mbox{\boldmath{$b$}}}\hat{\mbox{\boldmath{$b$}}})$
respectively, where
$\hat{\mbox{\boldmath{$b$}}}=\widetilde{\mbox{\boldmath{$B$}}}/\widetilde{B}$
and $\widetilde{\mbox{\boldmath{$B$}}}$ is the magnetic field in the frame
moving with $U$. $R_{xy}$ and $A_{xy}$ are computed for each species
separately, and then summed across species.

In Figure 10, we show the time evolution of the average $\alpha$ and its
components (left column) as well as the volume-averaged
$\langle\beta_{\|}\rangle$, $\langle\beta_{\perp}\rangle$, and the pressure
anisotropy $\langle(p_{\perp}-p_{\|})/p\rangle$ (right column) for the small-
box/high-frequency ratio case $2\times 2\lambda_{\mathrm{MRI}}^{2}$,
$\omega_{\rm C}/\Omega_{0}=240$ (top row), the medium-box/medium-frequency
ratio $8\times 8\lambda_{\mathrm{MRI}}^{2}$, $\omega_{\rm C}/\Omega_{0}=60$
(middle row), and the large-box/small-frequency ratio $32\times
32\lambda_{\mathrm{MRI}}^{2}$, $\omega_{\rm C}/\Omega_{0}=30$ (bottom row). In
all cases, $\beta_{0}=624$ and $v_{\mathrm{A},0}/c=0.01$. All simulations are
characterized by a rapid increase in all components of $\alpha$ during the
linear stage, with a peak at $\alpha\simeq 1$ right before the onset of
reconnection. At this point, different parameters determine strong differences
in the results: during the nonlinear stage, we observe that the volume-
averaged $M_{xy}$ can grow and remain close to 1 at early times in the small-
and medium-box cases; conversely, in the large-box case $M_{xy}$ rapidly
settles to $\sim 0.1$ after the onset of reconnection and slowly decreases to
$\sim 0.01$. In all cases, the anisotropic stress $A_{xy}$ achieves the
highest (volume-averaged) value $\sim 0.1$ at the onset of reconnection, and
then firmly settles around $0.02\mbox{--}0.05$ for the whole duration of the
nonlinear phase. The Reynolds stress is consistently the least important
component of $\alpha$ and remains much smaller than the other two during the
nonlinear phase. Overall, the $\alpha$-parameter closely follows the evolution
of $M_{xy}$ in the small- and medium-box cases, attaining values
$0.1\mbox{--}1$; in the large-box case, $\alpha\simeq 0.01\mbox{--}0.1$ and
the contribution of $A_{xy}$ is more significant, providing most of the
angular-momentum transport in the late nonlinear stage when the magnetic field
is dissipating via reconnection and turbulence. We note, however, that in the
large-box case $\omega_{\rm C}/\Omega_{0}$ is substantially smaller than in
the small-box case (respectively $\omega_{\rm C}/\Omega_{0}=30$ and
$\omega_{\rm C}/\Omega_{0}=240$); since the scale separation affects the
development of mirror modes (see section 4.2), pressure anisotropy may be
exaggerated in the small frequency-ratio case. Larger values of $\omega_{\rm
C}/\Omega_{0}$ would be needed to verify the trend of the anisotropic stress
in the large-box case, but this involves computational costs beyond our
current possibilities. We will investigate the matter more in detail in future
work.

These results align well with the discussion provided in the previous
sections. The MRI in the large-box case behaves qualitatively differently from
the small-box case: the system develops strong turbulence during the nonlinear
stage, promoting the continuous, more efficient dissipation of magnetic fields
and faster pressure isotropization. In small boxes, this dissipation is
primarily mediated by episodic, large-scale reconnection events, with phases
of strong magnetic-field amplification driving pressure anisotropy in between.
This suggests that, in a large box, $M_{xy}$ and $A_{xy}$ will be reduced in
comparison with a small-box case. The right-hand panels of Figure 10 confirm
this reasoning: pressure anisotropy is much more pronounced at the beginning
of the nonlinear stage in smaller boxes, while it remains rather limited in
the large-box case.

If $\alpha$ is invoked as a coefficient of turbulent angular-momentum
transport, then the large-box case where turbulence can develop is more
relevant to provide a first-principles measure of the effective turbulent
viscosity in a collisionless plasma. In this case, our 2D numerical
experiments indicate that the Maxwell stress provides the primary contribution
to $\alpha$, with the anisotropic stress achieving smaller values
$A_{xy}/M_{xy}\simeq 0.2$. This result differs from the findings of Kunz et
al. (2016), who reported $A_{xy}/M_{xy}\simeq 1$; however, we stress again
that 2D runs have the intrinsic limitation of lacking reconnection processes
in the dominant magnetic-field component $B_{y}$. In 3D, we may then expect
large differences with the results presented here: dissipating $B_{y}$ via
reconnection may imply that the Maxwell stress could achieve much smaller
values than those measured in 2D, possibly increasing the relative
contribution of pressure anisotropy in determining $\alpha$. On the other
hand, in the previous sections we have argued that reconnection along $y$ may
help reduce pressure anisotropy by increasing $p_{\|}$, hence decreasing
$A_{xy}$. We will discuss the dynamics of angular-momentum transport in 3D in
section 5.4.

###  4.7 Net flux vs. zero net flux

Figure 11: Evolution of the in-plane magnetic-energy variation for several 2D
pair-plasma simulations of the MRI with zero net magnetic flux. Left:
simulations with varying box size and fixed $\omega_{\rm C}/\Omega_{0}$,
$\beta_{0}$. Right: simulations with varying $\omega_{\rm C}/\Omega_{0}$ and
fixed box size and $\beta_{0}$.

We conclude our investigation of the 2D collisionless-MRI dynamics by briefly
considering the case of zero net flux. In contrast with the simulations
presented in the previous sections, here the initial magnetic field is set to
$\mbox{\boldmath{$B$}}(t=0)=-B_{z,0}\sin\left(2\pi
x/L_{x}\right)\hat{\mbox{\boldmath{$e$}}}_{z}$ such that no net flux is
present along $z$. This case was briefly discussed in Riquelme et al. (2012)
in 2D and was the focus of the only (to date) 3D hybrid study of the MRI (Kunz
et al. 2016).

In Figure 11 we plot the evolution of the in-plane magnetic energy for a few
simulations with exactly the same parameters employed in section 4.3 (aside
from the initial magnetic-field configuration). In particular, we compare
$B_{\mathrm{pol}}$ in several cases: first, by keeping a fixed $\omega_{\rm
C}/\Omega_{0}=30$, $\beta_{0}=624$, and $v_{\mathrm{A},0}/c=0.01$ (the latter
two measured with the maximum value of $B$), in larger and larger simulation
boxes $L_{z}/\lambda_{\mathrm{MRI}}=L_{x}/\lambda_{\mathrm{MRI}}=4,8,16,32$
(left panel); then, by fixing the box size
$L_{z}/\lambda_{\mathrm{MRI}}=L_{x}/\lambda_{\mathrm{MRI}}=8$,
$\beta_{0}=624$, and $v_{\mathrm{A},0}/c=0.01$, and increasing the frequency
ratio $\omega_{\rm C}/\Omega_{0}=6,15,30,60$ (right panel). From this
comparison, we observe that the general statements made for the nonzero net-
flux case hold here as well: the duration of the turbulent, saturated
nonlinear stage that follows the linear growth of the MRI increases both with
box size and with $\omega_{\rm C}/\Omega_{0}$. Over time, the turbulence
rapidly decays as magnetic fields grow and pressure anisotropy develops (not
shown), pushing the MRI wavelength to scales larger than the box size.
However, the zero net-flux case is also characterized by striking differences,
with the most prominent being a much more efficient breakup and mixing of
channel flows. This was also observed by Riquelme et al. (2012), who
attributed this behavior to the fact that the instability develops
nonuniformly in different regions of the domain; in particular, the most-
unstable MRI wavelength $\lambda_{\mathrm{MRI}}=2\pi
v_{\mathrm{A}}/\Omega_{0}$ is shorter where $B_{z}$ is weaker and longer where
$B_{z}$ is stronger, producing a spatial inhomogeneity in the dominant
instability wavevector. Channels forming inside regions of weaker/stronger
$B_{z}$ are therefore spatially nonuniform, and tend to be destroyed much more
quickly. Hence, differently from the net-flux case, the break-up of channels
formed during the linear stage is not solely mediated by the onset of tearing
modes (and concurrently of kink instabilities): the dynamics does not need to
“wait” for the magnetic field in the channels to grow strong enough (and for
current sheets to thin out enough) that reconnection can occur. Instead, the
intrinsic inhomogeneity of the background field mediates the mixing and
redistribution of magnetic energy. As a consequence, the saturation level
reached in all simulations with zero net flux is much smaller than in the
corresponding nonzero flux cases.

Figure 11 also shows important differences in how these simulations tend to
converge: while it appears that progressively better convergence is achieved
for increasing $\omega_{\rm C}/\Omega_{0}$ at a fixed box size, enlarging the
box size does not result in converging results. Larger boxes seemingly produce
a steady increase in the magnetic-field amplification, in contrast with the
clear convergence observed for the net-flux case. We also note that, since the
magnetic field here reverses its direction through the box, we observe no sign
of the differential heating of electrons vs. positrons that we measured in the
net-flux case (see section 4.5). Charges with opposite sign are indeed equally
likely to experience differential gyromotion and drift towards current sheets,
overall resulting in an almost perfectly balanced energization among electrons
and positrons.

In this work, we focus on nonzero net-flux simulations primarily because, in
any astrophysical accretion disk, the complete absence of vertical magnetic
flux is not expected; realistically, any portion of the disk will be threaded
by some magnetic field carrying a net flux, however small. Additionally, the
behavior of zero net-flux simulations is fundamentally different from the
expected development of a “pure” MRI: the former situation is more akin to a
dynamo problem, and does not feature the same qualitative evolution of the
nonzero net-flux case. Our results show that the case of zero net flux differs
enough from the case with nonzero magnetic flux to deserve a thorough,
dedicated exploration, which we will pursue in future work.

##  5 Three-dimensional simulations of the MRI in collisionless pair plasmas

Table 2: List of 3D simulations

Box size $(L_{x}\times L_{y}\times L_{z})$ | $\omega_{\mathrm{C}}/\Omega_{0}$ | $\beta_{0}$ | $\lambda_{\mathrm{MRI}}/\rho_{\mathrm{C}}$
---|---|---|---
$2\times 2\times 2\lambda_{\mathrm{MRI}}^{3}$ | 15 | 312 | 7.46
$4\times 4\times 2\lambda_{\mathrm{MRI}}^{3}$ | 15 | 312 | 7.46
$4\times 8\times 2\lambda_{\mathrm{MRI}}^{3}$ | 15 | 312 | 7.46

Note. — In all cases, $v_{\mathrm{A},0}/c\simeq 0.007$; the numerical
resolution is such that $\Delta x\simeq c/\omega_{\mathrm{p}}$, and we employ
27 particles per cell per species.

In this section we present 3D PIC simulations conducted with Zeltron, where we
have implemented the newly developed KSB-OA model (see section 2.3). In
general, we expect 3D runs to produce results more physically realistic (at
least in terms of the magnetic-field evolution and subsequent development of
turbulence) than corresponding 2D simulations, for the reasons discussed in
section 4.2. We will therefore select simulation parameters that strictly
enforce three physical constraints: i) a sufficiently subluminal velocity
offset at the $x$-boundaries, i.e. $v_{\mathrm{s}}(L_{x})=s\Omega_{0}L_{x}<c$;
ii) a vertical extent of the simulation box that does not exceed the disk
pressure scale height, i.e. $L_{z}/H\leq 1$; and iii) a sufficiently large
separation between macroscopic and microscopic scales (at least in the
beginning of each run), in terms of
$\lambda_{\mathrm{MRI}}/\rho_{\textrm{C}}$.

Due to the much higher computational cost of 3D runs, in this first study we
will explore a limited parameter range inspired by the results of 2D
simulations. We will show that 3D runs of relatively small size already
provide fundamental new insight into the physics of the collisionless MRI. The
3D simulations discussed in this work are summarized in Table 2; all runs are
initialized analogously to the 2D case, with a weak vertical magnetic field
determined from the initial simulation parameters. We always employ 27
particles per cell per species, and a numerical grid such that $\Delta x\simeq
c/\omega_{\mathrm{p}}$. In future work, we will conduct larger simulations
with better scale separation to achieve more realistic regimes.

In the following sections, we will first describe the physics of small- and
large-box 3D PIC simulations of the pair-plasma MRI, and how these relate to
corresponding 2D simulations. We will then analyze the 3D turbulence, particle
energization, and angular-momentum transport that characterize the 3D case.

###  5.1 Physics of the 3D collisionless MRI

Figures 12 and 13 show the results for two representative 3D simulations of
the collisionless MRI in pair plasmas, employing a small ($2\times 2\times
2\lambda_{\mathrm{MRI}}^{3}$) and a large ($4\times 8\times
2\lambda_{\mathrm{MRI}}^{3}$) box size, respectively. In these figures, we
show subsequent snapshots of the spatial distribution of $B_{x}$ and $B_{y}$
(left panels), the evolution of the change in magnetic energy (top-right
panel), and of the volume-averaged parallel and perpendicular plasma-$\beta$
(bottom-right panel). These simulations are initialized with parameters such
that $L_{z}=H$, $\lambda_{\mathrm{MRI}}/\rho_{\textrm{C}}\simeq 7.5$, and
$v_{\mathrm{s}}(L_{x})/c\simeq 0.125$ (in the small-box case) and
$v_{\mathrm{s}}(L_{x})/c\simeq 0.25$ (in the large-box case). In both cases,
imposing these conditions results in a frequency ratio $\omega_{\rm
C}/\Omega_{0}=15$ and a plasma temperature $\theta_{0}=1/128$ such that
$\beta_{0}=312$ (including both species), with an initial Alfvén speed
$v_{\mathrm{A},0}/c\simeq 0.007$.

Figure 12: Evolution of the pair-plasma MRI in a small 3D box ($2\times
2\times 2\lambda_{\mathrm{MRI}}^{3}$). Left: spatial distribution of the
radial ($B_{x}$) and toroidal ($B_{y}$) magnetic field at subsequent times
during the simulation. Right: evolution of the change in magnetic energy in
all three components of $B$ (top) and evolution of the volume-averaged
$\beta_{\|}$ and $\beta_{\perp}$ (bottom). Figure 13: As in Figure 12 but for
a box size of $4\times 8\times 2\lambda_{\mathrm{MRI}}^{3}$.

For the small-box case ($2\times 2\times 2\lambda_{\mathrm{MRI}}^{3}$) shown
in Figure 12, we observe that the evolution of the magnetic energy (top-right
panel) is initially qualitatively similar to that of 2D simulations of
comparable $L_{x},L_{z}$ (cf. the 2D case of Figure 2): from the initial
state, the field strength grows exponentially until saturation is reached. The
nonlinear stage, however, shows qualitative differences: we observe that the
magnetic energy in all components (dominated by $B_{y}$, which does not
saturate in 2D) remains approximately around the saturation level for $\sim 2$
orbital periods after the end of the linear stage, before starting to slowly
decrease. This behavior is due to the 3D geometry and is illustrated by
subsequent snapshots of the $B_{x}$ and $B_{y}$ distribution in the domain,
shown in Figure 12 (left panels). Around $t\simeq 3.7P_{0}$ (top-left panels)
we observe that channel flows form, during the linear stage, in the $xz$-plane
over length scales of $\lambda_{\mathrm{MRI}}$, similarly to the 2D dynamics.
However, in 3D these channels also span the entire $y$-extent of the
simulation domain. At the channel interfaces, reconnection starts dissipating
magnetic energy when current sheets become sufficiently thin and tearing modes
are excited (around $t\simeq 4.2P_{0}$, center-left panels). In 3D,
reconnection can take place also along $y$, thus dissipating energy stored in
$B_{y}$ (and not just in $B_{x}$ and $B_{z}$ as in 2D). The channel interfaces
are indeed heavily perturbed and magnetic reconnection can act along all
spatial directions (e.g. around $t\simeq 5.4P_{0}$, bottom-left panels). The
saturated energy level of the nonlinear stage we measure between $t\simeq
4P_{0}$ and $t\simeq 6P_{0}$ is maintained due to the balance achieved between
injection (via the MRI) and dissipation (via reconnection) of energy into the
system.

Starting from $t\simeq 6P_{0}$, the system slowly relaxes toward an end state
where the MRI slows down and then stops. This final state (not shown in the
left panels of Figure 12) is a two-channel configuration with lower $B_{y}$
energy. From $t\simeq 8P_{0}$ onward, we observe that the energy in $B_{x}$
and $B_{z}$ oscillates around approximately the same values; the spatial
distribution of the magnetic field tends to assume a “shear-wave”
configuration with coherent, large-scale structures aligned with the shear
flow appearing cyclically. No reconnection occurs at this point; the system
appears to maintain this state endlessly. The dynamics that characterizes this
state is different from that of the early nonlinear stage, and it is doubtful
whether it can even be associated with the MRI at all. Qualitatively, it
appears that MRI modes have migrated to wavelengths longer than the box size,
and that the 3D geometry is allowing a quasi-steady state that consists of
periodic oscillations. This dynamics is of no relevance for this work, but it
is interesting to note the difference between the 3D and 2D end states, with
the latter consisting of stable loops in the $xz$-plane.

The evolution of the volume-averaged $\beta_{\|}$ and $\beta_{\perp}$, shown
in Figure 12 (bottom right), closely follows that of the magnetic field: in
the early nonlinear stage (between $t\simeq 4P_{0}$ and $t\simeq 6P_{0}$), we
measure an increase
$\langle\beta_{\|}\rangle\simeq\langle\beta_{\perp}\rangle\sim
0.2\mbox{--}10$, with very brief periods where
$\langle\beta_{\perp}\rangle>\langle\beta_{\|}\rangle$ and subsequent “bursty”
reconnection events reestablishing isotropy periodically. At late times,
pressure anisotropy accumulates and pushes the MRI growth rate to smaller and
smaller values.

The small 3D run discussed above achieves magnetic-energy saturation in the
nonlinear stage, thus already improving over 2D runs of much larger size (in
$x$ and $z$). However, in such as small simulation we do not obtain a
developed turbulent state, consistently with the results of the sole 3D fully
kinetic study previously presented (Hoshino 2015). For this reason, we now
focus on the larger simulation shown in Figure 13 with box size $4\times
8\times 2\lambda_{\mathrm{MRI}}^{3}$. This case is chosen as representative of
a situation where the numerical domain size is rather limited in $z$ (in our
case, by the requirement $L_{z}\leq H$) — hence a small number of channel
modes are expected to grow during the linear stage — but the larger extent in
$x$ and $y$ now allows for additional instability mechanisms which cannot be
included in the previous $2\times 2\times 2\lambda_{\mathrm{MRI}}^{3}$ case.
For 2D simulations, in section 4.3.4 we have argued that boxes elongated in
$x$ promote the development of turbulence by allowing for drift-kink modes
that can disrupt the channels; here, we observe exactly this effect in 3D. The
top-right panel of Figure 13 shows that the system evolution now proceeds
through a linear stage where the magnetic field is amplified roughly 10 times
less than in the previous cubic-box case. The onset of magnetic reconnection
results in a nonlinear stage that is qualitatively different from that of the
small-box case: repeated cycles of period $\sim 2P_{0}$ are observed where the
magnetic energy in all components grows and is rapidly dissipated before
growing again. Each growth phase corresponds to the formation of macroscopic
channels, and each dissipation phase results in the complete breakup of said
channels, which lose coherence and decay into “turbulent” structures (whether
this can be truly called a turbulent state will be discussed in section 5.2).
This behavior is visible in the spatial distribution of $B_{x}$ and $B_{y}$ at
subsequent times in Figure 13 (left panels).

Reconnection in this elongated-box simulation is evidently much more active
than in the cubic-box case. We observe that the macroscopic current sheets
violently kink in both $x$ and $y$, eventually resulting in the emergence of
the turbulent state. The bottom-right panel of Figure 13 shows that the
stronger reconnection dynamics also heavily impacts the evolution of $\beta$:
because of the smaller magnetic-field amplification (with respect to the
cubic-box case), the nonlinear stage is characterized by a volume-averaged
$\langle\beta\rangle\sim 10$, which can grow up to $\sim 100$ during the
dissipation phase of the channel-breakup cycles. Pressure anisotropy is well
regulated, during the early nonlinear stage, by subsequent reconnection
events, which periodically reestablish isotropy; in the late nonlinear stage,
a finite (tens of percent), unquenched pressure anisotropy is retained,
marking the transition to the same end state described above for the small-box
simulation.

The numerical experiment reported in Figure 13, however limited in box size
and scale separation, provides a first, striking evidence that large 3D PIC
simulations of the collisionless MRI possess qualitative features analogous to
corresponding MHD simulations. In contrast to the 2D case, where important
physical processes are impeded by the reduced dimensionality, here we observe
an MRI behavior akin to that exhibited by 3D MHD runs, with the cyclic
creation (via the MRI) and destruction (via reconnection) of channel flows. In
the next subsections, we discuss more quantitative comparisons with 2D
simulations as well as the effect of choosing different aspect ratios in 3D
runs.

####  5.1.1 2D vs. 3D MRI dynamics

Figure 14: Left: evolution of the volume-averaged Alfvén speed (calculated
with the nonrelativistic expression) for several 3D PIC simulations of the
pair-plasma MRI, compared with a 2D simulation using the same initial physical
parameters. In the 3D runs, the constraint $v_{\rm A}^{\mathrm{NR}}/c<1$ is
much better respected. Right: evolution of the volume-averaged pressure
anisotropy, which influences the MRI growth rate, in the same runs. The 3D
runs develop smaller anisotropy over longer times.

Figure 14 shows the evolution of the volume-averaged Alfvén speed (left
panel), calculated with the nonrelativistic expression, and of the volume-
averaged pressure anisotropy $(p_{\perp}-p_{\|})/p$ (right panel) for the two
runs presented in the previous section (with box size $2\times 2\times
2\lambda_{\mathrm{MRI}}^{3}$ and $4\times 8\times
2\lambda_{\mathrm{MRI}}^{3}$, respectively), as well as for an analogous run
with box size $4\times 4\times 2\lambda_{\mathrm{MRI}}^{3}$. In the same
plots, these 3D runs are compared with a small-box 2D simulation (with box
size $2\times 2\lambda_{\mathrm{MRI}}^{2}$) employing the same physical
parameters. We observe that the 3D simulations display a qualitatively
different behavior overall.

As was also found in earlier work (Riquelme et al. 2012; Inchingolo et al.
2018), in the 2D run the average $v_{\textrm{A}}^{\mathrm{NR}}$ vastly exceeds
the speed of light during the nonlinear stage. For 2D, Inchingolo et al.
(2018) attributed this behavior to the limited size of the simulation box,
showing that for larger domains $\gtrsim 8\times 8\lambda_{\mathrm{MRI}}^{2}$
the constraint $v_{\mathrm{A}}^{\mathrm{NR}}/c<1$ was fulfilled over short
times. However, our 2D campaign (section 4) shows that in all cases
(irrespective of the box size) the system inevitably evolves such that
$v_{\mathrm{A}}^{\mathrm{NR}}/c>1$ over sufficiently long times. This occurs
because the out-of-plane $B_{y}$ fails to saturate in 2D, due to a lack of
dissipation processes (i.e. reconnection) along $y$ (see section 4.2). The
violation of the $v_{\mathrm{A}}^{\mathrm{NR}}/c<1$ condition breaks the
underlying nonrelativistic shearing-box assumption, and therefore puts into
question the physical validity of the 2D results. In addition, the lack of
saturation in 2D heavily contributes to the migration of MRI modes to scales
larger than the box size. Our 3D simulations shown in Figure 14 (left panel)
display a different behavior, with the magnetic-field amplification (and
therefore $v_{\mathrm{A}}^{\mathrm{NR}}/c$) reaching much smaller values. Even
a limited-size, cubic-box run already results in
$v_{\mathrm{A}}^{\mathrm{NR}}/c\sim 1$ (i.e. $\sim 10$ times smaller than in
the corresponding 2D case), and larger and larger 3D boxes progressively
improve this result. Our largest two cases, with box size $4\times 4\times
2\lambda_{\mathrm{MRI}}^{3}$ and $4\times 8\times
2\lambda_{\mathrm{MRI}}^{3}$, firmly maintain
$v_{\mathrm{A}}^{\mathrm{NR}}/c\lesssim 1$ for the whole duration of the
nonlinear stage.

The fact that in 3D MRI saturates at lower magnetic-field amplitudes implies
that unstable MRI wavelengths do not migrate to larger length scales as
quickly as in 2D. Thanks to the more efficient reconnection during the
nonlinear stage, in 3D the MRI is not aggressively suppressed by the unbound
growth of magnetic fields. This supports our statement on the significance of
3D simulations: capturing the MRI and reconnection dynamics along the toroidal
direction is essential to model a physically meaningful evolution of the
instability, which can then maintain a saturated state for longer times.

2D simulations also display a qualitatively different behavior in the
development of pressure anisotropy: Figure 14 (right panel) shows the
evolution of the volume-averaged $(p_{\perp}-p_{\|})/p$ in the same 2D and 3D
simulations. In 2D, the anisotropy rapidly increases during the nonlinear
stage, and progressively accumulates as the MRI continues to amplify $B_{y}$
without the dissipative effect of reconnection along $y$, up to a maximum
$\langle(p_{\perp}-p_{\|})/p\rangle\simeq 1$. In all 3D runs, this anisotropy
is limited to smaller values (up to a maximum
$\langle(p_{\perp}-p_{\|})/p\rangle\simeq 0.5$) for the whole duration of the
nonlinear phase, mildly increasing toward the end of the run. This is verified
both in the cubic-box case and in larger boxes, indicating that the additional
reconnection along $y$ plays a major role in constraining the development of
pressure anisotropy. This anisotropy, eventually, contributes to pushing MRI
wavelength to larger scales and to halt the MRI development, but over longer
timescales than in corresponding 2D cases.

####  5.1.2 Effect of the aspect ratio

The 2:4:1 aspect ratio employed for the simulation presented in Figure 13 has
been chosen by carefully analyzing a number of runs with varying
$L_{x},L_{y},L_{z}$. As illustrated in section 3, the physical domain size in
$z$ (i.e. $L_{z}/\lambda_{\mathrm{MRI}}$) heavily impacts the simulation cost,
when demanding that the macroscopic scale ordering $H/L_{z}\geq 1$ be
respected. In practice, choosing $L_{z}/\lambda_{\mathrm{MRI}}$ and
$\lambda_{\mathrm{MRI}}/\rho_{\mathrm{C}}$ determines the required number of
grid cells per $\lambda_{\mathrm{MRI}}$ and the number of time steps per
$P_{0}$. Then, the aspect ratios $L_{x}/L_{z}$ and $L_{y}/L_{z}$ are free
parameters, each impacting the overall simulation cost linearly (i.e.,
quadratically when combined). Ideally, we are free to choose a small vertical
extent, e.g. $L_{z}/\lambda_{\mathrm{MRI}}=2$, and then elongate the box in
$x$ and $y$ to better include the reconnection physics (as shown in sections
4.3.4 and 5.1). However, the $x$-extent of the box cannot be increased
indefinitely, since we demand that the condition $v_{\mathrm{s}}(L_{x})/c<1$
be satisfied at the radial boundaries. Our tests indicated that numerical
artifacts appear at the boundary when $L_{x}/L_{z}\geq 4$ (corresponding to
$v_{\mathrm{s}}(L_{x})/c\geq 0.5$). Hence, we choose $L_{x}/L_{z}=2$ in our 3D
simulations. Elongating the box in $y$ is less problematic, since no
constraints related to the nonrelativistic shearing-box model exist in the
toroidal direction111111Note however that curvature effects have been
neglected when constructing the shearing-box equations; $L_{y}/L_{z}$ can then
be larger than $L_{x}/L_{z}$, but must remain limited to a certain extent.. We
have chosen $L_{y}/L_{z}=4$ in our largest run.

In earlier sections, we have argued that boxes with larger $L_{x}/L_{z}$
produce better developed turbulence; Figure 14 indeed demonstrates that as we
enlarge a 1:1:1 box to 2:2:1 and then to 2:4:1, we can observe the emergence
of the characteristic cycles of creation and destruction of channel flows.
Reconnection acts more efficiently when $L_{x}/L_{z}>1$, resulting in dynamics
more similar to that expected from 3D MHD numerical experiments. It remains to
be assessed whether the effect of the aspect ratio changes when macroscopic
and microscopic length scales are better separated, and whether even larger
$L_{y}/L_{z}$ produce significant changes in the results. Due to the
exceptional computing costs involved for such simulations, we defer this
analysis to future work.

###  5.2 Spectra of 3D MRI-driven turbulence

Figure 15: Isotropic power spectrum of the poloidal (red line) and toroidal
(blue line) magnetic field for the 3D pair-plasma MRI simulation of Figure 13
($\omega_{\rm C}/\Omega_{0}=15$, $\beta_{0}=312$, and box size $4\times
8\times 2\lambda_{\mathrm{MRI}}^{3}$). The spectra are averaged over one
channel disruption/creation cycle during the nonlinear stage (between $t\simeq
9P_{0}$ and $t\simeq 11P_{0}$), showing a persisting peak around the most-
unstable MRI wavenumber $k_{\mathrm{MRI}}$. The spectra feature characteristic
slopes indicating an inertial range, and a spectral break in the vicinity of
the average Larmor-radius wavenumber $k_{\rho}$ measured at $t\simeq 10P_{0}$.

In Figure 15 we show the isotropic power spectrum of the poloidal
($B_{\mathrm{pol}}=\sqrt{B_{x}^{2}+B_{z}^{2}}$) and toroidal ($B_{y}$)
magnetic-field components for the large-box simulation discussed in Section
5.1 (with parameters $\omega_{\rm C}/\Omega_{0}=15$, $\beta_{0}=312$, and box
size $4\times 8\times 2\lambda_{\mathrm{MRI}}^{3}$; see Figure 13). This run
is characterized by an “episodic” behavior during the nonlinear stage, where
cycles of channel creation/destruction can be clearly identified. The
magnetic-field strength and configuration change drastically throughout each
cycle; for this reason, in Figure 15 we focus on the cycle occurring toward
the end of the nonlinear stage (between $t\simeq 9P_{0}$ and $t\simeq
11P_{0}$, see Figure 13), and time-average the power spectra throughout the
cycle.

The results show a persistent peak in the spectra around $k_{\mathrm{MRI}}$,
the wavenumber of the most-unstable MRI mode calculated with the initial
simulation parameters. This peak is indicative of the fact that the MRI is
active through the cycle, although magnetic-field dissipation via reconnection
is clearly dominant during the dissipation phase, when channels are being
disrupted and magnetic energy is converted into heat. Such dissipation results
in a turbulent cascade, which can be identified in the spectra via the
presence of an inertial range with characteristic slopes. At moderate
wavenumbers, the poloidal-field spectrum features a shallow slope which could
be comparable to $k^{-5/3}$, although the spectra are quite noisy (especially
in the inertial range) and a precise measurement of the spectral slopes is
hard to obtain. Regardless of the exact power-law index, we clearly observe
that a spectral break occurs around $k\simeq k_{\rho}$, with $k_{\rho}$ the
wavenumber corresponding to the average Larmor radius measured at $t\simeq
10P_{0}$ (i.e. at the midpoint in time through the cycle). In the sub-Larmor
range ($k>k_{\rho}$), the $B_{\mathrm{pol}}$ spectrum steepens to a slope
$\propto k^{-3.2}$; in Figure 15, a slope $\propto k^{-3}$ is also indicated
to guide the eye. The toroidal-field spectrum similarly features a shallow
slope consistent with $k^{-2}$ at $k<k_{\rho}$ (although this part of the
spectrum is also rather noisy), and a spectral break and steepening to
$\propto k^{-3.2}$ in the kinetic range.

These 3D spectra present similarities as well as substantial differences with
those shown in Figure 7 for the 2D case. Like in 2D, for the poloidal field we
observe an inertial range with a $\propto k^{-5/3}$ slope, indicative of a
turbulent MHD-type cascade driven by energy injection at large scales;
similarly to 2D, for $B_{\mathrm{pol}}$ we measure a spectral break and a
steepening in the sub-Larmor range, although in 3D this part of the spectra
features a slightly steeper $\propto k^{-3.2}$ slope. Especially at kinetic
scales, 2D and 3D spectra show substantial differences: the $B_{y}$ spectrum
in 3D presents a spectral break not observed in 2D, and the kinetic range is
here much more extended. In addition, we pointed out that in 2D the spectra
continuously evolve and, over time, the measured characteristic slopes do not
persist. We believe that this is an artifact of the reduced 2D dimensionality;
in 3D, we have verified that the spectra shown in Figure 15 consistently
appear over each channel creation/disruption cycle in the nonlinear stage as
long as the MRI is active.

Our 3D results are in good agreement, at least in the inertial range
($k<k_{\rho}$), with the 3D hybrid simulations of Kunz et al. (2016) as well
as with 3D MHD results (e.g. Walker et al. 2016). Slight differences with
these works arise in the kinetic range ($k>k_{\rho}$), where we measure a
$\propto k^{-3.2}$ slope instead of the $\propto k^{-3}$ reported by Kunz et
al. (2016), whose result was however related to the ion sub-Larmor scales in a
hybrid simulation, while here we are considering fully kinetic pair plasmas.
Future 3D PIC runs with electron-ion plasma will clarify whether this
discrepancy is related to our pair-plasma assumption.

###  5.3 Particle energization in 3D

Figure 16: Left: evolution of the average Lorentz factor for electrons,
positrons, and for the two species combined in the 3D pair-plasma MRI
simulation of Figure 13 ($\omega_{\rm C}/\Omega_{0}=15$, $\beta_{0}=312$, and
box size $4\times 8\times 2\lambda_{\mathrm{MRI}}^{3}$). As in 2D, electrons
heat significantly more than positrons. Right: evolution of the particle
energy distribution during the nonlinear stage ($t\simeq 4.7\mbox{--}15P_{0}$)
of the same simulation. Right after the onset of reconnection at $t\simeq
4.7P_{0}$, the distribution has developed a nonthermal tail following a power
law (the slopes of power laws with index 2.5 and 3 are shown to guide the
eye). The nonthermal population persists until the end of the run, while a
separated high-energy peak develops around $\gamma\sim 100$ (dotted line),
which corresponds to the maximum attainable energy due to the system size in
the $z$-direction.

In this section, we briefly analyze particle energization in the 3D pair-
plasma MRI. We focus on our largest simulation presented in Figure 13
($\omega_{\rm C}/\Omega_{0}=15$, $\beta_{0}=312$, and box size $4\times
8\times 2\lambda_{\mathrm{MRI}}^{3}$) and measure the evolution of the mean
Lorentz factor for electrons and positrons separately, as well as the average
value for both species. In addition, we compute the energy distribution
function at subsequent times during the simulation. The results are shown in
Figure 16.

Much like in the 2D case (section 4.5), we observe that the particle energy
constantly increases throughout the run, with the average Lorentz factor
reaching values $\sim 10\mbox{--}40$ at the end of the simulation (left
panel). Also similarly to our 2D simulations, we measure significant
differential heating between the two particle species, with electrons reaching
energies $\sim 4$ times larger than positrons on average. The energy
distribution functions measured for both species together (right panel) show
that right after the onset of reconnection at $t\simeq 4.7P_{0}$, there is a
significant nonthermal component in the energy spectra; this nonthermal tail
follows a power law that persists until the end of the run as the plasma
heats, while an additional high-energy peak progressively develops around
$\gamma\simeq 100$ (indicated by a dashed line in Figure 16). At the end of
the simulation ($t\simeq 15P_{0}$), particles have accumulated around the peak
energy, eroding the power-law tail of the spectrum.

These 3D results show both remarkable similarities and qualitative differences
with our 2D runs. First, the differential heating already found in 2D
simulations arises in 3D essentially in the same way, with electrons gaining
more energy than positrons; the mechanism driving this effect is the same as
in 2D simulations (on which we elaborate in appendix B). The plasma heating
proceeds, on average, very similarly to our large-box 2D runs (see Figure 8),
where developed Alfvénic turbulence drives a rather smooth increase in
particle kinetic energy. In 3D, the “episodic” nature of the channel
creation/disruption cycles manifests itself in the evolution of
$\langle\gamma\rangle$, where we observe periodic (small) increases
corresponding to the repeated onset of reconnection events. The particle
energy distribution also evolves similarly to the 2D case, at least in the
beginning of the nonlinear stage; like in 2D, the slope of the power-law tail
is $\sim 2.5\mbox{--}3$, persisting until the end of the run. A qualitative
difference with 2D simulations lies in the end state reached in 3D runs: the
high-energy peak shown in Figure 16 here precisely corresponds to the Lorentz
factor of particles with $\rho_{\mathrm{C}}\simeq L_{z}$ (calculated with the
mean magnetic-field strength at saturation). In 2D, we observed that a high-
energy peak develops due to the “quiet-loop” end state that emerges due to the
reduced dimensionality; here, it appears that the system is simply reaching a
state where particles accumulate around the highest attainable energy
according to the simulation-box size. This evolution is very similar to that
of forced-turbulence kinetic simulations, where particles accumulate, at late
times, around the maximum energy determined by the box size (Zhdankin et al.
2018b). In our runs, the accumulation of particles around the maximum energy
and the shift of the distribution to larger average energies progressively
reduce the extent of the power law. The development of the high-energy peak
also relates to (but is not a direct consequence of) the differential heating
discussed above, with electrons accumulating faster around the highest
attainable energy.

Both the differential heating and the emergence of the high-energy peak can be
ameliorated with larger boxes and/or better scale separation, which imply
higher computational costs; for example, we have conducted preliminary runs
with larger domain sizes, confirming that in such cases the high-energy peak
shifts to higher $\gamma$, proportionally to the extent of the smallest side
(usually $L_{z}$) of the simulation domain. Systematically studying the effect
of box size and scale separation on particle energization in 3D involves high
computational costs, and we thus defer it to future work.

###  5.4 Stresses and angular-momentum transport

Figure 17: Top left: evolution of $\alpha$ and of the volume-averaged Maxwell,
Reynolds, and anisotropic stress in the 3D pair-plasma MRI simulation of
Figure 13 ($\omega_{\rm C}/\Omega_{0}=15$, $\beta_{0}=312$, and box size
$4\times 8\times 2\lambda_{\mathrm{MRI}}^{3}$). Top right: comparison of
$\alpha$ for 2D and 3D simulations employing the same $\omega_{\rm
C}/\Omega_{0}$ and $\beta_{0}$ and different box sizes. Bottom left and right:
comparison of the volume-averaged Maxwell and anisotropic stress for the same
2D and 3D simulations.

To conclude our first study of 3D MRI simulations, we consider the evolution
of the $\alpha$-parameter and of the viscous stresses. Figure 17 (top left)
shows the measured values of $\alpha$ and of the volume-averaged $M_{xy}$,
$R_{xy}$, and $A_{xy}$ for the simulation presented in Figure 13 ($\omega_{\rm
C}/\Omega_{0}=15$, $\beta_{0}=312$, and box size $4\times 8\times
2\lambda_{\mathrm{MRI}}^{3}$). We clearly observe that, for this run, the
anisotropic stress dominates the overall angular-momentum transport; this
differs substantially from 2D cases, where $A_{xy}>M_{xy}$ was verified only
at late times (when the MRI is halting and pressure anisotropy accumulates).
The top-right, bottom-left, and bottom-right panels of Figure 17 show,
respectively, the evolution of $\alpha$ and of the volume-averaged $M_{xy}$
and $A_{xy}$ for the 2D and 3D cases already compared in section 5.1.1. In all
these cases $\omega_{\rm C}/\Omega_{0}=15$ and $\beta_{0}=312$, while the box
size and dimensionality are varied to compare the 2D and 3D dynamics as well
as the 3D small- and large-box evolution. We observe that, during the active
nonlinear phase, $\alpha$ is smaller for 3D runs than in the 2D case; a
further decrease in $\alpha$ occurs from small to large box sizes in 3D runs.
Comparing the most important stress components, we notice that $M_{xy}$
follows the same decreasing trend observed for $\alpha$, while $A_{xy}$
remains essentially unchanged across all runs. It is therefore evident that
the smaller value of $\alpha\simeq 0.1$ achieved in the large-box 3D case
(whereas $\alpha\geq 1$ in the active nonlinear phase of the 2D run) is solely
related to a decrease in the Maxwell stress.

The explanation for this behavior is straightforward: as we have argued in the
previous sections, large 3D boxes drive much more efficient reconnection
dynamics, which in turn results in a magnetic-field saturation at lower values
than both small 3D simulations and arbitrarily large 2D runs. The Maxwell
stress, on average, becomes correspondingly smaller and so does $\alpha$.
Since $A_{xy}$ remains somewhat constant across all these cases, our results
could in principle suggest that the main driving mechanism for angular-
momentum transport in large 3D simulations is represented by pressure
anisotropy (this result was also reported by Hoshino 2015). However, we
believe that the large value of $A_{xy}$ observed here may simply be an
artifact of the reduced scale separation employed in our 3D numerical
experiments. In particular, we measure $\langle A_{xy}\rangle/\langle
p\rangle\simeq 0.1$ in these runs, which is nearly twice as large as the value
found in 2D simulations with large boxes or large scale separation (compare
Figure 17). In such cases, pressure anisotropy is more efficiently limited by
developed turbulence and/or faster mirror modes; with the limited scale
separation and box size employed in our 3D runs, it is not unreasonable to
think that the measured anisotropy is thus exaggerated. 3D Braginskii-MHD and
hybrid simulations (Sharma et al. 2006; Kunz et al. 2016) have indeed shown
that $A_{xy}\sim M_{xy}$ is expected in fluid modeling, but that anisotropic
momentum transport does not typically dominate the overall value of $\alpha$
by the large margin we observe here. To assess the precise role and trend of
anisotropic stresses in 3D fully kinetic simulations, we will conduct large-
box studies with better scale separation in future work.

##  6 Discussion and conclusions

In this work, we have presented a comprehensive exploration of the
magnetorotational instability (MRI) in collisionless pair plasmas via fully
kinetic Particle-in-Cell (PIC) simulations. With a shearing-box setup
implemented in our relativistic PIC code Zeltron (Cerutti et al. 2013), we
have carried out a vast array of 2D runs, exploring an unprecedentedly large
parameter space. In particular, we have conducted very large-scale 2D
simulations with macroscopic-to-microscopic temporal-scale separation up to
$\sim 10$ times larger than previous works and with system size up to twice
larger than the largest simulations presented in literature. In addition, we
have carried out large-scale 3D PIC simulations of the MRI in pair plasmas,
achieving for the first time a global mesoscale dynamics akin to that observed
in MHD works. To study the axisymmetric MRI with PIC, we have resorted to the
established 2D shearing-coordinate approach employed in previous publications
(Riquelme et al. 2012; Inchingolo et al. 2018); for 3D nonaxisymmetric
simulations, we have developed and applied a novel “orbital-advection”
formulation of the shearing box that simplifies pre-existing methods (Hoshino
2013, 2015) and is less complicated to implement numerically.

An important result of our study is the demonstration that 2D (axisymmetric)
fully kinetic simulations of the MRI are qualitatively different from
corresponding 3D kinetic runs. We have pointed out for the first time such
differences, and provided a physical explanation for the behavior of the
kinetic MRI in 2D. We found that axisymmetric PIC simulations are
intrinsically limited in the physical processes they can model, and such a
limitation has important consequences for the overall MRI evolution. In
particular, the lack of reconnection along the out-of-plane (toroidal)
direction results in a lack of saturation of the strongest magnetic-field
component, $B_{y}$, with detrimental consequences for the development and
sustainment of turbulence. However, we also found that many important lessons
can still be learned from 2D runs which carry over to the 3D case: for
example, with our extensive parameter scans in 2D, we have determined the
minimal values of key physical parameters (in terms of the dimensionless scale
separation, system size, etc.) that need to be reached for meaningful studies.
Thus, 2D simulations give us guidance (starting point) for designing 3D runs.
Differences and similarities between the 2D and 3D cases are illustrated in
detail below.

By conducting large-scale 3D simulations of the MRI, for the first time we
have demonstrated that fully kinetic simulations can reproduce the mesoscale
dynamics typically observed in MHD: multi-dimensional, macroscopic channel
flows cyclically form and get disrupted by parasitic instabilities,
particularly via drift-kink modes. The channel disruption results in developed
turbulence, and several channel creation/growth/disruption cycles can be
observed during the nonlinear stage in large-scale 3D runs. This result also
aligns with previous 3D hybrid simulations (Kunz et al. 2016). Conversely,
channels developing in small-scale runs are not efficiently disrupted and no
turbulence develops, in agreement with the sole 3D fully kinetic study
presented in literature so far (Hoshino 2015), which employed a simulation box
of limited size.

In addition to studying the general MRI dynamics, we have explored in detail
several physical processes that develop concurrently with the linear and
nonlinear stage of the main instability. Our results can be summarized as
follows:

  * •

We have extensively described the typical nonlinear evolution of the MRI in
fully kinetic shearing-box simulations, explaining each phase of the evolution
and the differences between the small- and large-box cases, as well as between
the 2D and 3D cases. We have explained how small and large simulations differ
in producing developed turbulent states, and how the choice of physical
parameters impacts the overall evolution and duration of the nonlinear MRI. We
have also identified the role of pressure anisotropy in these simulations,
observing that a large $\Delta p=p_{\perp}-p_{\|}$ can develop and accumulate
throughout the nonlinear stage of the instability; this anisotropy is likely
exaggerated by our choice of modest scale separation (but still much larger
than that employed in previous studies) that necessarily characterizes PIC
simulations. A limited scale separation can result in inefficient mirror
modes, which grow much slower than the MRI, allowing for persistently large
pressure anisotropy. This anisotropy ultimately participates in pushing MRI
modes to wavelengths larger than the box size, halting the MRI dynamics over
long times.

  * •

By exploring a large parameter space, we have assessed the numerical
convergence of key quantities (e.g. the magnetic-field amplification) in 2D
simulations when the physical parameters are varied. We have studied the
effect of the separation between the macroscopic ($\Omega_{0}^{-1}$) and
microscopic ($\omega_{C}^{-1}$) temporal scales, $\omega_{\rm C}/\Omega_{0}$;
of the separation between box size and MRI scales $L/\lambda_{\mathrm{MRI}}$;
of the initial plasma-$\beta$; and of the box aspect ratio. Increasing
$\omega_{\rm C}/\Omega_{0}$ with a fixed box size generally results in a
longer duration of the nonlinear stage and in an increase of the magnetic
energy at saturation. We have shown that, at some large $\omega_{\rm
C}/\Omega_{0}$, the results (in particular, the saturated magnetic energy)
eventually converge; however, this converged nonlinear state is qualitatively
different for small and large boxes, and the precise value of $\omega_{\rm
C}/\Omega_{0}$ resulting in convergence depends on the box size. Similarly,
fixing $\omega_{\rm C}/\Omega_{0}$ and only increasing the box size produces
results that converge for sufficiently large box sizes; in large boxes, a
sustained-turbulence state can develop while ti does not in small simulations.
This partly agrees with previous 2D work by Inchingolo et al. (2018), who
focused their analysis on increasing box sizes while keeping $\omega_{\rm
C}/\Omega_{0}$ fixed. They concluded that convergence in the magnetic-field
amplification is attained, in 2D, for box sizes of at least $8\times
8\lambda_{\mathrm{MRI}}^{2}$, and that large-box simulations can maintain a
volume-averaged $v_{\mathrm{A}}/c<1$ (calculated with the nonrelativistic
expression) throughout the system’s evolution, respecting the underlying
assumptions of the nonrelativistic shearing box. We have instead demonstrated
that 2D simulations, over sufficiently long times, invariably develop an
average $v_{\mathrm{A}}/c>1$, owing to the absence of reconnection in the
$y$-direction. This does not occur, instead, in our 3D runs of sufficiently
large size; we conclude that 3D simulations are of key importance to obtain
physically valid results that respect all the underlying assumptions of the
shearing-box model. Considering the effect of the initial plasma-$\beta$, we
found that the results are practically unchanged in the range $\beta_{0}\simeq
0.1\mbox{--}2{,}500$: independent of the initial $\beta_{0}$, all simulations
end up with plasma beta more or less within the same universal final range, of
order $0.1\mbox{--}1$ for 2D runs. This result holds as long as the unstable
MRI modes can fit in the simulation box, i.e. under the constraint
$L_{z}>\lambda_{\mathrm{MRI}}$. Finally, we have observed that the geometry of
the simulation box can impact the development of the nonlinear MRI stage
dramatically: boxes with aspect ratio $L_{x}/L_{z}>1$ result in developed
turbulence much more easily than corresponding boxes with $L_{x}/L_{z}=1$,
both in 2D and 3D. We have observed that this is caused by additional drift-
kink modes that are impeded in small simulation domains, and that can promote
channel-disruption events in elongated boxes instead.

  * •

We have verified that sustained turbulence develops, in our simulations,
during the nonlinear MRI evolution. The isotropic power spectra of the
poloidal ($B_{\mathrm{pol}}=\sqrt{B_{x}^{2}+B_{z}^{2}}$) and toroidal
($B_{y}$) magnetic field during this phase show the presence of an inertial
range with characteristic power laws, indicative of turbulent activity. Both
in 2D and 3D, the $B_{\mathrm{pol}}$ spectrum features a shallow slope roughly
consistent with $\propto k^{-5/3}$ at large and intermediate spatial scales,
and a spectral break in the vicinity of the average Larmor-radius wavenumber
$k_{\rho}=\rho_{\mathrm{C}}^{-1}$. In the 3D case, at length scales below
$\rho_{\mathrm{C}}$ the $B_{\mathrm{pol}}$ spectrum steepens to a $\propto
k^{-3.2}$ slope and a clear kinetic range is present. The $B_{y}$ spectrum
appears to follow a $\propto k^{-3}$ slope at all scales in 2D; in 3D, $B_{y}$
instead follows a shallower slope compatible with $\propto k^{-2}$ in the
inertial range, and a $k^{-3.2}$ slope in the kinetic range. Our 2D results
are consistent with previous 2D studies (Inchingolo et al. 2018); our 3D
results for $B_{y}$, in the inertial range, are in agreement with MHD and
hybrid-kinetic simulations (Kunz et al. 2016; Walker et al. 2016). However,
the latter comparison is complicated by the difference in the underlying
models (our fully kinetic PIC vs. MHD or hybrid methods). Considering works
focusing on pair plasmas, a similarity can be drawn between our model and that
employed in Zhdankin et al. 2017, 2018a, since both cases consider a forced-
turbulence system evolution. However, in those works a $\propto k^{-4}$ slope
(or steeper) was found for $k>k_{\rho}$, which differs from our shallower
$k^{-3.2}$ result in the kinetic range. This measurement is also hard to
compare against analytic expectations, as no previous studies (to the best of
our knowledge) focused on the specific case of pair-plasma MRI-driven
turbulence we consider. For example, Loureiro & Boldyrev 2018 considered a
low-$\beta$ pair plasma in a tearing-mediated cascade, finding a $\propto
k^{-3}$ slope in the kinetic range; our conditions, however, are those of a
high-beta plasma and we have no basis to claim that our cascade is tearing-
mediated in nature. It is also interesting to note that solar-wind
observations commonly report a spectral slope $\propto
k^{-7/3}\mbox{--}k^{-8/3}$ (not extremely far from our $\propto k^{-3.2}$
result) at sub-ion-Larmor scales, consistent with gyrokinetic calculations
(e.g. Schekochihin et al. 2009). Whether a similarity can be drawn between
this and our case will be verified with future electron-ion simulations.

  * •

Concerning particle energization in 2D and 3D, we found that energy injection
in large runs is akin to that brought about by Alfvénic turbulence. In these
simulations turbulence is well developed, and particle heating proceeds
smoothly, with the system achieving a steady-state balance between turbulent
magnetic-energy dissipation and particle energization. In contrast, when
turbulence is not well developed (e.g. in small boxes), particle energization
occurs mostly during short “bursts” corresponding to large-scale reconnection
events. Moreover, for the first time we have reported that substantial
differential heating can occur between electrons and positrons in MRI
simulations. The physical mechanism behind this phenomenon is represented by
additional (related to the background differential rotation) drift forces that
affect opposite charges differently, causing a symmetry break in the
gyromotion of electrons and positrons in a uniform magnetic field. This effect
is unrealistically large in our simulations due to our choice of (necessarily)
limited scale separation $\omega_{\rm C}/\Omega_{0}$. We have elaborated on
how tuning the values of physical parameters can ameliorate this issue and its
implications for the interpretation of the results.

  * •

The energy distribution functions in our 2D and 3D simulations show the
presence of substantial nonthermal particle acceleration during the nonlinear
MRI stage. A power-law tail with index $\sim 2.5\mbox{--}3$ consistently
develops as the MRI transitions to sustained turbulence, both in 2D and in 3D.
Our results are in good agreement with previous works; however, here we have
demonstrated that in 2D this power-law state is transitory, and that the
overall evolution can produce very different energy distributions over time.
Inchingolo et al. (2018) carried out a similar analysis but focused on the
early nonlinear stage, where a power-law index $\sim 2$ can indeed be
realized; but such a state is clearly still evolving, and changes drastically
later on. At late times, the power-law part of the distributions progressively
disappears, and a high-energy peak develops. This is completely determined by
the artificial 2D end state (with characteristic “magnetic loops”), driven by
the reduced dimensionality of these runs. Riquelme et al. (2012) described a
similar late-time 2D dynamics, but attributed this evolution to the MRI; we
emphasize that such a state is in fact not representative of the response of
particles to the MRI, since the latter has slowed down (or completely stopped)
by the time when the quiescent end state has developed. Conversely, in our 3D
simulations the developed nonthermal features are maintained with the same
power-law index until the end of the run. However, the nonthermal tail is
progressively eroded by the accumulation of particles around the highest
energy attainable (which increases with box size). Hoshino (2015) observed a
similar effect in small-box 3D simulations; this result is also consistent
with forced-turbulence PIC simulations (Zhdankin et al. 2018b) and points to
the need for larger (and more expensive) 3D numerical experiments.

  * •

In our 2D and 3D simulations, we observed that viscous stresses (in
particular, the dominant $xy$-component of the stress tensor) can develop
during the nonlinear MRI stage. This can lead to efficient angular-momentum
transport: an effective (dimensionless) collisionless viscosity
$\alpha=(\langle M_{xy}\rangle+\langle A_{xy}\rangle+\langle
R_{xy}\rangle)/\langle p\rangle$ (Shakura & Sunyaev 1973) arises both in 2D
and 3D runs, mainly due to Maxwell and anisotropic stresses. In 2D, the
Maxwell stress $M_{xy}=-B_{x}B_{y}/(4\pi)$ is likely exaggerated due to a lack
of saturation in $B_{y}$ caused by the reduced dimensionality; in 3D large-
scale runs, where magnetic fields saturate at lower amplitudes, this stress
firmly settles on $\langle M_{xy}\rangle/\langle p\rangle\simeq
0.01\mbox{--}0.1$, consistently with previous hybrid simulations (Kunz et al.
2016). The anisotropic stress $A_{xy}=-(p_{\perp}-p_{\|})B_{x}B_{y}/B^{2}$,
instead, shows a more complex trend: when keeping the same $\omega_{\rm
C}/\Omega_{0}$, the average pressure anisotropy $(p_{\perp}-p_{\|})$ generally
decreases from 2D to 3D and from small to large boxes; but since the average
$B^{2}$ decreases as well, $\langle A_{xy}\rangle/\langle p\rangle$ maintains
roughly the same average value $\sim 0.1$ in 2D and 3D simulations. As a
result, we measure $A_{xy}>M_{xy}$ in our largest 3D runs, which is against
the general expectation $A_{xy}\lesssim M_{xy}$. We believe that this enhanced
anisotropic angular-momentum transport is related to our (necessarily) limited
scale separation and system size; with realistic parameters, efficient mirror
modes would rapidly quench pressure anisotropy, likely decreasing $A_{xy}$.
Hoshino (2015) measured a similarly enhanced angular-momentum transport due to
pressure anisotropy in small-box 3D simulations, which was likely due to the
same mechanisms we describe here.

Although we have addressed several key aspects of the kinetic physics of the
MRI, our study still presents several limitations. The most prominent open
question is whether 3D PIC simulations could produce results that completely
align with hybrid and MHD studies, and how the nonthermal particle
acceleration and angular-momentum transport depend on the box size and choice
of physical parameters. Here, we have conducted the largest (to date) 3D PIC
runs that successfully reproduce the global mesoscale MRI dynamics, but our
scale separation and box size remain rather limited. Larger simulations with
better separation would be needed, i) to completely characterize the
properties of particle energization, e.g. whether large-scale runs, where the
box size does not rapidly interfere with the nonthermal part of the spectrum,
present the same nonthermal features as those observed here; and ii) to
understand the general trend of viscous stresses when the scale separation is
substantially larger, possibly achieving the expected ordering $A_{xy}\lesssim
M_{xy}$ in simulations where pressure anisotropy is well regulated. Moreover,
the issue of substantial differential heating in 3D runs remains to be
addressed, as it can complicate the analysis of particle energization even in
pair-plasma simulations.

Despite these limitations, our work paves the way for several future studies
focused on the exploration of other properties of the MRI. For example, we
have reported here the stark contrast between simulations with and without net
magnetic flux; this is a well-known aspect of the MRI dynamics that has been
extensively studied with MHD shearing-box models (e.g. Gardiner & Stone 2005).
Several works have even pointed out that 3D simulations with and without net
magnetic flux behave substantially differently for boxes with different aspect
ratios (e.g. Bodo et al. 2008; Walker & Boldyrev 2017); moreover, in ideal MHD
it is generally observed that the $\alpha$-parameter representing angular-
momentum transport tends to 0 in the case of zero net magnetic flux, unless
explicit dissipation coefficients are included (e.g. Fromang & Papaloizou
2007; Fromang et al. 2007). In collisionless models (including ours), such
dissipative effects are naturally included (as shown e.g. by Kunz et al.
2016), and the phenomenology of angular-momentum transport can be studied from
first principles. All in all, we believe that a systematic PIC study of the
zero-net-flux case may be particularly important to understand the exact
mechanisms behind the phenomenology of MRI dynamo; we will pursue such a study
in the future. Another important analysis to be carried out in future work
concerns the initial magnetic-field configuration: here, we have limited
ourselves to the purely vertical-field case, but it is well known that
different field geometries (especially when an initial $B_{y}$ is present) can
result in qualitatively different global MRI evolution and turbulence
properties (e.g. Gardiner & Stone 2005; Goedbloed & Keppens 2022). This aspect
has not been considered in kinetic or hybrid works, but it would be important
to study generic field geometries, which are expected in realistic
astrophysical scenarios. Finally, we note that studying the kinetic MRI with
vertical stratification (as opposed to the case without stratification that we
considered) could in principle provide results that are more physically
relevant for modeling thick accretion disks (e.g. Hirabayashi & Hoshino 2017)
and resolve convergence issues that notably affect MHD studies (e.g. Regev &
Umurhan 2008; Bodo et al. 2014).

We conclude by remarking that the next step in our line of work concerns
electron-ion ($m_{i}\gg m_{e}$) simulations. Studying the case of large mass
ratio, and especially the energy partition between particle species, is
fundamental for the correct interpretation of data from current and future
observations targeting RIAFs (EHT Collaboration 2019). In these collisionless
astrophysical environments, electrons and ions are decoupled and a two-
temperature state can originate from turbulence and/or radiation-reaction
dynamics (Arzamasskiy et al. 2019; Zhdankin et al. 2019; Kawazura et al. 2020,
2021; Zhdankin et al. 2021). With our simulations, we will be able to study
the process of temperature decoupling (and possibly the dynamics of radiation)
during the MRI from first principles. Moreover, information from fully kinetic
simulations can be used as input for global fluid simulations of accretion
around compact objects employing subgrid models to account for the missing
microphysics (Ressler et al. 2015; Chael et al. 2018; Scepi et al. 2022).
Given the ever-growing interest in understanding the dynamics of collisionless
accretion disks around SMBHs, we believe that this work has the potential to
substantially impact our insight into the physics of such space environments.

## Acknowledgements

The authors would like to thank Mario Riquelme, Masahiro Hoshino, Giannandrea
Inchingolo, Andrea Mignone, Gianluigi Bodo, Jim Stone, Eliot Quataert, Anatoly
Spitkovsky, Matthew Kunz, Elias Most, Sasha Philippov, Nuno Loureiro, and
Nicolas Scepi for useful discussions and suggestions throughout the
development of this work.

The authors acknowledge support from a NASA ATP grant 80NSSC20K0545.
Additional support was provided by National Science Foundation grants AST
1806084 and AST 1903335 to the University of Colorado. F.B. was also partially
supported by a Junior PostDoctoral Fellowship (grant number 12ZW220N) from
Research Foundation – Flanders (FWO). Support for L.A. was provided by the
Institute for Advanced Study. Support for V.Z. at the Flatiron Institute is
provided by the Simons Foundation.

This work made use of substantial computational resources provided by several
entities. We acknowledge support by the VSC (Flemish Supercomputer Center),
funded by the Research Foundation – Flanders (FWO) and the Flemish Government
– department EWI. An award of computer time was provided by the Innovative and
Novel Computational Impact on Theory and Experiment (INCITE) program, using
resources (namely Theta) of the Argonne Leadership Computing Facility, which
is a DOE Office of Science User Facility supported under contract DE-
AC02-06CH11357. Additional computer time was provided by the Extreme Science
and Engineering Discovery Environment (XSEDE), which is supported by National
Science Foundation grant number ACI-1548562 (Towns et al., 2014), and by the
Frontera computing project (Stanzione et al., 2020) at the Texas Advanced
Computing Center (TACC, www.tacc.utexas.edu). We acknowledge TACC for
providing HPC resources on both Stampede2 and Frontera. Additional resources
(namely Pleiades) supporting this work were provided by the NASA High-End
Computing (HEC) Program through the NASA Advanced Supercomputing (NAS)
Division at Ames Research Center.

##  Appendix A Numerical methods for the kinetic shearing box

###  A.1 Shearing-coordinate simulations in 2D

For two-dimensional simulations, in this work we employ the KSB-SC
formulation, which is preferable since no complicated boundary conditions are
required (the simulation domain is periodic in all directions). Numerically,
the implementation of the KSB-SC system of equations in a PIC code requires
the solution of modified Maxwell’s equations and a peculiar particle-push
step. These modifications are however straightforward: for the electromagnetic
fields, the additional terms on the right-hand side of equations (13) and (14)
are linear and point-wise, and can therefore be treated explicitly with a
standard leapfrog algorithm by interpolating in time and space to maintain
second-order accuracy (see Riquelme et al. 2012). The particle-push step is
similarly formulated as a modification of a standard Boris scheme that
includes here extra (linear) terms. We refer to section A.2.2 for more details
on this modified Boris algorithm.

###  A.2 Orbital-advection simulations in 3D

For three-dimensional simulations, in this work we make use of the newly
developed KSB-OA framework (section 2.3). The numerical implementation of our
method requires solving modified Maxwell’s equations which generally include
implicit terms, as well as a modified Boris push (similar to that required for
2D runs, see section A.1 above). Our numerical approach also requires
shearing-periodic boundary conditions on electromagnetic fields, source terms,
and particles.

####  A.2.1 Solution of Maxwell’s equations

Our version of Maxwell’s equations (20)–(21) involves nonlocal, coupled terms
on the right-hand side; this complicates the numerical solution, since a
simple leapfrog scheme cannot be applied to explicitly push the
electromagnetic fields in time. To avoid costly matrix-inversion operations
(see e.g. Bacchini et al. 2019), in our first numerical implementation we
split the solution in two explicit steps. First, all terms on the right-hand
side of Maxwell’s equations that are not related to advection along $y$ are
discretized with a central-difference scheme in time and space. The resulting
system of discrete equations (excluding the advective terms) is iterated upon
with a simple predictor-corrector step, which is applied 3 times121212This
approach is essentially an iterative Crank-Nicolson scheme. It is well known
that this scheme is not asymptotically stable for an infinite number of
iterations (see e.g. Leiler & Rezzolla 2006); the scheme is, however, stable
for 3 iterations, which is our choice for all simulations. This represents a
good compromise between stability and accuracy in the solution.. At the end of
each predictor-corrector iteration, we can include the advection of $E$ and
$B$ along $y$ in several ways. In our numerical experiments, we have observed
that a simple central-difference scheme applied to these terms inevitably
results in instabilities developing at the shearing boundaries; we could
therefore opt for a more stable advection scheme among the many available in
literature. However, as a first approach we choose an even simpler strategy:
to model the advection of field quantities along $y$, we simply shift each
component of $E$ and $B$ in space by a distance $v_{\mathrm{s},y}(x)\Delta t$
(with appropriate signs). The shifted arrays are interpolated at the resulting
spatial locations, which approximates the advection step. This operation
introduces some numerical diffusion, but we observe that this is in fact
beneficial as it suppresses numerical instabilities at the shearing
boundaries. Our approach represents a cheap, less sophisticated alternative to
more classical schemes; while this leaves ample grounds for improvements, it
only requires simple spatial interpolation to model advection, and we observe
that this choice produces acceptable results for this first study. In future
works, we will implement more robust algorithms to accurately capture the
field advection and minimize the numerical diffusion introduced by our
approach. Finally, Gauss’s law (equation (23)) is enforced, in our approach,
with a Langdon-Marder correction (Marder 1987; Langdon 1991); the resulting
(small) errors on charge conservation are thus bounded in time.

####  A.2.2 Modified Boris push and particle-position update

To solve the momentum equation (40), we proceed as follows: first we define
$\mbox{\boldmath{$\epsilon$}}=q\Delta t\mbox{\boldmath{$E$}}/(2m)$,
$\mbox{\boldmath{$\tau$}}=q\Delta t\mbox{\boldmath{$B$}}/(2mc)$; we then
discretize the equation of motion between two consecutive time steps as

$\displaystyle\mbox{\boldmath{$u$}}^{n+1/2}-\mbox{\boldmath{$u$}}^{n-1/2}$
$\displaystyle=2\mbox{\boldmath{$\epsilon$}}+\frac{\mbox{\boldmath{$u$}}^{n+1/2}+\mbox{\boldmath{$u$}}^{n-1/2}}{\bar{\gamma}}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$\tau$}}$
(A1) $\displaystyle+\Delta
t\left(\mbox{\boldmath{$u$}}^{n+1/2}+\mbox{\boldmath{$u$}}^{n-1/2}\right)\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$\Omega$}}_{0}+\frac{s\Omega_{0}\Delta
t}{2}(u^{n+1/2}_{x}+u^{n-1/2}_{x})\hat{\mbox{\boldmath{$e$}}}_{y}$
$\displaystyle+s\Omega_{0}x^{n}\left(\frac{\mbox{\boldmath{$u$}}^{n+1/2}+\mbox{\boldmath{$u$}}^{n-1/2}}{\bar{\gamma}}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\epsilon$}}\right)\hat{\mbox{\boldmath{$e$}}}_{y},$

where $\bar{\gamma}$ is the Lorentz factor defined at an intermediate time
between $n+1/2$ and $n-1/2$. Now we recast this equation as

$\mbox{\boldmath{$u$}}^{+}-\mbox{\boldmath{$u$}}^{+}\,\mbox{\boldmath{$\times$}}\,\widetilde{\mbox{\boldmath{$\tau$}}}-\alpha_{1}u^{+}_{x}\hat{\mbox{\boldmath{$e$}}}_{y}-\alpha_{2}(\mbox{\boldmath{$u$}}^{+}\,\mbox{\boldmath{$\cdot$}}\,\widetilde{\mbox{\boldmath{$\epsilon$}}})\hat{\mbox{\boldmath{$e$}}}_{y}=\mbox{\boldmath{$u$}}^{-}+\mbox{\boldmath{$u$}}^{+}\,\mbox{\boldmath{$\times$}}\,\widetilde{\mbox{\boldmath{$\tau$}}}+\alpha_{1}u^{-}_{x}\hat{\mbox{\boldmath{$e$}}}_{y}+\alpha_{2}(\mbox{\boldmath{$u$}}^{-}\,\mbox{\boldmath{$\cdot$}}\,\widetilde{\mbox{\boldmath{$\epsilon$}}})\hat{\mbox{\boldmath{$e$}}}_{y},$
(A2)

where we have defined $\mbox{\boldmath{$u$}}^{\pm}=\mbox{\boldmath{$u$}}^{n\pm
1/2}\mp\mbox{\boldmath{$\epsilon$}}$,
$\widetilde{\mbox{\boldmath{$\tau$}}}=\mbox{\boldmath{$\tau$}}/\bar{\gamma}+\Delta
t\mbox{\boldmath{$\Omega$}}_{0}$,
$\widetilde{\mbox{\boldmath{$\epsilon$}}}=\mbox{\boldmath{$\epsilon$}}/\bar{\gamma}$,
$\alpha_{1}=s\Omega_{0}\Delta t/2$, $\alpha_{2}=s\Omega_{0}x^{n}$, and
$\mbox{\boldmath{$u$}}^{+}$ is the unknown to be solved for. Assuming that
$\bar{\gamma}$ is known (see next paragraph), equation (A2) can be solved
explicitly by recasting it as

$\textbf{A}\mbox{\boldmath{$u$}}^{+}=\mbox{\boldmath{$C$}},$ (A3)

where the right-hand side is the known term
$\mbox{\boldmath{$C$}}=\mbox{\boldmath{$u$}}^{-}+\mbox{\boldmath{$u$}}^{-}\,\mbox{\boldmath{$\times$}}\,\widetilde{\mbox{\boldmath{$\tau$}}}+\alpha_{1}u^{-}_{x}\hat{\mbox{\boldmath{$e$}}}_{y}+\alpha_{2}(\mbox{\boldmath{$u$}}^{-}\,\mbox{\boldmath{$\cdot$}}\,\widetilde{\mbox{\boldmath{$\epsilon$}}})\hat{\mbox{\boldmath{$e$}}}_{y}$,
and

$\textbf{A}=\begin{bmatrix}1&-\widetilde{\tau}_{z}&\widetilde{\tau}_{y}\\\
\widetilde{\tau}_{z}-\alpha_{1}-\alpha_{2}\widetilde{\epsilon}_{x}&1-\alpha_{2}\widetilde{\epsilon}_{y}&-\widetilde{\tau}_{x}-\alpha_{2}\widetilde{\epsilon}_{z}\\\
-\widetilde{\tau}_{y}&\widetilde{\tau}_{x}&1\end{bmatrix}.$ (A4)

The solution is then straightforward,

$\mbox{\boldmath{$u$}}^{+}=\textbf{A}^{-1}\mbox{\boldmath{$C$}},$ (A5)

with the inverse matrix

$\textbf{A}^{-1}=\frac{1}{D}\begin{bmatrix}1+\widetilde{\tau}_{x}^{2}+\alpha_{2}(\widetilde{\tau}_{x}\widetilde{\epsilon}_{z}-\widetilde{\epsilon}_{y})&\widetilde{\tau}_{z}+\widetilde{\tau}_{x}\widetilde{\tau}_{y}&-\widetilde{\tau}_{y}+\widetilde{\tau}_{x}\widetilde{\tau}_{z}+\alpha_{2}(\widetilde{\tau}_{y}\widetilde{\epsilon}_{y}+\widetilde{\tau}_{z}\widetilde{\epsilon}_{z})\\\
-\widetilde{\tau}_{z}+\widetilde{\tau}_{x}\widetilde{\tau}_{y}+\alpha_{1}+\alpha_{2}(\widetilde{\tau}_{y}\widetilde{\epsilon}_{z}+\widetilde{\epsilon}_{x})&1+\widetilde{\tau}_{y}^{2}&\widetilde{\tau}_{x}+\widetilde{\tau}_{y}\widetilde{\tau}_{z}-\alpha_{1}\widetilde{\tau}_{y}+\alpha_{2}(-\widetilde{\tau}_{y}\widetilde{\epsilon}_{x}+\widetilde{\epsilon}_{z})\\\
\widetilde{\tau}_{y}+\widetilde{\tau}_{x}\widetilde{\tau}_{z}-\alpha_{1}\widetilde{\tau}_{x}-\alpha_{2}(\widetilde{\tau}_{x}\widetilde{\epsilon}_{x}+\widetilde{\tau}_{y}\widetilde{\epsilon}_{y})&-\widetilde{\tau}_{x}+\widetilde{\tau}_{y}\widetilde{\tau}_{z}&1+\widetilde{\tau}_{z}^{2}-\alpha_{1}\widetilde{\tau}_{z}-\alpha_{2}(\widetilde{\tau}_{z}\widetilde{\epsilon}_{x}+\widetilde{\epsilon}_{y})\end{bmatrix},$
(A6)

where

$D=1+\widetilde{\tau}^{2}-\alpha_{1}(\widetilde{\tau}_{z}+\widetilde{\tau}_{x}\widetilde{\tau}_{y})-\alpha_{2}[\widetilde{\tau}_{y}(\widetilde{\mbox{\boldmath{$\tau$}}}\,\mbox{\boldmath{$\cdot$}}\,\widetilde{\mbox{\boldmath{$\epsilon$}}})+\widetilde{\epsilon}_{y}+\widetilde{\epsilon}_{x}\widetilde{\tau}_{z}-\widetilde{\epsilon}_{z}\widetilde{\tau}_{x}].$
(A7)

After $\mbox{\boldmath{$u$}}^{+}$ has been calculated, the particle velocity
at the next time step is simply

$\mbox{\boldmath{$u$}}^{n+1/2}=\mbox{\boldmath{$u$}}^{+}+\mbox{\boldmath{$\epsilon$}}.$
(A8)

The definition of $\bar{\gamma}$ is not unique and it is in fact the
peculiarity that generally distinguishes relativistic particle pushers (see
e.g. Ripperda et al. 2018). For simplicity, we can define the mid-point
Lorentz factor as

$\bar{\gamma}=\sqrt{1+\left[(u^{-})^{2}+2\alpha_{1}u^{-}_{x}u^{-}_{y}+2\alpha_{2}(\mbox{\boldmath{$u$}}^{-}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\epsilon$}})u^{-}_{y}/\bar{\gamma}\right]/c^{2}},$
(A9)

which is a choice similar to that used in the standard Boris pusher. This
definition is such that in the absence of electric fields, energy is conserved
exactly. We can solve the resulting depressed cubic equation for
$\bar{\gamma}$ with any available algorithm, e.g. the one described by Tiruneh
(2020): define $X=a_{1}/3$, $Y=-a_{2}/2$, where

$a_{1}=-1-\left[(u^{-})^{2}+2\alpha_{1}u^{-}_{x}u^{-}_{y}\right]/c^{2},$ (A10)
$a_{2}=-2\alpha_{2}(\mbox{\boldmath{$u$}}^{-}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\epsilon$}})u^{-}_{y}/c^{2},$
(A11)

such that $\bar{\gamma}^{3}+a_{1}\bar{\gamma}+a_{2}=0$. The discriminant
$Z=X^{3}+Y^{2}$ indicates that there is only one real root if $Z>0$,

$\bar{\gamma}=\sqrt[3]{Y+\sqrt{Z}}+\sqrt[3]{Y-\sqrt{Z}},$ (A12)

and three real roots (of which only one positive and $>1$) if $Z<0$,

$\bar{\gamma}_{j}=2\sqrt{-X}\cos\frac{\phi+2(j-1)\pi}{3},$ (A13)

where $j=1,2,3$ and $\phi=\cos^{-1}\left(Y/\sqrt{-X^{3}}\right)$.

When the new particle velocity is known, the position push can be carried out
in two steps. The position is first updated with the motion with respect to
the background, i.e.

$\mbox{\boldmath{$x$}}^{\prime}=\mbox{\boldmath{$x$}}^{n}+\Delta
t\frac{\mbox{\boldmath{$u$}}^{n+1/2}}{\gamma^{n+1/2}},$ (A14)

after which the current can be collected with the usual Villasenor-Buneman
charge deposition (Villasenor & Buneman 1992) using the displacement
$\mbox{\boldmath{$x$}}^{\prime}-\mbox{\boldmath{$x$}}^{n}$. Then, the position
is updated with the shearing-velocity part,

$\mbox{\boldmath{$x$}}^{n+1}=\mbox{\boldmath{$x$}}^{\prime}-s\Omega_{0}\frac{x^{\prime}+x^{n}}{2}\Delta
t\hat{\mbox{\boldmath{$e$}}}_{y},$ (A15)

which only modifies the $y$-coordinate, using the $x$-position at the midpoint
in time (note that $x^{n+1}\equiv x^{\prime}$). Depositing the current using
$\mbox{\boldmath{$x$}}^{\prime}-\mbox{\boldmath{$x$}}^{n}$ instead of
$\mbox{\boldmath{$x$}}^{n+1}-\mbox{\boldmath{$x$}}^{n}$ is essential to
recover $\mbox{\boldmath{$J$}}^{\prime}$ (in this paradigm, the current due to
velocity fluctuations with respect to the background) instead of $J$.

As noted in section 2.3, if particle velocities are relativistic the particle
position update above may result in superluminal motion, since it is not
guaranteed that
$(\mbox{\boldmath{$u$}}/\gamma+\mbox{\boldmath{$v$}}_{\mathrm{s}})^{2}<c^{2}$.
As an alternative strategy, we can employ the proper relativistic velocity
boost to update the particle position with a speed that is ensured to be
subluminal. Applying a central-difference scheme to equations (47)–(49), the
position update is then formulated as

$x^{n+1}=x^{n}+\Delta
t\frac{u_{x}^{n+1/2}}{[\gamma^{n+1/2}+u_{y}^{n+1/2}v_{\mathrm{s},y}^{n+1/2}]\gamma_{\mathrm{s}}^{n+1/2}},$
(A16) $y^{n+1}=y^{n}+\Delta
t\frac{u_{y}^{n+1/2}}{\gamma^{n+1/2}+u_{y}^{n+1/2}v_{\mathrm{s},y}^{n+1/2}}+\Delta
t\frac{\gamma^{n+1/2}v_{\mathrm{s},y}^{n+1/2}}{\gamma^{n+1/2}+u_{y}^{n+1/2}v_{\mathrm{s},y}^{n+1/2}},$
(A17) $z^{n+1}=z^{n}+\Delta
t\frac{u_{z}^{n+1/2}}{[\gamma^{n+1/2}+u_{y}^{n+1/2}v_{\mathrm{s},y}^{n+1/2}]\gamma_{\mathrm{s}}^{n+1/2}},$
(A18)

where $v_{\mathrm{s},y}^{n+1/2}=-s\Omega_{0}x^{n+1/2}$,
$\gamma_{\mathrm{s}}^{n+1/2}=1/\sqrt{1-(v_{\mathrm{s},y}^{n+1/2}/c)^{2}}$, and
$x^{n+1/2}=(x^{n+1}+x^{n})/2$. These equations are less trivial to employ,
especially since equation (A16) is now nonlinearly implicit in the unknown
$x^{n+1}$. Fortunately, the nonlinearity is weak and there are no numerical
pathologies on the right-hand side (e.g., the denominator is always $>1$). The
position update can then be carried out in the following steps:

  1. 1.

Solve the nonlinear equation (A16) for $x^{\prime}\equiv x^{n+1}$, obtaining
$v_{\mathrm{s},y}^{n+1/2},\gamma_{\mathrm{s}}^{n+1/2}$ in the process. We have
found that 3-4 fixed-point iterations with $x^{n}$ as initial guess typically
suffice in converging below an absolute tolerance of $10^{-10}$.

  2. 2.

Partially update the $y$-coordinate as

$y^{\prime}=y^{n}+\Delta
t\frac{u_{y}^{n+1/2}}{\gamma^{n+1/2}+u_{y}^{n+1/2}v_{\mathrm{s},y}^{n+1/2}}.$
(A19)

  3. 3.

Calculate $z^{\prime}\equiv z^{n+1}$ with equation (A18).

  4. 4.

Using $\mbox{\boldmath{$x$}}^{\prime}$, deposit the comoving-frame current
$\mbox{\boldmath{$J$}}^{\prime}$.

  5. 5.

Complete the update of the $y$-coordinate with the shearing-velocity part,

$y^{n+1}=y^{\prime}+\Delta
t\frac{\gamma^{n+1/2}v_{\mathrm{s},y}^{n+1/2}}{\gamma^{n+1/2}+u_{y}^{n+1/2}v_{\mathrm{s},y}^{n+1/2}}.$
(A20)

##  Appendix B Differential heating in kinetic shearing-box simulations of
pair plasmas

Kinetic shearing-box simulations are characterized by additional drift forces
in the particle equations of motion, related to the orbital motion of the
accretion disk. These forces are not charge-agnostic — that is, opposite
charges of equal mass experience difference forces acting upon their motion.
This causes a symmetry break in the dynamics of electrons and positrons; as a
result, in all our simulations (both 2D and 3D) we observe that the two
species in a pair plasma do not experience the same energization when solving
the shearing-box equations. This appears surprising if compared to standard
(non-shearing-box) PIC simulations of pair plasmas, where the two species
behave essentially in the same way, and the opposite sign in the charge only
serves to maintain global quasi-neutrality. Here, we prove that the
differential energization we observe is not a numerical artifact, but a
physical effect that characterizes shearing-box PIC simulations.

Figure 18: Left: trajectories of test particles (an electron and a positron)
moving in the $xy$-plane when a uniform and constant magnetic field
$\mbox{\boldmath{$B$}}=(0,0,B_{0})$ is present. The particle motion is
computed by solving the equations employed in the shearing-box paradigm.
Right: spatial distribution of the radial and toroidal magnetic field (top
row) after the formation of channel flows in a small-box, 2D MRI simulation of
size $2\times 2\lambda_{\mathrm{MRI}}^{2}$ (with $\omega_{\rm
C}/\Omega_{0}=120$, $\beta_{0}=624$). The corresponding spatial distribution
of particle drift velocity along $z$ (bottom row) shows that electrons are
pulled toward current sheets, and positrons are pushed away.

This effect stems from the kinetic shearing-box paradigm, which involves
additional forces in the particle equations of motion. In the KSB-SC system,
the momentum equation (18) contains the force term

$\mbox{\boldmath{$F$}}_{\mathrm{SB}}=2m\mbox{\boldmath{$u$}}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$\Omega$}}_{0}+ms\Omega_{0}u_{x}\hat{\mbox{\boldmath{$e$}}}_{y}=\begin{bmatrix}2m\Omega_{0}u_{y}\\\
m(s-2)\Omega_{0}u_{x}\\\ 0\end{bmatrix},$ (B1)

which is related to Coriolis and centrifugal forces in the comoving frame
where all quantities are measured (the same applies to the KSB-OA system).
This additional force causes a drift of the particle guiding center which is
absent in typical non-shearing-box kinetic simulations. In the nonrelativistic
limit, this drift has the usual expression,

$\mbox{\boldmath{$v$}}_{\mathrm{SB}}=\frac{\langle\mbox{\boldmath{$F$}}_{\mathrm{SB}}\rangle\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}}}{qB^{2}}c=\frac{mc}{qB^{2}}\begin{bmatrix}(s-2)\Omega_{0}\langle
u_{x}\rangle B_{z}\\\ -2\Omega_{0}\langle u_{y}\rangle B_{z}\\\
2\Omega_{0}\langle u_{y}\rangle B_{y}-(s-2)\Omega_{0}\langle u_{x}\rangle
B_{x}\end{bmatrix},$ (B2)

where $\langle...\rangle$ indicates averaging over the particle gyration.
Since the sign of the drift depends on the particle charge, in general
particles of equal mass but opposite charge will experience opposite drifts in
all spatial directions.

We now show that the $\mbox{\boldmath{$v$}}_{\mathrm{SB}}$ drift is directly
related to the dynamics leading to the measured differential energization.
First, we note that the initial conditions in our MRI simulations,
$\mbox{\boldmath{$B$}}=(0,0,B_{z,0})$, lead to a drift motion exclusively in
the $xy$-plane. Such motion is precisely the well-known epicyclic motion of
particles in orbit around massive objects (e.g. Balbus & Hawley 1998), which
is here modified by the presence of a vertical magnetic field. In Figure 18
(left panel), we show a portion of the trajectory of one electron and one
positron (with $m=|q|=1$) moving according to equation (18) in the presence of
a weak magnetic field $B_{z,0}$ such that $\omega_{\rm C}/\Omega_{0}=1$, with
initial velocity $\mbox{\boldmath{$v$}}=(0,0.1,0)$. We observe that the usual
particle gyromotion is here modified by a drift of the guiding center along
$\pm y$ depending on the particle charge $q=\mp 1$. Moreover, the gyromotion
is in itself different for the two particles, with the electron attaining
higher gyrofrequency and smaller gyroradius (and vice versa for the positron).
Although the guiding-center drift is purely in the toroidal direction in this
case, this differential gyromotion ought to manifest in our 2D ($x$-$z$)
simulations as well, given that the particle oscillation along $x$ occurs
differently for opposite charges. This difference in particle trajectories
implies that charge neutrality cannot be maintained at all times, and a net
current must arise.

While the differential gyromotion is present, at least in the beginning, in
all our simulations, it is arguably of minor importance in the overall system
evolution. Indeed, in section 4.5 we have shown that significant differential
heating primarily manifests during and after the linear stage, when the
magnetic-field geometry differs substantially from the initial (pre-
instability) phase. We find that relevant differential heating arises due to
the channel-flow configuration of $B_{x}$ and $B_{y}$ during the later stages:
$B_{z}$ here is heavily subdominant (see section 4.1), and the main drift
component of $\mbox{\boldmath{$v$}}_{\mathrm{SB}}$ is along $z$. In Figure 18
(right panels, top row) we plot the spatial distribution of $B_{x}$ and
$B_{y}$ (normalized over the initial field $B_{0}$) for a 2D, small-box
($2\times 2\lambda_{\mathrm{MRI}}^{2}$) simulation with parameters
$\omega_{\rm C}/\Omega_{0}=120$, $\beta_{0}=624$ at $t=3.3P_{0}$, when channel
flows are present and current sheets have formed at the channel interfaces. In
the same Figure (right panels, bottom row), we also plot the spatial
distribution of $v_{\mathrm{SB,z}}$ for electrons and positrons taken at the
same moment in time. We observe that, in the vicinity of current sheets, both
species are drifting along $z$; however, the particular magnetic-field
configuration causes this drift to be directed _toward_ current sheets in the
case of electrons, and _away_ from them in the case of positrons. In a
standard current-sheet (e.g. Harris-like) magnetic configuration, particles
drift towards current sheets due to the usual $E$ $\times$ $B$ motion
irrespective of their charge; this global behavior is still verified in our
case, but the addition of the $\mbox{\boldmath{$v$}}_{\mathrm{SB}}$ drift now
implies that charges of opposite sign drift toward current sheets at different
rates.

Since particles experience most of the energy gain via reconnection in current
sheets, this differential drift motion can explain the large difference in
particle energy observed between the two species. Electrons, which flow toward
current sheets faster, are capable of gaining more energy over shorter times;
this precisely corresponds to the faster electron energization rates reported
in section 4.5. Eventually, all available magnetic energy is distributed, via
reconnection, unequally among the pair-plasma species, with electrons
experiencing a larger energy gain.

The differential heating we report here for the first time (to the best of our
knowledge) is a genuine physical effect which, however, is arguably of very
minimal importance in realistic scenarios. Indeed, it can immediately be
observed that the drift speed $v_{\mathrm{SB}}$ is inversely proportional to
$\omega_{\rm C}/\Omega_{0}$ (and indeed, for the $\omega_{\rm
C}/\Omega_{0}=120$ case shown in Figure 18, $v_{\mathrm{SB}}/c\ll 1$). In
realistic environments where $\omega_{\rm C}/\Omega_{0}\sim 10^{7}$, this
drift would be completely negligible and therefore no differential heating
would be observed. In our simulations (as well as in any other presented in
literature), however, we are bound to employ much smaller values of
$\omega_{\rm C}/\Omega_{0}$ due to otherwise excessive computational costs.
Therefore, when measuring particle heating in simulations, this effect must be
taken into account; particularly in the context of future electron-ion
numerical experiments, where analyzing particle heating is of main interest
for observations of accreting compact objects. In such simulations, it will be
of major importance to quantify the “spurious” differential heating described
here, in order to correctly evaluate the energy partition between particle
species.

We conclude by reporting the differential heating measured in our 2D
simulation campaign presented in section 4. This is shown in Figure 19, where
the difference in electron and positron total energy $mc^{2}(\gamma-1)$ is
normalized over the instantaneous total energy of the system, in order to
measure the relative importance of differential heating with respect to the
global energetics. We observe that larger $\omega_{\rm C}/\Omega_{0}$ values
and larger boxes can ameliorate this issue substantially, keeping the
differential heating below 10%. While it is clear that a larger $\omega_{\rm
C}/\Omega_{0}$ implies smaller drift speeds and therefore smaller differences
in the rate at which particles flow towards current sheets, justifying the
positive effect of larger box sizes is less straightforward. Qualitatively, we
can consider that smaller boxes tend to create persistent large-scale channel
flows that are less likely to be destroyed by unstable modes and result in
developed turbulence (as illustrated in section 4.5); this large-scale,
ordered distribution of magnetic fields is the most favorable for particles to
drift undisturbed toward currents sheets at different rates. In large boxes
where turbulence develops, particles are less likely to experience such long-
lived, undisturbed motion, being instead continuously scattered across the
turbulent structures. Differential heating can still accumulate over long
times, reaching however much smaller levels.

Figure 19: Differential heating in 2D pair-plasma MRI simulations, measured as
the difference in kinetic energy between electrons and positrons normalized
over the instantaneous total energy of the system. Simulations with small box
size and variable $\omega_{\rm C}/\Omega_{0}$ (left), variable box size and
fixed $\omega_{\rm C}/\Omega_{0}$ (center), and large box size and variable
$\omega_{\rm C}/\Omega_{0}$ (right) are compared.

##  Appendix C Energy conservation in the shearing box

Figure 20: Energy conservation in a 2D pair-plasma MRI simulation with
parameters $\omega_{\rm C}/\Omega_{0}=30$, $\beta_{0}=624$, and box size
$8\times 8\lambda_{\mathrm{MRI}}^{2}$. Left: evolution of the total
electromagnetic and kinetic energy compared with the modified electromagnetic
and kinetic energy, where shearing-box source terms have been subtracted. The
modified total energy $\bar{\mathrm{E}}_{\mathrm{tot}}$ is a conserved
quantity in the absence of numerical errors. Right: evolution during the run
of the error on the conserved $\bar{\mathrm{E}}_{\mathrm{tot}}$ relative to
the instantaneous total energy $\mathrm{E}_{\mathrm{tot}}$ (the latter
includes source terms and is not conserved).

Considering for simplicity the 2D shearing box in shearing coordinates,
equations KSB-SC manifestly include source terms affecting the field and
particle evolution. These source terms translate into energy sources that are
not present in standard (non-shearing-box) PIC simulations. Evaluating energy
conservation in the presence of these sources is not straightforward.

To obtain a measure of energy conservation, we compute the energy sources
directly at runtime. For the fields, the energy source terms are obtained by
dotting Faraday’s and Ampère’s laws with $B$ and $E$, respectively:

$\mbox{\boldmath{$B$}}\,\mbox{\boldmath{$\cdot$}}\,\frac{\partial\mbox{\boldmath{$B$}}}{\partial
t}=\frac{1}{2}\frac{\partial B^{2}}{\partial t}=\mbox{(standard energy
terms)}-s\Omega_{0}B_{x}B_{y},$ (C1)
$\mbox{\boldmath{$E$}}\,\mbox{\boldmath{$\cdot$}}\,\frac{\partial\mbox{\boldmath{$E$}}}{\partial
t}=\frac{1}{2}\frac{\partial E^{2}}{\partial t}=\mbox{(standard energy
terms)}-s\Omega_{0}E_{x}E_{y},$ (C2)

where in addition to standard energy-exchange terms we now have extra source
terms stemming from the shearing-box paradigm. For the particles, we dot the
momentum equation with $\mbox{\boldmath{$u$}}/\gamma$,

$\frac{\mbox{\boldmath{$u$}}}{\gamma}\,\mbox{\boldmath{$\cdot$}}\,\frac{\textrm{d}\mbox{\boldmath{$u$}}}{\textrm{d}t}=\frac{\textrm{d}\gamma}{\textrm{d}t}=\mbox{(standard
energy terms)}+s\Omega_{0}\frac{u_{x}u_{y}}{\gamma},$ (C3)

which produces one extra energy source. All terms deriving from the shearing-
box framework are thus proportional to $\Omega_{0}$.

To quantify the importance of the shearing-box sources, we compute each source
term in the whole simulation box during a representative 2D pair-plasma MRI
simulation with parameters $\omega_{\rm C}/\Omega_{0}=30$, $\beta_{0}=624$,
and box size $8\times 8\lambda_{\mathrm{MRI}}^{2}$. The source terms are then
integrated in time, to obtain the net energy variation associated with these
terms. Subtracting the net contribution of each term from the electromagnetic
and particle energy produces a new measure of energy that only accounts for
standard energy-exchange channels (i.e. without external energy sources
related to the shearing box). In simulations, the total energies including the
source terms are given by

$B^{2}(t)=\sum_{g}B^{2}_{g}(t)\Delta V_{g},$ (C4)
$E^{2}(t)=\sum_{g}E^{2}_{g}(t)\Delta V_{g},$ (C5)

for the electromagnetic energy at each grid cell $g$ (with cell volume $\Delta
V_{g}$), and

$K(t)=\sum_{p}m_{p}c^{2}(\gamma_{p}(t)-1),$ (C6)

for the kinetic energy of all particles $p$. By subtracting the shearing-box
terms, the _modified_ field energies (i.e. without shearing-box sources) are

$\frac{1}{2}\bar{B}^{2}(t)\equiv\frac{1}{2}B^{2}(t)-\int_{0}^{t}dt\sum_{g}\Delta
V_{g}\left(-s\Omega_{0}B_{x,g}(t)B_{y,g}(t)\right),$ (C7)
$\frac{1}{2}\bar{E}^{2}(t)\equiv\frac{1}{2}E^{2}(t)-\int_{0}^{t}dt\sum_{g}\Delta
V_{g}\left(-s\Omega_{0}E_{x,g}(t)E_{y,g}(t)\right),$ (C8)

and the modified particle energy is

$\bar{K}(t)\equiv
K(t)-\int_{0}^{t}dt\sum_{p}s\Omega_{0}\frac{u_{x,p}(t)u_{y,p}(t)}{\gamma_{p}(t)}.$
(C9)

With these definitions, the modified energies are such that

$\partial_{t}\bar{\mathrm{E}}_{\mathrm{tot}}=\partial_{t}\left(\frac{\bar{B}^{2}+\bar{E}^{2}}{8\pi}+\bar{K}\right)=0,$
(C10)

i.e. the sum of the modified electromagnetic and kinetic energy integrated
over the simulation box remains constant in time (as in standard PIC
simulations), in the absence of numerical errors. Note that conversely, the
total energy including sources
$\mathrm{E}_{\mathrm{tot}}\equiv(B^{2}+E^{2})/(8\pi)+K$ is not constant: the
sources continuously inject energy into the system.

Figure 20 (left panel) shows the evolution in time of the modified
electromagnetic and kinetic energy (integrated over the box) as well as of the
sum of the two. The plot shows that the modified electromagnetic energy
becomes _negative_ : this is because the external energy input for the
electromagnetic fields is much larger than the initial energy, and subtracting
this energy source results in observing exclusively the standard process of
magnetic energy dissipation, e.g. via reconnection. Because of the external
source, more magnetic energy than initially available can be dissipated, and
the modified electromagnetic energy can fall below 0. The decrease in
electromagnetic energy corresponds almost exactly to an increase in particle
energy via dissipation processes. The total (electromagnetic plus kinetic)
modified energy is much smaller than its individual components, owing to the
fact that the modified electromagnetic energy becomes negative.

The right panel of Figure 20 shows a measure of energy conservation calculated
as
$|\bar{\mathrm{E}}_{\mathrm{tot}}(t)-\bar{\mathrm{E}}_{\mathrm{tot}}(0)|/\mathrm{E}_{\mathrm{tot}}(t)$.
This quantity is such that the numerator expresses the deviation from exact
energy conservation in terms of the absolute error on the conserved quantity
$\bar{\mathrm{E}}_{\mathrm{tot}}(t)$; normalizing over the instantaneous total
energy (including sources) serves to provide a measure of how important the
numerical errors are compared to the actual amount of energy in the system at
a certain time. We see that the error is at most of order $10^{-2}$, i.e.
numerical energy dissipation is negligible throughout a typical run.

## References

  * Alt & Kunz (2019) Alt, A., & Kunz, M. W. 2019, JPP, 85, 175850102
  * Arendt Jr. (1998) Arendt Jr., P. N. 1998, arxiv:astro-ph/9801194
  * Arzamasskiy et al. (2019) Arzamasskiy, L., Kunz, M. W., Chandran, B. D. G., & Quataert, E. 2019, ApJ, 879, 53
  * Bacchini et al. (2019) Bacchini, F., Amaya, J., & Lapenta, G. 2019, J. Phys: Conf. Ser., 1225, 012011
  * Bai & Stone (2013) Bai, X., & Stone, J. M. 2013, ApJ, 767, 1
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1
  * Ball et al. (2018) Ball, D., Sironi, L., & Feryal, O. 2018, ApJ, 862, 80
  * Bodo et al. (2014) Bodo, G., Cattaneo, F., Mignone, A., & Rossi, P. 2014, ApJL, 787, L13
  * Bodo et al. (2008) Bodo, G., Mignone, A., Cattaneo, F., Rossi, P., & Ferrari, A. 2008, A&A, 487,1
  * Cerri et al. (2021) Cerri, S. S., Arzamasskiy, L., & Kunz, M. W. 2021, ApJ, 916, 120
  * Cerutti et al. (2013) Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147
  * Chael et al. (2018) Chael, A., Rowan, M., Narayan, R., Johnson, M., & Sironi, L. 2018, MNRAS, 478, 4
  * Chandran et al. (2010) Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E., & Germaschewski, K. 2010, ApJ, 720, 503
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, PNAS, 46, 2
  * Chandrasekhar et al. (1958) Chandrasekhar, S., Kaufman, A. N., & Watson, K. M. 1958, RSPSA, 245, 435
  * Chen & Palmadesso (1984) Chen, J., & Palmadesso, P. 1984, The Physics of Fluids, 27, 1198
  * Comisso & Sironi (2018) Comisso, L., & Sironi, L. 2018, PRL, 121, 25
  * Comisso & Sironi (2019) —. 2019, ApJ, 886, 122
  * Dexter et al. (2021) Dexter, J., Scepi, N., & Begelman, M. C. 2021, ApJL, 919, L20
  * EHT Collaboration (2019) EHT Collaboration. 2019, ApJL, 875, L5
  * EHT Collaboration (2022) —. 2022, ApJL, 930, L16
  * Fromang & Papaloizou (2007) Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113
  * Fromang et al. (2007) Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123
  * Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, AIP Conference Proceedings, 784, 475
  * Goedbloed & Keppens (2022) Goedbloed, H., & Keppens, R. 2022, arXiv:2201.11551
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 2
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763
  * Goldreich & Tremaine (1978) Goldreich, P., & Tremaine, S. 1978, ApJ, 222, 850
  * Goodman & Xu (1994) Goodman, J., & Xu, G. 1994, ApJ, 432, 213
  * Gressel & Ziegler (2007) Gressel, O., & Ziegler, U. 2007, CPC, 176, 652
  * Guo et al. (2014) Guo, F., Li, H., Daughton, W., & Liu, Y. 2014, PRL, 113, 155005
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  * Hill (1878) Hill, G. W. 1878, American Journal of Mathematics, 1, 5
  * Hillas (1984) Hillas, A. M. 1984, Annu. Rev. Astron. Astr., 22, 425
  * Hirabayashi & Hoshino (2017) Hirabayashi, K., & Hoshino, M. 2017, ApJ, 842, 36
  * Hirai et al. (2018) Hirai, K., Katoh, Y., Terada, N., & Kawai, S. 2018, ApJ, 853, 174
  * Hoppock et al. (2018) Hoppock, I. W., Chandran, B. D. G., Klein, K. G., Mallet, A., & Verscharen, D. 2018, Journal of Plasma Physics, 84, 905840615
  * Hoshino (2013) Hoshino, M. 2013, ApJ, 773, 118
  * Hoshino (2015) —. 2015, PRL, 114, 061101
  * Hoshino (2020) —. 2020, ApJ, 900, 66
  * Howes (2010) Howes, G. G. 2010, MNRAS, 409, L104
  * Inchingolo et al. (2018) Inchingolo, G., Grismayer, T., Loureiro, N. F., Fonseca, R. A., & Silva, L. O. 2018, ApJ, 859,149
  * Isenberg (2004) Isenberg, P. A. 2004, JGRA, 109, A03101
  * Isenberg & Vasquez (2007) Isenberg, P. A., & Vasquez, B. J. 2007, ApJ, 668, 546
  * Kaisig (1989) Kaisig, M. 1989, A&A, 218, 1
  * Kawazura et al. (2019) Kawazura, Y., Barnes, M., & Schekochihin, A. A. 2019, PNAS, 116
  * Kawazura et al. (2021) Kawazura, Y., Schekochihin, A. A., Barnes, M., Dorland, W., & Balbus, S. A. 2021, arXiv:2110.12434
  * Kawazura et al. (2020) Kawazura, Y., Schekochihin, A. A., Barnes, M., et al. 2020, PRX, 10, 041050
  * Kennel & Engelmann (1966) Kennel, C. F., & Engelmann, F. 1966, Physics of Fluids, 9, 2377
  * Kunz et al. (2018) Kunz, M. W., Abel, I. G., Klein, K. G., & Schekochihin, A. A. 2018, Journal of Plasma Physics, 84, 715840201
  * Kunz et al. (2015) Kunz, M. W., Schekochihin, A. A., Chen, C. H. K., Abel, I. G., & Cowley, S. C. 2015, , 81, 325810501
  * Kunz et al. (2014a) Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014a, PRL, 112, 205003
  * Kunz et al. (2014b) Kunz, M. W., Stone, J. M., & Bai, X. 2014b, JCP, 259, 154
  * Kunz et al. (2016) Kunz, M. W., Stone, J. M., & Quataert, E. 2016, PRL, 117, 235101
  * Langdon (1991) Langdon, A. B. 1991, CPC, 70, 447
  * Latter et al. (2009) Latter, H. N., Lesaffre, P., & Balbus, S. A. 2009, MNRAS, 394, 715
  * Leiler & Rezzolla (2006) Leiler, G., & Rezzolla, L. 2006, PRD, 73, 044001
  * Ley et al. (2019) Ley, F., Riquelme, M., Sironi, L., Verscharen, D., & Sandoval, A. 2019, ApJ, 880, 100
  * Loureiro & Boldyrev (2018) Loureiro, N. F., & Boldyrev, S. 2018, ApJL, 866, L14
  * Loureiro & Boldyrev (2020) —. 2020, ApJ, 890, 55
  * Mallet et al. (2017) Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2017, MNRAS, 468, 4
  * Marder (1987) Marder, B. 1987, JCP, 68, 48
  * Masada & Sano (2008) Masada, Y., & Sano, T. 2008, ApJ, 689, 1234
  * Matsumoto & Tajima (1995) Matsumoto, R., & Tajima, T. 1995, ApJ, 445, 767
  * McChesney et al. (1987) McChesney, J. M., Stern, R. A., & Bellan, P. M. 1987, PRL, 59, 1436
  * Narayan et al. (1987) Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1
  * Parker (1958) Parker, E. N. 1958, PR, 109, 1874
  * Pessah & Goodman (2009) Pessah, M. E., & Goodman, J. 2009, ApJ, 698, L72
  * Pritchett et al. (1996) Pritchett, P. L., Coroniti, F. V., & Decyk, V. K. 1996, JGR, 101, A12
  * Quataert et al. (2002) Quataert, E., Dorland, W., & Hammett, G. W. 2002, ApJ, 577, 524
  * Regev & Umurhan (2008) Regev, O., & Umurhan, O. M. 2008, A&A, 481, 21
  * Ressler et al. (2015) Ressler, S. M., Tchekhovskoy, A., Quataert, E., Chandra, M., & Gammie, C. F. 2015, MNRAS, 454, 1848
  * Ripperda et al. (2018) Ripperda, B., Bacchini, F., Teunissen, J., et al. 2018, ApJS, 235, 21
  * Riquelme et al. (2012) Riquelme, M. A., Quataert, E., Sharma, P., & Spitkovsky, A. 2012, ApJ, 755, 50
  * Riquelme et al. (2015) Riquelme, M. A., Quataert, E., & Verscharen, D. 2015, ApJ, 800, 27
  * Riquelme et al. (2016) —. 2016, ApJ, 824, 123
  * Rosenbluth (1956) Rosenbluth, M. N. 1956, LANL Report LA-2030
  * Ryan et al. (2017) Ryan, B. M., Ressler, S. M., Dolence, J. C., et al. 2017, ApJL, 844, L24
  * Sagdeev & Shafranov (1960) Sagdeev, R. Z., & Shafranov, V. D. 1960, Soviet Physics JETP, 181
  * Sano & Inutsuka (2001) Sano, T., & Inutsuka, S. 2001, ApJ, 561, L179
  * Sano & Stone (2002) Sano, T., & Stone, J. M. 2002, ApJ, 577, 534
  * Scepi et al. (2022) Scepi, N., Dexter, J., & Begelman, M. C. 2022, MNRAS, 511, 3536
  * Schekochihin et al. (2009) Schekochihin, A. A., Cowley, S. C., Dorland, W., et al. 2009, ApJS, 182, 310
  * Schiff (1939) Schiff, L. I. 1939, PNAS, 25, 391
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Shapiro & Shevchenko (1964) Shapiro, V. D., & Shevchenko, V. I. 1964, JETP, 18, 1109
  * Sharma et al. (2006) Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. M. 2006, ApJ, 637, 952
  * Sironi & Spitkovsky (2014) Sironi, L., & Spitkovsky, A. 2014, ApJL, 783, 1
  * Squire et al. (2022) Squire, J., Meyrand, R., Kunz, M. W., et al. 2022, Nature Astronomy
  * Stanzione et al. (2020) Stanzione, D., West, J., Evans, R. T., et al. 2020, Frontera: The Evolution of Leadership Computing at the National Science Foundation (New York, NY, USA: Association for Computing Machinery), 106–111.  https://doi.org/10.1145/3311790.3396656
  * Stone & Gardiner (2010) Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142
  * Stone et al. (1996) Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656
  * Tiruneh (2020) Tiruneh, A. T. 2020, arXiv:2002.06976
  * Towns et al. (2014) Towns, J., Cockerill, T., Dahan, M., et al. 2014, Comput. Sci. Eng., 16, 62
  * Vedenov & Sagdeev (1958) Vedenov, A. A., & Sagdeev, R. Z. 1958, SPhD, 3, 278
  * Velikhov (1959) Velikhov, E. P. 1959, JETP, 36,1938
  * Villasenor & Buneman (1992) Villasenor, J., & Buneman, O. 1992, CPC, 69, 306
  * Walker & Boldyrev (2017) Walker, J., & Boldyrev, S. 2017, MNRAS, 470, 2653
  * Walker et al. (2016) Walker, J., Lesur, G., & Boldyrev, S. 2016, MNRAS, 457, L39
  * Werner & Uzdensky (2021) Werner, G. R., & Uzdensky, D. A. 2021, JPP, 87, 905870613
  * Werner et al. (2018) Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4
  * Werner et al. (2016) Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJL, 816, 1
  * Winarto & Kunz (2022) Winarto, H. W., & Kunz, M. W. 2022, JPP, 88, 905880210
  * Wong et al. (2020) Wong, K., Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2020, ApJL, 893, L7
  * Zenitani & Hoshino (2007) Zenitani, S., & Hoshino, M. 2007, ApJ, 670, 702
  * Zenitani & Hoshino (2008) Zenitani, S., & Hoshino, M. 2008, ApJ, 677, 1
  * Zhdankin et al. (2021) Zhdankin, V., Uzdensky, D. A., & Kunz, M. W. 2021, ApJ, 908, 71
  * Zhdankin et al. (2015) Zhdankin, V., Uzdensky, D. A., & PBoldyrev, S. 2015, PRL, 114, 065002
  * Zhdankin et al. (2013) Zhdankin, V., Uzdensky, D. A., Perez, J. C., & Boldyrev, S. 2013, ApJ, 771, 124
  * Zhdankin et al. (2018a) Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2018a, MNRAS, 474, 2
  * Zhdankin et al. (2018b) —. 2018b, ApJL, 867, L18
  * Zhdankin et al. (2019) —. 2019, PRL, 122, 055101
  * Zhdankin et al. (2017) Zhdankin, V., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2017, PRL, 118, 055103
  * Zweibel & Yamada (2009) Zweibel, E. G., & Yamada, M. 2009, Annu. Rev. Astron. Astr., 47, 1