# A general Monte Carlo sampling strategy for the adjoint Boltzmann equation in continuous energy

Paul Rovel, Coline Larmier, Davide Mancusi, Andrea Zoia

## ► To cite this version:

Paul Rovel, Coline Larmier, Davide Mancusi, Andrea Zoia. A general Monte Carlo sampling strategy for the adjoint Boltzmann equation in continuous energy. 2025. cea-05250148

HAL Id: cea-05250148

https://cea.hal.science/cea-05250148v1

Preprint submitted on 11 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
A general Monte Carlo sampling strategy for the adjoint

Boltzmann equation in continuous energy

Paul Rovel,*,a,b Coline Larmier,ª Davide Mancusi,ª and Andrea Zoiaª

<sup>a</sup>Université Paris-Saclay, CEA,
Service d'étude des réacteurs et de mathématiques appliquées
91191 Gif-sur-Yvette, France

b École nationale des ponts et chaussées, Institut Polytechnique de Paris,
77455 Champs-sur-Marne, France

*Email: paul.rovel@cea.fr

Number of pages: 67

Number of tables: 2

Number of figures: 17
Abstract

For neutron transport problems where the detector region is ‘smaller’ than the source region,
solving the adjoint Boltzmann equation by Monte Carlo methods might prove a more efficient
computational strategy, compared to solving the standard (forward) Boltzmann equation. Despite
intensive research efforts over the last decades, a robust and effective implementation of adjoint
Monte Carlo sampling strategies in continuous-energy configurations remains elusive, due to in-
herent issues encompassing numerical instabilities and explosive variance growth. In view of these
considerations, in this paper we investigate the reasons behind such behavior and propose a gen-
eral sampling strategy that can be tailored in order to achieve a satisfactory compromise between
computational complexity and resulting variance of the sought observables. We focus in particular
on the sampling of the collision events, in the context of neutron transport in infinite media, as
a stepping stone towards the treatment of the full adjoint Boltzmann equations by Monte Carlo
methods.

**Keywords** — Monte Carlo, Adjoint, Transport, Neutron, Variance reduction
# I. INTRODUCTION

Solving the adjoint Boltzmann equation by Monte Carlo methods represents a long-standing goal in the field of radiation shielding since the late 1960s [1, 2, 3], following the pioneering work by Maynard [4]. The main motivation is related to the fact that any detector response due to a particle source can be equivalently expressed as a scalar product involving solutions of either the direct or adjoint Boltzmann equations, thanks to the elegant reciprocity theorem [4]. In view of the fact that the adjoint Boltzmann equation can be identified with the equation for the particle importance with respect to a given detector [5, 6], it is in principle possible to conceive Monte Carlo methods to sample the histories of ‘importance particles’ and thus estimate the solutions of such equation.

This paves the way to estimating the sought response function by using an alternative strategy to the standard sampling of (forward) particle histories. In forward Monte Carlo simulations, neutron histories begin from the source, undergo flight and collision events, and accumulate their contributions in the detector region. Conversely, adjoint Monte Carlo simulations sample a reversed process: adjoint particles (sometimes called adjunctons [7]) begin from the detector, propagate backwards undergoing reversed flights and collisions, and finally accumulate their contributions in the source region [8].

The adjoint Monte Carlo approach might become computationally advantageous (as opposed to regular, forward sampling) whenever the detector region is ‘smaller’ than the source region, and in particular in the limit case of point detectors [8]. Indeed, in such cases forward Monte Carlo simulations will struggle to reach a statistically significant number of events contributing to the response, whereas adjoint Monte Carlo simulations would produce more histories contributing to the response and thus ensure a faster statistical convergence of the sought quantity.

The formalism of adjoint Monte Carlo methods was further developed starting from the 1970s, and several distinct algorithms were proposed [9, 10, 7, 11, 8, 12]. Extensive numerical investigations were carried out, especially in the context of continuous-energy neutron transport problems [13, 14, 15]: this revealed the occurrence of serious issues, encompassing numerical instabilities and explosive variance growth. The core of the difficulty lies in sampling the adjoint collision events, which are inherently different from their forward counterpart [14]. Adjoint flights are comparatively simpler, due to their similarity with the regular flights [14, 16]. Thus, despite the
intensive research efforts carried out in the past, the foundations of adjoint Monte Carlo methods
are still a topic of active investigation [17, 18, 16, 19, 20]. Furthermore, very few general-purpose
neutron-transport Monte Carlo codes implemented the adjoint sampling algorithms in continuous-
energy, including the full spectrum of nuclear reactions that are available in modern nuclear data
libraries [21, 22].

In view of these considerations, in this paper we investigate the statistical behavior of Monte
Carlo methods for adjoint particle transport, bearing in mind the application to radiation shielding
problems in continuous-energy settings. For the sake of simplicity, we will focus on stationary
fixed-source configurations in infinite media, where the only relevant sampling strategy concerns
the adjoint collision events.

This manuscript is organized as follows. First, in Section II we will briefly introduce the
forward and adjoint Boltzmann equations, and discuss the physical meaning of the forward and
adjoint particle densities. Then, in Section III we will focus on how the adjoint equations can
be solved by the Monte Carlo method, in the context of infinite media where the only phase-
space variable is the particle energy. We will revisit the physical interpretation of the adjoint
transport process, and elucidate some of its seemingly counterintuitive behaviors. Within this
broad framework, in Section IV we will propose several Monte Carlo sampling strategies for the
adjoint collision events, compare them to those previously proposed in the literature, and illustrate
their respective advantages and shortcomings. Special care will be given to techniques tailored to
reduce the variance of the sought response. In order to substantiate our findings, we will apply
these sampling strategies to the resolution of a few relevant adjoint problems in Section V, focusing
on the effects of different sampling strategies on the variance of the obtained results. Conclusions
will be finally drawn in Section VI. Additional technical details will be provided in a series of
Appendices.

## II. DIRECT AND ADJOINT TRANSPORT EQUATIONS

### II.A. Particle densities and the direct transport equations

In fixed-source time-independent particle-transport problems, the sought response function $R$ is typically a reaction rate or a dose. Generally speaking, quantity $R$ can be expressed in terms
of the directª collision density $\psi$ (i.e., the average number of collisions per unit of phase space) or the direct emission density $\chi$ (i.e., the average number of particles emitted per unit of phase space), namely

$$R = \int \eta_{\psi}(P) \psi(P) dP = \int \eta_{\chi}(P) \chi(P) dP, \qquad (1)$$

where $\eta_{\psi}$ and $\eta_{\chi}$ are the detector response functions associated to the collision and emission density, respectively. The shorthand notation $P = (\vec{r}, E, \vec{\Omega})$ denotes a phase-space point with position $\vec{r}$, energy $E$ and direction $\vec{\Omega}$. Unless stated otherwise, integration is performed over the entire phase space. In some cases, the angular particle flux $\phi$ is needed, which is related to $\psi$ by the algebraic relation $\phi(P) = \psi(P)/\Sigma_t(\vec{r}, E)$, the quantity $\Sigma_t(\vec{r}, E)$ expressing the total macroscopic cross section. As stated in the introduction, in this paper we focus on neutron transport, and all physical quantities will thus refer to neutron-matter interactions; however, the discussion in the following might apply more broadly to neutral particle transport, with minimal modifications.

The densities $\psi$ and $\chi$ satisfy the transport (Boltzmann) equation, which basically expresses a particle conservation law. An external particle source $S(P)$ is assigned, corresponding to the average number of particles externally added to the system around point $P$. Once born from the source, particles repeatedly undergo a series of successive flight and collision events. Flights correspond to particles streaming freely until undergoing a collision event, or leaving the viable domain. The flight kernel $T(P \to P')$ represents the average number of particles (per unit phase space volume) having their next collision around $P'$, for a single particle starting a flight at $P$. The kernel $T$ is related to $\chi$ and $\psi$ by

$$\psi(P') = \int T(P \to P')\chi(P)dP. \qquad (2)$$

Equation (2) states that the average number of particles entering collision around $P'$ is obtained by summing the particles being emitted at any point $P$ in phase space and undergoing their next collision around $P'$.

Upon collision at point $P$, the particle will be re-emitted at another point $P'$. The collision can result in zero, one or more particles being emitted for each incident particle. The collision kernel $C(P \to P')$ represents the average number of particles (per unit phase space volume) leaving

<footnote>ªWe will interchangeably use ‘forward’ and ‘direct’ in the following.</footnote>
a collision around $P'$, for a single particle entering a collision at $P$. The kernel $C$ is related to $\chi$, $\psi$ and $S$ by

$$ \chi(P') = \int C(P \to P')\psi(P)dP + S(P'). \quad (3) $$

Equation (3) states that the average number of particles being emitted around $P'$ is the sum of those coming from the external particle source $S(P')$ and those having a collision at any point $P$ and leaving the collision around $P'$.

Combining Eqs. (2) and (3), we obtain the (direct) Boltzmann equation for $\chi(P)$:

$$ \chi(P'') = \int C(P' \to P'') \int T(P \to P')\chi(P)dPdP' + S(P''), \quad (4) $$

which can be equivalently expressed in terms of $\psi(P)$ as

$$ \psi(P'') = \int T(P' \to P'') \int C(P \to P')\psi(P)dPdP' + \int T(P' \to P'')S(P')dP'. \quad (5) $$

The asymmetry of Eq. (4) and Eq. (5) with respect to the source term is due to $S$ specifically creating ‘emitted’ particles, and thus directly contributing to $\chi$ and not to $\psi$. However, one could imagine a special source that would create particles entering collision, in which case such a source would directly contribute to $\psi$ instead of $\chi$, and the situation would be inverted.

For the sake of conciseness, we define the flight operator

$$ \mathcal{T}f(P) = \int T(P' \to P)f(P')dP' \quad (6) $$

and the collision operator

$$ \mathcal{C}f(P) = \int C(P' \to P)f(P')dP', \quad (7) $$

for any sufficiently well-behaved $f(P)$. Then, under suitable assumptions, Eq. (5) can be solved for $\psi$ using a Neumann series, as

$$ \psi = [\mathcal{T} + \mathcal{T}\mathcal{C}\mathcal{T} + (\mathcal{T}\mathcal{C})^2\mathcal{T} + (\mathcal{T}\mathcal{C})^3\mathcal{T} + \dots] S = \sum_{n=0}^{\infty} (\mathcal{T}\mathcal{C})^n \mathcal{T}S = \mathcal{K}S, \quad (8) $$

which means that $\psi$ is obtained by successively applying the flight operator $\mathcal{T}$ and collision operator $\mathcal{C}$ to the source $S$. This whole procedure is encapsulated in the propagator, denoted $\mathcal{K}$, being
applied to the source S to obtain $\psi$. A similar solution exists for Eq. (4), yielding $\chi$ by repeated application of the operators $T$ and $C$.

## II.B. Importance and the adjoint transport equations

In addition to the direct particle densities, one can also determine the ‘importance’ of an event in phase space with respect to a given detector. More precisely, we introduce the collision importance $\psi^\dagger(P)$ as the expected contribution to the response $R$ for a neutron entering a collision at point $P$. Similarly, we introduce the emission importance $\chi^\dagger(P)$ as the average contribution to $R$ for a neutron being emitted at point $P$. In this work, we will follow the convention that the detector is expressed via the collision response function $\eta_\psi$; note, however, that one could introduce a broader definition of importance, relative to a more general response function.

The importance functions $\psi^\dagger(P)$ and $\chi^\dagger(P)$ carry dimensions of [unit of contribution] [particle]$^{-1}$ and are therefore not phase space densities, contrary to the particle densities $\psi(P)$ and $\chi(P)$. We will discuss in greater depth the consequences of this observation in Section III.D, in connection with the Monte Carlo interpretation of the direct and adjoint transport equations.

The quantities $\psi^\dagger(P)$ and $\chi^\dagger(P)$ are related to each other and satisfy conservation laws that are determined by examining the future contributions of a given particle to the detector. Let us consider the importance $\psi^\dagger(P)$ of a particle entering collision at $P$: the particle will directly contribute $\eta_\psi(P)$ to the response $R$ through the ongoing collision; the average number $C(P \to P')$ of particles that will emerge from the collision with outgoing coordinates $P'$ will each further contribute $\chi^\dagger(P')$ to $R$, by definition of the emission importance. We obtain then the law

$$ \psi^\dagger(P) = \int C(P \to P')\chi^\dagger(P')dP' + \eta_\psi(P), \qquad (9) $$

where we have integrated over all possible emission points $P'$. Let us now consider the importance $\chi^\dagger(P)$ of a particle emitted at $P$: the average number $T(P \to P')$ of particles ending a free flight at coordinate $P'$ will each have a collision there and thus contribute $\psi^\dagger(P')$ to $R$, by definition of the collision importance. We obtain then the law

$$ \chi^\dagger(P) = \int T(P \to P')\psi^\dagger(P')dP', \qquad (10) $$
where we have integrated over all possible collision points $P'$. Furthermore, it can be shown that the flux importance $\phi^\dagger(P)$, i.e. the importance of a particle in free-flight at $P$, is simply given by $\phi^\dagger(P) = \chi^\dagger(P)^b$. Similarly to what was done for the direct transport equation, we can combine Eqs. (9) and (10) into

$$ \psi^\dagger(P) = \int C(P \to P') \int T(P' \to P'') \psi^\dagger(P'') dP'' dP' + \eta_\psi(P), \quad (11) $$

or equivalently

$$ \chi^\dagger(P) = \int T(P \to P') \int C(P' \to P'') \chi^\dagger(P'') dP'' dP' + \int T(P \to P') \eta_\psi(P') dP'. \quad (12) $$

In Eq. (11), the detector response function $\eta_\psi(P)$ formally plays the role of a source term. Again, the lack of symmetry with respect to the source term $\eta_\psi$ is related to the fact that we have set a detector whose response function is associated to the collision density. If we had chosen a detector associated to the emission density, the situation would have been reversed.

Equations (11) and (12) for the importance involve the same kernels as Eqs. (4) and (5) for the particle densities, but the integrations are carried out ‘backward’. We formally introduce the adjoint flight operator

$$ T^\dagger f(P) = \int T(P \to P') f(P') dP' \quad (13) $$

and the adjoint collision operator

$$ C^\dagger f(P) = \int C(P \to P') f(P') dP', \quad (14) $$

for any sufficiently well-behaved $f(P)$, which formally corresponds to setting the adjoint kernels

$$ T^\dagger(P \to P') = T(P' \to P). \quad (15) $$

and

$$ C^\dagger(P \to P') = C(P' \to P). \quad (16) $$

<footnote>bThis holds true as neutrons here are only allowed to contribute when entering collision. If one were to use an emission-based detector $\eta_\chi$, the emission importance $\chi^\dagger$ would include the immediate payoff from $\eta_\chi$ while the importance of a particle in free-flight $\phi^\dagger$ would not, leading to the equation $\chi^\dagger(P) = \phi^\dagger(P) + \eta_\chi(P)$</footnote>
Therefore, under suitable assumptions, Eq. (12) can be expanded in Neumann series as

$$ \chi^{\dagger} = [T^{\dagger} + T^{\dagger}C^{\dagger}T^{\dagger} + T^{\dagger}(C^{\dagger}T^{\dagger})^{2} + T^{\dagger}(C^{\dagger}T^{\dagger})^{3} + \dots] \eta_{\psi} = \sum_{n=0}^{\infty} T^{\dagger}(C^{\dagger}T^{\dagger})^{n} \eta_{\psi} = K^{\dagger}\eta_{\psi}, \quad (17) $$

whose structure is very similar to Eq. (8): the emission importance $\chi^\dagger$ is obtained by repeatedly applying the adjoint flight operator $T^\dagger$ and the adjoint collision operator $C^\dagger$ to the detector response function $\eta_\psi$. This corresponds exactly to the application of the adjoint propagator $K^\dagger$ to the response function $\eta_\psi$. Due to the use of the adjoint propagator to obtain the importance $\chi^\dagger$, Eqs. (11) and (12) for the particle importance are usually called the adjoint transport equations, and the term $\eta_\psi$ is called the adjoint source.

Observe that the importance $\chi^\dagger(P)$ yields the expected contribution to the detector due to a single particle starting at $P$. The response $R$, which corresponds to the expected contribution of particles created by the source, can be then expressed as the importance per particle weighted by the source density $S(P)$. More formally, the emitted importance function $\chi^\dagger(P)$ and the particle density $\psi(P)$ are related to each other and to the direct and adjoint sources^c by the elegant property

$$ R = \int \eta_{\psi}(P) \psi(P) dP = \int \chi^{\dagger}(P) S(P) dP, \quad (18) $$

which stems from

$$ R = \langle \psi, \eta_{\psi} \rangle = \langle KS, \eta_{\psi} \rangle = \langle S, K^{\dagger} \eta_{\psi} \rangle = \langle S, \chi^{\dagger} \rangle. \quad (19) $$

Equation (18) takes the name of reciprocity theorem: the sought response $R$ can be equivalently computed by solving the direct equations (which determine the particle density $\psi$ in the phase space resulting from the source $S$) weighted by the detector response function $\eta_\psi$, or by solving the adjoint equations (which determine the importance function $\chi^\dagger$ in the phase space associated to the detector $\eta_\psi$) weighted by the source $S$.

<footnote>cNote that if we had used the emission-based response function $\eta_\chi$, the reciprocity theorem would have been expressed in terms of $\chi$, instead of $\psi$. Similarly, if we had used the first collision source $S_\psi(P) = \int T(P' \to P)S(P')dP'$ instead of $S$, then the reciprocity theorem would have been written in terms of the collisional importance $\psi^\dagger$.</footnote>
### III. SOLVING THE ADJOINT BOLTZMANN EQUATION VIA MONTE CARLO SAMPLING

Monte Carlo methods are customarily used to solve the direct transport equations, in view of the fact that they introduce almost none of the approximations (discretization, self-shielding models, ...) that are typically required by standard deterministic methods. This comes at the cost of a generally higher computational burden and an inherent statistical uncertainty on the obtained result, stemming from the stochastic nature of this approach. The basic Monte Carlo strategy consists in sampling the birth points of the particles based on the source $S$, and then simulating the random trajectories of the particles in the system by repeatedly sampling flight and collision events based on the kernels $T$ and $C$, respectively, until each history is terminated. Whenever a particle contributes to the detector, a tally is recorded using a given statistical estimator. Ensemble averages over all simulated histories will provide an unbiased estimate of the sought response $R$, provided that care is taken in selecting the estimator and the sampling laws for the histories.

Based on their similarity with respect to the direct Eq. (8) for the particle density, the adjoint Eq. (17) for the importance can be formally given a physical interpretation in terms of a ‘transport’ equation for ‘importance particles’, which have historically been named ‘adjunctons’ [10, 7]. According to this interpretation, adjunctons are born from the (direct) detector response function $\eta_{\psi}(P)$, which plays the role of the (adjoint) source. Then, adjunctons undergo ‘reversed’ flights and ‘reversed’ collision events according to the adjoint kernels $T^{\dagger}(P \rightarrow P')$ and $C^{\dagger}(P \rightarrow P')$, respectively, and finally contribute to the source $S$, which plays the role of the (adjoint) detector response function. We will make these statements more precise in the following section.

As briefly mentioned in the introduction, when the detector is sensitive to a relatively small volume in phase space, obtaining a sufficient number of contributing histories (and thus an acceptably small statistical uncertainty) might become a daunting task for direct Monte Carlo games. For configurations where the phase space volume of the source is much larger than that of the detector, more adjoint histories would contribute in the source than direct histories in the detector: it might be therefore tempting to use Monte Carlo in order to solve the adjoint equations. This is expected to lead to enhanced computational efficiency while preserving the lack of approximations. This strategy, despite its elegance and ostensible simplicity, leads unfortunately to several conceptual and practical issues, encompassing notably numerical instabilities, which make adjoint
Monte Carlo methods difficult in practice, as extensively pointed out by many authors over the last few decades [11, 15, 16, 19].

### III.A. Can the adjoint Boltzmann equation be interpreted as a transport equation?

Equation (8) shows that $\psi$ is obtained by successively applying the flight operator $\mathcal{T}$ and collision operator $\mathcal{C}$ to the source $S$. To sample a direct Monte Carlo game, we 'split' the source $S(P)$ into individual particles (called neutrons by analogy with the real physical process), each carrying its own portion of the total intensity of the source. These neutrons are then transported following $\mathcal{T}$ and $\mathcal{C}$, and accumulate their tallies along the way whenever they interact with the detector response $\eta_{\psi}$.

Historically, works concerned with adjoint transport have always purposefully tried to formulate the importance equation as a special form of transport equation. Equation (17) shows that $\chi^{\dagger}$ is obtained by successively applying the adjoint flight operator $\mathcal{T}^{\dagger}$ and adjoint collision operator $\mathcal{C}^{\dagger}$ to the response function $\eta_{\psi}$. By analogy with the direct problem, in order to sample an adjoint Monte Carlo game it has been suggested to 'split' the source $\eta_{\psi}(P)$ into individual 'importance particles', each carrying its own portion of the total intensity of the source $\eta_{\psi}$. Adjunctons can be then transported following $\mathcal{T}^{\dagger}$ and $\mathcal{C}^{\dagger}$, and accumulate their tallies whenever they interact with the detector response $S$ [10, 7]. By virtue of the reciprocity theorem, the tallies estimated in adjoint Monte Carlo games will indeed correspond to the sought result.

The strategy of building adjoint Monte Carlo games based on the close analogy with respect to the direct transport equations is extremely appealing in terms of building upon existing Monte Carlo codes developed for direct particle transport problems, with minimal or possibly no modifications [15]. However, a number of theoretical and algorithmic issues emerge from this approach. By inspection of Eqs. (8) and (17), it is natural to assume that the emission importance $\chi^{\dagger}$ is a collision density (and conversely that the collision importance $\psi^{\dagger}$ is an emission density). Related to this inversion, collision estimators (counting the contributions of particles entering collision) will estimate the emission importance $\chi^{\dagger}$, whereas emission estimators (counting the contributions of particles emerging from a collision) will estimate the collision importance $\psi^{\dagger}$; furthermore, the adjoint source creates adjunctons ready to perform a flight, although they carry 'collisional importance' in this state. As a consequence, in such adjoint Monte Carlo games, the collision importance
$\psi^\dagger$ must be estimated after emerging from the adjoint collision event, and not when entering the actual adjoint collision; similarly, the emission importance $\chi^\dagger$ must be estimated at the adjoint collision events, and not when starting an adjoint flight.

Another subtle point concerns the use of statistical weights in adjoint games. Because of the interpretation of $T^\dagger$ as a form of flight operator and of $C^\dagger$ as a modified collision operator, it has been shown that weight multiplication factors will appear even for non-multiplicative events, such as elastic scattering [19, 8, 17]. This issue, which will be discussed at length in the following, might lead to numerical instabilities, severely hindering the simulations or even preventing their convergence [19, 8].

To conclude, the adjoint Boltzmann equation can be interpreted as a transport equation and solved by Monte Carlo methods. This however comes at the cost of a few unusual behaviors for the transported particles.

In view of the previous considerations, we set out to conceive a general Monte Carlo sampling strategy for adjoint transport problems that can cope with the aforementioned problems and yet be simple enough to warrant implementation in general-purpose simulation codes. As a stepping stone towards this ambitious goal, in this paper we focus on the special case of the adjoint collision kernel, which physically corresponds to considering infinite homogeneous media. The case of the flight kernel is expected to involve fewer difficulties [14], and will be dealt with in a forthcoming paper.

### III.B. The adjoint collision kernel

Before examining how to sample according to the adjoint collision kernel, we will first recall a sampling strategy for the direct collision, which will serve as a reference. The direct collision kernel, as illustrated in Eq. (3), relates the collision density $\psi$ and the emission density $\chi$. Bearing a Monte Carlo interpretation in mind, it is useful to decompose the collision event into intermediate sampling steps [23]: the choice of the collided nuclide, denoted by an index $i$, the choice of the reaction channel $j$ for the collided nuclide $i$, and the choice of the outgoing states and the associated multiplicity. This is formalized by:

* the probability $N_i\sigma_{i,t}/\Sigma_t$ of interacting with a given nuclide $i$, where $N_i$ is the nuclide concentration, $\sigma_{i,t}$ is the total microscopic cross section of nuclide $i$, and $\Sigma_t$ is the total
macroscopic cross section of the considered mixture of nuclides in the traversed medium;

* the conditional probability $\sigma_{i,j}/\sigma_{i,t}$ of selecting the reaction channel $j$, given a collided nuclide $i$, where $\sigma_{i,j}$ is the microscopic cross section of reaction $j$ on nuclide $i$, and $\sigma_{i,t} = \sum_j \sigma_{i,j}$;

* the average multiplicity $\nu_{i,j}$ and the density $f_{i,j}$ of post-collision features for reaction channel $j$ on nuclide $i$.

As customary for stationary transport problems, the concentrations $N_i$ depend only on $\vec{r}$, the macroscopic cross section $\Sigma_t$ depends on $\vec{r}$ and energy $E$, the microscopic cross sections $\sigma_{i,j}$ and $\sigma_{i,t}$ depend only on energy $E$, the average multiplicity $\nu_{i,j}$ depends only on $E$ and the average collision distribution law $f_{i,j}$ depends only on the ingoing energy $E$, the outgoing energy $E'$ and the cosine of the deflection angle $\mu = \vec{\Omega} \cdot \vec{\Omega}'$. However, to keep the notation concise, we will sometimes omit to write the dependence entirely.

Correspondingly, the collision kernel will be decomposed as:

$$C(P \to P') = \sum_i \frac{N_i(\vec{r})\sigma_{i,t}(E)}{\Sigma_t(\vec{r}, E)} \sum_j \frac{\sigma_{i,j}(E)}{\sigma_{i,t}(E)} \times \nu_{i,j}(E) f_{i,j}(E \to E', \vec{\Omega} \cdot \vec{\Omega}'). \quad (20)$$

In a Monte Carlo game, each term in Eq. (20) is associated to a transition between states, as illustrated in Fig. 1, where we display all the intermediate steps visited by the particles undergoing a collision.

The scheme in Fig. 1 shows that there is a coefficient $N_i\sigma_{i,t}/\Sigma_t$ occurring in the collision kernel, between the collision density $\psi(P)$ at $P$, which collects all the particles entering collision at $P$, and the collision density at $P$ conditioned to nuclide $i$. The factors $N_i\sigma_{i,t}/\Sigma_t$ sum to one (for all possible nuclides), and can thus be used directly as a probability distribution for choosing nuclide $i$ when sampling the collision with Monte Carlo, in which case there will be no need for weight correction to preserve unbiasedness at this stage. The sampling of the reaction channel can be performed similarly, using $\sigma_{i,j}/\sigma_{i,t}$ as a distribution to choose reaction channel $j$. The law $f_{i,j}(E \to E', \vec{\Omega} \cdot \vec{\Omega}')$ is normalized, and can thus be used to sample the outgoing coordinates $E'$ and $\vec{\Omega}'$. The multiplicity $\nu_{i,j}$ is in general not equal to unit, but can be taken into account by applying a weight correction $\nu_{i,j}$ to the emitted particle, or by creating a random number of particles with average equal to $\nu_{i,j}$ (for sterile capture, with cross section $\sigma_{i,c}$, we have $\nu_{i,c} = 0$, and the particle is killed upon collision).
<figure><img src="image_2.png" /><figcaption>Fig. 1. The decomposed direct collision kernel $C$, describing all the intermediate states of the collision event, along with a possible sampling method. For the sake of illustration, the states are depicted as if they were discrete. The arrows represent the different possible transitions between states, their norm being the corresponding coefficient. When several options are possible, we have to sample the one towards which the particle is sent, among all possible states. Conversely, when many transitions lead to the same state, the value of such state is the sum of all contributions coming from all paths leading into it. Formulas in black represents the coefficients occurring in the collision kernel. Formulas in red represent a possible sampling distribution to be used by a Monte Carlo game, and formulas in purple represents the corresponding weight corrections to ensure unbiasedness. See text for detailed explanation.</figcaption></figure>

This sampling strategy for the direct collision event is not unique. For instance, implicit capture schemes inhibit the sampling of capture reaction channels, leading to $\nu_{i,c} = 0$, and correspondingly correct the outgoing particle weight by a factor of $1 - \sigma_{i,c}/\sigma_{i,t}$, to preserve the unbiasedness of the Monte Carlo game. In general, it is possible to use arbitrary distributions at each sampling step, as long as a weight correction is applied in order to ensure that the product of the probability of the event times the weight correction gives the corresponding term in the kernel, which preserves unbiasedness.

### III.B.1. Sampling of the adjoint collision kernel

The sampling of the adjoint collision event is more involved, since the adjunctons will visit the intermediary steps of the collision ‘backward’. The corresponding sampling scheme is illustrated in Fig. 2, where the sampling of the adjoint collision event is compared to that of the direct collision event.

In Fig. 2, the direct collision decomposition is further extended compared to the scheme in Fig. 1, in order to explicitly include the ‘emission by reaction $j$ on nuclide $i$’ and the ‘emission by nuclide $i$’ events as intermediary steps. Those extra steps are not usually useful for the sampling of the direct collision event, since a particle emitted at $P$ by $\{i, j\}$ will directly go to the state
<figure><img src="image_3.png" /><figcaption>Fig. 2. Decomposed direct and adjoint collision kernels, with the proposed sampling methods. Formulas in black represent the coefficients occurring in the collision kernel. Formulas in red represent a possible sampling distribution to be used by a Monte Carlo game, and formulas in purple represent the corresponding weight correction to ensure unbiasedness. The adjoint collision kernel features the same coefficients as the direct collision kernel, but the transitions between states are reversed. A general method, composed of arbitrary sampling distributions and corresponding weight corrections, is proposed to sample this adjoint collision. These distributions are stored in the form of adjoint distribution laws $f_{i,j}^{\dagger}$ and adjoint cross sections $\sigma_{i,j}^{\dagger}$, by analogy with the direct case. See text for detailed explanation.</figcaption></figure>

‘emitted at $P$ by nuclide $i$, and to the ‘emitted at $P'$ state, which justifies the unit coefficients in the kernel and the distribution probabilities. However, these states are useful when sampling the collision backward, since the choice of the nuclide $i$ and of reaction channel $j$ for the adjoint collision event is performed at these steps.

Similarly as for the sampling of the direct collision, the sampling of the adjoint collision involves the choice of a nuclide, a reaction channel, and finally the particle coordinates using the multiplicity and the post-collision features density. As illustrated in Fig. 2, the coefficients occurring in the adjoint collision kernel are the same as those in the direct collision kernel, and the only modification concerns the direction of the arrows corresponding to the transitions between
states.

The first choice to make to sample the adjoint collision is to select which nuclide and reaction channel emitted the particle. The arbitrary probability distributions needed to sample the nuclide $i$ and reaction $j$ are assigned by introducing (arbitrary and non-negative) adjoint microscopic cross sections $\sigma_{i,j}^{\dagger}$ for each nuclide $i$ and reaction channel $j$. This strategy was initially proposed in the pioneering work by Eriksson [1], and later adopted and expanded by many authors [8, 16, 7, 12, 9]. The total microscopic and macroscopic cross sections follow from the definitions $\sigma_{i,t}^{\dagger} = \sum_{j} \sigma_{i,j}^{\dagger}$ and $\Sigma_{t}^{\dagger} = \sum_{i} N_{i} \sigma_{i,t}^{\dagger}$. Building upon the adjoint cross sections, we can then formally establish the sought arbitrary probability $N_{i} \sigma_{i,t}^{\dagger} / \Sigma_{t}^{\dagger}$ for the choice of the nuclides, and $\sigma_{i,j}^{\dagger} / \sigma_{i,t}^{\dagger}$ for the choice of the reaction channels. Several distinct expressions for $\sigma_{i,j}^{\dagger}$ (and thus for the underlying arbitrary probability distributions) have been proposed in the literature: we will discuss the advantages and drawbacks of each approach in some depth in Section IV, in particular in terms of the resulting variance.

The scheme in Fig. 2 shows that all the possible transitions associated with the sampling of the nuclide and of the reaction channels come with the same unit coefficient in the kernel. In order to preserve the unbiasedness of the adjoint Monte Carlo game, a weight correction $\Sigma_t^\dagger/(N_i\sigma_{i,t}^\dagger)$ is applied after the sampling of the nuclide, and a weight correction $\sigma_{i,t}^\dagger/\sigma_{i,j}^\dagger$ is further applied after the sampling of the reaction channel, which ensures that the product of the probability distribution used with the weight correction matches the coefficient in the kernel, and thus the unbiasedness of the whole procedure [11].

The last choice concerns the selection of the entering state $P$. As illustrated in the scheme in Fig. 2, the coefficients $\nu_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}')$ occurring in the collision kernel are not normalized in $E, \vec{\Omega}$, and cannot be used for this purpose. We introduce then an arbitrary density $f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega})$ that will be used to sample the entering state. Unbiasedness will be ensured by adjusting the statistical weight of the sampled particle by the factor $\nu_{i,j} f_{i,j} / f_{i,j}^{\dagger}$. In the literature, some authors have further introduced an adjoint multiplicity $\nu_{i,j}^{\dagger}$, whose role is usually to carry all the weight corrections required by the adjoint collision kernel. This approach will not be pursued in this paper, but its advantages and drawbacks will be briefly discussed in Section III.B.2.

As illustrated in Fig. 2, in order to move to the state corresponding to the particle entering the collision at $P$ at the end of the sampling procedure, the coefficients are not unitary, although
no further sampling is required. The particle weight is thus multiplied by the factor $\sigma_{i,j}/\sigma_{i,t}$, and then by $N_i\sigma_{i,t}/\Sigma_t$.

### III.B.2. Other weight correction techniques

In our work, we have chosen to always correct the particle weight after sampling from a distribution that does not directly match the coefficient actually occurring in the collision kernel in Fig. 2. Determining precisely the steps at which all the weight corrections are applied along the adjoint collision sampling is key if one wants to estimate the importance at intermediary states, e.g. the importance of a particle entering a collision on nuclide i via reaction j. If one only wants to assess the importance before (emitted importance $\chi^\dagger$) or after (collision importance $\psi^\dagger$) the sampling of the adjoint collision event, then all the weight corrections can be factorized and applied at once.

In both cases, these techniques could be called ‘local’ weight corrections, because the bias introduced by using modified laws is corrected immediately after the sampled event, so that any score related to particle counting would be unbiased. The pioneering techniques for adjoint Monte Carlo games introduced by Irving in 1971 [10] and Hoogenboom in 1977 [8, 14] used instead a ‘global’ weight correction. This consists in adjusting the quantities actually estimated by a global factor: instead of computing $\psi^\dagger(P)$ and $\chi^\dagger(P)$, one evaluates $I(P) \times \psi^\dagger(P)$ and $I(P) \times \chi^\dagger(P)$, where $I(P)$ is an arbitrary positive function. The modified adjoint collision kernel for these I-weighted densities is:

$$ \tilde{C}^{\dagger}(P' \to P) = C^{\dagger}(P' \to P) \frac{I(P)}{I(P')}. \qquad (21) $$

These arbitrary modifications allow fitting the coefficients of the kernel to the actual distribution that we want to use in the adjoint Monte Carlo sampling to avoid weight correction. To ensure unbiased responses, one has to adapt the response function by multiplying it by $1/I(P)$, and to adapt the source term by multiplying it by $I(P)$.

For example, Hoogenboom has proposed two sampling strategies using $I(P) = \Sigma_t(\vec{r}, E)$ and $I(P) = \Sigma_t(\vec{r}, E)/E$ [8, 14]. These weight correction techniques have the same effect (in terms of explored trajectories, density of particles in the phase space, variance of the responses, and of course expected response) as the ‘local’ weight correction, provided that one uses the same
definitions of adjoint cross sections $\sigma_{i,j}^{\dagger}$ and distributions laws $f_{i,j}^{\dagger}$ d. The only difference is that the weights of the particles during their lives are not the same: differences can thus arise if some population control is enforced without adapting the thresholds.

The ‘global’ weight correction techniques also rely on the introduction of an adjoint multiplicity $\nu_{i,j}^{\dagger}$, which effectively acts as a weight correction [10, 8, 14]. The adjoint multiplicity carries all the weight correction terms that could not be removed via the ‘global’ weight correction discussed above. This leads to the following definition for $\nu_{i,j}^{\dagger}$:

$$ \nu_{i,j}^{\dagger}(\vec{r}, E', \vec{\Omega}', E, \vec{\Omega}) = \frac{\sigma_{i,j}(E)}{\Sigma_t(\vec{r}, E)} \times \frac{\nu_{i,j}(E) f_{i,j}(E \to E', \vec{\Omega} \cdot \vec{\Omega}')}{f_{i,j}^{\dagger}(E' \to E, \vec{\Omega} \cdot \vec{\Omega}')} \times \frac{\Sigma_t^{\dagger}(\vec{r}, E')}{\sigma_{i,j}^{\dagger}(E')} \times \frac{I(\vec{r}, E, \vec{\Omega})}{I(\vec{r}, E', \vec{\Omega}')} \quad (22) $$

This strategy is a perfectly valid way of taking into account the weight correction. However, in general it depends on both the entry point $P'$ and the exit point $P$. It would be desirable to define an adjoint multiplicity that depends only on entering coordinates, so that it looks similar to the classical multiplicity $\nu_{i,j}$ and can be easily pre-tabulated. For this purpose, one has to choose specific forms for $\sigma_{i,j}^{\dagger}$ and $f_{i,j}^{\dagger}$, and properly set the value of $I(P)$ so that all dependence on $\vec{r}, E, \vec{\Omega}$ and $\vec{\Omega}'$ cancels.

Using the global weight correction and the adjoint multiplicities, it might be feasible to run an adjoint calculation by changing only the nuclear data library, the sources and response function, without any modification to the Monte Carlo code itself. This has been successfully achieved for simple cases [15], where the adjoint data is not too hard to compute and manipulate, but this approach lacks the flexibility of the local weight correction technique.

The local weight correction proposed in this work can be seen as a generalization of the global weight correction, and offers many advantages. First, it can be computed and applied after selecting the exit point $P'$, thus allowing for more flexibility. Second, the local weight correction can help in compensating discretization errors that might occur in the production of adjoint data, as was first pointed out by De Matteis and Simonini in 1978 [11].

<footnote>dThe corresponding definitions will be given in Section IV.B.1, when treating the different possible sampling strategies.</footnote>
### III.C. Source and detectors

The sampling of adjoint transport problems further requires taking into account the source and the detectors.

#### III.C.1. Adjoint sources

The adjunctons are first sampled at the detector region, which plays the role of the adjoint source. The particles are created in the ‘emitted particle’ state (which means that they will go through the adjoint flight kernel first), although they represent collisional importance. The total norm of the adjoint source, i.e., the sum of all the weights of the sampled adjunctons, also needs to correspond to the total intensity of the detector response function, i.e. the integral of the response function in phase space. Thus, we have

$$S^{\dagger}(P) = \eta_{\psi}(P). \qquad (23)$$

In order to sample from $S^{\dagger}(P)$, we must take into account the fact that in general the adjoint source is not normalized: we will thus sample from an arbitrary distribution $\hat{S}^{\dagger}$, and adjust the statistical weight accordingly, to preserve the unbiasedness. The most common strategy is to set

$$\hat{S}^{\dagger}(P) = \frac{S^{\dagger}(P)}{\int S^{\dagger}(P')dP'}, \qquad (24)$$

which is usually properly normalized, and assign the particle a weight equal to the norm of the source, i.e. $\int S^{\dagger}(P')dP'$. Observe that $\eta_{\psi}$ can be of infinite norm. For example, if the sought response is the total number of collisions occurring in the whole phase space, then $\eta_{\psi} = 1$ everywhere, which results in an infinite norm whenever the total phase space volume is also infinite. In such cases, it is mandatory to sample adjunctons from an arbitrary distribution $\hat{S}^{\dagger}$, and then correct the weights of the source particles by a factor $S^{\dagger}(P)/\hat{S}^{\dagger}(P)$. Similarly as in direct Monte Carlo games, the choice of $\hat{S}^{\dagger}(P)$ will impact the final variance obtained on the sought physical quantity.
### III.C.2. Adjoint tallies

Scoring during an adjoint Monte Carlo game is quite similar to scoring in a direct game. However, there are a few issues to be addressed.

If the goal of the adjoint Monte Carlo game is to estimate the importance, one needs to first define which type of importance is sought. In order to estimate the importance of collided particles, one needs to count particles exiting the adjoint collision kernel (and not the flight kernel, as was the case in the direct simulation), as $\psi^\dagger$ is obtained after applying $C^\dagger$ (see Eq. (9)). Similarly, the emitted importance $\chi^\dagger$ would be estimated after the adjoint flight kernel (see Eq. (10)).

Most often, the sampling of adjoint Monte Carlo games is done in order to estimate the response $R$, in which case the adjoint score will be the direct source density function (source distribution times its norm), namely

$$ \eta_{\psi}^{\dagger}(P) = S(P). \qquad (25) $$

A possible collision-based estimator for Eq. (25) would be to accumulate $S(P)$ at each collision event.

### III.D. Extensive and intensive quantities

Formally, the direct transport equations and the importance equations are both Fredholm integral equations. However, as discussed in the previous section, their solution by Monte Carlo methods involves peculiar differences. In this section we will focus on the origins of this behavior.

In direct Monte Carlo games, the sought quantity is generally expressed in terms of a particle density, that is an intensive quantity. However, the transported quantity is a number of neutrons (i.e. an ‘integrated density’), carried by particles that represent a different number of neutrons depending on their statistical weight. This quantity is extensive: the total amount of neutrons in a given region is the sum of all the neutrons in the sub-regions composing that region: this additivity property ensures that the source can be split into ‘packets’ (the sampled particles) that are subsequently followed during their random trajectories. Additivity also applies to the tallies: the contributions of different particles to a detector are summed in order to infer the sought response. The number of particles is collected into detectors, and the result is divided by the volume of the detector to obtain an average particle density if needed.
In adjoint Monte Carlo games, conversely, the sought quantity is the importance, i.e. the average contribution to a detector for a particle of unit weight in a given phase space point. The importance is intensive and is not additive: the average importance in a given region is not the sum of the importance in all sub-regions, but their weighted average. As a consequence, the importance cannot be split into ‘packets’ carried by particles, and the importance tally cannot be collected by summing single contributions. Nonetheless, in order to solve adjoint problems by Monte Carlo methods, the adjoint equation is transformed into a transport equation for an extensive quantity. In practice, for adjoint transport problems, the quantity carried by the particles is rather an ‘integrated importance’, which is extensive: this is analogous to the case of direct Monte Carlo games. Adjunctons thus carry an amount of integrated importance representing an importance over a given volume of phase space, and have units of

$$ \frac{[\text{Detector response}]}{[\text{neutron}]} \times [\text{Phase space volume}]. \qquad (26) $$

Each adjuncton carries a product of importance × volume that can be distributed in a given portion of the phase space: the bigger the volume, the smaller the importance and vice versa. The integrated importance is collected by detectors, and the result is divided by the volume of the detector to obtain the average importance in that region.

Observe that the integrated importance is contravariant with respect to the units of phase space volume. The impacts of this fact will be detailed in Appendix A.

### III.D.1. Physical meaning of the weight corrections for the adjoint collision

The presence of all the weight corrections in the adjoint collision that are not present in the direct collision stem precisely from the intensive nature of importance.

In the Monte Carlo sampling strategy illustrated in Fig. 1 for the direct transport equation, the resulting branching structure for the choice of nuclides, reactions and post-collision features splits the number of particles undergoing a collision event into different reaction paths. Excluding the multiplicity $\nu_{i,j}$, all the coefficients appearing in Fig. 1 correspond to probability distributions, since the splitting distributes the total number of particles. It is therefore possible to directly sample from these distributions, without any weight correction factor. Once the complete path has been chosen, additivity ensures that all the contributions of the particles arriving at an exit
point $P'$, coming from any reaction path, can be summed without applying weight corrections. The only non-normalized step is the multiplicity, which requires either sampling multiple particles and preserving the average yield (this approach is usually called ‘analog’ Monte Carlo sampling), or applying a weight correction to the emitted particle (non-analog sampling).

For adjoint transport, four extra weight correction factors need to be applied, as shown in Fig. 2, namely $\Sigma_t^\dagger/(N_i\sigma_{i,t}^\dagger)$, $\sigma_{i,t}^\dagger/\sigma_{i,j}^\dagger$, $\sigma_{i,j}/\sigma_{i,t}$, and $N_i\sigma_{i,t}/\Sigma_t$, which have no equivalent in direct transport. Any reaction channel $j$ on nuclide $i$ that led to the emission of a neutron is associated to the same emission importance $\chi^\dagger(P')$: the importance is not split between the different paths, but all paths have the same importance as the original state (since the importance is an intensive quantity), which explains the unit coefficient occurring at this step in the collision kernel (see Fig. 2). In a Monte Carlo approach, we need to split the particle at the emission point into the different possible reaction channels: in order to obtain an unbiased result, we have then to either sample multiple particles, or multiply the weight of the particle.

After choosing a reaction channel, and an ingoing energy $E$ and direction $\Omega$, additional weight correction factors are required before concluding the sampling of the adjoint collision kernel. Let us define $\psi_{i,j}^\dagger(P)$ the importance of a neutron entering a collision on nuclide $i$ and reaction channel $j$ with coordinates $P$. The average number of neutrons actually sampling the pair $\{i, j\}$ is $(N_i\sigma_{i,t}/\Sigma_t) \times (\sigma_{i,j}/\sigma_{i,t})$. The importance of a particle entering a collision at $P$ is a then a weighted average, and the corresponding factor for the branch $\{i, j\}$ is precisely the fraction of the importance carried by the pair $\{i, j\}$:

$$ \psi^\dagger(P) = \sum_{\{i,j\}} \frac{N_i(\vec{r})\sigma_{i,t}(E)}{\Sigma_t(\vec{r}, E)} \times \frac{\sigma_{i,j}(E)}{\sigma_{i,t}(E)} \times \psi_{i,j}^\dagger(P) \quad (27) $$

Because of the underlying weighted average operation, we cannot simply sum all the adjunctons coming from all paths, but we also have to apply the corresponding weight corrections.

## IV. ADJOINT SAMPLING STRATEGIES

In the previous section, we presented a formalism yielding the possibility of sampling adjoint Monte Carlo games for the specific case of an infinite medium, where only the collision kernel needs to be considered, in addition to the source and the detector. For this purpose, we introduced
arbitrary adjoint cross sections $\sigma_{i,j}^{\dagger}$ and arbitrary adjoint distribution laws $f_{i,j}^{\dagger}$. Suitable weight corrections depending on the choice of $\sigma_{i,j}^{\dagger}$ and $f_{i,j}^{\dagger}$ have been set in order to ensure that the resulting game is unbiased. We will now discuss some specific strategies for adjoint Monte Carlo games, depending on the functional form chosen for $\sigma_{i,j}^{\dagger}$ and $f_{i,j}^{\dagger}$.

## IV.A. Sampling quality

A guiding principle of Monte Carlo games is that the effectiveness of sampling strategies can be assessed based on their impact on the calculation time as well as the variance on the sought response, which will be the focus of this section. As mentioned in Section III.B.1, the choice of the adjoint cross sections $\sigma_{i,j}^{\dagger}$ and adjoint distribution laws $f_{i,j}^{\dagger}$ that dictates the sampling distributions does not affect the unbiasedness of the average tallies (provided that the required weight corrections are applied), but has definitely an overall effect on the resulting variance of the sought response.

A naive approach to the sampling of the adjoint collision kernel could consist, for example, in sampling uniformly among the nuclides present in the mixture, and among the associated reaction channels<sup>e</sup>. However, it seems inappropriate to sample as many particles in low- and high-probability reaction channels. Generally speaking, such an approach would result in a broad dispersion of the weight corrections and yield a large variance in the sought tallies, since in most practical cases nuclide concentrations are not evenly distributed in the mixtures, and some reaction channels are similarly more likely than others. The most striking example is provided by the capture reaction: although a neutron ‘emitted’ by a capture event will virtually have the same importance as a neutron emitted by any other reaction at the same location in phase space, it is unwise to sample an adjuncton there, because the event of a neutron emerging from a capture is impossible, by definition of a capture event. In this case, a little anticipation would have shown that the null factor associated to the multiplicity of the capture reaction would kill the adjuncton anyway, so there is no interest in sampling the adjuncton for this reaction channel. A bad choice of the sampling scheme might even lead to an unacceptable amount of variance, making statistical convergence practically impossible.

Based on these considerations, in the following we will determine whether there exists a suitable sampling strategy for the adjoint kernel capable of effectively reducing the variance stemming

<footnote><sup>e</sup>The corresponding definition of adjoint cross sections is $\sigma_{i,j}^{\dagger} = 1/(N_i \times N_{\text{reac},i})$, with $N_{\text{reac},i}$ the number of different reaction channels in nuclide i, leading to $\sigma_{i,t}^{\dagger} = 1/N_i$.</footnote>
from said sampling.

IV.A.1. Insights from zero-variance forward games

Let us assume for the sake of simplicity that the response is associated to a single detector. Several variance-reduction techniques for the sampling of the collision events are widely used in direct Monte Carlo [23], e.g.:

* A particle being captured will usually not be able to contribute to the detector: it has thus zero importance, and should not be sampled; this is known as implicit capture.

* When sampling several reactions that have different yields $\nu_{i,j}$, it would make sense to modify the probability of sampling a reaction proportionally to the yield of each channel, as the combined importance of all the outgoing particles will be (to a first approximation) proportional to the number of outgoing particles: this is the idea at the root of ‘branchless’ collision sampling methods, which force a single particle to be emitted and adjust the weight accordingly.

Similarly, for the sampling of the flights in direct Monte Carlo one should favor particles going in the direction of the detector, as they will then have a higher importance: this has led to the development of strategies aimed at solving deep-penetration problems, like path-stretching or exponential transform [23].

More broadly, it has been shown that there exists a distinct way of sampling the particle trajectories for direct Monte Carlo games (including the source, and the flight and collision events) that yields the lowest possible variance of the sought tally, i.e. zero [24, 23]. Such optimal sampling strategy leading to zero variance depends on the chosen estimator and on specific implementation details. Assuming that the sampling of the collision events is performed using the branchless method, for the collision estimator one needs to sample the source event from the normalized density

$$ \tilde{S}(P) = S(P) \frac{\chi^{\dagger}(P)}{\int S(P) \chi^{\dagger}(P) dP}, \qquad (28) $$

the flight events from the normalized density

$$ \tilde{T}(P \to P') = T(P \to P') \frac{\psi^{\dagger}(P')}{\chi^{\dagger}(P)}, \qquad (29) $$
and the collision events from the normalized density

$$ \tilde{C}(P \to P') = C(P \to P') \frac{\chi^{\dagger}(P')}{\psi^{\dagger}(P)}. \qquad (30) $$

Equations (28) to (30) clearly show that in such games the source, flight and collision kernels must be sampled proportionally to the relative importance of the different outgoing states: thus, more important outgoing states will be favored, whereas less important will be disfavored. The statistical weights of the particles will be correspondingly adjusted in order for the sampling scheme to yield an unbiased result.

Such ‘zero-variance’ Monte Carlo games can be sampled in practice, but their application to practical problems is hindered by the fact that they would require a perfect knowledge of the solution of the adjoint Boltzmann equation, i.e. the importance function, in every point of the phase space. Since its knowledge is tantamount to the solution of the direct transport problem (which is precisely the equation that we are trying to solve in the first place), zero-variance games are not naively used as such, but provide the inspiration for approximate importance-based sampling strategies that would ideally yield zero variance in the limit of a perfectly known solution of the adjoint Boltzmann equation [24]. The most successful application is the use of an approximated adjoint function obtained by deterministic solvers to improve direct sampling, which goes under the name of Consistent Adjoint Driven Importance Sampling (CADIS), and has been shown to significantly reduce the variance of the tallies [25].

### IV.A.2. Implications for adjoint games

The use of the solution of the adjoint Boltzmann equation to conceive zero-variance direct Monte Carlo games suggests that a dual strategy might be adopted for adjoint Monte Carlo games, using the solutions $\psi$ and $\chi$ of the direct Boltzmann equation as ‘adjoint importance function’ to modify the adjoint sampling, via the adjoint cross sections and distribution laws. This stems from the fact that the adjoint of the adjoint of an operator is the original operator: the importance of an adjuncton at point $P$ relative to the source will be the collision density (or emission density) at that point $P$. Note however that, similarly as in the case of the direct Monte Carlo games, obtaining the exact solutions $\psi$ and $\chi$ to obtain zero-variance adjoint games is not possible in practice; nonetheless, we can use this ideal strategy as a guiding principle to build effective low-variance
adjoint games enabling a viable implementation.

This goal is achieved by suitably tailoring the adjoint cross sections $\sigma_{i,j}^{\dagger}$ and the adjoint distribution laws $f_{i,j}^{\dagger}$. We would like the adjoint cross sections at a point P to be proportional to the number of particles leaving the corresponding reaction at that point in a direct transport problem. Similarly, we would like $f_{i,j}^{\dagger}(P' \to P)$ to be proportional to the number of particles undergoing reaction j on nuclide i at P and leaving the collision at $P'$. This would ensure that the sampling of the adjoint collision kernel would yield zero-variance tallies, in the limit case where the adjoint importance function is perfectly known$^f$.

The knowledge of $\psi$ and $\chi$ can be either obtained by using a deterministic solver, or by simulating a direct game and estimating how many (direct) particles go through each reaction channel $\{i, j\}$ for each point in the phase space. The idea of coupled direct-adjoint Monte Carlo sampling was developed by Carter and McCormick [9] and implemented in a code named MCN-A [13]. One can then use this information to set the adjoint cross sections and adjoint distribution laws. This scheme can be considered as a sort of adjoint-CADIS method, leading to reduced variance (or even zero variance in the limit of perfect evaluation of the direct densities). Such calculations will be demonstrated in a separate work. Even when these detailed pieces of information are not available, a general knowledge of the transport problem under consideration (e.g. the concentrations of the nuclides and their cross section values at different energies) might provide sufficient hints to estimate the amount of neutrons emerging from each path, which can be then converted into sampling strategies for the adjoint games. In the next sections, we will examine different approaches that can be used to effectively build adjoint cross sections.

## IV.B. Strategies based on approximations of direct quantities

### IV.B.1. Flux-based sampling

We would like to build the adjoint cross sections $\sigma_{i,j}^{\dagger}(P')$ based on the average number $\chi_{i,j}(P')$ of particles exiting collision from $\{i, j\}$ at $P'$. The quantity $\chi_{i,j}(P')$ can be evaluated by summing the contributions by all possible entering energies and directions that emit particles at $P'$:

$$ \chi_{i,j}(\vec{r}, E', \vec{\Omega}') = \iint \varphi(\vec{r}, E, \vec{\Omega}) N_i(\vec{r}) \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \to E', \vec{\Omega}' \cdot \vec{\Omega}') dE d\vec{\Omega}. \quad (31) $$

<footnote>$^f$Provided that the source and the flight kernel are similarly subject to ideal importance-based sampling</footnote>
As shown in the scheme in Fig. 2, the probability of choosing the adjoint reaction path {$i, j$} is defined as $N_i \sigma_{i,j}^\dagger / \Sigma_t^\dagger$, and we want this probability to be proportional to $\chi_{i,j}$. This yields

$$ \sigma_{i,j}^{\dagger}(\vec{r}, E', \vec{\Omega}') \propto \frac{\chi_{i,j}(\vec{r}, E', \vec{\Omega}')}{N_i(\vec{r})} = \iint \varphi(\vec{r}, E, \vec{\Omega}) \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}') dE d\vec{\Omega}. \quad (32) $$

However, using this formula might prove unpractical, for several reasons. First, as previously mentioned in Section IV.A, perfectly knowing the angular flux $\varphi$ is not practically possible. Second, this formula gives adjoint cross sections that depend not only on the energy $E'$ (as customary for forward microscopic cross sections), but also on the direction $\vec{\Omega}'$ and the position $\vec{r}$. This dependence will make the discretization and storage of adjoint data much more difficult than those of their forward counterparts.

To overcome this issue, let us first assume that the angular flux $\varphi$ is roughly isotropic, and that its dependence on energy $E$ and position $\vec{r}$ can be factorized:

$$ \varphi(\vec{r}, E, \vec{\Omega}) \propto \varphi(\vec{r}, E) = g(E) \tilde{\varphi}(\vec{r}). \quad (33) $$

Under this assumption, we can write:

$$ \sigma_{i,j}^{\dagger}(\vec{r}, E') \propto \tilde{\varphi}(\vec{r}) \int g(E) \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \rightarrow E') dE. \quad (34) $$

We used the notation $f_{i,j}(E \rightarrow E') = \int f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}') d\vec{\Omega}'$. Note that we can exchange $\vec{\Omega}$ and $\vec{\Omega}'$ in the integration and keep no dependence on $\vec{\Omega}'$, since the distribution law $f_{i,j}$ only depends on the scalar product $\vec{\Omega} \cdot \vec{\Omega}'$ and not on both directions independently.

Since the probability of choosing the adjoint reaction path {$i, j$} implies ratios of adjoint cross sections, all evaluated at the same collision site $\vec{r}$, the dependence on $\tilde{\varphi}(\vec{r})$ can be dropped, and we define:

$$ \sigma_{i,j}^{\dagger}(E') = \int g(E) \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \rightarrow E') dE. \quad (35) $$

In order to construct the adjoint distribution law $f_{i,j}^\dagger$, we can follow a similar reasoning. We want $f_{i,j}^\dagger(P \rightarrow P')$ to be proportional to the number of particles going from $P$ to $P'$ via {$i, j$}:

$$ f_{i,j}^{\dagger}(\vec{r}, E' \rightarrow E, \vec{\Omega}', \vec{\Omega}) \propto \varphi(E, \vec{\Omega}', \vec{r}) N_i(\vec{r}) \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}'). \quad (36) $$
This formula has the same problems as Eq. (32). Using the simplifications on the flux, and evaluating the normalization constant$^g$, we have:

$$ 
\begin{aligned}
f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) &= \frac{g(E)N_i(\vec{r})\sigma_{i,j}(E)\nu_{i,j}(E)f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}')}{\iint g(E)N_i(\vec{r})\sigma_{i,j}(E)\nu_{i,j}(E)f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}')dEd\Omega} \\
&= \frac{g(E)\sigma_{i,j}(E)\nu_{i,j}(E)f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}')}{\sigma_{i,j}^{\dagger}(E')}. \qquad (37)
\end{aligned}
$$

The only dependence on the flux left in Eqs. (35) and (37) is through the function $g(E)$. Furthermore, note that the total norm of $g(E)$ is not important: multiplying $g(E)$ by a constant would not impact the sampling. We can then construct the adjoint cross sections by plugging in $g(E)$ the asymptotic energy dependence of the flux expected for the transport problem under consideration. By taking the simplest assumption $g(E) = 1$, i.e. a flat flux along the energy range, from Eqs. (35) and (37) we would obtain

$$ \sigma_{i,j}^{\dagger}(E') = \int \sigma_{i,j}(E) f_{i,j}(E \rightarrow E') \nu_{i,j}(E) dE \qquad (38) $$

and

$$ f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) = \frac{\sigma_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega}' \cdot \vec{\Omega}) \nu_{i,j}(E)}{\sigma_{i,j}^{\dagger}(E')} \qquad (39) $$

The definitions of $\sigma_{i,j}^{\dagger}$ and $f_{i,j}^{\dagger}$ presented in Eqs. (38) and (39) correspond to the ‘first sampling method’ for adjoint Monte Carlo games proposed by Hoogenboom and others [1, 8, 14]. Observe that in this case $\sigma_{i,j}^{\dagger}$ have the same dimensions as $\sigma_{i,j}$.

However, we might want to improve on this method using more accurate estimates concerning the direct particle transport. In a number of problems involving neutron slowing down, the flux generally takes a $1/E$ shape. Plugging in $g(E) = 1/E$ into Eqs. (35) and (37), we would obtain

$$ \sigma_{i,j}^{\dagger}(E') = \int \frac{\sigma_{i,j}(E) f_{i,j}(E \rightarrow E') \nu_{i,j}(E)}{E} dE \qquad (40) $$

and

$$ f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) = \frac{\sigma_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega}' \cdot \vec{\Omega}) \nu_{i,j}(E)}{E \sigma_{i,j}^{\dagger}(E')} \qquad (41) $$

<footnote>$^g$We need to sample directly from properly normalized $f_{i,j}^{\dagger}$.</footnote>
This assumption corresponds to the ‘second sampling method’ presented by Hoogenboom^h and is known to yield a lower variance than Eqs. (38) and (39) in most problems [8, 14]. The adjoint cross sections defined in Eq. (40) do not carry the physical units of cross sections anymore. However, recall that in adjoint Monte Carlo games the adjoint cross sections only occur in the form of ratios of partial cross sections to total cross sections, in order to define the probability of sampling a given nuclide and reaction pair. Thus, the physical units of the adjoint cross sections are of no practical importance.

To further improve the variance of the adjoint games, one could refine the weighting function $g(E)$ to take a shape tailored on the specific transport problem under examination. For example, we might add a Maxwellian component to the $1/E$ asymptotic slowing-down behavior to take into account the thermal distribution of neutrons, or add a Watt fission spectrum to better sample energies in fissile material. One could also produce different sets of adjoint cross sections and distribution laws generated with different weighting functions $g(E)$, and use the adapted data in portions of the problem where the shape of the direct flux fits best with the $g(E)$ used to generate the data. In this work, we focus on infinite spatially homogeneous media, so that this possibility is left as an avenue for future investigations.

### IV.B.2. Sampling based on other direct quantities

In the previous section, we used the direct angular flux $\phi$ to express the emission density (see Eq. (31)) used to define the adjoint cross section $\sigma_{i,j}^\dagger$ (see Eq. (32)) and the adjoint distribution laws $f_{i,j}^\dagger$ (see Eq. (36)). In principle, any other direct quantity can be used to build the adjoint cross sections, instead of the direct flux. For example, using the total reaction rate $\psi = \Sigma_t \phi$ would yield

$$ \sigma_{i,j}^{\dagger}(\vec{r}, E', \vec{\Omega}') \propto \iint \frac{\psi(\vec{r}, E, \vec{\Omega})}{\Sigma_t(\vec{r}, E)} \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}') dE d\vec{\Omega} \quad (42) $$

and

$$ f_{i,j}^{\dagger}(\vec{r}, E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) \propto \frac{\psi(\vec{r}, E, \vec{\Omega})}{\Sigma_t(\vec{r}, E)} N_i(\vec{r}) \sigma_{i,j}(E) \nu_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}'). \quad (43) $$

<footnote n="5">⁵Hoogenboom adds a reference energy $E_0$ as a factor in both definitions, to ensure homogeneity of adjoint cross sections. This has no practical impact on the adjoint Monte Carlo game behavior.</footnote>
Again, to avoid the aforementioned issues with the knowledge of $\varphi$, we suppose that the collision rate is isotropic and can be factorized as

$$ \psi(E, \vec{r}) = g(E) \tilde{\psi}(\vec{r}). \qquad (44) $$

Here we abusively used the same notation $g(E)$ for the energy dependence of the reaction rate $\psi$ as for the flux $\varphi$, since we anticipate that they will be roughly similar for our applications.

Using the same approximations as in the previous section, we obtain

$$ \sigma_{i,j}^{\dagger}(\vec{r}, E') = \int \frac{\sigma_{i,j}(E) f_{i,j}(E \rightarrow E') \nu_{i,j}(E) g(E)}{\Sigma_t(\vec{r}, E)} dE \qquad (45) $$

and

$$ f_{i,j}^{\dagger}(\vec{r}, E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) = \frac{\sigma_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega}' \cdot \vec{\Omega}) \nu_{i,j}(E) g(E)}{\Sigma_t(\vec{r}, E) \sigma_{i,j}^{\dagger}(\vec{r}, E')} \qquad (46) $$

The interest of these definitions is related to the fact that the approximation $\varphi(E) \propto 1/E$ is known to be poor in the presence of nuclides with resonant cross sections: in this case, the flux tends to dip within the resonances (with a $1/\Sigma_t$ behavior), whereas the reaction rate tends to preserve the asymptotic shape $\propto 1/E$ (at least in the limit of narrow resonances) [26]. Replacing then the approximation $\varphi(E) \propto 1/E$ by $\varphi(E)\Sigma_t \propto 1/E$ in the construction of the adjoint cross sections is expected to yield a more faithful energy dependence for the flux and thus lead to a lower variance in the sampling of the collision kernel.

Of course, this resilience comes at a price: the reaction rate sampling uses $\Sigma_t(\vec{r}, E)$ in the formula, which means that we have to create a set of adjoint nuclear data for each nuclide mixture occurring in the problem under examination. Furthermore, the use of $\Sigma_t$ implies that the sampling laws of the different nuclides can interfere with each other, because their cross sections occur in $\Sigma_t$. To overcome this dependence on mixtures and the inter-nuclide interference effects, one might want to generate adjoint nuclear data using the nuclide microscopic reaction rates $\sigma_{i,t}(E)$, instead of the total reaction rate. Following the same procedure, we would then obtain

$$ \sigma_{i,j}^{\dagger}(E') = \int \frac{\sigma_{i,j}(E)}{\sigma_{i,t}(E)} f_{i,j}(E \rightarrow E') \nu_{i,j}(E) dE \qquad (47) $$
and

$$f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) = \frac{\frac{\sigma_{i,j}(E)}{\sigma_{i,t}(E)} f_{i,j}(E \rightarrow E', \vec{\Omega}' \cdot \vec{\Omega}) \nu_{i,j} g(E)}{\sigma_{i,j}^{\dagger}(E')} \quad (48)$$

Finally, we could use the microscopic reaction rates $\sigma_{i,j}(E)$ as a weighting factor, which would lead to:

$$\sigma_{i,j}^{\dagger}(E') = \int f_{i,j}(E \rightarrow E') \nu_{i,j} g(E) dE \quad (49)$$

and

$$f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega}' \cdot \vec{\Omega}) = \frac{f_{i,j}(E \rightarrow E', \vec{\Omega}' \cdot \vec{\Omega}) \nu_{i,j} g(E)}{\sigma_{i,j}^{\dagger}(E')} \quad (50)$$

Formulas Eqs. (49) and (50) are much simpler than in all the previous strategies, insofar as they do not even depend on the direct cross sections. This comes at the cost that these adjoint cross sections would be less effective in reducing the variance of the collision sampling. Note in particular that the adjoint nuclear data will contain no information on the cross sections of the different reactions, and thus will be completely unable to anticipate which reactions are more likely to occur. This issue could be avoided by using different $g$ functions for different reaction paths, chosen properly to convey the information concerning the relative importance of the $\{i, j\}$ pairs.

### IV.C. Other sampling strategies

Using estimates stemming from direct problems as weighting factors for the adjoint nuclear data, although appealing for the aforementioned reasons, is by no means mandatory, and alternative strategies can be conceived: simpler formulas that would be easier to evaluate and to store, and yet approximate well enough the relative importance of the different paths, have been proposed in the literature.

Some authors have proposed for instance to take the adjoint cross sections equal to the direct cross sections [11, 12, 16], namely:

$$\sigma_{i,j}^{\dagger}(E') = \sigma_{i,j}(E'). \quad (51)$$

For the distribution laws, some authors [16, 12] use the flux-based sampling with $1/E$ shape (see Eq. (41)), while others [11] use a simpler version not taking into account the cross section and multiplicity

$$f_{i,j}^{\dagger}(E' \rightarrow E) = \frac{f_{i,j}(E \rightarrow E')}{E \int f_{i,j}(\tilde{E} \rightarrow E') / \tilde{E} d\tilde{E}}, \quad (52)$$
The main advantage of the strategy formulated in Eqs. (51) and (52) is that there is no need to store adjoint cross sections, and the distribution law is very simple (for some reactions, such as isotropic elastic scattering, it can even be determined analytically). The implicit hypothesis behind Eq. (51) is that the average number of particles emitted by a path {$i, j$} is roughly proportional to the cross section of this reaction {$i, j$} at the exiting energy. This could be very efficient for reactions with very regular cross sections, but might fail catastrophically in the case of threshold reactions that can emit neutrons below the threshold.

A representative example for this kind of behavior is the fission of $^{238}$U, whose cross section has an energy threshold at about one MeV (the cross section is orders of magnitude smaller for all energies under this so-called ‘threshold’ energy), but whose spectrum extends below that energy. Thus, an adjuncton entering collision right below the threshold energy should ideally see fission as a likely reaction (as it is a possible reaction path), but using the sampling technique in Eq. (51) the fission reaction would be sampled extremely rarely (as the direct fission cross section is very small under the energy threshold). Then, if this kind of reaction is triggered during the simulation, the weight correction will create a particle with a very large weight, and potentially a ‘shower’ of particles if splitting is used for population control.

Based on the scheme presented in Fig. 2, the total weight correction applied to an adjuncton going through the adjoint collision kernel is

$$ \frac{\sigma_{i,j}(E)}{\Sigma_t(\vec{r}, E)} \times \frac{\nu_{i,j}(E) f_{i,j}(E \rightarrow E', \vec{\Omega} \cdot \vec{\Omega}')}{f_{i,j}^{\dagger}(E' \rightarrow E, \vec{\Omega} \cdot \vec{\Omega}')} \times \frac{\Sigma_t^{\dagger}(\vec{r}, E')}{\sigma_{i,j}^{\dagger}(E')} \quad (53) $$

In the considered example the term $\sigma_{\text{fission}}(E)/\sigma_{\text{fission}}(E')$ is large, and leads to huge weight corrections. It has been suggested that these shower reactions can be handled by re-routing some adjunctons from a more regular reaction (such as elastic scattering) towards the problematic reaction, in order to ensure a minimum sampling of the shower reaction and thus avoid the weight explosion problem [12]. Unfortunately, This sort of manual re-routing undermines the simplicity and elegance of this method.

Some other strategies rely on the original relation (derived in Eq. (32))

$$ \sigma_{i,j}^{\dagger}(\vec{r}, E', \vec{\Omega}') \propto \frac{\chi_{i,j}(\vec{r}, E', \vec{\Omega}')}{N_i(\vec{r})} \quad (54) $$
and evaluate the $\chi_{i,j}(\vec{r}, E', \bar{\Omega}')$ term directly by running a direct Monte Carlo simulation and counting the number of particles exiting each reaction (this counting is based on a coarse mesh of the phase space). These counts are used as adjoint cross sections for the adjoint calculation [9, 13], which avoids the problem of computing adjoint cross sections (Eqs. (35), (38), (40), (45), (47) and (50)) by numerical integration. The accuracy of this approach will of course be limited by the number of particles in the direct simulation and by the mesh refinement. Doing the same counting to get an estimation of $f_{i,j}^\dagger$ is much more difficult (because it involves higher dimensional data), so these authors use ad-hoc laws depending on the reaction type to approximate the different reactions.

# V. NUMERICAL SIMULATION RESULTS

In order to substantiate the theoretical framework proposed in the previous sections, we carried out a series of numerical investigations on a few relevant benchmark configurations. For this purpose, we have developed a Monte Carlo mini-app capable of simulating continuous-energy direct and adjoint neutron transport in an infinite medium. The code takes into account the sampling of the (direct or adjoint) source, collision kernel, and detector. The only phase-space variable considered here will be energy. Nuclear data are not taken from an official library, but are synthetic, i.e. generated with given functional forms according to the specific needs of our numerical tests. Therefore, energy values used in our examples, although expressed in keV, do not carry physical significance.

## V.A. Purely elastic medium

### V.A.1. Collision importance

We begin by considering the evaluation of the collision importance $\psi^\dagger(P)$, i.e. the expected response for a particle having a collision at point $P$, for a purely elastic scattering medium with constant cross section. The adjoint source is a Dirac delta at $E = 1$ keV, and we will estimate the collision importance from 1 keV to 10 keV by adjoint Monte Carlo sampling. This is done by introducing a grid of 500 logarithmically placed bins, spanning the whole energy range, counting the number of particles exiting the adjoint collision kernel. The response collected by a detector (that corresponds to an 'integrated importance', as explained in Section III.D) is then divided by
the bin width to get the average importance over the considered bin. The results will always be
normalized to one emitted adjunction (i.e. the adjoint source normalization constant is one), so the
corresponding importance is related to a detector whose total response function has a unit norm.

The medium considered in this example has only one nuclide with elastic isotropic scattering
of constant cross section $\sigma_e = 1$. The atomic mass is chosen as $A = 12$, which yields $\alpha =
(A - 1)^2/(A + 1)^2 \approx 0.7$: in one collision, a neutron can lose up to 30% of its energy $^1$.

We simulated $10^6$ particles, separated into $10^3$ batches. For the sampling of the collision
events, we use the flux-based method with a constant weighting function $g(E) = 1$, which cor-
responds to Eqs. (38) and (39). A different sampling technique would also lead to an unbiased
estimate of the importance function $\psi^\dagger(P)$, but with a different variance. In addition to Monte
Carlo simulations, a simple deterministic solver has been developed to numerically determine the
importance $\psi^\dagger(P)$, for comparison.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Collision importance $\psi^\dagger(P)$ in a purely elastic medium with $A = 12$, for a point detector (adjoint source) at $E = 1$ keV. We compare Monte Carlo simulation results (with 3 standard deviations confidence intervals) to a numerical solution obtained using a deterministic solver. See text for detailed analysis.</figcaption></figure>

As illustrated in Fig. 3, the collision importance takes the shape of the inverse Placzek
transients, which is consistent with previous findings [16, 19]. The analytical solution provided in

<footnote n="1">¹In an isotropic elastic scattering reaction, the distribution of the outgoing energy of the neutron is uniform between $\alpha E$ and $E$, where $E$ is the incident energy [26].</footnote>
Ref. [16] coincides with the numerical results shown here.

The simulation results obtained here can be physically interpreted as follows. The importance of a particle entering a collision at some energy *E* depends on the likelihood for the particle to attain the detector (corresponding to the adjoint source, i.e. a Dirac delta at *E* = 1 keV). If the neutron is quite far from that energy, it will have to undergo several scattering reactions before reaching the sought energy, and the precise value of its current energy will not matter. This explains why the importance is asymptotically flat far from the source. At energies closer to 1 keV, the detector can be attained in one or two collisions, creating discontinuities in the importance and its derivatives of increasing order. In the curve shown in Fig. 3, the first discontinuity appears at *E*₀/*α*, between energies that are one collision away from detector and those that are two collisions away. At exactly *E* = 1 keV, the neutron is sure to score, leading to a Dirac of importance that is not shown in Fig. 3 (and will not be in subsequent importance graphs either) for readability. Finally, if the neutron is below the detector energy at *E* = 1 keV, its importance is zero, since up-scattering is forbidden and the particle cannot gain energy and contribute to the detector.

### V.A.2. Weight distribution

The importance curve estimated by Monte Carlo sampling does not convey much information regarding how the different sampling strategies work and how efficient they might be. One could, of course, compare the variances and global figures of merit of different strategies in different test cases, as previously done by many authors [11, 19]. This analysis, however, would be strongly case-dependent, and the variance and figures of merit obtained would not be applicable in real-world simulations (keep in mind that all nuclear data used here are synthetic, with arbitrary energy units).

In this work, we will instead focus on how different strategies lead to the sampling of different energies, and on the overall effect of the sampling strategies on the weights of the simulated particles. Studying the particle's weights is particularly relevant in adjoint calculations as weight correction is omnipresent in the sampling procedure. For this purpose, we have developed a series of tools for the statistical analysis of the adjoint Monte Carlo simulations, capable of extracting information concerning the collision events: incident and outgoing energy, weight correction factors applied during collision, etc. This data is then post-processed and visualized in order to better
grasp how adjoint Monte Carlo schemes behave. The behaviors of different sampling methods are
simple enough to remain likely valid in real-world applications, and can help us understand which
features of a simulation might pose convergence problems with a given sampling method.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Recorded weights of adjunctons during elastic collision, displayed as a function of collision energy. The test case is described in Section V.A. The different sub-figures correspond to using different weighting functions $g(E)$ to generate adjoint cross sections and distribution laws for adjoint Monte Carlo sampling. See text for detailed explanation.</figcaption></figure>

In Fig. 4, we display the weights of adjunctons at collision events, as a function of energy, for
a few choices of the function $g(E)$ for the adjoint Monte Carlo sampling. Both weights and energies
are recorded after the adjoint collision (hence the use of $E$ instead of $E'$ to express energy).

Figure 4(a) shows the case of flux-based collision sampling with $g(E) = 1$, which corresponds to Eqs. (38) and (39). Recall that this technique coincides with the first method used by Hoogenboom [8, 14], Irving [10] and Eriksson et al. [1]. The particles that undergo a collision have their weight multiplied by a constant factor, which has been computed in Ref. [19] and reads $\log(\alpha)/(\alpha - 1)$. Therefore, the weight distribution is composed of lines of equal weights, representing the particles that have undergone the same number of collisions after being sampled from the source. This creates a globally increasing trend for the weights, $\propto E$, which can be explained by observing that the importance tends to be asymptotically flat when the size of lethargy bins increase proportionally to $E$.
Since collisions can send adjunctons to different energies depending on the cosine of deflection, adjunctons that have undergone different amounts of collisions can end up at the same energy, effectively creating a variance in the weight of particles seen at a given energy. This effect depends on the 'spread' of the collision[^1]: if a collision event tends to induce more spread (for example for a lighter nuclide), the variance will get worse. In the limit case of $A = 1$, leading to $\alpha = 0$, the weight correction factor will even diverge to infinity, thus preventing the simulation from running [8].

This effect can be compensated for by using a more appropriate sampling. In such a uniform purely scattering medium with a Dirac source, we know that the actual direct flux would take an asymptotic shape of $1/E$. Thus, it would be more appropriate to use $g(E) = 1/E$ for the adjoint collision sampling, according to Eqs. (40) and (41). Recall that this choice corresponds to the second adjoint sampling method proposed by Hoogenboom [8, 14]. Compared to the $g(E) = 1$ method presented above, setting $g(E) = 1/E$ favors lower outgoing energies $E$ at each collision. The weight correction that ensures unbiasedness will correspondingly decrease the weight of particles going to low energies, and conversely. This yields the pattern displayed in Fig. 4(b).

Figure 4(b) shows that all the adjunctons at a given energy have a single weight, effectively ensuring zero variance in the distribution of population weights at a given point in phase space. Thus, using $g(E) = 1/E$ helps the adjoint Monte Carlo code in sampling the collision events in such a way that an adjuncton going to a specific energy will get a correct weight correction, and its weight will align with the weight of all the other adjunctons at that energy. This alignment property is extremely beneficial to the variance of any tallying that might be carried out on the sampled population.

The reason behind this behavior is that the $g(E) = 1/E$ weighting function takes the shape that would have the direct neutron flux in a purely scattering infinite medium with a direct source at $E = \infty$. In this respect, using this shape for $g(E)$ ensures an almost perfect zero-variance sampling of the adjoint collision kernel, as hinted in Section IV.A. The Monte Carlo game does not truly lead to zero variance, in that our adjoint detectors are not at infinite energy, and thus slight deviations are expected, yielding a finite, albeit small, variance[^k]. Nonetheless, we can see

<footnote>[^1]: The adjoint isotropic scattering can send an adjuncton of energy $E'$ into a range $E \in [E', E'/\alpha]$ [19].</footnote>

<footnote>[^k]: This notion of detector at infinity could be properly defined by measuring all particles crossing the upper energy limit $E_{max}$ allowed for the simulation, and setting a response function inversely proportional to the attained</footnote>
that taking $g(E) = 1/E$ has a stabilizing effect on the weights of the population (zero-variance of the weights at a given energy), leading to variance reduction for any detector that one might want to use in this energy range.

Independently of the type of sampling, the increasing trend of weight $\propto E$ is preserved. The fact that particles have weights that span orders of magnitude is not a problem in itself, as long as these weights are at different points in phase space. That is why, to measure the ‘weight distribution quality’ of the population at a point in phase space, we could use the relative variance of the weights at this point (or in a small region around such point). This notion will be clarified in Section V.D.2.

In this example, we did not test other sampling methods: due to the presence of a single nuclide with a single reaction channel, they all coincide.

## V.B. Absorption

We will now treat the case of absorption, in addition to the elastic scattering (a calculation with pure absorption would not be relevant, as it would yield zero importance everywhere). We use the simulation setup described in Section V.A (adjoint source at $E_{min} = 1$ keV, $10^3$ batches of $10^3$ particles each, and 500 energy bins), with an elastic cross section $\sigma_e = 1$ ($A = 12$); additionally, we set an absorption reaction of constant cross section $\sigma_a = 1$. For this configuration, we adopt $g(E) = 1$.

The obtained collisional importance $\psi^\dagger$ and recorded weights during collisions are shown in Fig. 5. The collisional importance $\psi^\dagger$ has inverse Placzek transients like in Fig. 3, but with an exponential decay for increasing $E$. The farther away from the detector (recall that $E_{min} = 1$ keV) a neutron is, the higher the number of collisions needed to reach the detector is. At each collision there is a probability of being absorbed (here, $\sigma_c/(\sigma_c + \sigma_e) = 1/3$): the probability of surviving the absorptions and reaching the detector decreases with increasing ‘distance’ of a neutron from the detector, which explains the damped behavior of the inverse Placzek transients.

Interestingly, during the adjoint simulation, the decay of the evaluated importance is not obtained by absorptions of adjunctons, but by a decay of their weights. As absorption has a zero multiplicity, the adjoint absorption has a zero adjoint cross section (cf Eq. (35)), and is never outgoing energy $E > E_{max}$. Using it in this case would indeed yield zero-variance.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Collisional importance $\psi^\dagger$ (left) and recorded weights of adjunctons during collision, displayed as a function of collision energy (right). Test case with elastic scattering and absorption (described in Section V.B). See text for detailed interpretation.</figcaption></figure>

sampled. The effect of absorption is taken into account by the weight corrections that occur while sampling the other adjoint reactions, in particular in two last weight correction factors (cf Fig. 2):

$$ \frac{N_i(\vec{r})\sigma_{i,t}(E)}{\Sigma_t(\vec{r}, E)} \times \frac{\sigma_{i,j}(E)}{\sigma_{i,t}(E)} = \frac{N_i(\vec{r})\sigma_{i,j}(E)}{\Sigma_t(\vec{r}, E)}. \qquad (55) $$

If one adds absorption to a nuclide, it will increase the $\Sigma_t$ term, and will decrease the weight correction for all the adjoint reaction paths {$i, j$} in the considered mixture.

This behavior is illustrated in Fig. 5, which shows that the weights of the adjunctons decrease at each collision. This effect counteracts the increasing weights observed in the case with no absorption (see Fig. 4(a)). For each value of A, there exists an amount of absorption $\sigma_a = \sigma_e((\alpha - 1)^{-1} \log(\alpha) - 1)$ that ensures that the increasing effect of adjoint scattering and the decreasing effect of absorption cancel out and the weight correction at each collision can be exactly equal to unity.
## V.C. Inelastic scattering

### V.C.1. Collision importance

We consider now a more involved configuration where an isotropic discrete-level inelastic scattering channel is present, in addition to the elastic scattering channel considered in Section V.A. For the inelastic scattering cross section, we take $\sigma_i = \sigma_e = 1$ for $E > \varepsilon$, and 0 below the threshold energy $\varepsilon$. The rest of the simulation (adjoint source, energy ranges, detector ranges) is kept identical to the specifications used in Section V.A. We take again $A = 12$. For illustration, a few examples of collision importance $\psi^\dagger(P)$ are provided in Fig. 6, by varying the mass parameter $A$ and the energy level $Q$ of the inelastic scattering, or equivalently the threshold energy$^1 \varepsilon = \frac{A}{A+1}Q$.

The obtained profiles of the collision importance displayed in Fig. 6 show a superposition between the importance shapes due to elastic scattering and those related to inelastic scattering, including the presence of a discontinuity due to the threshold energy $\varepsilon$.

### V.C.2. Weight distribution

Let us now investigate the impact of the inelastic reaction on the behavior of different collision sampling methods. The inelastic reaction allows for a much larger change in energy than elastic scattering. To emphasize this effect, we have purposefully moved the adjoint source farther away at $E_{min} = 10^{-4}$ keV, widened the considered energy range from $10^{-4}$ keV to $10^4$ keV, and set $\varepsilon = 1$ keV, so that the threshold lies right in the middle of the energy range (in log scale). We, once again, simulate $10^6$ particle histories, separated into $10^3$ batches. As done in Section V.A.2, several quantities of interest are recorded during collision events. The simulation results are displayed in Figs. 7 and 8. Similarly as in Fig. 4, the energy $E$ will denote the outgoing energy of the adjoint collision and $E'$ will denote the ingoing energy of the adjoint collision.

Figure 7 shows the measured particles weights as a function of the energy of the collision. Using the $g(E) = 1$ weighting function (corresponding to Fig. 7(a)), the population behaves similarly to the purely elastic scattering case (cf. Fig. 4(a)) for low energies, where the adjoint inelastic scattering does not play. However, at energies higher than $E = 10^{-3}$ keV, adjunctons start to undergo inelastic scattering and jump directly to energies higher than the threshold energy

<footnote>$^1$For a thorough definition of the inelastic scattering reaction on discrete level, and a derivation of its behavior, see for example Ref. [8].</footnote>
<figure><img src="image_7.png" /><figcaption>(c) $A = 12$ ($\alpha \approx 0.3$), $\varepsilon = 1.0$ keV</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) $A = 12$ ($\alpha \approx 0.3$), $\varepsilon = 0.1$ keV</figcaption></figure>

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

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

Fig. 6. Collision importance $\psi^\dagger$ corresponding to systems including both elastic and discrete level inelastic scattering, for different choices of physical parameters. See text for detailed problem description.
<figure><img src="image_11.png" /><figcaption>Fig. 7. Recorded weights of adjunctons during collision against collision energy, in a medium with both elastic scattering ($n, el$) and discrete level inelastic scattering ($n, n1$). Different reactions are displayed with different colors. Different sub-figures produced using different weighting function $g(E)$ used to generate adjoint cross section and distribution laws. See text for detailed explanation.</figcaption></figure>

$\epsilon = 1$ keV, skipping the entire range in between $10^{-3}$ keV and $\epsilon$. This introduces a very important weight dispersion after the threshold energy (recorded weights span more than four orders of magnitude). At high energies $E \gg \epsilon$, inelastic scattering tends to be equivalent to elastic scattering (at high energies, the energy deposited in the discrete energy level $Q$ is negligible compared to the energy exchanged by collision, and the distribution law of inelastic scattering converges towards the distribution law of elastic scattering), and the population behaves again as in the pure elastic case. However, as the weight distribution has been widened around the threshold energy $\epsilon$, this effect persists: the weights increase with the energy, effectively ‘shifting’ the weight distribution along a ‘weight $\propto E$’ line. Note that this ‘shift’ is multiplicative (between energies $E$ and $E'$ the weight distribution is multiplied by $E'/E$), so it is visible as a translation in log scale.

Using the $g(E) = 1/E$ sampling (corresponding to Fig. 7(b)), a similar behavior is observed, with a less pronounced variance of the weight distribution. At low energies, the weights align on the ‘weight $\propto E$’ line, as in Fig. 4(b). Getting closer to the threshold energy, particles of different weights start to appear at the same energy. This is coincidental with adjunctons emerging from the first inelastic reactions. Contrary to the $g(E) = 1$ case, there is no ‘hole’ in the energy range
sampled by the adjunctons. At high energies, we also see the ‘shift’ of the weight distribution
effect, albeit with a much narrower weight span (around two orders of magnitude).

<figure><img src="image_12.png" /><figcaption>Fig. 8. Entering energy $E'$ and exiting energy $E$ of adjoint collisions, in a medium with both elastic scattering ($n, el$) and discrete level inelastic scattering ($n, n1$). Different reactions are displayed with different colors. Different sub-figures correspond to using different weighting functions $g(E)$ to generate adjoint cross section and distribution laws. See text for detailed explanation.</figcaption></figure>

Let us now focus on the ‘hole’ in the energy range sampled by the $g(E) = 1$ method. To shed light on this behavior, in Fig. 8 we display the incident energy $E'$ and the outgoing energy $E$ of all recorded collision events, for both $g(E) = 1$ and $g(E) = 1/E$. For both sampling methods, once the inelastic scattering is sampled once, the energy of the adjuncton jumps to an outgoing energy $E$ larger than $\varepsilon$. For $g(E) = 1$, the inelastic scattering is sampled at low energies (it is active from the start at $E_{min} = 10^{-4}$ keV), effectively draining all adjunctons before they can reach $E \approx 10^{-2}$ keV. For $g(E) = 1/E$, the inelastic scattering starts to be sampled above $E \approx 10^{-3}$ keV, and some adjunctons can reach $E \approx \varepsilon$ without sampling the inelastic reaction, effectively filling the energy range where no adjunctons were found in the $g(E) = 1$ case.

The sampling of the different adjoint reactions is governed by the adjoint cross sections $\sigma_j^\dagger$ defined in Section III.B.1, which are shown in Fig. 9.

The effect of the weighting function $g(E)$ is apparent in the shape of the adjoint cross sections for the inelastic reaction. As shown in Fig. 9, for both weighting functions the cross section drops
<figure><img src="image_13.png" /><figcaption>Fig. 9. Adjoint cross sections $\sigma_j^\dagger$ for elastic scattering (n, el) and discrete level (with $\varepsilon = 1$ keV) inelastic scattering (n, n1), generated using flux-based sampling with $g(E) = 1$ and $g(E) = 1/E$, corresponding to Eq. (38) and Eq. (40). See text for detailed explanation.</figcaption></figure>

for energies under $10^{-2}$ keV, which corresponds to the quantity $\varepsilon/(A + 1)^2$, as pointed out by Hoogenboom [8]. However, for $g(E) = 1$, in the energy range between $E = \varepsilon/(A+1)^2 \simeq 10^{-2}$ keV and the threshold energy ($\varepsilon = 1$ keV), the adjoint cross sections of elastic and inelastic scattering are roughly equal. Cross sections having roughly the same size lead to adjoint inelastic scattering being sampled at low energies. For $g(E) = 1/E$, the adjoint cross section of elastic scattering follows a $1/E$ shape over the whole energy range, whereas the inelastic scattering cross section falls off for energies under the threshold energy: this leads to a significant difference in cross sections, and less inelastic scattering events sampled at low energies.

Using a $g(E) = 1/E$ sampling, high energies tend to be less favored. This effectively inhibits particles from gaining several orders of magnitude in energy, since the weighting function will be several orders of magnitude smaller there. Since the energy loss (or energy gain in the case of adjoint transport) in inelastic scattering is at least $\varepsilon$, the inelastic reaction will be disfavored compared to elastic scattering: the latter can lead to lower outgoing energies that will be comparatively favored by the $g(E) = 1/E$ weighting. The larger the gap between the energies attained by elastic and inelastic scattering, the larger will also be the gap between the cross sections.

Using $g(E) = 1$, on the other hand, the adjoint sampling has no 'energy preference' and no
reactions will be disfavored. Some paths that would be very unlikely in the corresponding direct game will be sampled by the adjoint game, and this leads to a larger weight dispersion (as shown in Fig. 7). Furthermore, there will be unexplored zones in the phase space, making importance tallying impossible in these regions.

Based on these considerations, the $g(E) = 1/E$ sampling should be favored for better convergence, although it still generates a large weight variance, especially in the zone around the threshold energy. We might therefore wonder whether better sampling strategies exist. As explained in Section IV.A, we might use a direct Monte Carlo calculation to determine $g(E)$ before sampling the adjoint game, which would provide the ‘best’ strategy for the adjoint collision kernel: it would yield zero variance for the adjoint detector corresponding to the direct source.

Using any particular source for the direct calculation might not be a good idea, since a zero-variance calculation for a specific detector might turn out to be detrimental for the variance of other detectors. In order to get a ‘global improvement’ (if any), similarly as in the case with only elastic scattering we seek a zero-variance game for the ‘detector placed at infinite energy’: this will, at least, have a stabilizing effect on the weights in the whole energy range, thus benefiting any potential detector that we might want to use. For this purpose, we run a direct calculation with a source placed at $E = 2 \times 10^4$ keV $> E_{max} = 10^4$ keV. Using a direct source far outside the sought energy range allows us to tally only the asymptotic flux, thus avoiding the Placzek transients. The flux is then tallied and used as weighting function $g(E) = \varphi(E)$ (using a deterministic solver to obtain the flux to use as weighting function is also possible). The simulation results corresponding to this strategy are shown in Fig. 10.

The direct flux $\varphi$ has the expected asymptotic $1/E$ shape, but with a dip and a discontinuity at the energy threshold $\varepsilon$. The discontinuity in the flux can be explained by the discontinuity in total cross section, which jumps from $\Sigma_t = 1$ below the threshold to $\Sigma_t = 2$ above it. The deviation of the computed flux from the $1/E$ approximation is responsible for the population dispersion that occurs in the game using $g(E) = 1/E$.

In Fig. 11 we illustrate the adjoint cross sections and $E'$ to $E$ diagram resulting from the ‘zero-variance at infinity sampling’. The ‘dip’ and the discontinuity in the weighting function that disfavors even more energies above $\varepsilon$ have the effect of further disfavoring the inelastic scattering from low $E'$, and help quench the sampling ‘blind spot’ entirely.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Direct flux to be used as weighting function g(E) for the adjoint calculation (left) and resulting weight distribution as a function of energy during the adjoint calculation (right). See text for detailed analysis.</figcaption></figure>

V.D. Impact of resonances

We address now a configuration where cross sections exhibit resonances, and assess their impact on flux-based sampling (using Eqs. (35) and (37)) and reaction-rate-based sampling (using Eqs. (45) and (46)).

When using flux-based sampling (for example the first method proposed by Hoogenboom [8]), the adjoint cross sections and adjoint collision distribution are:

$$
\sigma_{i,j}^{\dagger}(E') = \int \sigma_{i,j}(E) f_{i,j}(E \rightarrow E') \nu_{i,j}(E) dE. \qquad (56)
$$

$$
f_{i,j}^{\dagger}(E' \rightarrow E) = \frac{\sigma_{i,j}(E) f_{i,j}(E \rightarrow E') \nu_{i,j}(E)}{\sigma_{i,j}^{\dagger}(E')}, \quad (57)
$$

which follow from Eqs. (38) and (39) when angular dependence is neglected.

The formulas for σ†i,j and f†i,j make use of the direct cross sections σi,j. If σi,j contains a
resonance at an energy E*, the distribution law will also contain a peak at that energy (as it is
proportional to σi,j). As the distribution needs to be normalized, if the area under the peak is
large compared to the whole area in the viable energy range, a big resonance will not only lead to
a higher density of adjunctons exiting at resonant energy E*, but also a corresponding decrease of
(a) Adjoint cross sections for elastic (n, el) and dis-
crete level inelastic scattering (n, n1).

(b) Entering energy $E'$ and exiting energy $E$ of recorded adjoint collisions.

<figure><img src="image_15.png" /><figcaption>Fig. 11. Supplementary data regarding the ‘zero-variance at infinity’ run. Adjoint cross sections generated using Eq. (35) with $g(E) = \varphi(E)$. See text for analysis.</figcaption></figure>

all other available energies. The adjoint cross section involves an integral of the cross section: if the
resonance is large enough, this will also increase the value of the adjoint cross section for all energies
$E'$ where the resonant energy $E^*$ is attainable ($f_{i,j}(E^* \rightarrow E') \neq 0$). The resonance will ‘attract’
particles by increasing the adjoint cross section of the resonant nuclide (therefore increasing the
sampling probability of such nuclide) and concentrating the outgoing energy distribution around
the resonant energy $E^*$.

This ‘attractive behavior’ of resonances is due to the fact that the $g(E) = 1$ sampling method
supposes a constant flux: it does not take into account the self-shielding that will inevitably occur
within the resonance. Neglecting self-shielding will lead to an over-sampling of adjunctons in the
resonance, and conversely an under-sampling of the other energies and nuclides.

In order to examine these effects, we will examine the elastic scattering reaction on a nuclide
with mass parameter $\alpha$, having a resonant cross section at energy $E^{*m}$. We will impose an ideal
shape for the resonance: a step function of relative height $h$ (the cross section $\sigma$ is equal to a

<footnote><sup>m</sup>In most physical problems, scattering resonances are usually accompanied by absorption resonances, the latter will not be studied here.</footnote>
constant $\sigma_0$ outside of resonance and equal to $h\sigma_0$ in the resonance)$^n$ and full width $\Delta E$. The energy interval where the adjuncton can be sampled when entering the adjoint collision at an energy $E'$ is $[E', E'/\alpha]$, i.e. the interval of outgoing energies for the adjuncton has a width of $E'(1 - \alpha)/\alpha = \beta E'$, with $\beta = (1 - \alpha)/\alpha$. We will suppose that the resonance width is 'narrow', i.e. its width is small compared to the width of scattering at the resonant energy:

$$ \Delta E \ll \beta E^*, \qquad (58) $$

and $h \gg 1$, which means that the resonant cross section is 'large'. The setup is summarized in Fig. 12.

<figure><img src="image_16.png" /><figcaption>Fig. 12. Model of resonant cross section of the elastic scattering reaction. The cross section is modeled by a step function whose plateau value is $\sigma_0$ and resonant value is $h\sigma_0$. The resonance starts at $E^*$ and has a full width $\Delta E$. The width of the adjoint scattering energy interval around the resonant energy $E^*$ (all energies $E$ accessible from a fixed entering energy $E' \approx E^*$) is $\beta E' \approx \beta E^*$. We assume $\Delta E \ll \beta E^*$ and $h \gg 1$.</figcaption></figure>

In Appendix B, such a setup is studied, and a simplified model amenable to exact results is built. Based on this analysis, we obtain that:

* the quantity $x := h\Delta E/(\beta E^*)$, which is an ersatz of the size of the resonance, taking into account both the height and the width relative to the slowing down process, is key to characterizing the behavior of the resonance;

* whenever the resonance is in the range of outgoing energies, a fraction $x/(1+x)$ of adjunctons is 'attracted' by the resonance (cf Eq. (68));

<footnote>$^n$For all simulations, we will use $\sigma_0 = 1$.</footnote>
• adjunctons that scatter in the resonance see their weights multiplied by a factor $(1+x)/h$ (cf Eq. (71)); this factor is very small compared to unity, as $h \gg 1$ and $\Delta E/(\beta E^*) \ll 1$;

• adjunctons that scatter outside the resonance see their weights multiplied by a factor $1+x$ (cf Eq. (72)); this factor is bigger than unity;

• the resulting weight variance from one interaction with the resonance is $x$ (cf Eq. (73)).

### V.D.1. Weight distribution

Let us now consider a simple example. The resonance will be centered at $E^* = 1$ keV, and we will simulate adjoint Monte Carlo transport (using flux-based sampling) between $E_{min} = 10^{-1}$ keV and $E_{max} = 10^1$ keV, in order to assess the weight distribution of the adjunctons, before and after the resonance. The adjoint source is placed at $E_{min} = 10^{-1}$. The physical parameters are $\alpha = 0.5$ ($\beta = 1$), $\Delta E/E^* = 1/100$ and $h = 100$. For this calculation, we sample $10^3$ batches of $10^3$ particles each.

Since we want to focus on the effects of the resonance, the weight spreading induced by $g(E) = 1$ is to be avoided. We will therefore take $g(E) = 1/E$, which overcomes the spreading of the weight outside the resonance, i.e. where the cross section is constant.

As the resonant reaction is the only reaction in this configuration (and the medium is infinite), the norm of the cross section has no physical impact on the life of neutrons: they can only sample the elastic scattering, independently of the value of the cross section. Thus, the importance (not shown here) is the same as the importance for a constant cross section (cf. the inverse Placzek curve in Fig. 3).

The weight distribution of the collision as a function of energy is displayed in Fig. 13. At energies lower than the resonance energy, the weights display the 'weight $\propto E'$ behavior that we previously observed (see e.g. Fig. 4(b)). However, once the adjunctons get close to the resonance (i.e. when $E^*$ is in the range [$E'$, $\beta E'$] of outgoing energies), some adjunctons get 'attracted' by the resonance, and their weights decrease. This effect is responsible for the vertical structure at $E^* = 1$ keV in Fig. 13. This 'attractive' behavior of the resonance also leads to adjunctons scattering several times in the resonance, their weights decreasing at each collision. Conversely, the weights of the particles avoiding the resonance will increase at each collision. This explains the lines parallel to $E$ on top of the main line in Fig. 13.
<figure><img src="image_17.png" /><figcaption>Fig. 13. Recorded weights of adjunctons during collisions, displayed as a function of collision energy. Infinite medium with isotropic elastic scattering, whose cross section has a single resonance at $E^* = 1$ keV. We used flux-based sampling with $g(E) = 1/E$ (Eqs. (40) and (41)). See the text for a detailed description of the simulated problem and an interpretation.</figcaption></figure>

Given the numerical values chosen for our example, we have $x = 1$: according to our simplified model, approximately half of the population will be ‘attracted’ by the resonance for each interaction with the resonance (cf. Eq. (68)), i.e. whenever the resonance lies in the accessible range of outgoing energies. Our simulation results confirm that the weight deviation factor for the particles that avoid the resonance is $1+x = 2$ (cf Eq. (72)), and the weight deviation factor for the particles entering the resonance is $(1+x)/h = 1/50$ (cf Eq. (71)).

At energies larger than the resonance energy, the particle weights display again a $\propto E$ trend. We recover the ‘shift of the weight distribution’ behavior, previously observed in Fig. 7.

### V.D.2. Variance of the weights

In order to quantify the impact of resonances, we assessed the variance resulting from the dispersion of the weights. Instead of considering a single detector (which would be case-dependent), we considered the relative variance of the weights (rescaled by the square of the average weight) at a given location in the phase space. This proxy indicator measures the ‘relative spreading’ of the weights, which is an ‘intrinsic’ property of the adjunctonº population at a point: it does

<footnote>ºThis technique could also be used with direct particles.</footnote>
not depend on the possible choice of a detector, and it also does not depend on the amount of
simulated particles. The relation between the variance of the weights and the variance of a specific
detector is not straightforward. However, it seems intuitive that increasing the dispersion of weights
will lead to increased variance for detectors that collect the contributions of this population. In
Appendix C, we illustrate the effect of relative variance of the weights on the relative variance of
detectors. Focusing on minimizing the weight variance over the whole range of energies seems to
be a good heuristic approach in view of minimizing the variance of any detector that one might
want to use.

In this paper, the variance of the weights will be tallied in two ways. The simplest is to
create an array of ‘detectors’ spanning the whole energy range: these ‘detectors’ will tally the
variance of the weights of the particles that collide within their respective bins. The tallying bins
need to be small enough, so that the weight variation will be entirely explained by the ‘intrinsic’
weight dispersion of the particles, and not by the size of the detector (measuring different average
weights at different positions in phase space). Conversely, the smaller the bin size, and the less
accurate will be the tallying of the variance, as fewer particles will actually collide in it. Note that
these ‘detectors’ are not collecting physical quantities: using different adjoint cross sections $\sigma_{i,j}^{\dagger}$
and distribution laws $f_{i,j}^{\dagger}$, e.g., would yield the same ‘physical’ information, but different weight
distributions.

An alternative strategy is based on tallying the weights of particles as they ‘cross’ a specific
energy surface. This requires to monitor each adjoint collision event of the adjunctons, which
is then attributed to all the bins in between their incident and outgoing energy. This indicator
converges faster, as the entire population will contribute to ‘detectors’ of this type, but is more
complicated to implement and less generalizable to more complex problems. An issue of this
method is that one cannot compare the weights of particles going to different energies $E$, since the
average estimated weight depends on this energy (with a weight $\propto E$ relation). To overcome this
issue, such a detector rescales the estimated weight by an educated guess$^P$ of the expected weight
at the energy $E$ before tallying the relative variance: this ensures that different energies $E$ will
have the same average ‘relative weight’, and the relative variance measured is consistent. Such
surface detectors will always yield a value of the intrinsic variance of the weights at a position,

<footnote>$^P$Using a weight $\propto E$ guess for the average weight at energy $E$ is usually correct in most problems. Using an incorrect function can only lead to an increased measured weight variance.</footnote>
averaged over all the ingoing and outgoing energies of the adjunctons.

These two tallying methods have been implemented with 100 log-spaced bins for the ‘bin-wise’ tallying and 100 log-spaced energy surfaces for the ‘surface-wise’ tallying, covering the total energy range. Measured relative variances are reported in Fig. 14.

<figure><img src="image_18.png" /><figcaption>Fig. 14. Variance of the recorded weights of adjoint collisions as a function of energy. Variance of weights measured using two ‘detector’ types. Infinite medium with a resonant elastic scattering reaction. A single resonance is located at $E^* = 1$ keV. See text for detailed explanation.</figcaption></figure>

The relative variance of the weights is around zero before the resonance, undergoes a transition around the resonance, and stabilizes around a value of about 11 after the resonance. This global behavior is measured by both variance tallies, and is in agreement with the intuitive explanation that, in this Monte Carlo game, the flux-based sampling ‘introduces’ no variance when the cross section is constant, but spreads the weights at the resonance. The transition phase around the resonance energy is different for the two tallying methods, since they take into account particles in different states. The variance before the resonance is very small, but not null. Numerical errors in the sampling process lead to a relative weight variance around $10^{-8}$, which is the value tallied by the ‘surface-wise’ estimator. The ‘bin-wise’ estimator is sensitive to the error due to the size of the bins. Here the relative size of the bin is $\Delta E_{bin}/E_{bin} \approx 5\%$, leading to different points in phase space being measured by the same ‘detector’: this yields an expected ‘noise’ on the relative
variance<fnref n="9" /> of about $2 \times 10^{-4}$, which is coherent with the tallied value of about $2 \times 10^{-4}$.

Using the aforementioned techniques to tally the relative variance of the weights resulting from the resonance, a parametric study was conducted to check the dependence of the generated weight variance on the resonance parameters. This analysis was done by varying the scattering width $\beta$, the resonance width $\Delta E$, and the resonance relative height $h$. The results are summarized in Table I.

<table><thead><tr><th>β</th><th>ΔE/E*</th><th>h</th><th>x</th><th>V<sub>r</sub></th></tr></thead><tbody><tr><td>0.3</td><td>1/10000</td><td>3000</td><td>1</td><td>5</td></tr><tr><td>0.3</td><td>1/1000</td><td>300</td><td>1</td><td>6</td></tr><tr><td>0.3</td><td>1/100</td><td>30</td><td>1</td><td>5</td></tr><tr><td>0.3</td><td>1/30</td><td>9</td><td>1</td><td>3</td></tr></tbody></table>

a Varying $\Delta E$ and $h$ at fixed $x$

<table><thead><tr><th>β</th><th>ΔE/E*</th><th>h</th><th>x</th><th>V<sub>r</sub></th></tr></thead><tbody><tr><td>0.01</td><td>1/10000</td><td>100</td><td>1</td><td>4</td></tr><tr><td>0.1</td><td>1/1000</td><td>100</td><td>1</td><td>4</td></tr><tr><td>0.2</td><td>1/500</td><td>100</td><td>1</td><td>5</td></tr><tr><td>0.5</td><td>1/200</td><td>100</td><td>1</td><td>7</td></tr><tr><td>1</td><td>1/100</td><td>100</td><td>1</td><td>11</td></tr></tbody></table>

b Varying $\Delta E$ and $\beta$ at fixed $x$

<table><thead><tr><th>β</th><th>ΔE/E*</th><th>h</th><th>x</th><th>V<sub>r</sub></th></tr></thead><tbody><tr><td>0.3</td><td>1/10000</td><td>30</td><td>0.01</td><td>0.03</td></tr><tr><td>0.3</td><td>1/1000</td><td>30</td><td>0.1</td><td>0.3</td></tr><tr><td>0.3</td><td>1/100</td><td>30</td><td>1</td><td>5</td></tr><tr><td>0.3</td><td>1/75</td><td>30</td><td>1.33</td><td>10</td></tr><tr><td>0.3</td><td>1/50</td><td>30</td><td>2</td><td>19</td></tr><tr><td>0.3</td><td>1/10</td><td>30</td><td>10</td><td>~250</td></tr></tbody></table>

c Varying $x$ at fixed $\beta$

TABLE I

Relative variance of the weights of the population tallied after a resonance, depending on different resonance parameters.

The simplified theoretical model developed in Appendix B suggested that the weight variance should depend on $x$. The findings in Table I confirm indeed this strong dependence on $x$, the variance being roughly proportional to $x$ for small values. A rapid explosion of the variance for $x > 1$ is observed: this behavior is not predicted by the model, since it does not take into account the possibility of several collisions occurring around the resonance.

When varying other factors by keeping $x$ fixed, the resulting variance has only mild variations, which is in agreement with the theoretical model. We restricted ourselves to $\beta < 1$, as the theoretical model supposed $\beta \ll 1$ (heavy nuclide hypothesis). Furthermore, the case $\beta > 1$ is not

<footnote n="9">⁹Since we know that weights are proportional to the energy, with $\Delta E_{bin}/E_{bin} \approx 5\%$, the weights will follow a uniform distribution of relative width $\sim 5\%$. The variance of a uniform distribution is $V = (a-b)^2/12$, where $a$ and $b$ are the bounds of the probability density function. This leads to a relative variance of $V_r = (5\%)^2/12 \approx 2 \times 10^{-4}$.</footnote>
realistic, as it would correspond to very light nuclides ($\beta > 1 \Leftrightarrow A < (2 + \sqrt{2})/(2 - \sqrt{2}) \approx 5.8$),
which have few or no resonances.

### V.D.3. Impact of several resonances

The problem posed by resonances when using flux-based sampling can be further exacerbated when several resonances are present and interact to yield a large relative weight variance. This can be illustrated with a simulation involving three well-separated resonances. For this configuration, the simulation domain will range from the adjoint source energy $E_{min} = 10^{-1}$ keV to $E_{max} = 10^3$ keV, and three resonances will be placed at $E_1^* = 0.5$ keV, $E_2^* = 5$ keV and $E_3^* = 50$ keV, respectively. Each resonance will have a relative width of $\Delta E/E^* = 1/100$, and a relative height of $h = 50$. The considered nuclide is such that $\beta = 1$. For all resonances, we thus have $x = 0.5$. We will use $10^4$ batches of $10^4$ particles each: a large number of particles is needed in order to ensure the convergence of the estimators for the variance of the weights.

The weight distribution as a function of energy is displayed in Fig. 15(a). Outside of resonances, the weights of the population grow proportionally to $E$, this behavior is punctually disrupted by the resonances. Between resonances, the adjunctons undergo a sufficient number of collisions to properly mix$^r$ and recovers the ‘outside of resonance’, proportional to $E$ weight growth.

The relative variance of the weights is shown in Fig. 15(b) as a function of the energy of the particles (once again measured by 100 log-spaced ‘surface-based’ variance detectors). The variance of the weights stabilizes at three values after each resonance, which corresponds to the stable phase in between resonances where the distribution of the weights is simply ‘shifted’ along a weight $\propto E$ line. After the first resonance, a relative variance around 2.9 is tallied; after the second resonance, a relative variance around 13 is tallied; finally, after the third resonance, a relative variance around 60 is tallied.

This behavior is due to the fact that each variance contributes approximately the same multiplicative factor to the final relative variance of the weights. In order to better grasp this phenomenon, we will model the three interactions. We are only interested in deviations from the

<footnote>$^r$As the resonances are separated by a tenfold increase in energy, the adjunctons have to undergo at least three scattering reactions between two resonances. This ensures that the sampling of the following resonance can be considered independent of the outcome of the sampling of the previous resonance.</footnote>
<figure><img src="image_19.png" /><figcaption>Fig. 15. Recorded weights and their variance as a function of the outgoing collisional energy $E$. Variance of weights measured using energy surface based ‘detector’ type. Infinite medium with a resonant elastic scattering reaction. Three resonances are located at $E_1^* = 0.5$ keV, $E_2^* = 5$ keV and $E_3^* = 50$ keV. See text for detailed explanation.</figcaption></figure>

initial weight $\propto E$ line, so we define weight deviation factors as factors applied to the ‘weight deviation’ defined as:

$$ \text{weight deviation} = \frac{\text{weight}}{E/E_{min}}. \qquad (59) $$

Let $W_1, W_2$ and $W_3$ be the random variables describing the weight deviation factor applied by the interaction with each resonance, and $W_l$ be the total weight deviation after $l$ resonances. The resonances being identical, we suppose that the $W_k$ random variables follow the same law. Each resonance being separated by a decade of energy, representing at least 3 scattering reactions, we suppose furthermore that the sampling of each resonance is independent of the sampling of the others. As the growth proportional to $E$, seen outside of resonances does not impact the weight deviation, $W_l$ after $l$ resonances is the product of all resonances weight deviations:

$$ W_l = \prod_{k=1}^{l} W_k, \qquad (60) $$
with $W_k$ being independent and identically distributed (iid) variables. Thus, we have<fnref n="8" />

$$ \frac{\mathbb{E}[W_l^2]}{\mathbb{E}[W_l]^2} = \prod_{k=1}^{l} \frac{\mathbb{E}[W_k^2]}{\mathbb{E}[W_k]^2} = \left( \frac{\mathbb{E}[W_1^2]}{\mathbb{E}[W_1]^2} \right)^l, \qquad (61) $$

and in terms of relative variances we finally have:

$$ \mathbb{V}_r[\mathcal{W}_l] = (\mathbb{V}_r[\mathcal{W}_1] + 1)^l - 1. \qquad (62) $$

In our example, one resonance creates a relative variance of about $\mathbb{V}_r[\mathcal{W}_1] = 2.9$, tallied after the first resonance, so we have $\mathbb{V}_r[\mathcal{W}_1] = 2.9$, $\mathbb{V}_r[\mathcal{W}_2] = (2.9 + 1)^2 - 1 = 14.2 \approx 13.5$, and $\mathbb{V}_r[\mathcal{W}_3] = (2.9)^3 - 1 = 58.3 \approx 60$, which supports our simulation findings.

The multiplicative behavior of the relative variance of the weights is problematic when the particles may undergo several variance-inducing events: the variance of the weights will diverge after a few collisions, severely hindering the statistical convergence of the Monte Carlo calculation. A representative example of this behavior would be a slowing-down problem with nuclides having resonant cross sections: all the resonances will interact and broadly spread out the population weights. To alleviate this problem, one could enforce population control (such as roulette and splitting), which counteracts the spreading of the weights. This might at least prevent the multiplicative behavior emphasized above, by ensuring that the population entering each resonance has a controlled weight variance (but might result in longer calculation times if it increases the total population size). This analysis is left for future work.

To summarize our findings, using flux-based sampling with a simple weighting function ($g(E) = 1$ or $g(E) = 1/E$) does not properly take self-shielding into account and potentially leads to a catastrophic behavior. One way of overcoming this issue is to use instead reaction-rate-based sampling, or equivalently use a $1/\Sigma_t$ shape in the $g(E)$ function, which would account for self-shielding. In the case of a single resonant reaction, this cancels the cross section dependence of the adjoint cross sections and of the adjoint distribution laws, and thus prevents the 'attractive' behavior of resonances. By applying these techniques, we would recover the behavior observed in the cases using pure elastic scattering with constant cross sections (presented in Section V.A).

<footnote n="8">⁸If X and Y are two independent random variables, then $\mathbb{E}[(XY)] = \mathbb{E}[X]\mathbb{E}[Y]$, and $\mathbb{E}[(XY)^2] = \mathbb{E}[X^2]\mathbb{E}[Y^2]$.</footnote>
VI. CONCLUSIONS

In this work, we have investigated the features of adjoint Monte Carlo methods for particle
transport. We have revisited various sampling techniques that have been proposed in the litera-
ture to solve the adjoint Boltzmann equation. Focusing specifically on the adjoint collision kernel
for continuous-energy neutron transport problems, we have shown that all the methods proposed
in the literature actually belong to a broader set of sampling strategies. Within this theoretical
framework, we have elucidated the underling assumptions for each technique, which is key to under-
standing their relative strengths and shortcomings. We have proposed a general strategy to tailor
the sampling schemes, that consists in starting from zero-variance schemes and making simplifying
assumptions in order to achieve a reasonable compromise between computational complexity and
resulting variance of the sought response.

Then, in order to further substantiate our findings, we have illustrated the statistical behav-
ior of some of these techniques in simple infinite-medium benchmark problems involving elastic
and inelastic scattering. These tests allowed us to pinpoint situations where a given set of sampling
schemes might show unwanted behavior, leading to excessive variance creation. We have deter-
mined how the variance is related to the evolution of the particle population and the associated
statistical weights. We were thus able to quantify the impact of each step of the simulation on
the increase of the overall variance, depending on the chosen sampling strategy. Special empha-
sis has been given to the case of resonances occurring in the cross sections, which are commonly
met in real-world neutron transport problems, and are an extremely challenging issue for adjoint
sampling.

The results obtained in this paper provide a stepping stone towards the ambitious goal of
establishing a full Monte Carlo sampling strategy for the adjoint Boltzmann equation. For this
purpose, further work is needed in several directions. A first requirement concerns the processing,
storage and sampling of the adjoint cross sections and post-collision distributions stemming from
modern nuclear data libraries, based on the theoretical formulas discussed here. This endeavor,
which will necessarily be reaction-type dependent and involve specific numerical methods, will be
the object of a follow-up paper. Furthermore, an additional avenue for research concerns the inves-
tigation of efficient sampling strategies for spatially heterogeneous media, which demand in turn
the sampling of the adjoint flights. As mentioned earlier, this treatment should be more straight-
forward to implement, since the adjoint flight events can be sampled very similarly to their direct counterparts [14, 16]. Once these two objectives are met, the analysis of how variance-reduction techniques behave in realistic configurations will enable the selection of the most promising methods, in view of their implementation in production Monte Carlo codes.

## A. IMPACT OF PHASE SPACE UNITS ON ADJUNCTONS

In Section III.D we observed that the ‘integrated importance’, i.e. the quantity being carried by adjunctons, is contravariant with respect to the units of phase space volume. If we choose to measure phase space with different units (for example we measure energy in MeV instead of keV), then the weight transported by the particles needs to be divided by the corresponding factor (here, $10^3$). This normalization expresses the fact that the measured volume of the detector would also be divided by the same factor, and the correction ensures that the importance stays the same. This behavior can be interpreted as attributing dimensions to the weight transported by the adjunctons.

This dependence is generally not an issue, since units are usually set before running the simulations: no special weight correction is thus needed during sampling. However, the dependence on units still has an impact that can be illustrated with a simple example. Let us assume that a particle at an energy $E$ makes a transition to an energy $E'$, as illustrated in Fig. 16. Both $E$ and $E'$ are continuous states. The corresponding probability density in the direct case is $f(E \rightarrow E')dE'$, whereas in the adjoint case we have $f(E \rightarrow E')dE$.

<figure><img src="image_20.png" /><figcaption>Fig. 16. Example of a transition from state $E$ to state $E'$.</figcaption></figure>

Observe that the density $f(E \rightarrow E')$ is the same for both cases, whereas the differential element is different. We can actually interpret the term $dE/dE'$ in the adjoint kernel precisely as the correction factor needed to transform the integrated importance evaluated at $E'$ into an integrated importance evaluated at $E$. If we needed to measure $E$ and $E'$ in different units, then the correction factor would be carried by the $dE/dE'$ term.
The effect of the differential element ratio might seem superfluous, since $E$ and $E'$ are generally expressed in the same units. However, we need to keep in mind that one might want to provide the distribution in other units, for reasons due to easier storage and ease of sampling. In the case of elastic scattering, for instance, the density $f(E \to E')$ is usually stored in the form of a deflection angle distribution $g(\mu, E)$ for the deflection cosine $\mu$. Imagine that we draw the deflection cosine $\mu$ of the adjuncton from $g^\dagger(\mu, E')$: the required weight correction factor might be assumed to be $g(\mu, E)/g^\dagger(\mu, E')$, by analogy with the previous example concerning energies. However, the decomposed scheme displayed in Fig. 17 shows that this is not the case.

<figure><img src="image_21.png" /><figcaption>Fig. 17. Example of decomposition of a transition from $E$ to $E'$ using $\mu$ as an intermediary variable.</figcaption></figure>

The change of coordinates from $(E, \mu)$ to $(E', \mu)$, which had no effect for a direct particle, introduces now a correction factor. Thus, if we use the distribution $g^\dagger(\mu, E')d\mu$ to draw the $\mu$ value, the total correction weight becomes

$$ \frac{g(\mu, E)}{g^\dagger(\mu, E')} \times \left. \frac{\partial E}{\partial E'} \right|_{\mu=\text{const}}. \qquad (63) $$

This example shows that we always need to use the distributions in the laboratory frame to apply weight corrections, or otherwise take into account the proper correction factors.

## B. MODELING THE BEHAVIOR OF ADJUNCTONS CROSSING A RESONANCE

Consider the case of a single elastic scattering reaction with a resonance at energy $E^*$, as described in Fig. 12. We will model the adjuncton crossing a resonance as a binary event: the adjuncton is either sampled in the resonance, or avoids the resonance (its energy is sampled outside resonant energies). We will derive the probability of each event, and its corresponding weight correction, in view of obtaining an estimate of the variance of the weights generated by this
interaction with the resonance.

First, let us derive a sampling strategy that induces no weight variance. For example a ‘zero-variance for the detector at infinite energy’ strategy as introduced in Section V.A.2, ensures all the weights remain on a line with no weight variance.

Since the resonant reaction is the only reaction in this configuration (and the medium is infinite), the norm of the cross section has no physical impact on the life of neutrons (they can only sample the elastic scattering channel, independently of the value of the cross section). This means that the actual asymptotic reaction rate that one would obtain in a direct calculation is the same as if no resonance was present, i.e. a $1/E$ shape. To simplify the analysis, let us suppose that the resonant nucleus has a large mass, leading to $\beta \ll 1$. For $E$ between $E'$ and $(1 + \beta)E'$ (the range of accessible energies via one adjoint scattering event), we can then approximate

$$ \psi(E) \approx \text{constant}, \qquad (64) $$

and

$$ f(E \to E') = \frac{1}{(1 - \alpha)E} \approx \text{constant}. \qquad (65) $$

Sampling proportionally to the importance (in our case the importance of adjunctons is the asymptotic reaction rate) leads to ‘zero variance for a detector at infinite energy’ (cf. Section IV.A.1 and Eq. (30)). Sampling the outgoing energies uniformly in $[E', \beta E']$ is then the ‘zero variance at infinity’ solution, leading to no weight dispersion. This corresponds to a probability

$$ p_{zv} = \frac{\Delta E}{\beta E^*} \qquad (66) $$

of sampling the outgoing energy in the resonance, and

$$ q_{zv} = 1 - p_{zv} \approx 1 \qquad (67) $$

outside the resonance. Here we used the ‘narrow resonance’ hypothesis stated in Eq. (58), leading to $p_{zv} \ll 1$.

Now, let us derive these probabilities using the standard flux-based sampling (Eq. (57)). We
will sample E within the resonance with probability

$$p = \frac{\int_{E^*}^{E^*+\Delta E} \sigma(E) f(E \rightarrow E') dE}{\int_{E'}^{\beta E'} \sigma(E) f(E \rightarrow E') dE} = \frac{hp_{zv}}{(1-p_{zv}) + hp_{zv}} \approx \frac{x}{1+x}, \qquad (68)$$

where we have introduced

$$x := \frac{h\Delta E}{\beta E^*}, \qquad (69)$$

and we will sample E outside the resonance with probability

$$q = 1 - p \approx \frac{1}{1 + x}. \qquad (70)$$

The deviation from the zero-variance probabilities leads to an inversely proportional deviation in weights. Compared to the zero-variance weight correction, the weights of the particles sampled in the resonance will be multiplied by a factor

$$a = \frac{p_{zv}}{p} = \frac{1+x}{h}, \qquad (71)$$

and the weights of the particles sampled outside the resonance will be multiplied by a factor

$$b = \frac{q_{zv}}{q} = 1 + x, \qquad (72)$$

which are required to ensure that the sampling is unbiased.

The weight deviation being described by a binary variable, we can evaluate its variance using

$$V[\text{weights after resonance}] = (b-a)^2 pq = \left(1+x - \frac{1+x}{h}\right)^2 \frac{x}{(1+x)^2} \approx x, \qquad (73)$$

where we used $h \gg 1$. According to this model, the resulting variance only depends on $x$, which is an ersatz of the size of the resonance, taking into account both the height and the width relative to the slowing down process.

The formula for the variance in Eq. (73) corresponds to a single interaction within the resonance. In practice, the adjunctons might undergo more than one interaction (as observed experimentally in Section V.D): this would lead to weights being multiplied several times, and to
an even larger observed variance. As we supposed $A \gg 1$ to obtain $\beta \ll 1$, we can expect small discrepancies when studying resonant light nuclides (which anyway in most cases have fewer or no resonances).

### C. FROM RELATIVE VARIANCE OF THE WEIGHTS TO DETECTOR VARIANCE

In Section V, the relative variance of the adjunction weights was used as a proxy to characterize the statistical convergence. To illustrate how the variance of the statistical weights can influence the convergence of physical detectors, tests with different trial detectors were run with different amounts of relative weight variance, and the obtained relative variance was compared to that of the detectors.

For this purpose, we consider a configuration involving purely elastic media (isotropic elastic scattering) with $A = 12$, simulated in the range $E_{min} = 0.1$ keV to $E_{max} = 10$ keV. The adjoint source is a point source at $E_{min} = 0.1$ keV. We simulate $10^6$ histories separated into $10^3$ batches. The weight variance is generated by a single resonance in the cross section at $E^* = 1$ keV, and the shape of the resonance is set as described in Fig. 12, with $h = 30$ and variable $\Delta E$ (to obtain different values of weight variance). To ensure that the variance of the weights is only generated by the resonance, we use the flux-based sampling with $g(E) = 1/E$, corresponding to Eqs. (40) and (41).

The obtained relative weight variance is measured by a 'surface-wise' variance tallying (as described in Section V.D.2) at $E = 6$ keV. Several detectors are placed in the energy range to test the impact of the relative weight variance on their statistical convergence:

* A reference detector is placed before the resonance (and should therefore not be impacted by the generated weight variance), counting collisions in the energy interval $[0.3 \text{ keV}, 0.31 \text{ keV}]$.

* Detector #1 is a collision-based detector in the interval $[3 \text{ keV}, 3.1 \text{ keV}]$.

* Detector #2 is a collision-based detector in the interval $[6 \text{ keV}, 6.2 \text{ keV}]$. This detector has the same relative size compared to Detector #1 and no weight variance is added in between so it should yield the same convergence behavior.
* Detector #3 is a collision-based detector in the interval [8 keV, 8.01 keV]. The interval is narrower than the previous detectors, so we expect a larger variance.

* Detector #4 is a special 'detector at infinity', defined as in Section V.A.2. This detector counts all adjunctons whose collision energy crosses the upper energy boundary of the simulation $E_{max} = 10$ keV. Before tallying, the particle weight is divided by its outgoing energy $E$. This corresponds to using the detector response function:

$$h(E', E) = \begin{cases} 1/E & \text{if } E' \le E_{max} < E \\ 0 & \text{otherwise} \end{cases}$$

<table><caption>TABLE II</caption><thead><tr><th>ΔE/E*</th><th>V<sub>r</sub>[weights]</th><th>Ref. Det.</th><th>Det. #1</th><th>Det. #2</th><th>Det. #3</th><th>Det. #4</th></tr></thead><tbody><tr><td>0</td><td>1.6 × 10<sup>-8</sup></td><td>4.7</td><td>4.3</td><td>3.9</td><td>140</td><td>1.6 × 10<sup>-8</sup></td></tr><tr><td>1/200</td><td>1.2</td><td>4.6</td><td>12</td><td>11</td><td>260</td><td>1.2</td></tr><tr><td>1/100</td><td>3.4</td><td>4.4</td><td>23</td><td>23</td><td>530</td><td>3.4</td></tr><tr><td>1/50</td><td>13</td><td>4.6</td><td>81</td><td>75</td><td>1600</td><td>13</td></tr><tr><td>1/30</td><td>50</td><td>4.5</td><td>240</td><td>190</td><td>3000</td><td>51</td></tr></tbody></table>

Relative variance of trial detectors depending on the relative variance of the weights of the adjunctons. $ΔE/E^*$ represents the relative width of the resonance used to generate the weight dispersion. $V_r[weights]$ is the relative variance of the weights measured after the resonance. The five remaining columns are the unitary relative variance ($V_{ur}$), that is the relative variance of the trial detectors multiplied by the total number of particles simulated (giving a theoretical variance for one particle). See text for detailed description of the simulation and for interpretation.

In Table II, the relative variances of all the trial detectors are reported. In the first case, where no resonance is included, we can see that no weight variance is created (the tally ~ $10^{-8}$ is due to discretization errors in the preparation of the adjoint cross sections and distribution laws). As expected, detectors #1 and 2 yield similar results, detector #3 has a larger variance, and detector #4 has a very low variance, since it is built to ideally have a zero variance.

When the resonance size increases, the generated weight variance correspondingly increases. As expected, the reference detector tallying adjunctons before the resonance is not affected. All other detectors see their relative variance increase roughly proportionally to $(1+V_r[weight])$, except for detector #4, which is built for zero variance.
REFERENCES

[1] B. ERIKSSON, C. JOHANSSON, M. LEIMDORFER, and M. H. KALOS, "Monte Carlo Integration of the Adjoint Neutron Transport Equation," *Nuclear Science and Engineering*, **37**, 3, 410 (1969); 10.13182/NSE69-A19116., URL https://doi.org/10.13182/NSE69-A19116.

[2] M. H. KALOS, "Monte Carlo Integration of the Adjoint Gamma-Ray Transport Equation," *Nuclear Science and Engineering*, **33**, 3, 284 (1968); 10.13182/NSE68-A19235., URL https://doi.org/10.13182/NSE68-A19235.

[3] L. B. LEVITT and J. SPANIER, "A New Non-Multigroup Adjoint Monte Carlo Technique," *Nuclear Science and Engineering*, **37**, 2, 278 (1969); 10.13182/NSE69-A20688., URL https://doi.org/10.13182/NSE69-A20688.

[4] C. W. MAYNARD, "An Application of the Reciprocity Theorem to the Acceleration of Monte Carlo Calculations," *Nuclear Science and Engineering*, **10**, 2, 97 (1961); 10.13182/NSE61-A25945., URL https://doi.org/10.13182/NSE61-A25945.

[5] A. M. WEINBERG and E. P. WIGNER, *The physical theory of neutron chain reactors*, University of Chicago Press (1959)URL http://archive.org/details/physicaltheoryof0000wein.

[6] J. LEWINS, *Importance, the Adjoint Function: The Physical Basis of the Variational and Perturbation Theory in Transport and Diffusion Problems*, Elsevier Science & Technology (1965).

[7] A. DE MATTEIS, "Phenomenological interpretation of the adjoint neutron transport equation," *Meccanica*, **9**, 3, 162 (1974); 10.1007/BF02317707., URL https://doi.org/10.1007/BF02317707.

[8] J. E. HOOGENBOOM, "Adjoint Monte Carlo methods in neutron transport calculations," PhD Thesis, Delft University Press (1977)URL https://www.researchgate.net/profile/Jan-Hoogenboom/publication/34734817_Adjoint_Monte_Carlo_methods_in_neutron_transportCalculations/links/02bfe50f6e1f50afac000000/Adjoint-Monte-Carlo-methods-in-neutron-transport-calculations.pdf.
[9] L. L. CARTER and N. J. MCCORMICK, “Coupled Sampling with the Monte Carlo Method in Neutron Transport Calculations,” *Nuclear Science and Engineering*, **39**, 3, 296 (1970); 10.13182/NSE70-A19991., URL https://doi.org/10.13182/NSE70-A19991.

[10] D. C. IRVING, “The adjoint Boltzmann equation and its simulation by Monte Carlo,” *Nuclear Engineering and Design*, **15**, 273 (1971); 10.1016/0029-5493(71)90069-0., URL https://www.sciencedirect.com/science/article/pii/0029549371900690.

[11] A. D. MATTEIS and R. SIMONINI, “A New Monte Carlo Approach to the Adjoint Boltzmann Equation,” *Nuclear Science and Engineering*, **65**, 1, 93 (1978); 10.13182/NSE78-A27129., URL https://doi.org/10.13182/NSE78-A27129.

[12] R. J. BRISSENDEN, “Continuous energy adjoint Monte Carlo,” *Progress in Nuclear Energy*, **24**, 1, 129 (1990); 10.1016/0149-1970(90)90030-9., URL https://www.sciencedirect.com/science/article/pii/0149197090900309.

[13] L. CARTER, “MCNA: a computer program to solve the adjoint neutron transport equation by coupled sampling with the Monte Carlo method.” LA-4488, 4090972, Los Alamos National Lab.(LANL) (1970)URL http://www.osti.gov/servlets/purl/4090972/.

[14] J. E. HOOGENBOOM, “A Practical Adjoint Monte Carlo Technique for Fixed-Source and Eigenfunction Neutron Transport Problems,” *Nuclear Science and Engineering*, **79**, 4, 357 (1981); 10.13182/NSE81-A21387., URL https://doi.org/10.13182/NSE81-A21387.

[15] J. E. HOOGENBOOM, “Continuous energy adjoint Monte Carlo in MCNP with minor code extensions,” *M and C'99-Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Applications*, 1912-1921, Madrid (1999)URL https://inis.iaea.org/records/4qp6h-z4m25.

[16] P. SARACCO, S. DULLA, and P. RAVETTO, “The adjoint neutron transport equation and the statistical approach for its solution,” *The European Physical Journal Plus*, **131**, 11, 412 (2016); 10.1140/epjp/i2016-16412-0., URL http://link.springer.com/10.1140/epjp/i2016-16412-0.
[17] J. E. HOOGENBOOM, "Methodology of Continuous-Energy Adjoint Monte Carlo for Neutron, Photon, and Coupled Neutron-Photon Transport," *Nuclear Science and Engineering*, **143**, 2, 99 (2003); 10.13182/NSE03-A2322., URL https://doi.org/10.13182/NSE03-A2322.

[18] C. M. DIOP, O. PETIT, C. JOUANNE, and M. COSTE-DELCLAUX, "Adjoint Monte Carlo neutron transport using cross-section probability table representation," *Annals of Nuclear Energy*, **37**, 9, 1186 (2010); 10.1016/j.anucene.2010.04.018., URL https://www.sciencedirect.com/science/article/pii/S0306454910001544.

[19] V. VITALI, S. DULLA, P. RAVETTO, and A. ZOIA, "Comparison of Monte Carlo methods for adjoint neutron transport," *The European Physical Journal Plus*, **133**, 8, 317 (2018); 10.1140/epjp/i2018-12132-9., URL https://doi.org/10.1140/epjp/i2018-12132-9.

[20] S. OCCHIPINTI, "Stochastic methods for the solution of the adjoint transport equation," Master's Thesis, Politecnico di Torino (2025) URL https://webthesis.biblio.polito.it/secure/34956/1/tesi.pdf.

[21] M. GRIMSTONE, "Extension of the MCBEND Monte Carlo code to perform adjoint calculations using point energy data," *ANS Radiation Protection and Shielding Division Topical Conference*, 143-150, Nashville (1998) URL https://answerssoftwareservice.com/resource/pdfs/log033.pdf.

[22] T. SHUTTLEWORTH, M. GRIMSTONE, and S. CHUCAS, "Application of Acceleration Techniques in MCBEND," *Journal of Nuclear Science and Technology*, **37**, sup1, 406 (2000); 10.1080/00223131.2000.10874916., URL https://doi.org/10.1080/00223131.2000.10874916.

[23] I. LUX and L. KOBLINGER, *Monte Carlo Particle Transport Methods: Neutron and Photon Calculations*, CRC Press (2018) URL https://www.taylorfrancis.com/books/9781351083287.

[24] J. E. HOOGENBOOM, "Zero-Variance Monte Carlo Schemes Revisited," *Nuclear Science and Engineering*, **160**, 1, 1 (2008); 10.13182/NSE160-01., URL https://doi.org/10.13182/NSE160-01.
[25] J. C. WAGNER and A. HAGHIGHAT, "Automated Variance Reduction of Monte Carlo Shielding Calculations Using the Discrete Ordinates Adjoint Function," *Nuclear Science and Engineering*, **128**, 2, 186 (1998); 10.13182/NSE98-2., URL https://www.tandfonline.com/doi/full/10.13182/NSE98-2.

[26] G. I. BELL and S. GLASSTONE, "Nuclear Reactor Theory," TID-25606, United States Atomic Energy Commission (1970) URL https://www.osti.gov/biblio/4074688.