Adaptive time-stepping Monte Carlo integration of Coulomb collisions

TL;DR

This paper introduces a robust Monte Carlo tool for simulating Coulomb collisions in relativistic plasmas, utilizing adaptive time stepping to enhance efficiency and accuracy in diverse plasma physics scenarios.

Contribution

It presents a new adaptive time-stepping Monte Carlo method for Coulomb collisions that supports relativistic effects and can be integrated into existing plasma simulation tools.

Findings

Operator converges to analytical results

Adaptive time step improves computational efficiency

Careless implementation can cause errors

Abstract

We report an accessible and robust tool for evaluating the effects of Coulomb collisions on a test particle in a plasma that obeys Maxwell-J\"uttner statistics. The implementation is based on the Beliaev-Budker collision integral which allows both the test particle and the background plasma to be relativistic. The integration method supports adaptive time stepping, which is shown to greatly improve the computational efficiency. The Monte Carlo method is implemented for both the three-dimensional particle momentum space and the five-dimensional guiding center phase space. Detailed description is provided for both the physics and implementation of the operator. The focus is in adaptive integration of stochastic differential equations, which is an overlooked aspect among existing Monte Carlo implementations of Coulomb collision operators. We verify that our operator converges to known analytical results and demonstrate that careless implementation of the adaptive time step can lead to severely erroneous results. The operator is provided as a self-contained Fortran 95 module and can be included into existing orbit-following tools that trace either the full Larmor motion or the guiding center dynamics. The adaptive time-stepping algorithm is expected to be useful in situations where the collision frequencies vary greatly over the course of a simulation. Examples include the slowing-down of fusion products or other fast ions, and the Dreicer generation of runaway electrons as well as the generation of fast ions or electrons with ion or electron cyclotron resonance heating.

Equations

$$\frac{\partial f_a}{\partial t}=\sum _b{C}_{ab}[f_a,f_b].$$

$$C_{ab}[f_a,f_b]=\frac{{\Gamma }_{ab}}{2m_a}\frac{\partial }{\partial \mathbf{u}}\cdot {\int }_{{\mathbb{R}}^3}d\bar{\mathbf{u}}{\mathbf{U}}_{\mathrm{BB}}\cdot \left(\frac{{\bar{f}}_b}{m_a}\frac{\partial f_a}{\partial \mathbf{u}}-\frac{f_a}{m_b}\frac{\partial {\bar{f}}_b}{\partial \bar{\mathbf{u}}}\right),$$

$${\mathbf{U}}_{\mathrm{BB}}=\frac{r^2}{\bar{\gamma }\gamma w^3}\left(w^2\mathbf{I}-\mathbf{u}\mathbf{u}-\bar{\mathbf{u}}\bar{\mathbf{u}}+r(\mathbf{u}\bar{\mathbf{u}}+\bar{\mathbf{u}}\mathbf{u})\right),$$

$${\mathbf{U}}_{\mathrm{L}}=\frac{1}{|\mathbf{u}-\bar{\mathbf{u}}|}\left(\mathbf{I}-\frac{(\mathbf{u}-\bar{\mathbf{u}})(\mathbf{u}-\bar{\mathbf{u}})}{|\mathbf{u}-\bar{\mathbf{u}}{|}^2}\right).$$

$$\frac{\partial f_a}{\partial t}=-\frac{\partial }{\partial \mathbf{u}}\cdot \left({\mathbf{K}}_a{f}_a\right)+\frac{\partial }{\partial \mathbf{u}}\frac{\partial }{\partial \mathbf{u}}:\left({\mathbf{D}}_a{f}_a\right),$$

$${\mathbf{D}}_{ab}[f_b]$$

$${\mathbf{K}}_{ab}[f_b]$$

$$f_b(\mathbf{u})=\frac{n_b{e}^{-\sqrt{1+u^2}/{\Theta }_b}}{4\pi {\Theta }_b{\mathrm{K}}_2(1/{\Theta }_b)}.$$

$${\mathbf{D}}_{ab}$$

$${\mathbf{K}}_{ab}$$

$$K_{ab}$$

$$D_{ab,\parallel }$$

$$D_{ab,\perp }$$

$${\mu }_0$$

$${\mu }_1$$

$${\mu }_2$$

$$L_0$$

$$L_1$$

$$d\mathbf{u}(t)={\mathbf{K}}_a(u(t),t)dt+{\mathbf{\sigma }}_a(u(t),t)\cdot d\mathbf{W},$$

$${\mathbf{\sigma }}_a{\mathbf{\sigma }}_a^{⊺}=2{\mathbf{D}}_a.$$

$${\mathbf{\sigma }}_a=\sqrt{2D_{a,\parallel }(u(t))}\hat{\mathbf{u}}\hat{\mathbf{u}}+\sqrt{2D_{a,\perp }(u(t))}\left(\mathbf{I}-\hat{\mathbf{u}}\hat{\mathbf{u}}\right).$$

$$\frac{\partial \cdot }{\partial \mathbf{u}}=\{\mathbf{x},\cdot \},$$

$$\frac{\partial }{\partial \mathbf{u}}\cdot {\mathbf{D}}_{ab}=\frac{\partial D_{ab,\parallel }}{\partial u}\hat{\mathbf{u}},$$

$$\frac{\partial f_a}{\partial t}=-\frac{\partial }{\partial \mathbf{u}}\cdot \left({\mathbf{Q}}_a{f}_a\right)+\frac{\partial }{\partial \mathbf{u}}\cdot \left({\mathbf{D}}_a\cdot \frac{\partial f_a}{\partial \mathbf{u}}\right),$$

$${\mathbf{Q}}_a(\mathbf{u})$$

$$Q_{ab}(u)$$

$$\frac{\partial f_a}{\partial t}=-\{x_i,Q_{a,i}{f}_a\}+\{x_i,D_{a,ij}\{x_j,f_a\}\},$$

$$\frac{\partial F_a}{\partial t}=-\frac{1}{\mathcal{J}}\frac{\partial }{\partial \mathbf{Z}}\cdot \left(\mathcal{J}{\mathcal{Q}}_a{F}_a+\mathcal{J}{\mathcal{D}}_a\cdot \frac{\partial F_a}{\partial \mathbf{Z}}\right),$$

$${\mathcal{Q}}_{a,\alpha }$$

$${\mathcal{D}}_{a,\alpha \beta }$$

$${\mathbf{\Delta }}_{\alpha }=\{{\mathcal{T}}^{-1}\mathbf{X},Z_{\alpha }{\}}_{gc},$$

$$\frac{\partial F_a}{\partial t}=-\frac{1}{\mathcal{J}}\frac{\partial }{\partial \mathbf{Z}}\cdot \left(\mathcal{J}{\mathcal{K}}_a{F}_a\right)+\frac{1}{\mathcal{J}}\frac{\partial }{\partial \mathbf{Z}}\frac{\partial }{\partial \mathbf{Z}}:\left(\mathcal{J}{\mathcal{D}}_a{F}_a\right),$$

$${\mathcal{K}}_a={\mathcal{Q}}_a+\frac{1}{\mathcal{J}}\frac{\partial }{\partial \mathbf{Z}}\cdot \left(\mathcal{J}{\mathcal{D}}_a\right).$$

$${\mathcal{D}}_a=\sum _b{\mathcal{D}}_{ab,\mathbf{X}}({\mathbf{I}}_{\mathbf{X}}-\hat{\mathrm{b}}\hat{\mathrm{b}})+\sum _b{\mathcal{D}}_{ab,u}\hat{u}\hat{u}+\sum _b{\mathcal{D}}_{ab,\xi }\hat{\xi }\hat{\xi },$$

$${\mathcal{Q}}_{ab,u}=-\frac{{\Gamma }_{ab}{n}_b}{m_a^2{c}^3}\frac{1}{u^2}\left(\frac{m_a}{m_b}{\mu }_1\right).$$

Full text

Detailed description is provided for both the physics and implementation of the operator. The focus is in adaptive integration of stochastic differential equations, which is an overlooked aspect among existing Monte Carlo implementations of Coulomb collision operators. We verify that our operator converges to known analytical results and demonstrate that careless implementation of the adaptive time step can lead to severely erroneous results.

The operator is provided as a self-contained Fortran 95 module and can be included into existing orbit-following tools that trace either the full Larmor motion or the guiding center dynamics. The adaptive time-stepping algorithm is expected to be useful in situations where the collision frequencies vary greatly over the course of a simulation. Examples include the slowing-down of fusion products or other fast ions, and the Dreicer generation of runaway electrons as well as the generation of fast ions or electrons with ion or electron cyclotron resonance heating.

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: AMCC – (A)daptive (M)onte-(C)arlo (C)oulomb collisions

Licensing provisions: LGPL-2.1

Programming language: Fortran 95

Nature of problem: Test-particle tracing is a common feat within existing fusion applications. While efficient adaptive methods exist for integrating the incompressible Hamiltonian flow, the effects of Coulomb collisions are commonly implemented with far less sophisticated algorithms.

Solution method: The relativistic Fokker-Planck equation for test-particles in Maxwell-Jüttner background plasmas is converted into a stochastic differential equation. The stochastic differential equation is solved using adaptive Monte Carlo techniques. Methods to evaluate the effect of Coulomb collisions for both the three-dimensional particle momentum space and the five-dimensional reduced guiding center phase space are included.

*Additional comments including Restrictions and Unusual features:

*The package includes optionality to evaluate the relativistic Fokker-Planck coefficients, a feature useful for constructing accurate orbit averaged collision operators. The package also provides explicit one-step sympletic integrator for the relativistic Lorentz force that can be used for tracing test-particles in given electromagnetic backgrounds.

1 Introduction

Standard Runge-Kutta methods developed for ordinary differential equations (ODEs), are based on first converting the ODE into an integral equation, and then discretizing the integral. Stochastic differential equations (SDEs), often encountered when stochastic processes such as collisions are of interest, can equivalently be transformed into integral equations, with the exception that the resulting integral is non-Riemannian [1]. The discretization rules used for ODEs then no longer apply and, instead, either Itô or Stratonovich calculus must be adopted to obtain numerical methods for solving SDEs.

Perhaps the most straight-forward discretization method for SDEs is the Euler-Maruyama method. While simple to implement, the Euler-Maruyama method suffers from constant time-step requirement in the sense that one is not allowed to simply recompute the step if the estimated integration error turns out too large, as this might lead to a convergence to a wrong solution [2]. The reason is that while Euler-Maruyama has weak order of convergence 1.0, its strong order of convergence is only 0.5. For adaptive integration of SDEs, a method must have strong convergence of at least 1.0 to guarantee that it converges to a correct solution. A common application, where time steps are tested and rejected when needed, is test-particle tracing where the Hamiltonian motion of the particle is solved with adaptive Runge-Kutta methods. Therefore, to treat collisions with a manner consistent with the Hamiltonian motion, a higher order scheme than the Euler-Maruyama must be used.

Further, the test-particle collision frequencies in a plasma can vary greatly during the course of following the particle trajectory: energetic particles start as nearly collisionless and slowly drift toward the more collisional bulk population, while some particles from the bulk may accelerate to high energies. Therefore, adaptive time-stepping methods should be adopted also for the stochastic contribution to the particle motion, to reduce both the convergence errors and computational costs.

We thus see fit to introduce an adaptive time-stepping algorithm that would allow incorporating the Coulomb scattering into tracing of either the particle or guiding-center dynamics. Our Monte Carlo implementation evaluates the momentum change of a test-particle when it collides with a Maxwell-Jüttner background plasma, consisting of electrons and possibly multiple ion species. The operator is based on the Beliaev-Budker collision integral [3] and, therefore, is applicable even if either or both the test particle and the background plasma populations are relativistic. Advanced Monte Carlo Coulomb collision operators have been developed quite recently [4, 5], but the operator developed here differs from those in that it is relativistic, adaptive, and also applicable in guiding center dynamics.

The paper is organized as follows. In section 2, we begin by discussing how the Monte Carlo operator is obtained from the Beliaev-Budker collision integral by finding the related SDE. The corresponding guiding-center operator is presented in section 3. Discretization of the relativistic particle collision operator is first done with the Euler-Maruyama method, in section 4, where we also show how to compute the collision coefficients efficiently as these are also needed for the adaptive time-stepping. The adaptive scheme is described in detail in section 5 for both particle and guiding center picture. The collision operators are verified in section 6, where we also compare the adaptive method with the fixed time step scheme, and confirm that the adaptive method is both faster and more accurate than the commonly used Euler-Maruyama method.

2 Test-particle Fokker-Planck equation

In the limit of binary collisions, the collisional evolution of the particle distribution function of species $a$, interacting with species $b$, is determined by

[TABLE]

In a plasma, the collisions are dominated by the small angle scattering events, so that the collision operator can be written in the Landau-Fokker-Planck form [6]

[TABLE]

where $\mathbf{U}_{\mathrm{BB}}$ is the relativistic Beliaev-Budker tensor. The coordinate $\bm{u}=\bm{p}/m_{a}c$ is particle $a$ momentum normalized to the rest mass $m_{a}$ and speed of light $c$, while $\bm{\bar{u}}=\bm{\bar{p}}/m_{b}c$ denotes the same for particle of species $b$. The quantities with an overbar are evaluated at $\bm{\bar{u}}$. The operator also depends on the species charges, $q_{a}$ and $q_{b}$, through the coefficient $\Gamma_{ab}=q^{2}_{a}q^{2}_{b}\ln\Lambda_{ab}/(4\pi\varepsilon_{0}^{2})$ where $\varepsilon_{0}$ is the vacuum permittivity, and $\ln\Lambda$ is the Coulomb logarithm. The Coulomb logarithm describes the ratio of minimum and maximum impact parameters $\ln\Lambda=\ln r_{\mathrm{max}}/r_{\mathrm{min}}$ and as such it indicates by which factor small angle scattering dominates the large angle scattering. Natural choice for the maximum impact parameter is the Debye length , $\lambda_{D}=\sqrt{\varepsilon_{0}\sum_{b}\left(T_{b}/n_{b}q_{b}^{2}\right)}$, where $n$ is density and $T$ is temperature. Depending which one is smaller, the minimum impact parameter is determined either via classical electron radius, $r_{\mathrm{cl}}=q_{a}q_{b}/(4\pi\varepsilon_{0}m_{r}\tilde{v}^{2})$, where $m_{r}=m_{a}m_{b}/(m_{a}+m_{b})$ is the reduced mass and $\tilde{v}=\left<|\bm{v}_{a}-\bm{v}_{b}|\right>$ is the mean relative velocity, or by a quantum mechanical limit, $r_{\mathrm{qm}}=\hbar/(2m_{r}\tilde{v})$, where $\hbar$ is the reduced Planck constant. In fusion plasmas, $\ln\Lambda$ has values in the range of 10 – 20, but the uncertainties in estimating $r_{\mathrm{min}}$ means that the Coulomb logarithm is only accurate to within $1/\ln\Lambda$.

Returning to the Beliaev-Budker tensor, we find that it has a rather complicated expression [3]

[TABLE]

with $r=\gamma\bar{\gamma}-\bm{u}\cdot\bm{\bar{u}}$, $w=\sqrt{r^{2}-1}$, and $\gamma=\sqrt{1+u^{2}}$ is the Lorentz factor. In the non-relativistic limit, $c\rightarrow\infty$, the Beliaev-Budker tensor reduces to the better known Landau tensor [7]

[TABLE]

From now on, we will focus solely on the relativistic expressions which, of course, are valid also in the non-relativistic regime. In order to obtain the test particle collision operator, we first write the Beliaev-Budker collision integral, Eq. (2), in an explicit Fokker-Planck form

[TABLE]

where the vector $\bm{K}_{a}=\sum_{b}\bm{K}_{ab}[f_{b}]$ and the tensor $\mathbf{D}_{a}=\sum_{b}\mathbf{D}_{ab}[f_{b}]$ are summations of the species-wise coefficients $\bm{K}_{ab}[f_{b}]$ and $\mathbf{D}_{ab}[f_{b}]$, and functionals of the distribution $f_{b}$. The expressions for the species-wise coefficients are

[TABLE]

From now on, we will assume that the distributions $f_{b}$ are Maxwell-Jüttner distributions

[TABLE]

Here $\mathrm{K}_{\nu}(x)$ (not to be confused with the coefficient $\bm{K}_{ab}$) is the $\nu^{\mathrm{th}}$ order modified Bessel function of the second kind, and $\Theta_{b}=T_{b}/m_{b}c^{2}$ is the normalized temperature. In this case, the species-wise diffusion tensor $\mathbf{D}_{ab}$ and the force $\bm{K}_{ab}$ become isotropic [7]

[TABLE]

where $\bm{\hat{u}}\equiv\bm{u}/u$ is the unit vector parallel to $\bm{u}$. The coefficients $K_{ab}$, $D_{ab,\parallel}$, and $D_{ab,\perp}$ are defined in terms of three special functions $\mu_{0}(u;\Theta_{b})$, $\mu_{1}(u;\Theta_{b})$, and $\mu_{2}(u;\Theta_{b})$, and they are given by the following expressions [8]

[TABLE]

The coefficients are illustrated in Fig. 1 as a function of $u$.

The special functions $\mu_{0}(u;\Theta)$, $\mu_{1}(u;\Theta)$, and $\mu_{2}(u;\Theta)$ are given by

[TABLE]

where the functions $L_{0}(u;\Theta)$ and $L_{1}(u;\Theta)$ are

[TABLE]

Our notation here differs from the Ref. [8] as we have multiplied both the numerator and denominator in Eqs. (14) - (16) by $e^{1/\Theta}$ to avoid floating point errors when $\Theta$ is small.

The Monte Carlo simulation of a given test particle distribution is based on the connection between the Fokker-Planck equation and stochastic differential equations. If the distribution function $f_{a}(\bm{u},t)$ satisfies equation (5), then an individual sample particle from $f_{a}(\bm{u},t)$ obeys the following stochastic differential equation, known as the Langevin equation, of Itô kind [9]

[TABLE]

where the rank-2 tensor $\boldsymbol{\sigma}_{a}$ satisfies the condition

[TABLE]

Here $d\bm{W}$ is a differential of a vector-valued, uncorrelated standard Wiener processes $\bm{W}\sim{\cal N}(\bm{0},t\mathbf{I})$, with ${\cal N}$ being the standard multivariate normal distribution. Since the diffusion tensor is diagonal, the decomposition Eq. (20) is easy to accomplish, and we find

[TABLE]

3 Guiding-center test-particle operator

In many applications the rapid oscillation of a charged particle in plane perpendicular to a magnetic field is of little interest. If the magnetic field is slowly varying, one can resort to guiding center formalism which omits this gyro motion and depends only on the gyro-averaged quantities, thus reducing the 6D particle phase space into a 5D guiding center phase space.

Rigorous transformation of the particle phase-space Fokker-Planck equation (5) into guiding-center phase-space was first carried out in Ref. [10]. The transformation is based on the observation that the noncanonical particle phase-space Poisson bracket could be used to express the momentum-space derivatives according to

[TABLE]

where $\mathbf{x}$ is the spatial coordinate. This fact can be put to use by first observing that the isotropic test-particle diffusion tensor satisfies

[TABLE]

so that the particle phase-space collision operator can be written in a form

[TABLE]

where the modified friction coefficient is defined by

[TABLE]

Now, using the Poisson brackets, one obtains

[TABLE]

where the subscripts $i,j$ denote the Cartesian indices for the vector and tensor components, and summation over repeated indices is assumed.

The guiding-center transformation then follows via Lie-transform (see Ref. [10] for details), and the final form of the collision operator for the gyroangle-independent guiding-center distribution function $F$ becomes

[TABLE]

where $\mathbf{Z}$ are guiding center phase-space coordinates and $\mathcal{J}$ is the transformation Jacobian. The guiding center friction and diffusion coefficients are

[TABLE]

where ${\cal T}^{-1}$ is the guiding-center push-forward, $\langle\cdot\rangle$ denotes a gyroaverage, and $\bm{\Delta}_{\alpha}$ are the so-called projection vectors, defined according to

[TABLE]

with $\{{\cal F},{\cal G}\}_{gc}$ being the guiding center Poisson bracket. Here indices $\alpha,\beta$ denote the guiding center coordinates. Writing the guiding-center Fokker-Planck equation in a form similar to Eq. (5), we find

[TABLE]

where the drift coefficient is

[TABLE]

The guiding center formalism rests on the assumption that the magnetic moment is invariant. This invariance can be used to reduce the number of equations of motion to four by choosing the magnetic moment as one coordinate. However, the diffusion tensor is not diagonal in this case [11] but, as we saw in the last section, a diagonal basis is desired when considering the numerical implementation. Fortunately, there exists a suitable set of coordinates where the diffusion tensor is (almost) diagonal. This basis is $\mathbf{Z}=(\mathbf{X},u,\xi)$ where $\mathbf{X}$ is the guiding center location, $u$ is the magnitude of the normalized momentum, and pitch is $\xi=\hat{\boldsymbol{\mathrm{u}}}\cdot\hat{\boldsymbol{\mathrm{b}}}$, where $\hat{\boldsymbol{\mathrm{b}}}$ is the unit vector parallel to the magnetic field $\mathbf{B}(\mathrm{X})$. In these curvilinear coordinates the Jacobian is $\mathcal{J}=m_{a}Bu$, and the diffusion tensor becomes

[TABLE]

with diagonal matrix $\mathbf{I}_{\mathbf{X}}$ having non-zero elements only in coordinates $\mathbf{X}$. The friction coefficient has only one component

[TABLE]

The guiding center collision operator has few notable differences to the particle collision operator. First, the collisions now cause also spatial diffusion, with a diffusion coefficient [11]

[TABLE]

where $\Omega=q_{a}B/m_{a}$ is the gyrofrequency. In uniform magnetic field this corresponds to classical diffusion. Second, the momentum magnitude and direction now have separate coefficients for the diffusion,

[TABLE]

and for the drift

[TABLE]

where

[TABLE]

is the pitch collision frequency. Third, both $\nu_{ab}$ and ${\cal K}_{ab,u}$ diverge at $u=0$, which makes the particle collision operator more attractive when simulating thermal particles. Also the guiding-center specific coefficients $\mathcal{D}_{ab,\mathrm{X}}$, ${\cal K}_{ab,u}$, and $\nu_{ab}$ are illustrated in Fig. 1. Note that ${\cal K}_{ab,\mathrm{X}}=0$, and that we have implicitly assumed uniform magnetic field when writing down the coefficients above, see Ref. [11] for details.

The guiding center collision operator we have is in curvilinear coordinates so obtaining the corresponding Langevin equation is not as trivial as it was in the particle picture. Details on the derivation of the Monte-Carlo operator are found in A, and here we only show the result

[TABLE]

where $\boldsymbol{\Sigma}_{a}$ is again easily obtained from the decomposition $(1/2)\boldsymbol{\Sigma}_{a}\boldsymbol{\Sigma}_{a}^{\intercal}=\boldsymbol{\mathcal{D}}_{a}$. Written explicitly, the guiding center collision operator is a set of equations

[TABLE]

with $\bm{W}_{\mathrm{X}}$, $W_{u}$, and $W_{\xi}$ being independent Wiener processes.

4 Monte Carlo algorithm

The equation (19) does not have a known analytical solution and, therefore, numerical methods are required. Given initial condition $\bm{u}(t_{0})$, the numerical approximation for $\bm{u}(t)$ is obtained by discretizing the time coordinate as $t_{k+1}=t_{k}+\Delta t$, $k=0,\dots,n$, and evaluating $\bm{u}(t_{k+1})$ at each step using the solution of the previous step as an initial condition. The evaluation is commonly done with the Euler-Maruyama method which, when leaving out the species subscript $a$ and switching to Einstein notation for brevity, reads

[TABLE]

where the index $i=1,2,3$ represents the coordinates in Cartesian basis. In this basis, the matrix $\mathbf{G}$ reads

[TABLE]

where $\delta_{ij}$ is the Dirac delta. The discretized differentials of Wiener processes are drawn from a normal distribution $\Delta W_{i}=W_{i}(t_{k+1})-W_{i}(t_{k})\sim{\cal N}(0,\Delta t)$.

Implementing Euler-Maruyama method is straightforward, and all the struggle is in evaluating the coefficients $\mathbf{K}$ and $\mathbf{D}$ that, through the special functions $\mu_{0}$, $\mu_{1}$, and $\mu_{2}$, depend on the integrals $L_{0}$ and $L_{1}$. These integrals (Eqs. (17) and (18)) cannot be solved analytically and, therefore, we evaluate them via the adaptive Simpson’s method. The adaptive Simpson’s method divides the integration interval into subintervals until the difference in the resulting numerical approximation between the successive divisions is less than a given tolerance. However, most of the contribution to $L_{0}$ and $L_{1}$ comes from small values of $u$, and the adaptive Simpson’s method converges to an incorrect value if the upper limit for $u$ is large while $\Theta$ is small. This can be avoided by noting that both integrands equal to unity at $u=0$ and, as the integrands decay rapidly, the integrals have practically constant values beyond a certain point. Therefore, we can make the integration robust by defining a cut-off limit for $u$ as

[TABLE]

where the accuracy is controlled with parameter $\epsilon$, and then applying the adaptive Simpson’s method separately on intervals $[0,u_{c}]$ and $[u_{c},u]$ (when $u_{c}<u$).

However, evaluating the integrals this way at every integration time step would be very inefficient. Therefore, we calculate and tabulate the $L_{0}$ and $L_{1}$ values for a wide range of $u$ and $\Theta$ values before proceeding to solve Eq. (46). From the cut-off limit, Eq. (48), we can deduce that the integrals are sensitive to $u$ and $\Theta$ only near the curve $u=\sqrt{\Theta^{2}+2\Theta}$, but not elsewhere as Fig. 2 illustrates. This means that $L_{0}$ and $L_{1}$ do not have to be tabulated with high resolution and, therefore, the tabulated values can be used to efficiently compute the coefficients $\mathbf{K}$ and $\mathbf{D}$ at every time step.

5 Integration with an adaptive time step

The collision operator introduced in the last section, Eq. (46), was based on the Euler-Maruyama method and, therefore, is not applicable for an adaptive time step scheme. An adaptive integration is desirable since especially the pitch collision frequency has a strong dependency on test particle momentum (recall Fig. 1). For example, in a proton-electron plasma with $T=10$ keV, $n=10^{20}$ m*-3*, a thermal electron has a pitch collision frequency in the order of $\nu\approx 10^{3}$ $s^{-1}$, while, for a high-energy electron with $E_{\mathrm{kin}}=1$ MeV, the frequency is two orders of magnitude less $\nu\approx 10$ $s^{-1}$.

The simplest method which is suitable for the adaptive scheme is the Milstein method which has both weak and strong convergence of 1.0. When the evolution of an $N$-dimensional stochastic variable $\bm{Y}$ is given by the following Itô form

[TABLE]

the Euler-Maruyama method is acquired by simple discretization: $d\bm{Y}\rightarrow\Delta\bm{Y}$, $dt\rightarrow\Delta t$, and $d\bm{W}\rightarrow\Delta\bm{W}$. The Milstein method

[TABLE]

has a form similar to the Euler-Maruyama, with the difference being the additional last term, which contains the double Itô integral

[TABLE]

In the adaptive scheme, the integration time step is varied so that local truncation error stays within user-defined boundaries. For our purposes, we follow the scheme proposed in Ref. [12], which was later refined in Ref. [13], in which two different error controls are used. The first error control, defined in terms of the deterministic component $\mathbf{p}$, is

[TABLE]

where $\varepsilon_{\mathrm{abs,i}}$ is the tolerated error for the variable $Y_{i}$. The second error control is defined in terms of the stochastic component $\mathbf{g}$,

[TABLE]

Alternatively, the second error control can be defined in terms of both $\mathbf{p}$ and $\mathbf{g}$

[TABLE]

which is useful when $\mathbf{g}^{\prime}$ is zero or expensive to compute. Here $J$ is the Stratonovich integral, $J_{i}=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{t}\circ dW_{i}dt$, which has the value [13]

[TABLE]

where $\omega_{i}\sim\mathcal{N}(0,\Delta t)$. Approximating $|\omega_{i}|\approx\sqrt{\Delta t}$, we get $|J_{i}|=(1/2)\Delta t(|\Delta W_{i}|+\sqrt{\Delta t/3})$. One can see that the error controls given by Eqs. (53) and (54) are of the order $\mathcal{O}((\Delta t)^{3/2})$, while Eq. (52) has higher order $\mathcal{O}((\Delta t)^{2})$. The latter then seems unnecessary. However, this is not the case in the limit of weak diffusion, where $\mathbf{g}\approx 0$, and the deterministic component $\mathbf{p}$ dominates.

5.1 Particle operator

The Milstein method for the particle operator is obtained from by first substituting the tensor $\mathbf{G}$, Eq. (47), in the extra term in Eq. (50). The next task would be to discretize the double integral $\mathbf{I}$, but this turns out to be fatal to our approach. The diagonal components of $\mathbf{I}$ can be shown to be exactly $I_{ii}=(\Delta W_{i})^{2}-\Delta t$, but the non-diagonal ones have to be approximated [14, 15]. There are ways to do this but they come with a hefty price: the strong order of converge would be reduced to 0.5. Since we chose the Milstein method for the exact reason that its strong order is 1.0, this is obviously unacceptable. On the other hand, approximations that would maintain the strong order 1.0 are complicated and expensive to compute, reducing the efficiency gains of the adaptive integration, so we seek a better way to implement the operator.

The fact that $\mathbf{G}$ is diagonal in $(u_{\parallel},u_{\perp})$ basis provides us with a suitable path. With only diagonal components included, i.e. the noise being commutative, the Milstein method reads

[TABLE]

so that the collision operator becomes

[TABLE]

where $(\hat{\boldsymbol{\mathrm{u}}},\hat{\boldsymbol{\mathrm{\perp}}}_{1},\hat{\boldsymbol{\mathrm{\perp}}}_{2})$ form an orthogonal basis, and prime denotes partial derivative with respect to $u$. Equation (57) has formally reduced to the Euler-Maruyama method as the noise is additive, i.e., $\partial D_{\perp}/\partial u_{\perp}=0$, so this is a special case where even the Euler-Maruyama method has strong convergence of 1.0. The error terms, Eqs. (52) and (53), are now solely determined by Eq. (5.1):

[TABLE]

where we have chosen Eq. (53) as the second error estimate instead of Eq. (54). This is because $K^{\prime}$ can be zero (recall Fig. 1), in which case both error estimates, Eqs. (52) and (54), would yield zero values.

The Milstein method requires additional computation of $D_{\parallel}^{\prime}$, while $K^{\prime}$ is needed for the error estimates. Both can be evaluated analytically ($D_{\perp}^{\prime}$ is included for completeness sake):

[TABLE]

where the derivatives of the special functions are

[TABLE]

One can note that $\mu_{1}^{\prime}=(u/\gamma)\mu_{2}$.

5.2 Guiding center operator

With Milstein discretization, the guiding center collision operator, Eqs. (43) - (45), becomes

[TABLE]

Equation (69) was reduced to the Euler-Maruyama form as we have assumed uniform magnetic field, i.e., $(\partial/\partial\mathrm{X})(\cdot)=0$. In addition to these equations, there are also boundary conditions. The particle pitch is limited to the interval $\left[-1,1\right]$ which can be enforced with a reflecting boundary condition,

[TABLE]

applied if $|\xi(t_{k+1})|>1$. Similarly, $u$ cannot have negative values, so a reflecting boundary condition should be set at $u=0$. However, as both $\nu$ and $\mathcal{K}_{u}$ diverge at $u=0$ (recall Fig. 1), setting the reflecting boundary condition to a small positive value but still below the thermal momentum value $u\approx\sqrt{2\Theta}$, e.g. at $u=0.05\sqrt{2\Theta}$, ensures that the time step in the adaptive scheme does not become extremely small, and that the method is stable if a fixed time step is used.

The error estimates from Eqs. (52) - (54) are

[TABLE]

Note that we have two tolerances: $\varepsilon_{\mathrm{abs,u}}$ for the momentum and $\varepsilon_{\mathrm{abs,\xi}}$ for the pitch which we set $\varepsilon_{\mathrm{abs,\xi}}=\varepsilon_{\mathrm{tol}}$ as the pitch values are bounded in an interval. Strictly speaking, the error tolerance Eq. (74), should be defined in terms of $\mathcal{K}_{u}$, but using $\mathcal{Q}_{u}$ we only have to evaluate the derivative

[TABLE]

$\mathcal{Q}_{u}$ is the dominant term in $\mathcal{K}_{u}$ outside the diffusion dominated regime so this alteration is justifiable. For the $\xi$ error term in Eq. (5.2), we had to choose diffusion error estimate Eq. (54) instead because the estimate derived from Eq. (53)

[TABLE]

diverges when $|\xi|=1$. No separate error limit was set for the spatial coordinate $\mathrm{X}$.

5.3 Optimal time step and Brownian bridge

Now that the criteria for time step rejection are established, the next task is to choose an optimal time step to minimize the number of rejections. A good guess for the initial step is $\Delta t_{\mathrm{init}}=\epsilon_{\mathrm{tol}}^{3/2}/\nu$, but the presence of $|\Delta\bm{W}|$ in the error estimates complicates finding the optimal time step. The simplest scheme, known as the Brownian tree, is based on halving and doubling the current time step. However, this scheme is far from optimal and, therefore, we choose to implement the algorithm, described in detail in Ref [12], where the next time step depends on the value of the error estimates. In the regime of weak diffusion, $\varepsilon_{\mathrm{drift}}>\varepsilon_{\mathrm{diff}}$, we could treat the collision operator as an ODE, and choose the next time step as

[TABLE]

where $\beta<1$ is a safety factor for which we set value $\beta=0.9$. However, $\varepsilon_{\mathrm{drift}}>\varepsilon_{\mathrm{diff}}$ could also be due to extraordinary small $|\Delta W|$, not because drift dominates, so the next time step is chosen as

[TABLE]

where

[TABLE]

Here $\Delta W_{\mathrm{opt}}=\beta\varepsilon_{\mathrm{diff}}^{-1/3}|\Delta W|$ is the estimate for the “optimal” value of Wiener process. In the diffusion dominated regime, $\varepsilon_{\mathrm{diff}}>\varepsilon_{\mathrm{drift}}$, we again determine the next step iteratively from

[TABLE]

where $n$ is determined by the condition

[TABLE]

where $l_{\mathrm{max}}=2$ if the current step was rejected, $l_{\mathrm{max}}=4$ if $\Delta W/\sqrt{\Delta t}<2$, i.e., the current Wiener process value was not an outlier, and otherwise $l_{\mathrm{max}}=6$.

Whenever new Wiener processes are generated, these must always be stored – even when the time step they are associated with is rejected or they are generated for the sole purpose of determining the next time step. Only when integration has reached time $t$ can processes $\bm{W}(t^{\prime})$, $t^{\prime}<t$, be discarded. The reason for this is that the realized values condition the distribution of the Wiener processes, so discarding them can lead to a bias if the discard mechanism is systematic. A systematic mechanism can arise, e.g., from error estimate Eq. (53) as discarding Wiener processes would lead to over-representation of small values of $\Delta W$ when $\mathbf{g}$ and $\mathbf{g}^{\prime}$ are large.

The bias is avoided by introducing the conditioned probability distribution known as the Brownian bridge. First, let $\bm{W}_{t_{-}}$ and $\bm{W}_{t_{+}}$ be adjacent realized Wiener processes with $t_{-}<t_{+}$. When $t_{-}<t<t_{+}$, the process $\bm{W}_{t}$ is no longer normally distributed as $\mathcal{N}(\bm{W}_{t_{-}},(t-t_{-})\mathbf{I})$, but follows a different normal distribution where the mean is

[TABLE]

and the variance

[TABLE]

One can observe that the variance of $\bm{W}_{t}$ has its maxima at the center of the interval $[t_{-},t_{+}]$ while the expected values follow a straight line from $\bm{W}_{t_{-}}$ to $\bm{W}_{t_{+}}$. Only when $\bm{W}_{t_{+}}$ for $t_{+}>t$ does not exist, $\bm{W}_{t}\sim\mathcal{N}(\bm{W}_{t_{-}},(t-t_{-})\mathbf{I})$.

6 Verification and benchmark

To summarize, we have now developed following collision operators: the fixed step Euler-Maruyama method in the particle phase space, Eq. (46), the adaptive Milstein method in the particle phase space, Eqs. (57) – (59), and the adaptive Milstein method in the guiding center phase space, Eqs. (69) – (71). From now on, we refer to these as FEP, AMP, and AMG, respectively. These operators should yield equivalent results which preferably are the same as those obtained analytically – a topic we investigate here.

Our first task is to verify that a given test particle population relaxes to Maxwell-Jüttner distribution, Eq. (8). In equilibrium, the magnitude of the momentum is distributed as

[TABLE]

and the pitch is distributed as $\xi\sim\mathcal{U}(-1,1)$, where $\mathcal{U}$ is the uniform distribution. According to Eq. (85), the test particle mean momentum should converge to $0.56$ and variance to $3.62$ when considering a test case where both background and the test particle population consist of (relativistic) electrons with $\Theta=1\times 10^{-1}$. Likewise mean pitch should converge to [math] and variance to $1/3$. From Fig. 3 we see that this is the case: All operators converge to the equilibrium values at the same rate and, thus, are verified in this regard.

The above test did not verify the spatial collision operator, Eq. (69), which, in a uniform magnetic field, should correspond to the classical diffusion given by the coefficient

[TABLE]

where $\rho$ is the Larmor radius, $\nu$ is the pitch collision frequency, and $f(u)$ is the momentum distribution function. The spatial diffusion coefficient $D_{B}$ can be estimated with a Monte Carlo method as $D_{B}=\sum_{j=1}^{n}D_{j}/n$, where $D_{j}$ are test particle diffusion coefficients, and $n$ is either the number of test particles or, for a single particle, the number of time steps. It can be shown that the ratio $D_{j}/D_{X}$ obeys $\chi^{2}$ distribution with 1 degree of freedom [16]. With the AMG operator, $D_{j}$ is easy to calculate as $D_{j}=(\Delta X_{j})^{2}/2\Delta t_{j}$, where $\Delta X$ is the change in guiding center position, along some predefined direction, during time $\Delta t$. With the AMP or FEP operator, the collisions only affect particle momentum and, hence, do not lead to spatial diffusion, unless the particle collision operator is coupled with the Lorentz force. For numerical evaluation of the Lorentz force, we use the energy conserving scheme [17].

Using a test case where the test particles are thermalized, and evaluating diffusion coefficient with different operators, we obtain the test particle diffusion coefficient distributions shown in Fig. 4. We thus confirm that all operators yield results matching well to the analytical result calculated directly from Eq. (86).

We have not yet shown that the rate at which test population relaxes to equilibrium is correct, or that the operators are valid also in the regime of low diffusivity, where the deterministic parts in the collision operators dominate. To tackle both issues, we perform a slowing-down simulation using the same background plasma but this time with test particles having initially $u=5$. In the slowing-down simulation, fast test particles are simulated until they reach a certain momentum below which they can be considered part of the thermal population. Here the common figure of merit is the slowing-down time, i.e., the time it takes for a particle to thermalize.

There are ways to analytically estimate the slowing-down time [18], but here we derive a simple estimate from the theory of stochastic processes. The slowing-down process is analogous to the problem of finding the so-called first passage time. Consider a stochastic 1D process defined by the Langevin equation (49) where coefficients $p$ and $g$ are constant. If this process has an initial value $Y(t_{0})=Y_{0}$, then the time, $\tau$, it takes to reach value $Y_{0}+\alpha$ for the first time, obeys an inverse Gaussian distribution $\tau\sim IG(\alpha/p,\alpha^{2}/g^{2})$. Here the 1D process in question is the guiding center momentum equation (5.2), $u_{0}$ is the particle initial momentum and $u_{\mathrm{min}}$ is the momentum value below which a particle can be considered thermalized. Setting $u_{0}=5$ and $u_{min}=1$, we have $\alpha=-4$. Furthermore, setting the background temperature to $\Theta=1\times 10^{-2}$, the coefficients $Q$ and $D_{\parallel}$ are approximately constant on the interval $[u_{min},u_{0}]$ (recall Fig. 1). Figure 5 shows the slowing down distribution estimated this way overlapping with the distributions obtained numerically with the different collision operators. We can therefore conclude that the operators yield a correct relaxation rate and are valid also in the regime of low diffusionality.

Having verified the collision operators, it is time to benchmark them to see what can be gained by using the adaptive time step. To this end, we introduce one additional operator, FMG, which is the guiding center AMG operator but with a fixed time step. The test case consists of simulating $u=5$ electrons in $\Theta=1\times 10^{-2}$ plasma until they slow down below the energy corresponding to the background temperature. We measure how the mean slowing-down time converges when decreasing the time step in the fixed schemes or the error tolerance in the adaptive ones. Here we noticed that the operators FEP and AMP converged to 0.670 s while FMG and AMG converged to 0.673 s. This difference could originate from the guiding center transformation, but it is insignificant as the Coulomb logarithm is only accurate to within $1/\ln\Lambda$. The rate of convergence for all operators is shown in Fig. 6, where the error is plotted with respect to the elapsed cpu time. In all cases, the slope of the fitted curves is approximately -1, which confirms that the operators have a weak order of convergence 1.0 as expected. The benchmark shows that the adaptive method reduces the computational time by a factor of 10 in the particle picture, and by a factor of 3 in the guiding center picture. The guiding center operators, both fixed and adaptive, are more efficient compared to the corresponding particle operators even though the guiding center operator has two more variables to be solved for.

Now that we have verified and benchmarked the collision operators, one might wonder was it necessary to go trough all the trouble with the Milstein method and Brownian bridge? Would the same results be achieved with using the adaptive scheme with the Euler-Maruyama method and by omitting Brownian bridge when rejecting time steps? The answer is no as figure 7 clearly shows. Using the Euler-Maruyama method and omitting Brownian bridge leads to a distribution that is strongly peaked and slightly biased to lower $u$ values. When using the Milstein method but still omitting Brownian bridge, the bias is no longer present but the distribution remains peaked. The peak is exactly what we would expect from earlier discussion of a biasing mechanism: the small values of $\Delta W$ are over represented as $D_{\parallel}$ is large, which leads to drift term being too dominant. This in turn drives markers towards the peak where $K$ changes sign. When using the Euler-Maruyama method with Brownian bridge, the peak disappears but the distribution is biased, which confirms that the Euler-Maruyama method is unsuitable for adaptive time-stepping. The difference to the analytical result is not large but it could be more significant in more complex cases than our test case. Note that the error estimate Eq. (62) is for the Milstein method so it cannot be used when using the Euler-Maruyama method adaptively. Instead, the “extra” term that separates the Euler-Maruyama and Milstein methods becomes the error estimate: $\varepsilon_{\mathrm{diff}}=|D_{\parallel}^{\prime}[(\Delta W_{3})^{2}-\Delta t]|/2\varepsilon_{\mathrm{abs}}$. Therefore, using Euler-Maruyama adaptively still requires computation of $D_{\parallel}^{\prime}$, so no computational benefits are gained when using it instead of the Milstein method.

7 Summary and conclusion

We have developed a robust and computationally efficient collision operator for simulating test particle Coulomb collisions with a background obeying Maxwell-Jüttner statistics. It features an adaptive time-step integration scheme and is based on the Milstein method that has both weak and strong order of convergence of 1.0. The collision operator is fully relativistic and can operate either in the 3D particle momentum space or in the 5D guiding center phase space. The operator was verified by showing that it converges to the known analytical results. The adaptive scheme decreased computation time by a factor of 10 (particle momentum operator) or 3 (guiding center operator) in comparison to the fixed time step integration, when simulating slowing down of fast particles during which the collision frequency changes significantly. The collision operator developed here is implemented in the accompanied code package which is intended to serve as a library for other codes featuring Coulomb collisions.

Acknowledgements

This work is part of the EUROFUSION Enabling Research project ER15-CEA-09. The work was partially funded by the Academy of Finland projects No. 259675 and No. 298126, and the work of E. Hirvijoki was supported by US DOE Contract DE-AC02-09-CH11466.

Appendix A Langevin equation in curvilinear coordinates

The relation between the partial differential Fokker-Planck equation and the stochastic differential Langevin equation is widely known [9]. However, proofs of this relation usually consider a system in Cartesian coordinates, but the guiding center Fokker-Planck equation (32) is in curvilinear coordinates. Therefore, here we derive the relation in curvilinear coordinates to show that the guiding center Langevin equation indeed is Eq. (42).

Let $\mathbf{z}$ be a stochastic process given by the Langevin equation

[TABLE]

where we assume $\bm{W}$ are independent Wiener processes. Applying Itô’s Lemma on an arbitrary function $a(\mathbf{z})$, we get

[TABLE]

A Wiener process has the property $E[d\bm{W}]=0$ and, therefore, taking the expectation value on both sides and formally dividing by $dt$ yields

[TABLE]

When coordinates $\mathbf{z}$ form a curvilinear $n$-dimensional system, the differential volume element becomes $d\mathbf{z}=Jdz_{1}dz_{2}\cdots dz_{n}$, where $J$ is the Jacobian. Assuming that the values of $\mathbf{z}$ obey a probability distribution $f(\mathbf{z},t)$, the left-hand side of Eq. (89) becomes

[TABLE]

where the integration is over the phase space $\Omega$. Writing the first term on the right-hand side of Eq. (89) explicitly

[TABLE]

and integrating by parts results in

[TABLE]

where the first term vanishes when we assume $f(\mathbf{z},t)\rightarrow 0$ when $\mathbf{z}\rightarrow\infty$. Now, repeating the above manipulation twice also for the second term on the right-hand side, the equation (89) can be written as

[TABLE]

Since $a$ is an arbitrary function, the final step is to define $\mathbf{G}\equiv(1/2)\mathbf{g}\mathbf{g}^{\intercal}$ and divide by $J$ to obtain

[TABLE]

which is the Fokker-Planck equation in curvilinear coordinates.