Randomized quasi-Monte Carlo simulation of fast-ion thermalization

The quasi-Monte Carlo method is a well-established method for improving the statistical convergence for problems solved with the standard Monte Carlo method. It is commonly used for improving high-dimensional integration, but there is growing research on the applicability of the method for simulation of diffusion processes. Simulation of diffusion is achieved by moving particles in the Lagrangian framework according to the associated stochastic differential equation (SDE). Examples in fusion plasma physics are the simulation of neutral beam injection [1, 2] and radio frequency heating [3]. In the simulation, two sources of error appear; the time discretization error from the numerical scheme of the SDE and the statistical error from the finite number of particles describing the density. The standard Monte Carlo method gives a statistical convergence of αN^β with β = −1/2. To improve this value one often tries to reduce the value of α with some sort of variance reduction method. Common methods are the importance sampling and control-variate methods. In the plasma community a variance reduction concept known as the δf-method is commonly used, which is a family of different methods based on combinations of the importance sampling and control-variate methods [4–9]. Instead of reducing α, the noise can be reduced by improving the order of convergence, β. This can be achieved by replacing the pseudo-random numbers in the Monte Carlo method with deterministic numbers with the low-discrepancy property. An early application of low-discrepancy sequences is the quiet start method in particle-in-cell simulations [10, 11] where quasi-random points are used for sampling the initial particle distribution.

The asymptotic convergence of the quasi-Monte Carlo method is,

$$\begin{equation} \mathcal{O}(N^{-1}\log(N)^s), \end{equation} \tag{ 1 }$$

where s is the number of dimensions. For modest dimensions, the log(N)^s term can be ignored. The convergence is dependent on the number of dimensions, forcing a greater number of particles to be used as the number of dimensions increase. For integration problems over the unit hypercube we require at least a $N \geqslant \exp (s)$ number of particles, (from ∂_NN⁻¹log(N)^s = 0), for (1) to converge. The number of particles required explodes as the number of dimensions increase. This suggests that the quasi-Monte Carlo method is only efficient for small dimensions and suffers from the 'curse of dimensionality'. Fortunately it has been shown in [12] and references therein, that the quasi-Monte Carlo method has better convergence in practice than theory predicts. This can partly be explained by the fact that a large class of integrands have a so called low effective dimension; the function to be integrated has most of its variation located in few dimensions.

The purpose of this paper is to test the applicability of the quasi-Monte Carlo method for simulation of the fast-ion thermalization process of a neutral beam injection scenario. In this paper, the Brownian bridge method and the method used in [13], called 'the sorting and mixing method' herein, is tested for the unscrambled and scrambled Faure sequence, which is a quasi-random sequence in the family of (t,s)-sequences. In section 2 we will briefly introduce the Fokker–Planck equation and its connection to SDEs. In section 3 the quasi-Monte Carlo method is introduced where the concept of low-discrepancy sequence, randomization and effective dimensions are briefly touched upon. The sorting and mixing method and the Brownian bridge method are presented in sections 3.6 and 3.7. In section 4, we derive the SDE for the fast-ion thermalization process, which is followed by the simulation results in section 5 and ends with conclusions.

Consider the Fokker–Planck equation for particle interactions in s coordinates x = x₁,...,x_s,

$$\begin{eqnarray} &&\begin{array}{l} \displaystyle \frac{\partial f}{\partial \tau} = L(f), \quad f : (\mathbf{y},\tau) \in \mathbb{R}^s\times[0,T] \mapsto \mathbb{R},\\[12pt] \displaystyle L= - \sum_i^s \frac{\partial}{ \partial y_i} A_i(\mathbf{y}) +\frac{1}{2} \sum_{i,j}^{s}\frac{\partial^2}{ \partial y_i \partial y_j} B_{ij}(\mathbf{y}), \end{array} \end{eqnarray} \tag{ 2 }$$

with initial condition f(y,0) = δ(x − y) where A_i and $B_{ij} = \sum _{l=1}^s \sigma _{il} \sigma _{jl}$ are the drift vector and diffusion tensor respectively, both independent of time. The 'characteristics' of (2) are described by an Itô SDE

$$\begin{equation} \mathrm{d}\mathbf{X}(t) = A(\mathbf{X}(t))\mathrm{ d}t + \sigma(\mathbf{X}(t))\mathrm{ d}\mathbf{W}(t),\quad \mathbf{X}_0 = \mathbf{x}, \end{equation} \tag{ 3 }$$

where dW(t) is the Wiener process also known as Brownian motion, with normal distributed components having zero mean and dt variance. Many realizations of the SDE give a distribution of particles, which is the solution of (2). A common time-discretization method of the above SDE is the Euler–Maruyama scheme [14]. Let i = {1,...,I} and define the end time by T = IΔt and assume $X \in \mathbb {R}$; then the Euler–Maruyama scheme of (3) is,

$$\begin{equation} X_{i+1} = X_i + A(X_i) (t_{i+1} - t_i )+ \sigma(X_i)( W_{i+1} -W_i), \end{equation} \tag{ 4 }$$

where ΔW_i = W_i+1 − W_i are normally distributed pseudo-random numbers with zero mean and (t_i+1 − t_i)-variance. In the following section we will discuss how the pseudo-random numbers can be replaced with quasi-random sequences.

The naive approach of replacing the pseudo-random numbers with quasi-random numbers will not work for simulation of SDEs since the quasi-random numbers are deterministic and the correlation will introduce artificial drift of the the particles. However, it is possible to use the quasi-random points, if combined with methods that break the correlation.When using the quasi-Monte Carlo method for SDEs the term dimension has a different meaning than the physical number of dimensions. The number of dimensions of a process, X(t), which solves a certain SDE is the number of time steps times the number of physical dimensions. From here onward we let s denote the physical number of dimensions. This is most easily explained by the following example.

3.1. Example

Consider a one-dimensional process X(t) satisfying the scalar SDE dX = σ(X) dW with initial condition X₀ = x. We would like to simulate this SDE for three time steps using the discretization of Euler–Maruyama. Unrolling the SDE we obtain,

$$\begin{eqnarray*} X_0 &=& x,\\ X_1 &=& X_0 + \sigma(X_0)\sqrt{\Delta t} \xi_1 = X_1 (X_0, \xi_1),\\ X_2 &=& X_1(X_0, \xi_1) + \sigma(X_1(X_0, \xi_1))\sqrt{\Delta t} \xi_2 = X_2(X_0, \xi_1, \xi_2),\\ X_3 &=& X_2(X_0, \xi_1, \xi_2) + \sigma(X_2(X_0, \xi_1, \xi_2))\sqrt{\Delta t} \xi_3 = X_3(X_0, \xi_1, \xi_2, \xi_3), \end{eqnarray*}$$

where ξ_i∈N(0,1), are normally distributed random numbers with zero mean and unit variance. Clearly after three time steps the process X is dependent on three random numbers. Now assume that we are interested in calculating an expected value of the process 〈g(X(t₃))〉 for a known function g(·). We know that the normally distributed random numbers ξ_i are related to uniformly distributed numbers z_i∈U(0,1) over the unit interval [0, 1) by the inverse cumulative normal distribution ξ = Φ⁻¹(z). Using this in the above sequence we obtain,

$$\begin{eqnarray*} X_0 &=& x,\\ X_1 &=& X_0 + \sigma(X_0)\sqrt{\Delta t} \Phi^{-1}(z_1)= X_1 (X_0, z_1),\\ X_2 &=& X_1(X_0, z_1) + \sigma(X_1(X_0, z_1))\sqrt{\Delta t} \Phi^{-1}(z_2) = X_2(X_0, z_1, z_2),\\ X_3 &=& X_2(X_0, z_1, z_2) + \sigma(X_2(X_0, z_1, z_2))\sqrt{\Delta t} \Phi^{-1}(z_3) = X_3(X_0, z_1, z_2, z_3). \end{eqnarray*}$$

The expected value of g(X₃) is now a four-dimensional integration,

$$\begin{equation} \langle g(X_3)\rangle = \int_{\mathbb{R}}\int_{[0,1)^3}\delta(x-y) g(x_3(z_1, z_2, z_3;y) )\mathrm{d}z_1\, \mathrm{d}z_2\,\mathrm{d}z_3\,\mathrm{d}y. \end{equation} \noindent \tag{ 5 }$$

From this example we clearly see how the number of dimensions is dependent on the number of physical dimensions and the number of time steps. An estimate of the expected value is obtained by sampling the hypercube with a finite number of particles,

$$\begin{equation} \langle g(X_3)\rangle \approx \frac{1}{N} \sum_{j=1}^N g(x_3(z_1^j, z_2^j, z_3^j) ). \end{equation} \tag{ 6 }$$

If the hypercube is sampled with quasi-random points, we obtain a more accurate estimate of the expected value than if we would use pseudo-random points. This is because the quasi-random points are more uniformly distributed than pseudo-random numbers and do not form clusters, which is a common effect for pseudo-random numbers.

3.2. Low discrepancy sequences

Two common classes of quasi-random sequences are (digital) nets/sequences and lattice rules [15]. In this paper we have used the Faure sequence [15], which is a (0,s)-sequence in base b, where b is the smallest prime number greater than the number of dimensions s. Similar to pseudo-random numbers the law of large numbers must hold for these sequences. A deterministic version of the law of large numbers is provided by the Koksma–Hlawka inequality, which gives a bound on the quadrature error,

$$\begin{equation} \left| \int_{[0,1)^s} g(\mathbf{z}) \mathrm{d}z_1\,\mathrm{d}z_2\cdots \mathrm{d}z_s - \frac{1}{N} \sum_{j=1}^ N g(\mathbf{Z}^j) \right| \leqslant D_N^* V_{\mathrm{HK}}(g), \end{equation} \tag{ 7 }$$

where V_HK(g) is bounded variation of g in the sense of Hardy and Krause [15]. The exact definition of the Hardy and Krause variation contains higher order derivatives of g. Therefore the smoothness of g determines the size of V_HK and is small for smooth functions. D*_N is called the discrete star discrepancy and it measures how uniformly the samples are distributed over the s-dimensional unit cube. A sequence with a small D*_N is called a low-discrepancy sequence. The discrete star discrepancy is defined as,

$$\begin{equation} D^*_N(Z_1,\ldots, Z_N) := \sup \left| \frac{1}{N}\sum_{j=1}^N \mathbb{I}_{0 \leqslant Z_j < p}- \mathrm{vol}([0, p) )\right|, \end{equation} \tag{ 8 }$$

where $\mathbb {I}$ is the indicator function and is 1 if Z_j is inside the box [0,p)^s and zero otherwise. The second term, vol measures the volume of the s-dimensional box [0,p)^s. Equation (8) measures the maximum difference between the number of points in the box and the volume of the box over all possible boxes with one corner in the origin. An illustration of the star discrepancy in two dimensions is given in figure 1.

Zoom In Zoom Out Reset image size

3.3. Randomization

Scrambling is a randomization method for obtaining new uncorrelated (digital) nets or sequences generated from a common mother net/sequence and can partly break the correlation between particles. There are many different randomization methods, the simplest method is so called linear scrambling also known as the generalized Faure for scrambled Faure sequences, [16, 17]. The idea of scrambling is to combine the benefit of calculating sample variance, as in the standard Monte Carlo method, with the improved convergence of the quasi-Monte Carlo method. Scrambling shuffles the position of the points and is achieved by randomizing the digits in the b-adic expansion with a matrix vector multiplication of a matrix with random elements. In short, the nth s-dimensional quasi-random point is generated by,

$$\begin{eqnarray*} &&X_n = [\Phi_{b}(L^{(1)}P^0 \mathbf{n} ), \Phi_{b}(L^{(2)}P^1 \mathbf{n} ), \ldots, \Phi_{b}( L^{(s)}P^{s-1}\mathbf{n} )] \end{eqnarray*}$$

where $\Phi _{b}(\mathbf {n}) = \frac {n_0}{b} + \frac {n_1}{b^2} + \cdots + \frac {n_m}{b^m}$, is the so called radical inverse. The vector n = [n_m,...,n₀] is the b-adic representation of the scalar n, (n = n_mb^m + ··· + n₁b + n₀). The matrix P^s is the sth power of the Pascal matrix modulo b where the (k,l)-element of P is equal to $\left({{l-1}\atop {k-1}}\right){\mathrm { mod}}\,b$. The matrices, L^(s) are independent lower triangular matrices with diagonal elements selected randomly from {1,...,b − 1} and the other elements chosen randomly from {0,...,b − 1}. For a more in depth treatment on randomization we refer the reader to [15, 18–20] (see figure 2).

Zoom In Zoom Out Reset image size

3.4. Effective dimension

High-dimensional problems depend primarily on a few important dimensions, thus reducing the 'effective dimensionality'. This property was analyzed in [21, 22] in terms of the truncated ANOVA decomposition f_d of a function f. The effective dimension can be defined as the lowest possible dimension of f_d for which,

$$\begin{equation} \mathrm{Var}(f_{\mathrm{d}}) \geqslant 0.99 \mathrm{Var}(f). \end{equation} \tag{ 9 }$$

We refer the reader to [23] for a detailed description on effective dimensions.

3.5. Quasi-Monte Carlo for diffusion

The quasi-Monte Carlo method for simulation of diffusion can be used in two ways. In the first concept, the Brownian motion is simulated using an s × I-dimensional sequence. Here I is the number of time steps. The other alternative is to use a 2s-dimensional sequence. Successive realizations of the Brownian motion are obtained by shifting the index pointer to the sequence with N, e.g. use the N first points for the first time-step and use the [N + 1,2N] points for the second time-step and so forth. What is important for both methods is that the realizations should be uncorrelated. In the beginning of this section, we mentioned the importance of breaking the correlation and pointed out that scrambling can partly break the correlation. This is only true in high dimensions. The difference in the auto-correlation between elements of 50 one-dimensional scrambled Faure points and the auto-correlation between elements of one 50-dimensional scrambled Faure point is illustrated in figures 3(a) and (b). The figures show that correlation is very strong between points in the one-dimensional case but not very strong in the 50-dimensional case, which is of the same order as for a standard pseudo-random generator. Scrambling is apparently not sufficient for breaking correlation in the low-dimensional case without using more advanced techniques. In [13], a method for resolving this issue was suggested. It is based on sorting and mixing the particles in a clever way such that decorrelation is ensured, and is the topic of the next section.

Zoom In Zoom Out Reset image size

3.6. The sorting and mixing method

The idea of the sorting and mixing method is to break the correlation by sorting the particles along each dimension. We will give a simplified presentation of the method. For an in depth treatment of the concept we refer the reader to [13]. The foremost drawback of this method is that it requires at least a b^s number of particles being simulated where b is the least prime number greater than 2s. In six dimensions this would give 13⁶ ≈ 4 × 10⁶ number of particles. For higher dimensions this value will explode.

The first step of the sorting and mixing method is to treat the drift and diffusion part of (3) separately,

$$\begin{equation} X_{i+1/2} = X_i + A(X_i) (t_{i+1} - t_i), \end{equation} \tag{ 10 }$$

$$\begin{equation} X_{i+1} = X_{i+1/2} + \sigma(X_i)( W_{i+1} -W_i ). \end{equation} \tag{ 11 }$$

The drift (10) is evaluated with the standard forward Euler method. For the diffusion (11) part we need to generate quasi-random numbers for the Wiener process. Let N = b^m be the number of particles to be simulated. Here b is the least prime number greater than 2s, e.g. for s = 2⇒b = 5, and let d₁,...,d_s be integers greater than zero such that m = d₁ + ··· + d_s. For a two-dimensional diffusion process a minimum of N = 5² particles are needed.

• Sorting. First sort the particles in ascending order of the magnitude of the first coordinate of the particles into b^d₁ groups. Secondly sort the particles in each group in ascending order of the magnitude of the second coordinate into b^d₂ subgroups. For a three-dimensional problem sort the particles in each b^d₂ subgroup in ascending order of the magnitude of the third coordinate of the particles into b^d₃ sub-subgroups. Continue sorting if dimensions are greater than three. This is a computationally demanding operation requiring $\mathcal {O}(N \log N)$ operations (with quick-sort) for each sort.
• Mixing. Consider a vector y = (y₁,...,y_2s) from a (0,2s)-net in base b and define two selection functions P' and P'' by P'y: = (y₁,...,y_s) and P''y: = (y_s+1,...,y_2s). Use the first s dimension of the (0,2s)-net for mixing. Let a(y): = (⌊ b^d₁y₁ ⌋,...,⌊ b^d_sy_s) ⌋) where the floor function, ⌊ · ⌋, returns the greatest integer less than or equal to the argument. Define the one-to-one map between array indices using the first s variables of the (0,2s)-net, j → a(P'y_iN+j). The map is from [0,N) into [0,b^d₁) × ··· × [0,b^d_s), e.g. in two dimensions we map the vector [0,N) to the matrix [0,b^d₁) × [0,b^d₂). The particle with array index a(P'y_iN+j) is pushed by the quasi-random points, P''y_iN+j.

The sorting and mixing procedures are illustrated in figures 4(a) and (b) for s = 2. It was proved in [13] that this construction has the following bound:

$$\begin{eqnarray} &&\fl D_N^*(X_i^j)\leqslant D_N^*(X_0^j) + b^{d_1 + \cdots + d_{s-1} + \lfloor d_s/2 \rfloor}\sum_{k=0}^{i-1} D_N^*(Y^k) + i \left( \frac{1}{b^{d_1}} + \cdots + \frac{1}{b^{d_s -1}} + \frac{1}{b^{\lfloor d_s/2 \rfloor}} \right), \end{eqnarray} \tag{ 12 }$$

where 0 ⩽ j < N and D*_N(Y^k) is the star discrepancy of the (0,2s)-net, Y^k in base b. Since N is defined as the product of b^d_i, the second term converges in the limit of infinite number of particles. Note that for a finite N the bound increases with the number of time steps, i. To keep the error fixed, more particles are needed as the simulation time increases.

Zoom In Zoom Out Reset image size

3.7. The Brownian bridge method

In this section we will turn our attention to the s × I-dimensional case. If many dimensions are used, the number of effective dimensions become important. It has been argued in [12, 21], that the midpoint Brownian bridge method effectively reduces the effective dimensions. The general Brownian bridge formula is defined as

$$\begin{equation} W_{t_j} = (1-\lambda)w_{t_i} + \lambda w_{t_k} + \left( \lambda (1- \lambda ) (k-i) \Delta t \right)^{1/2} Z, \end{equation} \tag{ 13 }$$

where λ = (j − i)/(k − i) for t_i < t_j ⩽ t_k and Z is a normally distributed random number. The variables w_ti,w_tk are previous realization of the Wiener process. The midpoint Brownian bridge formula is obtained from the special case where the time axis is populated by successive midpoint splits. First, generate the start- w₀ and end-point w_T . Then draw the midpoint and continue splitting, see figure 5. The generated order is: w₀,w_T ,w_T/2,w_T/4,w_3T/4,...,w_T−1. One interesting note is that the midpoint Brownian bridge construction reduces the variance of W_tj by a factor of 2 for each split,

$$\begin{eqnarray*} \mathrm{Var}[W_{t_{k/2}}] &=& 1/4 k \Delta t, \quad t_0 < t_{k/2} < t_{k},\\ \mathrm{Var}[W_{t_{k/4}}] &=& 1/4 k/2 \Delta t , \quad t_0 < t_{k/4} < t_{k/2}.\\ &\vdots& \end{eqnarray*}$$

The success of the Brownian bridge method can be explained by the way it utilizes low-discrepancy sequences. An important property of low-discrepancy sequences is that they are typically more uniformly distributed in lower dimensions than in higher. This is efficiently used by the Brownian bridge construction, since the first few quasi-random points describe the complete Brownian path on a coarse level. On the finer scale the reduction of variance in the Brownian bridge construction, reduces the impact of the degradation in the low-discrepancy sequence. It should be noted that the Brownian bridge method does not always perform better than the standard Monte Carlo method. For example in [24] it is shown that the performance of the Brownian bridge construction is not consistent and can actually perform worse than the standard Monte Carlo method for certain types of integrands. Thus the Brownian bridge method is problem dependent and has the potential to generate many papers.

Zoom In Zoom Out Reset image size

In this section we present an SDE for a simple model of fast-ion thermalization in fusion plasmas given by the Fokker–Planck equation with the Coulomb collision operator of Spitzer [25].

$$\begin{eqnarray} &&\fl \frac{\partial f}{\partial t} = \!-\!\frac{1}{v^2}\frac{\partial }{\partial v}\left[ v^2 \langle \Delta v_{\Vert} \rangle + \frac{v}{2}\langle (\Delta v_{\perp})^2\rangle \right]f + \frac{1}{2v^2}\frac{\partial^2}{\partial v^2} \left[ v^2 \langle (\Delta v_{\Vert})^2\rangle\right]f + \frac{1}{4 v^2} \frac{\partial }{\partial \xi} (1\!-\!\xi^2)\frac{\partial }{\partial \xi}[ \langle (\Delta v_{\perp})^2\rangle ] f. \end{eqnarray} \tag{ 14 }$$

The above equation is not of the form of (2) due to the Jacobian v², therefore the drift A and diffusion coefficients σ are calculated from the moments of the the Coulomb operator,

$$\begin{eqnarray*} &&\fl C(\cdot) = -\frac{1}{v^2}\frac{\partial }{\partial v}\left[ v^2 \langle \Delta v_{\Vert} \rangle + \frac{v}{2}\langle (\Delta v_{\perp})^2\rangle \right] + \frac{1}{2v^2}\frac{\partial^2}{\partial v^2} \left[ v^2 \langle (\Delta v_{\Vert})^2\rangle\right] + \frac{1}{4 v^2} \frac{\partial }{\partial \xi} (1-\xi^2)\frac{\partial }{\partial \xi}[ \langle (\Delta v_{\perp})^2\rangle ]. \end{eqnarray*}$$

The first two moments of the above operator are given by,

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \langle v \rangle = \alpha(v), \end{equation} \tag{ 15 }$$

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \langle \xi \rangle = - \frac{1}{2 v^2 } \xi \gamma(v), \end{equation} \tag{ 16 }$$

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \langle v^2 \rangle = 2 v \alpha(v) +\beta(v), \end{equation} \tag{ 17 }$$

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \langle \xi^2 \rangle = \frac{1}{2v^2} (1-3\xi^2)\gamma(v), \end{equation} \tag{ 18 }$$

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \langle v\xi \rangle = \xi \langle \Delta v_{\Vert} \rangle , \end{equation} \tag{ 19 }$$

where we have introduced the variables, following [25], $\alpha (v) = \langle \Delta v_{\Vert} \rangle + \frac {1}{2v}\langle (\Delta v_\perp )^2\rangle$, β(v) = 〈(Δv_∥)²〉, γ(v) = 〈(Δv_⊥)²〉. The explicit expressions for the moments are,

$$\begin{eqnarray*} &&\begin{array}{l} \displaystyle \langle \Delta v_{\Vert} \rangle = \sum_f C_f v_f^{-2}[ 1 + m/m_f]G( v/v_f), \\[12pt] \displaystyle \langle (\Delta v_{\parallel})^2\rangle = \sum_f \frac{C_f}{v} G(v/v_f), \\[12pt] \displaystyle \langle (\Delta v_{\perp})^2\rangle = \sum_f \frac{C_f}{v} [ \mathrm{erf}(v/v_f) - G(v/v_f) ], \end{array} \end{eqnarray*}$$

where v_f is the thermal velocity, C_f is given in [25] and f is the index over the species of the field ions. The function G is given by,

$$\begin{eqnarray*} &&G(x) = \frac { 1}{2x^2}\left(\mathrm{erf}(x) - 2x \pi^{-1/2} \,\mathrm{e}^{-x^2} \right). \end{eqnarray*}$$

Armed with the equations above, the covariance matrix is calculated from the following general formula,

$$\begin{equation} \sigma^2_{ij} = \frac{\mathrm{d}}{\mathrm{d}t} \langle X_i X_j\rangle - X_i \frac{\mathrm{d}}{\mathrm{d}t}\langle X_j \rangle- X_j \frac{\mathrm{d}}{\mathrm{d}t} \langle X_i\rangle . \end{equation} \tag{ 20 }$$

Inserting (15)–(19) in the equation above we obtain the following:

$$\begin{equation} \begin{array}{l} \displaystyle \sigma^2_{11} = \beta(v), \\[12pt] \displaystyle \sigma^2_{12} = \sigma^2_{21} = 0,\\[12pt] \displaystyle \sigma^2_{22} = \frac{\gamma(v)}{2v^2}( 1 - \xi^2) . \end{array} \end{equation} \tag{ 21 }$$

Since the nonzero elements of the covariance matrix are on the matrix diagonal, the diffusion coefficient σ is immediately obtained. The elements of the drift vector are given by,

$$\begin{equation} \begin{array}{l} \displaystyle A_1 = \alpha(v), \\[12pt] \displaystyle A_2 = - \frac{1}{2 v^2 } \xi \gamma(v). \end{array} \end{equation} \tag{ 22 }$$

The drift and diffusion coefficients are in v,ξ coordinates. Using the Itô formula for E = g(V ) = mV²/2 we obtain an SDE in energy and pitch-angle,

$$\begin{equation} \begin{array}{l} \displaystyle \mathrm{d}E =\left( \alpha m^{1/2} (2E)^{1/2} + m/2 \beta \right)\mathrm{d} t + \left( 2mE \beta \right)^{1/2} \mathrm{d}W,\\[12pt] \displaystyle \mathrm{d}\xi = - \frac{m\gamma }{4E}\xi \mathrm{d}t + \left(\frac{m \gamma}{4E}( 1 - \xi^2) \right)^{1/2} \mathrm{d}W. \end{array} \end{equation} \tag{ 23 }$$

This system is solved with the Euler–Maruyama scheme,

$$\begin{equation} \begin{array}{l} \displaystyle E_{n+1} = E_n + \left( \alpha(E_n) m^{1/2} (2E_n)^{1/2} +m/2\beta(E_n) \right)\Delta t + \left( 2m E_n \beta \right)^{1/2} \sqrt{\Delta t} Z_1, \\[12pt] \displaystyle \xi_{n+1} = \xi_n - \frac{m \gamma(E_n)}{4 E_n } \xi_n \Delta t + \left( \frac{m \gamma(E_n)}{4E_n } \left(1 - (\xi_n)^2\right) \right)^{1/2} \sqrt{\Delta t} Z_2, \end{array} \end{equation} \tag{ 24 }$$

where Z₁,Z₂ are normal distributed random numbers with zero mean and unit variance.

This section concerns modeling the thermalization of fast ions from neutral beam injection in fusion devices. In practice these models often aim at determining the heat and momentum transfer from the fast ions to the thermal plasma components, along with the fast ion current drive. For this purpose it is often sufficient to solve a linearized Fokker–Planck equation for the fast ions, in which the thermal components of the plasma are considered stationary, or quasi-stationary [25]. The separation of fast and thermal ions can be done in different ways, e.g. in the test particle Monte Carlo codes NUBEAM [26] and ASCOT [27] where the test particles are dropped when they reach 3/2 of the thermal energy. In this work the separation is achieved by simulating the test particles for half a collisional slowing down time against the electrons. Thus, the energy distributions calculated here differ slightly from those obtained with NUBEAM and ASCOT in the thermal energy range, but not for the supra-thermal energies. Furthermore, by assuming the background plasma to be stationary the equation becomes autonomous. This allows us to use the technique used in [27], where the steady state distribution function is proportional to the time that a set of test particles spend in each volume element.

Four methods have been tested, the standard Monte Carlo, the sorting and mixing method, the Brownian bridge method and a naive method where the pseudo-random numbers are replaced with Faure numbers, scrambled and unscrambled without any extra modification. Measures of the run times of the simulations are not considered since the implementation of the Faure quasi-random generator has not been optimized in contrast to the pseudo-random generator, which is highly optimized in Matlab^TM. However, the sorting mixing method is slower than the Brownian bridge method since it requires sorting the particles each time step. The plasma is assumed to consist of protons and electrons, both with densities 3 × 10¹⁹ 1 m⁻³ and temperatures 4 keV. The beam ions are protons injected with 100 keV energy and a pitch angle ξ = −0.8, where ξ = v_∥/v, v is the speed and v_∥ is the velocity component along the magnetic field. The fast protons will collide primarily with the electron above 59 keV and with the thermal protons below. The total simulation time is half a slowing down time and is split into 2⁷ time steps. The parameters for the sorting and mixing method are, d₁ = {1,2,2,2,2}, d₂ = {1,1,2,3,4} and b = 5 for N = b^(d₁+d₂). The Faure sequence was scrambled to a depth of 30-digits. The tests was conducted in Matlab^TM version R2010b on a twelve core Intel machine @ 2.67 GHz with 62 GB of RAM. Since the entire particle trajectory was saved, a massive amount of memory was required with a maximum of 5⁹ × 2¹⁴ = 32 × 10⁹ particles. The exact solution for this test case is unknown. Therefore we have measured the convergence in terms of sample variance on a batch of 20 runs, each simulation with different initial seed of the randomly scrambled Faure sequence. Two different sample-means were considered, the parallel velocity $\langle g(E, \xi )\rangle = \langle v_{\Vert}\rangle = \langle v \xi \rangle =\langle \xi \sqrt {2E/m}\rangle$, which is related to the fast-ion current drive and the energy squared 〈g(E)〉 = 〈E²〉. The batch-mean value is calculated from,

$$\begin{equation} \hat{\mu} = \frac{1}{M N} \sum^M_{j=1} \sum^N_{i=1} g(E_i^j, \xi_i^j) \end{equation}\noindent \tag{ 25 }$$

and the variance of the batch-mean is given by,

$$\begin{equation} \hat{\sigma}^2 = \frac{1}{M-1} \sum_{j=1}^M ( \hat{\mu}_j - \hat{\mu})^2. \end{equation} \tag{ 26 }$$

From these expressions we have calculated the 99% confidence interval for the batch-mean from the Student-distribution with M − 1 degrees of freedom,

$$\begin{equation}\hat{\mu} \pm t_{99\%, M-1} \sqrt{\hat{\sigma}^2/M}, \end{equation}\noindent \tag{ 27 }$$

where t_99%,M−1 = 2.86. Simulation results for the parallel velocity for the test cases are given in table 1 and for the energy squared test case in table 2 .

The convergence of the normalized root mean square error, $\mathrm {n.r.m.s} = \sqrt {\hat {\sigma }^2}/( \hat {\mu }\sqrt {M} )$ for the two test cases is illustrated in figures 6(a) and (b). In these figures, BB denotes the Brownian bridge method and SM the method of sorting and mixing. The histograms of the stationary distributions for the different methods are plotted in figures 7(a)–10(b). A reference distribution, obtained from a standard Monte Carlo simulation, is plotted in figure 7(a). In figure 7(b), the distribution of the naive method with unscrambled Faure points is plotted. The method does not converge as expected since the correlation in the Faure sequence introduces spurious drifts of the particles. The result is slightly improved if combined with the Brownian bridge method, see figure 8(a), but the distribution does not agree with the reference figure 7(a). The sorting and mixing method is the only method with a theoretical proof of the convergence for unscrambled Faure points. The distribution for the sorting and mixing method is plotted in figure 8(b) and agrees with the standard Monte Carlo distribution. The distribution from a simulation with the scrambled Faure sequence is plotted in figure 9. From this figure we can see that the scrambling reduces the error from the correlation. Even though the distribution of the naive scrambled Faure method is in agreement with the Monte Carlo reference distribution, the moments are not in agreement for the case 〈vξ〉 as seen in table 1. The distribution for the scrambled Faure with the Brownian bridge method is plotted in figure 10(a) and it is in agreement with the Monte Carlo reference. The distribution for the sorting and mixing method is plotted in figure 10(b) and agrees with the Monte Carlo reference. The sorting and mixing method and the Brownian bridge method perform much better than standard Monte Carlo for 2⁷ time steps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1. Very long-time simulation

In order to test the applicability for long-time particle simulations we have estimated the slope of the convergence with a least-squares fit for different number of time steps. The result is plotted in figure 11. The results indicate that the sorting and mixing method performs worse than the standard Monte Carlo above 400 time steps when measuring parallel velocity and remains very low for 〈E²〉 for all tested values. The convergence of the Brownian bridge method degrades almost linearly for both measured quantities, but the convergence is better than the standard Monte Carlo even for a thousand time steps. A rough extrapolation of the slope for the Brownian bridge method indicates that it will perform worse than the standard Monte Carlo at about 1700 time steps.

Zoom In Zoom Out Reset image size

In summary, we have tested the methods for a very large number of time steps, I = 2¹⁴ = 16384. The computational results are given in tables 3 and 4. For this value the number of dimensions of the hypercube is 2¹⁵ for the Brownian bridge method. The theoretical number of particles required to uniformly cover the hypercube is astronomical, therefore we can not expect the method to perform well. The sorting and mixing method cannot perform well either because the bound of the discrepancy (12) increases with the number of time steps, the sum in (12) is of the order of b^{d₁+⌊ d₂/2 ⌋} × I × D*(Y⁰) ≈ 10⁷ × D*(Y⁰). The Brownian bridge method only gives an acceptable estimate of 〈E²〉 for N = 5⁶ but has the same confidence interval as the standard Monte Carlo, table 4. The performance is slightly better for 〈vξ〉 where the value is estimated correctly except for N = 5² but the confidence interval is smaller than the standard Monte Carlo at N = 5⁶, table 3. The sorting and mixing method does not converge for the parallel velocity case, as seen in table 3 and it gives an inconsistent estimate for the energy squared case, table 4.

Table 3. Computational results for 〈v_||〉 with the Brownian bridge method and the sorting and mixing method at dt = 2⁻¹⁵.

N	Monte Carlo	Brownian bridge	Sorting and mixing
52	−2.617±0.096	−1.859±0.174	−2.716±0.003
53	−2.604±0.025	−2.478±0.120	−2.681±0.004
54	−2.590±0.011	−2.545±0.034	−2.635±0.005
55	−2.588±0.005	−2.577±0.008	−2.600±0.005
56	−2.588±0.003	−2.585±0.001

Table 4. Computational results for 〈E²〉 with the Brownian bridge method and the sorting and mixing method at dt = 2⁻¹⁵.

N	Monte Carlo	Brownian bridge	Sorting and mixing
52	5.171±0.184	2.112±0.208	4.958±0.014
53	5.128±0.081	4.161±0.358	5.019±0.009
54	5.123±0.040	4.882±0.077	5.082±0.006
55	5.128±0.018	5.066±0.026	5.120±0.002
56	5.123±0.007	5.113±0.007

We have tested the applicability of the quasi-Monte Carlo method on fast-ion thermalization using a simplified neutral beam injection model. The method of sorting and mixing and the Brownian bridge method are much better than the standard Monte Carlo method with a faster convergence and more accurate estimate up to a thousand time steps for simulation of fast-ion thermalization. When very long simulations are required these methods can fail to give accurate estimates and they may not converge. For modest number of time steps (2⁷) the value of 〈vξ〉 converges as $\mathcal {O}(N^{-1})$ for the Brownian bridge method and as $\mathcal {O}(N^{-0.6})$ for the sorting and mixing method, while for 〈E²〉 both methods converge as $\mathcal {O}(N^{-1})$. The sorting and mixing method is more accurate than the Brownian bridge method but is also more computationally demanding due to the sorting stage and requires simultaneous simulation of an ensemble of particles. The sorting and mixing method is therefore not suitable for higher-dimensional problems s > 3 . The Brownian bridge method has a similar convergence to the sorting and mixing method and can be used for evolution of single particles, but the method can converge to the wrong value for a large number of time steps with few particles and depends strongly on the measured quantity. The Brownian bridge method and the sorting and mixing method combined with scrambled Faure are very efficient within their working domain.