Eliminating poor statistics in Monte-Carlo simulations of fast-ion losses to plasma-facing components and detectors

TL;DR

This paper introduces a Backward Monte Carlo method for simulating fast-ion losses in plasma devices, which overcomes the poor statistical issues of traditional Forward Monte Carlo approaches, leading to more accurate and efficient predictions.

Contribution

The paper presents the BMC algorithm tailored for fast-ion loss simulations, detailing its implementation and advantages over the conventional FMC method.

Findings

BMC eliminates poor statistics in small detector domains.

The method provides more accurate loss signal simulations.

Numerical examples demonstrate improved efficiency.

Abstract

With Wendelstein 7-X now up and running, and the construction of ITER proceeding, predicting fast-ion losses to sensitive plasma-facing components and detectors is gaining significant interest. A common recipe to perform such studies is to push a large population of marker particles along their equations of motion, the trajectories randomized with Monte Carlo operators accounting for Coulomb collisions, and to record possible intersections of the marker trajectories with synthetic detectors or areas of interest in the first wall. While straightforward to implement and easy to parallelize, this Forward Monte Carlo (FMC) approach tends to suffer from poor statistics and error estimation as the detector domain is often small: it is difficult to guess how to set up the initial weights and locations of the markers for them to remain representative of the source distribution, yet record enough hits to the detector for good statistics. As an alternative, the FMC method can be replaced with a so-called Backward Monte Carlo (BMC) algorithm. Instead of starting with a given initial marker population, one starts from the end condition at the detector and records how the hit probability evolves backwards in time. The scheme eliminates the statistics issue present in the FMC scheme and may provide more accurate and efficient simulations of fast-ion loss signals. The purpose of this paper is to explain the BMC recipe in the fast-ion setting and to discuss the associated nuances, especially how to negate artificial diffusion. For illustration purposes, our numerical example considers a 1-D stochastic harmonic oscillator as a mock-up of a charged particle.

Figures

Equations

$$\mathbb{E}[Z_{t+\Delta t}=z_{i-1}|Z_t=z_i]$$

$$\mathbb{E}[Z_{t+\Delta t}=z_i|Z_t=z_i]$$

$$\mathbb{E}[Z_{t+\Delta t}=z_{i+1}|Z_t=z_i]$$

$$\Phi (z_i,m)=\mathbb{E}[Z_{M\Delta t}=y|Z_{m\Delta t}=z_i],$$

$$\Phi (x,m-1)$$

$$=\sum _{z_i\in \Omega }\mathbb{E}[Z_{M\Delta t}=y|Z_{m\Delta t}=z_i]\mathbb{E}[Z_{m\Delta t}=z_i|Z_{(m-1)\Delta t}=x]$$

$$=\sum _{z_i\in \Omega }\Phi (z_i,m)\mathbb{E}[Z_{m\Delta t}=z_i|Z_{(m-1)\Delta t}=x]$$

$$\equiv \mathbb{E}[\Phi (Z_{m\Delta t},m)|Z_{(m-1)\Delta t}=x],$$

$$d{\mathbf{X}}_t$$

$$d{\mathbf{V}}_t$$

$$U_{\alpha }^{\mathbf{x}}(\mathbf{z},t)$$

$$U_{\alpha }^{\mathbf{v}}(\mathbf{z},t)$$

$$\mathbf{z}(t)\mapsto {\phi }_{U_{\alpha }}^{s,t}(\mathbf{z}(t))$$

$$\displaystyle\frac{d\varphi^{s,t}_{U_{\alpha}}(\bm{z})}{ds}\Big{|}_{s=0}$$

$${\mu }_{\alpha }^i(\mathbf{x},\mathbf{v},t)$$

$${\sigma }_{\alpha }^{ij}(\mathbf{x},\mathbf{v},t)$$

$${\nu }_{\alpha }$$

$$D_{\alpha \parallel }$$

$$D_{\alpha \perp }$$

$$\mathbf{Z}\mapsto {\phi }_{U_{\alpha }}^{\Delta t,t}(\mathbf{Z}),$$

$$(\mathbf{X},\mathbf{V})\mapsto {\phi }_C^{\Delta t,t}(\mathbf{X},\mathbf{V})=(\mathbf{X},\mathbf{V}+\mathbf{\mu }\left(\mathbf{Z},t\right)\Delta t+\mathbf{\sigma }\left(\mathbf{Z},t\right)\cdot \mathbf{\xi }\sqrt{\Delta t}),$$

$$I(\Omega ,t)=\int S({\mathbf{z}}^{\prime })\mathbb{E}[P({\mathbf{Z}}_{\tau };\Omega ,t)|{\mathbf{Z}}_0={\mathbf{z}}^{\prime }]d{\mathbf{z}}^{\prime }.$$

$$\Phi (\mathbf{z},m)=\mathbb{E}[P({\mathbf{Z}}_{\tau };\Omega ,M\Delta t)|{\mathbf{Z}}_{m\Delta t}=\mathbf{z}],\tau \in [m\Delta t,M\Delta t],$$

$$\Phi (\mathbf{z},m-1)$$

$$=\int \mathbb{E}[P({\mathbf{Z}}_{\tau };\Omega ,M\Delta t)|{\mathbf{Z}}_{m\Delta t}={\mathbf{z}}^{\prime }]\mathbb{E}[{\mathbf{Z}}_{m\Delta t}={\mathbf{z}}^{\prime }|{\mathbf{Z}}_{(m-1)\Delta t}=\mathbf{z}]d{\mathbf{z}}^{\prime }$$

$$=\int \Phi ({\mathbf{z}}^{\prime },m)\mathbb{E}[{\mathbf{Z}}_{m\Delta t}={\mathbf{z}}^{\prime }|{\mathbf{Z}}_{(m-1)\Delta t}=\mathbf{z}]d{\mathbf{z}}^{\prime }$$

$$\equiv \mathbb{E}[\Phi ({\mathbf{Z}}_{m\Delta t},m)|{\mathbf{Z}}_{(m-1)\Delta t}=\mathbf{z}].$$

$$\mathbf{z}\mapsto {\phi }^{\Delta t,(m-1)\Delta t}(\mathbf{z},\mathbf{\xi }).$$

$$\Phi (\mathbf{z},m-1)={\int }_{{\mathbb{R}}^3}\Phi ({\phi }^{\Delta t,(m-1)\Delta t}(\mathbf{z},\mathbf{\xi }),m)\frac{\exp [-{\mathbf{\xi }}^2/2]}{(2\pi {)}^{3/2}}d\mathbf{\xi }.$$

$$d{X}_t$$

$$d{V}_t$$

$$(x,v)\mapsto {\phi }_H^{\Delta t}(x,v)$$

$$(x,v)\mapsto {\phi }_S^{\Delta t}(x,v)$$

$$(x,v)\mapsto {\phi }_S^{\Delta t}\circ \underset{n}{\underbrace{{\phi }_H^{\Delta t/n}\circ \cdots \circ {\phi }_H^{\Delta t/n}(x,v)}}\equiv {\phi }^{\Delta t}(x,v,\xi ),$$

Full text

Introduction – To get a feeling for Backward Monte Carlo (BMC) modeling, we start with a one-dimensional random process $Z_{t}\in\{z_{i}\}_{i=1}^{N}$ that has discrete outcomes depending only on the previous value of the process. We shall assume the process $Z_{t}$ to satisfy the following "go left – stay put – go right" jump probabilities on the $z$-axis

[TABLE]

Then we ask: Given any $Z_{0}=x\in\{z_{i}\}_{i=1}^{N}$, at what probability does $Z_{M\Delta t}=y\in\Omega$ happen? In mathematical terms, this is equal to computing the expectation $\mathbb{E}[Z_{M\Delta t}=y|Z_{0}=x]$. The Forward Monte Carlo (FMC) recipe to solve this problem involves populating each $z_{i}$ with multiple marker processes, simulating each marker according to the jump probabilities given above, and, after $M$ time steps, checking whether the marker has the value $y$. The FMC algorithm is straightforward to implement on a computer and trivial to parallelize as each marker process is independent of each other. The method might, however, become computationally expensive if a strict error estimate is needed – the error is proportional to $1/\sqrt{N_{i}}$ with $N_{i}$ the number of markers launched from site $z_{i}$ – or if the set $\{z_{i}\}_{i=1}^{N}$ is very large.

There is, of course, a simpler solution to the problem. Let us define the function

[TABLE]

and note that $\Phi(x,0)=\mathbb{E}[Z_{M\Delta t}=y|Z_{0}=x]$ would give us the expectation value we are chasing and that $\Phi(x,M)=\mathbf{1}_{y}(x)$ is the indicator function. Using conditional probabilities, we compute a recursive expression for $\Phi(x,m)$, namely

[TABLE]

valid for any $0<m\leq M$. Updating the function $\Phi(x,m)$, by letting $m$ run from $M$ down to $1$, constitutes to computing how the probability to reach the end condition branches out "backwards" in time or, put more precisely, how the probability depends on the time given for the markers to reach the end condition when this time gradually increases. For example, let us choose $M=3$, an end condition $Z_{3\Delta t}=2$, and the set $\{z_{i}\}_{i=1}^{4}=\{1,2,3,4\}$. The resulting branching is illustrated below, with bolded numbers denoting the probability (or $\Phi(z,m)$) of each state to reach the end condition, during each of the three rounds of iterations.

z$$t$$0\Delta t$$1\Delta t$$2\Delta t$$3\Delta t$$1$$2$$3$$4[math]$1$[math][math]\bf{\tfrac{1}{2}}$$\bf{\tfrac{1}{3}}$$\bf{\tfrac{1}{3}}$$\bf{0}

z$$t$$0\Delta t$$1\Delta t$$2\Delta t$$3\Delta t$$1$$2$$3$$4[math]$1$[math][math]\tfrac{1}{2}$$\tfrac{1}{3}$$\tfrac{1}{3}[math]\bf{\tfrac{5}{12}}$$\bf{\tfrac{7}{18}}$$\bf{\tfrac{2}{9}}$$\bf{\tfrac{1}{6}}

z$$t$$0\Delta t$$1\Delta t$$2\Delta t$$3\Delta t$$1$$2$$3$$4[math]$1$[math][math]\tfrac{1}{2}$$\tfrac{1}{3}$$\tfrac{1}{3}[math]\tfrac{5}{12}$$\tfrac{7}{18}$$\tfrac{2}{9}$$\tfrac{1}{6}$$\bf{\tfrac{29}{72}}$$\bf{\tfrac{37}{108}}$$\bf{\tfrac{7}{27}}$$\bf{\tfrac{7}{36}}

Obviously, there is no uncertainty in determining $\Phi(x,m)$ but, where the FMC method is trivially parallel, the BMC algorithm does require a gather operation of the information $\Phi(z_{i},m)$ after each time step. Fortunately, for physical processes, the jump probability $\mathbb{E}[Z_{t+\Delta t}=z_{j}|Z_{t}=z_{i}]$ is often peaked and effectively limited to a narrow range. Hence, even if the set $\{z_{i}\}_{i=1}^{N}$ grew significantly, domain decomposition methods would likely be effective.

Before moving on, we make an additional note. In the example described above, we have used a jump probability that depends only on the previous value of the process. In the iterative expression (Eliminating poor statistics in Monte-Carlo simulations of fast-ion losses to plasma-facing components and detectors), the expectation $\mathbb{E}[Z_{m\Delta t}=y|Z_{(m-1)\Delta t}=x]$ could, however, also depend on the time index $m$, accounting, e.g., for temporal changes in conditions the jump probabilities depend on. With this introduction to BMC calculation, we shall apply the principle to model propabilities of fast ions in tokamak or stellarator plasmas to reach preset phase-space domains such as specific wall tiles or particle detectors.

Traditional modeling of fast-ion-loss signals: – Arguably, the Vlasov equation, supplemented with the Landau collision operator, is a sufficient model to describe the distribution function of charged particles in many magnetized plasmas, such as those in tokamaks and stellarators. While the Landau operator is non-linear, in case of fast ions that form a minority population in typical plasmas, the so-called field-particle distribution in the collision operator can be approximated to be a Maxwellian and the self-collisions neglected. The kinetic equation for fast ions can then be assumed linear.

The resulting, so-called test-particle kinetic equation is a Fokker-Planck equation with a well-known connection to continuous stochastic processes. In a given background plasma, an individual marker particle with coordinates $\bm{Z}_{t}=(\bm{X}_{t},\bm{V}_{t})$, sampled from the distribution function of species $\alpha$, is subject to the following stochastic differential equation (SDE) of Itô kind

[TABLE]

The vector field $U_{\alpha}(\bm{z},t)=(U_{\alpha}^{\bm{x}}(\bm{z},t),U_{\alpha}^{\bm{v}}(\bm{z},t))$ has the components

[TABLE]

with $\bm{E}$ and $\bm{B}$ the electric and magnetic fields and $e_{\alpha}$ and $m_{\alpha}$ the charge and mass of species $\alpha$, and represents time-dependent Hamiltonian flow

[TABLE]

The vector $\bm{\mu}_{\alpha}$ and the matrix $\bm{\sigma}_{\alpha}$ appearing in the SDE for $\bm{Z}_{t}$ originate from the test-particle collision operator, with Cartesian components given by

[TABLE]

and the coefficients $\nu_{\alpha}$, $D_{\alpha\parallel}$, and $D_{\alpha\perp}$ being

[TABLE]

The argument $y_{\beta}$, at which the function $\Psi(y)=\pi^{-1/2}\exp(-y^{2})+(y+1/(2y))\text{erf}(y)$ and its derivatives are to be evaluated, is the particle speed normalized to the thermal velocity, i.e., $y_{\beta}=\sqrt{m_{\beta}|\bm{v}|^{2}/2T_{\beta}}$, with $n_{\beta}(\bm{x},t)$ and $T_{\beta}(\bm{x},t)$ the density and temperature of species $\beta$. The coefficient $c_{\alpha\beta}=e_{\alpha}^{2}e_{\beta}^{2}\ln\Lambda_{\alpha\beta}/(8\pi\varepsilon_{0}^{2})$ involves the Coulomb logarithm $\ln\Lambda_{\alpha\beta}$ and, finally, $\bm{W}_{t}\sim\mathcal{N}(\bm{0},t\mathbf{I})$ is the standard multivariate Wiener process. In essence, $\nu_{\alpha}$ acts as a friction coefficient, slowing the particle down, whereas $D_{\alpha\parallel}$ results in velocity diffusion along the particle velocity and $D_{\alpha\perp}$ in the direction perpendicular to it.

The above set of equations and expressions describes how a test particle with phase-space coordinates $\bm{Z}_{t}=(\bm{X}_{t},\bm{V}_{t})$ behaves in a plasma, including the fast Larmor gyration around the magnetic field line. One could write an analogous set of equations describing the guiding-center motion of test particles and, also, further approximate the expressions for $\bm{\mu}$ and $\bm{\sigma}$ to reflect the speed of the fast ions with respect to a thermal plasma.

Out in the field, deep in the trenches, several numerical codes Kramer et al. (2013); White et al. (2010); Akers et al. (2012); Hirvijoki et al. (2014); Pfefferlé et al. (2014); Shinohara et al. (2016); Reusch et al. (2014); Clauser et al. (2019); Carbajal et al. (2017); Xu et al. (2018) use the equations above – or their guiding-center counterparts – to perform FMC simulations of various fast ion (or electron) signals such as the energy flux onto certain wall tiles or the distribution of fast ions recorded by a fast-ion-loss detector (FILD). In practice, each code implements the following algorithm: (i) Assign a marker particle with initial coordinates $\bm{Z}_{0}=(\bm{X}_{0},\bm{V}_{0})$. (ii) Iterate the phase-space position of the marker by advancing the Hamiltonian flow along the vector field $U_{\alpha}$ with the map

[TABLE]

and by updating the marker velocity to account for Coulomb collisions with the Euler-Maryama scheme

[TABLE]

where $\bm{\xi}\sim\mathcal{N}(\bm{0},\mathbf{I})$ are random numbers obeying standard multivariate normal distribution. (iii) If the trajectory intersects a wall tile or a diagnostic, represented by a phase-space domain $\Omega$, record a hit, accumulate the recorded quantity, and terminate the marker if necessary. In this process, a large number of marker particles is used to populate the phase-space, and often only a few hit the detector, resulting in poor statistics.

BMC approach: – Fast-ion signals could be equivalently computed from the expression

[TABLE]

Here $I$ is the signal strength, $\Omega$ is the domain from which the signal is produced, and $t$ is the time interval over which the marker particles are allowed to evolve from their initial state and possibly intersect the domain $\Omega$. The signal is then computed as a sum of the source particles, represented by $S(\bm{z})$, weighted with the probability $\mathbb{E}[P(\bm{Z}_{\tau};\Omega,t)|\bm{Z}_{0}=\bm{z}]$ for a marker from $\bm{z}$ to reach the domain $\Omega$ during the interval $[0,t]$. The expression $P(\bm{Z}_{\tau};\Omega,t)$ is a statement of the end condition, and could read "a marker $\bm{Z}_{\tau}$ has intersected the domain $\Omega$ with $\tau\in[0,t]$". Obviously, there is no Monte Carlo sampling error in the signal $I$ unless evaluation of the expectation $\mathbb{E}[P(\bm{Z}_{\tau};\Omega,t)|\bm{Z}_{0}=\bm{z}]$ introduces such.

As we saw in the introductory part, an expectation value to reach an end condition can be computed in a deterministic fashion using the BMC algorithm: it was demonstrated explicitly for a process with discrete outcomes. For a process with continuous outcomes, such as the trajectory of fast ions resulting from the equations (17) and (18), the approach is almost exactly the same. The only significant difference is that, in computing the expectation value, the sum over different outcomes has to be computed as an integral over the associated probability measure. In our case this means that each realization of the fast ion trajectory $\bm{Z}_{m\Delta t}$ should be interpreted as a function of the random variable $\bm{\xi}\sim\mathcal{N}(0,\mathbf{I})$ and the expectation ultimately taken with respect to the probability density of $\bm{\xi}$.

If we now define the function

[TABLE]

where $P(\bm{Z}_{\tau};\Omega,M\Delta t)$ refers to the end condition and $m$ labels the time step index, we may yet again build upon the conditionality of the probabilities and compute a recursion formula similar to that argued in the introduction:

[TABLE]

To evaluate the expectation value, we rely on the fact that each realization of the fast-ion trajectory from point $\bm{Z}_{(m-1)\Delta t}$ to $\bm{Z}_{m\Delta t}$, computed as a composition of the maps (17) and (18), is a time-dependent random map which we denote by

[TABLE]

Taking this into account, the recursive formula (Eliminating poor statistics in Monte-Carlo simulations of fast-ion losses to plasma-facing components and detectors) can be expressed as

[TABLE]

In computing fast ion signals, one may thus use this recursive integral expression to evaluate the necessary expectation value, without introducing Monte Carlo sampling noise. As of yet, the only source of error in the formula (23) results from the integration of the stochastic differential equations to construct the fast-ion trajectory. As the map $\bm{z}\mapsto\varphi^{\Delta t,(m-1)\Delta t}(\bm{z},\bm{\xi})$ is shared with the FMC scheme, the BMC expression (23) as such is no less accurate than the process of following markers in the FMC scheme.

Although the BMC formula might not immediately ring a bell in the fast-ion community, it has been put to good use elsewhere in the past. The methodology has helped to estimate, e.g., runaway-electron generation rates in a tokamak Zhang and del-Castillo-Negrete (2017), to price options in finance Bormetti et al. (2018), and even to craft a Langevin approach to multi-scale modeling Hirvijoki (2018) and a fluid-kinetic framework for self-consistent runaway-electron modeling Hirvijoki et al. (2018). Furthermore, the BMC algorithm shares a close ideology with the adjoint methods for Fokker-Planck equations which have been used to study, e.g., current drive Karney and Fisch (1986) and runaway electrons Liu et al. (2016).

An example: – We will demonstrate the numerical implementation of the expression (23) using a 1-D stochastic harmonic oscillator as a mock-up as it displays the most essential characteristics of charged-particle motion in a plasma: the circular Hamiltonian motion is akin to the orbit topologies of charged particles with their gyro, bounce, and passing motion, and the stochastic addition to the velocity mimics the effects of Coulomb collisions, resulting in transport across the isosurfaces of energy. More importantly, this simplified model allows us to illustrate the effect of artificial diffusion which might arise if the marker evolution is dominated by Hamiltonian dynamics – this is the case for fast ions in tokamak or stellarator plasmas – and proper care is not taken while advancing the marker trajectories.

Our model stochastic differential equations read

[TABLE]

and, with $\sigma$ set to zero, the marker characteristics would draw circles centered around the point $(x,v)=(0,0)$. With $\sigma$ being non-zero, we expect deviations from the circles and, were we to draw a rectangular box centered at $(x,v)=(0,0)$, eventual intersections of the marker trajectories with the box. To construct the map for moving markers around, we use the symplectic Euler for the Hamiltonian part and the Euler-Maryama scheme for the stochastic part, providing us with the explicit maps

[TABLE]

We will then combine these maps into the following composition

[TABLE]

where the number of symplectic Euler substeps, $n$, can be chosen freely. The actual time step $\Delta t$ will be chosen so that the condition $\sigma^{2}\Delta t\ll 1$ is satisfied.

The reasoning behind multiple, optional Hamiltonian substeps is to minimize artificial numerical diffusion when $\sigma$ is small but yet finite: In computing the recursive update (23), the function $\Phi$ has to be presented in some finite form, typically in a mesh. As the map $\bm{z}\mapsto\varphi^{\Delta t,(m-1)\Delta t}(\bm{z},\bm{\xi})$ will rarely result in locations that coincide with the mesh points, interpolation of $\Phi$ will be needed. This interpolation will introduce artificial diffusion in a similar way as is encountered in semi-Lagrangian schemes for advection equations. By choosing the time step $\Delta t$ to allow for maximum effect of physical diffusion to take place, the effect of artificial diffusion can be negated. If $\sigma$ is small, a large time step can be allowed, and hence sub steps become necessary for correctly evaluating the Hamiltonian contribution.

In our demonstration simulation, we choose a rectangular box for the computational domain, with corners at $(x,v)=(-2,-2),(2,-2),(2,2),(-2,2)$. The boundaries of the box will be indicated with a purple contour. The vertical lines at $x=\pm 2$ can be seen to represent the wall surface of a tokamak or a stellarator and, to model a FILD, we consider segments of the vertical lines, with each segment naturally corresponding to a different velocity regime, mimicking the different energy channels of a real FILD. The function $\Phi(x,v)$ is represented on a fixed 2-D rectilinear mesh, with 51 mesh points in both $x$ and $v$ directions, and linear interpolation applied in between the node points. The integral in (23), in our mock-up case a 1-D Gaussian integral, is computed with a Gauss-Hermite quadrature. The characteristic time scale in the problem, the Hamiltonian orbit time, is $2\pi$ for all closed Hamiltonian orbits.

In figure 1, we see how $\Phi(x,v)$, the probability for markers to reach the red segment in given time, branches out backwards in time, reaching further out with more time given for markers to reach the red segment. One should note that the domain of the red segment would be entirely inaccessible for closed Hamiltonian orbits without diffusion. The example thus demonstrates that the BMC algorithm is able to capture probabilities of particles to reach domains that have a very narrow window of entrance, such as the FILD in a tokamak or stellarator environment. For parameters in this simulation, we chose $\sigma=0.1$, $\Delta t=0.01/\sigma^{2}=1.0$, and $n=12$. Figure 2, on the other hand, illustrates the effect of choosing the time step for the recursive iteration correctly. Too short a time step results in artificial diffusion similarly as in semi-Lagrangian schemes for advection equations and is clearly visible in the right frame of the figure, in the form of extra smearing of the function $\Phi(x,v)$.

Summary and concluding thoughts – In this paper, we have depicted a scheme to eliminate statistical errors from Monte Carlo simulations of fast-ion-loss signals. This is a problem that plagues, e.g., estimations of FILD signals in tokamak and stellarator plasmas as very few of the vast amount of marker particles used in the forward Monte Carlo simulations reach the final target domain. Our results with the mock-up, the 1-D stochastic harmonic oscillator, are encouraging in that the scheme is practical and easy to program but also reveal that further studies should perhaps be conducted to experiment with mesh generation techniques before plunging into the deep end of the pool, the fast ions in tokamak and stellarator geometry. One can, e.g., interpret from the figures that a rectilinear mesh might not offer the best fit to circular Hamiltonian orbits and sparse-grid techniques could result in savings in both memory and idling of computing nodes in a large-scale implementation. Finally, we mention that the Backward Monte Carlo scheme is likely to scale up favourably in a parallel computing environment. After all, the time stepping is explicit in nature and communication between different nodes could be optimized with proper domain decomposition mimicking the trajectories of the Hamiltonian orbits.