A fully implicit, scalable, conservative nonlinear relativistic Fokker-Planck 0D-2P solver for runaway electron

TL;DR

This paper introduces a fully implicit, scalable, and conservative solver for the relativistic nonlinear Fokker-Planck equation, improving accuracy and efficiency in modeling runaway electrons in tokamak disruptions.

Contribution

It develops a novel implicit solver with exact conservation, positivity preservation, and scalable parallel algorithms for the relativistic Fokker-Planck equation including advanced multigrid techniques.

Findings

Accurately reproduces runaway tail dynamics under strong electric fields.

Ensures exact conservation and positivity in numerical solutions.

Demonstrates efficiency and scalability through numerical tests.

Abstract

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.

Tables

Table 1. Table 1: Parallel and algorithmic scaling tests: Thermalization of a random distribution in domain with p ∥ ∈ ( − 10 , 10 ) subscript 𝑝 parallel-to 10 10 p_{\parallel}\in(-10,10) , p ⟂ ∈ ( 0 , 10 ) subscript 𝑝 perpendicular-to 0 10 p_{\perp}\in(0,10) , and Δ ​ t = 1 . Δ 𝑡 1 \Delta t=1. The results are averaged over 2 time steps with nonlinear relative convergence tolerance of r N ​ L = 10 − 5 subscript 𝑟 𝑁 𝐿 superscript 10 5 r_{NL}=10^{-5} (considering more time steps is not useful, as the solution has already settled into a MJ distribution).

$N_{\parallel }$	$N_{⟂}$	$np$	NLI per $\Delta t$	WCT (sec) per $\Delta t$	WCT/NLI	$\Delta t/\Delta t_{\mathrm{explicit}}$	PLI
128	64	4	4.5	4.85	1.07	82	7.8
256	128	16	4	6.0	1.5	329	9.14
512	256	64	4	7.4	1.85	1318	9.475
1024	512	256	5	11.1	2.2	5270	9.088
2048	1024	1024	6.5	15.8	2.43	21083	8.969
4096	2048	4096	12	30.6	2.55	84331	9.0

Equations

$${\partial }_t{f}_e+{\partial }_{\vec{p}}\cdot \left[(\vec{E}+{\vec{F}}_S)f_e\right]=\sum _{\beta =i,e}C(f_{\beta },f_e),$$

$${\tau }_{ee}^{relativistic}=\frac{4\pi {ϵ}_0^2{m}_e^2{c}^3}{q_e^4{n}_e\ln {\Lambda }_{ee}},$$

$$C(f_{\beta },f_e)={\partial }_{\vec{p}}\cdot \left[{\overline{\overline{D}}}_{\beta }\cdot {\partial }_{\vec{p}}{f}_e-\frac{m_e}{m_{\beta }}\overrightarrow{F_{\beta }}{f}_e\right],$$

$${\overline{\overline{D}}}_e$$

$$\overrightarrow{F_e}$$

$$\overline{\overline{L}}\psi$$

$$\vec{K}\psi$$

$$\overline{\overline{P}}$$

$$[L+1]h_0$$

$$[L-3]h_1$$

$$[L-3]h_2$$

$$L{g}_0$$

$$L{g}_1$$

$$L\psi =\overline{\overline{P}}:\frac{{\partial }^2\psi }{\partial \vec{p}\partial \vec{p}}+3\vec{p}\cdot \frac{\partial \psi }{\partial \vec{p}}.$$

$$h_0$$

$$h_1$$

$$h_2$$

$$g_0$$

$$g_1$$

$$D_{i,\parallel \parallel }=\frac{Z_{\mathrm{eff}}}{2v}\frac{p_{\perp }^2}{p^2},D_{i,\parallel \perp }=D_{i,\perp \parallel }=-\frac{Z_{\mathrm{eff}}}{2v}\frac{p_{\perp }{p}_{\parallel }}{p^2},D_{i,\perp \perp }=\frac{Z_{\mathrm{eff}}}{2v}\frac{p_{\parallel }^2}{p^2},{\vec{F}}_i=\vec{0},$$

$$\frac{1}{v}\approx \frac{1}{\sqrt{v^2+v_{cut}^2}},$$

$$F_{S,\perp }=-S\frac{p_{\perp }}{\gamma }(1+p_{\perp }^2),F_{S,\parallel }=-S\frac{p_{\parallel }}{\gamma }{p}_{\perp }^2,$$

$$C(f_e,f_e)={\partial }_{\vec{p}}\cdot ({\overline{\overline{D}}}_e{\nabla }_p{f}_e-\overrightarrow{F_e}{f}_e)\approx {\delta }_{\vec{p}}\cdot ({\vec{R}}_D-{\vec{R}}_F)={\delta }_{\vec{p}}\cdot \vec{R},$$

$${\left({\delta }_{\vec{p}}\cdot \vec{R}\right)}_{j,k}=\left(\frac{R_{\parallel ,j+1/2,k}-R_{\parallel ,j-1/2,k}}{\Delta p_{\parallel ,j}}+\frac{p_{\perp ,k+1/2}{R}_{\perp ,j,k+1/2}-p_{\perp ,k-1/2}{R}_{\perp ,j,k-1/2}}{p_{\perp ,k}\Delta p_{\perp ,k}}\right).$$

$$R_{D,\parallel ,j+\frac{1}{2},k}={\left(D_{\parallel \parallel }{\partial }_{p_{\parallel }}{f}_e+D_{\parallel \perp }{\partial }_{p_{\perp }}{f}_e\right)}_{j+\frac{1}{2},k},R_{D,\perp ,j,k+\frac{1}{2}}={\left(D_{\perp \parallel }{\partial }_{p_{\parallel }}{f}_e+D_{\perp \perp }{\partial }_{p_{\perp }}{f}_e\right)}_{j,k+\frac{1}{2}},$$

$$R_{F,\parallel ,j+\frac{1}{2},k}=F_{\parallel ,j+\frac{1}{2},k}{f}_{e,j+\frac{1}{2},k},R_{F,\perp ,j,k+\frac{1}{2}}=F_{\perp ,j+\frac{1}{2},k}{f}_{e,j,k+\frac{1}{2}}.$$

$$⟨1,C(f_e,f_e){⟩}_p$$

$$⟨p_{\parallel },C(f_e,f_e){⟩}_p$$

$$⟨\gamma ,C(f_e,f_e){⟩}_p$$

$$⟨A,B{⟩}_p^D\approx 2\pi \sum _{j=1}^{N_{\parallel }}\sum _{k=1}^{N_{\perp }}{A}_{j,k}{B}_{j,k}{p}_{\perp ,k}\Delta p_{\parallel ,j}\Delta p_{\perp ,k},$$

$$\eta =1+{\eta }_0+{\eta }_1(p_{\parallel }-{\bar{p}}_{\parallel }),$$

$$C^D(f_e,f_e)={\delta }_{\vec{p}}\cdot (\overrightarrow{\eta R_D}-{\vec{R}}_F),$$

$$⟨p_{\parallel },C^D(f_e,f_e){⟩}_p^D$$

$$⟨\gamma ,C^D(f_e,f_e){⟩}_p^D$$

$$\left[\begin{array}[]{cc}\langle\gamma\delta_{\vec{{p}}}\cdot\vec{{R}}\rangle_{p}^{D}&\langle\gamma(p_{\parallel}-\bar{{p}}_{\parallel})\delta_{\vec{{p}}}\cdot\vec{{R}}_{D}\rangle_{p}^{D}\\ \langle p_{\parallel}\delta_{\vec{{p}}}\cdot\vec{{R}}_{D}\rangle_{p}^{D}&\langle p_{\parallel}(p_{\parallel}-\bar{{p}}_{\parallel})\delta_{\vec{{p}}}\cdot\vec{{R}}_{D}\rangle_{p}^{D}\end{array}\right]\left[\begin{array}[]{c}\eta_{0}\\ \eta_{1}\end{array}\right]=\left[\begin{array}[]{c}\langle\gamma\delta_{\vec{{p}}}\cdot(\vec{{R}}_{F}-\vec{{R}}_{D})\rangle_{p}^{D}\\ \langle p_{\parallel}\delta_{\vec{{p}}}\cdot(\vec{{R}}_{F}-\vec{{R}}_{D})\rangle_{p}^{D}\end{array}\right],$$

$${\delta }_t{f}_e^n=C^D(f_e^n,f_e^n)+\underset{\mathcal{E}(f_e^n)}{\underbrace{\tilde{C}(f_{i,}^n{f}_e^n)-{\delta }_{\vec{p}}.\left[({\vec{F}}_S+\vec{E})f_e^n\right]}},$$

$${\delta }_t{f}_e^n=\sum _{i=0,1,\mathrm{2...}}\frac{b_i{f}_e^{n-i}}{\Delta t},$$

Full text

Los Alamos National Laboratory, New Mexico, USA

1 Introduction

Relativistic Coulomb collisions are modeled using an extended version of the Landau-Fokker-Planck collision operator [1]. Similarly to its non-relativistic counterpart [2], the operator assumes small-angle collisions, is well-posed, and features strict conservation of total particle number, total momentum, and total energy [3]. 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 [4] and CQL3D [5, 6] 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 [7, 8]. 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 [9]. 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 [1] or in differential form [10]. 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 $\mathcal{{O}}(N)$ solver, where $N$ is the number of grid points, as the integral form produces an $\mathcal{{O}}(N^{2})$ scaling when symmetry is preserved in the integral operators (which is required to achieve strict conservation properties [11, 12]). To obtain a nearly optimal scaling ($\mathcal{{O}}(N^{\alpha}\log N)$, with $\alpha\gtrsim 1$) with the differential formulation, we propose a multigrid-preconditioned GMRES [13] solver for the potentials along with an adaptive treatment for evaluating potential far-field boundary conditions [which greatly decreases their computational complexity from $\mathcal{O}(N^{3/2})$ to $\mathcal{O}(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 [14]. 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-Jüttner distribution as a steady-state solution to our system). Earlier algorithmic approaches proposed for this system are either linearly implicit (e.g., NORSE [4]), or lack strict conservation properties (e.g., CQL3D [5, 6]). 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 $\mathcal{{O}}(N^{1.1}\log N).$ Our conservation and preconditioning strategies follow closely those proposed in Ref. [9].

The rest of the paper is organized as follows. In $\S$2, 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 $\S$3, 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 $\S$4, we briefly describe our fully implicit nonlinear solver using an Anderson Acceleration scheme. In $\S$5, we discuss the numerical results that demonstrate the correctness of our implementation, and finally in $\S$6 we list the conclusions and scope for future work.

2 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 $\beta$,

[TABLE]

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

[TABLE]

$\vec{{p}}$ is the momentum vector normalized with $m_{e}c$, $m_{e}$ is the electron mass, $c$ is the speed of light, $q_{e}$ is the electron charge, $\vec{{E}}$ is the electric field normalized with the critical value for runaway electron generation [8], $E_{c}=n_{e}q_{e}^{3}\ln\Lambda_{ee}/4\pi\epsilon_{0}^{2}m_{e}c^{2}$, $n_{e}$ is the electron number density, $\epsilon_{0}$ is the electrical permittivity, $\ln\Lambda_{ee}$ is the Coulomb logarithm, $\vec{{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,

[TABLE]

where $\overline{{\overline{{D}}}}_{\beta}$ represents the collisional diffusion tensor coefficients and $\vec{{F}}_{\beta}$ represents the collisional friction vector coefficients (computed based on the appropriate background species $f_{\beta}$). Though Eq. (1) 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_{\parallel},p_{\perp}$), with the subscripts $\parallel$ and $\perp$ referring to directions parallel and perpendicular to the magnetic field, respectively (see Fig. 1). 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}>>\gamma m_{e}$ with $\gamma=\sqrt{{1+p^{2}}}$ the Lorentz factor).

2.1 Electron-electron collisions

The collisional coefficients, $\overline{{\overline{{D}}}}_{\beta}$ and $\vec{{F}}_{\beta}$, for electrons in Eq. (2) are expressed in terms of the Braams-Karney potentials [10]. 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 $\mathcal{{O}}(N\log N)$. The collisional coefficients are given by [10]:

[TABLE]

where the operators $\overline{{\overline{{L}}}}$, $\vec{{K}}$, and $\overline{{\overline{{P}}}}$ are defined as:

[TABLE]

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

[TABLE]

and then the $g$ potentials by solving:

[TABLE]

Here, the operator $L$ is defined as:

[TABLE]

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. (5,6) [1]:

[TABLE]

where $r=\gamma\gamma^{\prime}-\vec{{p}}\cdot\vec{{p^{\prime}}}$. Note that the integral kernels of $h_{0}$ and $g_{0}$ are singular when $r\rightarrow 1\,(\vec{{p}}\rightarrow\vec{{p^{\prime}}}),$ which require a specialized numerical treatment in terms of elliptic integrals for accuracy and efficiency (see $\S$3.5 and App. B). Also, we have devised an efficient adaptive algorithm to fill ghost cells that prevents these boundary integrals from leading to an $\mathcal{{O}}(N^{3/2})$ scaling of the computational complexity (see also $\S$3.5).

2.2 Modeling external effects

We consider several external effects, including an imposed electric field, $\vec{{E}}=(E_{\parallel},0)$, ions, and synchrotron radiation, $\vec{{F}}_{S}$.

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

[TABLE]

where $p^{2}=p_{\perp}^{2}+p_{\parallel}^{2}$, $v$ is the velocity magnitude (normalized with $c$), and $Z_{\mathrm{{eff}}}=\sum n_{i}Z_{i}^{2}/\sum n_{i}Z_{i}$ is the effective ion-charge state ($n_{i}$ and $Z_{i}$ refers to ion densities and charges). For a quasi-neutral plasma, $\sum n_{i}Z_{i}=n_{e}.$ Note that the electron-ion collision operator becomes singular at the origin $v\rightarrow 0$. We mollify this singularity by reformulating the singular part, as:

[TABLE]

where $v_{cut}=p_{cut}/\sqrt{{1+p_{cut}^{2}}}$ is the velocity cut-off, with $p_{cut}=2\Delta p$. Note this approximation of the singular term in the cylindrical space introduces a finite but small amount of heating as $p\rightarrow 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 [15, 16]. The reaction term has finite friction coefficients, given by:

[TABLE]

where $S=\tau_{ee}^{relativistic}/\tau_{r}$ relates the time scale of the synchrotron-radiation damping, $\tau_{r}$, to that of electron-electron relativistic collisions, $\tau_{ee}^{relativistic}$.

3 Algorithm

3.1 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,

[TABLE]

where $\vec{{R}}_{D}$ and $\vec{{R}}_{F}$ are the diffusion and friction fluxes, respectively, and $\delta_{\vec{{p}}}\cdot$ denotes the discrete form of the divergence operator, which in cylindrical-momentum space is written as:

[TABLE]

Fluxes at cell faces are given by:

[TABLE]

[TABLE]

The potential operator $L$ (Eq. 7) 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. (5)-(6), we apply far-field Dirichlet boundary conditions using Eqs. (8) as discussed in $\S$3.5. The collisional coefficients, $\overline{{\overline{{D}}}}_{\beta}$ and $\vec{{F}}_{\beta}$, 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. (9), synchrotron damping effects, Eq. (10), 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.

3.2 Discrete conservation strategy for the e-e collision operator

The relativistic electron-electron collision operator conserves the total particle number, momentum, $\vec{{p}}=\gamma\vec{{v}}$, and energy $E=\gamma$, as the moments of the collision operator satisfy:

[TABLE]

where $\langle a,b\rangle_{p}=\int_{p}ab2\pi p_{\perp}\mathrm{{d}}p_{\parallel}\mathrm{{d}}p_{\perp}$. Discretely, these inner products may be approximated by a mid-point quadrature rule as:

[TABLE]

where the superscript $D$ refers to the discrete representation of the summation operator, and $\Delta p_{\parallel,j}$ and $\Delta p_{\perp,k}$ are the width and height of a rectangular cell located at $(j,k).$

In general, the relationships shown in Eqs. (13-15) will not be satisfied due to numerical errors. Discrete particle number conservation (Eq. 13) is trivially satisfied by setting the normal component of diffusion and friction fluxes to zero at the boundary. However, enforcing Eqs. (14,15) is more challenging. A recent study [9] 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

[TABLE]

where the magnitudes of $\eta_{0}$ and $\eta_{1}$ are expected to be of the order of truncation error, and $\bar{{p}}_{\parallel}=\langle f_{e},p_{\parallel}\rangle_{p}/\langle 1,f_{e}\rangle_{p}$ is the mean momentum. Thus, the discrete collision operator is of the form,

[TABLE]

Integrating over the cylindrical-momentum domain,

[TABLE]

we obtain a system of two equations,

[TABLE]

for unknowns $[\eta_{0},\>\eta_{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).

3.3 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, $\mathcal{{O}}(\tau_{ee}^{relativistic})$. For typical bulk temperatures $\Theta=T/m_{e}c^{2}\sim 10^{-4}$, this implies a time-scale separation of six orders of magnitude between thermal and relativistic time scales (as $\tau_{ee}^{thermal}=\Theta^{3/2}\tau_{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, $\mathcal{{E}}$, on electron evolution can be written as:

[TABLE]

where the superscript $D$ represents the appropriate discrete form defined in $\S$3.2, and $\tilde{{C}}$ represents the Lorentz operator (see Eq. 9). For a general implicit backward time discretization scheme at time step $n$, we have,

[TABLE]

where constants, $b_{i}$, satisfy $\sum_{i}b_{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. (18) with $\vec{{c}}=(1,p_{\parallel},\gamma$), and averaging over the momentum space, we obtain:

[TABLE]

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.

3.4 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 [17] to construct the associated fluxes. For the diagonal components of the tensor diffusion term, $\overline{{\overline{{D}}}}\cdot\nabla_{p}f_{e}\,|_{\parallel\parallel}$ and $\overline{{\overline{{D}}}}\cdot\nabla_{p}f_{e}\,|_{\perp\perp}$, we employ a standard second-order discretization:

[TABLE]

[TABLE]

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. [18]:

[TABLE]

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.

3.5 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. (8). For the $h_{1},h_{2},g_{1}$ relativistic potentials, we use a trapezoidal-rule numerical integration with 24 discrete points in the $\phi$ angle. However, the kernels in $g_{0},h_{0}$ become singular in the limit of $p^{\prime}\rightarrow p\implies r\rightarrow 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:

[TABLE]

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

However, even after these reductions, evaluating potentials at all ghost points in the boundary remains expensive. There are approximately $\mathcal{{O}}(N^{1/2})$ ghost-cell boundary points, each point requiring $\mathcal{{O}}(N)$ integrals when using Eqs. (8). This makes the potential boundary evaluations scale poorly with the number of mesh points, $N$ [i.e., $\mathcal{{O}}(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. 2. 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, $\phi_{S}$, and the actual value of the potential integral at the targeted points $\phi_{I}$:

[TABLE]

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:

[TABLE]

Figure 3 illustrates the adaptive knots generated with the adaptive spline algorithm for the $g_{0}$ potential for a Maxwell-Jüttner distribution of $\Theta=10^{-2}$ in a mesh of $N_{\parallel}=2048$ and $N_{\perp}=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. 21:

[TABLE]

Here, $n_{s}$ is the order of the spine. For instance, for a fourth-order spline, this result predicts $N_{p}\sim 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. $\S$5.2.1 show that $N_{p}\sim\mathcal{O}(N^{0.13})$ for a fourth-order spline.

4 Nonlinear solver

The spatial and temporal discretization techniques prescribed in $\S$3 lead to a coupled nonlinear system of equations, which requires an iterative nonlinear solver for the distribution function. We use an Anderson Acceleration scheme [19] to converge iteratively the system, which we briefly summarize next.

Given a fixed point map based Picard iteration,

[TABLE]

where the superscript $k$ denotes the iteration step, Anderson Acceleration scheme [20] accelerates the convergence of the Picard iteration by using the history of past nonlinear solutions via:

[TABLE]

where in this study $m_{k}=\min(5,k).$ The coefficients $\alpha_{i}^{k}$ are determined via an optimization procedure that minimizes,

[TABLE]

subject to $\sum_{i=0}^{m_{k}}\alpha_{i}^{k}=1$.

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

[TABLE]

with $P^{k}$ the preconditioner, $\delta f^{k}$ the nonlinear increment, and $\mathcal{{R}}^{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}=(\partial\mathcal{{R}}/\partial f_{e})^{k}$, then Eq. (24) 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. (5,6). These are inverted for each nonlinear iteration at flux-assembly time along with the computation of conservation constraints $\eta_{0}$ and $\eta_{1}$, see Eq. 17. The nonlinear elimination of the residuals associated with the potentials and conservation constraints follows from previous studies [21, 9], and enables a conservative, optimal $\mathcal{{O}}(N\log N)$ solver when the Poisson operators are inverted optimally and scalably. Here, the linear potential equations are solved using a multigrid-preconditioned GMRES [22] 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. (24) is obtained by Picard linearization of the potentials and subsequent discretization of the full system,

[TABLE]

where $\mathcal{{E}}$ is a linear operator representing the net external effects on electrons, see Eq. (18). 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}\delta f^{k}=-\mathcal{{R}}^{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,

[TABLE]

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.

5 Results

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

5.1 Verification

We first verify conservation properties of the equilibrium Maxwell-Jüttner 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 [3, 23, 4]. Finally, we verify our algorithm with recent calculations of runaway dynamics in the $t\rightarrow\infty$ limit [8, 15, 16]. Verification results were obtained using the second-order BDF2 time-stepping scheme except when explicitly stated (see $\S$3.3 for details on the time-stepping scheme).

5.1.1 Preservation of stationary and boosted Maxwell-Jüttner distributions.

The computational domain is uniform with $N_{\parallel}=256$ and $N_{\perp}=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. 4$(a)$ illustrates a static Maxwell-Jüttner (MJ):

[TABLE]

in log scale with normalized temperature $\Theta=T/m_{e}c^{2}=1$. In Eq. 25, $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. 4$(b)$ demonstrates the evolution of the relative errors in the number density, momentum and energy for 200 $\tau_{ee}$. The relative errors are measured as,

[TABLE]

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:

[TABLE]

where the subscript $b$ denotes the values in the boosted frame. Fig. 5 illustrates a boosted MJ equilibrium distribution with $\Theta_{b}=0.15$, in a frame boosted by $p_{b}=2$ and with $\gamma_{b}=\sqrt{{1+p_{b}^{2}}}.$ The relative errors are shown for $10\tau_{ee}^{relativistic}\sim 100\tau_{ee}^{thermal}$, demonstrating identical behavior as in the stationary MJ case.

5.1.2 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-Jüttner.

Figure 6$(a)$ illustrates the collisional relaxation of two MJ distributions boosted by 2 units in opposite directions. Fig. 6$(b)$ depicts the relative errors in number density, momentum, and energy. After an initial transient stage $t\in(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:

[TABLE]

where $\Theta=1$ and $J$ is a random number function with a uniform distribution in the range $[0,1]$, see Fig. 7. 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}$.

5.1.3 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 [3]. To measure conductivity, we apply a small electric field, $\hat{{E}}_{\parallel}=10^{-3}E_{D}$ where $\hat{{E}}_{\parallel}$ is the parallel electric field, $E_{D}=E_{c}/\Theta$ is the Dreicer field, and $E_{c}$ the Connor-Hastie critical electric field [8]. The electron distribution is initialized using the Maxwell-Jüttner distribution at various temperatures $\Theta$ and the normalized electrical conductivity is computed as:

[TABLE]

where $v_{\parallel}=p_{\parallel}/\gamma$. A small electric field deforms the Maxwellian slightly to produce a net electron flow in the positive $p_{\parallel}$ direction. To prevent numerical overflow, for $\Theta\apprle 10^{-3}$ the initial distribution is defined using a non-relativistic Maxwellian. Fig. 8$(a)$ depicts $\bar{\sigma}/Z_{\mathrm{{eff}}}$ at various temperatures. The numerical simulation results (circles) are in excellent agreement with the analytical results (lines) from Ref. [3]. 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, $\Theta=\langle f_{e}(\vec{{p}},t)p^{2}/2\gamma\rangle_{p}$, are larger (but close) to their initial value. The electrical conductivity, $\bar{{\sigma}}$, in Eq. (29) is computed using the quasi-steady-state temperature.

Figure 8$(b)$ illustrates the time evolution of electrical conductivity for various electric-field strengths in the non-relativistic limit. Results of NORSE [4] (circles) and Weng et al. [23] (squares) are also shown. The electron distribution is initialized with a Maxwellian corresponding to an initial temperature of $\Theta=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 $\hat{{E}}_{\parallel}=0.01E_{D},$ we have better agreement with NORSE than Weng et al.. Ref. [4] hypothesizes that the observed deviations between NORSE and Weng et al. in the small $\hat{{E}}/E_{D}$ limit may be due to numerical heating in Weng et al.. Our results also suggest the same.

5.1.4 Reproducing runaway-electron tail dynamics

To verify runaway dynamics with existing linear test-particle studies [8, 16], we performed two linearized numerical simulations where we keep the collisional cofficients fixed in time to those of a Maxwell-Jüttner distribution with $n_{e}=1$ and $\Theta=0.01$. In the first simulation, we applied an electric field 2.25 times the critical value, i.e $E_{\parallel}=2.25.$ This causes some electrons to overcome the frictional force and accelerate to high speeds. Fig. 9 demonstrates the evolution of the runaway tail at $42000\,\tau_{ee}^{thermal}$ collision times. The asymptotic slope of the runaway tail as predicted by Connor-Hastie [8] is:

[TABLE]

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. (10). Because we are interested in the steady-state as $t\rightarrow\infty$, time marching is performed efficiently with a BDF1 time stepping scheme (rather than BDF2). Fig. 10 shows the electron distribution function at $t=777$ (i.e., $\approx 777,000\,\tau_{ee}^{thermal}$) with a damping coefficient of $S=0.1$. Other parameters are $E_{\parallel}=2.25,\Theta=0.01,$ and $Z_{\mathrm{{eff}}}=1$. Ref. [15] 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}\approx 18$, to form a second maximum. The location of the second maximum is in good agreement with Ref. [16]. 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. [4].

5.2 Solver performance

5.2.1 Algorithmic and parallel scalability.

Table 2 lists weak parallel scalability results for the thermalization of a random electron distribution, see $f_{e}^{rand}$ in Eq. (28), with $\Theta=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 $\Delta 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. (21) and a fourth-order spline is used for piece-wise interpolation. The explicit time step is calculated as:

[TABLE]

Table 2 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 2 lists parallel and algorithmic scalability results for the case of collisional relaxation of two Maxwell-Jüttner 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. (27). We observe good parallel (WCT/NLI) and algorithmic scalability (NLI) up to 4096 processors with $\Delta t/\Delta t_{\mathrm{{explicit}}}\sim 460$ for the high resolution case of $4096\times 2048$.

Fig. 11$(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 2) and the boosted MJ relaxation (blue line, Table 2). We observe excellent parallel scalability in both cases. The expected scaling for parallel multigrid-based solvers is WCT/NLI $\sim\mathcal{O}(\log N)$, resulting in an overall algorithmic scaling of $\mathcal{{O}}(N\log{N})$. Instead, we find WCT/NLI $\sim\mathcal{O}(N^{0.1}\log{N})$. The additional factor of $N^{0.1}$ originates in the growth of the number of spline knots with $N$, as predicted by Eq. 22 and demonstrated in Fig. 11$(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. 22 (which assumed uniformly spaced knots), and much more efficient than the naive scaling of $\mathcal{O}(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 $\mathcal{O}(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.

5.2.2 Spatial and temporal accuracy

Figure 12$(\emph{{a}}$) illustrates the spatial accuracy of the proposed scheme with the boosted MJ collisional relaxation problem described in $\S$5.2.1. The ‘exact’ electron distribution $f_{e}^{exact}$ is obtained at $t=0.15$ (before steady-state is reached) with $N_{\parallel}=1024$, $N_{\perp}=512$ and $\Delta t=0.01$, while the three data points (solid dots) correspond to coarser grids of $128\times 64$, $256\times 128$ and $512\times 256$. The blue dashed line corresponds to a second-order error scaling. The $\ell_{2}-$norm of the error between the ‘exact’ and numerical electron distributions is computed as:

[TABLE]

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

Figure 12$(\emph{{b}}$) illustrates the temporal accuracy of the implementation. In this case, the ’exact’ electron distribution $f_{e}^{exact}$ is obtained in a $256\times 128$ grid using the BDF2 time advancement scheme with $\Delta t=5\times 10^{-4}$, see description in $\S$3.3. The four data points correspond to larger time steps of $\Delta t=10^{-4},2\times 10^{-4},4\times 10^{-4},$ and $10^{-3}$. The proposed implementation is confirmed to be second-order accurate in time.

6 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 $\mathcal{O}(N^{0.1})$. Positivity of the distribution function is ensured using a continuum-based reformulation approach [18] combined with robust positivity-preserving discretizations schemes [17]. 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 $\mathcal{{O}}(N^{1.1}\log{N})$. The $\log{N}$ 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 [8, 15, 16]. 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.

Appendix A Discretization of operators in potential equations and collisional

coefficients

A1. Linear potential equations.

The potential operator $L$ consists of Hessian and advective terms,

[TABLE]

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

[TABLE]

where,

[TABLE]

The remaining Hessian piece is computed as:

[TABLE]

where

[TABLE]

[TABLE]

The advective piece is computed as:

[TABLE]

Note the cell faced values of $\psi$ are found by linear averaging across cell centered values, for example $\psi_{j+1/2,k}=0.5(\psi_{j,k}+\psi_{j+1,k})$ and $\psi_{j,k+1/2}=0.5(\psi_{j,k}+\psi_{j,k+1})$.

A2. Collisional coefficients.

Once the potentials are determined, the friction coefficients are evaluated using Eq. (4). The components of $\vec{{K}}\psi$ at the cell center are defined as:

[TABLE]

The cell face values of $\psi$ are found by taking the average of cell-centered values. A similar discretization approach is used when evaluating the diffusion coefficient, Eq. (3) .

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_{\parallel\perp}\partial\ln{f}/\partial p_{\perp}$ and $D_{\perp\parallel}\partial\ln{f}/\partial p_{\parallel}$, see Eq (19). The momentum-space derivatives of $\ln{f}$ at cell centers are evaluated by averaging the cell vertex values, for example:

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\epsilon_{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.

Appendix B Solution of singular integrals in relativistic potentials

B1. Solution of first singular integral

We seek a solution of the integral:

[TABLE]

with:

[TABLE]

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

[TABLE]

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

[TABLE]

This gives:

[TABLE]

Solution of $I_{2}$ integral

We begin with the integral $I_{2}$. We follow Abramowitz & Stegun [24], and consider the polynomials:

[TABLE]

These polynomials have real roots $\pm 1$, $-\frac{a^{2}+1}{b^{2}}$, $-\frac{a^{2}-1}{b^{2}}$. Because of Eq. (33), it is apparent that the last two roots are $\leq-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:

[TABLE]

Seeking a zero discriminant for the quadratic form in $t$ gives the following value for $\lambda$:

[TABLE]

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

[TABLE]

Also, it is clear that

[TABLE]

and that:

[TABLE]

Since the discriminant for Eq. (34) vanishes, it follows that the roots of $Q_{1}-\lambda Q_{2}$ are perfect squares and are given by:

[TABLE]

Therefore:

[TABLE]

At this point, it is useful to point out a few properties of the roots $t_{\pm}$ in Eq. (38). Firstly, from Eq. (36) it follows that:

[TABLE]

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

[TABLE]

which can be used to prove that:

[TABLE]

(needed for later) as follows:

[TABLE]

which can be shown to be true when noting that:

[TABLE]

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

[TABLE]

which will be important later.

Eqs (39, 40) can be solved for $Q_{1}$ and $Q_{2}$ as follows:

[TABLE]

Here:

[TABLE]

Note that:

• •

$a_{2+}>a_{2-}$ (from Eq. 36).

• •

$a_{1+}/a_{1-}=t_{-}/t_{+}<1$ (from Eq. 41).

From the expressions of $Q_{1}$, $Q_{2}$, one can write:

[TABLE]

Following Ref. [24], we postulate the change of variables:

[TABLE]

Hence:

[TABLE]

Here:

[TABLE]

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

[TABLE]

to find:

[TABLE]

Here, we have used the surprising property that:

[TABLE]

which can be demonstrated by using the definition of $t_{\pm}$ (Eq. 38) and $t_{\pm}^{2}$ (Eq. 42).

Solution of $I_{1}$ integral

The solution of the integral $I_{1}$ follows a similar development, except now:

[TABLE]

With these definitions, it can be shown that the discriminant of the combination $Q_{1}-\lambda Q_{2}$ is exactly the same, and therefore so are the solutions $\lambda_{\pm}$. However the roots in $t$ now have opposite signs:

[TABLE]

and the factorization of $Q_{1,2}$ reads:

[TABLE]

From the expressions of $Q_{1}$, $Q_{2}$, one can write:

[TABLE]

Following Ref. [24], we postulate the change of variables:

[TABLE]

When postulating this change of variables, we have taken into account the fact that $t\leq 1<t_{+}$ (Eq. 45), and that $t_{-}<1$ (Eq. 43). It follows that:

[TABLE]

where:

[TABLE]

As before, one can readily prove that:

[TABLE]

and therefore $w_{1}=e$. There results:

[TABLE]

Solution of total integral $I$

When combining these two solutions, we find:

[TABLE]

which can be written in terms of the complete elliptic integral of the second kind as [24]:

[TABLE]

where in the last step we have used the definition of $d$ (Eq. 47), and where:

[TABLE]

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. [25] to express the following elliptic integral,

[TABLE]

in terms of Legendre’s elliptic functions. Recall that

[TABLE]

The integral can thus be expressed as,

[TABLE]

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

[TABLE]

Removing the odd terms in the radicand, we obtain,

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

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. [25].

[TABLE]

[TABLE]

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:

[TABLE]

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

[TABLE]

Putting the above expression into $H_{1}$ and $H_{2}$ , the contributions of the even terms may be expressed as,

[TABLE]

[TABLE]

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

[TABLE]

where $\Pi$ 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,

[TABLE]

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

[TABLE]