# MODEL FREE COLLISION AGGREGATION FOR THE COMPUTATION OF ESCAPE DISTRIBUTIONS

LAETITIA LAGUZET<fnref n="1" />* AND GABRIEL TURINICI<fnref n="2" />, *

**Abstract.** Motivated by a heat radiative transport equation, we consider a particle undergoing collisions in a space-time domain and propose a method to sample its escape time, space and direction from the domain. The first step of the procedure is an estimation of how many elementary collisions is safe to take before chances of exiting the domain are too high; then these collisions are aggregated into a single movement. The method does not use any model nor any particular regime of parameters. We give theoretical results both under the normal approximation and without it and test the method on some benchmarks from the literature. The results confirm the theoretical predictions and show that the proposal is an efficient method to sample the escape distribution of the particle.

**Mathematics Subject Classification.** 65C05, 80M31, 35Q49.

Received March 26, 2024. Accepted September 19, 2024.

## 1. INTRODUCTION

Particle simulations offering insights into complex chemical systems at the molecular level and can help elucidate reaction mechanisms, predict thermodynamic properties, and explore molecular assembly; such methods have been successfully applied in various fields such as catalysis [1], atmospheric modelling [2, 3], radiation transport [4], etc. We will focus on the integro-differential transport equation:

$$ \frac{1}{c}\partial_t u(t, x, \omega) + \omega \cdot \nabla u(t, x, \omega) + (\sigma_a(t, x) + \sigma_s(t, x))u(t, x, \omega) = \sigma_s(t, x)\langle u \rangle(t, x), \quad (1) $$

with time variable $t \in [0, T]$, position variable $x \in \mathcal{X} \subset \mathbb{R}^d$ (here $d \ge 1$ is the spatial dimension), the angle of propagation $\omega \in \mathcal{S}^d$ (unit sphere in $\mathbb{R}^d$) and $\langle u \rangle(t, x) = \frac{\int_{\mathcal{S}^d} u(t, x, \omega') d\omega'}{\int_{\mathcal{S}^d} 1 \cdot d\omega'}$ the angular average of $u$ on $\mathcal{S}^d$. The model represents heat radiative transfer equations that can be used for both photons and neutrons (we use the former in the numerical results).

The model also comes with the constant speed $c$; for instance for photons this will be the speed of light. The absorption opacity $\sigma_a$ and the scattering opacity $\sigma_s$ are (known) functions of $x$ and $t$ that describe the collision dynamics of the particles and more precisely the time to next collision, see Section 2.1 for details.

*Keywords and phrases.* Heat radiative transfer equation, Monte Carlo simulation, collision transport equations.

<footnote n="1">¹ CEA-DAM-DIF, F-91297 Arpajon, France.</footnote>

<footnote n="2">² CEREMADE, Université Paris Dauphine – PSL, 75116 Paris, France.</footnote>

<footnote>*Corresponding author: Gabriel.Turinici@dauphine.fr</footnote>
The spatial domain $\mathcal{X}$ is usually a mesh simulation cell and various approximations are invoked to compute relevant quantities and manage the transition of particles from one mesh cell to another. We will not discuss this but refer to [5, 6] and related literature. The main focus of this paper will be on how to compute the evolution of one particle from the initial time $t_0 = 0$ and initial position $x_0 \in \mathcal{X}$ to the moment when it exits the domain $\mathcal{X} \times [0, T]$, i.e., either reaches the spatial boundary $\partial\mathcal{X}$ or consumes all available time $T$.

We are concerned here with “Monte Carlo” approaches that regard (1) as the time-evolving probability density of a stochastic process, see [7]. When the parameter $\sigma_s$ is large, many collisions occur before final time $T$ because the average time to next collision is $\frac{1}{\sigma_s}$. This is the so-called “diffusion regime” [8] and approximation methods exist to exploit this remark, in particular the Random Walk (RW) methods [4, 9]. On the contrary, when $\sigma_s$ is small, the particle will not undergo many collisions before exiting the mesh cell and the ballistic regime is important. In between, there are situations where the diffusion limit is not valid but the number of collisions is still important and requires extensive numerical simulations. Like in the diffusion limit, one would like to somehow accelerate this computation by replacing a large sum of independent collisions with some aggregate step, without resorting to diffusive approximations. So we focus in this paper on a model free method to aggregate many collisions into a single displacement without affecting the escape distribution of the particle. The simulation of the particle’s trajectory stops when either time ends (at $t = T$) or the particle reaches the spatial boundary $\partial\mathcal{X}$. Note that a single aggregated step will probably not suffice to end the simulation for the particle so several such movements will probably be used.

As a technical circumstance, we will invoke the discrete ordinate method, denoted $S_N$, which consists in discretizing the angular direction variable i.e., replacing $S^d$ by a set of discrete directions $\mathcal{D}$.

The outline of the paper is the following: we describe the collision dynamics of one particle in Section 2.1. Then in Section 2.2 we describe the general principle of the proposed method; in Sections 2.3 and 2.4 we give several theoretical results. The procedure is then tested in Section 3.

## 2. PROCEDURE AND ASSOCIATED THEORETICAL INSIGHTS

### 2.1. Description of the collision dynamics

We present briefly the dynamical setting. For graphical convenience this problem is presented in 1D but it transcribes without difficulty to the multi-dimensional case.

For further simplification we will restrict to a situation where the direction of the particle is either +1 or -1; this is the so-called $S_2$ discrete-ordinates method (see [10], Sect. 16.3, p. 502, which relates to the Schuster-Schwarzschild equations ([10], Sect. 14.3, p. 456)). We will denote $\mathcal{D}_2 = \{-1, 1\}$ and in general

$$ \mathcal{D}_N = \left\{ \cos\left(\frac{2\pi j}{N}\right), j = 0, \dots, N - 1 \right\}. \qquad (2) $$

Note that, for dimensions $d$ higher than 1, the direction of the particle is an element of the unit sphere $S^d$. We consider a 1D particle in the spatial domain $\mathcal{X} := [0, L]$ starting from $x(0) = x_0 \in \mathcal{X}$ at $t_0 = 0$. We are also given a maximum time $T$. The particles evolve as follows: the collision counter is set to $s = 1$ and a direction $a_s \in \mathcal{D}$ is chosen at random uniformly from $\mathcal{D}$. Here $\mathcal{D} = \mathcal{D}_N$ and unless stated otherwise we set $N = 2$ but our most general results in Section 2.5 consider arbitrary values of $N$.

A time $\tau_s$ is sampled from an exponential law of mean $1/\sigma$ (this will be denoted $\tau_s \sim \text{Exp}(\sigma)$). We define $t_s = t_{s-1} + \tau_s$. The particle moves on a straight line in the direction $a_s$ during the time $\tau_s$ at constant speed $v$. Thus for any $t \in [t_{s-1}, t_s]$ : $x(t) = x(t_{s-1}) + v(t - t_{s-1})a_s$.

Next, the collision counter $s$ is incremented $s \to s + 1$ and the process repeats until either $t_{s*} \ge T$ or $x_{t_{s*}} \notin \mathcal{X}$. The precise space-time coordinates when the particle touches the first time the boundary of the domain $\mathcal{X} \times [0, T]$ are computed, i.e., for our simple 1D case the smallest $t^* \in [t_{s*-1}, t_{s*}]$ such that $x(t^*) = 0$ or $x(t^*) = L$ or $t^* = T$.
<figure><img src="image_1.png" /><figcaption>FIGURE 1. An illustration of the time dynamics and escape space-time distribution for a particle starting at $x_0$, $D = D_2$. The abscissas represent the position of the particle and the ordinates the time. For each remarkable point we give its $x, t$ coordinates, for instance $A$ has $x = 0$ and $t = 0$ and represents the left extremity of the segment at $t = 0$ while $C$ is the left extremity at the final time. The time to next collision is an exponential random variable of average $1/\sigma$ i.e., distributed $\text{Exp}(\sigma)$. The joint escape distribution $\mathcal{E}(\sigma, \mathcal{X}, T, x_0)$ collects the location $x^*$, time $t^*$ and direction $a^*$ when the boundary $x = 0, x = L$ or $t = T$ is reached for the first time (i.e., particle hits $AB, BC$ or $CD$). The simulation is stopped when this happens. The colored histograms are artist views of each conditional escape distributions: the magenta/orange is the histogram of the escape time $t^*$ given that particle escaped through the left/right while the green is the histogram of escape position $x^*$ for particles that did not escape before time $T$ or, equivalently, the escape is due to time $T$ being totally consumed.</figcaption></figure>

The quantity of interest is the distribution of the escape space-time, more precisely the joint distribution of the escape position $x^* := x(t^*)$, escape direction $a^* \in D$ and escape time $t^*$ when the boundary is reached. This triplet $(x^*, t^*, a^*)$ is a random variable whose distribution $\mathcal{E}(\sigma, \mathcal{X}, T, x_0)$ depends only on $\sigma, \mathcal{X}, T$ and $x_0$.

An illustration is given in Figure 1 for $D = D_2$. Note that the escape space-times values are elements of $\mathbb{R}^{2\text{dim}(\mathcal{X})+1} = \mathbb{R}^3$ thus the support of the distribution $\mathcal{E}(\sigma, \mathcal{X}, T, x_0)$ is included in $\mathbb{R}^3$; but there are many restrictions on this support, for instance $t \in [0, T]$ and so on, leaving the support to be included in (see Fig. 1 for notations):

$$ (AB \times \{-1\}) \cup (BC \times \{-1, 1\}) \cup (CD \times \{1\}). \qquad (3) $$

For instance, the explanation of the first element $AB \times \{-1\}$ is that if the particle exists through its left border the exit direction will point to the left. Or in $D$ there is only one direction that points to the left which is $-1$; same for $CD \times \{1\}$.

The colored areas in Figure 1 are a general illustration, not corresponding to any specific parameters, of the following three conditional distributions:

- the left area (magenta in color figure) is the distribution of the escape time $t^*$ at which $x(t^*) = 0$, conditioned by the fact that the particle touched $x = 0$ before $T$ and before touching $x = L$;

- the top area (green in color figure) is the distribution of the escape position $x(t^*)$ conditioned by escaping because time $T$ was reached before reaching $x = 0$ or $x = L$, i.e., $t^* = T$;

- the right area (orange in color figure) is the equivalent to the magenta area when border is first reached for $x = L$.
## 2.2. Idea and first estimations

The dynamics described in Section 2.1 is used to sample from the distribution $\mathcal{E}$. In particular the trajectory before the exit time i.e., $x(t)$ for $t < t^*$, is not useful and not used. We can thus imagine a way to accelerate the computation by “skipping” these intermediary steps. For instance, when the diffusion parameter $\sigma$ is very large many collisions will occur before particle exits the space-time domain and in this case a random walk approximation could be valid. We do not want to use this kind of approximation here but remain as close as possible to the collisional dynamics.

Recall that the collision time $t_n$ is the sum of i.i.d Exp($\sigma$) random variables $\tau_j: t_n = \sum_{j=1}^n \tau_j$; the law of $t_n$ is a Gamma distribution of parameters $n$ and $\sigma$.

The position $X_n = x(t_n)$ is such that $X_n = x_0 + \sum_{j=1}^n v a_j \tau_j$, its law is a mixture of sums of two Gamma distributed random variables (one for each value in $\mathcal{D}$); we will make this precise latter. In any case we will show that for any $n$ one can sample directly and exactly from the joint law $(X_n, t_n)$. In this way we can advance the time by $t_n$ units and replace $n$ individual collisions by only one sample from this joint law.

But the question that arises is the following: what is the value of $n$ such that, with high certainty, we can make $n$ steps without exiting the space-time domain $\mathcal{X} \times [0, T]$? Given some tolerance $\epsilon > 0$ the overall goal of the estimations in following sections is to find how large $n$ can be while still ensuring that:

$$ \mathbb{P}[\{\exists t \le t_n : x(t) \notin \mathcal{X}\} \cup \{t_n \ge T\}] \le \epsilon. \quad (4) $$

**Remark 1.** The event $\{\exists t \le t_n : x(t) \notin \mathcal{X}\} \cup \{t_n \ge T\}$ whose probability is evaluated above concerns trajectories where either $t_n$ is larger than $T$ or trajectory already exited the domain $\mathcal{X}$ before $t_n$. The complementary of this event are trajectories $(x(t), t)$ that remained in $\mathcal{X} \times [0, T]$ for all $t \le t_n$. For trajectories in the complementary of $\{\exists t \le t_n : x(t) \notin \mathcal{X}\} \cup \{t_n \ge T\}$ we do not need to sample $x(t_1), \dots, x(t_{n-1})$ because, with high probability, no exit occurred in $[0, t_n]$. If $\epsilon$ is small then with probability at least $1 - \epsilon$ one can thus only sample $t_n, x(t_n)$ (and the direction at time $t_n$). This is the goal of this estimation, to know how many intermediary steps $x(t_1), \dots, x(t_{n-1})$ can be safely omitted. Note that $x(t_n)$ is still sampled exactly; the acceleration provided by our procedure comes from the fact that we do not need to sample $x(t_1), \dots, x(t_{n-1})$. Then the procedure computes, starting from $x(t_n)$ another estimation for how many steps can be aggregated from this new position with very low risk of exiting the domain $\mathcal{X} \times [0, T - t_n]$ and so on.

Note that when $\epsilon$ is small this does not imply $n$ is also small, cf. for instance estimation (5) below where $\epsilon = 10^{-9}$ and $n$ can be as large as $T\sigma + 18 - 6\sqrt{T\sigma} + 9$ (which is large provided $T\sigma$ is large, see Tab. 2 for examples).

## 2.3. Estimates for $t_n$

### 2.3.1. Normal approximation for $t_n$

The first approach is to use a normal (Gaussian) approximation. This idea is close to the Random Walk regime (cf. Introduction). For instance we know that $t_n$ has mean $n/\sigma$ and variance $n/\sigma^2$. For very large values of $n$, $t_n$ will behave as a normal variable with same mean and variance, i.e., $\frac{\sigma}{\sqrt{n}}(t_n - n/\sigma)$ is “close” to a standard normal. So, we will write $\mathbb{P}[t_n \ge T] \le \epsilon$ is the same as $\mathbb{P}[\frac{\sigma}{\sqrt{n}}(t_n - n/\sigma) \ge \frac{\sigma}{\sqrt{n}}(T - n/\sigma)] \le \epsilon$ which, if the normal approximation holds, will be true when $\frac{\sigma}{\sqrt{n}}(T - n/\sigma)$ is larger than the $1 - \epsilon$ quantile of the normal distribution. To simplify things we take as small error $\epsilon = 10^{-9}$ corresponding to slightly more than 6 standard deviations (that we conservatively set to 6). In this case, with high probability, $t_n$ will not be larger than $T$ if $\frac{\sigma}{\sqrt{n}}(T - n/\sigma) > 6$ or, equivalently, $n + 6\sqrt{n} \le \sigma T$. So the aggregation rule becomes :

$$ n \le T\sigma + 18 - 6\sqrt{T\sigma} + 9. \quad (5) $$

Note that here $T\sigma$ is the mean number of collisions required to reach the final time $T$ (each collision “consumes” in average $1/\sigma$ time units).
### 2.3.2. Exact tail approximation for $t_n$

We can also give a more precise, but slightly less convenient, estimation for $n$ coming from the tail estimates for the gamma distribution.

**Proposition 2.** With previous notations, choosing $n \le T\sigma/2$:

$$ \mathbb{P}[t_n \ge T] \le \frac{z^{z/2}e^{-z}}{\Gamma(1+z/2)}, \quad \text{where } z = T\sigma. \qquad (6) $$

In particular:

$$ \mathbb{P}[t_n \ge T] \le 10^{-9}, \quad \text{if } n \le T\sigma/2, n \ge 50. \qquad (7) $$

*Proof.* We are interested in $\mathbb{P}[t_n \ge T] = \mathbb{P}[\sigma t_n \ge T\sigma]$. Note that $\sigma t_n$ is a gamma random variable with parameters shape = $n$ and scale = 1 and density $f(y) = \frac{1}{(n-1)!}y^{n-1}e^{-y}$. We need to give a bound for $\mathbb{P}[\sigma t_n \ge T\sigma] = \int_{T\sigma}^{\infty} \frac{1}{(n-1)!}y^{n-1}e^{-y} dy$. By integration by parts for general $z > 0$:

$$ \int_z^\infty \frac{y^{n-1}}{(n-1)!} e^{-y} dy = \left. \frac{y^n}{n!} e^{-y} \right|_{y=z}^\infty + \int_z^\infty \frac{y^n}{n!} e^{-y} dy \ge \frac{z}{n} \int_z^\infty \frac{y^{n-1}}{(n-1)!} e^{-y} dy - \frac{z^n}{n!} e^{-z}. \qquad (8) $$

We used the inequality $\int_z^\infty \frac{y^n}{n!} e^{-y} dy \ge \frac{z}{n} \int_z^\infty \frac{y^{n-1}}{(n-1)!} e^{-y} dy$ that is true because $y \ge z$ on the domain of integration. The term $\int_z^\infty \frac{y^{n-1}}{(n-1)!} e^{-y} dy$ appears now in both sides and thus one can write $(\frac{z}{n} - 1) \int_z^\infty \frac{y^{n-1}}{(n-1)!} e^{-y} dy \le \frac{z^n}{n!} e^{-z}$; dividing by $(\frac{z}{n} - 1)$ it follows that, for $z > n$:

$$ \int_z^\infty \frac{y^{n-1}}{(n-1)!} e^{-y} dy \le \frac{n z^n e^{-z}}{n! (z-n)} = \frac{n z^n e^{-z}}{\Gamma(n+1)(z-n)}. \qquad (9) $$

Let us take now $n = z/2$; the right hand side equals $\frac{z^{z/2}e^{-z}}{\Gamma(1+z/2)}$. The function $z \mapsto \ln(\frac{z^{z/2}e^{-z}}{\Gamma(1+z/2)}) = \frac{z}{2}\ln(z) - z - \ln(\Gamma(1+z/2))$ has derivative $\frac{\ln(z)-1-\psi(1+z/2)}{2} = \frac{\ln(z/2)-\psi(1+z/2)+\ln(2/e)}{2} \le 0$ we have used the notation for the digamma function $\psi(z) = \frac{\Gamma'(z)}{\Gamma(z)}$, inequality $\ln(u) \le \psi(u+1)$ for all $u > 0$ and $\ln(2/e) < 0$. This shows that the error term $\frac{z^{z/2}e^{-z}}{\Gamma(1+z/2)}$ is decreasing for $z \ge 2$. In particular its value at $z = 100$ is $2.44 \times 10^{-10} < 10^{-9}$. Conclusion (7) follows because $n \mapsto \mathbb{P}[t_n \ge T]$ is obviously increasing with $n$ and the value at $n = T\sigma/2$ is less than $10^{-9}$. $\square$

### 2.3.3. Summary for $t_n$

To summarize, in order to satisfy $\mathbb{P}[t_n \ge T] < 10^{-9}$ we have two possible choices

- the rigorous, conservative estimate from (7) with choice $n \le T\sigma/2$ as soon as $n \ge 50$;

- the normal approximation (5) resulting in the bound: $n \le T\sigma + 18 - 6\sqrt{T\sigma} + 9$. In practice we require $n$ to be larger than some $n_N$ that we set to $n_N := 300$ in order to ensure that the normal approximation is in the asymptotic regime.

The normal approximation provides larger (thus less restrictive) values for $n$ but its quality is not precisely quantified. On the other hand, the conservative estimate (7) has a known error bound and works even if only one hundred average collisions are left in the time interval. If no other parameters enter into the decision, the number of such “aggregated collisions” required to reach final time $T$ is $O(1)$ for the normal approximation and $O(\log(T\sigma))$ for the conservative estimate (each step halves the time “left”). Both give very good results for the numerical regimes we are interested in.
## 2.4. Estimates for the spatial boundary

We now inquire about estimating the number of collisions that can be aggregated without reaching the spatial boundary. Define $D$ to be the distance from $x_0$ to the boundary of $\mathcal{X}$ i.e., $D = \min\{x_0, L - x_0\}$. In general, if the spatial domain is $\mathcal{X}$ we set $D = \text{dist}(x_0, \partial\mathcal{X})$. Note that $\mathbb{E}[X_n - x_0] = 0$, $\mathbb{V}[X_n - x_0] = nv^2/\sigma^2$. But, this still does not tell us if some other $X_k$ for $k \le n$ did not already exited through the spatial domain boundaries $\{0\}$ or $\{L\}$. We can write:

$$ \mathbb{P}[\{\exists t \le t_n : x(t) \notin \mathcal{X}\}] \le \mathbb{P}\left[\left\{\max_{s \le n} |X_s - x_0| \ge D\right\}\right] = \mathbb{P}\left[\left\{\max_{s \le n} \left|\sum_{j=1}^s a_j \tau_j\right| \ge D/v\right\}\right]. \quad (10) $$

## 2.5. Spatial boundary treatment, no Gaussian approximation, general N

We will consider a general case when $\mathcal{D}$ is not necessary $\{-1, 1\}$. We will assume:

$$ \text{Hypothesis on directions (HD) : } \forall j \ge 0 : \mathbb{E}[a_j] = 0, \mathbb{V}[a_j] = \varsigma^2, \text{ the law of } a_j \text{ is symmetric.} \quad (11) $$

In practice this can restrict for instance the number of directions to be even for the $S_N$ model (to ensure symmetry). Note that for $\mathcal{D} = \mathcal{D}_N: \varsigma = \sqrt{N/2}$. We will need the following result.

**Proposition 3.** Consider $Y_n = \sum_{j=1}^n Z_j$ with $Z_j$ symmetric i.i.d. variables such that for some $\lambda > 0$: $\mathbb{E}[e^{\lambda Z_1}] < \infty$. Then for any $C \ge 0$

$$ \mathbb{P}\left[\max_{j \le n} |Y_j| \ge C\right] \le 4\mathbb{P}[Y_n \ge C] = 2\mathbb{P}(|Y_n| \ge C] \quad (12) $$

$$ \mathbb{P}\left[\max_{j \le n} |Y_j| \ge C\right] \le 2e^{-\lambda C} (\mathbb{E}[e^{\lambda Z_1}])^n \quad (13) $$

In particular $\mathbb{P}[\max_{j \le n} |Y_j| \ge C] \le \epsilon$ as soon as:

$$ n \le \frac{\lambda C + \ln(\epsilon) - \ln(2)}{\ln(\mathbb{E}[e^{\lambda Z_1}])}. \quad (14) $$

**Remark 4.** The estimation (12) reminds of the Kolmogorov's inequality that would read:

$$ \mathbb{P}\left[\max_{j \le n} |Y_j| \ge C\right] \le \frac{1}{C^2} \mathbb{V}[Y_n] = \frac{n}{C^2} \mathbb{V}[Z_1]. \quad (15) $$

Such an inequality is not useful because it would constraint $n$ to not be larger than $\epsilon C^2/\mathbb{V}[Z_1]$ which is very disappointing when $\epsilon \to 0$. On the other hand, Doob's inequality will be invoked in the proof of the upper bound (13) and the estimation (14) where $\epsilon$ only appears through its logarithm.

**Proof.** **Proof of inequality (12):** the general idea is that $Y_n$ can be thought close, by Donsker's reflection principle, to a Brownian motion; for a Brownian motion the reflection principle related the maximum deviation before time $t_n$ with the value at time $t_n$, i.e., the estimation (12) is true with equality. So we follow the proof of the Brownian motion reflection principle. Note first that

$$ \begin{aligned} \mathbb{P}\left[\max_{j \le n} |Y_j| \ge C\right] &= \mathbb{P}\left[\left\{\max_{j \le n} Y_j \ge C\right\} \cup \left\{\min_{j \le n} Y_j \le -C\right\}\right] \le \mathbb{P}\left[\left\{\max_{j \le n} Y_j \ge C\right\}\right] + \mathbb{P}\left[\left\{\min_{j \le n} Y_j \le -C\right\}\right] \\ &= \mathbb{P}\left[\left\{\max_{j \le n} Y_j \ge C\right\}\right] + \mathbb{P}\left[\left\{C \le -\min_{j \le n} Y_j\right\}\right] = \mathbb{P}\left[\left\{\max_{j \le n} Y_j \ge C\right\}\right] + \mathbb{P}\left[\left\{C \le \max_{j \le n} Y_j\right\}\right] \end{aligned} $$
$$ = 2\mathbb{P}\left[\left\{\max_{j \le n} Y_j \ge C\right\}\right], \qquad (16) $$

where we used the symmetry of $Y_n$. By symmetry we also obtain

$$ \mathbb{P}[|Y_n| \ge C] = 2\mathbb{P}[Y_n \ge C], \qquad (17) $$

so it is enough to show that

$$ \mathbb{P}\left[\max_{j \le n} Y_j \ge C\right] \le 2\mathbb{P}[Y_n \ge C]. \qquad (18) $$

Denote $u$ the stopping time to be the first $s$ such that $Y_s \ge C$ or $n$ if no such $s$ exists; define $\mathcal{F}_u$ to be its associated sigma-algebra. Then:

$$ \begin{aligned} \mathbb{P}\left[\max_{j \le n} Y_j \ge C\right] &= \mathbb{P}\left[\max_{j \le n} Y_j \ge C, Y_n \ge C\right] + \mathbb{P}\left[\max_{j \le n} Y_j \ge C, Y_n < C\right] \le \mathbb{P}[Y_n \ge C] \\ &\quad + \mathbb{P}\left[\max_{j \le n} Y_j \ge C, Y_n - Y_u < 0\right] \\ &= \mathbb{P}[Y_n \ge C] + \mathbb{E}\left[\mathbb{1}_{\max_{j \le n} Y_j \ge C, Y_n - Y_u < 0}\right] = \mathbb{P}[Y_n \ge C] + \mathbb{E}\left[\mathbb{1}_{\max_{j \le n} Y_j \ge C, Y_n - Y_u < 0} | \mathcal{F}_u\right] \\ &= \mathbb{P}[Y_n \ge C] + \mathbb{E}\left[\mathbb{1}_{\max_{j \le n} Y_j \ge C} \mathbb{E}\left[\mathbb{1}_{Y_n - Y_u < 0}\right]\right] \\ &= \mathbb{P}[Y_n \ge C] + \mathbb{P}\left[\max_{j \le n} Y_j \ge C\right] \cdot \mathbb{P}[Y_n - Y_u < 0] \le \mathbb{P}[Y_n \ge C] + \frac{1}{2}\mathbb{P}\left[\max_{j \le n} Y_j \ge C\right], \end{aligned} \qquad (19) $$

where we used the fact that $Y_n - Y_u$ is symmetric thus $\mathbb{P}[Y_n - Y_u < 0] \le 1/2$. The relation (18) follows by subtracting $\frac{1}{2}\mathbb{P}[\max_{j \le n} Y_j \ge C]$ from the first and last terms of the formula above.

**Proof of inequality (13):** we already proved $\mathbb{P}[\max_{j \le n} |Y_j| \ge C] \le 2\mathbb{P}[\{\max_{j \le n} Y_j \ge C\}] = 2\mathbb{P}[\{\max_{j \le n} e^{\lambda Y_j} \ge e^{\lambda C}\}]$. Since for $\lambda > 0$, the function $x \mapsto \exp(\lambda x)$ is convex, $e^{\lambda Y_n}$ is a sub-martingale. By Doob's inequality,

$$ \mathbb{P}\left[\left\{\max_{j \le n} e^{\lambda Y_j} \ge e^{\lambda C}\right\}\right] \le e^{-\lambda C} \mathbb{E}[e^{\lambda Y_n}] = e^{-\lambda C} \prod_{j=1}^{n} \mathbb{E}[e^{\lambda Z_j}] = e^{-\lambda C} (\mathbb{E}[e^{\lambda Z_1}])^n. \qquad (20) $$

**Proof of inequality (14):** it results from (13) by taking the logarithm in $2e^{-\lambda C} (\mathbb{E}[e^{\lambda Z_1}])^n \le \epsilon$. $\square$

**Corollary 5.** Let $\epsilon > 0$. Then, with the previous notations, when $\mathcal{D} = \mathcal{D}_N$ for some even value of $N \ge 2$:

$$ \mathbb{P}[\{\exists t \le t_n : x(t) \notin \mathcal{X}\}] \le \epsilon \text{ if } n \le \frac{\frac{\sigma D}{v} + \sqrt{2} \ln(\epsilon/2)}{\sqrt{2} \ln(2)}. \qquad (21) $$

*Proof.* We use (10) and inequality (14) from Proposition 3. Here $Z_j = a_j \tau_j$, $C = D/v$. We write for $\lambda < \sigma$ and $\mathcal{D} = \mathcal{D}_2$:

$$ \begin{aligned} \mathbb{E}[e^{\lambda Z_1}] &= \mathbb{E}[e^{\lambda a_1 \tau_1}] = \frac{1}{2}\mathbb{E}[e^{\lambda \tau_1} + e^{-\lambda \tau_1}] \\ &= \frac{\sigma}{2} \int_0^\infty e^{-\sigma\tau} (e^{\lambda\tau} + e^{-\lambda\tau}) d\tau = \frac{\sigma}{2} \left( \frac{1}{\sigma-\lambda} + \frac{1}{\sigma+\lambda} \right) = \frac{\sigma^2}{\sigma^2 - \lambda^2}. \end{aligned} \qquad (22) $$

Take now $\lambda = \sigma/\sqrt{2}$. The term in (14) becomes $\frac{\sigma D/(v\sqrt{2})+\ln(\epsilon/2)}{\ln(2)}$.
For $\mathcal{D} = \mathcal{D}_N$ with $N > 2$ (even) denote $d_j = \cos(\frac{2\pi j}{N})$. Then

$$ \mathbb{E}[e^{\lambda Z_1}] = \mathbb{E}[e^{\lambda a_1 \tau_1}] = \frac{1}{N} \sum_{j=0}^{N-1} \mathbb{E}[e^{\lambda d_j \tau_1}] = \frac{\sigma}{N} \sum_{j=0}^{N-1} \int_0^{\infty} e^{-\sigma \tau} e^{\lambda d_j \tau} d\tau = \frac{1}{N} \sum_{j=0}^{N-1} \frac{\sigma}{\sigma - d_j \lambda}. \quad (23) $$

In the sum over $j$, for even values of $N$, values can be regrouped by associating some value $d$ with the value $-d$. There are $N/2$ such couples and each of them will contribute $\frac{\sigma}{\sigma-d\lambda} + \frac{\sigma}{\sigma+d\lambda} = \frac{\sigma^2}{\sigma^2-d^2\lambda^2} \le \frac{\sigma^2}{\sigma^2-\lambda^2}$. We conclude as in the case $N = 2$. $\square$

This corollary can be used to know how many elementary collisions can be aggregated while keeping the chances to reach the frontier very small. In practice we will improve this bound but it has the merit to show that $\sigma D/v$ should be large enough with respect to $\ln(\epsilon)$. Conservative values for $\epsilon$ will be of order $10^{-9}$ which is a chance in a billion to be wrong and assume that particle still stays inside the domain when in reality it has exited. This would give $\sqrt{2}\ln(\epsilon/2) = -30.28$ which leads, using $\sqrt{2}\ln(\epsilon/2) \sim 0.98 < 1$ to:

$$ n \le \frac{\sigma D}{v} - 31. \quad (24) $$

This is a more useful value than $\frac{2}{N} \times 10^{-9} \times (\frac{\sigma D}{v})^2$ from Remark 4. Of course, which one is larger depends on the precise value of $\frac{\sigma D}{v}$ but in the regimes where this is of interest to us $\frac{\sigma D}{v}$ can be quite large. From a qualitative point of view, combining the two behaviors is even better; the goal would be to have an estimate that contains $(\frac{\sigma D}{v})^2$ with a very weak dependence on $\epsilon$, for instance logarithmic or even weaker. This will be provided in the context of the normal approximation presented in Section 2.6 and using Hoeffding inequality in the next result.

**Proposition 6.** Let $\epsilon > 0$. Then, with the previous notations, for even values of $N$:

$$ \mathbb{P}[\{\exists t \le t_n : x(t) \notin \mathcal{X}\}] \le \epsilon \text{ if } 2n \left[ \ln(n) + \ln\left(\frac{2}{\epsilon}\right) \right]^2 \ln\left(\frac{8}{\epsilon}\right) \le \frac{D^2 \sigma^2}{v^2}. \quad (25) $$

*Proof.* We invoke (12) from Proposition 3 for $Z_j = a_j \tau_j$, $C = D/v$. We choose some level to be defined $\xi > 0$ and recall that for $U \simeq \text{Exp}(1)$: $\mathbb{P}[U \ge \xi] = e^{-\xi}$ and thus $\mathbb{P}[\tau_i \ge \xi/\sigma] = e^{-\xi}$. We can replace each $\tau_i$ by $\tilde{\tau}_i = \min\{\tau_i, \xi/\sigma\}$. Note that $\mathbb{P}[\tilde{\tau}_i = \tau_i] = 1 - e^{-\xi}$ and $\mathbb{P}[\tilde{\tau}_i = \tau_i, \forall i \le n] \ge 1 - ne^{-\xi}$.

Consider the modified trajectory $\tilde{x}(t)$ constructed as $x(t)$ but with $\tilde{\tau}_i$ instead of $\tau_i$. Previous estimation informs that $\mathbb{P}[\exists t \le t_n : \tilde{x}(t) \ne x(t)] \le ne^{-\xi}$. We will use the following Hoeffding inequality: if $\chi_1, \dots, \chi_n$ are i.i.d. random variables such that $\alpha_i \le \chi_i \le \beta_i$ almost surely, then $\mathbb{P}[|\sum_{i=1}^n (\chi_i - \mathbb{E}[\chi_i])| \ge m] \le 2 \exp(-\frac{2m^2}{\sum_{i=1}^n (\beta_i - \alpha_i)^2})$. We use this Hoeffding inequality for the variables $\tilde{Z}_j = a_j \tilde{\tau}_j$ which have values inside $[-\xi/\sigma, \xi/\sigma]$ and average to zero:

$$ \begin{aligned} \mathbb{P}[\{\exists t \le t_n : x(t) \notin \mathcal{X}\}] &\le ne^{-\xi} + \mathbb{P}[\{\exists t \le t_n : \tilde{x}(t) \notin \mathcal{X}\}] \le ne^{-\xi} + 2\mathbb{P}\left[\left|\sum_{j=1}^n \tilde{Z}_j\right| \ge \frac{D}{v}\right] \\ &\le ne^{-\xi} + 4e^{-\frac{2D^2\sigma^2}{4nv^2\xi^2}}. \end{aligned} \quad (26) $$

We now choose $\xi$ such that $ne^{-\xi} = \epsilon/2$ i.e., $\xi = \ln(n) + \ln(2/\epsilon)$. Then, it is enough to find $n$ such that $4e^{-\frac{2D^2\sigma^2}{4nv^2\xi^2}} \le \epsilon/2$ to conclude. Replacing the value of $\xi$ we obtain the required estimation. $\square$

This result is to be used when $\frac{D\sigma}{v}$ is large. For example, the bound in equation (24) is better for values of $\frac{D\sigma}{v}$ up to $\simeq 5 \times 10^4$. On the contrary, for $\frac{D\sigma}{v} \simeq 10^6$, estimation (25) gives a "n" one order of magnitude larger than that of (24); for $\frac{D\sigma}{v} \simeq 10^7$ more than two orders of magnitude are gained.
<table><caption>TABLE 1. Summary of estimations of the number of steps. We used that $\zeta^2 = N/2$ for $\mathcal{D} = \mathcal{D}_N$.</caption><thead><tr><th>Situation</th><th>Time</th><th>Condition</th><th>Spatial</th><th>Condition</th></tr></thead><tbody><tr><td>Normal approximation</td><td>Equation (5):<br>$n \le T\sigma + 18 - 6\sqrt{T\sigma} + 9$</td><td>$n \ge n_N := 300$<br>(empirical)</td><td>Equation (28): $n \le \frac{1}{25N} \left(\frac{D\sigma}{v}\right)^2$</td><td>$n \ge n_N := 300$<br>(empirical)</td></tr><tr><td>Conservative (probability $\ge 1 - 10^{-9}$)</td><td>Equation (7): $n \le T\sigma/2$</td><td>$n \ge 50$</td><td>Equations (21) and (24): $n \le \frac{\sigma D}{v^2} - 31$<br>Equation (25): $n[\ln(n) + 21.42]^2 \le \frac{1}{45.61} \frac{D^2\sigma^2}{v^2}$</td><td>$n \ge 1$</td></tr></tbody></table>

## 2.6. Spatial boundary treatment, normal approximation

We consider now the normal approximation. Then in this case $\sum_{j=1}^n a_j \tau_j$ looks like a Brownian motion and, in view of the estimation (12) if we accept and error $\epsilon$, we will set $n$ such that:

$$2\mathbb{P}\left[\left|\sum_{j=1}^{n} a_j \tau_j\right| \geq D/v\right] \leq \epsilon. \qquad (27)$$

But $\sum_{j=1}^n a_j \tau_j$ has mean 0 and variance $n\zeta^2/\sigma^2$ so the probability to be larger than $D/v$ will be small when $\frac{D\sigma}{v\zeta\sqrt{n}}$ is large enough (of the order of 6), giving the aggregation rule:

$$n \le \frac{1}{50} \left( \frac{D\sigma}{v\zeta} \right)^2. \qquad (28)$$

As before, this rule will be used when $n$ is large enough, we take $n$ larger than some $n_N$ that is set to $n_N := 300$. Note that the dependence with respect to $D\sigma$ is quadratic and not linear as in Corollary 21 (recall that $\zeta = 1$ when $\mathcal{D} = \{-1, 1\}$). The dependence on the tolerance $\epsilon$ is very weak, because the $1 - \epsilon/2$ quantile of the normal law increase very slowly when $\epsilon \to 0$.

## 2.7. Summary of the aggregation rules

We summarize the estimations obtained in Table 1.

## 2.8. Sampling once aggregation level $n$ is given

Previous considerations helped to choose a $n$ such that $n$ collisions later the particle has a very high probability to be still in the space-time domain $\mathcal{X} \times [0, T]$. But we still need to place the particle somewhere i.e., we need to explain how to sample from the distribution of the position and time coordinates of the particle after $n$ collision steps.

As soon as $n > 1$ the direction of the particle is very easy to sample: it is a value from $\mathcal{D}$ sampled uniformly. For the time $t_n$ we take a gamma variable with shape $n$ and scale $1/\sigma$ (average $= n/\sigma$). What about the position $x(t_n)$? We know that $x(t_n) = x_0 + \sum_{j=1}^n a_j \tau_j$. To compute this efficiently we sample from the multinomial distribution $\mathcal{M}(1/N, \dots, 1/N)$ with $N$ events and obtain positive integers $n_j$ that sum up to $n$ with the convention that $n_j$ will represent the number of times the value $d_j \in \mathcal{D}$ has been taken by some $a_\ell$. So finally the term $\sum_{j=1}^n a_j \tau_j$ will be a sum of $N$ gamma distributions, with $j$-th term being of parameters $n_j$ and $d_j/\sigma$. The resulting procedure is described in the Algorithm 1.
**Algorithm 1.** Aggregating collision escape sampling algorithm, normal approximation.

**Inputs:** no. directions $N$, spatial size $L$, initial position $x_0$, total time $T_f$, collision parameter $\sigma$, speed $v$, approximation type apr $\in$ {"normal", "notnormal"}.

**Outputs:** space-time exit point $x^*, a^*, t^*$, number of collisions $c$

1: **procedure** SINGLE COLLISION($N, L, x, t, T_f, \sigma, v$)

2: set $\mathcal{X} = [0, L]$

3: sample $\tau \sim \text{Exp}(1/\sigma)$, $a$ uniform from $\mathcal{D}$

4: **if** $t + \tau < T_f$ and $x + va\tau \in \mathcal{X}$ **then**

5: return $(t + \tau, x + va\tau, a$, "not escaped")

6: **else**

7: compute time to escape when either $T_f$, or $\partial\mathcal{X}$ are reached: $\Delta t^* = \min\{\tau, T_f - t, \frac{\text{dist}(x,\partial\mathcal{X})}{va}\}$

8: return $(t + \Delta t^*, x + va\Delta t^*, a$, "escaped")

9: **end if**

10: **end procedure**

11: **procedure** ESCAPE TIME SAMPLING($N, L, x_0, T_t, \sigma, v$)

12: Set $t = 0, x = x_0$, collision counter $c = 0$, flag="not escaped",

13: **if** apr=="normal" **then** $n_{\mathcal{N}} = 300$

14: **else** $n_{\mathcal{N}} = 50$

15: **end if**

16: compute $d_k = \cos(2\pi j/N)$, $j = 0, \dots, N - 1$.

17: **while** flag="not escaped" **do**

18: set $D := \min\{x, L - x\}, T = T_t - t$.

19: **if** apr=="normal" **then** $n_c = \min\{T\sigma + 18 - 6\sqrt{T\sigma + 9}, \frac{1}{25N}\left(\frac{D\sigma}{v\zeta}\right)^2\}$

20: **else** $n_c = \max\{n \le T\sigma/2 | n \le \frac{\sigma D}{2} - 31 \text{ or } n[\ln(n) + 21.42]^2 \le \frac{1}{45.61} \frac{D^2\sigma^2}{v^2}\}$

21: **end if**

22: **if** $n_c \ge n_{\mathcal{N}}$ **then**:

23: sample $\tau \sim \Gamma(n, 1)$, $n_1, \dots, n_N \sim \mathcal{M}(1/N, \dots, 1/N)$, $\tau_\ell \sim \Gamma(n_\ell, 1)$, $\ell = 1, \dots, N$

24: **if** $t + \tau/\sigma > T$ or $x + \sum_{\ell=1}^N d_\ell \tau_\ell/\sigma \notin ]0, L[$ **then**

25: print "error in" apr "approximation"

26: $(t, x, a, \text{flag}) = \text{SINGLE COLLISION}(N, L, x, t, t + T, \sigma, v)$.

27: $c = c + 1$

28: **else**

29: set $t = t + \tau/\sigma$, $x = x + \sum_{\ell=1}^N d_\ell \tau_\ell/\sigma$

30: $c = c + n_c$

31: **end if**

32: **else**

33: $(t, x, a, \text{flag}) = \text{SINGLE COLLISION}(N, L, x, t, t + T, \sigma, v)$.

34: $c = c + 1$

35: **end if**

36: **end while**

37: return $(t, x, a)$, $c$

38: **end procedure**

# 3. NUMERICAL RESULTS

## 3.1. One dimension two directions

We start with the 1D case where the number of directions is $N = 2$, i.e., the $S_2$ model. This model is interesting in itself and not necessarily as a discretization of the situation when $\mathcal{D}$ is the unit sphere (see the introduction).

For the numerical tests, we set the segment length to $L = 0.01$ cm, final time to $T = 4 \times 10^4$ fs, speed $v = 3.0 \times 10^{-5}$ fs/cm, and $\sigma \in \{10^{-2}, 10, 10^4, 10^6\}$ (unit is cm$^{-1}$). The values have been chosen to see all
<table><caption>TABLE 2. Values of $T\sigma$ and $D\sigma/v$ for simulations in Section 3.1.</caption><thead><tr><th></th><th>σ = 10<sup>-2</sup></th><th>σ = 10<sup>1</sup></th><th>σ = 10<sup>4</sup></th><th>σ = 10<sup>6</sup></th></tr></thead><tbody><tr><th>Tσ</th><td>4 × 10<sup>2</sup></td><td>4 × 10<sup>5</sup></td><td>4 × 10<sup>7</sup></td><td>4 × 10<sup>10</sup></td></tr><tr><th>Dσ/v</th><td><sup>10</sup>&frasl;<sub>6</sub></td><td><sup>10<sup>4</sup></sup>&frasl;<sub>6</sub></td><td><sup>10<sup>7</sup></sup>&frasl;<sub>6</sub></td><td><sup>10<sup>11</sup></sup>&frasl;<sub>6</sub></td></tr></tbody></table>

regimes (see also results below): for $\sigma = 10^{-2}$ there are very few collisions while for $10^6$ there are a large number of them. Then the sampling inside these bounds is to obtain an almost log-uniform $T\sigma$ grid (mean number of collisions). The initial point is in the middle of the spatial interval ($x_0 = L/2$) thus $D = L/2$. We plot in Figure 2 some examples of trajectories and in Figure 3 the resulting escape laws. The values of the parameters $T\sigma$ and $D\sigma/v$ at the initial time are given in Table 2.

For $\sigma = 10^{-2}$ (first row of Figs. 2 and 3) the particle has very few collisions (here 2) and exists the spatial domain well in advance of the final time $T_f$. The aggregation mechanism is never activated. The escape distributions concentrate on the borders $x = 0$ and $x = L$; the middle column depicts the escape law conditional to having escaped because final time is reached; here this column has no mass at all.

For $\sigma = 10^1$ (second row of Figs. 2 and 3) the computation can be accelerated by aggregating several collisions, up to $\simeq 4 \times 10^4$ for an average of 399918/187 = 2138 collisions per step. The particular trajectory depicted here exists because final time is reached but $\simeq 14\%$ of trajectories exit because the spatial border is touched first.

For $\sigma = 10^4$ and $\sigma = 10^6$ (last two rows of Figs. 2 and 3) the particle undergoes an important number of collisions. The proposed procedure turns out to be very useful with a number of collisions per aggregated step being $\frac{399987252}{266} = 1.503 \times 10^6$ and $\frac{39999999269}{118} = 3.390 \times 10^8$ respectively. The procedure reaches thus acceleration factors of up to $10^8$ and the numerical resolution is very expensive without it. All particles exit because the time is up, none exists through the spatial borders: the particle behavior is that of a random walk around the initial point.

## 3.2. One dimension, $N$ directions

We take now $N > 2$ directions, i.e., the $S_N$ model. Such a model can for instance be used as a discretization of the situation when $D$ is the whole unit sphere. We will not test the same things as before but instead investigate the two possible sources of error: the fact that $N$ is finite (so $D_N$ is not a perfect representation of the unit sphere) and the fact that we used the normal approximation. Moreover we will take the most difficult test case which is the situation of a long time and a moderate value of $\sigma$, neither in the ballistic (small $\sigma$) nor in the diffusion (large $\sigma$) regime: $T = 2 \times 10^4$ fs, $\sigma = 1$. The initial point is in the middle of the spatial interval ($x_0 = L/2$).

We consider several test cases: parameters $N = 10$ can be or $N = 100$; the $n_N$ parameter (that decides when the normal approximation is to be used) can be $n_N = 50$, $n_N = 300$ or $n_N = 1000$. In each case we compute $10^4$ escape points $t, x, a$. The nominal values are: $n_N = 300$ cf. Table 2 and $N = 100$. The results for these nominal parameters are given in Figures 4 and 5. It is seen in Figure 4 that in this case all escape sides are populated, i.e., the particle can escape either through $x = 0$ or $t = T$ or $x = L$. As expected the situation is symmetric and this is confirmed by the fact that the probability to exit though $x = 0$ is the same as that for $x = L$ (up to precision $10^{-2}$). The Figure 5 presents an example of trajectory. It is seen that the aggregation procedure has been effective because it reduced the total number of steps as can be seen comparing the total number of collisions with the actual number of steps taken. As expected, when the trajectory (the leftmost plot) approaches the boundary the collisions are treated one by one and when the trajectory is close to the middle of the interval circa 600 collisions are aggregated together.

We make now variations around the nominal parameters. The results for $n_N = 1000$ and $N = 100$ are given in Figures 6 and 7. Here $n_N$ is large enough so that all collisions are treated alone, no aggregation is used so the result is the one that we could have obtained without any aggregation procedure, cf. Figure 7 middle where the
<figure><img src="image_2.png" /><figcaption>FIGURE 2. Examples of trajectories for the 1D $S_2$ test case with $L = 0.01$ cm, $T_f = 4 \times 10^4$ fs, $v = 3.0 \times 10^{-5}$ fs/cm, $\sigma \in \{10^{-2}, 10, 10^4, 10^6\}$ (cm$^{-1}$), $x_0 = L/2$, initial direction +1 (oriented towards right). Each row is a single trajectory, each trajectory corresponding to a different value of $\sigma$. First column plots the particle evolution $x(t)$ with ordinate axis being the time and abscissas the location in the $[0, L]$ segment. Second column displays the histogram of the $n_c$ values (cf. Algorithm 1). The third column gives information on the total steps taken and total number of collisions treated. The abscissas are the number of collisions and the ordinate axis is time in femtoseconds. The particle in the first row has an exit direction $a^* = +1$ because no collision takes place. The particles in the last three rows have an exit direction $a^*$ that will be drawn at random from $\{-1, 1\}$ because in all cases exit is due to the fact that time was consumed.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIGURE 3. The escape laws corresponding to Figure 2. Same convention for the signification of the rows is used. For the middle column the abscissas are time values expressed in femtoseconds; for all other plots the abscissas are *x* values in the segment [0, *L*]. Each column is a conditional distribution: in the first column are drawn histograms of *t*\* for escape points with *x*(*t*\*) = 0 (escape through the left extremity of the domain), second column draws histograms of *x*\* = *x*(*T*) for escape points with *t* = *T* (escape because time is consumed) and third histograms of *t*\* for escape points with *x*(*t*\*) = *L* (escape through the right extremity of the domain). The probability of escape for each alternative is given above the plot. When a histogram is void there is no escape for the alternative the histogram is supposed to represent. For instance, in the third row, first column one expects to see the histogram of escape time *t*\* conditioned by the fact that escape occurred through the left side i.e., *x*(*t*\*) = 0 (before time *T* and before reaching the right side). But this never happens because in this case *σ* is too large, many collisions occur and the particle does not have the time to reach the left side (the right side neither in fact) before time *T*.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 4. Results for the setting in Section 3.2 for $n_N = 300$, $N = 100$. The histograms of the escape distributions $t, x, a$ are plotted as follows: first column are histograms of the time values $t^*$, second column the positions $x^*$ and third the directions $a^*$. Each row is a conditional distribution: the first row are escape points with $x^*(t^*) = 0$ (escape through left side), second row with $t^* = T$ (escape due to total time $T$ being consumed) and third with $x^*(t^*) = L$ (escape through the right side).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 5. One example trajectory corresponding to results in Figure 4, i.e., $n_N = 300$, $N = 100$. Left plot: the trajectory $x(t)$ of the particle, $y$ axis is the time. Exit direction is $a^* = -1$. Middle plot: the distribution of the number of collisions that have been aggregated. Right plot: the number of collisions with time in the $y$ axis. Many are equal to 1 but some go as high as 600–800.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 6. Results for the setting in Section 3.2 for $n_N = 1000$, $N = 100$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 7. One example trajectory corresponding to results in Figure 6, i.e., $n_N = 1000$, $N = 100$. Signification of the plots is as in Figure 5.</figcaption></figure>

histogram is a Dirac mass in 1. On the other hand Figure 6 shows that the escape distributions are very similar to that in Figure 4 which shows that the aggregation procedure obtains comparable quality at lower costs. We have also tested $n_N = 50$ and $N = 100$ and the results (not shown here) are the same as in Figure 4 showing that our estimate $n_N = 300$ is actually a conservative one.

As a final test we will lower the parameter $N$ to inquire whether the escape distributions are sensitive to the discretization of the number of directions in $\mathcal{D}$. We take $N = 10$, $n_N = 300$, results are in Figures 8 and 9. Of course, the distribution in the third column of Figure 8 is very discrete (only 10 directions are possible and some of them coincide). But the reassuring result is that the other two columns look very much like that in
<figure><img src="image_8.png" /><figcaption>FIGURE 8. Results for the setting in Section 3.2 for $n_N = 300$, $N = 10$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 9. One example trajectory corresponding to results in Figure 8, i.e., $n_N = 300$, $N = 10$. Exit direction is $a^* = -1$. Signification of the plots is as in Figure 5.</figcaption></figure>

Figure 4 which shows that the discrete nature of $\mathcal{D}$ does not seem to play an important role in the shape of the escape distributions.

4. CONCLUSIONS

We presented a method that accelerates the sampling of the escape times, position and directions for particles undergoing collisions separated by exponentially long times. The procedure works by aggregating several collisions into a single step. The advantage of the method is that it does not uses any model or Brownian walk approximation and therefore can treat in the same way a large range of collision parameters including those where the Brownian walk approximation (as used in [4, 9]) may not work well. The procedure operates by esti-
mating conservatively the number of collisions that can safely be made before escaping the time-space domain; once this number estimated, the resulting position after those collisions is sampled exactly. The empirical results show that the method shows promising results for substantially diminishing the numerical cost while retaining excellent quality for the escape distribution.

The method was demonstrated in one dimension. Extensions to multi-dimensional domains or more general forms are left for future works. We have reasons to be optimistic because, in the worse case scenario, when the domain is a rectangle in high dimensions one can work by projecting on each of the dimensions and the setting will be that of a 1D segment; the number of aggregate collisions will be the minimum over each dimension; when the domain is not a rectangle one can consider the largest rectangle included in the domain. Finding even more efficient approaches in high dimensions is left for future work.

## REFERENCES

[1] M.I. Cabrera, O.M. Alfano and A.E. Cassano, Absorption and scattering coefficients of titanium dioxide particulate suspensions in water. *J. Phys. Chem.* **100** (1996) 20043-20050.

[2] G. Li, C. Rosenthal and H. Rabitz, High dimensional model representations. *J. Phys. Chem. A* **105** (2001) 7765-7777.

[3] J. Shorter, P. Ip and H. Rabitz, Radiation transport simulation by means of a fully equivalent operational model. *Geophys. Res. Lett.* **27** (2000) 3485-3488.

[4] J.A. Fleck and E.H. Canfield, A random walk procedure for improving the computational efficiency of the implicit Monte Carlo method for nonlinear radiation transport. *J. Comput. Phys.* **54** (1984) 508-523.

[5] J.A. Fleck and J.D. Cummings, An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport. *J. Comput. Phys.* **8** (1971) 313-342.

[6] E.D. Brooks III, Symbolic implicit Monte Carlo. *J. Comput. Phys.* **83** (1989) 433-446.

[7] B. Lapeyre, É. Pardoux and R. Sentis, Introduction to Monte-Carlo Methods for Transport and Diffusion Equations. Oxford University Press (2003).

[8] E.W. Larsen, Diffusion theory as an asymptotic limit of transport theory for nearly critical systems with small mean free paths. *Ann. Nucl. Energy* **7** (1980) 249-255.

[9] J. Giorla and R. Sentis, A random walk method for solving radiative transfer equations. *J. Comput. Phys.* **70** (1987) 145-165.

[10] M.F. Modest and S. Mazumder, Radiative Heat Transfer. Academic Press (2021).

<figure><img src="image_10.png" /></figure>

**Please help to maintain this journal in open access!**

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.