# Coupling Fluid Plasma and Kinetic Neutral Models using Correlated Monte
Carlo Methods

G.J. Parker, M.V. Umansky, B.D. Dudson  Stanford University, Department of
Mathematics, Stanford CA 94305, USA  Lawrence Livermore National Lab,
Livermore, CA 94550, USA

###### Abstract

While boundary plasmas in present day tokamaks generally fall in a fluid
regime, neutral species near the boundary often require kinetic models due to
long mean-free-paths compared to characteristic spatial scales in the region.
Monte-Carlo (MC) methods provide a complete, high-fidelity approach to solving
kinetic models, and must be coupled to fluid plasma models to simulate the
full plasma-neutrals system. The statistical nature of MC methods, however,
prevents convergence of coupled fluid-kinetic simulations to an exact self-
consistent steady-state. Moreover, this forces the use of explicit methods
that can suffer from numerical errors and require huge computational
resources.

Correlated Monte-Carlo (CMC) methods are expected to alleviate these issues,
but have historically enjoyed only mixed success. Here, a fully implicit
method for coupled plasma-neutral systems is demonstrated in 1D using the
UEDGE plasma code and a homemade CMC code. In particular, it is shown that
ensuring the CMC method is a differentiable function of the background plasma
is sufficient to employ a Jacobian-Free Newton-Krylov solver for implicit time
steps. The convergence of the implicit coupling method is explored and
compared with explicit coupling and uncorrelated methods. It is shown that
ensuring differentiability by controlling random seeds in the MC is sufficient
to achieve convergence, and that the use of implicit time-stepping methods has
the potential for improved stability and runtimes over explicit coupling
methods.

###### keywords:

fusion, boundary plasma, numerical algorithms, Monte-Carlo methods

##  1 Introduction

High fidelity modeling for coupled systems of plasma and neutral species is
essential in the design and study of divertor regions in fusion devices, where
strong interactions occur [1, 2]. Fluid models provide accurate approximations
for the plasma behavior, and robust approaches using finite-volume methods
have now been employed with great success for several decades [3]. In cases
where the mean-free paths of neutral particles in edge plasmas are long,
neutral species exhibit non fluid-like behavior and are more accurately
modeled using the kinetic Boltzman equation [4, 5]. Monte-Carlo (MC) methods
[1, 6, 7, 8] are by far the most developed among the several known approaches
[9, 10] to solving the kinetic equation.

The stochastic nature of Monte-Carlo methods however, results in unfavorable
convergence properties when coupling to fluid plasma models. Simulations
typically converge to stochastic equilibria, for which robust convergence
criteria may be difficult to determine, and which are known to exhibit bias
compared to true solutions [1, 11, 12]. More importantly, the stochastic
nature of Monte-Carlo methods typically forces the use of explicit coupling
methods. By standard numerical stability theory [13, 14], such explicit
methods can result in large numerical errors and require excessively small
timesteps to achieve desired accuracy. As a result, resolving long timescale
(1s) effects to necessary accuracies using explicit methods can require
prohibitively long runtimes or heavy computational resources. Implicit methods
enjoy much more desirable stability properties and, in theory, have the
potential to be more accurate by an exponential factor for the same size
timestep, but are impossible when coupling to standard Monte-Carlo methods.

It has been known to experts for many years that correlated Monte-Carlo (CMC)
methods have the potential to address these problems, but attempts to realize
this have had only mixed success. To the authors’ knowledge, correlated
methods have only been successfully demonstrated in one spatial dimension, and
only using explicit coupling for which the advantages of correlation were
shown to be negligible [11]. The purpose of this study is to demonstrate a
fully implicit coupling scheme in one spatial dimension with robust stability
and convergence properties employing correlated methods.

Section 2 explains the coupled model studied in detail, and reviews the
difference between explicit and implicit coupling methods. Section 3 offers a
new perspective on Correlated Monte-Carlo methods that is useful in the
development of implicit coupling algorithms. Section 4 introduces the UEDGE
plasma code, a homemade Monte-Carlo code using the JULIA programming language
and the methods of coupling. Section 5 displays the results of coupled
simulations, and Sections 6 and 7 provide the discussion and conclusion
respectively.

##  2 Numerical Methods

The model investigated couples a set of fluid-plasma evolution equations to a
chosen equation modeling neutral species with two basic interaction: (i)
ionization, and (ii) charge exchange. Let $p(x)$ and $n(x)$ denote the plasma
and neutral state vectors, so that e.g. $p=(n_{i},n_{e},T_{i},T_{e},u_{i})$ is
a vector-valued function combining the ion and electron densities,
temperatures, and velocities. The coupled model leads to a system of equations
with the schematic form

$\displaystyle\partial_{t}p+P(p)$ $\displaystyle=$ $\displaystyle S_{1}(p,n)$
(1) $\displaystyle N(n)$ $\displaystyle=$ $\displaystyle S_{2}(p)$ (2)

wherein $P,N$ are the plasma and neutral evolution equations respectively, and
$S_{i}$ comprise the interaction terms. We make the additional assumption that
the neutral species’ relaxation time is on a significantly faster timescale
than the plasma timescale – in this situation it is only necessary to consider
the time-independent version of $N(n)$. It is additionally assumed that
recombination is negligible, and that the source of neutrals from the wall is
independent of the neutral distribution, i.e. $S_{2}$ is independent of $n$.

While we focus here on two specific plasma models encapsulated by a simple
non-linear diffusion equation and the UEDGE model [15, 16], we emphasize that
the techniques described are equally applicable to any chosen pairs of plasma
and neutral models, or more generally to any chosen pair of continuum (e.g.
finite volume) and particle-based (e.g. Monte Carlo) simulations. In
particular, including additional physics effects such as recombination of ions
or additional neutral species should constitute no fundamental change.

###  2.1 Explicit vs. Implicit Coupling

Suppose that known techniques may be applied to time evolve the plasma for a
given fixed neutral distribution, and vice versa. That is, assume that one has
known numerical algorithms implementing the time evolution operator $U_{\Delta
t}(p,S)$ for plasma with timestep $\Delta t$ and fixed source term $S$, and
$N^{-1}$ the inverse operator of the time-independent neutral equation. The
system (1–2) may then be coupled using either explicit or implicit methods.
Note that this terminology refers to the method used to couple the system,
while the individual solvers for $U_{\Delta t},N^{-1}$ may themselves include
an internal choice of an explicit or implicit method.

####  2.1.1 Explicit Coupling

Given a coupled state $(p_{j},n_{j})$ at time $t_{j}=j\Delta t$, an explicit
method advances the state by setting

$\displaystyle p_{j+1}$ $\displaystyle=$ $\displaystyle U_{\Delta t}(p_{j}\ ,\
S_{1}(p_{j},n_{j}))$ (3) $\displaystyle n_{j+1}$ $\displaystyle=$
$\displaystyle N^{-1}(S_{2}(p_{j+1})).$ (4)

The method is explicit in the sense that the states at timestep $j+1$ are
defined entirely using quantities calculated in the previous steps. The above
demonstrates a standard first-order splitting method, but higher order or
Strang splitting methods would be equally applicable. The above time-step is
iterated until a sufficient equilibrium is reached.

Explicit methods have several key advantages and disadvantages. Advantages
include that they are relatively easy to implement and parallelize, and can be
implemented treating the existing individual solvers or codebases for
$U_{\Delta t},N^{-1}$ as “black boxes”. The one main disadvantage is that
explicit coupling methods display the same quite severe numerical instability
as explicit time evolution for a single equation: small errors may compound
exponentially as they propagate forward in time. To guarantee a desired
accuracy on a time interval $[0,T]$, the time-step must be taken to be of size
$O(T^{-1}e^{-\lambda T})$, where $\lambda$ is the largest eigenvalue of the
(linearization of) the time-independent equation [13, Ch. 7]. The second
disadvantage is that explicit coupling schemes rely on the assumption that the
interaction terms are relatively weak, and so for strongly coupled systems
convergence is not guaranteed. Despite these significant shortcomings,
explicit methods are currently the only method available for coupling kinetic
models using Monte-Carlo methods. Contrasting the above with implicit methods
makes the reason for this clear.

####  2.1.2 Implicit Coupling

A first-order implicit method advances the state by $\Delta t$ by taking
$(p_{j+1},n_{j+1})$ as the unique solution of the non-linear system:

$\displaystyle(I+\Delta tP)p_{j+1}+\Delta tS_{1}(p_{j+1},n_{j+1})$
$\displaystyle=$ $\displaystyle p_{j}$ (5) $\displaystyle n_{j+1}$
$\displaystyle=$ $\displaystyle N^{-1}(S_{2}(p_{j+1}),$ (6)

where $P$ now denotes the time-independent plasma equation. Substituting (6)
into (5) results in an equation purely for $p_{j+1}$ defining the evolution of
the plasma. For $\Delta t$ sufficiently small, the equation is well-
conditioned and may be solved using a suitable implementation of Newton’s
method.

Implicit methods enjoy far better stability properties than explicit ones, and
in practice the time-step may usually be taken as large as the Newton solver
will allow without impacting accuracy. The seemingly innocuous disadvantage of
implicit methods is that they require the functions $N,P,S_{i}$ to be
differentiable functions of the plasma and neutral states. This is the
requirement that precludes the use of Monte-Carlo methods (which are a priori
not even constant for repeated evaluation at the same input). From a
mathematical perspective, the key result of this study is to demonstrate that
correlated Monte Carlo methods can be made sufficiently differentiable as a
function of parameters to apply Newton solvers, and thus implicit methods.

##  3 A New View of Correlated MC

When the model $N(n)$ for the neutral species is the kinetic Boltzman
equation, Monte-Carlo methods are the most widely-used technique for finding
the equilibrium distribution. The standard viewpoint of a Monte-Carlo method
(see e.g. [6, 7, 8, 1]) is the following: particle testflights are run using
(pseudo) random numbers to generate particle trajectories, and each testflight
gives a sample from the equilibrium particle distribution. Thus with
sufficiently many flights, the tallied output converges to the solution of
$N(n)=S_{2}(p)$.

A correlated Monte-Carlo method does the same, but fixes the random numbers
generated for each individual testflight to make the output reproducible. The
result is a sample from the same equilibrium distribution, and the
distribution is now obtained by averaging over increasingly large numbers of
possibilities for the fixed random seeds.

###  3.1 Correlated MC Constructs the Inverse Operator

When considering a Monte-Carlo method as a particular step in a larger
algorithm as in (5–6), it is advantageous adopt a more abstract viewpoint.
Rather than working at the level of equations as explained above, let us
consider the situation at the level of operators. Suppose that the neutral
state is described by a vector $n(x)\in\mathbb{R}^{m}$ on a discretized
spatial grid of size $m$. Thus the Boltzman equation $N_{p}(n)=N(n)-S_{2}(p)$
may be considered as a matrix $N_{p}\in M(\mathbb{R}^{m},\mathbb{R}^{m})$ in
the space of operators on the state space. In our case, this operator is
linear, but the same framework applies equally well to a space of nonlinear
operators including neutral-neutral interactions.

The existence of a physical equilibrium distribution dictates that the
equation should be always be solvable for reasonable plasma backgrounds $p$,
i.e. that $N^{-1}_{p}\in M(\mathbb{R}^{m},\mathbb{R}^{m})$ exists. One may
combine all the random choices required for Monte-Carlo flights into a joint
probability distribution $\mathbb{X}$. Then, at an operator level, a
correlated Monte-Carlo method is a (parameterized) operator-valued probability
distribution

$\displaystyle P\times\mathbb{X}$ $\displaystyle\longrightarrow$
$\displaystyle M(\mathbb{R}^{m},\mathbb{R}^{m})$ (7) $\displaystyle(p,\xi)$
$\displaystyle\mapsto$ $\displaystyle N_{p,\xi}^{-1}.$ (8)

That associates to any pair of a background plasma state $p$ and joint draw of
pseudorandom numbers $\xi$ a particular approximate inverse to $N_{p}$. The
true inverse is then obtained as $N_{p}^{-1}=\int
N^{-1}_{p,\xi}d\mathbb{X}(\xi)$, i.e. averaging over increasingly large
samples of the joint probability distribution $\mathbb{X}$ gives better
approximations to the true inverse. For algorithms that require multiple
applications of the operator $N^{-1}_{p}$ (e.g. explicit or implicit
coupling), it can be advantageous to fix a particular choice and carry out the
algorithm to completion, rather than generate a new inverse each time
$N_{p}^{-1}$ is called within the algorithm.

####  3.1.1 Smooth Dependence of MC

The space $M(\mathbb{R}^{m},\mathbb{R}^{m})$ where the approximate inverse
operators are valued is a vector space (more generally a manifold for non-
linear operators), as is the space of plasma parameters $p$. It therefore
makes sense to ask whether, for fixed $\xi\in\mathbb{X}$, the operator inverse
$N_{p,\xi}^{-1}$ is a continuous, differentiable, or smooth function of the
plasma state $p$. Note that without correlation, continuity certainly fails —
even for fixed $p$ each call to an uncorrelated MC method produces a different
result. If differentiability is correctly imposed for fixed $\xi$, this
function may be called repeatedly as a step in larger algorithms relying on
the linearized approximation of a function, in particular various version of
Newton’s method.

In practice, there are at least two challenges in making a MC output (i.e.
$N_{p,\xi}^{-1}$) depend differentiably on $p$. These include:

  1. (i)

De-synchronization of trajectories

  2. (ii)

Discreteness of estimators

If one naively fixes the random seeds, particle flights can de-synchronize for
infinitesimal changes in background parameters. To elaborate, the cut-off
determining whether a particular collision is a charge-exchange or ionization
is a discrete transition; thus a tiny variation in background parameters
causing a fixed random seed to land on the other side of this cut-off may
cause the particle flight to proceed for more steps, thereby shifting the
seeds for subsequent flights. Even if this shift is remedied, the estimator
used to tally the output will change by a discrete transition when the
particle flight elongates, which creates a point of discontinuity with respect
to the parameters.

To fix the first of these issues, it is necessary to identify the joint
probability distribution $\mathbb{X}$ and have a reliable method of generating
joint draws so that these shifts do not occur, as well as imposing continuity
of the tallies. In practice, the former means pre-generating random numbers
for all possible branchings of the particle trajectory before the flight
begins. To fix the second, one must use more advanced estimators with better
continuity properties. In this study, the attenuated absorption method [6] was
used with a track length estimator [17, Eq. 2.35]. If one prefers a
collisional estimator, then finite-size particle methods provide continuity
properties. In general, it is advantageous to identify any source of
discreteness within the MC simulation and replace this by a continuous
quantity.

Figure 1: The differentiability of the correlated Monte-Carlo method (blue)
compared to the discontinuous uncorrelated method (orange). The derivative
displayed is the partial derivative at a particular volume cell with respect
to varying the plasma background $p$ at equilibrium in the direction of a
gaussian perturbation $dp$.

When smoothness is imposed correctly, the derivative

$\frac{\partial N_{p,\xi}^{-1}}{\partial p}\approx\frac{N^{-1}_{p+\Delta
p,\xi}-N^{-1}_{p,\xi}}{\Delta p}$

exists, and can be calculated through finite-difference approximation schemes
as we use here, or potentially by auto-differentiation methods.

##  4 Simulation Setup

The system of equations (1–2) is discretized on a 1-dimensional grid, assuming
both rotational and poloidal symmetry. The grid spacing decreases toward the
domain wall to better resolve edge effects, and upstream boundary conditions
are imposed on the inner boundary.

###  4.1 Simulation Codes

Two plasma models are implemented: first, a simplified model denoted $P_{0}$
for troubleshooting. This model assumes constant velocity, and temperatures,
and that the ion densities satisfies a second-order non-linear diffusion
model. The second model, denoted $P_{1}$ is given by the UEDGE code, described
in detail in [15, 16]. The code implements finite-volume methods for a non-
linear system of coupled fluid and transport equations in the variables
$p=(n_{i},T_{i},T_{e},u_{i})$ of ion density, temperature, and velocity, and
electron temperature $t_{e}$.

The Correlated Monte-Carlo methods are implemented using a simplified,
homemade code written in JULIA. This code includes 1-dimensional capabilities
following large codebases such as DEGAS2 [6, 17]. Testflights are run using
the attenuated absorption method, and tallied with a track length estimator.

###  4.2 Coupling Methods

Although the UEDGE code includes a sophisticated and optimized JFNK solver for
time-stepping, we here opt for a more accessible setup using an external
solver. Thus the above codes are used simply to call their respective
functions $P_{i}$ and $N_{p,\xi}^{-1}$. The coupling schemes (3–4) and (5–6)
are implemented in a Python overhead using SciPy’s JFNK and finite-difference
methods. The workflow of the implicit solver is depicted in Figure 2 below.

Figure 2: The code workflow for fully implicit coupling. The derivative of the
Monte-Carlo method is called repeatedly to generate the Krylov subspace for
the Jacobian during Newton iteration.

###  4.3 Simulation Parameters

Simulations are run beginning from a uniform ion density of $10^{19}$
particles/m2, with uniform starting temperatures $10$ eV, and velocity
$10^{4}$ m/s toward the divertor. The a fixed value $R$ is chosen for the wall
recycling coefficient (generally $R=1$ is taken). Dirichlet boundary
conditions are imposed on the upstream plasma boundary, and Robbin boundary
conditions are imposed at the plate. Neutrals are assumed to have perfect
reflection at walls.

##  5 Simulation results

###  5.1 Correlated vs. Uncorrelated Explicit Methods

The convergence properties of the coupled system (1–2) was investigated using
the explicit coupling scheme (3–4) for the UEDGE plasma model with both
correlated and uncorrelated Monte-Carlo Methods. The rate of convergence is
displayed in Figure 3 on a logarithmic scale as a function of time.

Figure 3:  (Top) Comparison of residual difference from converged solution for
correlated vs uncorrelated explicit coupling on a log scale as a function of
time (Bottom) The average of the stochastic equilibrium at two selected volume
elements compared to the converged correlated solution for a variety of fixed
random seeds.

With uncorrelated MC, the coupled system reaches a stochastic equilibrium
after only a few orders of magnitude of convergence (Figure 3a). The magnitude
of the stochastic oscillations and the amount of convergence displayed is
proportional to $\sqrt{M}$ for $M$ the number of flights. In contrast, the
correlated method produces convergence to the machine floating-point precision
(or set tolerance of the plasma step), regardless of the number of flights.
For the uncorrelated method, true equilibria are obtained by averaging the
stochastic oscillations, whereas for correlated methods by averaging over
different seeds (Figure 3b). These averages are known to each also include a
bias arising from the non-linearity of the plasma model, but this is a second-
order (i.e. $O(M)$) effect and becomes negligible for large particle numbers
[11].

###  5.2 Explicit vs. Implicit Correlated Methods

To begin, implicit methods were tested using a simplified plasma model given
described by the non-linear diffusion equation

$\frac{\partial p}{\partial t}-\nabla(D(p)\nabla p)=0$

where $D(p)=cp^{n}$ is the diffusion constant. The converged equilibria for
the explicit and implicit coupling schemes (3–4) and (5–6) are compared below
in Figure 4.

Figure 4:  (Top) Time evolution of plasma density at a single cell $x=3.8$ for
implicit and explicit coupling with various timesteps $dt$. (Bottom) residual
difference from converged solution as a function of machine runtime for both
explicit and implicit methods plotted on log-log scale.

The implicit coupling method is found to obtain an equilibrium which is
independent of the time-step used up to machine accuracy, while the explicit
coupling method produces an equilibrium which depends on the time-step, and
approaches that of the implicit method as $dt\to 0$ (Figure 4a).

Figure 4b compares the runtimes of the two methods to attain a certain
accuracy. The explicit coupling method is found to require exponentially more
computation to reach the same levels of accuracy as the implicit method.

###  5.3 Implicit Coupling to UEDGE

The fully implicit coupling scheme was also tested for the UEDGE plasma model
in one spatial dimension. Convergence to machine floating-point accuracy was
achieved for a variety of time steps, up to approximately $\sim 0.1$ ms. The
results displayed in Figure 4 are reproducible using UEDGE, though the trends
persist only for several orders of magnitude smaller than the starting
residual, after which long time-scale effects in the UEDGE model dominate. The
exploration of the capabilities of this fully implicit model using UEDGE will
be the subject of future work.

##  6 Discussion

The use of correlated Monte-Carlo methods in coupled Plasma-neutral systems
was investigated. Previous results [11] demonstrating convergence to machine
precision were reproduced using the UEDGE plasma code and a homemade Monte-
Carlo code.

Implicit methods are demonstrated to have significant advantage over explicit
coupling schemes. From a naive perspective, it is not so clear that the
implicit method should be advantageous; indeed, to do a single implicit time-
step requires calling the Monte-Carlo code $O(10^{3})$ times to differentiate,
generate the Krylov subspace, and perform Newton iteration. In contrast each
time-step of the explicit scheme requires only one invocation of the Monte-
Carlo solver. The advantage of the implicit scheme appears when considering
the number of calls to achieve a particular accuracy or a particular
simulation time. Standard numerical theory shows explicit methods can result
in numerical error proportional to $dt\cdot Te^{\lambda T}$ where $dt$ is the
timestep, $T$ is the time to convergence, and $\lambda$ is the (maximum in $t$
of the) largest eigenvalue of the Jacobian, requiring the number of timesteps
to increase exponentially as a function of desired accuracy. In contrast,
implicit methods have accuracy independent of stepsize provided stability
conditions are met, allowing timesteps to be taken several orders of magnitude
larger to model effects on longer timescales without significantly increasing
runtime. Testing both coupling schemes confirms that the fully implicit
coupling scheme indeed realizes these improvements, despite the apparent
difficulty in performing a single time-step.

The use of fully implicit methods presents several avenues for further study.
The robust stability and convergence properties of these methods make them
ideal candidates for generalizing to study the convergence of two and three-
dimensional coupled systems using correlated methods, which have eluded
previous efforts. It is also an interesting prospect to explore whether the
same differentiability properties of the Monte-Carlo output may be obtained
using existing functions in established codebases such as DEGAS2 [8] and
EIRENE [7]. Finally, there are many avenues for optimizing and improving the
performance of the method. In particular, one interesting prospect is to
investigate the use of auto-differentiation methods to generate the Jacobian
of the Monte-Carlo output.

##  7 Conclusions

The methods investigated here display significant potential advantages over
the current coupling schemes of fluid plasma and kinetic neutral models.
First, a shift in viewpoint of Monte-Carlo methods is proposed, to one in
which the code is considered as a “black box” producing an approximate inverse
operator, which allows the use and analysis of more sophisticated (in
particular implicit) coupling methods. Then, by imposing smoothness of the
correlated Monte-Carlo output with respect to background plasma parameters, we
show that the Monte-Carlo code may be differentiated using finite-difference
schemes making it viable for using in solving algorithms based on Newton’s
method, in particular implicit time steps using Krylov based methods.
Converged equilibria using fully implicit coupled simulations are obtained for
a variety of parameters. Comparisons of runtimes demonstrate that implicit
methods enjoy an exponential speedup in machine time over explicit methods to
obtain the same accuracy. Moreover, the equilibria obtained are shown to be
independent of the timestep to machine accuracy, allowing the use of long
timesteps to investigate effects on different timescales without significantly
increasing runtime.

##  8 Acknowledgments

G.P. is supported by an NSF Mathematical Sciences Postdoctoral Research
Fellowship (Award No. 2303102). Performed in part by LLNL under Contract DE-
AC52-07NA27344. G.P. wishes to acknowledge Andreas Holm and George Wilkie for
helpful discussions.

## References

  * Krasheninnikov  _et al._ [2020] S. Krasheninnikov, A. Smolyakov, and A. Kukushkin,  _On the edge of magnetic fusion devices_ (Springer, 2020).
  * Militello [2022] F. Militello,  _Boundary Plasma Physics: An Accessible Guide to Transport, Detachment, and Divertor Design_ , Vol. 123 (Springer Nature, 2022).
  * Wiesen  _et al._ [2015] S. Wiesen, D. Reiter, V. Kotov, M. Baelmans, W. Dekeyser, A. Kukushkin, S. Lisgo, R. Pitts, V. Rozhansky, G. Saibene, _et al._ , Journal of nuclear materials 463, 480 (2015).
  * Van Uytven  _et al._ [2022] W. Van Uytven, W. Dekeyser, M. Blommaert, S. Carli, and M. Baelmans, Nuclear Fusion 62, 086023 (2022).
  * Blommaert  _et al._ [2019] M. Blommaert, N. Horsten, P. Börner, and W. Dekeyser, Nuclear Materials and Energy 19, 28 (2019).
  * Heifetz [1986] D. Heifetz, in  _Physics of Plasma-Wall Interactions in Controlled Fusion_ (Springer, 1986) pp. 695–771.
  * Reiter  _et al._ [2005] D. Reiter, M. Baelmans, and P. Boerner, Fusion science and technology 47, 172 (2005).
  * Stotler and Karney [1994] D. Stotler and C. Karney, Contributions to Plasma Physics 34, 392 (1994).
  * Hakim  _et al._ [2020] A. Hakim, M. Francisquez, J. Juno, and G. W. Hammett, Journal of Plasma Physics 86, 905860403 (2020).
  * Wersal [2017] C. Wersal, (2017).
  * Ghoos  _et al._ [2016] K. Ghoos, W. Dekeyser, G. Samaey, P. Börner, and M. Baelmans, Journal of computational physics 322, 162 (2016).
  * Boeyaert  _et al._ [2022] D. Boeyaert, S. Carli, K. Ghoos, W. Dekeyser, S. Wiesen, and M. Baelmans, Nuclear Fusion 63, 016005 (2022).
  * LeVeque [2007] R. J. LeVeque,  _Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems_ (SIAM, 2007).
  * Roache [1976] P. J. Roache, Computational Fluid Dynamics (1976).
  * Rognlien  _et al._ [1992] T. Rognlien, J. Milovich, M. Rensink, and G. Porter, J. Nuc. Mat. 196, 347 (1992).
  * Rognlien  _et al._ [1999] T. D. Rognlien, D. D. Ryutov, N. Mattor, and G. D. Porter, Physics of Plasmas 6, 1851 (1999).
  * Stotler  _et al._ [2022] D. Stotler, C. Karney, G. Wilkie, R. Kanzleiter, and S. Jaishankar, (2022).