A fully implicit, scalable, conservative nonlinear relativistic Fokker-Planck 0D-2P solver for runaway electron
Don Daniel William T. Taitano Luis Chacón





A fully implicit, scalable, conservative nonlinear relativistic Fokker-Planck
0D-2P solver for runaway electrons 
    Don Danielcorresponding author email: dond@lanl.gov,
William T. Taitano, Luis Chacn
    
=================================================================================================================





Los Alamos National Laboratory, New Mexico, USA



Upon application of a sufficiently strong electric field, electrons
break away from thermal equilibrium and approach relativistic speeds.
These highly energetic ‘runaway’ electrons
( MeV) play a significant role in tokamak disruption
physics, and therefore their accurate understanding is essential to
develop reliable mitigation strategies. For this purpose, we have
developed a fully implicit solver for the 0D-2P (i.e., including two
momenta coordinates) relativistic nonlinear Fokker-Planck equation
(rFP). As in earlier implicit rFP studies (NORSE, CQL3D), electron-ion
interactions are modeled using the Lorentz operator, and synchrotron
damping using the Abraham-Lorentz-Dirac reaction term. However, our
implementation improves on these earlier studies by 1) ensuring exact
conservation properties for electron collisions, 2) strictly preserving
positivity, and 3) being scalable algorithmically and in parallel.
Key to our proposed approach is an efficient multigrid preconditioner
for the linearized rFP equation, a multigrid elliptic solver for the
Braams-Karney potentials [Braams and Karney, Phys. Rev. Lett.
59, 16 (1987)], and a novel adaptive technique to determine
the associated boundary values. We verify the accuracy and efficiency
of the proposed scheme with numerical results ranging from small electric-field
electrical conductivity measurements to the accurate reproduction
of runaway tail dynamics when strong electric fields are applied.




§ INTRODUCTION


Relativistic Coulomb collisions are modeled using an extended version
of the Landau-Fokker-Planck collision operator <cit.>. Similarly
to its non-relativistic counterpart <cit.>, the operator
assumes small-angle collisions, is well-posed, and features strict
conservation of total particle number, total momentum, and total energy
<cit.>. However, its accurate numerical solution is difficult
because of its integro-differential formulation, which introduces
scalability and discretization challenges. In this study, we propose
a finite-difference-based conservative, parallel, fully implicit solver
for the 0D-2P relativistic Fokker-Planck (rFP) electron-electron collision
operator. This work builds and improves on earlier rFP algorithms
as implemented in the NORSE <cit.> and CQL3D <cit.>
codes.

The solver proposed in this study is primarily designed to simulate
runaway electrons produced by a large loop voltage. In tokamaks, the
loop voltage induced during disruptions can produce a large amount
of runaway electrons, which may severely damage plasma facing materials
<cit.>. The generated runaway current is also affected
by secondary mechanisms such as energy transfers from the primary
runaway electron current to the thermal electrons through knock-on
(large-angle) collisions. Understanding these nonlinear mechanisms
may be essential to develop either avoidance or mitigation strategies
for runaway electrons in tokamaks. 

A solver designed to capture runaway-electron dynamics benefits from
certain features. For example, capturing small-amplitude tails necessitates
strict positivity preservation. Runaway-electron generation time may
be large: a sizable runaway tail length may take hundreds of electron-electron
thermal collision-time scales to develop. Therefore, an implicit solver
that can step over stiff thermal collision-time scales is desirable.
The ability to use large time steps also demands that the scheme be
asymptotic preserving, which in turn requires enforcing strict conservation
properties <cit.>. It is also essential that the solver be
optimal and scales with the number of mesh points, as resolving small-scale
features may require fine grids while fitting long tails may require
large domains. The relativistic Fokker-Planck operator can be expressed
either in integral form <cit.> or in differential form <cit.>.
We employ the differential form, in which the collisional coefficients
are expressed in terms of relativistic potentials. This form is more
conducive to an optimal 𝒪(N) solver, where N is
the number of grid points, as the integral form produces an 𝒪(N^2)
scaling when symmetry is preserved in the integral operators (which
is required to achieve strict conservation properties <cit.>).
To obtain a nearly optimal scaling (𝒪(N^αlog N),
with α≳1) with the differential formulation, we propose
a multigrid-preconditioned GMRES <cit.> solver for
the potentials along with an adaptive treatment for evaluating potential
far-field boundary conditions [which greatly decreases their computational
complexity from 𝒪(N^3/2) to 𝒪(N^1.1)].
We note that it is in principle possible to improve the scaling resulting
from the integral form of the collision operator by the use of optimal
integral methods such as fast multipole methods <cit.>. However,
such methods can break the numerical symmetry of the integrals, resulting
in the loss of strict conservation properties. In this regard, fast
integral methods would not improve on the proposed optimal algorithm
for the differential formulation, and would require similar strategies
to recover the strict conservation properties of the collision operator.

With regard to time-stepping, we propose a conservative, fully implicit
nonlinear scheme, which as a result is asymptotic preserving (i.e.,
it captures the Maxwell-Jttner distribution as a steady-state solution
to our system). Earlier algorithmic approaches proposed for this system
are either linearly implicit (e.g., NORSE <cit.>), or lack
strict conservation properties (e.g., CQL3D <cit.>).
The implicit solver proposed in this study satisfies discrete conservation
properties, is preconditioned for optimal algorithmic performance,
and is scalable in parallel. As a result, our algorithm scales as
𝒪(N^1.1log N). Our conservation and preconditioning
strategies follow closely those proposed in Ref. <cit.>.

The rest of the paper is organized as follows. In <ref>,
we discuss the full relativistic electron-electron operator, the Lorentz
operator for electron-ion interactions, and the Abraham-Lorentz-Dirac
reaction term for modeling losses due to synchrotron damping. Then,
in <ref>, we discuss the algorithmic aspects
with regard to the discrete conservation strategy, positivity preservation,
and our near-optimal strategy for determination of the potentials.
In <ref>, we briefly describe
our fully implicit nonlinear solver using an Anderson Acceleration
scheme. In <ref>, we discuss the numerical
results that demonstrate the correctness of our implementation, and
finally in <ref> we list the conclusions
and scope for future work. 



§ FORMULATION


We model a homogeneous quasi-neutral plasma. We evolve the electron
species with the relativistic Fokker-Planck equation for the electron
distribution function, f_e, in the presence of background species
β, 


    ∂_tf_e+∂_p⃗·[(E⃗+F⃗_S)f_e]=∑_β=i,eC(f_β,f_e),

where t is time normalized with the relativistic electron collision
time, 

    τ_ee^relativistic=4πϵ_0^2m_e^2c^3/q_e^4n_elnΛ_ee,

p⃗ is the momentum vector normalized with m_ec, m_e
is the electron mass, c is the speed of light, q_e is the
electron charge, E⃗ is the electric field normalized with
the critical value for runaway electron generation <cit.>, E_c=n_eq_e^3lnΛ_ee/4πϵ_0^2m_ec^2,
n_e is the electron number density, ϵ_0 is the electrical
permittivity, lnΛ_ee is the Coulomb logarithm, F⃗_S
refers to the electron friction coefficients associated with synchrotron
radiation damping effects (defined in detail later), and C is the
collision operator given by, 


    C(f_β,f_e)=∂_p⃗·[D_β·∂_p⃗f_e-m_e/m_βF⃗_⃗β⃗f_e],

where D_β represents the collisional
diffusion tensor coefficients and F⃗_β represents
the collisional friction vector coefficients (computed based on the
appropriate background species f_β). Though Eq. (<ref>)
in principle may be used for multiple species, here we only consider
the evolution of electrons interacting with themselves, ions and external
electric fields. 

The distribution function is described in a two-dimensional cylindrical
domain (p_∥,p_⊥), with the subscripts ∥
and ⊥ referring to directions parallel and perpendicular to
the magnetic field, respectively (see Fig. <ref>). The
azimuthal direction is ignored because the distribution is axisymmetric.
The electron-electron interactions are described using the full form
of the collision operator, while the electron-ion interaction is modeled
with the Lorentz operator (which assumes the ions to be cold and infinitely
massive, m_i>>γ m_e with γ=√(1+p^2) the
Lorentz factor). 






 §.§ Electron-electron collisions


The collisional coefficients, D_β
and F⃗_β, for electrons in Eq. (<ref>)
are expressed in terms of the Braams-Karney potentials <cit.>.
These potentials are obtained by inverting a set of elliptic equations.
In this study, the elliptic solves are performed optimally with parallel
multigrid-preconditioned GMRES techniques, with a scaling of 𝒪(Nlog N).
The collisional coefficients are given by <cit.>:


    D_e    =     -4π/n_βγ^-1[L+P]h_1+4γ^-1[L-P]h_2,
    F⃗_⃗e⃗    =     -4π/n_βγ^-1K⃗(g_0-2g_1),

where the operators L, K⃗,
and P are defined as:

    Lψ    =    P·∂^2ψ/∂p⃗∂p⃗·P-P(p⃗·∂ψ/∂p⃗),
    K⃗ψ    =    P·∂ψ/∂p⃗,
    P    =    I+p⃗p⃗.

To obtain the transport coefficients, we first compute the h potentials
by solving the partial differential equations, 


    [L+1]h_0    =     f_e,
    
    [L-3]h_1    =     h_0,
    
    [L-3]h_2    =     h_1,

and then the g potentials by solving: 


    Lg_0    =     f_e,
    
    Lg_1    =     g_0.

Here, the operator L is defined as:

    L ψ=P:∂^2ψ/∂p⃗∂p⃗+3p⃗·∂ψ/∂p⃗.

To solve these linear potential equations, we require far-field boundary
conditions. They are determined from the Green's function solution
of the elliptic equations, Eqs. (<ref>,<ref>)
<cit.>:


    h_0    =     -1/4π∫(r^2-1)^-1/2f_β(p⃗'⃗)/γ'd^3p⃗'⃗,
    
    h_1    =     -1/8π∫√(r^2-1)f_β(p⃗'⃗)/γ'd^3p⃗'⃗,
    
    h_2    =     -1/32π∫(rcosh^-1r-√(r^2-1))f_β(p⃗'⃗)/γ'd^3p⃗'⃗,
    
    g_0    =     -1/4π∫ r(r^2-1)^-1/2f_β(p⃗'⃗)/γ'd^3p⃗'⃗,
    
    g_1    =     -1/8π∫cosh^-1rf_β(p⃗'⃗)/γ'd^3p⃗'⃗,

where r=γγ'-p⃗·p⃗'⃗. Note that the integral
kernels of h_0 and g_0 are singular when r→1 (p⃗→p⃗'⃗),
which require a specialized numerical treatment in terms of elliptic
integrals for accuracy and efficiency (see <ref>
and App. <ref>). Also, we have devised
an efficient adaptive algorithm to fill ghost cells that prevents
these boundary integrals from leading to an 𝒪(N^3/2)
scaling of the computational complexity (see also <ref>).



 §.§ Modeling external effects


We consider several external effects, including an imposed electric
field, E⃗=(E_∥,0), ions, and synchrotron radiation,
F⃗_S. 

Electron-ion scattering is modeled with the Lorentz or pitch-angle
scattering operator <cit.>, which
assumes ions are cold and infinitely massive. The operator causes
scattering of the electrons in the pitch angle (arccos[p_∥/p])
and, in this simplified form, it preserves kinetic energy. It has
finite diffusion coefficients and zero friction coefficients, given
by:


    D_i,∥∥=Z_eff/2vp_⊥^2/p^2,     D_i,∥⊥=D_i,⊥∥=-Z_eff/2vp_⊥p_∥/p^2,     D_i,⊥⊥=Z_eff/2vp_∥^2/p^2,    F⃗_i=0⃗,

where p^2=p_⊥^2+p_∥^2, v is the velocity
magnitude (normalized with c), and Z_eff=∑ n_iZ_i^2/∑ n_iZ_i
is the effective ion-charge state (n_i and Z_i refers to
ion densities and charges). For a quasi-neutral plasma, ∑ n_iZ_i=n_e.
Note that the electron-ion collision operator becomes singular at
the origin v→0. We mollify this singularity by reformulating
the singular part, as:


    1/v≈1/√(v^2+v_cut^2),

where v_cut=p_cut/√(1+p_cut^2) is the velocity cut-off,
with p_cut=2Δ p. Note this approximation of the singular
term in the cylindrical space introduces a finite but small amount
of heating as p→0. 

Finally, we consider synchrotron radiation, which results in loss
of momentum for the electrons. We model this with the Abraham-Lorentz-Dirac
reaction term <cit.>. The reaction term
has finite friction coefficients, given by: 


    F_S,⊥=-Sp_⊥/γ(1+p_⊥^2),     F_S,∥=-Sp_∥/γp_⊥^2,

where S=τ_ee^relativistic/τ_r relates the time scale
of the synchrotron-radiation damping, τ_r, to that of electron-electron
relativistic collisions, τ_ee^relativistic.



§ ALGORITHM




 §.§ General discretization strategy


We employ a conservative finite-difference scheme. The distribution
is evaluated at cell centers, while the friction and diffusion fluxes
are evaluated at cell faces. Recall the electron-electron collision
operator is the divergence of a collisional flux,


    C(f_e,f_e)=∂_p⃗·(D_e∇_pf_e-F⃗_⃗e⃗f_e)≈δ_p⃗·(R⃗_D-R⃗_F)=δ_p⃗·R⃗,

where R⃗_D and R⃗_F are the diffusion and
friction fluxes, respectively, and δ_p⃗· denotes
the discrete form of the divergence operator, which in cylindrical-momentum
space is written as:


    (δ_p⃗·R⃗)_j,k=(R_∥,j+1/2,k-R_∥,j-1/2,k/Δ p_∥,j+p_⊥,k+1/2R_⊥,j,k+1/2-p_⊥,k-1/2R_⊥,j,k-1/2/p_⊥,kΔ p_⊥,k).

Fluxes at cell faces are given by:


    R_D,∥,j+1/2,k=(D_∥∥∂_p_∥f_e+D_∥⊥∂_p_⊥f_e)_j+1/2,k,   R_D,⊥,j,k+1/2=(D_⊥∥∂_p_∥f_e+D_⊥⊥∂_p_⊥f_e)_j,k+1/2,



    R_F,∥,j+1/2,k=F_∥,j+1/2,kf_e,j+1/2,k,   R_F,⊥,j,k+1/2=F_⊥,j+1/2,kf_e,j,k+1/2.


The potential operator L (Eq. <ref>) is discretized
using central differences (see App. A1 for details). The potentials
are evaluated at cell centers and their boundary conditions are specified
at ghost cells. For the potential Eqs. (<ref>)-(<ref>),
we apply far-field Dirichlet boundary conditions using Eqs. (<ref>)
as discussed in <ref>. The collisional
coefficients, D_β and F⃗_β,
are evaluated at cell centers using the computed potentials (see App.
A2 for discretization details). The collisional coefficients at the
ghost cells are evaluated by linearly extrapolating the values from
adjacent cell-centered values. Note that the ghost cells also store
distribution and coefficient data for cross-processor communication
using an MPI framework. 

The coefficients for external effects (such as scattering due to ion
interactions, Eq. (<ref>), synchrotron
damping effects, Eq. (<ref>), and electric
field acceleration terms) are evaluated at cell-centers. Where needed,
values at cell faces are found by averaging two adjacent cell-centered
values within the computational domain. As is typical in kinetic simulations,
the outer domain boundaries are selected such that the distribution
is sufficiently small there.



 §.§ Discrete conservation strategy for the e-e collision operator




The relativistic electron-electron collision operator conserves the
total particle number, momentum, p⃗=γv⃗, and
energy E=γ, as the moments of the collision operator satisfy:


    ⟨1,C(f_e,f_e)⟩_p    =     0,
    ⟨ p_∥,C(f_e,f_e)⟩_p    =     0,
    ⟨γ,C(f_e,f_e)⟩_p    =     0,

where ⟨ a,b⟩_p=∫_pab2π p_⊥dp_∥dp_⊥.
Discretely, these inner products may be approximated by a mid-point
quadrature rule as: 


    ⟨ A,B⟩_p^D≈2π∑_j=1^N_∥∑_k=1^N_⊥A_j,kB_j,kp_⊥,kΔ p_∥,jΔ p_⊥,k,

where the superscript D refers to the discrete representation of
the summation operator, and Δ p_∥,j and Δ p_⊥,k
are the width and height of a rectangular cell located at (j,k).

In general, the relationships shown in Eqs. (<ref>-<ref>)
will not be satisfied due to numerical errors. Discrete particle number
conservation (Eq. <ref>) is trivially satisfied
by setting the normal component of diffusion and friction fluxes to
zero at the boundary. However, enforcing Eqs. (<ref>,<ref>)
is more challenging. A recent study <cit.> enforced these
conservation properties discretely by redistributing the numerical
errors via discrete nonlinear constraints. We employ a similar methodology
here. Firstly, we multiply the diffusion flux by a factor


    η=1+η_0+η_1(p_∥-p̅_∥),

where the magnitudes of η_0 and η_1 are expected to
be of the order of truncation error, and p̅_∥=⟨ f_e,p_∥⟩_p/⟨1,f_e⟩_p
is the mean momentum. Thus, the discrete collision operator is of
the form,


    C^D(f_e,f_e)=δ_p⃗·(η⃗R⃗_⃗D⃗-R⃗_F),

Integrating over the cylindrical-momentum domain, 


    ⟨ p_∥,C^D(f_e,f_e)⟩_p^D    =     0,
    ⟨γ,C^D(f_e,f_e)⟩_p^D    =     0,

we obtain a system of two equations,


    [[                ⟨γδ_p⃗·R⃗⟩_p^D    ⟨γ(p_∥-p̅_∥)δ_p⃗·R⃗_D⟩_p^D;           ⟨ p_∥δ_p⃗·R⃗_D⟩_p^D ⟨ p_∥(p_∥-p̅_∥)δ_p⃗·R⃗_D⟩_p^D ]][[ η_0; η_1 ]]=[[    ⟨γδ_p⃗·(R⃗_F-R⃗_D)⟩_p^D; ⟨ p_∥δ_p⃗·(R⃗_F-R⃗_D)⟩_p^D ]],

for unknowns [η_0, η_1], which can be inverted straightforwardly.
This strategy conserves momentum and energy at the discrete level
for electron-electron collisions. Note that, because we assume the
ions to be cold and infinitely massive, there are no conservation
properties associated with electron-ion collisions (i.e. discrete
representation may result in a net energy and momentum loss).



 §.§ Time stepping strategy


A huge separation in time scales exists in runaway-electron dynamics.
Long time-scales are of the order of the relativistic collision times,
𝒪(τ_ee^relativistic). For typical bulk temperatures
Θ=T/m_ec^2∼10^-4, this implies a time-scale separation
of six orders of magnitude between thermal and relativistic time scales
(as τ_ee^thermal=Θ^3/2τ_ee^relativistic).
Stepping over fast time scales demands a fully implicit temporal scheme
with strict conservation and positivity preservation properties. We
describe our approach next.

The discrete system of equations representing the effects of electron-electron
collisional interactions C and external effects, ℰ,
on electron evolution can be written as:


    δ_tf_e^n=C^D(f_e^n,f_e^n)+C̃(f_i,^nf_e^n)-δ_p⃗.[(F⃗_S+E⃗)f_e^n]_ℰ(f_e^n),

where the superscript D represents the appropriate discrete form
defined in <ref>, and
C̃ represents the Lorentz operator (see Eq. <ref>).
For a general implicit backward time discretization scheme at time
step n, we have,


    δ_tf_e^n=∑_i=0,1,2...b_if_e^n-i/Δ t,

where constants, b_i, satisfy ∑_ib_i=0. We use both
first-order (Euler, BDF1) and second-order (BDF2) schemes for time
advancing. For BDF1, b_0=-1 and b_1=1 and for BDF2 with
constant time steps, b_o=3/2, b_1=-2, b_2=1/2. The
coefficients can be generalized for non-uniform time steps. 

Multiplying Eq. (<ref>) with c⃗=(1,p_∥,γ),
and averaging over the momentum space, we obtain: 


    ∑_i=0,1,2...b_i⟨c⃗,f_e^n-i⟩_p^D/Δ t=⟨c⃗,C^D(f_e^n,f_e^n)⟩_p^D+⟨c⃗,ℰ(f_e^n)⟩_p^D.

Because of the discrete conservation properties of the electron-electron
collision operator, the first term in the right hand side vanishes.
Therefore, any change in the total momentum or energy of electrons
can only be due to external effects such as ion-electron collisions,
synchrotron radiation, and electric field acceleration. 



 §.§ Positivity-preserving strategy


Positivity-preserving schemes are essential to capture small-amplitude
runaway tails. Our strategy is to leverage the structure of the differential
operators (advection-diffusion), and use existing positivity-preserving
discretizations for these terms.

For all advective terms in the relativistic kinetic equation, we use
the positivity-preserving SMART flux limiter <cit.> to construct
the associated fluxes. For the diagonal components of the tensor diffusion
term, D·∇_pf_e |_∥∥
and D·∇_pf_e |_⊥⊥,
we employ a standard second-order discretization:


    (D·∇_pf_e )_∥∥,j+1/2,k=(D_∥∥∂_p_∥f_e)_j+1/2,k≈D_∥∥,j+1,k+D_∥∥,j,k/2f_e,j+1,k-f_e,j,k/Δ p_∥,j+1/2,



    (D·∇_pf_e )_⊥⊥,j,k+1/2=(D_⊥⊥∂_p_⊥f_e)_j,k+1/2≈D_⊥⊥,j,k+1+D_⊥⊥,j,k/2f_e,j,k+1-f_e,j,k/Δ p_⊥,k+1/2,

which is numerically well-posed (does not feature a null space and
features a maximum principle). However, unless care is taken, the
off-diagonal diffusion tensor terms do not feature a discrete maximum
principle, resulting in loss of boundedness. To address this issue,
we reformulate the off-diagonal components as effective friction forces
as proposed in Ref. <cit.>:


    (D·∇_pf_e )_∥⊥    =     D_∥⊥∂_p_⊥f_e=f_eD_∥⊥∂_p_⊥ln f_e_F_∥^eff=f_eF_∥^eff=R_F,∥^eff,
    
    (D·∇_pf_e )_⊥∥    =f_e   D_⊥∥∂_p∥ln f_e_F_⊥^eff=f_eF_⊥^eff=R_F,⊥^eff.

Once formulated as advective terms, we use flux-limiting advective
schemes (similar to the collisional friction terms) to calculate the
effective flux. Discretization details can be found in App. A3. 



 §.§ Strategy for evaluating boundary conditions of collision potentials 


The boundary conditions for the relativistic potential equations for
h and g are found using the integral formulations in Eq. (<ref>).
For the h_1,h_2,g_1 relativistic potentials, we use a trapezoidal-rule
numerical integration with 24 discrete points in the ϕ angle.
However, the kernels in g_0,h_0 become singular in the limit
of p'→ p r→1, complicating a direct
numerical integration. These complexities can be eliminated by reformulating
these integrals in terms of complete elliptic integrals. We begin
by noting that, because the distribution is axisymmetric, the 3D momentum-space
integration of the Green's function can be rewritten as a 2D momentum
space integration over the PDF and a 1D azimuthal angle integration
as:


    h_β,0    =     -1/4π∫f_β(p'_∥,p'_⊥)/γ'p_⊥dp'_∥dp'_⊥∫_ϕ1/(r^2-1)^1/2dϕ_I=-1/4π∫f_β(p'_∥,p'_⊥)/γ'I(p_∥,p_⊥,p'_∥,p'_⊥)p_⊥dp'_∥dp'_⊥,
    
    g_β,0    =     -1/4π∫f_β(p'_∥,p'_⊥)/γ'p_⊥dp'_∥dp'_⊥∫_ϕr/(r^2-1)^1/2dϕ_H=-1/4π∫f_β(p'_∥,p'_⊥)/γ'H(p_∥,p_⊥,p'_∥,p'_⊥)p_⊥dp'_∥dp'_⊥,

The segregated integrals I and H are then written in terms of
complete integrals of the first and third kind (see App. <ref>).


However, even after these reductions, evaluating potentials at all
ghost points in the boundary remains expensive. There are approximately
𝒪(N^1/2) ghost-cell boundary points, each point requiring
𝒪(N) integrals when using Eqs. (<ref>).
This makes the potential boundary evaluations scale poorly with the
number of mesh points, N [i.e., 𝒪(N^3/2)].
To ameliorate the scaling for the boundary-condition treatment, we
adaptively select a small number of boundary points for the potential
evaluations to match a given accuracy, with the remaining ghost points
found by interpolation using a high-order spline. The adaptive algorithm
to find the minimum number of spline knots needed for a given tolerance
is outlined in Algorithm 1, and illustrated in Fig. <ref>.
We begin with a set of uniformly distributed ghost points at the boundary,
for example, four points (black crosses in the first row), where we
evaluate the values of the potential integral. We fit a cubic (or
higher order) spline through these values (blue crosses in the second
row). New knots are then created by bisection (black crosses in the
third row), where integrals are again evaluated. The absolute error
is then computed as the difference between the value given by the
spline interpolation, ϕ_S, and the actual value of the potential
integral at the targeted points ϕ_I:


    a_b=|ϕ_I-ϕ_S|.

Intervals delimited by the set of knots that do not satisfy the prescribed
tolerance (e.g., red knot in the fourth row) are bisected further.
This process is continued until a spline fit of the desired accuracy
is obtained. 



To ensure the spline error is commensurate with other sources of error
in the algorithm, in practice the absolute tolerance criterion is
chosen to be a function of the momentum mesh spacing as: 


    a_b∼0.05Δ p_⊥Δ p_∥.


Figure <ref> illustrates the adaptive knots
generated with the adaptive spline algorithm for the g_0 potential
for a Maxwell-Jttner distribution of Θ=10^-2 in a mesh
of N_∥=2048 and N_⊥=512. The red dots along
the left, top, and right boundaries point to the location of the spline
knots generated using Algorithm 1. The contour of g_0 is also
illustrated to demonstrate the variation of g_0 at points close
and far away from the distribution. The algorithm generates more spline
knots where the function varies significantly. At the far right boundary,
the points are few and equally spaced, as the function variation is
small. A clear benefit of the adaptive spline approach can be seen
at the top boundary, where it is determined that only 41 functional
evaluations are needed for an accurate estimate of the potential along
the entire boundary, which spans a total of 2048 mesh points. 

Assuming equi-spaced knots, we can get an estimate of how the number
of splining knots, N_p, scales with the total degrees of freedom,
N, by comparing the spline error with the tolerance in Eq. <ref>:


    1/N∼(1/N_p)^n_s+1 N_p∝ N^1/(n_s+1).

Here, n_s is the order of the spine. For instance, for a fourth-order
spline, this result predicts N_p∼ N^0.2. However, we expect
this estimate to be very conservative, and it does not take into the
account the adaptive distribution of the knots. Numerical experiments
in Sec. <ref> show that N_p∼𝒪(N^0.13)
for a fourth-order spline. 



§ NONLINEAR SOLVER 


The spatial and temporal discretization techniques prescribed in <ref>
lead to a coupled nonlinear system of equations, which requires an
iterative nonlinear solver for the distribution function. We use an
Anderson Acceleration scheme <cit.> to converge iteratively
the system, which we briefly summarize next. 

Given a fixed point map based Picard iteration,


    f^k+1=G(f^k),

where the superscript k denotes the iteration step, Anderson Acceleration
scheme <cit.> accelerates the convergence of the
Picard iteration by using the history of past nonlinear solutions
via: 


    f^k+1=∑_i=0^m_kα_i^kG(f^k-m_k+i)_f^k-m_k+i+1,

where in this study m_k=min(5,k). The coefficients α_i^k
are determined via an optimization procedure that minimizes, 


    ‖∑_i=0^m_kα_i^k(G(f^k-m_k+i)-f^k-m_k+i)‖ ,

 subject to ∑_i=0^m_kα_i^k=1.

To enable preconditioning of the Anderson iteration, our fixed map
is based on a quasi-Newton iteration, where:


    f^k+1=G(f^k)=f^k+δ f^k=f^k-(P^k)^-1ℛ^k,

with P^k the preconditioner, δ f^k the nonlinear increment,
and ℛ^k the nonlinear residual. Given an electron
distribution, f_e, the residual for the nonlinear system is evaluated
as outlined in Algorithm 2. Note that if P is the Jacobian, i.e.
P^k=(∂ℛ/∂ f_e)^k, then Eq. (<ref>)
becomes a Newton iteration. 

The residual contribution from electron-electron collisions requires
the solutions of five potentials, which require inversions of the
linear equations in Eqs. (<ref>,<ref>).
These are inverted for each nonlinear iteration at flux-assembly time
along with the computation of conservation constraints η_0
and η_1, see Eq. <ref>. The nonlinear elimination
of the residuals associated with the potentials and conservation constraints
follows from previous studies <cit.>, and enables
a conservative, optimal 𝒪(Nlog N) solver when the
Poisson operators are inverted optimally and scalably. Here, the linear
potential equations are solved using a multigrid-preconditioned GMRES
<cit.> solver. The multigrid preconditioner features 1 V cycle
with 4 passes of damped Jacobi (damping factor of 0.7), along with
agglomeration for restriction and a second-order prolongation. At
the beginning of the solve, the five potentials are solved using a
tighter relative tolerance criteria of 10^-8 and then followed
by a looser relative-tolerance criteria of 10^-5-10^-7 during
each nonlinear solve, depending on the problem. 

The preconditioner in Eq. (<ref>) is obtained by
Picard linearization of the potentials and subsequent discretization
of the full system,


    P^kδ f=δ_tδ f-C(f_e^k-1,δ f)-ℰ(δ f),

where ℰ is a linear operator representing the net external
effects on electrons, see Eq. (<ref>). The transport
coefficients in the electron-electron collision operator, C, are
Picard-linearized and computed at the previous nonlinear iteration,
k. All advective terms in the preconditioner are discretized using
a linear upwinding scheme. During each nonlinear step k, the linear
system P^kδ f^k=-ℛ^k is solved with one
multigrid V-cycle and 3 passes of damped Jacobi (with damping constant
0.7). We use agglomeration for restriction and second-order prolongation. 

The nonlinear iteration ends when the desired relative nonlinear residual
convergence ratio r_NL is reached, 


    r_NL=ℛ^k/ℛ^k=0.

Cases with large disparities in signal amplitudes, for example a Maxwellian
thermal bulk along with runaway tail, may require a tighter convergence
ratio, r_NL=10^-7, to capture accurately the small-amplitude
tail. In contrast, a single deforming electron thermal bulk can use
a significantly looser nonlinear convergence, r_NL=10^-4, for
accurate results.






§ RESULTS


We begin this section with some verification studies, and finish it
with scalability and accuracy studies to assess the performance of
the algorithm.



 §.§ Verification


We first verify conservation properties of the equilibrium Maxwell-Jttner
distributions either at rest or moving with a mean momentum, p_b.
Note that all computations are performed in the stationary reference
frame. Then, we verify conservation properties during collisional
relaxation, and also benchmark the calculation of electrical conductivity
under the action of both weak and strong electric fields with previous
studies <cit.>. Finally, we verify our algorithm
with recent calculations of runaway dynamics in the t→∞
limit <cit.>. Verification results
were obtained using the second-order BDF2 time-stepping scheme except
when explicitly stated (see <ref>
for details on the time-stepping scheme).



  §.§.§ Preservation of stationary and boosted Maxwell-Jttner distributions. 


The computational domain is uniform with N_∥=256 and
N_⊥=128. The nonlinear residual, r_NL, is converged
to a relative tolerance of 10^-4 unless otherwise specified.
The discrete conservation properties are satisfied to nonlinear tolerance
and are independent of the time step used. The electron number density
is normalized, n_e=1. The domain is chosen such that the distribution
function is sufficiently small at boundaries. The entire domain is
shown in the figures illustrating the distribution function. Fig.
<ref>(a) illustrates a static Maxwell-Jttner
(MJ):


    f_e^MJ=n_e/4πΘ K_2(1/Θ)exp[-γ(p)/Θ],

in log scale with normalized temperature Θ=T/m_ec^2=1.
In Eq. <ref>, K_2 is the modified Bessel function of
the second kind. We have confirmed that the distribution retains its
initial shape for the whole simulation. To illustrate this, Fig. <ref>(b)
demonstrates the evolution of the relative errors in the number density,
momentum and energy for 200 τ_ee. The relative errors are
measured as, 


    relative error=| g(t)-g(0)|/g(0).

where g is either the number density, momentum or energy. The figure
shows that the relative errors in number density are one part in 10^11,
while errors in relativistic momentum and energy remain smaller than
one part in 10^8. 


For a boosted (translated) MJ, the equilibrium distribution appears
deformed and is given by: 


    f_e^BMJ=n_e/4πΘ_bγ_bK_2(1/Θ_b)exp[-γ_bγ-p_bp_∥/Θ_b],

where the subscript b denotes the values in the boosted frame.
Fig. <ref> illustrates a boosted MJ equilibrium distribution
with Θ_b=0.15, in a frame boosted by p_b=2 and with
γ_b=√(1+p_b^2). The relative errors are shown
for 10τ_ee^relativistic∼100τ_ee^thermal, demonstrating
identical behavior as in the stationary MJ case. 






  §.§.§ Conservation properties during collisional relaxation dynamics


To explore collisional relaxation dynamics, we consider two cases,
one which features an initial configuration of two boosted MJ distributions,
and the other which features a randomized initial distribution. Simulations
have been run till the distributions have relaxed to a single Maxwell-Jttner.

Figure <ref>(a) illustrates the collisional relaxation
of two MJ distributions boosted by 2 units in opposite directions.
Fig. <ref>(b) depicts the relative errors in number
density, momentum, and energy. After an initial transient stage t∈(0,300),
the relative errors in momentum and number density remain small and
bounded in time. 


To demonstrate that the discrete conservation strategy works in more
complicated cases, we consider the thermalization of a random distribution,
of the form:


    f_e^rand=𝒫/⟨𝒫,1⟩_p,    where    𝒫(p_∥,p_⊥)=J(p_∥,p_⊥)/4πΘ K_2(1/Θ)exp[-γ(p)/Θ],

where Θ=1 and J is a random number function with a uniform
distribution in the range [0,1], see Fig. <ref>.
The presence of large gradients in the distribution and small tails
makes this an excellent problem to test discrete conservation errors
and positivity preservation. For a nonlinear relative tolerance of
10^-4, the relative errors in momentum are larger than in previous
cases, 1 part in 10^4. Tightening the relative nonlinear tolerance
to 10^-6 results in a commensurate decrease of the errors to
1 part in 10^6. 






  §.§.§ Electrical conductivity in weak and strong electric fields


We consider next the case where collisional friction balances an externally
imposed electric field, leading to finite electrical conductivity.
We verify the code for a wide range of initial temperatures with electrical
conductivity results provided by Braams and Karney <cit.>. To
measure conductivity, we apply a small electric field, Ê_∥=10^-3E_D
where Ê_∥ is the parallel electric field, E_D=E_c/Θ
is the Dreicer field, and E_c the Connor-Hastie critical electric
field <cit.>. The electron distribution is initialized using
the Maxwell-Jttner distribution at various temperatures Θ
and the normalized electrical conductivity is computed as:


    σ̅=Z_eff/n_eq_eΘ^3/2ĵ/E_∥,        ĵ=-n_eq_ev_∥,

where v_∥=p_∥/γ. A small electric field
deforms the Maxwellian slightly to produce a net electron flow in
the positive p_∥ direction. To prevent numerical overflow,
for Θ10^-3 the initial distribution is defined using
a non-relativistic Maxwellian. Fig. <ref>(a) depicts
σ̅/Z_eff at various temperatures. The numerical
simulation results (circles) are in excellent agreement with the analytical
results (lines) from Ref. <cit.>. The electrical conductivity
measurements are made after the simulation reaches a quasi-steady-state
after an initial transient. Because of the applied electric field,
the plasma slowly heats up, and the quasi-steady-state temperatures,
Θ=⟨ f_e(p⃗,t)p^2/2γ⟩_p, are
larger (but close) to their initial value. The electrical conductivity,
σ̅, in Eq. (<ref>) is computed using the
quasi-steady-state temperature.



Figure <ref>(b) illustrates the time evolution
of electrical conductivity for various electric-field strengths in
the non-relativistic limit. Results of NORSE <cit.> (circles)
and Weng et al. <cit.> (squares) are also shown. The electron
distribution is initialized with a Maxwellian corresponding to an
initial temperature of Θ=10^-4 and n_e=1. Note that
the Weng et al. study uses the nonrelativistic Fokker-Planck operator.
For all values of electric field, we have good agreement with earlier
studies. For the case of Ê_∥=0.01E_D, we have
better agreement with NORSE than Weng et al.. Ref. <cit.> hypothesizes
that the observed deviations between NORSE and Weng et al. in the
small Ê/E_D limit may be due to numerical heating in
Weng et al.. Our results also suggest the same. 



  §.§.§ Reproducing runaway-electron tail dynamics


To verify runaway dynamics with existing linear test-particle studies
<cit.>, we performed two linearized numerical simulations
where we keep the collisional cofficients fixed in time to those of
a Maxwell-Jttner distribution with n_e=1 and Θ=0.01.
In the first simulation, we applied an electric field 2.25 times the
critical value, i.e E_∥=2.25. This causes some electrons
to overcome the frictional force and accelerate to high speeds. Fig.
<ref> demonstrates the evolution of the runaway
tail at 42000 τ_ee^thermal collision times. The asymptotic
slope of the runaway tail as predicted by Connor-Hastie <cit.>
is:


    f_e^tail∝1/p_∥exp(-(E_∥+1)p_⊥^2/2(1+Z_eff)p_∥).

As can be seen in the figure, the runaway tail produced by the algorithm
is in excellent agreement with the asymptotic theoretical results.


In the second simulation, we verify the runaway electron dynamics
in the presence of the synchrotron radiation damping term in Eq. (<ref>).
Because we are interested in the steady-state as t→∞,
time marching is performed efficiently with a BDF1 time stepping scheme
(rather than BDF2). Fig. <ref> shows the electron
distribution function at t=777 (i.e., ≈777,000 τ_ee^thermal)
with a damping coefficient of S=0.1. Other parameters are E_∥=2.25,Θ=0.01,
and Z_eff=1. Ref. <cit.> performed
a linearized initial value problem and describes the evolution of
runaway electrons in the momentum space as a two-step phenomenon,
beginning with the relative fast formation of a long runaway tail
and a much slower rise of the bump to a steady-state solution. We
find similar behavior here with the electrons accumulating in the
momentum space around p_0≈18, to form a second maximum.
The location of the second maximum is in good agreement with Ref.
<cit.>. Note that, when collisional coefficients are evolved
nonlinearly, we see heating of the bulk (not shown) leading to slide-away
effects as previously reported in Ref. <cit.>. 






 §.§ Solver performance




  §.§.§ Algorithmic and parallel scalability.


Table <ref> lists weak parallel
scalability results for the thermalization of a random electron distribution,
see f_e^rand in Eq. (<ref>), with Θ=1.
For simplicity, we employ the BDF1 time-stepping scheme for scalability
tests. We report the wall clock time (WCT) per time step, the number
of nonlinear iterations (NLI) per implicit time step Δ t,
the ratio WCT/NLI (which is an indirect measure of communication costs),
and the ratio between implicit and explicit time steps (which is a
measure of numerical stiffness). The potential linear iterations (PLI)
is the average number of GMRES iterations required for each potential
solve. The absolute tolerances are set to 10^-8. For the initial
solve, the potentials are converged tightly to a relative tolerance
of 10^-8, with a looser tolerance of 10^-6 used for subsequent
nonlinear iterations. The adaptive spline tolerance is defined by
Eq. (<ref>) and a fourth-order spline is used
for piece-wise interpolation. The explicit time step is calculated
as:


    Δ t_explicit=0.25min_β=i,e(Δp_∥^2/max(D_β,∥∥),Δp_⊥^2/max(D_β,⊥⊥),Δp_⊥/max(A_β,⊥),Δp_∥/max(A_β,∥)).

Table <ref> demonstrates excellent
parallel scalability (WCT/NLI) up to 4096 processors. When increasing
the number of processors while keeping the problem size per processor
constant, we are effectively increasing the resolution of the problem,
thus making the problem harder to solve for a fixed time step (as
evidenced by the increasing implicit-to-explicit timestep ratio).
This manifests in a very mild growth of the number of nonlinear iterations
(NLI) as we increase the number of cores (the iteration count increases
only by a factor of 3 over a three-order-of-magnitude increase of
the problem size).Table <ref> lists
parallel and algorithmic scalability results for the case of collisional
relaxation of two Maxwell-Jttner distributions, boosted by one momentum
unit in opposite directions and a normalized temperature of 10^-1
in their frames of reference, see f_e^BMJ in Eq. (<ref>).
We observe good parallel (WCT/NLI) and algorithmic scalability (NLI)
up to 4096 processors with Δ t/Δ t_explicit∼460
for the high resolution case of 4096×2048. 

Fig. <ref>(a) illustrates weak scaling results of the
wall clock time per nonlinear iteration (WCT/NLI) vs. the number of
cores for the random electron thermalization (red line, Table <ref>)
and the boosted MJ relaxation (blue line, Table <ref>).
We observe excellent parallel scalability in both cases. The expected
scaling for parallel multigrid-based solvers is WCT/NLI ∼𝒪(log N),
resulting in an overall algorithmic scaling of 𝒪(NlogN).
Instead, we find WCT/NLI ∼𝒪(N^0.1logN). The additional
factor of N^0.1 originates in the growth of the number of spline
knots with N, as predicted by Eq. <ref> and demonstrated
in Fig. <ref>(b). In the figure, we show that the cumulative
number of spline knots (for all potentials) increases as O(N^0.13)
for a fourth-order spline, which is more benign than the predicted
one in Eq. <ref> (which assumed uniformly spaced
knots), and much more efficient than the naive scaling of 𝒪(N^3/2).
We have also confirmed that the scaling of number of spline knots
with N depends on the spline order, with lower orders resulting
in a larger exponent. For instance, a cubic spline for the same tests
results in N_p scaling as O(N^0.17) (results not shown).
It is interesting to note that, from the figure, the 𝒪(N^0.1)
scaling seems to disappear at large core count (large problem sizes),
which we speculate is due to the problem becoming too well resolved
by the mesh.








  §.§.§ Spatial and temporal accuracy 




Figure <ref>(a) illustrates the spatial
accuracy of the proposed scheme with the boosted MJ collisional relaxation
problem described in <ref>.
The `exact' electron distribution f_e^exact is obtained at
t=0.15 (before steady-state is reached) with N_∥=1024,
N_⊥=512 and Δ t=0.01, while the three data points
(solid dots) correspond to coarser grids of 128×64, 256×128
and 512×256. The blue dashed line corresponds to a second-order
error scaling. The ℓ_2-norm of the error between the `exact'
and numerical electron distributions is computed as: 


    f_e-f_e^exact_2=(∑_j=1^N_∥∑_k=1^N_⊥(f_e,j,k-f_e,j,k^exact)^22π p_⊥,kΔ p_∥,jΔ p_⊥,k)^0.5.

We confirm that the proposed implementation is second-order accurate
in space. 

Figure <ref>(b) illustrates the temporal
accuracy of the implementation. In this case, the 'exact' electron
distribution f_e^exact is obtained in a 256×128 grid
using the BDF2 time advancement scheme with Δ t=5×10^-4,
see description in <ref>. The
four data points correspond to larger time steps of Δ t=10^-4,2×10^-4,4×10^-4,
and 10^-3. The proposed implementation is confirmed to be second-order
accurate in time.



§ CONCLUSIONS


We have developed a fully implicit, nearly optimal, relativistic nonlinear
Fokker-Plank algorithm with strict conservation properties.. We consider
a 0D2P cylindrical momentum-space representation. The solver employs
the differential form of the Fokker-Planck equation, which requires
the solution of five relativistic potentials in momentum space to
obtain the collisional coefficients. Singularities in the potential
integral formulations are resolved by expressing them in terms of
complete elliptic integrals of the second and third kind. To ensure
a benign scaling of the potential solves with the total number of
mesh points N, we employ a multigrid-preconditioned GMRES solver,
and have developed an adaptive spline technique for finding far-field
boundary conditions for the potentials. The adaptive spline technique
results in a small additional exponent in the algorithmic scaling
of 𝒪(N^0.1). Positivity of the distribution function
is ensured using a continuum-based reformulation approach <cit.>
combined with robust positivity-preserving discretizations schemes
<cit.>. Using an Anderson Acceleration fixed-point iteration
scheme for our nonlinear solves, also preconditioned with multigrid
techniques, we obtain an algorithm that overall scales as 𝒪(N^1.1logN).
The logN contribution is due to the parallel multigrid techniques
employed, and the N^0.1 contribution is from the proposed adaptive
spline technique. We have demonstrated second-order accuracy in both
space and time, and characterized the performance of our parallel
implementation. We have verified our solver by comparing with previous
results for electrical conductivity measurements in the weak and strong
electric field limits. We have demonstrated the accuracy of conserved
quantities in electron-electron collisions, with small relative errors,
in number density, relativistic momentum, and energy. In addition,
we have examined runaway dynamics and verified it by comparing to
known results <cit.>. In future work,
we will extend this method to study inhomogeneous plasmas by considering
the spatial dependence of the electron distribution function.



§ ACKNOWLEDGEMENTS


The authors thank Z. Guo for help in verifying the algorithm and E.
Hirvijoki for insightful discussions on the properties of the relativistic
operator. The authors also thank C. McDevitt and X. Tang for useful
inputs during the course of this project. This work was supported
by the US Department of Energy through the Los Alamos National Laboratory.
Los Alamos National Laboratory is operated by Triad National Security,
LLC, for the National Nuclear Security Administration of U.S. Department
of Energy (Contract No. 89233218CNA000001). This research used resources
provided by the Los Alamos National Laboratory Institutional Computing
Program.





§ DISCRETIZATION OF OPERATORS IN POTENTIAL EQUATIONS AND COLLISIONAL
COEFFICIENTS




 §.§ A1. Linear potential equations. 


The potential operator L consists of Hessian and advective terms,


    L ψ=(I+p⃗p⃗):∂^2ψ/∂p⃗∂p⃗+3p⃗.∂ψ/∂p⃗.

The terms are discretized using central differencing. The Laplacian
piece is computed as:


    (I:∂^2ψ/∂p⃗∂p⃗)_j,k=(∂^2ψ/∂ p_∥^2+∂^2ψ/∂ p_⊥^2)_j,k=X_j+1/2,k-X_j-1/2,k/Δ p_∥,j+p_⊥,k+1/2Y_j,k+1/2-p_⊥,k-1/2Y_j,k-1/2/p_⊥,kΔ p_⊥,k,

where,


    X_j+1/2,k=ψ_j+1,k-ψ_j,k/Δ p_∥,j+1/2,         Y_j,k+1/2=(ψ_j,k+1-ψ_j,k)/Δ p_⊥,k+1/2.

The remaining Hessian piece is computed as:


    (p⃗p⃗:∂^2ψ/∂p⃗∂p⃗)_j,k    =    [[ p_∥,jp_∥,j p_⊥,kp_∥,j; p_∥,jp_⊥,k p_⊥,kp_⊥,k ]]:[[    ∂^2ψ/∂ p_∥^2 ∂^2ψ/∂ p_∥∂ p_⊥; ∂^2ψ/∂ p_⊥∂ p_∥    ∂^2ψ/∂ p_⊥^2 ]]_j,k
        =     p_∥,jp_∥,jQ_j+1/2,k-Q_j-1/2,k/Δ p_∥,j+p_⊥,kp_⊥,kR_j,k+1/2-R_j,k-1/2/Δ p_⊥,k
        +     2p_⊥,kp_∥,jT_j,k+1/2-T_j,k-1/2/Δ p_⊥,k,

where


    Q_j+1/2,k=ψ_j+1,k-ψ_j,k/Δ p_∥,j+1/2,   R_j,k+1/2=ψ_j,k+1-ψ_j,k/Δ p_⊥,k+1/2,



    T_j,k+1/2=1/2(ψ_j+1/2,k+1-ψ_j-1/2,k+1/Δ p_∥,j+ψ_j+1/2,k-ψ_j-1/2,k/Δ p_∥,j).

The advective piece is computed as:


    (3p⃗.∂ψ/∂p⃗)_j,k=3(p_∥,jψ_j+1/2,k-ψ_j-1/2,k/Δ p_∥,j+p_⊥,kψ_j,k+1/2-ψ_j,k-1/2/Δ p_⊥,k).

Note the cell faced values of ψ are found by linear averaging
across cell centered values, for example ψ_j+1/2,k=0.5(ψ_j,k+ψ_j+1,k)
and ψ_j,k+1/2=0.5(ψ_j,k+ψ_j,k+1). 



 §.§ A2. Collisional coefficients. 


Once the potentials are determined, the friction coefficients are
evaluated using Eq. (<ref>). The components of K⃗ψ
at the cell center are defined as:


    (Kψ)_j,k    =    (I+p⃗p⃗)_j,k·(∂ψ/∂p⃗)_j,k
        =    [[ (1+p_∥,jp_∥,j)(ψ_j+1/2,k-ψ_j-1/2,k/Δ p_∥,j)+p_∥,jp_⊥,k(ψ_j,k+1/2-ψ_j,k-1/2/Δ p_⊥,k); (1+p_⊥,kp_⊥,k)(ψ_j,k+1/2-ψ_j,k-1/2/Δ p_⊥,k)+p_∥,jp_⊥,k(ψ_j+1/2,k-ψ_j-1/2,k/Δ p_∥,j) ]]

The cell face values of ψ are found by taking the average of
cell-centered values. A similar discretization approach is used when
evaluating the diffusion coefficient, Eq. (<ref>) .



 §.§ A3. Reformulated off-diagonal tensor diffusion terms (effective friction
coefficients). 


The off-diagonal diffusion coefficients are expressed as effective
friction coefficients of the form D_∥⊥∂lnf/∂ p_⊥
and D_⊥∥∂lnf/∂ p_∥, see
Eq (<ref>). The momentum-space derivatives of lnf at
cell centers are evaluated by averaging the cell vertex values, for
example:


    (∂lnf/∂ p_⊥)_j,k=1/4((∂lnf/∂ p_⊥)_j+1/2,k+1/2+(∂lnf/∂ p_⊥)_j-1/2,k+1/2+(∂lnf/∂ p_⊥)_j-1/2,k-1/2+(∂lnf/∂ p_⊥)_j+1/2,k-1/2),

where the cell vertex value is obtained by averaging over adjacent
face-centered values:


    (∂lnf/∂ p_⊥)_j+1/2,k+1/2=1/2(ln(|f_j,k+1|+ϵ_l)-ln(| f_j,k|+ϵ_l)/Δ p_⊥,k+1/2+ln(|f_j+1,k+1|+ϵ_l)-ln(|f_j+1,k|+ϵ_l)/Δ p_⊥,k+1/2),



    (∂lnf/∂ p_⊥)_j+1/2,k-1/2=1/2(ln(|f_j,k|+ϵ_l)-ln(| f_j,k-1|+ϵ_l)/Δ p_⊥,k-1/2+ln(|f_j+1,k|+ϵ_l)-ln(|f_j+1,k-1|+ϵ_l)/Δ p_⊥,k-1/2),



    (∂lnf/∂ p_⊥)_j-1/2,k+1/2=1/2(ln(|f_j-1,k+1|+ϵ_l)-ln(| f_j-1,k|+ϵ_l)/Δ p_⊥,k+1/2+ln(|f_j,k+1|+ϵ_l)-ln(|f_j,k|+ϵ_l)/Δ p_⊥,k+1/2),



    (∂lnf/∂ p_⊥)_j-1/2,k-1/2=1/2(ln(|f_j-1,k|+ϵ_l)-ln(| f_j-1,k-1|+ϵ_l)/Δ p_⊥,k-1/2+ln(|f_j,k|+ϵ_l)-ln(|f_j,k-1|+ϵ_l)/Δ p_⊥,k-1/2),

where ϵ_l=10^-30 is added to mollify singularities.
Once computed at the cell centers, the friction coefficients at the
cell faces are found by linear averaging. 



§ SOLUTION OF SINGULAR INTEGRALS IN RELATIVISTIC POTENTIALS






 §.§ B1. Solution of first singular integral


We seek a solution of the integral:

    I=∫_0^2πdϕ/√(r^2-1),

with:

    r=√((1+p^2)(1+(p')^2))-𝐩·𝐩'=√((1+p^2)(1+(p')^2))-p_∥p'_∥_a^2-p_⊥p'_⊥_b^2cosϕ=a^2-b^2cosΦ.

We consider the case of b^2>0. Note r^2-1=(r+1)(r-1). Since
r≥1, it follows that:

    a^2≥ b^2+1.


To begin, we consider the change of variable t=cosϕ. We consider
the following cases:

    t=cosϕ    ,    ϕ∈[0,π/2],ϕ∈[3π/2,2π] ; dϕ=-dt/√(1-t^2)
    
    t=-cosϕ    ,    ϕ∈[π/2,3π/2] ; dϕ=dt/√(1-t^2).

This gives:

    I=2∫_0^1dt/√((1-t^2)(a^2+1-b^2t)(a^2-1-b^2t))_I_1+2∫_0^1dt/√((1-t^2)(a^2+1+b^2t)(a^2-1+b^2t))_I_2.





  §.§.§ Solution of I_2 integral


We begin with the integral I_2. We follow Abramowitz & Stegun
<cit.>, and consider the polynomials:

    Q_1    =     1-t^2,
    
    Q_2    =     (a^2+1+b^2t)(a^2-1+b^2t).

These polynomials have real roots ±1, -a^2+1/b^2,
-a^2-1/b^2. Because of Eq. (<ref>),
it is apparent that the last two roots are ≤-1, and hence Q_1
and Q_2 do not have nested roots. In this case, one can consider
the transformation to the canonical forms of the elliptic integrals
by constructing the polynomial:

    Q_1-λ Q_2=-(1+λ b^4)t^2-2b^2a^2λ t+1-λ(a^4-1).

Seeking a zero discriminant for the quadratic form in t gives the
following value for λ:

    b^4a^4λ^2    =     (λ(a^4-1)-1)(1+λ b^4)⇒λ^2b^4-λ(a^4-b^4-1)+1=0
       ⇒   λ_±=(a^4-b^4-1)±√((a^4-b^4-1)^2-4b^4)/2b^4.

Note that these roots are real and semi-positive, since, by Eq. (<ref>):

    (a^4-b^4-1)≥2b^2.

Also, it is clear that 

    λ_+>λ_->0,

and that:

    λ_+λ_-=1/b^4.

Since the discriminant for Eq. (<ref>) vanishes,
it follows that the roots of Q_1-λ Q_2 are perfect squares
and are given by:

    t=-t_±  ;  t_±=λ_±b^2a^2/1+λ_±b^4.

Therefore:

    Q_1-λ_+Q_2    =     -(1+λ b^4)(t+t_+)^2,
    
    Q_1-λ_-Q_2    =     -(1+λ b^4)(t+t_-)^2.


At this point, it is useful to point out a few properties of the roots
t_± in Eq. (<ref>). Firstly, from Eq. (<ref>)
it follows that:

    t_+>t_-.

Secondly, from the polynomial in Eq. (<ref>)
and the properties of the quadratic equations, we can write:

    t_±^2=λ_±(a^4-1)-1/1+λ_±b^4,

which can be used to prove that: 

    t_-≤1

(needed for later) as follows:

    t_-^2≤1⇔λ_-(a^4-b^4-1)_≥2b^2<2⇔λ_-≤1/b^2,

which can be shown to be true when noting that:

    (a^4-b^4-1)^2-4b^4=(a^4-b^4-1-2b^2)(a^4-b^4-1+2b^2)≥(a^4-b^4-1-2b^2)^2.

The inequality follows from Eq. (<ref>). Finally, using
Eq. (<ref>) and the inequality above, we can
also readily prove that:

    t_+=λ_+b^2a^2/1+λ_+b^4=a^2/b^21/λ_-+1≥a^2/b^2+1≥1,

which will be important later.

Eqs (<ref>, <ref>) can be solved for
Q_1 and Q_2 as follows:

    Q_2    =     a_2+(t+t_+)^2-a_2-(t+t_-)^2,
    
    Q_1    =     a_1+(t+t_+)^2-a_1-(t+t_-)^2.

Here:

    a_2±=1+λ_±b^4/λ_+-λ_-  ;  a_1±=λ_∓(1+λ_±b^4)/λ_+-λ_-.

Note that:


  * a_2+>a_2- (from Eq. <ref>).

  * a_1+/a_1-=t_-/t_+<1 (from Eq. <ref>).

From the expressions of Q_1, Q_2, one can write:

    Q_1Q_2=(t+t_+)^4[a_2+-a_2-(t+t_-)^2/(t+t_+)^2][a_1+-a_1-(t+t_-)^2/(t+t_+)^2].

Following Ref. <cit.>, we postulate the change of variables:

    w=(t+t_-)/(t+t_+)⇒ dw=t_+-t_-/(t_++t)^2dt.

Hence:

    I_2=∫_0^1dt/√(Q_1Q_2)=1/t_+-t_-∫_w_0^w_1dw/√([a_2+-a_2-w^2][a_1+-a_1-w^2]).

Here:

    w_0=t_-/t_+<1  ;  w_1=1+t_-/1+t_+  ;  w_0<w_1<1.

The result in Eq. (<ref>) can be written as a canonical
elliptic integral by considering:

    a_1+/a_1-=t_-/t_+=w_0=e^2<1  ;  a_2+/a_2-=1+λ_+b^4/1+λ_-b^4=d^2>1,

to find:

    I_2=∫_0^1dt/√(Q_1Q_2)=1/(t_+-t_-)√(a_2-a_1-)∫_e^2^edw/√([d^2-w^2][e^2-w^2]).

Here, we have used the surprising property that:

    w_1^2=(1+t_-/1+t_+)^2=t_-/t_+=e^2⇒ w_1=e,

which can be demonstrated by using the definition of t_± (Eq.
<ref>) and t_±^2 (Eq. <ref>). 



  §.§.§ Solution of I_1 integral


The solution of the integral I_1 follows a similar development,
except now:

    Q_1    =     1-t^2,
    
    Q_2    =     (a^2+1-b^2t)(a^2-1-b^2t).

With these definitions, it can be shown that the discriminant of the
combination Q_1-λ Q_2 is exactly the same, and therefore
so are the solutions λ_±. However the roots in t now
have opposite signs:

    t=t_±  ;  t_±=λ_±b^2a^2/1+λ_±b^4,

and the factorization of Q_1,2 reads:

    Q_2    =     a_2+(t-t_+)^2-a_2-(t-t_-)^2,
    
    Q_1    =     a_1+(t-t_+)^2-a_1-(t-t_-)^2.

From the expressions of Q_1, Q_2, one can write:

    Q_1Q_2=(t-t_+)^4[a_2+-a_2-(t-t_-)^2/(t-t_+)^2][a_1+-a_1-(t-t_-)^2/(t-t_+)^2].

Following Ref. <cit.>, we postulate the change of variables:

    w=(t-t_-)/(t_+-t)⇒ dw=t_+-t_-/(t_+-t)^2dt.

When postulating this change of variables, we have taken into account
the fact that t≤1<t_+ (Eq. <ref>), and that
t_-<1 (Eq. <ref>). It follows that:

    I_1=∫_0^1dt/√(Q_1Q_2)=1/t_+-t_-∫_w_0^w_1dw/√([a_2+-a_2-w^2][a_1+-a_1-w^2]),

where:

    w_0=-t_-/t_+=-e^2<0  ;  w_1=1-t_-/t_+-1>0.

As before, one can readily prove that:

    w_1^2=(1-t_-/t_+-1)^2=t_-/t_+=e^2,

and therefore w_1=e. There results:

    I_1=∫_0^1dt/√(Q_1Q_2)=1/(t_+-t_-)√(a_2-a_1-)∫_-e^2^edw/√([d^2-w^2][e^2-w^2]).





  §.§.§ Solution of total integral I


When combining these two solutions, we find:

    I     =     2(I_1+I_2)=2/(t_+-t_-)√(a_2-a_1-)[∫_e^2^edw/√([d^2-w^2][e^2-w^2])+∫_-e^2^edw/√([d^2-w^2][e^2-w^2])]
        =    2/(t_+-t_-)√(a_2-a_1-)[∫_e^2^e+∫_0^e+∫_-e^2^0]=2/(t_+-t_-)√(a_2-a_1-)[∫_e^2^e+∫_0^e+∫_0^e^2]
        =    4/(t_+-t_-)√(a_2-a_1-)∫_0^edw/√([d^2-w^2][e^2-w^2]),

which can be written in terms of the complete elliptic integral of
the second kind as <cit.>:

    I=4K(m)/d(t_+-t_-)√(a_2-a_1-)=4K(m)/(t_+-t_-)√(a_2+a_1-),

where in the last step we have used the definition of d (Eq. <ref>),
and where:

    m=e^2/d^2.





 §.§ B2. Solution of second singular integral


In the previous section, we determined the root structure of the radicand
and removed the odd terms in the radicand. We employ this approach
and also use ideas from Ref. <cit.> to express the following
elliptic integral, 


    H=∫_o^2πr dϕ/√(r^2-1),

 in terms of Legendre's elliptic functions. Recall that


    r=a^2-b^2cosϕ.

The integral can thus be expressed as,


    H     =     2∫_0^1(a^2-b^2t)dt/√((1-t^2)(a^2+1-b^2t)(a^2-1-b^2t))+2∫_0^1(a^2+b^2t)dt/√((1-t^2)(a^2+1+b^2t)(a^2-1+b^2t)).

We can regroup as we know from B1 the solution when the numerator
is unity, 


    H=a^2I+2∫_0^1-b^2t dt/√((1-t^2)(a^2+1-b^2t)(a^2-1-b^2t))_H_1+2∫_0^1b^2t dt/√((1-t^2)(a^2+1+b^2t)(a^2-1+b^2t))_H_2.

Removing the odd terms in the radicand, we obtain, 


    H_1=∫_0^1-b^2t dt/√(Q_1Q_2)=-b^2/(t_+-t_-)√(a_2-a_1-)∫_-e^2^e(wt_++t_-)/(1+w)dw/√([d^2-w^2][e^2-w^2])



    =-b^2/(t_+-t_-)√(a_2-a_1-)∫_-e^2^eR_1(w)dw/√([d^2-w^2][e^2-w^2]),

and


    H_2=∫_0^1b^2t dt/√(Q_1Q_2)=b^2/(t_+-t_-)√(a_2-a_1-)∫_e^2^e(-wt_++t_-)/(-1+w)dw/√([d^2-w^2][e^2-w^2])



    =b^2/(t_+-t_-)√(a_2-a_1-)∫_e^2^eR_2(w)dw/√([d^2-w^2][e^2-w^2]).

The rational functions of w, R_1 and R_2, can be expressed
in terms of odd and even functions. This is because the odd term can
be simplified into elementary functions via trigonometric substitutions,
see Ref. <cit.>. 


    R_1(w)=w/1-w^2(t_+-t_-)+t_--w^2t_+/1-w^2



    R_2(w)=w/1-w^2(t_+-t_-)-t_--w^2t_+/1-w^2

However, in our case we observe that these terms cancel each other.
Examining the odd terms in R_1 and R_2, and adding together
their contribution to H, we find:


    H_1^odd+H_2^odd=-∫_-e^2^e+∫_e^2^e=-∫_-e^2^e^2=0,

as the odd function is asymmetric about the origin. The even term
can be further factorized into:


    t_--w^2t_+/1-w^2=t_+-t_+-t_-/1-w^2.

Putting the above expression into H_1 and H_2 , the contributions
of the even terms may be expressed as, 


    H_1+H_2=-b^2t_+/(t_+-t_-)√(a_2-a_1-)(∫_-e^2^edw/√([d^2-w^2][e^2-w^2])+∫_e^2^edw/√([d^2-w^2][e^2-w^2]))



    +b^2/√(a_2-a_1-)(∫_-e^2^edw/(1-w^2)√([d^2-w^2][e^2-w^2])+∫_e^2^edw/(1-w^2)√([d^2-w^2][e^2-w^2]))

Simplifying and regrouping, we obtain the expression for H, see Eq.
(<ref>) as,


    H=(a^2-b^2t_+)I+4b^2/√(a_2+a_1-)Π(e^2,m)

where Π is the complete elliptic integral of the third kind,
with e^2<1. This formula has been verified numerically. Also
in the above, we made use of the following step which was described
earlier in the previous section. It is as follows,


    ∫_-e^2^e+∫_e^2^e=∫_-e^2^0+2∫_0^e-∫_0^e^2

Then substituting t=w/e to get the final form (<ref>).
The complete elliptic integral of the third kind is given by,


    Π(e^2,m)=∫_0^1dt/(1-e^2t)√((1-t^2)(1-mt^2)).


unsrt


