Sensitivity Analysis in Coupled Monte Carlo Radiation Transport Simulations

Abstract

Sensitivity analysis methods have found extensive use in nuclear criticality safety applications for understanding the impact of uncertain nuclear data on eigenvalue estimates. Significant uncertainty exists in nuclear data and nuclear physics models for photon and electron transport applications, and the goal of this work is to explore whether recently developed adjoint-based, first-order generalized perturbation theory reaction rate sensitivity methods can be extended to coupled Monte Carlo radiation transport simulations. This paper presents a rigorous theoretical derivation for this extended sensitivity analysis method, which is then implemented in a one-dimensional test Monte Carlo code. The adjoint-based sensitivity coefficients are found to agree well with reference direct perturbation and deterministic SENSMG sensitivity coefficients for a simple test problem.

Keywords:

• Monte Carlo
• CLUTCH
• sensitivity analysis

I. INTRODUCTION

Sensitivity analysis methods have found extensive use in nuclear criticality safety applications, for example, by allowing analysts to propagate the impact of uncertainty in evaluated neutron data, by determining the degree of similarity between benchmark experiments and a target validation application, and by using the results of integral benchmark measurements to calibrate neutron cross-section data. Significant uncertainty exists in nuclear data and nuclear physics models for photon and electron transport, and the goal of this work is to explore the potential for sensitivity analysis in coupled radiation transport applications. While previous work developed sensitivity analysis methods in deterministic neutron-gamma radiation transport simulations¹ and recent work developed generalized perturbation theory (GPT) reaction rate sensitivity methods in Monte Carlo neutron transport simulations,² this work explores whether GPT sensitivity analysis methods can be extended to coupled Monte Carlo radiation transport simulations. Although the authors have an interest in ultimately exploring sensitivity and uncertainty analysis for photon and electron transport applications, the work presented here is implemented in a particle-agnostic, one-dimensional (1-D) Monte Carlo code and can therefore be applied to any coupled transport application.

As shown in Eq. (1)(1) $S_{R,{\mathrm{\Sigma }}_x}=\frac{\delta R/R}{\delta {\mathrm{\Sigma }}_x/{\mathrm{\Sigma }}_x}.$(1) , sensitivity coefficients predict the impact of changes to or uncertainties in the nuclear data parameter ${\mathrm{\Sigma }}_x$ on some integral response of interest $R$:

(1) $S_{R,{\mathrm{\Sigma }}_x}=\frac{\delta R/R}{\delta {\mathrm{\Sigma }}_x/{\mathrm{\Sigma }}_x}.$(1)

Sensitivity coefficients for eigenvalue responses (i.e., $R=k_{eff}$) have historically been computed using first-order perturbation theory, which relies on determining some representation of the adjoint flux, i.e., the particle’s “importance.” The derivation of first-order perturbation theory is well-documented by various sources, and it has recently been applied to compute sensitivity coefficients in continuous-energy, eigenvalue simulations.^3,4

Recent work by Perfetti and Rearden^2,5 developed a methodology for estimating sensitivity coefficients for responses that are ratios of reaction rates using continuous-energy Monte Carlo neutron transport methods; i.e.,

(2) $R=\frac{({\mathrm{\Sigma }}_1\varphi )}{({\mathrm{\Sigma }}_2\varphi )}.$(2)

Algorithms for computing these GPT sensitivity estimates use a combination of the Contributon-Linked eigenvalue sensitivity/Uncertainty estimation via Track length importance CHaracterization (CLUTCH) methodology³ and a modified version of the Iterated Fission Probability (IFP) method,⁴ which is needed only for sensitivity analysis of eigenvalue problems.

This work has demonstrated proof of principle for extending these GPT sensitivity analysis methods to coupled Monte Carlo radiation transport simulations. Theoretical developments are presented that provide a methodology for extending these methods in these coupled simulations, and a simple 1-D Monte Carlo test code was developed to achieve proof of principle for the new methodology.

II. CONTRIBUTON THEORY FOR COUPLED RADIATION TRANSPORT SIMULATIONS

We will now develop an expression for estimating the impact of perturbations to transport operators on linear response functionals in coupled radiation transport simulations. This derivation will assume photon-electron transport but is in fact applicable to any coupled radiation transport simulation. Consider, for example, a photon flux $\mathrm{\Phi }$ and an electron flux $\psi$. A response $R_{tot}$ is a combination of arbitrary responses generated by the photon and electron fluxes integrated over some phase-space:

(3) $R_{tot}=⟨R\mathrm{\Phi }⟩+⟨r\psi ⟩.$(3)

These responses in Eq. (3)(3) $R_{tot}=⟨R\mathrm{\Phi }⟩+⟨r\psi ⟩.$(3) are arbitrary but are here considered to be linear functionals, i.e., the product of flux and some cross section ${\mathrm{\Sigma }}_R$:

(4a) $R_{pho.}=⟨R\mathrm{\Phi }⟩={\mathrm{\Sigma }}_{R,pho.}\mathrm{\Phi }$(4a)

and

(4b) $R_{ele.}=⟨r\psi ⟩={\mathrm{\Sigma }}_{R,ele.}\psi .$(4b)

By introducing a perturbation to all system properties (e.g., $\mathrm{\Phi }{ }^{\mathrm{\prime }}$ = $\mathrm{\Phi }$ + $\delta \mathrm{\Phi }$, etc.), ignoring higher-order terms, and subtracting Eq. (3)(3) $R_{tot}=⟨R\mathrm{\Phi }⟩+⟨r\psi ⟩.$(3) , the change in the total response that occurs in response to perturbations or uncertainties in system parameters is described by

(5) $\delta R_{tot}=⟨\delta R\mathrm{\Phi }⟩+⟨\delta r\psi ⟩+⟨R\delta \mathrm{\Phi }⟩+⟨r\delta \psi ⟩.$(5)

The $⟨\delta R\mathrm{\Phi }⟩$ and $⟨\delta r\psi ⟩$ terms in Eq. (5)(5) $\delta R_{tot}=⟨\delta R\mathrm{\Phi }⟩+⟨\delta r\psi ⟩+⟨R\delta \mathrm{\Phi }⟩+⟨r\delta \psi ⟩.$(5) represent the “direct effect” of data perturbations on the response, in other words, the impact of perturbing the nuclear data that is directly contained in the responses of interest from Eq. (4)(4a) $R_{pho.}=⟨R\mathrm{\Phi }⟩={\mathrm{\Sigma }}_{R,pho.}\mathrm{\Phi }$(4a) . The terms containing the flux perturbations in Eq. (5)(5) $\delta R_{tot}=⟨\delta R\mathrm{\Phi }⟩+⟨\delta r\psi ⟩+⟨R\delta \mathrm{\Phi }⟩+⟨r\delta \psi ⟩.$(5) represent the “indirect effect” that the data perturbations will have on the response by changing the flux encountered by the response.² In Eq. (6)(6) $\frac{\delta R_{tot}}{R_{tot}}=\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }{\mathrm{\Sigma }}_X}\delta {\mathrm{\Sigma }}_X+\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\mathrm{\Phi }}\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }{\mathrm{\Sigma }}_X}\delta {\mathrm{\Sigma }}_X+\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\psi }\frac{\mathrm{\partial }\psi }{\mathrm{\partial }{\mathrm{\Sigma }}_X}\delta {\mathrm{\Sigma }}_X.$(6) , the direct effect terms from Eq. (5)(5) $\delta R_{tot}=⟨\delta R\mathrm{\Phi }⟩+⟨\delta r\psi ⟩+⟨R\delta \mathrm{\Phi }⟩+⟨r\delta \psi ⟩.$(5) are combined into a single term, and all perturbations are recast as derivatives with respect to an arbitrary nuclear data parameter ${\mathrm{\Sigma }}_X$ and are multiplied by a perturbation $\delta {\mathrm{\Sigma }}_X$ to ${\mathrm{\Sigma }}_X$:

(6) $\frac{\delta R_{tot}}{R_{tot}}=\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }{\mathrm{\Sigma }}_X}\delta {\mathrm{\Sigma }}_X+\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\mathrm{\Phi }}\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }{\mathrm{\Sigma }}_X}\delta {\mathrm{\Sigma }}_X+\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\psi }\frac{\mathrm{\partial }\psi }{\mathrm{\partial }{\mathrm{\Sigma }}_X}\delta {\mathrm{\Sigma }}_X.$(6)

The direct effect sensitivity terms are generally simple to calculate; therefore, this work will focus on developing expressions to estimate the indirect effect sensitivities.

The transport equations governing the balance of photon and electron fluxes in a hypothetical system are given as follows:

Photons:

(7a) $T\mathrm{\Phi }=Q_{total}=Q+P\psi ;$(7a)

Electrons:

(7b) $t\psi =q_{total}=q+p\mathrm{\Phi }.$(7b)

This system assumes vacuum boundary conditions. These equations assume that photons create electrons via the $p\mathrm{\Phi }$ operator, that electrons create photons via the $P\psi$ operator, and that the system is subject to an external source of photons and electrons ($Q$ and $q$, respectively). In these and future equations, uppercase letters will generally denote photon transport equation terms, and lowercase letters will denote electron transport equation terms.

Adjoint transport equations will now be defined and manipulated to calculate the indirect effect of data perturbations from Eqs. (4)(4a) $R_{pho.}=⟨R\mathrm{\Phi }⟩={\mathrm{\Sigma }}_{R,pho.}\mathrm{\Phi }$(4a) and (5)(5) $\delta R_{tot}=⟨\delta R\mathrm{\Phi }⟩+⟨\delta r\psi ⟩+⟨R\delta \mathrm{\Phi }⟩+⟨r\delta \psi ⟩.$(5) . The adjoint transport equations for the coupled photon and electron fluxes can be shown to be

(8a) $T^{†}{\mathrm{\Phi }}^{†}=S^{†}+p^{†}{\psi }^{†}$(8a)

and

(8b) $t^{†}{\psi }^{†}=s^{†}+P^{†}{\mathrm{\Phi }}^{†}.$(8b)

The adjoint source terms $S^{†}$ and $s^{†}$ can be defined in a way that enables computing the indirect effect terms:

(9a) $S^{†}\equiv \frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\mathrm{\Phi }}={\mathrm{\Sigma }}_{R,pho.}$(9a)

and

(9b) $s^{†}\equiv \frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\psi }={\mathrm{\Sigma }}_{R,ele}.$(9b)

Multiplying Eqs. (8a)(8a) $T^{†}{\mathrm{\Phi }}^{†}=S^{†}+p^{†}{\psi }^{†}$(8a) and (8b)(8b) $t^{†}{\psi }^{†}=s^{†}+P^{†}{\mathrm{\Phi }}^{†}.$(8b) by perturbations in the photon and electron fluxes, respectively, and integrating over all phase-space give

(10a) $\delta \mathrm{\Phi }{T}^{†}{\mathrm{\Phi }}^{†}=\delta \mathrm{\Phi }{S}^{†}+\delta \mathrm{\Phi }{p}^{†}{\psi }^{†}$(10a)

and

(10b) $\delta \psi t^{†}{\psi }^{†}=\delta \psi s^{†}+\delta \psi P^{†}{\mathrm{\Phi }}^{†}.$(10b)

Applying the property of the adjoint (e.g., ${\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }=\delta \mathrm{\Phi }{T}^{†}{\mathrm{\Phi }}^{†}$) several times to Eq. (10)(10a) $\delta \mathrm{\Phi }{T}^{†}{\mathrm{\Phi }}^{†}=\delta \mathrm{\Phi }{S}^{†}+\delta \mathrm{\Phi }{p}^{†}{\psi }^{†}$(10a) gives Eq. (11)(11a) ${\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }=\delta \mathrm{\Phi }{S}^{†}+{\psi }^{†}p\delta \mathrm{\Phi }$(11a) :

(11a) ${\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }=\delta \mathrm{\Phi }{S}^{†}+{\psi }^{†}p\delta \mathrm{\Phi }$(11a)

and

(11b) ${\psi }^{†}t\delta \psi =\delta \psi s^{†}+{\mathrm{\Phi }}^{†}P\delta \psi .$(11b)

It should be noted that applying the adjoint property introduces a bilinear concomitant term that disappears upon applying vacuum boundary conditions.⁶ Equation (11)(11a) ${\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }=\delta \mathrm{\Phi }{S}^{†}+{\psi }^{†}p\delta \mathrm{\Phi }$(11a) will be revisited momentarily, but first, we will allow perturbations to occur to all terms in Eqs. (7a)(7a) $T\mathrm{\Phi }=Q_{total}=Q+P\psi ;$(7a) and (7b)(7b) $t\psi =q_{total}=q+p\mathrm{\Phi }.$(7b) ; ignore all higher-order ${\delta }^2$ terms; multiply the two equations by ${\mathrm{\Phi }}^{†}$ and ${\psi }^{†}$, respectively; and take the inner product to give

(12a) ${\mathrm{\Phi }}^{†}\delta T\mathrm{\Phi }+{\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }={\mathrm{\Phi }}^{†}\delta Q+{\mathrm{\Phi }}^{†}\delta P\psi +{\mathrm{\Phi }}^{†}P\delta \psi$(12a)

and

(12b) ${\psi }^{†}\delta t\psi +{\psi }^{†}t\delta \psi ={\psi }^{†}\delta q+{\psi }^{†}\delta p\mathrm{\Phi }+{\psi }^{†}p\delta \mathrm{\Phi }.$(12b)

Upon examination, one notices that Eq. (11a)(11a) ${\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }=\delta \mathrm{\Phi }{S}^{†}+{\psi }^{†}p\delta \mathrm{\Phi }$(11a) contains one term in both Eqs. (12a)(12a) ${\mathrm{\Phi }}^{†}\delta T\mathrm{\Phi }+{\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }={\mathrm{\Phi }}^{†}\delta Q+{\mathrm{\Phi }}^{†}\delta P\psi +{\mathrm{\Phi }}^{†}P\delta \psi$(12a) and (12b)(12b) ${\psi }^{†}\delta t\psi +{\psi }^{†}t\delta \psi ={\psi }^{†}\delta q+{\psi }^{†}\delta p\mathrm{\Phi }+{\psi }^{†}p\delta \mathrm{\Phi }.$(12b) and that Eq. (11b)(11b) ${\psi }^{†}t\delta \psi =\delta \psi s^{†}+{\mathrm{\Phi }}^{†}P\delta \psi .$(11b) contains one term in both Eqs. (12a)(12a) ${\mathrm{\Phi }}^{†}\delta T\mathrm{\Phi }+{\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }={\mathrm{\Phi }}^{†}\delta Q+{\mathrm{\Phi }}^{†}\delta P\psi +{\mathrm{\Phi }}^{†}P\delta \psi$(12a) and (12b)(12b) ${\psi }^{†}\delta t\psi +{\psi }^{†}t\delta \psi ={\psi }^{†}\delta q+{\psi }^{†}\delta p\mathrm{\Phi }+{\psi }^{†}p\delta \mathrm{\Phi }.$(12b) . These terms can be made to cancel by summing Eqs. (12a)(12a) ${\mathrm{\Phi }}^{†}\delta T\mathrm{\Phi }+{\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }={\mathrm{\Phi }}^{†}\delta Q+{\mathrm{\Phi }}^{†}\delta P\psi +{\mathrm{\Phi }}^{†}P\delta \psi$(12a) and (12b)(12b) ${\psi }^{†}\delta t\psi +{\psi }^{†}t\delta \psi ={\psi }^{†}\delta q+{\psi }^{†}\delta p\mathrm{\Phi }+{\psi }^{†}p\delta \mathrm{\Phi }.$(12b) and by then subtracting Eqs. (11a)(11a) ${\mathrm{\Phi }}^{†}T\delta \mathrm{\Phi }=\delta \mathrm{\Phi }{S}^{†}+{\psi }^{†}p\delta \mathrm{\Phi }$(11a) and (11b)(11b) ${\psi }^{†}t\delta \psi =\delta \psi s^{†}+{\mathrm{\Phi }}^{†}P\delta \psi .$(11b) , giving an updated form for Eq. (9)(7b) $t\psi =q_{total}=q+p\mathrm{\Phi }.$(7b) :

$\begin{array}{rl}& \frac{\delta R_{tot}^{Indirect}}{R_{tot}}=\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\mathrm{\Phi }}\delta \mathrm{\Phi }+\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\psi }\delta \psi =\delta \mathrm{\Phi }{S}^{†}+\delta \psi s^{†}\\ & ={\mathrm{\Phi }}^{†}\delta Q-{\mathrm{\Phi }}^{†}\delta T\mathrm{\Phi }+{\mathrm{\Phi }}^{†}\delta P\psi +{\psi }^{†}\delta q-{\psi }^{†}\delta t\psi \\ & +{\psi }^{†}\delta p\mathrm{\Phi }.(13)\end{array}$

This expression allows for the estimation of the indirect effect if the photon and electron adjoint fluxes are known.

II.A. Contributon Theory Formulation

We will now develop a method for computing the adjoint fluxes in Eq. (13)$\begin{array}{rl}& \frac{\delta R_{tot}^{Indirect}}{R_{tot}}=\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\mathrm{\Phi }}\delta \mathrm{\Phi }+\frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\psi }\delta \psi =\delta \mathrm{\Phi }{S}^{†}+\delta \psi s^{†}\\ & ={\mathrm{\Phi }}^{†}\delta Q-{\mathrm{\Phi }}^{†}\delta T\mathrm{\Phi }+{\mathrm{\Phi }}^{†}\delta P\psi +{\psi }^{†}\delta q-{\psi }^{†}\delta t\psi \\ & +{\psi }^{†}\delta p\mathrm{\Phi }.(13)\end{array}$. This method is based on Contributon Theory, which was developed originally by Williams⁷ and adapted to continuous-energy Monte Carlo simulations by Perfetti et al. via the CLUTCH method.³ To estimate the adjoint flux for eigenvalue sensitivity coefficient estimates, Contributon Theory first takes all sources of neutrons $Q$ and equates it to the system’s fission source weighted by the system eigenvalue $\lambda F\varphi$:

(14) $Q=\lambda F\varphi .$(14)

Next, we multiply both sides by the adjoint flux and take the inner product:

(15) $Q{\varphi }^{†}=\lambda {\varphi }^{†}F\varphi .$(15)

Next, we will reconsider our neutron source to represent one single particle traveling in phase-space $\tau$, i.e., $Q=Q_o\delta \left(\tau \right)$, to yield the following Green’s function:

(16) ${\varphi }^{†}\left(\tau \right)=\frac{\lambda }{Q_o}{\varphi }^{†}\left(r\right)F\left(r\right)\varphi \left(\tau \to r\right).$(16)

Equation (16)(16) ${\varphi }^{†}\left(\tau \right)=\frac{\lambda }{Q_o}{\varphi }^{†}\left(r\right)F\left(r\right)\varphi \left(\tau \to r\right).$(16) states that the importance of a particle traveling in some phase-space is determined by the number of fission neutrons it creates at location $r$ after leaving $\tau$, weighted by the importance of these daughter fission neutrons. In essence, a particle’s importance is determined by how much response it and any daughter particles create in the future.

We will now adapt the existing Contributon Theory interpretation for determining the importance of individual particles to coupled transport simulations. Weighting Eqs. (7a)(7a) $T\mathrm{\Phi }=Q_{total}=Q+P\psi ;$(7a) and (7b)(7b) $t\psi =q_{total}=q+p\mathrm{\Phi }.$(7b) by ${\mathrm{\Phi }}^{†}$ and ${\psi }^{†}$, respectively, and taking the inner product give

(17a) ${\mathrm{\Phi }}^{†}T\mathrm{\Phi }={\mathrm{\Phi }}^{†}{Q}_{total}$(17a)

and

(17b) ${\psi }^{†}t\psi ={\psi }^{†}{q}_{total}.$(17b)

Next, the adjoint balance equations in Eqs. (8a)(8a) $T^{†}{\mathrm{\Phi }}^{†}=S^{†}+p^{†}{\psi }^{†}$(8a) and (8b)(8b) $t^{†}{\psi }^{†}=s^{†}+P^{†}{\mathrm{\Phi }}^{†}.$(8b) are weighted by $\mathrm{\Phi }$ and $\psi$, respectively; the inner product is taken; and the property of the adjoint is applied. Noting that the left-hand sides of Eqs. (8a)(8a) $T^{†}{\mathrm{\Phi }}^{†}=S^{†}+p^{†}{\psi }^{†}$(8a) and (8b)(8b) $t^{†}{\psi }^{†}=s^{†}+P^{†}{\mathrm{\Phi }}^{†}.$(8b) and those of Eqs. (17a)(17a) ${\mathrm{\Phi }}^{†}T\mathrm{\Phi }={\mathrm{\Phi }}^{†}{Q}_{total}$(17a) and (17b)(17b) ${\psi }^{†}t\psi ={\psi }^{†}{q}_{total}.$(17b) are equivalent gives

(18a) ${\mathrm{\Phi }}^{†}{Q}_{total}=\mathrm{\Phi }{S}^{†}+{\psi }^{†}p\mathrm{\Phi }$(18a)

and

(18b) ${\psi }^{†}{q}_{total}=\psi s^{†}+{\mathrm{\Phi }}^{†}P\psi .$(18b)

We again consider a photon or electron source that consists of one single particle traveling in phase-space; i.e., $Q_{total}\left(\tau \right)=Q_o\delta \left(\tau \right)$ or $q_{total}\left(\tau \right)=q_o\delta \left(\tau \right)$. It should be noted that this source is either one photon or one electron, but not a combination of the two. This assumption allows one to compute the adjoint flux for photons and electrons using the following Green’s function interpretation:

(19a) ${\mathrm{\Phi }}^{†}\left(\tau \right)=\mathrm{\Phi }\left(\tau \to r\right)S^{†}\left(r\right)+{\psi }^{†}\left(r\right)p\left(r\right)\mathrm{\Phi }\left(\tau \to r\right)$(19a)

and

(19b) ${\psi }^{†}\left(\tau \right)=\psi \left(\tau \to r\right)s^{†}\left(r\right)+{\mathrm{\Phi }}^{†}\left(r\right)P\left(r\right)\psi \left(\tau \to r\right).$(19b)

Similarly to Eq. (16)(16) ${\varphi }^{†}\left(\tau \right)=\frac{\lambda }{Q_o}{\varphi }^{†}\left(r\right)F\left(r\right)\varphi \left(\tau \to r\right).$(16) , Eqs. (19a)(19a) ${\mathrm{\Phi }}^{†}\left(\tau \right)=\mathrm{\Phi }\left(\tau \to r\right)S^{†}\left(r\right)+{\psi }^{†}\left(r\right)p\left(r\right)\mathrm{\Phi }\left(\tau \to r\right)$(19a) and (19b)(19b) ${\psi }^{†}\left(\tau \right)=\psi \left(\tau \to r\right)s^{†}\left(r\right)+{\mathrm{\Phi }}^{†}\left(r\right)P\left(r\right)\psi \left(\tau \to r\right).$(19b) show that a particle’s importance is determined by how much both that particle and its secondary particle progeny contribute to some response of interest. The definitions for $S^{†}$ and $s^{†}$ in Eqs. (9a)(9a) $S^{†}\equiv \frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\mathrm{\Phi }}={\mathrm{\Sigma }}_{R,pho.}$(9a) and (9b)(9b) $s^{†}\equiv \frac{1}{R_{tot}}\frac{\mathrm{\partial }{R}_{tot}}{\mathrm{\partial }\psi }={\mathrm{\Sigma }}_{R,ele}.$(9b) allow the first inner products on the right-hand side of Eqs. (19a)(19a) ${\mathrm{\Phi }}^{†}\left(\tau \right)=\mathrm{\Phi }\left(\tau \to r\right)S^{†}\left(r\right)+{\psi }^{†}\left(r\right)p\left(r\right)\mathrm{\Phi }\left(\tau \to r\right)$(19a) and (19b)(19b) ${\psi }^{†}\left(\tau \right)=\psi \left(\tau \to r\right)s^{†}\left(r\right)+{\mathrm{\Phi }}^{†}\left(r\right)P\left(r\right)\psi \left(\tau \to r\right).$(19b) to be computed by tallying the contribution of the original particle in phase-space $\tau$ to the responses of interest after it is emitted. In practice, this means that the first component to a particle’s importance is equal to the response ($R_{pho.}$ or $R_{ele.}$) that is generated by the particle from the time it leaves the event until its death by weighting this future flux by $S^{†}$ or $s^{†}$. The second inner products on the right-hand side of Eqs. (19a)(19a) ${\mathrm{\Phi }}^{†}\left(\tau \right)=\mathrm{\Phi }\left(\tau \to r\right)S^{†}\left(r\right)+{\psi }^{†}\left(r\right)p\left(r\right)\mathrm{\Phi }\left(\tau \to r\right)$(19a) and (19b)(19b) ${\psi }^{†}\left(\tau \right)=\psi \left(\tau \to r\right)s^{†}\left(r\right)+{\mathrm{\Phi }}^{†}\left(r\right)P\left(r\right)\psi \left(\tau \to r\right).$(19b) rely on computing the importance of any secondary electrons or photons that are created by the primary particles by weighting their future fluxes by $S^{†}$ or $s^{†}$. This process is then repeated as necessary for tertiary particles, quaternary particles, etc., until all importance generated by daughter particles has been accounted for.

III. PROOF OF PRINCIPLE

A simple Monte Carlo test code was developed to demonstrate proof of principle for the proposed CLUTCH methodology of computing GPT reaction rate sensitivity coefficients in coupled radiation transport simulations. The Monte Carlo code examined the sensitivity of responses in a simple 1-D slab geometry that occurred in response to an incident flux of radiation—here assumed to represent photons and electrons. The results presented here build on previously published, initial results.⁸ These particles were assumed to have an initial direction that was normal to the surface of the slab and scattered isotropically once inside the slab. This slab was divided into ten 1-cm-thick regions, each with identical cross sections and density, for an overall slab thickness of 10 cm. These simulations assumed a three-group energy format and used artificial nuclear data for the simulation cross sections, emission probabilities, scattering kinematics, and secondary particle production probabilities. As mentioned in Eqs. (7a)(7a) $T\mathrm{\Phi }=Q_{total}=Q+P\psi ;$(7a) and (7b)(7b) $t\psi =q_{total}=q+p\mathrm{\Phi }.$(7b) , it was assumed that photons could create daughter electrons and that electrons could make daughter photons.

As shown in Table I, the electron absorption cross sections were approximately three times as large as photon absorption cross sections in hopes of emulating true physics; in truth, electron absorption cross sections are often much more than three times larger than photon absorption cross sections, but these cross sections were chosen to be much closer in magnitude to increase the degree of coupling between the two types of radiation and to make any defects in the proposed method more readily apparent. Table I also gives the distribution and spectra of the photons and electrons incident on the slab. Table II gives the photon and electron group-to-group scattering probabilities. Table III gives the probability of an absorption reaction producing a secondary particle; note that these secondary particle production probabilities were chosen to be large enough to sufficiently couple the behavior of photons and electrons. In Tables I, II, and III, Group 1 is intended to serve as the highest-energy group, and Group 3 is the lowest. Figure 1 plots the photon and electron flux solutions throughout the slab. This problem is solved deterministically by coupling the particle types in the scattering matrix.

TABLE I Sample Problem Cross Sections and Source Spectra

	Photons (80% of Particle Source)				Electrons (20% of Particle Source)
	Total	Scatter	Absorption	Source Spectrum	Total	Scatter	Absorption	Source Spectrum
Group 1	1.00	0.70	0.30	0.70	2.50	1.50	1.00	1.00
Group 2	2.00	1.00	1.00	0.20	4.00	1.00	3.00	0.00
Group 3	5.00	1.00	4.00	0.10	6.00	1.00	5.00	0.00

TABLE II Sample Problem Scattering Kinematics

Photon Outgoing Group	Photon Incident Group			Electron Outgoing Group	Electron Incident Group
	Group 1	Group 2	Group 3		Group 1	Group 2	Group 3
Group 1	0.90	0.10	0.00	Group 1	0.60	0.10	0.05
Group 2	0.05	0.40	0.20	Group 2	0.30	0.50	0.15
Group 3	0.05	0.50	0.80	Group 3	0.10	0.40	0.80

TABLE III Sample Problem Secondary Production Probabilities

Electron Birth Group	Photon Incident Group			Photon Birth Group	Electron Incident Group
	Group 1	Group 2	Group 3		Group 1	Group 2	Group 3
Group 1	0.25	0.10	0.00	Group 1	0.05	0.00	0.00
Group 2	0.15	0.20	0.25	Group 2	0.15	0.15	0.00
Group 3	0.15	0.30	0.25	Group 3	0.20	0.15	0.10

Fig. 1. Photon and electron fluxes in a sample 1-D slab.

The CLUTCH methodology for estimating GPT response sensitivities was implemented in the test Monte Carlo code and was used to estimate the sensitivity of GPT responses to the various nuclear data in this simulation. These simulations assumed the GPT responses of interest to be particle absorption rates in each of the 10 regions; i.e.,

(20) $R={\mathrm{\Sigma }}_{Absorb}\varphi { }_{Celli}.$(20)

Responses and response sensitivities were computed for each of the 10 cell regions for either photon or electron absorption rates. In practice, applications are more likely to compute the sensitivity of responses that sum both photon and electron reaction rates, but this work instead examined the sensitivity of each individual particle reaction rate to uncertain nuclear data. Examining these responses individually allows for more straightforward debugging and more apparently reveals any defects in the method’s theory or implementation, as a method that incorrectly computes the sensitivity of the individual responses will not succeed in computing the sensitivity of the sum of two responses.

Direct perturbation simulations were used to obtain reference sensitivity coefficient estimates. For the sake of time, these direct perturbation calculations were performed for only the photon cross-section data. Select perturbations were performed for select electron data and successfully confirmed the accuracy of these electron data sensitivity coefficients; however, resolving the stochastic noise in the electron data perturbations generally required simulating significantly more particle histories due to the smaller magnitude of these electron sensitivities. The electron sensitivities share the exact same methodology and implementation as the photon sensitivities; they differ only due to the magnitude of the arbitrary, artificial nuclear data used in these simulations. Therefore, it is reasonable to assume that the electron sensitivities are accurate if the photon sensitivities can be found to be accurate.

The SENSMG sensitivity analysis wrapper code enables the PARTISN deterministic transport code to estimate sensitivity coefficients.⁹ Although SENSMG was originally intended for sensitivity analysis of neutron-gamma coupled transport simulations, the test problem’s physics are sufficiently particle-agnostic to allow PARTISN to simulate the test problem. Figure 1 includes plots of the SENSMG photon and electron flux solutions for the test problem, and all future sensitivity plots include SENSMG sensitivity coefficient estimates; these SENSMG results are presented as solid squares in each of these plots.

IV. RESULTS

Figure 2 compares the CLUTCH, direct perturbation, and SENSMG estimates for the sensitivity of photon absorption reaction rates throughout the 1-D slab to photon Group 1 total, absorption, and scattering cross sections. The CLUTCH sensitivity estimates are plotted using solid lines with red error bars representing the 1σ stochastic uncertainty intervals, and the direct perturbation estimates are plotted using dashed lines with orange error bars. The CLUTCH error bars are often not visible in this and later sensitivity plots because they are generally very small in magnitude. Shaded area plots in Fig. 2 depict the difference between the CLUTCH and the direct perturbation sensitivity estimates in terms of the number of effective standard deviations. Comparing sensitivity coefficients using the number of standard deviations of difference is generally more effective than comparing sensitivities using relative error because sensitivities that happen to be close to zero produce very large relative errors, even if the sensitivity methods and implementations are accurate.

Fig. 2. Sensitivity of the photon absorption reaction rate to photon Group 1 total, absorption, and scattering cross sections.

The sensitivity coefficients in Fig. 2 are initially positive in magnitude near the front of the slab, as increasing cross-section magnitudes near the front of the slab will stop particles from escaping these front cells. The sensitivities grow larger and more negative deeper into the shield, as increasing the magnitude of these cross sections throughout the entire slab will strongly depress the flux that reaches these regions.

Figure 2 indicates that the CLUTCH, direct perturbation, and SENSMG sensitivity coefficients agree very well. There is virtually no observable difference between the sensitivity coefficients near the front of the slab, and the sensitivities only begin to differ deep into the slab as the uncertainty in the direct perturbation sensitivities increases due to stochastic noise. Sensitivity coefficients can generally be considered inaccurate if they differ by more than 2σ or 3σ, and all sensitivities in Fig. 2 are within approximately 1σ of agreement.

Figure 3 shows the sensitivity of the electron absorption rate to photon Group 1 total, absorption, and scattering cross sections. The CLUTCH and direct perturbation sensitivity coefficients generally agree well in Fig. 3, but some significant disagreement (−4.42$\sigma$) exists in Cell 5. This disagreement is presumably due to inadequate convergence in the direct perturbation sensitivity coefficients and their uncertainty estimates, as the CLUTCH and SENSMG sensitivity coefficients agree well for Cell 5.

Fig. 3. Sensitivity of the electron absorption reaction rate to photon Group 1 total, absorption, and scattering cross sections.

Obtaining high-fidelity direct perturbation reference sensitivities can be difficult: Nuclear data perturbations must be large enough to be resolved beyond stochastic noise, but very large perturbations may introduce inaccuracy due to higher-order effects. The optimal size for nuclear direct nuclear data perturbations depends on the magnitude of the sensitivity being verified, and verifying sensitivities individually can lead to a significant number of direct perturbation calculations. This study performed roughly 200 direct perturbation calculations and verified groups of sensitivity coefficients together; verifying each sensitivity coefficient individually would have required approximately 3000 direct perturbation calculations. Therefore, these reference direct perturbation sensitivity coefficients contain some small degree of inaccuracy. As will be discussed later, this inaccuracy was more pronounced for direct perturbation sensitivity coefficients that are smaller in magnitude and thus more difficult to resolve.

Errors in sensitivity methods are usually abundantly clear when prototyping sensitivity methods—differences between flawed sensitivity coefficients and reference direct perturbation coefficients are generally on the order of 10s, 100s, or 1000s of standard deviations, not in the range of 4$\sigma$ or 5$\sigma$. Furthermore, these differences grow as simulations are performed with longer run times and thus smaller stochastic errors. The differences observed here between the CLUTCH and direct perturbation sensitivity estimates differ by a somewhat large, yet still plausible, degree and did not grow after longer CLUTCH simulations were performed. Furthermore, the direct perturbation sensitivity coefficients in Fig. 3 behave suspiciously, oscillating above and below the steadier CLUTCH sensitivity estimates. The CLUTCH and SENSMG sensitivity coefficients exhibit consistent, strong agreement, even when the direct perturbation sensitivity coefficients contain significant stochastic noise.

Figure 4 compares the sensitivity of the photon absorption reaction rate to photon Group 2 cross-section data. Like the Group 1 sensitivities, the Group 2 sensitivities are initially positive and become negative farther into the slab, but the Group 2 sensitivities do not continue becoming more and more negative; in fact, the opposite occurs: They become less negative and closer to a sensitivity of zero. Similar behavior is seen for the sensitivity of the electron absorption reaction rate to photon Group 2 cross sections in Fig. 5. This behavior likely occurs because Group 2 absorption cross sections are three times greater than the Group 1 absorption cross sections, meaning that any Group 2 particles present in the deeper regions of the slab are almost entirely produced either by downscatter of Group 1 particles or are secondary radiation. Thus, Group 2 particles do not produce large, negative sensitivities due to self-shielding because very few Group 2 particles that are incident on the surface of the slab penetrate deeply into the slab.

Fig. 4. Sensitivity of the photon absorption reaction rate to photon Group 2 total, absorption, and scattering cross sections.

Fig. 5. Sensitivity of the electron absorption reaction rate to photon Group 2 total, absorption, and scattering cross sections.

The CLUTCH, direct perturbation, and SENSMG sensitivity coefficients exhibited good agreement for the Group 2 cross-section sensitivities. Most of the sensitivities were within 1σ or 2σ, with several sensitivities differing by more than 4σ. The large error bars in Figs. 4 and 5 indicate that a significant amount of stochastic noise was present in the direct perturbation Group 2 cross-section sensitivities. These sensitivities were much smaller in magnitude than the Group 1 cross-section sensitivities; therefore, it was much more difficult to obtain accurate and well-resolved direct perturbation sensitivity estimates. Nonetheless, the higher-resolution direct perturbation sensitivities at the front of the shield, the generally good agreement between the CLUTCH and direct perturbation sensitivities deeper into the shield, and the consistent strong agreement with SENSMG sensitivities suggest that the CLUTCH sensitivity estimates are accurate.

The trends observed for the Group 2 cross-section sensitivity coefficients, including the significant impact of stochastic noise, were also present in the Group 3 cross-section sensitivity estimates, as shown in Figs. 6 and 7. Stochastic noise affected these sensitivities to a greater degree; in fact, the CLUTCH sensitivities, which barely produced visible error bars for the Group 1 or 2 sensitivity estimates, produced visible error bars for the Group 3 sensitivity estimates. Nonetheless, these sensitivities appeared to agree well through the entire slab.

Fig. 6. Sensitivity of the photon absorption reaction rate to photon Group 3 total, absorption, and scattering cross sections.

Fig. 7. Sensitivity of the electron absorption reaction rate to photon Group 3 total, absorption, and scattering cross sections.

Figures 8 and 9 show the sensitivity of the photon and electron absorption reaction rates, respectively, to photon-to-electron secondary production cross sections. Incident particles are unlikely to reach the deeper regions of the slab unattenuated; therefore, the absorption rate in these deeper regions is attributable primarily to the secondary radiation production. This self-shielding causes the secondary particle production sensitivities to be high deeper into the shield and lower at the front of the shield where unattenuated incident radiation is still abundant. Indeed, Fig. 9 shows that the electron absorption rate in the deepest regions of the shield is almost 100% sensitive to the production of secondary electrons by photons; in other words, electrons are present in these regions only because they are progeny of more penetrating photons. Figure 8 shows that the production of secondary electrons by photons also contributes significantly to the photon absorption rate in deeper regions of the shield, as these secondary electrons can affect the photon absorption rate by producing tertiary photons.

Fig. 8. Sensitivity of the photon absorption reaction rate to photon → electron secondary production cross sections.

Fig. 9. Sensitivity of the electron absorption reaction rate to photon → electron secondary production cross sections.

The CLUTCH sensitivity coefficients for secondary particle production cross sections once again agree well with the reference direct perturbation and SENSMG sensitivity coefficients, and stochastic noise again impedes verifying the accuracy of these sensitivities in the deepest regions of the slab.

Figures 10 and 11 show the sensitivity of the photon and electron absorption reaction rates, respectively, to the spectrum of incident photon radiation that is incident on the slab. These figures plot the “unconstrained” spectral sensitivity coefficients because sensitivity method debugging is generally more straightforward and effective when examining unconstrained sensitivity coefficients.¹⁰ Because distributions must be normalized, increasing the probability that radiation will be emitted at one energy requires a corresponding decrease in the probability that radiation will be emitted at other energies. This relationship between different emission probabilities within the same probability distribution necessitates relating the sensitivity coefficients for emission probabilities within that distribution by “constraining” them; however, it is easier to verify the accuracy of the unconstrained sensitivities because they do not contain the noise introduced by the sensitivity constraining process.

Fig. 10. Sensitivity of the photon absorption reaction rate to the photon source spectrum.

Fig. 11. Sensitivity of the electron absorption reaction rate to the photon source spectrum.

As observed for the secondary particle production sensitivity coefficients in Figs. 8 and 9, the spectral sensitivity coefficients for Group 1 photons grow in magnitude in the deeper regions of the slab; in fact, these spectral sensitivities approach a value of 1.00. This again is because Group 1 photons are significantly more likely to penetrate the deeper regions of the shield and to contribute to the photon and (through secondary particle production) the electron absorption rates.

The CLUTCH, direct perturbation, and SENSMG sensitivity coefficients again exhibit strong agreement, and stochastic noise again affects the sensitivity comparison in deeper regions of the shield. This noise is apparent in the direct perturbation sensitivities for Cell 9 in both Figs. 10 and 11: Spectral sensitivities can generally reach a maximum value of 1.00, which corresponds to a case when a response is 100% attributed to that type and energy of radiation, yet the direct perturbation sensitivities for Cell 9 in Figs. 10 and 11 exceed a value of 1.00 due to stochastic noise.

V. CONCLUSIONS

This work has outlined an extension of GPT sensitivity analysis and Contributon Theory that enables adjoint-based sensitivity coefficient calculations in coupled particle Monte Carlo transport simulations. This new sensitivity methodology was implemented into a simple 1-D test Monte Carlo code, and direct perturbation and deterministic SENSMG calculations were used to confirm the accuracy of the new sensitivity methodology. In general, the CLUTCH and direct perturbation sensitivities agreed well, and the few instances of nonnegligible disagreement were attributed to stochastic noise in the direct perturbation sensitivities. Future studies could further resolve these direct perturbation reference sensitivities and remove the few instances where sensitivities disagreed by more than 3σ, but the strong agreement between the CLUTCH and the SENSMG results and the authors’ previous experience from benchmarking sensitivity methods suggest that the coupled CLUTCH sensitivity methodology produces accurate sensitivity coefficient estimates.

Future work includes implementing this methodology in a production-level coupled transport Monte Carlo code, such as the Integrated Tiger Series (ITS) code¹¹ and using this methodology to quantify the impact of uncertain photon and electron nuclear data on photon and electron radiation transport applications.

Acknowledgments

The employee owns all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan.

The authors are grateful to Jeffrey A. Favorite, Los Alamos National Laboratory, for performing the SENSMG sensitivity calculations presented in this work.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE).