UBER v1.0: a universal kinetic equation solver for radiation belts

Recent proceedings in radiation belt studies have proposed new requirements for numerical methods to solve the kinetic equations involved. In this article, we present a numerical solver that can solve the general form of the radiation belt Fokker–Planck equation and Boltzmann equation in arbitrarily provided coordinate systems and with user-specified boundary geometry, boundary conditions, and equation terms. The solver is based upon the mathematical theory of stochastic differential equations, whose computational accuracy and efficiency are greatly enhanced by specially designed adaptive algorithms and a variance reduction technique. The versatility and robustness of the solver are exhibited in four example problems. The solver applies to a wide spectrum of radiation belt modeling problems, including the ones featuring non-diffusive particle transport such as that arising from nonlinear wave–particle interactions.

1 Introduction

In the space plasma environment, radiation belts refer to torus-shaped regions surrounding Earth and other magnetized planets that are filled with highly energetic charged particles trapped in the planetary magnetic field. Since their discovery (Van Allen and Frank, 1959; Vernov et al., 1959), radiation belts have been the focus of intense research due to the innumerable unknowns concerning their extremely dynamic behavior and their damaging effects on spacecraft (e.g., Baker, 2000; Welling, 2010). During slowly changing conditions, radiation belt particles undergo three types of periodic motion: gyration about field lines, bounce along field lines, and drift about the planet. With each periodic motion there is a corresponding adiabatic invariant, defined through the Hamiltonian action integral (e.g., Goldstein, 1980, chap. 10), that is only violated when the conditions are undergoing change on timescales shorter than the period. A widely adopted method to study the dynamics of radiation belts is to solve a kinetic equation describing the evolution of particle phase space density. In quasi-linear theory, this kinetic equation is usually a Fokker–Planck equation that takes the general covariant form (Schulz, 1991)

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial Q^{\mathit{\alpha }}}}\left(G{D}^{\mathit{\alpha }\mathit{\beta }}{\displaystyle \frac{\partial \overline{f}}{\partial Q^{\mathit{\beta }}}}\right)-{\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial Q^{\mathit{\alpha }}}}\left(G{h}^{\mathit{\alpha }}\overline{f}\right)+S\overline{f}+v,\end{array}$$

where $\overline{f}$ is the phase-averaged phase space density, $G=\text{det}\left(\frac{\partial J^{\mathit{\alpha }}}{\partial Q^{\mathit{\beta }}}\right)$ is the Jacobian determinant for the transformation from canonical action variables J^α ($\mathit{\alpha }=\mathrm{1},\mathrm{2},\mathrm{3}$) to the generalized coordinates Q^α, and D^αβ, h^α, S, and v are coefficients of the equation. Summation on repeated Greek indices is implied throughout this paper. In different radiation belts, the number of terms emerging on the right-hand side of Eq. (1) and their respective physical backgrounds may be different. For the Earth's outer radiation belt, the second and the fourth terms are usually missing; the first term represents diffusion caused by wave–particle interactions, and the third term is often a loss characterized by the particle lifetime (e.g., Li and Hudson, 2019, and the reference therein). In the low-altitude inner radiation belt where wave–particle interactions are not as significant, the first and second terms are often provided by the diffusion and dynamic friction caused by inter-particle Coulomb collisions (e.g., Selesnick, 2012), and the fourth term may be a source from cosmic ray albedo neutron decay (CRAND; e.g., Selesnick, 2015; Li et al., 2017). For radiation belts of the gas giants, all terms could be present (e.g., Horne et al., 2008; Lorenzato et al., 2012). The first two terms may be attributable to both wave–particle interactions and inter-particle collisions and, in addition, synchrotron radiation, which is negligible in Earth's radiation belts, bleeds energy for the ultra-relativistic electrons, and thus also contributes to the second term (e.g., Bolton et al., 2002, 2004). The third term could represent the moon-sweeping loss, and the fourth term may come from moon volcanic activities as a plasma source (e.g., Nénon et al., 2017).

In some circumstances, the dependence of phase space density on certain phases φ^ι could be discerned, and the radiation belt kinetic equation takes the form

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+{\dot{\mathit{\phi }}}^{\mathit{\iota }}{\displaystyle \frac{\partial \overline{f}}{\partial {\mathit{\phi }}^{\mathit{\iota }}}}={\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial Q^{\mathit{\mu }}}}\left(G{D}^{\mathit{\mu }\mathit{\nu }}{\displaystyle \frac{\partial \overline{f}}{\partial Q^{\mathit{\nu }}}}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial Q^{\mathit{\mu }}}}\left(G{h}^{\mathit{\mu }}\overline{f}\right)+S\overline{f}+v,\end{array}\end{array}$$

where a dot over φ^ι indicates its time derivative. The phase space density $\overline{f}$ here is only averaged over the phases varying faster than φ^ι, and the mutually exclusive indices ι and μ together form the complete set of α. The most common situation is perhaps the dependence of $\overline{f}$ on the drift phase φ³, which in the Earth's outer radiation belt may be caused by the wave activity dependence on magnetic local time (e.g., Shprits et al., 2009) and in the inner belt by the longitudinal variation of drift shell altitude (e.g., Tu et al., 2010; Xiang et al., 2019). With the spatial derivative term on the left-hand side, Eq. (2) appears as a Boltzmann equation for $\overline{f}$, but by Hamiltonian mechanics, the conjugating term ${\dot{Q}}^{\mathit{\iota }}(\partial \overline{f}/\partial Q^{\mathit{\iota }})$ should have also appeared on the left-hand side. Its absence is due to the fact that, for particles in the radiation belt energy range, the drift-phase-dependent electric potential energy is usually negligible in the unperturbed particle Hamiltonian so that φ³ becomes a cyclic variable. Appendix A provides a more comprehensive explanation of this equation, and we will return to the general case in which φ³ is not cyclic in the “Conclusions and discussion” section.

Various numerical models have been built to solve a specific form of either Eq. (1) or (2) (e.g., Beutier et al., 1995; Selesnick et al., 2003; Tao et al., 2008; Albert et al., 2009; Subbotin et al., 2010; Tu et al., 2013; Wang et al., 2017; Xiang et al., 2020, to name a few). Though their underlying numerical schemes might not restrict, these existing models have in practice been implemented with hard-coded choices of coordinates, relatively simple boundary geometry, and a roughly fixed number of equation terms; therefore, each model is applicable to a specific set of problems. This situation could become quite inconvenient when adiabatic invariants of particle motion are used as coordinates of phase space to model radiation belt dynamics, as promoted by Schulz (1996). The reasons are twofold: first, due to their vast range of magnitude and dramatically varying resolution, adiabatic invariant coordinates often require some kind of rescaling and transformation (e.g., Zheng et al., 2014), specific to the problem, to be computationally efficient; second, boundary geometry becomes complicated in adiabatic invariant coordinates, which could be challenging for finite-difference methods and led Subbotin and Shprits (2012) to seek new coordinates from combinations of the adiabatic invariants. However, the use of adiabatic invariant coordinates is crucial for some compelling problems in the radiation belts: for example, the mechanisms of storm-time electron loss in which adiabatic modulations due to magnetic field configuration change must be separated from non-adiabatic processes (Kim and Chan, 1997; Turner and Ukhorskiy, 2020) and the relative significance of Earthward diffusion versus CRAND as a possible inner belt electron source with the drift shell splitting effect contributing (Cunningham et al., 2018). It is the purpose of this article to present a numerical code, named UBER (for “universal Boltzmann equation solver”), that solves Eqs. (1) and (2) in an arbitrarily user-specified coordinate system up to three dimensions, with great freedom in specifying boundary geometry and boundary conditions and with various combinations of equation terms. Therefore, it is expected that the solver can be applied to a wide spectrum of radiation belt modeling problems. More importantly, the freedom of specifying equation terms implies that, in an asymptotic manner, UBER can even solve the integro-differential kinetic equations arising from non-diffusive particle transport, such as that formulated in Artemyev et al. (2018) for nonlinear wave–particle interactions, and thereby provides a viable means to incorporate non-diffusive transport into global radiation belt modeling.

The underlying mathematical theory of the solver is stochastic differential equation (SDE) theory. The SDE method had been utilized by Tao et al. (2008), Selesnick et al. (2013), and Zheng et al. (2014) in their modeling of radiation belts. The method is grid-free and enjoys unparalleled advantages in dealing with cross-diffusion components and complicated boundary geometry (e.g., Zheng et al., 2016), but it is meanwhile notorious for low efficiency ascribed to its Monte Carlo nature. In this article, we also describe specially designed numerical techniques that have enhanced the computational speed of the SDE method by an order of magnitude, thus making the solver much more affordable for large-scale simulations. Four example problems with distinct physical backgrounds are provided in this article to demonstrate the abilities and versatility of the solver.

2 Mathematical theory

The kinetic equations given as Eqs. (1) and (2) above are parabolic partial differential equations (PDEs). Written in the Kolmogorov backward form (see below), a parabolic PDE corresponds to a multidimensional SDE that describes the motion of an Itô stochastic process whose certain functional expectation satisfies the PDE; the PDE can then be solved by calculating path integrals of the corresponding stochastic process (e.g., Freidlin, 1985; Øksendal, 1998).

Let us consider the following partial differential problem composed of a Kolmogorov backward equation and a general set of initial and boundary conditions:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}f=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}^{\mathit{\alpha }\mathit{\beta }}(t,\mathit{x}){\partial }_{\mathit{\alpha }}{\partial }_{\mathit{\beta }}f+b^{\mathit{\alpha }}(t,\mathit{x}){\partial }_{\mathit{\alpha }}f+c(t,\mathit{x})f\\ & +u(t,\mathit{x}),\end{array}\text{(4)}& {\displaystyle }f(\mathrm{0},\mathit{x})=g_{\mathrm{0}}\left(\mathit{x}\right),\mathit{x}\in \overline{\mathrm{\Omega }},\text{(5)}& {\displaystyle }f(t,\mathit{x})=g_{\mathrm{1}}(t,\mathit{x}),\mathit{x}\in {\partial }_{\mathrm{1}}\mathrm{\Omega },\text{(6)}& {\displaystyle }\widehat{\mathit{\gamma }}\left(\mathit{x}\right)\cdot \mathrm{\nabla }f-\mathit{\lambda }(t,\mathit{x})f=\mathrm{0},\mathit{x}\in \partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega },\end{array}$$

where ∂_t and ∂_α are shorthand for the partial differentials with respect to t and the αth coordinate, respectively. In Eqs. (4)–(6), $\overline{\mathrm{\Omega }}$ denotes the closure of the domain and ∂Ω its boundary. In particular, ∂₁Ω represents the boundary pieces of the first type (Dirichlet) boundary condition, and $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$ indicate the boundary pieces excluding those in ∂₁Ω, which are of the second (Neumann, λ≡0) or the third type (Robin, λ≠0) boundary conditions. The unit vector $\widehat{\mathit{\gamma }}$ points into $\overline{\mathrm{\Omega }}$ and is not tangent to the local boundary.

The mathematical theory of SDEs establishes a relation between Eqs. (3)–(6) and the Itô stochastic process, whose spatial positions are denoted by the random variable X_s in Ω that obeys the reflected SDE:

$$\begin{array}{}\text{(7)}& \mathrm{d}{\mathit{X}}_s=\mathit{b}(t-s,{\mathit{X}}_s)\mathrm{d}s+\mathit{\sigma }(t-s,{\mathit{X}}_s)\cdot \mathrm{d}{\mathit{W}}_s+\widehat{\mathit{\gamma }}\left({\mathit{X}}_s\right)\mathrm{d}{k}_s,\end{array}$$

where the dot product on the right-hand side is between a rank-2 tensor and a vector, and the parameter s runs from 0 to t so that the stochastic process retrogrades in time from t to 0. The first term on the right-hand side describes the ballistic part of its motion. The second term describes the stochastic part, with the coefficient tensor σ satisfying $\mathit{\sigma }\cdot {\mathit{\sigma }}^T=\mathbf{a}$ (whose components are a^αβ). Note that this condition does not uniquely determine σ, but all satisfying σ tensors are equivalent (Levi's theorem; Freidlin, 1985; Zheng et al., 2014). W_s is a vector Wiener process of the same dimensions as X_s, with each dimension an independent Gaussian stochastic variable that has zero mean and variance s. The third term describes reflection of the stochastic process in the direction given by $\widehat{\mathit{\gamma }}$ on the boundary $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$, and k_s is a monotonic stochastic variable that only increases when the stochastic process is on that boundary to force X_s to stay in $\overline{\mathrm{\Omega }}$. k_s can thus be considered a measure of the time that the stochastic process spent on $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$ and hence has the name local time. The Itô process stops either in $\overline{\mathrm{\Omega }}$ when s=t or on ∂₁Ω at $s=\mathit{\tau }<t$.

A formal solution of the problem in Eqs. (3)–(6) is given by the Feynman–Kac formula (e.g., Kac, 1949; Øksendal, 1998; Klebaner, 2005):

$$\begin{array}{}\text{(8)}& f(t,\mathit{x})=\mathbb{E}\left[{\mathcal{F}}^{t,\mathit{x}}\right[{\mathit{X}}_s\left]\right],\end{array}$$

in which 𝔼 is the expectation operator, and ℱ^t,x[X_s] is a functional of the stochastic path X_s started from t and x that has the expression

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\mathcal{F}}^{t,\mathit{x}}\left[{\mathit{X}}_s\right]\\ & ={\mathbb{I}}_{\mathit{\tau }\ge t}{g}_{\mathrm{0}}\left({\mathit{X}}_t\right)\exp \left[\underset{\mathrm{0}}{\overset{t}{\int }}c(t-s,{\mathit{X}}_s)\mathrm{d}s-\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }(t-s,{\mathit{X}}_s)\mathrm{d}{k}_s\right]\\ & +{\mathbb{I}}_{\mathit{\tau }<t}{g}_{\mathrm{1}}(t-\mathit{\tau },{\mathit{X}}_{\mathit{\tau }})\exp \left[\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}c(t-s,{\mathit{X}}_s)\mathrm{d}s-\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathit{\lambda }(t-s,{\mathit{X}}_s)\mathrm{d}{k}_s\right]\\ & +\underset{\mathrm{0}}{\overset{t\wedge \mathit{\tau }}{\int }}u(t-s,{\mathit{X}}_s)\exp \left[\underset{\mathrm{0}}{\overset{s}{\int }}c(t-r,{\mathit{X}}_r)\mathrm{d}r-\underset{\mathrm{0}}{\overset{s}{\int }}\mathit{\lambda }(t-r,{\mathit{X}}_r)\mathrm{d}{k}_r\right]\mathrm{d}s,\end{array}\end{array}$$

where the symbol 𝕀_τ≥t is equal to 1 when τ≥t, which means the stochastic process has stopped in $\overline{\mathrm{\Omega }}$ before it had a chance to reach ∂₁Ω, and zero otherwise; t∧τ means the smaller between the two. Physically, the functional ℱ^t,x[X_s] is a propagator of contribution carried along the stochastic path from either the initial condition or the first type boundary condition to the point of solution, and the exponential functions indicate how this contribution enhances or decays along this path.

To formally solve the Fokker–Planck equation (Eq. 1) by the Feynman–Kac formula (Eq. 8), it remains to transform the equation together with its proper initial and boundary conditions into the form of Eqs. (3)–(6). To this end, directly expanding Eq. (1) and collecting terms with the same differentiation order yields its Kolmogorov backward form:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\partial }_t\overline{f}=& D^{\mathit{\alpha }\mathit{\beta }}{\partial }_{\mathit{\alpha }}{\partial }_{\mathit{\beta }}\overline{f}+\left[({\partial }_{\mathit{\beta }}{D}^{\mathit{\alpha }\mathit{\beta }}+D^{\mathit{\alpha }\mathit{\beta }}{\partial }_{\mathit{\beta }}\ln G)-h^{\mathit{\alpha }}\right]{\partial }_{\mathit{\alpha }}\overline{f}\\ & +\left[S-({\partial }_{\mathit{\alpha }}{h}^{\mathit{\alpha }}+h^{\mathit{\alpha }}{\partial }_{\mathit{\alpha }}\ln G)\right]\overline{f}+v.\end{array}\end{array}$$

Comparing Eq. (10) with Eq. (3) and taking Q^α equivalent to x^α, we thus have the correspondences of coefficients:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}\mathbf{a}=\mathrm{2}\mathbf{D},\\ \mathit{b}=\mathrm{\nabla }\cdot \mathbf{D}-\mathit{h},\\ c=S-\mathrm{\nabla }\cdot \mathit{h},\\ u=v,\end{array}\right.\end{array}$$

where in curvilinear coordinates, the divergence operator on a tensor field Γ is

$$\begin{array}{}\text{(12)}& \mathrm{\nabla }\cdot \mathbf{\Gamma }={\partial }_{\mathit{\alpha }}{\mathrm{\Gamma }}^{\mathit{\alpha }\mathrm{\dots }}+{\mathrm{\Gamma }}^{\mathit{\alpha }\mathrm{\dots }}{\partial }_{\mathit{\alpha }}\ln G,\end{array}$$

in which the dots stand for all other indices irrelevant to the operation, and the terms ∂_αln G come from summation of the Christoffel symbols in a covariant derivative (e.g., Mathews and Walker, 1970, chap. 15). It is worth remarking that −h appears in the expression for b so that the Itô process travels against the advection velocity. This is indeed the case since it is time-backwards. Also, from the expression for c, divergence of the advection serves as a loss of phase space density.

Initial and boundary conditions to Eq. (1) are transformed as follows. For the initial condition and the first type boundary condition, values of $\overline{f}(t,\mathit{x})$ are specified just as in Eqs. (4) and (5). For a flux boundary condition of the form $\mathrm{\Phi }=g_{\mathrm{2}}(t,\mathit{x})\overline{f}$, we note that the outward flux Φ across a boundary is given by $(\widehat{\mathit{n}}\cdot \mathbf{D}\cdot \mathrm{\nabla }\overline{f}-\widehat{\mathit{n}}\cdot \mathit{h}\overline{f})$, with $\widehat{\mathit{n}}$ the unit inward normal vector of $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$. Therefore, the corresponding boundary condition is

$$\begin{array}{}\text{(13)}& \widehat{\mathit{n}}\cdot \mathbf{D}\cdot \mathrm{\nabla }\overline{f}-(\widehat{\mathit{n}}\cdot \mathit{h}+g_{\mathrm{2}})\overline{f}=\mathrm{0}.\end{array}$$

Comparing Eq. (13) with Eq. (6), we determine that

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}\widehat{\mathit{\gamma }}={\displaystyle \frac{\widehat{\mathit{n}}\cdot \mathbf{D}}{|\widehat{\mathit{n}}\cdot \mathbf{D}|}},\\ \mathit{\lambda }={\displaystyle \frac{\widehat{\mathit{n}}\cdot \mathit{h}+g_{\mathrm{2}}}{|\widehat{\mathit{n}}\cdot \mathbf{D}|}}.\end{array}\right.\end{array}$$

Although the SDE (Eq. 7) does not prevent σ, and hence D, from being zero, the expressions in Eq. (14) do become singular for vanishing D on $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$. In the region where D vanishes, Eq. (3) is no longer parabolic but degenerates to an advection equation (a first-order PDE), for which imposing a Neumann or Robin boundary condition is over-determinant. In this case, we invoke on the boundary minimal diffusion in the eigen-direction of $\widehat{\mathit{n}}$ so that $\widehat{\mathit{\gamma }}=\widehat{\mathit{n}}$, and let $g_{\mathrm{2}}\equiv -\widehat{\mathit{n}}\cdot \mathit{h}$ so that λ=0, which means the advective flow is free to cross the boundary. The situation that D is finite but $|\widehat{\mathit{n}}\cdot \mathbf{D}|$ vanishes is considered pathological to our problem.

Up to this point, we have transformed the Fokker–Planck equation (Eq. 1) and its initial and boundary conditions to the problem in Eqs. (3)–(6) and gathered all expressions in Eqs. (11) and (14) for the constructing components of the SDE (Eq. 7) as well as the functional (Eq. 9). In order to solve the Boltzmann equation (Eq. 2), it suffices for us to just transform the equation into the form of Eq. (1). To this end, we expand the phase space by concatenating the coordinates Q^μ and φ^ι so that $x^{\mathit{\alpha }}=\mathit{\{}{Q}^{\mathit{\mu }},{\mathit{\phi }}^{\mathit{\iota }}\mathit{\}}$ (recall that $\mathit{\alpha }=\mathit{\{}\mathit{\mu },\mathit{\iota }\mathit{\}}$) and introduce the new coefficients ${\stackrel{\mathrm{̃}}{D}}^{\mathit{\alpha }\mathit{\beta }}$, ${\stackrel{\mathrm{̃}}{h}}^{\mathit{\alpha }}$, $\stackrel{\mathrm{̃}}{S}$, and $\stackrel{\mathrm{̃}}{v}$ that satisfy the following conditions:

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{l}{\stackrel{\mathrm{̃}}{D}}^{\mathit{\mu }\mathit{\nu }}=D^{\mathit{\mu }\mathit{\nu }},{\stackrel{\mathrm{̃}}{D}}^{\mathit{\alpha }\mathit{\iota }}=\mathrm{0},\\ {\stackrel{\mathrm{̃}}{h}}^{\mathit{\alpha }}=\mathit{\{}{h}^{\mathit{\mu }},{\dot{\mathit{\phi }}}^{\mathit{\iota }}\mathit{\}},\\ \stackrel{\mathrm{̃}}{S}=S+{\partial }_{\mathit{\iota }}{\dot{\mathit{\phi }}}^{\mathit{\iota }},\\ \stackrel{\mathrm{̃}}{v}=v.\end{array}\right.\end{array}$$

It can be verified that Eq. (1) in the new x^α coordinates with the new coefficients given by Eq. (15) transforms into Eq. (2) after replacing x^α by Q^μ and φ^ι. The transformation (Eq. 15) essentially treats φ^ι as new dimensions of the stochastic motion, except that the stochastic part of the motion in these dimensions is identically zero. A new type of boundary condition might emerge for problems involving Eq. (2), which is the periodic boundary condition for the phases φ^ι. From the viewpoint of stochastic motion, though, such periodicity is not really a boundary but rather a topology of Ω. The treatment of periodic boundary condition will be exemplified in the third problem in Sect. 4 below.

Table 1 User input items to the UBER code.

A set of the boundary-related items for each piece of boundary.

To summarize this section, the above mathematical theory allows us to fully define a PDE problem involving Eq. (1) or (2) in an arbitrary coordinate system given the input functions and equations as listed in Table 1, which can be either analytical or numerical in the UBER code. The equation terms may be freely turned off by setting their corresponding coefficients to zero. The number of boundary pieces is totally up to choice, which can even be zero to put the boundary at infinity. The boundary geometry may be time-variable for boundary pieces in ∂₁Ω in the UBER code, but it must be fixed for those in $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$. Solutions of this problem are obtained once we find a way to evaluate the functional in Eq. (9) for a realization of a stochastic path and to estimate the expectation of the functional. These numerical techniques are the subject of the next section.

3 Numerical techniques

We give an outline of the algorithms used by the UBER code in this section, with an emphasis on the techniques that improve both its accuracy and efficiency. Lower-level numerical techniques, such as the generation of pseudo-random variables, linear algebraic operations, and parallelized computation, are based on the works presented in Zheng (2015). The general idea for numerically implementing the SDE method is as follows: (i) for a given spatiotemporal position (t,x) for which an equation solution is wanted, a number of stochastic paths starting from this common position are simulated; (ii) for each stochastic path, its functional value is evaluated by the path integrals as in Eq. (9); and (iii) from these sampled functional values, their expectation is estimated, and this gives the solution at (t,x). Therefore, the SDE method is essentially a Monte Carlo method. It does not rely on a computational grid and is able to solve the problem locally. However, on many occasions it is still worth obtaining global solutions on a grid so that the solutions at time stamp T_i may be used as the initial condition for the solutions at T_i+1, analogous to the idea of layer methods (e.g., Tao et al., 2009). In this way, the stochastic processes need only to be simulated for a short duration of $t=T_{i+\mathrm{1}}-T_i$ to obtain the new solutions, for which the calculation of functional expectation would converge much faster than those simulated for the full length $t=T_{i+\mathrm{1}}$. The only operation on this grid would be interpolation and possibly extrapolation; therefore, unlike in layer methods, the grid can be irregular and thus allows for the use of sophisticated interpolation libraries on irregular grids, although in the present version only a nonuniform Cartesian grid (in the given coordinate system) is implemented with user-specified nodes.

Integration of the SDE (Eq. 7) employs the Euler–Maruyama scheme that is order 1 for weak convergence problems such as ours, meaning that when only the statistical distribution of stochastic paths matters but not the individual path, the expectation of the schematic error is proportional to the first power of the time step size (Kloeden and Platen, 1992). To further reduce the schematic error, an adaptive time step size is used in UBER. It can be shown that (e.g., Zheng, 2015) the root mean square (rms) distance an Itô stochastic process travels in infinitesimal time ds is

$$\begin{array}{}\text{(16)}& \mathrm{d}{\overline{X}}_s=\sqrt{\text{tr}\left(\mathbf{a}\right)\mathrm{d}s}.\end{array}$$

Numerically, the first-order contribution from $\left|\mathit{b}\mathrm{\Delta }s\right|$ cannot be neglected due to the finite Δs. Therefore, we prescribe a desired rms spatial step size $\mathrm{\Delta }{\overline{X}}_s$, which is sufficiently small compared to the size of Ω and any scale length of the equation coefficients, and then choose the smaller Δs inferred from either Eq. (16) or $\mathrm{\Delta }{\overline{X}}_s=\left|\mathit{b}\mathrm{\Delta }s\right|$ at every step of integration as the adaptive step size. This scheme evidently reduces to a simple adaptive Euler scheme for integrating ordinary differential equations when a approaches zero.

Oblique reflection of the stochastic process on $\partial \mathrm{\Omega }\backslash {\partial }_{\mathrm{1}}\mathrm{\Omega }$ and the calculation of dk_s follow the projected half-space algorithm presented in Gobet (2001), which is also order 1 in the weak convergence sense. The idea is that, for an exact half-space boundary, k_s can be proven to share the same probabilistic distribution with a composite stochastic variable involving W_s, coefficients of the SDE, the normal vector $\widehat{\mathit{n}}$, and an independent exponential random variable with parameter (2s)⁻¹ (Lépingle, 1995); therefore, k_s can be explicitly calculated by these known quantities. For general smooth boundary geometry, an additional contribution to dk_s may also come from possible projection along the $\widehat{\mathit{\gamma }}$ direction needed to keep the stochastic process within domain. With dk_s obtained and the SDE (Eq. 7) integrated, the functional (Eq. 9) can be readily evaluated by an ordinary numerical integration technique implemented along the realized stochastic path.

Expectation of the functionals can be estimated, in principle, from an arithmetic mean of a number N of sampled stochastic path integrals. The error of this estimation, $\mathit{ϵ}=\left|\mathbb{E}\right[{\mathcal{F}}^{t,\mathit{x}}]-\langle {\stackrel{\mathrm{̃}}{\mathcal{F}}}^{t,\mathit{x}}\rangle |$, where a tilde is used to indicate a numerical realization in this section and 〈⋯〉 indicates averaging over samples, can be estimated by dividing the simulation of stochastic processes into batches (Zheng, 2015). Although the probabilistic distribution of individual ${\stackrel{\mathrm{̃}}{\mathcal{F}}}^{t,\mathit{x}}$ is generally far from normal and largely unknown, that of the batch-wise mean of ${\stackrel{\mathrm{̃}}{\mathcal{F}}}^{t,\mathit{x}}$ approaches a Gaussian for a large enough sample number per batch due to the central limit theorem, and a confidence interval can thereby be calculated for the batch-wise means using the Student's t distribution (e.g., Kloeden and Platen, 1992). We use this confidence interval as an approximation to ϵ. In this way, UBER adaptively stops simulating more batches of stochastic processes when the estimated error meets a prescribed tolerance.

In typical radiation belt problems, the functional values from various stochastic paths may differ by orders of magnitude, and hence their contributions to the arithmetic mean also differ by orders of magnitude, whereas their computational efforts are of the same order. Therefore, straightforward calculation of their arithmetic mean could result in extremely slow convergence with N and squander computational power. To reduce statistical variance in this procedure, a process-splitting technique is developed based on the idea of importance sampling, i.e., to make “denser” sampling in more important “regions”. In conventional Monte Carlo methods, the “region” is an “area” in a parameter space, and importance sampling effectively splits one sample point therein that would have made a huge contribution to the calculation into many sample points nearby, while weights of these samples are reduced accordingly to keep the probabilistic distribution of samples unbiased (e.g., Press et al., 1992). But unlike conventional Monte Carlo methods, the samples in the SDE method are paths which belong to a functional space. To still implement this idea, we split the stochastic path when it is projected to contribute a large functional value.

Figure 1 Schematic illustration of process splittings in a t⊗ℝ² space. A stochastic process travels backward in time from point P and splits into two at point Q, where its projected functional value is found to be sufficiently large (see text for exact meaning). One child process splits again at point R, where its projected functional value is found to be even larger. The independent child processes would eventually stop either in $\overline{\mathrm{\Omega }}$ as at points A and B or on ∂₁Ω as at point C.

Figure 1 gives an illustration of this technique in a t⊗ℝ² space. As a stochastic path being integrated from point P, the functional value of the entire path (from s=0 to s=t) is continuously predicted based on the partial path that has been realized. This projected functional value is compared to the value of some quantile (e.g., the 80th percentile) statistically derived from all previously completed stochastic paths starting from the same position. When at some place Q the projected functional value falls above this quantile, the stochastic process is deemed to make a significant contribution to the arithmetic mean. It is then split into a number of child processes at Q, and each child process traces down an independent path thereafter. These child paths, together with their common parent path segment PQ, hence constitute “nearby samples” in the functional space. This procedure can be further iterated if the projected functional value later falls into an even higher quantile (e.g., the 90th percentile), as shown at R. After all procedures finish, the eventual result is a tree structure of stochastic paths rooted at P. For the illustration in Fig. 1, the actual functional value of the path PQA will be weighted by $\mathrm{1}/\mathrm{2}$, and those of PQRB and PQRC will be weighted by $\mathrm{1}/\mathrm{4}$, when calculating their contributions to the mean. In the UBER code, a practical choice for the number of children at each splitting is 4, and that for the upper limit of offspring generations is 3, so that a stochastic process can be split into a maximum of 4³=64 processes. Effects of the process-splitting technique are studied in the first problem in the next section.

It still remains to find a method to project the functional value of a stochastic path when it is only partially realized. For this purpose, we insert a break point at $s=s^{\prime }\in (\mathrm{0},t)$ to the integrations in Eq. (9) and see how it transforms. We simplify the situation by only considering the stochastic processes stopping in $\overline{\mathrm{\Omega }}$ for the moment and denote the following functional integrals:

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathcal{U}}_{\mathrm{0}}^t=\exp \left(\underset{\mathrm{0}}{\overset{t}{\int }}c\mathrm{d}s-\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\mathrm{d}{k}_s\right),\text{(18)}& {\displaystyle }{\mathcal{V}}_{\mathrm{0}}^t=\underset{\mathrm{0}}{\overset{t}{\int }}u\exp \left(\underset{\mathrm{0}}{\overset{s}{\int }}c\mathrm{d}r-\underset{\mathrm{0}}{\overset{s}{\int }}\mathit{\lambda }\mathrm{d}{k}_r\right)\mathrm{d}s,\end{array}$$

in which the integrand functions c, λ, and u are as those in Eq. (9). Then, the functional ℱ^t,x with the above presumptions and notations is transformed as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\mathcal{F}}^{t,\mathit{x}}& =g_{\mathrm{0}}\left({\mathit{X}}_t\right){\mathcal{U}}_{\mathrm{0}}^t+{\mathcal{V}}_{\mathrm{0}}^t\\ & =g_{\mathrm{0}}\left({\mathit{X}}_t\right){\mathcal{U}}_{\mathrm{0}}^{s^{\prime }}{\mathcal{U}}_{s^{\prime }}^t+\left({\mathcal{V}}_{\mathrm{0}}^{s^{\prime }}+{\mathcal{U}}_{\mathrm{0}}^{s^{\prime }}{\mathcal{V}}_{s^{\prime }}^t\right)\\ & =\left[g_{\mathrm{0}}\left({\mathit{X}}_t\right){\mathcal{U}}_{s^{\prime }}^t+{\mathcal{V}}_{s^{\prime }}^t\right]{\mathcal{U}}_{\mathrm{0}}^{s^{\prime }}+{\mathcal{V}}_{\mathrm{0}}^{s^{\prime }}\\ & ={\mathcal{F}}^{t-s^{\prime },{\mathit{X}}_{s^{\prime }}}{\mathcal{U}}_{\mathrm{0}}^{s^{\prime }}+{\mathcal{V}}_{\mathrm{0}}^{s^{\prime }},\end{array}\end{array}$$

where ${\mathcal{F}}^{t-s^{\prime },{\mathit{X}}_{s^{\prime }}}$ is the functional for a stochastic process that starts from the break point $(t-s^{\prime },{\mathit{X}}_{s^{\prime }})$ and continues until s = t.

Suppose that a partial path has been realized up to $s=s^{\prime }$. From it we can readily evaluate ${\mathcal{U}}_{\mathrm{0}}^{s^{\prime }}$ and ${\mathcal{V}}_{\mathrm{0}}^{s^{\prime }}$ in Eq. (19) and therefore need an estimated ${\overline{\mathcal{F}}}^{t-s^{\prime },{\mathit{X}}_{s^{\prime }}}$ to project the functional value ${\overline{\mathcal{F}}}^{t,\mathit{x}}$, where a bar is put over all unrealized entities. Specifically, we would need these three estimates: ${\overline{\mathit{X}}}_t$, ${\overline{\mathcal{U}}}_{s^{\prime }}^t$, and ${\overline{\mathcal{V}}}_{s^{\prime }}^t$. In principle, a good estimation of ${\overline{\mathit{X}}}_t$ is given by integrating along the streamline of the $\mathit{b}(t-s,\mathit{x})$ field through ${\mathit{X}}_{s^{\prime }}$ until s=t, which means projecting for ${\overline{\mathit{X}}}_t$ along the ballistic trajectory of motion while ignoring all the stochasticity since the Wiener process has zero mean. However, this integration is not much cheaper than the realization of ${\mathcal{F}}^{t-s^{\prime },{\mathit{X}}_{s^{\prime }}}$ itself and is thus unaffordable. In anticipation that the total time length t would not be too large, especially when using a solution grid, and that $\mathit{b}(t-s,\mathit{x})$ would not vary drastically in this time interval, mapping ${\overline{\mathit{X}}}_t$ along the constant vector $\mathit{b}(t-s^{\prime },{\mathit{X}}_{s^{\prime }})$ is a good enough but much cheaper approximation. If ${\overline{\mathit{X}}}_t$ is mapped out of $\overline{\mathrm{\Omega }}$ so that $g_{\mathrm{0}}\left({\overline{\mathit{X}}}_t\right)$ is unable to be evaluated, the particular stochastic process is then disabled from splitting.

The functional values ${\overline{\mathcal{U}}}_{s^{\prime }}^t$ and ${\overline{\mathcal{V}}}_{s^{\prime }}^t$ are estimated by assuming that, for all possible stochastic paths belonging to the same solution point, there are mean functions $\overline{c}$, $\overline{\mathit{\lambda }}$, and $\overline{u}$ that are independent of time and that the mean local time is proportional to the total time length of the stochastic process so that ${\overline{k}}_s=\overline{k}s$ with $\overline{k}$ the proportionality constant. Under these assumptions, ${\overline{\mathcal{U}}}_{s^{\prime }}^t$ and ${\overline{\mathcal{V}}}_{s^{\prime }}^t$ can be expressed by

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\overline{\mathcal{U}}}_{s^{\prime }}^t=& \exp \left(\overline{c}\underset{s^{\prime }}{\overset{t}{\int }}\mathrm{d}s-\overline{\mathit{\lambda }}\underset{s^{\prime }}{\overset{t}{\int }}\mathrm{d}{k}_s\right)\\ =& \exp \left[(\overline{c}t-\overline{\mathit{\lambda }}\overline{k}t)-(\overline{c}{s}^{\prime }-\overline{\mathit{\lambda }}\overline{k}{s}^{\prime })\right]\\ =& {\displaystyle \frac{{\overline{\mathcal{U}}}_{\mathrm{0}}^t}{\exp \left({\displaystyle \frac{s^{\prime }}{t}}\ln {\overline{\mathcal{U}}}_{\mathrm{0}}^t\right)}},\end{array}\end{array}$$

and

$$\begin{array}{}\text{(21)}& {\overline{\mathcal{V}}}_{s^{\prime }}^t={\displaystyle \frac{{\overline{\mathcal{V}}}_{\mathrm{0}}^t}{{\overline{\mathcal{U}}}_{\mathrm{0}}^t-\mathrm{1}}}\left[{\overline{\mathcal{U}}}_{\mathrm{0}}^t-\exp \left({\displaystyle \frac{s^{\prime }}{t}}\ln {\overline{\mathcal{U}}}_{\mathrm{0}}^t\right)\right],\end{array}$$

if ${\overline{\mathcal{U}}}_{\mathrm{0}}^t\ne \mathrm{1}$, or by

$$\begin{array}{}\text{(22)}& {\displaystyle }{\overline{\mathcal{U}}}_{s^{\prime }}^t=\mathrm{1},\text{(23)}& {\displaystyle }{\overline{\mathcal{V}}}_{s^{\prime }}^t={\overline{\mathcal{V}}}_{\mathrm{0}}^t\left(\mathrm{1}-{\displaystyle \frac{s^{\prime }}{t}}\right),\end{array}$$

if ${\overline{\mathcal{U}}}_{\mathrm{0}}^t=\mathrm{1}$. The values of ${\overline{\mathcal{U}}}_{\mathrm{0}}^t$ and ${\overline{\mathcal{V}}}_{\mathrm{0}}^t$ can be well estimated, respectively, by the medians of ${\stackrel{\mathrm{̃}}{\mathcal{U}}}_{\mathrm{0}}^t$ and ${\stackrel{\mathrm{̃}}{\mathcal{V}}}_{\mathrm{0}}^t$ that are obtained from all previously completed stochastic paths. Medians are preferred to means here because the probabilistic distributions of these functionals are usually very skewed and heavy-tailed. This projection mechanism would become statistically more accurate with more stochastic processes having been simulated.

4 Example problems

Four example problems are provided in this section. In the first problem, we solve a Fokker–Planck equation with two source terms, one proportional to the unknown function and the other independent of the unknown function, in both a spherical coordinate system and a Cartesian coordinate system. Effects of the process-splitting technique are analyzed in this example. In the second problem, an advection-dominated Fokker–Planck equation is considered. We further show that, even for a pure advection equation, the UBER code still gives the correct solutions, although it is not designed for such an equation and may not be the most efficient method. Code behavior with advection equations is further studied in the third problem in which the treatment of periodic boundary condition is also illustrated. In the last problem, we simulate the Earth's inner radiation belt by solving its Boltzmann equation involving realistic pitch-angle diffusion and CRAND source.

4.1 Problem 1: neutron generation and diffusion in nuclear material

In this problem, we consider the diffusion and generation of neutrons in a spherical nuclear material at detonation, with an initially injected Gaussian neutron distribution from a small source at the center and a neutron-reflecting coat that allows only one-half of the surface neutrons to escape. In a spherical coordinate system, the equation, initial condition, and boundary conditions are (Serber, 1992):

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial f}{\partial t}}={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial r}}\left(r^{\mathrm{2}}D{\displaystyle \frac{\partial f}{\partial r}}\right)+Sf+v\left(r\right),\text{(25)}& {\displaystyle }f(\mathrm{0},r)=\exp \left(-{\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{0.02}}}\right),\text{(26)}& {\displaystyle \frac{\partial f}{\partial r}}{|}_{r=\mathrm{0}}=\mathrm{0},\text{(27)}& {\displaystyle }\left(D{\displaystyle \frac{\partial f}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}f\right){|}_{r=\mathrm{1}}=\mathrm{0},\end{array}$$

where f is neutron density, the constant diffusion coefficient D=0.1, the constant rate of neutron generation from chain reaction S=2.5, and $v\left(r\right)={\mathrm{10}}^{-\mathrm{6}}/(\mathrm{1}+r)$ characterizes a weak source of neutrons spontaneously emitted in the material. The values and functional forms of these coefficients are solely designed for demonstration purposes and are not meant to be experimentally accurate.

Figure 2 (a) UBER and finite-difference solutions (dashed line) to the problem in Eqs. (24)–(27). The UBER 1D solutions (circles) are obtained in a one-dimensional spherical coordinate system, and the UBER 3D solutions (triangles) are obtained in three-dimensional Cartesian coordinates along a sphere radius. (b) Left y axis: the relative errors of the UBER 3D solutions at T=0.05, respectively obtained with the same total number of stochastic processes (2048 per batch) but different upper limits of offspring generations (ν) in the process-splitting technique. ν=0 means the process splitting is turned off. Right y axis: the percentage of stochastic processes that have undergone splitting for ν=4. (c) The reduction of relative errors with an increasing number of stochastic processes at the slowest-converging solution point (r=0.76) for ν=0 (dashed line and squares) and 4 (solid line and triangles). Colors denote different numerical experiments.

UBER solutions are obtained at four time stamps and are compared with those from a staggered-grid finite-difference method (e.g., Ames, 2014), as shown in Fig. 2a. A turning point is observed in the solutions at T=0.05, which marks the transition of the dominating neutron source from chain reaction at high background density to spontaneous emission at low density. As time goes by, the effect of the spontaneous emission is overwhelmed by the fast-growing chain reaction. Even though the solutions span 8 orders of magnitude, the UBER results are virtually identical to the finite-difference ones, and statistical fluctuation, which is a typical feature in Monte Carlo methods, is not observed in these solutions due to the adaptive algorithms and the variance reduction technique.

To demonstrate UBER's ability in multiple dimensions with a complicated boundary geometry, the same problem is also solved in a three-dimensional Cartesian coordinate system along a sphere radius. In this coordinate system, the diffusion coefficient becomes a 3 × 3 matrix with each diagonal component equal to D, and the boundary condition in Eq. (27) is applied to the only boundary that is a sphere with unit radius. The solutions are over-plotted in Fig. 2a. Consistency between the one-dimensional and the three-dimensional results is quite evident.

To analyze the effects of the process-splitting technique, we repeated the three-dimensional solutions at T=0.05, but with a fixed number of stochastic processes (2048 samples per batch, 200 batches) for each solution point and with various upper limits of the offspring generations ν; ν=0 indicates that the process-splitting technique is disabled. For a fixed number of samples, the relative error of a solution is proportional to the square root of the variance of sampled functional values and determines how fast the calculation of expectation converges. The relative errors as functions of r are plotted against the left y axis of Fig. 2b, and each curve is in fact formed by the medians from eight independent and identical numerical experiments to be more statistically representative. In the range $\mathrm{0.4}<r<\mathrm{0.9}$, the relative errors are consistently reduced with higher offspring generations. At the slowest-converging point r=0.76, the process-splitting technique with a maximum of four offspring generations could reduce the relative error by an order of magnitude compared to that without splitting. For this curve (ν=4), the percentages of stochastic processes that have undergone splitting are plotted as a shaded area against the right y axis. For r<0.4, the relative errors are small and computational convergence is fast enough, so process splitting is automatically suppressed by the code to achieve an optimal speed. If the relative errors are large, usually a small fraction of split stochastic processes can be rather effective.

To further reveal the behavior of the process-splitting technique, Fig. 2c plots how the relative error reduces with an increasing number of samples (N) in the Monte Carlo procedure for the solution point at r=0.76. There are eight independent and identical numerical experiments, respectively, for ν=0 and 4, and each line represents the results from one numerical experiment. The general trend is that the relative error reduces linearly in a log–log scale plot, resembling its dependence on $N^{-\mathrm{1}/\mathrm{2}}$. However, without process splitting, the relative error often jumps up sharply due to the occurrence of a very low-probability sample that made a very large contribution, which severely slows down the computational convergence. With process splitting, such jumps are largely avoided, and on average, the code uses just a little more than $\mathrm{1}/\mathrm{100}$ of the samples without process splitting to achieve the same relative error of 0.1.

Table 2 Normalized wall clock time versus maximum offspring generations (ν) for the UBER 3D solutions at T=0.05.

^* Median value from eight independent numerical tests.

In practical UBER usage, solutions are achieved with a prescribed tolerance of relative error and an adaptive number of samples. Therefore, fast convergence with the process-splitting technique could save a significant amount of computational effort even with its extra computational burden. Table 2 lists the normalized wall clock time consumed by UBER for obtaining the solution curve in three dimensions at T=0.05 with a relative error tolerance of 0.1 and a range of maximum offspring generations in process splitting. Again, each of these numbers is the median from eight independent and identical numerical experiments. With ν=3 and 4, the code is nearly an order of magnitude faster than without process splitting. The same wall clock time in these two cases indicates that the faster convergence with more offspring generations starts to be traded off by the computational overhead associated with more complicated splitting, and therefore further increasing ν would not be optimal.

Figure 3 (a) UBER (circles) and staggered-grid finite-difference (dashed line) solutions to the problem in Eqs. (28)–(31). (b) UBER (circles) and Lax–Wendroff (dashed line) solutions to the same problem but with zero diffusion.

4.2 Problem 2: magnetized plasma evolution under instability

In the second problem, we consider a Fokker–Planck equation for the pitch-angle distribution of a magnetized plasma (e.g., Dendy, 1990). Suppose that the electrons are initially in a sin ²(x) background pitch-angle distribution with x as the pitch angle. An electron beam is injected into the system centered at pitch angle x=0.4. In addition to pitch-angle diffusion, the injected beam excites some kind of plasma instability that kinetically transports the distribution toward the $\mathit{\pi }/\mathrm{2}$ pitch angle. The equation, initial condition, and boundary conditions are written as

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\partial f}{\partial t}}={\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial x}}\left[GD(t,x){\displaystyle \frac{\partial f}{\partial x}}\right]-{\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial x}}\left[Gh\left(x\right)f\right],\text{(29)}& {\displaystyle }f(\mathrm{0},x)={\sin }^{\mathrm{2}}\left(x\right)+\exp \left[-{\displaystyle \frac{(x-\mathrm{0.4}{)}^{\mathrm{2}}}{\mathrm{0.02}}}\right],\text{(30)}& {\displaystyle }f{|}_{x=\mathrm{0.05}}=\mathrm{0},\text{(31)}& {\displaystyle \frac{\partial f}{\partial x}}{|}_{x=\mathit{\pi }/\mathrm{2}}=\mathrm{0},\end{array}$$

where f is the electron distribution function, the Jacobian determinant G=sin (x), the diffusion coefficient $D(t,x)=(\mathrm{1}/\mathrm{2}{\mathit{\pi }}^{\mathrm{2}})\text{erf}(t/\mathrm{2})[\mathrm{1}+{\sin }^{\mathrm{2}}(\mathrm{2}x\left)\right]$, and the advection coefficient h(x)=cos (x). Note that, in most of the x range, the advection coefficient is about an order of magnitude larger in value than the diffusion coefficient. Equation (30) indicates a loss cone at pitch angle x=0.05. UBER solutions for this problem are plotted in Fig. 3a as circles and are in excellent agreement with those from the staggered-grid finite-difference method. In these solutions, the beam evolves toward $x=\mathit{\pi }/\mathrm{2}$ because of the kinetic advection. As the system relaxes, the beam eventually merges into the background, and a final stable distribution is then approached.

Equation (28) degenerates to a continuity equation if pitch-angle diffusion is turned off by setting D(t,x) to zero. Even for such a pure advection problem, UBER can still obtain accurate and robust solutions compared to the widely used Lax–Wendroff method (e.g., Ames, 2014), as shown in Fig. 3b. Before T=2, an advection of the beam toward $x=\mathit{\pi }/\mathrm{2}$ is seen in the solutions without dispersion, and UBER results are almost identical to the Lax–Wendroff ones. The system, however, is unstable due to the positive advection velocity at $x<\mathit{\pi }/\mathrm{2}$ and the zero advection velocity at $x=\mathit{\pi }/\mathrm{2}$ so that the electron distribution will be piled up near $x=\mathit{\pi }/\mathrm{2}$ and ultimately evolve into a singularity. For this reason, the Lax–Wendroff method begins to fail at T=2 by generating unphysical negative solutions near x=1.3 and will be divergent henceforth; this problem would require careful choices of additional flux limiters to overcome. UBER nonetheless gives the correct results that still resolve the peak height and position of the beam.

Figure 4 (a) Initial condition (solid black line) and UBER solutions of Eqs. (32) and (33) for a constant drift frequency after 1 round (T=1, blue dashed line) and 20 rounds (T=20, red dotted line) of drift around Earth. (b) Initial condition and UBER solutions of the same problem but with an accelerating drift frequency after 1 round ($T=\sqrt{\mathrm{2}}$) and 20 rounds ($T=\mathrm{2}\sqrt{\mathrm{10}}$) of drift around Earth.

4.3 Problem 3: particle adiabatic drift around Earth

In this problem, we study the UBER code behavior when solving the advection equation resulting from Eq. (2) but without any non-adiabatic process or source and loss, which describes the adiabatic drift of radiation belt particles around Earth. To be specific, the equation and boundary condition are

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+\dot{\mathit{\phi }}{\displaystyle \frac{\partial \overline{f}}{\partial \mathit{\phi }}}=\mathrm{0},\text{(33)}& {\displaystyle }\overline{f}{|}_{\mathit{\phi }=\mathrm{0}}=\overline{f}{|}_{\mathit{\phi }=\mathrm{2}\mathit{\pi }},\end{array}$$

and an idealized initial condition g₀(φ) consisting of a step function, a triangle, and a semicircle is given as that in Fig. 4a, which can induce unphysical negative solutions to finite-difference methods due to its discontinuities and infinite gradients (e.g., Aseev et al., 2016). Equation (33) specifies the periodic boundary condition for the drift phase φ. In the UBER code, the periodic boundary condition is not really considered a boundary condition; rather, it is dealt with by extending the computational domain to include multiple periods so that the Itô stochastic processes would not move out of the domain within the given time duration, except for stopping on other first type boundaries. In this specific problem, the time stamp for obtaining solutions is every one-half drift period; therefore, the computational domain is extended for one extra period of φ from 0 to −2π since the stochastic processes retrograde in time. However, solutions are only sought in the right half of the domain for φ between 0 and 2π at each time stamp, and after that, they are copied to the left half to form the entire initial condition for the next time stamp.

By the method of characteristics, Eq. (32) is trivially solved for a constant drift frequency $\dot{\mathit{\phi }}=\mathrm{2}\mathit{\pi }$, and the solution $\overline{f}(t,\mathit{\phi })=g_{\mathrm{0}}(\mathit{\phi }-\mathrm{2}\mathit{\pi }t)$ preserves its shape and returns to its initial position after each period of t=1. As shown in Fig. 4a, this is exactly what the UBER solutions do. Unlike finite-difference methods which often introduce numerical dispersion to advection equations, the UBER solutions can be exact because, without the stochastic terms, the solver of the SDE (Eq. 7) actually integrates along the characteristic curves of the advection equation, and for constant advection velocity, these characteristic curves are straight lines for which the Euler scheme is exact.

The Euler scheme is well known to overshoot when integrating along a curve that consistently curls to one side, and this would cause error in UBER solutions. To show this, we replace the drift frequency by $\dot{\mathit{\phi }}=\mathrm{2}\mathit{\pi }t$ in Eq. (32), for which the characteristic curves are now a family of parabolas whose equations are $\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}=\mathit{\pi }{t}^{\mathrm{2}}$, with φ₀ the parameter. Note that the drift is no longer periodic but forever accelerating. UBER solutions after the first round ($T=\sqrt{\mathrm{2}}$) and the 20th round ($T=\mathrm{2}\sqrt{\mathrm{10}}$) of drift are plotted in Fig. 4b. Compared with the initial condition, the UBER solutions are slightly displaced due to the overshoot, and the solution shapes are very slightly distorted, too. For radiation belt applications, though, any frequency ${\dot{\mathit{\phi }}}^{\mathit{\iota }}$ must meanwhile be an approximate periodic function of φ^ι, and hence the corresponding characteristic curves would at most wobble rather than curling, and the Euler scheme overshoot largely cancels rather than accumulating. Therefore, we would not expect the UBER solution errors in those applications to be worse than that illustrated in this case.

4.4 Problem 4: Earth's inner radiation belt simulation

In the last problem, we demonstrate UBER's ability to solve a radiation belt Boltzmann equation by performing an inner radiation belt simulation involving both the stably trapped (out of the drift loss cone) and the quasi-trapped (in the drift loss cone) electron populations. Inspired by Xiang et al. (2020), we consider the 304 keV electrons at McIlwain's L_M=1.25, which are subject to pitch-angle scattering caused by Coulomb collisions with upper-atmospheric neutrals and ionospheric ions and electrons. The equation, initial condition, and boundary conditions are

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+\dot{\mathit{\phi }}{\displaystyle \frac{\partial \overline{f}}{\partial \mathit{\phi }}}={\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial {\mathit{\alpha }}_{\mathrm{0}}}}\left(G{D}_{\mathit{\alpha }\mathit{\alpha }}{\displaystyle \frac{\partial \overline{f}}{\partial {\mathit{\alpha }}_{\mathrm{0}}}}\right)+{\displaystyle \frac{S_{\mathrm{e}}}{p^{\mathrm{2}}}},\text{(35)}& {\displaystyle }\overline{f}(\mathrm{0},\mathit{\phi },{\mathit{\alpha }}_{\mathrm{0}})=\mathrm{0},\text{(36)}& {\displaystyle }\overline{f}{|}_{\mathit{\phi }=\mathrm{0}}=\overline{f}{|}_{\mathit{\phi }=\mathrm{2}\mathit{\pi }},\text{(37)}& {\displaystyle }\overline{f}{|}_{{\mathit{\alpha }}_{\mathrm{0}}={\mathit{\alpha }}_{\mathrm{L}}}=\mathrm{0},\text{(38)}& {\displaystyle \frac{\partial \overline{f}}{\partial {\mathit{\alpha }}_{\mathrm{0}}}}{|}_{{\mathit{\alpha }}_{\mathrm{0}}=\mathit{\pi }/\mathrm{2}}=\mathrm{0}.\end{array}$$

Figure 5 (a) Bounce-averaged pitch-angle diffusion coefficient D_αα (s⁻¹) for 304 keV electrons. The blank area is in the bounce loss cone. (b) CRAND electron source rate $S_{\mathrm{e}}/p^{\mathrm{2}}$ (c³ ${\mathrm{cm}}^{-\mathrm{3}}{\mathrm{MeV}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) for 304 keV electrons. The black line plots the variation of dipole L-shell versus geomagnetic longitude against the right y axis, corresponding to the McIlwain's L_M=1.25. (c) Calculated electron fluxes (${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}{\mathrm{MeV}}^{-\mathrm{1}}$) at T = 2 h. (d) Calculated electron fluxes (${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}{\mathrm{MeV}}^{-\mathrm{1}}$) at T = 10 h.

In Eq. (34), drift frequency,

$$\begin{array}{}\text{(39)}& \dot{\mathit{\phi }}={\displaystyle \frac{\mathrm{3}cL{R}_{\mathrm{E}}}{e{\mathit{\mu }}_{\mathrm{E}}}}{\displaystyle \frac{p^{\mathrm{2}}}{m_{\mathrm{e}}}}{\displaystyle \frac{D\left(\sin {\mathit{\alpha }}_{\mathrm{0}}\right)}{T\left(\sin {\mathit{\alpha }}_{\mathrm{0}}\right)}},\end{array}$$

is evaluated using dipole field approximation (Schulz, 1991), in which c is the speed of light in a vacuum, L is dipole L-shell, R_E is the radius of Earth, e is the elementary charge, μ_E is the magnetic moment of Earth's intrinsic dipole field, m_e is electron mass, p is electron momentum, α₀ is electron equatorial pitch angle, and the functions D(sin α₀) and T(sin α₀) are bounce motion integrals in the dipole field that are given in Schulz (1991, pp. 205–210). For simplicity, we ignore the dependence of drift frequency on φ so that the drift phase becomes equivalent to geomagnetic longitude. The Jacobian determinant G=T(sin α₀)sin (2α₀). The bounce-averaged pitch-angle diffusion rate is empirically given by

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{D}_{\mathit{\alpha }\mathit{\alpha }}=& {\mathrm{10}}^{-\mathrm{5}}\exp \left\{\mathrm{92.55}\left[{\cos }^{\mathrm{4}}{\mathit{\alpha }}_{\mathrm{0}}-{\cos }^{\mathrm{4}}{\mathit{\alpha }}_{\mathrm{L}}\left(\mathit{\phi }\right)\right]\right\}\\ & +{\mathrm{10}}^{-\mathrm{9}}\left({\mathrm{s}}^{-\mathrm{1}}\right),\end{array}\end{array}$$

which features quantitative resemblance to that calculated by realistic atmosphere and ionosphere models in Xiang et al. (2020). In this expression, α_L(φ) is the bounce loss cone angle dependent on geomagnetic longitude that is determined by drift shell tracing in the International Geomagnetic Reference Field (IGRF; Finlay et al., 2010). D_αα as a function of φ and α₀ is plotted in Fig. 5a: it is only significant near the bounce loss cone and in the South Atlantic Anomaly (SAA) centered at about 20^∘ of geomagnetic longitude due to the closer proximity of the drift shell to the upper atmosphere in these regions. The CRAND source rate $S_{\mathrm{e}}/p^{\mathrm{2}}$ is approximated by (Lenchek et al., 1961; Selesnick, 2015)

$$\begin{array}{}\text{(41)}& \begin{array}{rl}& {\displaystyle \frac{S_{\mathrm{e}}}{p^{\mathrm{2}}}}\approx \mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{12}}{\displaystyle \frac{(E_{max}-E{)}^{\mathrm{2}}}{L^{\mathrm{2.7}}\sin {\mathit{\alpha }}_{\mathrm{0}}}}\\ & \left(c^{\mathrm{3}}{\mathrm{cm}}^{-\mathrm{3}}{\mathrm{MeV}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right),\end{array}\end{array}$$

where E_max is the maximum kinetic energy (782 keV) available to electrons from neutron β decay and E is the electron kinetic energy in question, with both being measured in units of the electron rest energy (511 keV), and L is the dipole L-shell, which is a variable dependent on geomagnetic longitude due to multipoles of the Earth's magnetic field. Figure 5b plots the CRAND source rate and the dipole L-shell values corresponding to McIlwain's L_M=1.25 obtained from drift shell tracing in IGRF, which vary from less than 1.2 in the SAA to above 1.3 near 180^∘ of geomagnetic longitude.

The simulation is performed with an initially empty radiation belt as indicated by Eq. (35), and electrons are gradually generated by the CRAND source and meanwhile lost to the bounce loss cone. Figure 5c and d show the solution electron fluxes calculated from $\overline{j}=\overline{f}{p}^{\mathrm{2}}$ after 2 h (about one drift period) and 10 h, respectively. The characteristic west–east electron flux gradient is formed for the quasi-trapped population (${\mathit{\alpha }}_{\mathrm{0}}<\mathrm{60}{}^{\circ }$) within the first 2 h and changes very little over time because the SAA sweeps these electrons out every drift period. Weak pitch-angle diffusion of electron fluxes from the quasi-trapped population toward the stably trapped population can be observed at T = 2 h when the stably trapped fluxes are still low due to the stronger source rate in the quasi-trapped region. At T = 10 h, the direction of the pitch-angle diffusion is reversed. Even with atmospheric loss, the CRAND source is strong enough to continuously contribute to the trapped electron fluxes, which are increased by 1 order of magnitude in 8 h.

5 Conclusions and discussion

In conclusion, we have built a numerical solver for the general form of kinetic equations that appear in radiation belt studies. Based on the SDE method, the solver is coded to work in arbitrarily provided coordinate systems of up to three dimensions with user-specified boundary geometry, boundary conditions, and equation terms. We have also designed adaptive algorithms and a variance reduction technique for the SDE method, which enhanced its computational speed by 1 order of magnitude in our test. The example problems in this article demonstrated the solver's versatility and robustness in dealing with a range of problems that might each require a different solver in other methods. The solver, named UBER, has been programmed into a Fortran library that can be easily incorporated with other more complicated space physics models.

Figure 6 Left y axis: UBER code speedup, defined as the ratio between wall clock time of serial and parallel executions, as a function of parallel threads for solving Problem 1 in spherical coordinates. Each data point (cross) is the median from eight independent and identical numerical experiments, with the associated error bar covering the range of those eight experiments. The dashed line gives a fit to the data points using Amdahl's law with parameter p=0.85 (see text for its meaning). The right y axis marks the actual wall clock time in seconds for the numerical experiments on a Linux server with a maximum of 86 CPUs (the rightmost data point) at 2.1 GHz of floating-point operation frequency per processor.

The strengths of the SDE method lie in its abilities to solve problems that are difficult for other methods, but not in its speed. Even with the presented improvements, the SDE method is in general nowhere near the finite-difference counterparts in terms of speed due to its Monte Carlo nature. However, Monte Carlo methods are perfectly parallelized for no communication among parallel tasks; therefore, the UBER code gains handsome speedup from parallel computation. Figure 6 displays the UBER code speedup, measured as the ratio between wall clock time of serial and parallel executions, against different numbers of parallel threads employed in solving Problem 1 in spherical coordinates. For a fixed amount of computation, speedup is limited by the portion of work that cannot be parallelized; this is described by Amdahl's law (inset formula in Fig. 6; Amdahl, 1967), which fits the numerical experiment data with parameter p=0.85, meaning that effectively 85 % of the total work has been parallelized. With 86 threads, the code approaches its theoretical maximum speedup of 6.7 for this problem, which uses about 620 s of wall clock time on a Linux server. The referenced finite-difference code uses 65 s in serial execution to solve the problem on the same server. In other words, the UBER code, in its ideal parallelization, is an order of magnitude slower than serial finite-difference methods, and it is also roughly the case for Problem 4 (Zheng Xiang, personal communication, 2021) for which the UBER code spends 22 min to simulate one drift period. This conclusion, however, is rather qualitative because in these comparisons different methods have used quite different grids. Finite-difference methods generally require a much finer grid than the SDE method in order to deal with short-wavelength components of the solution, and time spacing is restricted accordingly by the Courant–Friedrichs–Lewy condition (e.g., Press et al., 1992). On the other hand, the SDE method affords a grid that is only concentrated in interesting areas as long as global interpolation can be warranted, and the time stamps for dumping solutions are up to user choice. Moreover, if only local solutions are sought over a short time period, the SDE method could totally disregard the uninteresting solutions and might be more efficient. For radiation belt simulations in particular, finite-difference models usually rely on the operator splitting technique (e.g., Press et al., 1992) to solve tensorial diffusion in three dimensions, which involves frequent interpolation between two sets of grids and hence greatly hinders speed, whereas problem dimensionality has little effect on the SDE method.

Several other forms of radiation belt kinetic equations should also be solvable by the method presented in this article. In formulating the Boltzmann equation (Eq. 2), we have assumed that the unperturbed particle Hamiltonian H₀ is independent of phases of particle motion. For lower-energy ring current particles, the convective electric field potential energy is not negligible in their Hamiltonian, and therefore H₀ would be dependent on the drift phase. As such, expanding the Poisson bracket $[\overline{f},H_{\mathrm{0}}]$ on the left-hand side of the Boltzmann equation will result in additional terms involving partial differentials with respect to the generalized momenta Q^ι (see Appendix A). For a radiation belt model including ring current particles, the general form of the Boltzmann equation will be

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+{\dot{\mathit{\phi }}}^{\mathit{\iota }}{\displaystyle \frac{\partial \overline{f}}{\partial {\mathit{\phi }}^{\mathit{\iota }}}}+{\dot{Q}}^{\mathit{\iota }}{\displaystyle \frac{\partial \overline{f}}{\partial Q^{\mathit{\iota }}}}=\mathrm{\dots },\end{array}$$

in which the omitted right-hand side is exactly the same as that of Eq. (2). The Boltzmann equations of the so-called four-dimensional radiation belt models, such as the CIMI model (Fok et al., 2014), the VERB-4D model (Aseev et al., 2016), and the K2 MHD-particle model (Elkington et al., 2019), are of this type. Similar to the treatment of Eq. (2), Eq. (42) can be obtained from Eq. (1) by introducing the new coordinates $x^{\mathit{\xi }}=\mathit{\{}{Q}^{\mathit{\mu }},Q^{\mathit{\iota }},{\mathit{\phi }}^{\mathit{\iota }}\mathit{\}}\equiv \mathit{\{}{x}^{\mathit{\mu }},x^{\mathit{\kappa }}\mathit{\}}$, which enlarges the index set from α to ξ, and performing the following transformation of equation coefficients:

$$\begin{array}{}\text{(43)}& \left\{\begin{array}{l}{\stackrel{\mathrm{̃}}{D}}^{\mathit{\mu }\mathit{\nu }}=D^{\mathit{\mu }\mathit{\nu }},{\stackrel{\mathrm{̃}}{D}}^{\mathit{\xi }\mathit{\kappa }}=\mathrm{0},\\ {\stackrel{\mathrm{̃}}{h}}^{\mathit{\xi }}=\mathit{\{}{h}^{\mathit{\mu }},{\dot{x}}^{\mathit{\kappa }}\mathit{\}},\\ \stackrel{\mathrm{̃}}{S}=S+{\partial }_{\mathit{\kappa }}{\dot{x}}^{\mathit{\kappa }}+{\dot{Q}}^{\mathit{\iota }}{\partial }_{\mathit{\iota }}\ln G,\\ \stackrel{\mathrm{̃}}{v}=v.\end{array}\right.\end{array}$$

Therefore, the Boltzmann equation (Eq. 42) can also be solved by the method presented in this article in principle. However, such four-dimensional simulations are beyond the current scope of the UBER code since it is only coded for up to three dimensions in space.

Nonlinear evolution of phase space density occurs when the particle scatterings are not only small-scale but also large-scale, usually as a result of trapping by intense plasma waves (e.g., Bortnik et al., 2008; Albert et al., 2013). In this case, the right-hand side of the kinetic equation must include terms of nonlocal transport of phase space density by these large-scale scatterings, and the equation is formulated as (Artemyev et al., 2016; Zheng et al., 2019)

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{\displaystyle \frac{\partial \overline{f}}{\partial t}}=& {\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial Q^{\mathit{\alpha }}}}\left(G{D}^{\mathit{\alpha }\mathit{\beta }}{\displaystyle \frac{\partial \overline{f}}{\partial Q^{\mathit{\beta }}}}\right)-{\displaystyle \frac{\mathrm{1}}{G}}{\displaystyle \frac{\partial }{\partial Q^{\mathit{\alpha }}}}\left(G{h}^{\mathit{\alpha }}\overline{f}\right)\\ & -\left(\int P_{Q\to \stackrel{\mathrm{̃}}{Q}}\stackrel{\mathrm{̃}}{G}\mathrm{d}{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }}\right)\overline{f}+\int P_{\stackrel{\mathrm{̃}}{Q}\to Q}\stackrel{\mathrm{̃}}{f}\stackrel{\mathrm{̃}}{G}\mathrm{d}{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }},\end{array}\end{array}$$

in which $\stackrel{\mathrm{̃}}{f}$ is a shorthand for the function $\overline{f}(t,{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }})$, and $\stackrel{\mathrm{̃}}{G}$ is the Jacobian determinant evaluated at ${\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }}$. With nonlinear wave–particle interactions, the phase bunching effect gives rise to advection characterized by the coefficients h^α. The function $P_{Q\to \stackrel{\mathrm{̃}}{Q}}$ is the trapping probability density per unit time from Q^α to ${\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }}$; that is, particles are trapped by the wave field at Q^α and subsequently escape from trapping at ${\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }}$. This is considered a known function which can be evaluated from single particle behaviors by either perturbation theory of Hamiltonian mechanics (e.g., Artemyev et al., 2016) or test-particle simulations (e.g., Vainchtein et al., 2018). Note that, since the unknown function is contained in the last integral term, Eq. (44) is an integro-differential equation. However, formal similarity between Eq. (44) and the Fokker–Planck equation (Eq. 1) suggests that an asymptotic solution of Eq. (44) may be achieved by Taylor expanding $\stackrel{\mathrm{̃}}{f}$ as

$$\begin{array}{}\text{(45)}& \stackrel{\mathrm{̃}}{f}={\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}+{\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}^{\prime }t+\mathrm{\dots },\end{array}$$

where ${\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}=\overline{f}(\mathrm{0},{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }})$, and ${\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}^{\prime }$ indicates its time-derivative function evaluated at t=0. When applying the SDE method with a solution grid, the functions ${\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}$ and ${\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}^{\prime }$ can be obtained from solutions of previous time stamps. Then, by defining the following coefficients

$$\begin{array}{}\text{(46)}& {\displaystyle }S\left(Q^{\mathit{\alpha }}\right)=-\int P_{Q\to \stackrel{\mathrm{̃}}{Q}}\stackrel{\mathrm{̃}}{G}\mathrm{d}{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }},\text{(47)}& {\displaystyle }\begin{array}{rl}& v(t,Q^{\mathit{\alpha }})=\int P_{\stackrel{\mathrm{̃}}{Q}\to Q}{\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}\stackrel{\mathrm{̃}}{G}\mathrm{d}{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }}\\ & +t\int P_{\stackrel{\mathrm{̃}}{Q}\to Q}{\stackrel{\mathrm{̃}}{f}}_{\mathrm{0}}^{\prime }\stackrel{\mathrm{̃}}{G}\mathrm{d}{\stackrel{\mathrm{̃}}{Q}}^{\mathit{\alpha }}+\mathrm{\dots },\end{array}\end{array}$$

which are now known functions, Eq. (44) is transformed into the form of Eq. (1) and is readily solvable by the UBER code. In this way, simulations of nonlinear wave–particle interactions in radiation belts could hence be unified with well-developed simulations in quasi-linear theory.

Appendix A: A formal derivation of radiation belt kinetic equations

We consider, for simplicity, a hypothetical radiation belt whose particle motion has two well-separated periods which define two pairs of action-angle variables $\mathit{\{}{\mathit{\phi }}^{\mathrm{1}},{\mathit{\phi }}^{\mathrm{2}},J^{\mathrm{1}},J^{\mathrm{2}}\mathit{\}}$. We assume that the phase angle φ¹ changes much faster than φ² and hence call {φ¹,J¹} the fast variables and {φ²,J²} the slow variables. We further assume, for a moment, that the Hamiltonian of particle motion,

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& H({\mathit{\phi }}^{\mathrm{1}},{\mathit{\phi }}^{\mathrm{2}},J^{\mathrm{1}},J^{\mathrm{2}},t)=H_{\mathrm{0}}({\mathit{\phi }}^{\mathrm{2}},J^{\mathrm{1}},J^{\mathrm{2}},t)\\ & +\mathit{\delta }H({\mathit{\phi }}^{\mathrm{1}},J^{\mathrm{1}},J^{\mathrm{2}},t),\end{array}\end{array}$$

is constituted of an unperturbed part H₀ that depends on the slow phase and a perturbation δH that is caused by electromagnetic forces whose variation timescale is shorter than the periodicity $\mathrm{2}\mathit{\pi }/{\dot{\mathit{\phi }}}^{\mathrm{1}}$. Apparently, δH is a periodic function of φ¹. Upon averaging over φ¹, the perturbation cancels out so that

$$\begin{array}{}\text{(A2)}& \langle \mathit{\delta }H\rangle \equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathit{\delta }H\mathrm{d}{\mathit{\phi }}^{\mathrm{1}}=\mathrm{0}.\end{array}$$

These presumptions allow the particle phase space density,

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& f({\mathit{\phi }}^{\mathrm{1}},{\mathit{\phi }}^{\mathrm{2}},J^{\mathrm{1}},J^{\mathrm{2}},t)=\overline{f}({\mathit{\phi }}^{\mathrm{2}},J^{\mathrm{1}},J^{\mathrm{2}},t)\\ & +\mathit{\delta }f({\mathit{\phi }}^{\mathrm{1}},J^{\mathrm{1}},J^{\mathrm{2}},t)\end{array}\end{array}$$

to be so decomposed into a fast-phase-averaged part $\overline{f}\equiv \langle f\rangle$ and a perturbation δf, also periodic in φ¹, which has 〈δf〉=0 by definition.

For a collisionless plasma, evolution of f is governed by the Vlasov equation:

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\partial f}{\partial t}}+{\dot{\mathit{\phi }}}^{\mathit{\alpha }}{\displaystyle \frac{\partial f}{\partial {\mathit{\phi }}^{\mathit{\alpha }}}}+{\dot{J}}^{\mathit{\alpha }}{\displaystyle \frac{\partial f}{\partial J^{\mathit{\alpha }}}}=\mathrm{0},\end{array}$$

in which the index $\mathit{\alpha }=\mathrm{1},\mathrm{2}$. Expressing ${\dot{\mathit{\phi }}}^{\mathit{\alpha }}$ and ${\dot{J}}^{\mathit{\alpha }}$ by Hamilton's canonical equations, the Vlasov equation can be expanded in light of Eqs. (A1) and (A3). When averaging the expanded equation over φ¹, all terms to the first order in perturbation vanish due to either their null phase average or periodicity in φ¹, and the remaining terms form the equation

$$\begin{array}{}\text{(A5')}& \begin{array}{rl}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+{\displaystyle \frac{\partial H_{\mathrm{0}}}{\partial J^{\mathrm{2}}}}{\displaystyle \frac{\partial \overline{f}}{\partial {\mathit{\phi }}^{\mathrm{2}}}}-{\displaystyle \frac{\partial H_{\mathrm{0}}}{\partial {\mathit{\phi }}^{\mathrm{2}}}}{\displaystyle \frac{\partial \overline{f}}{\partial J^{\mathrm{2}}}}=\\ & -〈{\displaystyle \frac{\partial \mathit{\delta }H}{\partial J^{\mathrm{1}}}}{\displaystyle \frac{\partial \mathit{\delta }f}{\partial {\mathit{\phi }}^{\mathrm{1}}}}-{\displaystyle \frac{\partial \mathit{\delta }H}{\partial {\mathit{\phi }}^{\mathrm{1}}}}{\displaystyle \frac{\partial \mathit{\delta }f}{\partial J^{\mathrm{1}}}}〉,\end{array}\end{array}$$

organized into Poisson brackets with respect to the canonical coordinates {φ^α,J^α},

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+[\overline{f},H_{\mathrm{0}}]=-〈[\mathit{\delta }f,\mathit{\delta }H]〉.\end{array}$$

In fact, this equation form is more neatly derived from Liouville's theorem (e.g., Goldstein, 1980, chap. 9), which says

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\mathrm{d}f}{\mathrm{d}t}}={\displaystyle \frac{\partial f}{\partial t}}+[f,H]=\mathrm{0}.\end{array}$$

Phase averaging Eq. (A6) over φ¹ and noting that $\langle [\mathit{\delta }f,H_{\mathrm{0}}]\rangle =\langle [\overline{f},\mathit{\delta }H]\rangle =\mathrm{0}$ directly gives Eq. (A5).

Equation (A5) appears in the form of a Boltzmann equation for the phase-averaged phase space density: the left-hand side describes evolution of the unperturbed system in the slow variables, whereas the right-hand side, involving only the perturbed quantities and the fast variables, serves the role of a collision integral. Indeed, it can be viewed as “collisions” between particles and the perturbing electromagnetic field. In this regard, we symbolically denote the right-hand side ${\displaystyle }\langle ({\displaystyle \frac{\partial f}{\partial t}}{)}_{\mathrm{w}}\rangle$ in analogy to that caused by real collisions, with the subscript designating wave–particle interaction.

Equation (A5) is closed when its right-hand side can be expressed in terms of $\overline{f}$ under certain approximations. If both δH and δf are small compared to their unperturbed counterparts, δf can be directly solved from the linearized Vlasov equation retaining only the first-order terms in expansion and, after mathematical transformations, gives the expression (e.g., Lerche, 1968; Kaufman, 1972)

$$\begin{array}{}\text{(A7)}& 〈{\left({\displaystyle \frac{\partial f}{\partial t}}\right)}_{\mathrm{w}}〉={\displaystyle \frac{\partial }{\partial J^{\mathrm{1}}}}\left(D_{\mathrm{w}}{\displaystyle \frac{\partial \overline{f}}{\partial J^{\mathrm{1}}}}\right),\end{array}$$

where the coefficient D_w is a functional of δH. The corresponding theory is called the quasi-linear theory. However, when the perturbing electromagnetic wave is sufficiently coherent, δf may become large even if δH remains small. In this situation, ${\displaystyle }\langle ({\displaystyle \frac{\partial f}{\partial t}}{)}_{\mathrm{w}}\rangle$ is estimated by considering particle phase trajectories near the resonance point (e.g., Artemyev et al., 2016). The result would then contain corrections to Eq. (A7), which are due to particles trapped in phase with the wave, whose formulation in the current setup could be inferred from that of Eq. (44) in the body; we hereby do not elaborate.

Taking account of collisions in the plasma would introduce to Eq. (A4) a collision term not describable by the single-particle Hamiltonian so that the transport equation becomes the Boltzmann equation:

$$\begin{array}{}\text{(A8)}& {\displaystyle \frac{\partial f}{\partial t}}+{\dot{\mathit{\phi }}}^{\mathit{\alpha }}{\displaystyle \frac{\partial f}{\partial {\mathit{\phi }}^{\mathit{\alpha }}}}+{\dot{J}}^{\mathit{\alpha }}{\displaystyle \frac{\partial f}{\partial J^{\mathit{\alpha }}}}={\left({\displaystyle \frac{\partial f}{\partial t}}\right)}_{\mathrm{c}}.\end{array}$$

Following the same treatment from Eqs. (A4) to (A5), Eq. (A8) leads to the fast-phase-averaged equation

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\partial \overline{f}}{\partial t}}+[\overline{f},H_{\mathrm{0}}]=〈{\left({\displaystyle \frac{\partial f}{\partial t}}\right)}_{\mathrm{w}}〉+〈{\left({\displaystyle \frac{\partial f}{\partial t}}\right)}_{\mathrm{c}}〉.\end{array}$$

For Coulomb collisions, small-angle scatterings at large impact parameters dominate due to the long range of Coulomb force, and consequently the phase-averaged collision integral can be expanded into a Fokker–Planck form in the generalized momenta that are changed by the collisions (Lifshitz and Pitaevskii, 1981, chaps. 2 and 4), i.e.,

$$\begin{array}{}\text{(A10)}& 〈{\left({\displaystyle \frac{\partial f}{\partial t}}\right)}_{\mathrm{c}}〉={\displaystyle \frac{\partial }{\partial J^{\mathrm{1}}}}\left(D_{\mathrm{c}}{\displaystyle \frac{\partial \overline{f}}{\partial J^{\mathrm{1}}}}\right)-{\displaystyle \frac{\partial }{\partial J^{\mathrm{1}}}}\left(h_{\mathrm{c}}\overline{f}\right),\end{array}$$

which usually only involves the fast momentum J¹ on timescales much shorter than $\mathrm{2}\mathit{\pi }/{\dot{\mathit{\phi }}}^{\mathrm{2}}$. The transport coefficients D_c and h_c are determined from the particle species and their collision cross sections. We note again that, in the phase-averaged kinetic equation (Eq. A9), slow and fast variables are separated onto each side of the equation.

Neglecting the source and loss terms, the quasi-linear kinetic equations in the body of this paper could all be recovered from Eqs. (A7), (A9), and (A10), which are already in the same form as Eq. (42). If there are no slow variables, the Poisson bracket on the left-hand side of Eq. (A9) vanishes, and the equation reduces to the Fokker–Planck equation (Eq. 1). If there are slow variables but φ² is cyclic to the unperturbed Hamiltonian, the left-hand side of Eq. (A9) would then contain the first two terms shown in Eq. (A5'), which is in the form of Eq. (2). In this case, the dependence of $\overline{f}$ on φ² is introduced by means other than the Hamiltonian, such as the φ²-dependent boundary geometry, boundary conditions, or collision terms.