Backward Monte-Carlo applied to muon transport

TL;DR

This paper introduces a backward Monte-Carlo method for muon transport, implemented in the PUMAS library, demonstrating high accuracy and potential advantages for muography imaging applications.

Contribution

It presents a novel backward Monte-Carlo technique for muon transport, offering greater versatility and comparable accuracy to forward schemes in muography.

Findings

Backward Monte-Carlo agrees with forward schemes within 1% accuracy.

The PUMAS library successfully implements the backward Monte-Carlo method.

The approach enhances muography imaging capabilities.

Abstract

We discuss a backward Monte-Carlo technique for muon transport problem, with emphasis on its application in muography. Backward Monte-Carlo allows exclusive sampling of a final state by reversing the simulation flow. In practice it can be made analogous to an adjoint Monte-Carlo, though it is more versatile for muon transport. A backward Monte-Carlo was implemented as a dedicated muon transport library: PUMAS. It is shown for case studies relevant for muography imaging that the implementations of forward and backward Monte-Carlo schemes agree to better than 1%.

Tables

Table 1. Table 1: Selection probability for the toy case.

code	sampling	scheme	$p_s$ (%)	${\sigma }_{p_s}$ (%)
MUM	forward	hybrid	52.750	0.016
PUMAS	forward	hybrid	58.097	0.017
PUMAS	backward	hybrid	58.142	0.037
Geant4	forward	detailed	48.632	0.016
PUMAS	forward	detailed	49.008	0.017
PUMAS	backward	detailed	49.492	0.034

Equations

$$\mathcal{B}(ϵ)=-\frac{\partial \mathcal{E}}{\partial x}({\mathcal{E}}^{-1}(ϵ))$$

$$P_{\lambda ,\mathcal{B}}({ϵ}_f;{ϵ}_i)=\exp (-{\int }_{{ϵ}_f}^{{ϵ}_i}\frac{dϵ}{\lambda (ϵ)\mathcal{B}(ϵ)})$$

$$\begin{array}{rl}{\rho }_f({\mathbf{s}}_f)=& \int {\tau }_{i,f}({\mathbf{s}}_f;{\mathbf{s}}_i){\rho }_i({\mathbf{s}}_i)d{\mathbf{s}}_i\simeq \frac{1}{N}\sum _{k=1}^N{\omega }_{i,k}{\rho }_i({\mathbf{s}}_{i,k})\\ {\mathbf{s}}_{i,k}=& g^{-1}({\mathbf{s}}_f;x_k)\\ {\omega }_{i,k}=& \det (J_{g^{-1},{\mathbf{s}}_f}({\mathbf{s}}_f;x_k))\end{array}$$

$${ϵ}_f=P_{\lambda ,\beta }^{-1}(P_{\lambda ,\beta }({ϵ}_i)+\ln (x))$$

$${ϵ}_i=P_{\lambda ,\beta }^{-1}(P_{\lambda ,\beta }({ϵ}_f)-\ln (x))$$

$$\left|\frac{\partial {ϵ}_i}{\partial {ϵ}_f}\right|=\frac{\lambda ({ϵ}_i)\beta ({ϵ}_i)}{\lambda ({ϵ}_f)\beta ({ϵ}_f)}$$

$$\begin{array}{rl}{\rho }_n({\mathbf{s}}_n)=& \int {\tau }_{n-1,n}d{\mathbf{s}}_{n-1}\int {\tau }_{n-2,n-1}d{\mathbf{s}}_{n-2}...\int {\tau }_{0,1}{\rho }_0({\mathbf{s}}_0)d{\mathbf{s}}_0\\ \simeq & \frac{1}{N}\sum _j{\rho }_0({\mathbf{s}}_{0,j})w_{0,j}\\ s_{0,j}=& g_0^{-1}\circ ...g_{n-2}^{-1}\circ g_{n-1}^{-1}({\mathbf{s}}_n;x_{n-1,j})\\ w_{0,j}=& \prod _{k=1}^n\det (J_{g_{k-1}^{-1},{\mathbf{s}}_k})=\left|\frac{\partial {\mathbf{s}}_{0,j}}{\partial {\mathbf{s}}_n}\right|\end{array}$$

$$\begin{array}{rl}{\rho }_f({\mathbf{s}}_f)=& \int {\tau }_{i,f}({\mathbf{s}}_f;{\mathbf{s}}_i){\rho }_i({\mathbf{s}}_i)d{\mathbf{s}}_i\simeq \frac{1}{N}\sum _{k=1}^N{\rho }_i({\mathbf{s}}_{i,k}){\omega }_{i,k}\\ {\mathbf{s}}_{i,k}=& g_b^{-1}({\mathbf{s}}_f;x_k)\\ {\omega }_{i,k}=& \frac{{\tau }_{i,f}({\mathbf{s}}_f;{\mathbf{s}}_{i,k})}{{\tau }_{i,f}^{(b)}({\mathbf{s}}_f;{\mathbf{s}}_{i,k})}{J}_{g_b^{-1},{\mathbf{s}}_f}({\mathbf{s}}_f;x_k)\end{array}$$

$${\tau }_{i,f}({\mathbf{s}}_f;{\mathbf{s}}_i)=\sum _{\ell =1}^m{p}_{\ell }({\mathbf{s}}_i){\tau }_{i,f}^{(\ell )}({\mathbf{s}}_f;{\mathbf{s}}_i)$$

$$\begin{array}{rl}{\rho }_f({\mathbf{s}}_f)=& \int {\tau }_{i,f}({\mathbf{s}}_f;{\mathbf{s}}_i){\rho }_i({\mathbf{s}}_i)d{\mathbf{s}}_i\simeq \frac{1}{N}\sum _{k=1}^N{\rho }_i({\mathbf{s}}_{i,k}){\omega }_{i,k}\\ {\mathbf{s}}_{i,k}=& g_{{\ell }_0,b}^{-1}({\mathbf{s}}_f;x_{{\ell }_0,k})\\ {\omega }_{i,k}=& \frac{p_{{\ell }_0}({\mathbf{s}}_{i,k}){\tau }_{i,f}^{({\ell }_0)}({\mathbf{s}}_f;{\mathbf{s}}_{i,k})}{p_{{\ell }_0,b}({\mathbf{s}}_f){\tau }_{i,f}^{({\ell }_0,b)}({\mathbf{s}}_f;{\mathbf{s}}_{i,k})}{J}_{g_{{\ell }_0,b}^{-1},{\mathbf{s}}_f}({\mathbf{s}}_f;x_{{\ell }_0,k})\end{array}$$

$${\omega }_{i,f}=\frac{\rho ({\vec{r}}_f)}{\rho ({\vec{r}}_i)}\frac{\lambda ({ϵ}_i)}{\lambda ({ϵ}_f)}\frac{\beta ({ϵ}_i)}{\beta ({ϵ}_f)}$$

$${\rho }_S(ϵ,\vec{r},\vec{u})={\varphi }_S(ϵ,\vec{r},\vec{u})\frac{\rho (\vec{r})}{\lambda (ϵ)}$$

$${\omega }_{i,f}=\left|\frac{\partial {ϵ}_i}{\partial {ϵ}_f}\right|$$

$${\omega }_{j+1,j}^{(T)}=\frac{\beta ({ϵ}_{j+1})}{\beta ({ϵ}_j)}$$

$${\omega }_{j+1,j}^{(V)}=\frac{\lambda ({ϵ}_j)}{\lambda ({ϵ}_{j+1})}\left|\frac{\partial {ϵ}_{j+1}}{\partial {ϵ}_j}\right|{\omega }_b({\mathbf{s}}_{j+1};{\mathbf{s}}_j)$$

$${\varphi }_f({ϵ}_f,{\vec{r}}_f,{\vec{u}}_f)\simeq \frac{1}{N}\sum _{k=1}^N{w}_{i,k}{\varphi }_i({ϵ}_{i,k},{\vec{r}}_{i,k},{\vec{u}}_{i,k})$$

$${\int }_{\mathcal{D}}{\varphi }_{f}d{\mathbf{s}}_f\simeq \frac{1}{N}\sum _{k=1}^N{w}_{i,k}\frac{{\varphi }_i({\mathbf{s}}_{i,k})}{{\varphi }_b({\mathbf{s}}_{f,k})}$$

$$p_{\ell }({ϵ}_i)=\frac{\lambda ({ϵ}_i)}{{\lambda }_{\ell }({ϵ}_i)}$$

$$F_{\sigma ,h_b}^{†}({ϵ}_i;{ϵ}_f)={\int }_{{ϵ}_f}^{{ϵ}_i}\frac{\partial \sigma }{\partial {ϵ}_f}({ϵ}_f;{ϵ}_i)h_b({ϵ}_i,{ϵ}_f)d{ϵ}_i$$

$$P_{\sigma ,h_b}^{†}({ϵ}_i;{ϵ}_f)=F_{\sigma ,h_b}^{†}({ϵ}_i;{ϵ}_f)/F_{\sigma ,h_b}^{†}(+\infty ;{ϵ}_f)$$

$${\omega }_{i,f}^{(V)}=\frac{\lambda ({ϵ}_f)}{{\lambda }^{†}({ϵ}_f)}\frac{1}{h_b({ϵ}_i,{ϵ}_f)}$$

$${\omega }_{j+1,j}^{(T)}=\frac{⟨\beta ⟩({ϵ}_{j+1})}{⟨\beta ⟩({ϵ}_j)}$$

$$p_0=\frac{\ln (\frac{1}{{\nu }_{\text{CEL}}})}{\ln (\frac{1}{{\nu }_{\text{CEL}}})+\ln \left(\frac{{ϵ}_0}{{ϵ}_0-{ϵ}_f}\right)}\theta ({ϵ}_0-{ϵ}_f)$$

$${\rho }_{\mathrm{DEL}}({ϵ}_f;{ϵ}_i)\propto \frac{1}{({ϵ}_i-{ϵ}_f{)}^{{\alpha }_{\mathcal{E}}}}$$

$${\mu }_n(\mathbb{I})=\frac{1}{N}\sum _{i\in \mathbb{I}}{\omega }_i^n$$

$$\begin{array}{rl}{\rho }_Y=& {\rho }_X/\sqrt{\det (J_{g^{-1},X}^T{J}_{g^{-1},X})}\\ =& {\tau }_{i,f}\sqrt{\frac{\det (J_{g,X}^T{J}_{g,X})}{\det (g_{g^{-1},X}^T{J}_{g^{-1},X})}}\end{array}$$

$$g^{-1}(g({\mathbf{s}}_i;X);X)={\mathbf{s}}_i$$

$$\det (J_{g^{-1},{\mathbf{s}}_f}{)}^2=\frac{\det (J_{g^{-1},X}^T{J}_{g^{-1},X})}{\det (J_{g,X}^T{J}_{g,X})}$$

$$\begin{array}{rl}\frac{1}{N}\sum _{k=1}^N{\omega }_{i,k}{\rho }_i({\mathbf{s}}_{i,k})\simeq & \int {\rho }_i({\mathbf{s}}_i)\det (J_{g^{-1},{\mathbf{s}}_f}){\rho }_Y({\mathbf{s}}_i;{\mathbf{s}}_f)d{\mathbf{s}}_i\\ \simeq & \int {\tau }_{i,f}({\mathbf{s}}_f;{\mathbf{s}}_i){\rho }_i({\mathbf{s}}_i)d{\mathbf{s}}_i={\rho }_f({\mathbf{s}}_f)\end{array}$$

$$\omega_{f,i,\mathcal{E}}=\left|\frac{\partial\epsilon_{f}}{\partial\epsilon_{i}}\right|\ \mathrm{and}\ \omega_{f,i,G}=\left|\begin{array}[]{cc}\frac{\partial\vec{r}_{f}}{\partial\vec{r}_{i}}&\frac{\partial\vec{r}_{f}}{\partial\vec{u}_{i}}\\ \frac{\partial\vec{u}_{f}}{\partial\vec{r}_{i}}&\frac{\partial\vec{u}_{f}}{\partial\vec{u}_{i}}\end{array}\right|$$

Full text

††journal: Computer Physics Communications

1 Introduction

In the present section we start by giving the motivation for this work and by setting the context. Section 2 is dedicated to a general discussion of the Backward Monte-Carlo (BMC) technique as applied to the muon transport problem. Then, a specific implementation is detailed in section 3. Finally, validation results can be found in section 4.

1.1 Motivation

This work has been largely motivated by the muography imaging. Muography is a radiography technique for large scale, $0.1-2$ km, dense structures, e.g. pyramids [1, 2], small volcanoes [3, 4, 5, 6, 7], etc …. The measurement of the transmitted flux of atmospheric muons through a target allows inferring the integrated density along lines of sight. A precise density measurement requires dedicated Monte-Carlo (MC) computations. These computations can be very inefficient since in a classical MC one has to sample muons produced over the whole atmosphere but going through a tiny $\simeq 1\text{m}^{2}$ detection plane. This is especially problematic for low energy muons that can scatter away a lot from their initial directions at production. Although those non ballistic muons are not of direct interest for muography they contribute -as well as other particles- to the background [8, 9]. Therefore it is desirable to account for them by simulation, if they cannot be rejected experimentally. A successful strategy in such cases is to reverse the simulation flow using Adjoint Monte-Carlo (AMC) techniques. This guarantees to sample only useful events by starting from the detector and going backwards to the source.

So far there has not been any implementation of AMC for muography. The closest work is the one of Desorgher et al. [10] with an implementation in Geant4 [11] for space radiation problems, i.e. $e^{\pm}$, $p$ and $\gamma$. It would have been possible to extend the work of Desorgher et al. to muons. But a dedicated code that would allow specific optimisations for muography was thought to be more efficient. Following this path triggered an alternative scheme for the backward sampling of muons, which is detailed herein.

The starting idea of the BMC is to directly reverse the simulation algorithmic step by step. For example, given a random procedure to generate a final state from an initial state and some random variates, one simply reverses it to express the initial state as function of the final state. When the inverse exists this is a correct procedure up to a weight factor for the conservation of the probability density, as detailed in lemma 1 of section 2. Applying this method to the muon transport problem one gets an algorithm very close to the AMC of Desorger et al. though more versatile. Since it does not rely on an adjoint formulation we call this method backward Monte-Carlo. The term reverse Monte-Carlo also found in the literature would have been misleading to us because this technique does not solve any inverse problem. Indeed, backward sampling a bunch of scattered particles will not re-focus them on the initial source. Instead, BMC deals with the exact same problem as the forward case but it allows to generate only useful events, i.e. final states that are actually observed.

1.2 Definitions

For the present discussion a MC state is simplified as a set $\mathbf{s}=\{\epsilon,\vec{r},\vec{u}\}$ describing the components of a stationary flux of relativistic muons of kinetic energy $\epsilon$, at the space coordinate $\vec{r}$ and with momentum’s direction $\vec{u}$. The muons interactions with matter are split into a continuous component describing collective processes, e.g. multiple scattering and Continuous Energy Loss (CEL), and one or more discrete processes leading to catastrophic interactions. The following is assumed throughout this work:

Hypothesis 1** (H1)**

the work of all external forces has a null sum at every space coordinate. 2. 2.

the medium density is constant or atomic cross-sections do not depend on the medium density. 3. 3.

the differential cross-sections are invariant by rotation, e.g. unpolarised.

In practice, atomic cross-sections do actually depend on the medium density, due for example to the density effect in ionisation loss. However, in some circumstances this might be negligible, e.g. when considering a porous material with minerals of constant density and pores of air. Then, the bulk density variations actually reflect porosity variations while the energy loss is dominated by the minerals of constant density.

Note that in a realistic case the simulated geometry actually consists of several media with different compositions and densities, e.g. rocks and air in muography. However, let us first consider a single medium and then generalise to a succession of different media, provided that they all satisfy to hypothesis 1.

Under hypothesis 1, the continuous energy loss depends only on the column density along the path, $x=\int{\rho ds}$. Therefore, the CEL could be generated a priori, parametrised only by the column density, although this is not a realistic MC implementation. Thus, it is conceptually convenient to model it as following:

Definition 1** (Continuous energy loss)**

A random function $\mathcal{E}(x)$ of the column density $x$ generates a continuous energy loss from $\epsilon_{0}$ to $\epsilon_{1}$, with $\epsilon_{0}>\epsilon_{1}\geq 0$, if and only if there exists $x_{1}>0$ such that :

$\mathcal{E}(0)=\epsilon_{0}$* and $\mathcal{E}(x_{1})=\epsilon_{1}$.* 2. 2.

$\mathcal{E}$* is a strictly decreasing function over $\mathbb{I}=[0,x_{1}]$.* 3. 3.

$\mathcal{E}$* is continuous and differentiable over $\mathbb{I}$.*

Following, $\mathcal{E}^{-1}(\epsilon)$ the inverse of $\mathcal{E}$ over $\mathbb{I}$ exists. The continuous energy loss random function, $\mathcal{B}$, is defined over $[\epsilon_{1},\epsilon_{0}]$ as:

[TABLE]

where it must be understood that $\mathcal{E}$, $\partial\mathcal{E}/\partial x$ and $\mathcal{E}^{-1}$ refer to the same random variate.

Let us point out that we follow the usual convention that a random variable or a random function is written with a capital letter while a given outcome is written with the corresponding lower case letter.

The previous definition for CEL is purely mathematical. It states that CEL is indeed sustained by a continuous variation of the energy along the particle trajectory and that there is a bijective mapping between kinetic energy and column density variations. It is up to the implementation to ensure that the corresponding CEL actually observes the Physics. Given definition 1, one can associate to any MC particle a unique CEL occurrence, $\beta(\epsilon)$, fully pre-generated at the start according to $\mathcal{B}$. Again this is not a realistic MC procedure to implement but it is conceptually convenient and asymptotically exact. A realistic MC procedure would actually compute the CEL on the flight and might implement discontinuous effects during its stepping procedure as an approximation.

For discrete processes, on the contrary the kinetic energy and direction vary non continuously at some space coordinates, i.e. at vertices where interactions occur. In the following, let us call Discrete Energy Loss (DEL) these processes. Let us also write $\lambda>0$ the total interaction length for discrete processes. Under hypothesis 1, the interaction length depends only on the kinetic energy. As a result it is noteworthy that the survival probability to DEL does not depend on the geometric details of the CEL, but only on the kinetic energy, as:

[TABLE]

Note that $P_{\lambda,\mathcal{B}}$ is a random function again, derived from $\mathcal{B}$. This result is well known in the case of a deterministic CEL, e.g. see Appendix B of Sokalski et al. [12]. It is also true for a stochastic CEL provided that it satisfies definition 1.

In addition, in order to simplify the present discussion it is convenient to make the following assumption:

Hypothesis 2** (H2)**

The total macroscopic cross-section $\sigma=1/\lambda$ is strictly positive for all values of the kinetic energy in $\mathbb{R}_{+}=[0,+\infty[$.

Note that hypothesis 2 is not rigorously true in a practical implementation, as for muon transport. However, a simple way to patch this is to regularise the total cross-section by adding a do nothing sub-process, as illustrated in section 3. But, this only shifts the problem to the sub-process sampling. Under hypothesis 2, $P_{\lambda,\mathcal{B}}$ can be inverted over $\mathbb{R}_{+}$ with respect to either $\epsilon_{i}$ or $\epsilon_{f}$. Consequently the forward and backward sampling procedures are symmetrical.

2 Backward Monte-Carlo

The following is a general discussion of the BMC technique. It provides a collection of lemmas allowing to build an end to end backward Monte-Carlo simulation of high energy muons. In order not to dilute the main results all proofs have been relegated to the appendices.

2.1 Formulation

The key element of the BMC method is given by lemma 1. It states that if one can invert the MC process for generating the final state, $\mathbf{s}_{f}$, from the initial state, $\mathbf{s}_{i}$, the inverse process provides an exclusive MC estimate of the final state density, after applying a Jacobian weight factor $|\partial\mathbf{s}_{i}/\partial\mathbf{s}_{f}|$.

Lemma 1

Let $\tau_{i,f}(\mathbf{s}_{f};\mathbf{s}_{i})$ be the probability density function, with respect to $\mathbf{s}_{f}$, of a transition probability from an initial state $\mathbf{s}_{i}$ to a final state $\mathbf{s}_{f}$. Let $\mathbf{s}_{f}=g(\mathbf{s}_{i};X)$ be a Monte-Carlo process for generating $\mathbf{s}_{f}$ distributed as $\tau_{i,f}$, given $\mathbf{s}_{i}$ and a set of random variables $X$ that do not depend on $\mathbf{s}_{i}$. If $g^{-1}$ the inverse of $g$ with respect to $\mathbf{s}_{i}$ exists, then for any probability density function $\rho_{i}$ for the initial state, a backward Monte-Carlo estimate of the probability density function $\rho_{f}$ of the final state is:

[TABLE]

where $x_{k}$ is a random variate drawn from $X$ and with $J_{g^{-1},\mathbf{s}_{f}}$, informally $\partial\mathbf{s}_{i}/\partial\mathbf{s}_{f}$, the Jacobian matrix of $g^{-1}$ with respect to $\mathbf{s}_{f}$.

Note that lemma 1 is still valid when the state transition is a deterministic process, i.e. $X=\varnothing$. In this case the initial states $\mathbf{s}_{i,k}$ and their weights $\omega_{i,k}$ are all equal and one simply recovers the classical conservation law for a change of variable in a probability density function (pdf), i.e. $\rho_{f}\mathrm{d}\mathbf{s}_{f}=\rho_{i}\mathrm{d}\mathbf{s}_{i}$. The non obvious result stated here is that the same can be done event by event treating the underlying random variates, $x_{k}$, as constants. The proof is given in A.

Let us illustrate the method with a simple case useful for muon transport. Let us consider a unidimensional propagation with a deterministic CEL, $\beta$, e.g. as implemented in MUM [12]. Considering a muon with initial kinetic energy $\epsilon_{i}$, under hypothesis 2, the final kinetic energy $\epsilon_{f}$ at which the next stochastic interaction would occur can be randomised as:

[TABLE]

with $x$ a random variate drawn from a uniform distribution over $[0,1]$ and where $P_{\lambda,\beta}(\epsilon)\equiv P_{\lambda,\beta}(\epsilon;0)$ was given in equation (2). Inverting equation (4) for $\epsilon_{i}$ is straightforward, yielding:

[TABLE]

from which the Jacobian weight factor can be computed directly as:

[TABLE]

Hence, in order to sample the final density for a given kinetic energy $\epsilon_{f}$ one simply needs to draw initial values $\epsilon_{i}$ from equation (5) and to average the corresponding initial densities, $\rho_{i}$, applying the weight factor of equation (6). Note that the forward MC algorithm would instead sample $\epsilon_{i}$ from $\rho_{i}$ and then average the corresponding transition probabilities, $\tau_{i,f}$.

Given only the previous example, there would not be a significant difference between the forward and backward sampling. However, backward sampling becomes of interest when considering a succession of $n$ transition processes from state $\mathbf{s}_{0}$ to state $\mathbf{s}_{n}$. Let us write $\tau_{j,j+1}$ the density probability for a transition from $\mathbf{s}_{j}$ to $\mathbf{s}_{j+1}$ and $\mathbf{s}_{j+1}=g_{j}(\mathbf{s}_{j};X_{j})$ the corresponding MC process. Then, applying lemma 1 recursively, the BMC estimate of the final state density $\rho_{n}(\mathbf{s}_{n})$ is given as:

[TABLE]

with $\rho_{0}$ the initial density. Equation (7) states that in order to estimate the final state density for a given state $\mathbf{s}_{n}$ one applies, in reverse order, the inverse MC processes at each step as well as the corresponding Jacobian weight factors. By doing so, one gets a MC initial state that is used to sample the initial state density. In comparison, if only a given final state $\mathbf{s}_{n}$ is of interest, a forward MC estimate can be very inefficient since $\mathbf{s}_{n-1}$ might deviate significantly from the desired final state $\mathbf{s}_{n}$ after $n-1$ steps.

2.2 Biasing

In practice MC processes can be complicated to invert. Hopefully, as for the forward MC, one can substitute the true transition distribution with an approximate but more compliant one. This leads to the following biasing variant of lemma 1.

Corollary 1

Let $\tau_{i,f}^{(b)}(\mathbf{s}_{f};\mathbf{s}_{i})$ be a probability density function for $\mathbf{s}_{f}$ such that $g_{b}^{-1}$, the inverse of $\mathbf{s}_{f}=g_{b}(\mathbf{s}_{i};X)$ with respect to $\mathbf{s}_{i}$, exists. Then, a Monte-Carlo estimate of $\rho_{f}(\mathbf{s}_{f})$ is:

[TABLE]

provided that $\tau_{i,f}^{(b)}>0$ for any $\mathbf{s}_{i}$ where $\tau_{i,f}>0$.

The proof is almost identical to the previous proof of lemma 1, hence it is not detailed in this paper. Corollary 1 allows to postulate an a priori distribution for the transition probability, easier to invert, but correcting a posteriori for the true distribution value. Let us point out that the weight factor conveniently involves the Jacobian of the bias MC process, $g_{b}$, not the true one, $g$. Thus, combined with biasing, BMC becomes simple to implement.

A case of particular interest is a mixture pdf where a same final state could be reached through several components. Let us write such a pdf as:

[TABLE]

with $p_{\ell}\in[0,1]$ the weights of the mixture and $\sum_{\ell}p_{\ell}=1$. In a forward sampling one would first draw out a component, let’s say $\ell_{0}$, according to the probabilities $p_{\ell}$. Then, one generates the final state from $\mathbf{s}_{f}=g_{\ell_{0}}(\mathbf{s}_{i};X_{\ell_{0}})$. When the mixture weights depend on the initial state this procedure cannot be inverted directly. A first solution is to use an approximation to the mixture, but invertible, and rely on corollary 1 to draw the initial state. Then one can proceed as in the forward case and draw the mixture component. This can be rough though in some cases. Instead, in a backward sampling one might have a good a priori estimate of the mixture weights even though the initial state is not known. If so, corollary 2 hereafter provides an efficient bias method for the backward sampling of the mixture. Note that one might use the exact mixture components if they are easily inverted. Again the proof is similar to the one of lemma 1.

Corollary 2

Let $\tau_{i,f}(\mathbf{s}_{f};\mathbf{s}_{i})$ be a mixture pdf for $\mathbf{s}_{f}$ given by equation (9). Let us write $\tau_{i,f}^{(\ell,b)}(\mathbf{s}_{f};\mathbf{s}_{i})$ a set of $m$ pdf for $\mathbf{s}_{f}$ such that $g_{\ell,b}^{-1}$, the inverse of $\mathbf{s}_{f}=g_{\ell,b}(\mathbf{s}_{i};X_{\ell})$ with respect to $\mathbf{s}_{i}$, exists. In addition, let us write $p_{\ell,b}(\mathbf{s}_{f})$ an estimate of the probability to select the $\ell^{\text{th}}$ component in a backward sampling. Then, provided that $p_{\ell,b}\tau_{i,f}^{(\ell,b)}>0$ for any $\mathbf{s}_{i}$ where $p_{\ell}\tau^{(\ell)}_{i,f}>0$, a Monte-Carlo estimate of $\rho_{f}(\mathbf{s}_{f})$ is:

[TABLE]

where the index $\ell_{0}(k)$ is drawn according to the probabilities $p_{\ell,b}(\mathbf{s}_{f})$ and where $x_{\ell_{0},k}$ is a random variate drawn from $X_{\ell_{0}}$.

2.3 Transport

Let us consider more in detail the inversion of the MC process for muon transport. Let us consider the generic case where the initial state, $\mathbf{s}_{i}$, is at a vertex right after a DEL occured and the final state, $\mathbf{s}_{f}$, is at the next interaction vertex, before the following DEL occurs. In the forward MC procedure the muon transport between these two interaction vertices is approximated by a succession of steps over which the muon energy and direction are varied due to CEL, multiple scattering or an external magnetic field. Performing an exact inversion of this stepping process would be very inefficient, because e.g. the length of each transport step usually depends on its initial state. However the initial state is not known a priori in a backward simulation. Furthermore, the stepping part of the simulation implements a discrete approximation for a continuous transport process. Therefore, it is advantageous to define instead an approximate backward stepping process between the interaction vertices. Provided that the forward stepping is properly designed it can actually be used in backward mode with little to no modification. A possible implementation is discussed in section 3. Presently let us consider only the Jacobian weight factor corresponding to the backward transport between the two interaction vertices. This Jacobian weight is not straightforward to derive because it depends a priori on the discretisation procedure. The asymptotic result in the continuous limit is however simple and it is given hereafter, in lemma 2. A proof is provided in B. Note that in the particular case where there is no CEL, i.e. $\beta=0$, lemma 2 remains valid setting $\beta(\epsilon_{i})/\beta(\epsilon_{f})=1$ in equation (11).

Lemma 2

Let us consider a Monte-Carlo transport process satisfying to hypotheses 1 and 2 with total interaction length $\lambda>0$ and continuous energy loss $\mathcal{B}$, both in units of column density. Let $\mathbf{s}=\{\epsilon,\vec{r},\vec{u}\}$ be a Monte-Carlo state and $\rho(\vec{r})$ the variable medium density. Then, the backward Monte-Carlo weight for the transport between two interaction vertices does not depend on the details of the propagation but only on the end state properties, as:

[TABLE]

with $\epsilon_{i,f}$ and $\vec{r}_{i,f}$ the kinetic energy and corresponding vertex position after the previous interaction or before the next interaction. $\beta(\epsilon)$ is the randomised continuous energy loss between the two vertices.

2.4 Boundaries

Let us now consider the case where the MC particle hits a medium boundary or a user specified limit during its backward travel to the previous interaction vertex. Applying the backward weight given by lemma 2 would not be correct in this case. The same holds when the MC particle starts a backward propagation from a boundary or a final state, since the final position is no more an interaction vertex. A proper way to handle these boundary effects is to consider that one actually simulates a flux of particles, not a density. Thus, the initial or final state densities are defined on a surface rather than over a volume. Lemma 3, hereafter, provides a simple and generic way to convert a surface boundary to a volume boundary or reciprocally. Note that the result is still correct in the particular case of a null CEL. A proof is given in C.

Lemma 3

Let us consider a Monte-Carlo transport process with total interaction length $\lambda>0$ and continuous energy loss $\mathcal{B}$, both in units of column density. Let $\mathbf{s}=\{\epsilon,\vec{r},\vec{u}\}$ be a Monte-Carlo state and $\phi_{S}$ the state density on an interface. Then, under hypotheses 1 and 2, one can arbitrarily start or stop the particle transport on the interface using an equivalent volume state density, $\rho_{S}$, as:

[TABLE]

with $\rho(\vec{r})$ the variable medium density. The corresponding transport weight is given by lemma 2 substituting the initial or final state kinetic energy and position by the values at the interface crossing.

2.5 Vertex sampling

At an interaction vertex the MC particle’s energy and direction change discontinuously while its position remains unchanged. Let us call vertex transform the corresponding process in order to distinguish it from the transport. Inverting the vertex transform depends strongly on the nature of the process involved thus we will not give a general procedure here. Some examples are given in section 3. If there is no straightforward direct inversion it is preferable to rely on corollary 1.

However, as a general result, let us point out that the backward weight at the vertex transform depends only on the energies of the initial and final states. Therefore, there is no need to compute the Jacobian for the change of direction. A formal statement is given by lemma 4 hereafter and a proof in D.

Lemma 4

Let there be a discrete vertex transform from state $\mathbf{s}_{i}=\{\epsilon_{i},\vec{r},\vec{u}_{i}\}$ to state $\mathbf{s}_{f}=\{\epsilon_{f},\vec{r},\vec{u}_{f}\}$. Let the corresponding differential cross-section be invariant by rotation. Then, the backward weight of the vertex transform, $\omega_{i,f}$, depends only on the initial and final kinetic energies as:

[TABLE]

2.6 Flux sampling

A generic backward sampling procedure for a stationary flux of MC particles can be derived by combining the previous results. Indeed, when both initial and final conditions are specified on a surface the backward weights given by lemmas 2 and 3 nicely simplify out. Thus the final recipe is analogous to the AMC procedure.

Corollary 3

Let there be a Monte-Carlo transport process with interaction length $\lambda>0$ and continuous energy loss $\mathcal{B}$, both in units of column density. Let us assume that the vertex transforms can be backward sampled, directly or using a bias distribution. Let us write $\omega_{b}$ the corresponding bias weight if any. Under hypotheses 1 and 2, a backward Monte-Carlo estimate of the flux of a final state $\mathbf{s}_{f}=\{\epsilon_{f},\vec{r}_{f},\vec{u}_{f}\}$ is given by backward propagating $N$ particles from $\mathbf{s}_{f}$ up to an initial flux boundary $\phi_{i}(\epsilon_{i},\vec{r}_{i},\vec{u}_{i})$, with the following set of rules:

The candidate interaction vertices are drawn according to equation (5) or any equivalent procedure. 2. 2.

Any backward transport from the kinetic energy $\epsilon_{j}$ to $\epsilon_{j+1}\geq\epsilon_{j}$ generates an effective backward weight:

[TABLE]

whatever the start or end condition at $\mathbf{s}_{j}$ and $\mathbf{s}_{j+1}$ are, where $\beta(\epsilon)$ is the realised CEL on the path. 3. 3.

At any interaction vertex the effective backward weight for the interaction is:

[TABLE]

with $\mathbf{s}_{j}$ and $\mathbf{s}_{j+1}$ the states after and before the interaction, such that $\epsilon_{j+1}\geq\epsilon_{j}$. Note that if a bias procedure is used for the backward sampling of the interaction, the Jacobian factor $|\partial\epsilon_{j+1}/\partial\epsilon_{j}|$ refers to the bias process. Otherwise one must set $\omega_{b}=1$.

Let $\omega_{i,k}$ be the total backward weight given by the previous rules for the $k^{\text{th}}$ Monte-Carlo event, and $\mathbf{s}_{i,k}=\{\epsilon_{i,k},\vec{r}_{i,k},\vec{u}_{i,k}\}$ be the corresponding initial state. Then a backward estimate of the stationary flux of the final state is:

[TABLE]

Another obvious corollary concerns the BMC integration over the final states of a stationary flux. It is stated hereafter as corollary 4. Note that the integration might concern only a subset of attributes, for example only the kinetic energy thus yielding the integral flux instead of the rate.

Corollary 4

Let us consider the conditions of corollary 3 with a bias distribution, $\phi_{b}(\mathbf{s}_{f})$, for the final state over a domain $\mathcal{D}$. Then, a backward Monte-Carlo estimate of the rate of the final flux over $\mathcal{D}$ is given as:

[TABLE]

with $\omega_{i,k}$ computed according to corollary 3 and where the Monte-Carlo final states $\mathbf{s}_{f,k}$ are generated over $\mathcal{D}$ according to $\phi_{b}$.

2.7 Multiple discrete processes

So far it was convenient to consider the DEL as a single process with a total interaction length $\lambda$. Following this scheme, in the forward case whenever an interaction occurs one must first randomise the sub-process responsible for the interaction. This is done according to the occurrence probability, $p_{\ell}$, as:

[TABLE]

where $\lambda_{\ell}$ is the interaction length for the $\ell^{\text{th}}$ process, with $1/\lambda=\sum_{\ell=1}^{m}{1/\lambda_{\ell}}$, and where $\epsilon_{i}$ is the kinetic energy before the interaction. In the backward case $\epsilon_{i}$ is not known a priori. Hence one needs to rely on a bias for the sub-process selection. For example, one can rely on corollary 2 and draw the process according to $p_{\ell}(\epsilon_{f})$, with $\epsilon_{f}$ the kinetic energy after interaction. This will be a good selection procedure provided that $p(\epsilon_{i})\simeq p(\epsilon_{f})$ which is likely to be correct for muon transport since $\epsilon_{f}\simeq\epsilon_{i}$ for most DEL occurrences. A natural way to perform this biasing, in close analogy with the forward MC, is described in lemma 5 hereafter. A proof is given in E.

Lemma 5

Let there be a Monte-Carlo transport process with $m$ independent discrete energy loss sub-processes with interaction lengths $\lambda_{\ell=1,...,m}>0$ and a continuous energy loss $\mathcal{B}$, all in units of column density. Let us write $\omega_{b}^{(\ell)}$ the bias weight for the randomisation of the vertex transform for process $\ell$, if any biasing is used. Under hypothesis 1 and assuming that $\sigma_{\ell}=1/\lambda_{\ell}>0$ for all sub-processes, the following backward procedure is a correct simultaneous sampling of the previous hard interaction vertex and of its corresponding sub-process:

For all $\ell$, draw $\epsilon_{i,\ell}$ according to equation (5) but substituting $\lambda$ with $\lambda_{\ell}$. 2. 2.

Select $\epsilon_{i,\ell_{0}}=\min(\epsilon_{i,\ell})$ as the kinetic energy of the previous vertex, after interaction, and consider the corresponding sub-process, $\ell_{0}$, as the interacting one. 3. 3.

At any interaction vertex, the vertex weight, $\omega^{(V)}_{j+1,j}$, for this backward sampling procedure is given by equation (15) but substituting $\lambda$ with $\lambda_{\ell_{0}}$ and $\omega_{b}$ with $\omega_{b}^{(\ell_{0})}$.

Note that lemma 5 requires the macroscopic cross-sections of all sub-processes to be not null, $\sigma_{\ell}>0$, not only the total one as stated in hypothesis 2. This is not rigorously true for a practical implementation, as discussed in section 3 where another strategy is detailed. For example, the cross-section for a given sub-process could be set to zero below a threshold value of the kinetic energy, $\epsilon_{0}$. In this case the present bias procedure would yield a null probability to originate from the corresponding sub-process for all particles with $\epsilon_{f}<\epsilon_{0}$. This is obviously incorrect since a particle could have $\epsilon_{i}>\epsilon_{0}$ before the interaction but $\epsilon_{f}<\epsilon_{0}$ after.

2.8 Adjoint Monte-Carlo

Although derived using different methods, the BMC procedure given by corollaries 3, 4 and lemma 5 and the AMC derived by Desorgher et al. [10] are almost identical in practice. The only difference concerns the randomisation of the vertex transform and the corresponding weight, $\omega^{(V)}_{j+1,j}$. Using the adjoint process as an approximation of the inverse process is a particular biasing procedure for backward randomising the vertex transform. This particular biasing procedure leads directly to the results of Desorgher et al. [10] as shown formally in lemma 6, for which a proof is provided in F. Note however that for this method to be efficient one must postulate a relevant initial state density, as a bias factor $h_{b}(\epsilon_{f};\epsilon_{i})$, when computing the adjoint cross-section. Desorgher et al. set $h_{b}(\epsilon_{f};\epsilon_{i})=\epsilon_{f}/\epsilon_{i}$ which is also good for muon transport since the Differential Cross-Section (DCS) are soft, i.e. it is unlikely that the initial state has a much larger kinetic energy than the final one. In addition one must set an a priori high cut kinetic energy, $\epsilon_{\text{max}}$, above which no muon might exist. In this case the backward procedure is more flexible and directly well conditioned. If no simple direct inverse exists, one can bias with an analytically convenient approximation of the forward DCS, easy to invert. Thus, there is no restriction on the primary muon energy. In addition, there is no need to duplicate the pre-computations and tabulations required for the adjoint cross-sections. Note also that for efficiency, the adjoint path lengths of all sub-processes are rescaled to the forward ones in Desorgher’s et al. implementation. This is not correct when the forward cross-section of a sub-process is null, which happens in practice, for example below an energy threshold. As a result the corresponding inverse processes are arbitrarily suppressed. The AMC method however offers a conceptually straightforward inversion scheme since it is analogue to the forward case. But this comes at the price of obfuscating the flexibility of the backward scheme.

Lemma 6

Let us consider a discrete energy loss process with interaction length $\lambda$. Let us define the following integral:

[TABLE]

where $\sigma=1/\lambda$ is the macroscopic cross-section and with $h_{b}>0$ a unit less bias function depending on $\epsilon_{i}$ and $\epsilon_{f}$, e.g. $h_{b}(\epsilon_{i},\epsilon_{f})=\epsilon_{f}/\epsilon_{i}$. If $\lambda^{\dagger}=1/F_{\sigma,h_{b}}^{\dagger}(+\infty;\epsilon_{f})>0$, then:

[TABLE]

is a valid bias cumulative distribution function for the inverse vertex transform, i.e. for drawing the initial kinetic energy, $\epsilon_{i}$, given the final one, $\epsilon_{f}$. The corresponding vertex weight, defined in corollary 3 or lemma 5, is:

[TABLE]

3 Implementation

The BMC procedure discussed in section 2 has been implemented as a library initially dedicated to the transport of high energy muons used in muography. Later it was extended to taus as well. The library is named Semi-Analytical MUons Propagation (PUMAS) where the acronym reads backwards. For portability the code is written in C99 with the standard library as sole dependency. It is available from GitHub [13] under the terms of the GNU LGPLv3 licence. PUMAS was designed with a multi-threaded usage in mind, though it does not rely on a specific multi-threading implementation. The bulk of its used memory is read-only data relative to material properties and physics tabulations. This shared memory is dynamically preallocated by the user with a critical initialisation call. Afterwards, each simulation thread implements its own context and consumes a negligible amount of memory. The memory footprint of PUMAS is indeed small, $300$ kB for the library code plus $100$ to $200$ kB per loaded material. Additionally, calls to the standard library require a few extra MB. Thus, PUMAS parallelisation is very likely to be only CPU-bound on modern hardware. With a proper parallel framework, e.g. using HPCSim [14], a linear scaling of the processing speed with the number of available CPU cores can be achieved as was shown by Schweitzer [15].

In the following, we start with a brief overview of the forward MC procedure that has been implemented in PUMAS. Then, we discuss the modifications for the backward case. Note that for the sake of simplicity we limit the discussion to muons and muography, although PUMAS can also handle taus.

3.1 Energy loss

A key property for muography is the muon average energy loss, $\left<dE/dx\right>$, including both CEL and DEL. Within the Continuously Slowing Down Approximation (CSDA), the average energy loss governs the observed rate of transmitted muons through rocks in most practical cases, i.e. below $\simeq 1$ km rock depth. Tabulations of muon’s $\left<dE/dx\right>$ are provided online by the Particle Data Group (PDG) [16, 17] for a very broad range of materials. Therefore, PUMAS was designed to take these tabulations as numerical input. The PDG tables provide the $\left<dE/dx\right>$ split between the four energy loss processes relevant for muons, i.e. ionisation, pair production, bremsstrahlung and photonuclear interactions. Using these tabulated data, PUMAS allows three levels of detail for the simulation of the energy loss, configurable per simulation context as follows.

CSDA : the energy loss is fully deterministic with only CEL, i.e. $\beta=\left<dE/dx\right>$ and $\lambda=+\infty$. 2. 2.

hybrid : the energy loss is split between a deterministic CEL, $\beta<\left<dE/dx\right>$, and stochastic catastrophic DEL, $0<\lambda<+\infty$. 3. 3.

detailed : the energy loss is split as previously but the ionisation component of the CEL also has small stochastic fluctuations.

For the hybrid and detailed cases the splitting between CEL and DEL is fixed at a relative cut value on the transferred energy of $\nu_{CEL}=5\%$. This might seem high but it was shown by Sokalski et al. [12] that values as high as $10$ % can reasonably be used for high energy muons. The value set in PUMAS was cross-checked to be a good compromise between speed and accuracy for muography applications. As a result the hybrid simulation scheme à la MUM is very efficient.

For the simulation of DEL, the hybrid and detailed cases require explicit DCS for each process. The bremsstrahlung DCS has been implemented following Groom et al. [18], i.e. the same as was used for computing the PDG tables. However, for pair production Geant4’s implementation [11, 19] is used instead, since it is numerically much simpler and provides identical results. Also, for photonuclear interactions the DRSS DCS [20] is used rather than the Bezrukov and Bugaev DCS [21], used also by Groom et al.. Indeed, the latter cross-section underestimates the energy loss at PeV energies and above. Though the impact is mild for muography, it matters when considering ultra relativistic taus. Finally, for ionisation it is enough to consider only close interactions for DEL. The corresponding knock-on electron cross-section, including a high energy radiative correction, is modelled according to MUM [12].

At PUMAS initialisation the average energy loss of all DEL, $\left<dE/dx\right>_{DEL}$, is computed by numerical integration of the DCS. Then, the average value of the CEL is rescaled as $\left<\beta\right>=\left<dE/dx\right>_{PDG}-\left<dE/dx\right>_{DEL}$, where $\left<dE/dx\right>_{PDG}$ is the average energy loss read from the PDG tables. Note however that these tabulations have been initially computed only up to $10^{5}$ GeV and are not accurate above even though higher values are reported. Indeed, radiative losses are linearly extrapolated above $10^{5}$ GeV and hence underestimated. In rocks the underestimation is of the order of $\simeq 5\%$ at $10^{5}$ GeV. Of course alternative $\left<dE/dx\right>$ tables might be used as input as long as they have a format similar to the PDG ones. In particular, PUMAS includes an executable that generates energy loss tables, in the PDG format, according to its own DCS. In addition, since the initialisation of PUMAS can take a few seconds the library can also dump and reload a binary summary of all its static data, hence allowing to skip the initial pre-computations.

The detailed case has been implemented following PENELOPE [22, 23]. The CEL is fluctuated around its mean value, $\left<\beta\right>$, using a truncated Gaussian distribution for large steps, a uniform distribution for intermediate steps and a mixture pdf for short step lengths, as given in section 4.2 of PENELOPE’s 2001 manual [24]. In the continuous limit, the straggling parameter $\Omega^{2}(\epsilon)$ is given by the energy loss variance in a thick absorber, see e.g [25, 19]. In order to improve the accuracy for large steps and for consistency with the backward case, $\frac{1}{2}(\Omega^{2}(\epsilon_{i})+\Omega^{2}(\epsilon_{f}^{*}))$ is used instead of $\Omega^{2}(\epsilon_{i})$. Thus, one performs a trapeze discrete integral of the straggling along the step. In previous expressions, $\epsilon_{i}$ stands for the kinetic energy at step’s start and $\epsilon_{f}^{*}$ is the kinetic energy at step’s end in the ${\it hybrid}$ case, i.e. for mean CEL $\left<\beta\right>$. Note that the effect of this soft straggling is small : for GeV muons it leads to a $\sim 5$ % fluctuation on the muon range.

For the transport, the four inelastic DEL are treated as a single process with total interaction length $\lambda_{\text{in}}$. In the hybrid case the kinetic energy at the next interaction is randomised a priori from equation (4). Then, the corresponding column density can be derived in a unique way. In the detailed case this is not efficient. Instead the rejection method described in section 4.3.1 of PENELOPE’s manual is used at each transport step. At an interaction vertex, a sub-process is selected and the muon final kinetic energy, $\epsilon_{f}$, must be randomised according to the corresponding differential cross-section. Above $10$ GeV, considering the high cut value $\nu_{CEL}$ used, the DCS are monotone decreasing functions of the energy transfer. Thus, a Zigourat algorithm [26] is used to efficiently sample $\epsilon_{f}$ by rejection sampling from pre-computed envelopes. In addition, in order to speed up their evaluation, the DCS are initially fitted by polynomials using a log-log mapping. The order of the polynomials was set so that an agreement better than $0.3$ % is achieved over the range of use. Only in the case of ionisation the exact DCS for the knock-on electrons is computed on the flight, since its expression is already very light CPU-wise. Below $10$ GeV the sampling of $\epsilon_{f}$ is biased using a power law and the Monte-Carlo weights are corrected accordingly. Consequently, even in the forward case, PUMAS is not rigorously an analog Monte-Carlo.

3.2 Transverse transport

Transverse Transport (TT) can be enabled or disabled per simulation context. This allows useful optimisations since the impact of TT on muography is negligible for high energy muons, above $\simeq 100$ GeV, while the simulation of TT is very CPU consuming. For TT, in addition to the four energy loss processes one has to consider the Coulomb scattering off screened nuclei. It is the dominant source of TT though its contribution to the energy loss is negligible. An accurate estimation of the Coulomb differential cross-section is given by Salvat [27], for protons, using the Eikonal approximation. The application to muons is straightforward. However, the Eikonal computation is too CPU intensive to be directly used in MC. Therefore, Salvat relies on detailed tabulations. Instead, in PUMAS the Coulomb scattering has been implemented following Salvat [27] but substituting an approximate result to the Eikonal estimate, i.e. a Wentzel distribution. The atomic screening parameter is taken from Kuraev et al. [28]. It includes high energy Coulomb corrections to Moliere’s first estimate. In addition, the finite nucleus size, relevant for high energy muons, is taken into account by a nuclear form factor. Following ELSEPA [29], the nucleus is modelled by Helm’s uniform-uniform distribution. This allows analytical integration while being accurate enough for the present purpose. The final result is a fair trade off between speed, memory usage and accuracy. It yields the correct tail behaviour while being accurate to a few percent on the standard deviation, for relativistic muons.

As for the energy loss, the TT is split between a continuous multiple scattering component and catastrophic discrete events. For Coulomb interactions the split is done by cutting on the deflection angle following Fernández-Varea et al. [30]. The corresponding interaction length, $\lambda_{\text{el}}$, is tabulated at initialisation. Hence, a discrete Coulomb event is handled as a second process distinct from the DEL. The deflection angle at the vertex of a Coulomb event is randomised with the inverse cumulative distribution function (cdf) method. The implemented process has an analytical expression for the cdf. The inversion is done by a bracketing algorithm, initialised with an approximate solution. In the hybrid and detailed energy loss cases the deflection angle at the vertex of a DEL event is randomised following Geant4’s implementation. Finally, the multiple scattering process is simulated using a truncated Gaussian distribution of width $\sigma_{\theta}(\epsilon)$. The dominant contribution arises from Coulomb scattering. The soft part of photonuclear interactions is also included in $\sigma_{\theta}(\epsilon)$. Other soft processes are neglected. For accuracy and consistency with the backward case, the multiple scattering width is integrated over MC transport steps, as in the case of the energy straggling. Hence, $\frac{1}{2}(\sigma^{2}_{\theta}(\epsilon_{i})+\sigma^{2}_{\theta}(\epsilon_{f}))$ is used instead of $\sigma^{2}_{\theta}(\epsilon_{i})$ or $\sigma^{2}_{\theta}(\epsilon_{f})$.

3.3 Backward tweaks

For deterministic CEL, a convenient backward transport process is to simply use the forward process with the substitution $\vec{u}\rightarrow-\vec{u}$, i.e. the propagation direction is inverted. Since the forward transport procedure was designed accordingly, the accuracy of the corresponding backward transport is similar to the forward one. This is why it was important to integrate the multiple scattering over the transport steps. Using only the initial or final kinetic energy of the step to randomise the multiple scattering would lead to an asymmetric result and less accurate forward or backward multiple scattering. In the detailed case, for small fluctuations a good approximation for the straggling of the energy loss is to use the forward algorithm by exchanging $\epsilon_{i}$ and $\epsilon_{f}$. This is obviously wrong since one uses $\epsilon_{i}^{*}$ instead of $\epsilon_{i}$ which is not known in a backward step. But in case of large steps where the relative straggling is small, or if several transport steps are summed, the effect is negligible. The main problem is actually to correctly estimate the backward weight. When considering small fluctuations of the energy loss, a good approximation for the backward transport weight between two vertices, boundaries, etc … is given as:

[TABLE]

Provided that the stepping is long enough, the fluctuations of the energy loss at end points are negligible with respect to the bulk variation. As shown in section 4 this is accurate enough for most cases considered in muography, resulting in discrepancies on the muon flux smaller than $1$ %. Thus, in practice the same backward weight is used in the hybrid and detailed cases. But in the first case $\epsilon_{j+1}$ is deterministic, given $\epsilon_{j}$, while in the second it is stochastic.

The backward sampling of DEL or of discrete Coulomb events is straightforward. Note that in the detailed case the rejection method used in PENELOPE and PUMAS is identical to its inverse. A difficulty arises however from the fact that for numerical convenience the total inelastic cross-section, $1/\lambda_{\text{in}}$, is set to zero below a threshold value $\epsilon_{0}>0$, corresponding to the smallest non null tabulated value. In addition, even for $\epsilon>\epsilon_{0}$, for some sub-processes the cross-section is also set to zero for small kinetic energy values. For example, the parametrisations for the computation of the photonuclear cross-section become unstable at low energies. It follows that the backward sampling method discussed in section 2 would fail for inelastic processes. Consequently, some possible cases occuring during the forward sampling would be missed in the backward sampling.

Instead, the following procedure has been implemented, relying on corollaries 1 and 2. The total inelastic cross-section is regularised by adding a virtual do nothing sub-process for $\epsilon\leq\epsilon_{0}$ with a large interaction length, set equal to the last non null tabulated value. The backward sampling of DEL at vertices proceeds in three steps, as following.

First, a backward do nothing process might occur with a probability $p_{0}$, set as:

[TABLE]

If a do nothing process is drawn nothing occurs at the vertex, except that the MC weight is updated by multiplying it by $1/p_{0}$. Otherwise, an inelastic event occurs and an additional weight factor $1/(1-p_{0})$ has to be applied. 2. 2.

Then, a bias power law distribution is assumed for the discrete energy transfer at the vertex, independently of the yet unknown inelastic sub-process, as:

[TABLE]

Using corollary 1 the initial kinetic energy, $\epsilon_{i}$, is backward sampled given $\epsilon_{i}\geq\min(\epsilon_{0},\epsilon_{f}/(1-\nu_{\text{CEL}}))$. 3. 3.

Finally, the inelastic sub-process is drawn according to the differential cross-sections from $\epsilon_{i}$ to $\epsilon_{f}$, instead of using the integrated ones at $\epsilon_{i}$.

If TT is enabled, given both $\epsilon_{i}$, $\epsilon_{f}$ and the sub-process, the backward randomisation of the deflection angle is identical to the forward case. Note that (i) is a direct application of corollary 2. The probability law for $p_{0}$ is arbitrary. It was set by assuming a $1/\epsilon_{f}$ distribution for events with $\epsilon_{f}<\epsilon_{0}$ originating from inelastic processes. Indeed, particles starting with a kinetic energy well below $\epsilon_{0}$ are unlikely to originate directly from an inelastic DEL since these processes are soft.

Drawing the initial kinetic energy, $\epsilon_{i}$, before the sub-process ensures that no possibility is left aside. For relativistic muons, a value of $\alpha_{\mathcal{E}}=1.4$ was found to give good results, i.e. low sample variance for the MC. Randomising over the differential cross-sections in (iii) is not mandatory. However it ensures that the selected sub-process has both a non null interaction length for $\epsilon_{i}$ and that it indeed allows a transition from $\epsilon_{i}$ to $\epsilon_{f}$. As a result one does not backward generate useless non-physical events. The final MC weight, corrected for the biasing procedure, is guaranteed to be not null.

4 Validation

Let us now consider here some benchmark cases in order to validate the BMC implementation described in section 3. Two simulation schemes are considered. The first one, named hybrid hereafter, is a simplified 1-dimensional propagation with a deterministic CEL. It can be compared to MUM[12]. The second scheme, named detailed, implements a stochastic CEL and transverse transport as discussed previously. It can be compared to Geant4 [11], for what concerns muon transport. Comparisons with MUM and Geant4 are provided in order to show that the present implementation is indeed realistic. But the main goal is actually to cross-check the BMC and forward implementations used in PUMAS.

The comparisons have been done using the tagged version 0.8 of PUMAS. The propagation medium was set to standard rock [31] with a uniform density of $2.65\ \mathrm{g/cm^{3}}$, since it is a common material for most high energy muon transport codes. In addition the energy loss table has been self generated using PUMAS.

In order to clarify the discussion let us define formaly a few criteria used for the forecoming comparisons. Let us write $\mu_{n}(\mathbb{I})$ the nth sample moment for a subset $\mathbb{I}$ of MC events belonging to a given population, $\mathcal{P}$.

[TABLE]

where $N$ is the total sample size and $\omega_{i}$ the MC weight of the $i^{\mathrm{th}}$ event. A MC estimate of the probability to belong to $\mathcal{P}$ is given by the $1^{\mathrm{st}}$ sample moment, i.e. $\mathrm{Prob}[\mathcal{P}]\sim\mu_{1}$. Let us call systematic error the difference on $\mathrm{Prob}[\mathcal{P}]$ obtained with two different simulations, in the limit of an infinite sample size. An estimate of the statistical uncertainty on $\mathrm{Prob}[\mathcal{P}]$ is $\sigma_{\mathcal{P}}=\sqrt{(\mu_{2}-\mu_{1}^{2})/(N-1)}$. Let us call $\sigma_{\mathcal{P}}$ the MC* error* and the ratio $\sigma_{\mathcal{P}}/\mu_{1}$ the relative MC error.

4.1 Low energy case

To start with, let us consider the following toy geometry in order to cross-check the BMC implementation for muons in ionisation regime, i.e. well below the muon critical energy in standard rock, $E_{c}=693$ GeV.

The propagation medium is a box of standard rock of $10$ m thickness along the $z$-axis and infinite extension along the $x$ and $y$ directions. The initial muon flux, on the $z=0$ m face, has a kinetic energy uniformly distributed over $\epsilon_{i}\in[5,6]$ GeV and direction $u_{i,z}=\vec{u}_{i}\cdot\vec{e}_{z}\in[0.9,1]$ uniformly distributed in solid angle as well.

Then, one is interested in computing the outgoing flux on the $z=10$ m face. In particular, let us consider the distribution of the final state kinetic energy, $\epsilon_{f}$, of the muon direction, $u_{f,z}$, and of its transverse deflection, $|x_{f}-x_{i}|$. In addition let us select only outgoing muons with kinetic energy $\epsilon_{f}\in[0.1,1]$ GeV and direction $u_{f,z}\in[0.9,1]$. Let us write $p_{s}$ the corresponding probability that the outgoing muons fulfill the selection criteria.

The selection probabilities obtained with PUMAS for different settings are given in table 1. In the hybrid case the agreement between the forward and backward simulations is better than the statistical uncertainty of $0.04$ %. In the detailed case, due to the use of an approximate backward CEL process, there is a slight discrepancy of $1$ %, arising from low energy transport steps. Switching on and off the transverse transport does not change the discrepancy. Increasing the initial kinetic energy however decreases the forward/backward discrepancy since the impact of small CEL decreases. The selection probability obtained with PUMAS is close to the one obtained using Geant4, though a statistically significant ($16~{}\sigma$) deficit of $0.8$ % is visible. MUM result however should only be considered qualitatively. MUM was not designed to be used at such low energies. Below $10$ GeV the mean CEL computed by MUM starts to deviate significantly from the PDG or Geant4 tabulations.

Figures 2 to 3 show the distribution of the selected muons. For PUMAS there is an excellent agreement on the shapes obtained with forward and backward sampling. The pdfs agree within $1$ %. When comparing Geant4 results to PUMAS ones the strongest discrepancy concerns the final state kinetic energy, $\epsilon_{f}$. Geant4 yields a $3$ % larger average value of $\epsilon_{f}$ while having a slightly lower value for $p_{s}$. The missing events are mainly located in the forward region, for small deflections. These differences are likely due to the use of a higher cut value in PUMAS, $\nu_{\text{CEL}}=5$ %, for the splitting of the energy loss into discrete and continuous components. As a result, in this particular case PUMAS underestimates a few low energy catastrophic events that are wrongly redistributed into the CEL. The global impact on $p_{s}$ is however mild and decreases as the muon initial energy increases. Hence, the high value selected for $\nu_{\text{CEL}}$ is of little impact for muography applications while providing an impressive CPU gain as shown in section 4.2.

4.2 High energy case

Let us now cross-check the high energy behaviour with a toy case taken from muography, i.e. the rate of atmospheric muons surviving after propagation through a given column density of standard rock. Let the initial spectrum of atmospheric muons, $\phi_{i}$, be given by Gaisser’s parametrisation [32, 16] for an elevation of $10$ deg and an initial kinetic energy $\epsilon_{i}\geq 0.1$ GeV. Note that this model is not very accurate since the parametrisation used for $\phi_{i}$ neglects the Earth curvature and the muon decay, which is not correct below $10$ deg of elevation nor below $0.1$ GeV. However it is precise enough for the present purpose and has the advantage of being a widely used model. Let us write $\Phi(s)$ the rate of remnant muons with kinetic energy $\epsilon_{f}\geq 0.1$ GeV after a path length $s$ in the rock. Note that in the detailed case the path length is integrated along the muon trajectory. Hence, at low energies it differs from the rock thickness by a detour factor (0.98 in the worst case [18]). Considering the path length however allows for a direct comparison of the hybrid and detailed cases.

The muon rate as function of the path length was computed by MC varying the path length from $1$ mm to $10$ km with a logarithmic stepping. A log-uniform bias distribution was used for the initial or final state kinetic energy. For MUM and PUMAS the initial kinetic energy was limited to $\epsilon_{i}\leq 10^{8}$ GeV. However, Geant4 does not provide a model for photonuclear interactions at such high energies. Hence, we had to limit $\epsilon_{i}\leq 10^{6}$ GeV in this case. For the hybrid case $10^{8}$ events were simulated per value of path length, resulting in a negligible MC error. For the detailed case this number of events was limited to $10^{7}$.

The estimated rate is shown on figure 4. It spans $20$ orders of magnitude. Above $1$ km of standard rock the rate falls drastically as the range of high energy muons becomes limited by radiative processes. Hence, the MC estimate of the rate becomes very inefficient even with a log biasing. A more detailed comparison is shown on figure 5 where the PUMAS forward hybrid case is taken as a reference for the comparison.

Considering only PUMAS the forward and backward results are in perfect agreement for the hybrid case, up to statistical uncertainties of $0.1$ %. For the detailed case a $0.3$ % difference is observed with respect to the hybrid case. The MC errors are mostly negligible when compared to the observed difference, which can be considered as a systematic error. This is consistent with discretisation errors. Indeed, the hybrid scheme relies on a tabulation of the total number of interactions lengths, i.e. the integral in equation (2), while the detailed case uses a tabulation of the total cross-section instead. Both are consistent only up to discrete integration errors.

Below a few hundred meters of path length, Geant4 and PUMAS agree within $\simeq 0.5$ % systematic effects. MUM however exhibits larger deviations, $\simeq 2$ %, over the whole energy range investigated. This could be due to a simplified CEL, or to differences in the interpolation of the cross-sections for DELs. Above a few kilometers of path length Geant4 starts to deviate significantly from PUMAS and MUM. This is consistent with uncertainties on the energy loss due to photonuclear interactions. Actually, the photonuclear cross-section is the main source of systematics on the muon range at those energies, see e.g. Sokalski et al. [33]. In addition, in this comparison Geant4 suffers from a bias due to the initial kinetic energy being limited to $10^{6}$ GeV, resulting in an underestimated rate for very large path lengths.

Finally let us compare the MC integration performances. Figure 6 shows the MC relative error for the different simulation codes. For PUMAS, the detailed and hybrid cases lead to identical errors up to a factor $\sqrt{10}$ due to the different number of events. This was expected since the detailed and hybrid cases are compared for a same path length. Thus, the details of the transverse transport play no role. In the forward case, PUMAS and MUM accuracies are identical. However, Geant4 is slightly more efficient than PUMAS below $1$ km path lengths. This is because the muons propagated with Geant4 have a lower maximal value for their initial kinetic energy.

Considering the backward case, the MC error is systematically better than the forward one. The improvement is marginal though in this case, up to a factor $2$ between $10$ m and $1$ km path lengths. But since the present purpose is to validate the BMC implementation it is convenient to consider cases where the forward and backward cases are evenly efficient.

To complete this performance comparison, let us take into account the CPU cost. Let us introduce $\tau_{100}$ as the average CPU time required to reach a $1$ % MC error on the muon rate, pruned from any initialisation time. Of course the results depend on the hardware on which the code runs, so this is only indicative. For MUM and PUMAS the same PC was used, an Intel Xeon E5-$4620$ at $2.2$ GHz having $64$ cores with $128$ GB total memory. For Geant4, since the required CPU time was significantly higher the code was run on the CPU cluster of CC-IN2P3. It was however cross-checked that the average CPU performances of both systems agree within 10% for Geant4.

A summary of the achieved performance is shown on figure 7. Below $3$ km path lengths the best performance is achieved with PUMAS in backward hybrid mode. At larger depths the forward hybrid mode becomes more efficient however, e.g. by $50$ % at $10$ km. This could be due to a non optimal backward sampling of DEL or to an interplay between the MC error and the CPU cost of events. Nevertheless, path lengths larger than $3$ km are not of practical interest for muography because the transmitted flux becomes ridiculously small. Note also that PUMAS forward hybrid and MUM have similar performances, within $10$ %, when crossing more than hundred meters of rock.

Achieving such performances while being accurate enough is not straightforward. It requires a balance between float and double usage as well as efficient approximations for the differential cross-sections.

When considering the detailed simulations, PUMAS, even in the forward mode, is more efficient by at least 1 order of magnitude with respect to Geant4. In addition, for PUMAS, a significant gain is observed by using the backward mode, a factor of $10$ at large depths with respect to the forward one. In the detailed case the bottleneck is the simulation of the TT. In the forward mode a lot of CPU is wasted in the detailed TT of low energy events that would actually die without ever reaching large depths. Using the backward mode and switching between detailed mode at low energies and hybrid mode at high energies, e.g. above $100$ GeV, the CPU time can be decreased by $2$ to $3$ orders of magnitude with respect to Geant4.

5 Conclusion

BMC provides a simple way to invert the simulation flow in a muon transport MC. The elementary step processes are reverted one by one, when a trivial inverse exists, or using a bias strategy when it does not. This is a correct procedure provided that a Jacobian weight factor is applied at each step, in order to preserve unitarity. Backward sampling can provide drastic CPU gains by simulating only events that can be observed. In practice, it is similar to an adjoint MC but more flexible. A proof of principle has been provided using a dedicated library implemented for muography imaging: PUMAS. The agreement between the forward and backward sampling result is better than the $0.1$ % MC error at high energies. At low energy, due to an approximate backward procedure for the fluctuations in the ionisation loss, systematics of $1$ % can be observed. However, this systematic is washed out when considering that muography imaging uses atmospheric muons with a broad energy spectrum.

Appendix A Proof of lemma 1

Let us write $Y=g^{-1}(\mathbf{s}_{f};X)$ the random variable derived from $X$ and $g^{-1}$. The transformed variable has a density $\rho_{Y}$ given as:

[TABLE]

where $J_{g,X}$ is the Jacobian matrix with respect to $X$. In the last line of equation (26) the density $\rho_{X}$ has been substituted considering that $\mathbf{s}_{f}=g(\mathbf{s}_{i};X)$. Furthermore, since $g^{-1}$ is the inverse of $g$ with respect to $\mathbf{s}_{f}$ one has:

[TABLE]

Taking the gradient with respect to $X$ one finds after a few manipulations:

[TABLE]

Substituting back into equation (26) one gets the very simple result that $\rho_{Y}=\tau_{i,f}/\det(J_{g^{-1},\mathbf{s}_{f}})$, i.e. $Y$ is distributed as $\mathbf{s}_{f}$ up to a Jacobian factor $|\partial\mathbf{s}_{i}/\partial\mathbf{s}_{f}|$. With a Physicist’s notation this writes as a conservation law, as $\tau_{i,f}d\mathbf{s}_{f}=\tau_{f,i}d\mathbf{s}_{i}$, with $\tau_{f,i}\equiv\rho_{Y}$. Thus:

[TABLE]

where in the first line the classical MC result for weighted particles is used, with $\omega_{i,j}=\det(J_{g^{-1},\mathbf{s}_{f}})$, and in the second $\rho_{Y}$ has been substituted. This completes the proof.

Appendix B Proof of lemma 2

Let us compute the backward weight from the forward process. It goes as the inverse of the Jacobian determinant of the transform $\mathbf{s}_{i}\rightarrow\mathbf{s}_{f}$, i.e. $\omega_{i,f}=1/\omega_{f,i}=1/\det(J_{\mathbf{s}_{f},\mathbf{s}_{i}})$. Let us write $\beta(\epsilon)$ the pre-generated CEL. Then, the randomisation of the final kinetic energy, $\epsilon_{f}$, is given by equation (4). It does not depend on $\vec{r}_{i}$ nor on $\vec{u}_{i}$, therefore the forward Jacobian determinant can be factorised as $\omega_{f,i}=\omega_{i,f,\mathcal{E}}\omega_{i,f,G}$ with:

[TABLE]

The kinetic weight factor $\omega_{i,f,\mathcal{E}}$ was computed previously and is given by equation (6). The second weight factor is geometric and depends a priori on the transport from $(\vec{r}_{i},\vec{u}_{i})$ to $(\vec{r}_{f},\vec{u}_{f})$. However, it should not depend on the details of the stepping procedure, up to discretisation errors. The proof of lemma 2 requires considering the continuous limit. Therefore, let us consider the following stepping procedure with $n$ transport steps, from $s_{0}=s_{i}$ to $s_{n}=s_{f}$, where the step length is given by a uniform column density increment, as:

[TABLE]

with $x(\epsilon_{f};\epsilon_{i})$ the randomised column density between the two vertices. The stepping rule from state $\mathbf{s}_{j}$ to $\mathbf{s}_{j+1}$ can be written generically as:

[TABLE]

with $\Delta d_{j}=d_{j+1}-d_{j}$ and where $\vec{\Omega}_{j}$ generates an infinitesimal rotation, e.g. due to multiple scattering or an external magnetic field, depending on $\vec{r}_{j}$, $\vec{u}_{j}$ and $\epsilon_{j}$. The quantities $\Delta d_{j}$ and $\vec{\Omega}_{j}$ are both of order $\mathcal{O}(\frac{1}{n})$. At same order, the Jacobian for the step writes as a block matrix, as:

[TABLE]

where $I_{3}$ is the identity matrix of rank 3 and with $\otimes$, the outer product. At order $\mathcal{O}(\frac{1}{n})$ the lower right block, $D_{j}$, is invertible and its inverse writes:

[TABLE]

Following, the Schur complement of $M_{j,j+1}$ with respect to $D_{j}$ is simply given as:

[TABLE]

where we have neglected terms of order $\mathcal{O}(\frac{1}{n^{2}})$ or lower. The matrices $M_{j,j+1}$ and $D_{j}(M_{j,j+1}/D_{j})$ have the same determinant. Therefore let us write:

[TABLE]

By construction, the matrices $M_{i,f}^{(n)}$ converge to a matrix whose determinant is the geometric weight factor $J_{i,f,G}$. Furthermore, according to the matrix determinant lemma, one has:

[TABLE]

since $\nabla_{u}\cdot(\vec{u}_{j}\times\vec{\Omega}_{j})=0$. Similarly, one finds:

[TABLE]

where the last lines falls from the expression of the step length given by equation (28). Collecting results, one finally has:

[TABLE]

which completes the proof.

Appendix C Proof of lemma 3

Let us first consider the initial state boundary case with no transverse transport, i.e. 1 dimensional propagation. Then, a MC state can be summarised as $\mathbf{s}=\{\epsilon,s\}$ with $s$ the curvilinear abscissae along a straight line trajectory such that increasing values of $s$ correspond to a forward propagation. Let $s_{I}$ be the start position and $\phi_{S}(\epsilon_{I})$ be the known initial flux. Although the initial state density is not specified before $s_{I}$ one can consider an equivalent case, with an initial state density $\rho_{i}(\epsilon_{i},s_{i})$ for $s_{i}\leq s_{I}$, and yielding the same transported spectrum, $\phi_{S}$, through the interface. For the sake of simplicity let us set the arbitrary medium density before the interface to $\rho(s_{I})$. A procedure for generating the initial state density is as following. Let $\beta(\epsilon)$ be the pre-generated CEL and $\epsilon(x)$ its generating function. Then, let us first draw the kinetic energy $\epsilon_{I}$ at the interface crossing according to $\phi_{S}(\epsilon_{I})$ and read the corresponding column density, $x_{I}$, from $\epsilon(x_{I})=\epsilon_{I}$. Secondly, the initial position is drawn from an a priori arbitrary distribution $\rho_{d}(s_{i};\epsilon_{I})$. Finally, the initial kinetic energy is set according to the generated CEL, as $\epsilon_{i}=\epsilon(x_{I}-\rho_{I}(s_{I}-s_{i}))$. The corresponding initial state density goes as:

[TABLE]

By construction this generation procedure would yield the MC flux $\phi_{S}$ at the start interface if there were no interactions, only continuous energy loss. To be consistent, in a forward MC, events that interact before the start interface can be discarded and their weight set to zero. Thus, in order to yield the correct flux at the interface the density given by equation (37) needs to be divided by the survival probability up to the interface, $P_{\lambda,\beta}(\epsilon_{I};\epsilon_{i})$, given in equation (2). With this generation recipe the forward and backward procedures are identical at the interface crossing. Thus, a backward transported particle can be sampled before the initial interface according to lemma 2 and using the equivalent initial state density, $\rho_{i}$.

Let us know consider a particular case where the initial positions $s_{i}$ are distributed according to the interaction probability density, as:

[TABLE]

where one should recall that $\epsilon_{i}$ is actually a function of $s_{i}$. Following, the initial state density has a simple expression, as:

[TABLE]

where $\omega_{i,I}$ is the BMC weight factor for the transport from the initial position to the interface, given by lemma 2. Thus, in a BMC step a particle crossing the initial state interface samples the equivalent initial density, $\rho_{i}(\epsilon_{i},s_{i})$, with a weight given as:

[TABLE]

Since this weight does not depend on the initial state properties $\{\epsilon_{i},s_{i}\}$, it is useless to perform the backward transport before the interface. Everything goes as if the MC particle would have originated from the interface with an equivalent initial state density as stated in lemma 3.

With a similar reasoning one can prove the same result for the final state density. Given, a flux $\phi_{S}(\epsilon_{F})$ through an interface at $s_{F}$, one can build an equivalent final state density, $\rho_{f}(\epsilon_{f},s_{f})$, such that:

[TABLE]

Once again, everything goes as if the MC particle would have interacted on the interface with an equivalent final state density as stated in lemma 3.

The direct proof for arbitrary MC trajectories, when there is transverse transport, is not obvious because building a realistic equivalent state density would a priori depend on the 3 dimensional topology now. However, an equivalent generation procedure does not necessarily require a realistic geometry. It only needs to yield the correct flux at the interface while conforming to the Physics of the transport. In order to simplify the discussion let us only consider the case of an initial state interface. For any initial state on the interface, $\mathbf{s}_{I}=\{\epsilon_{I},\vec{r}_{I},\vec{u}_{I}\}$, let us attach an initial state space of infinite extension connected at $\vec{r}_{I}$ with a uniform density $\rho(\vec{r}_{I})$. The material is the same than the one after the interface, such that the Physics of the transport is unchanged. Since the initial state space is invariant by space translation and rotation an initial state is summarised as $\{\epsilon_{i},s_{i}\}$ where $s_{i}\leq s_{I}$ is the curvilinear abscissae along the MC trajectory, whose details do not need to be simulated. Following, an initial state is generated as previously by drawing a state $\mathbf{s}_{I}=\{\epsilon_{I},\vec{r}_{I},\vec{u}_{I}\}$ on the interface and a curvilinear distance $d$ over $\rho_{d}(s_{i};\epsilon_{I})$ from which the initial kinetic energy, $\epsilon_{i}$ is computed, according to the precomputed continuous energy loss. The particle is transported in the initial state space until it reaches $s_{I}=s_{i}+d$ where it enters the true geometry at $\vec{r}_{I}$ with direction $\vec{u}_{I}$. The corresponding backward transport is straightforward. When the backward state reaches the interface it continues propagating in the same material with a uniform density and no boundaries. With this procedure and the same reasoning as for the previous unidimensional case, one obviously recovers the results of equation (40).

Appendix D Proof of lemma 4

Let us consider the following two steps procedure for simulating the vertex transform. First the discrete energy loss, $q=\epsilon_{i}-\epsilon_{f}$ is randomised from the marginal macroscopic cross-section $\partial\sigma/\partial q$. Second, $\vec{u}_{f}$ is generated given $\vec{u}_{i}$, $\epsilon_{i}$ and $\epsilon_{f}$, e.g. according to the kinematics or from a conditional distribution. Following this procedure, the randomisation of $\epsilon_{f}$ does not depend on $\vec{u}_{i}$, since the differential cross-sections are invariant by rotation. Thus the backward weight can be factorised as:

[TABLE]

The transform $\vec{u}_{i}\rightarrow\vec{u}_{f}$ is a pure rotation, $\mathcal{R}$. Again, since the differential cross-sections are invariant by rotation it does not depend on $\vec{u}_{i}$. It depends only on $\epsilon_{i}$, $\epsilon_{f}$ and some random numbers. Therefore its Jacobian determinant is $\det(\mathcal{R})=1$, which completes the proof.

Appendix E Proof of lemma 5

Let us write $\mathbf{s}^{*}_{i}$ the initial state at the interaction vertex but before the interaction occurs, $\mathbf{s}_{i}$ the state right after and $\mathbf{s}_{f}$ the final state at the next vertex, before interaction. Let us first prove that the proposed procedure yields the correct total interaction length $1/\lambda=1/\lambda_{1}+...1/\lambda_{m}$ while selecting the interacting process with the bias probability $p_{\ell}(\epsilon_{i})=\lambda(\epsilon_{i})/\lambda_{\ell}(\epsilon_{i})$, instead of $p_{\ell}(\epsilon_{i}^{*})$. An inductive procedure is used for the proof. Let us initialise the procedure by considering the case of two independent processes. Let us write $\epsilon_{i,1}$ and $\epsilon_{i,2}$ the random variables with complementary cdf $P_{\lambda_{1},\beta}(\epsilon_{i};\epsilon_{f})$ and $P_{\lambda_{2},\beta}(\epsilon_{i};\epsilon_{f})$. The corresponding pdf write:

[TABLE]

Let us now consider the random variable $\mathcal{E}_{i}=\min(\mathcal{E}_{i,1},\mathcal{E}_{i,2})$ and the plane $(\epsilon_{i,1},\epsilon_{i,2})$ for the joint occurence of $\mathcal{E}_{i,1}$ and $\mathcal{E}_{i,2}$. The event $\{\mathcal{E}_{i}\geq\epsilon_{i}\}$ requires both $\{\mathcal{E}_{i,1}\geq\epsilon_{i}\}$ and $\{\mathcal{E}_{i,2}\geq\epsilon_{i}\}$ which are independent events. Thus:

[TABLE]

where $1/(\lambda_{1}\oplus\lambda_{2})=1/\lambda_{1}+1/\lambda_{2}$. Therefore $\mathcal{E}_{i}$ has the same cdf than a single process that would be randomised with the equivalent interaction length $\lambda=\lambda_{1}\oplus\lambda_{2}$. Furthermore, for a given occurrence $\epsilon_{i}$ of $\mathcal{E}_{i}$ the probability $p_{1}(\epsilon_{f};\epsilon_{i})$ that the $1^{\text{st}}$ process was the minimiser of $\mathcal{E}_{i}$ writes:

[TABLE]

which is indeed the expected bias probability.

Let us now consider an inductive step and assume that the randomisation procedure is correct for $m-1\geq 2$ processes. Then, let us write $\mathcal{E}_{i}^{(m-1)}=\min(\mathcal{E}_{i,1},...,\mathcal{E}_{i,m-1})$. Following from the previous assumption, $\mathcal{E}_{i}^{(m-1)}$ is distributed as a single random variable with an equivalent interaction length $\lambda^{(m-1)}=\lambda_{1}\oplus...\oplus\lambda_{m-1}$ and with occurrence probability for the $\ell^{\text{th}}$ sub-process given as $p_{\ell}^{(m-1)}(\epsilon_{i})=\lambda^{(m-1)}(\epsilon_{i})/\lambda_{\ell}(\epsilon_{i})$. In addition, since the procedure is in particular correct for $m=2$, the random variable $\mathcal{E}_{i}=\min(\mathcal{E}_{i}^{(m-1)},\mathcal{E}_{i,m})$ is also distributed as a single random variable with equivalent interaction length $\lambda^{(m)}=\lambda^{(m-1)}\oplus\lambda_{m}$ and selection probabilities $p_{\ell<m}^{(m)}=\lambda^{(m)}/\lambda^{(m-1)}p_{\ell}^{(m-1)}$ or $p_{m}^{(m)}(\epsilon_{i})=\lambda^{(m)}/\lambda^{(m-1)}$. Developing one finds $\lambda^{(m)}=\lambda$ and $p_{\ell}^{(m)}(\epsilon_{i})=\lambda(\epsilon_{i})/\lambda_{\ell}(\epsilon_{i})$, which completes the inductive demonstration.

Finally, at each interaction vertex one needs to correct the backward weight from the bias factor due to the sub-process selection, as given by corollary 2. The corresponding bias weight writes:

[TABLE]

where $\omega^{(j)}_{b}$ is any additional bias due to the vertex transform of the sub-process. This is equivalent to substituting $\lambda$ by $\lambda_{j}$ and $\omega_{b}$ by $\omega_{b}^{(j)}$ in equation (15), which completes the proof.

Appendix F Proof of lemma 6

Let us consider the forward interaction process at the vertex and let us explicitly write $\tau_{i,f}$ the pdf of the transform from initial state $\mathbf{s}_{i}$ to final state $\mathbf{s}_{f}$, as:

[TABLE]

where $\theta$ is the Heaviside step function. Let us write $\tau^{\dagger}$ the pdf of the adjoint process, given as:

[TABLE]

The corresponding cdf is given by $P_{\sigma,h_{b}}$, defined in equation (20). It provides a possible randomisation process, $g^{\dagger,-1}$, for the initial kinetic energy, biased with respect to the direct inverse, $g^{-1}$ as defined in lemma 1. With a similar reasoning than for the proof of lemma 1 one can show that the pdf of the corresponding biased forward process, $g^{\dagger}$, writes:

[TABLE]

Thus, the bias factor, $\omega_{b}$, in equation (15) goes as:

[TABLE]

Injecting this result in equation (15) yields equation (21) which completes the proof.