Fourth-Order Adjoint Sensitivity Analysis of an OECD/NEA Reactor Physics Benchmark: I. Mathematical Expressions and CPU-Time Comparisons for Computing 1^st-, 2^nd- and 3^rd-Order Sensitivities

Dan Gabriel Cacuci^*, Ruixian Fang

Center for Nuclear Science and Energy, Department of Mechanical Engineering, University of South Carolina, Columbia, SC, USA.

DOI: 10.4236/ajcm.2021.112009                

Abstract

This work extends to fourth-order previously published work on developing the adjoint sensitivity and uncertainty analysis of the numerical model of a polyethylene-reflected plutonium (acronym: PERP) OECD/NEA reactor physics benchmark. The PERP benchmark comprises 7477 imprecisely known (uncertain) model parameters which have nonzero values. These parameters are as follows: 180 microscopic total cross sections; 7101 microscopic scattering sections; 60 microscopic fission cross sections; 60 parameters that characterize the average number of neutrons per fission; 60 parameters that characterize the fission spectrum; 10 parameters that characterize the fission source; and 6 parameters that characterize the isotope number densities. Previous works have used the adjoint sensitivity analysis methodology to compute exactly and efficiently all of the 7477 first-order and 27,956,503 second-order sensitivities of the PERP benchmark’s leakage response to all of the benchmark’s uncertain parameters. These works showed that largest response sensitivities involve the total microscopic cross sections, which motivated the recent computation of all of the (180)³ third-order sensitivities of the PERP leakage response with respect to these total microscopic cross sections. It turned out that some of these 3^rd-order cross sections were far larger than the corresponding 2^nd-order ones, thereby having the largest impact on the uncertainties induced in the PERP benchmark’s response. This finding has motivated the development of the original 4^th-order formulas presented in this work, which are valid not only for the PERP benchmark but can also be used for computing the 4^th-order sensitivities of response of any nuclear system involving fissionable material and internal or external neutron sources. Subsequent works will use the adjoint-based mathematical expressions obtained in this work to compute exactly and efficiently the numerical values of the largest fourth-order sensitivities of the PERP benchmark’s response to the total microscopic cross section and use them for a pioneering fourth-order uncertainty analysis of the PERP benchmark’s response.

Keywords

Polyethylene-Reflected Plutonium Sphere Benchmark, 1^st-Order, 2^nd-Order, 3^rd-Order Forward Sensitivity Analysis, Finite-Differences, 1^st-Order, 2^nd-Order, 3^rd-Order Adjoint Sensitivity Analysis, Cross Sections

1. Introduction

Until recently, only the first-order sensitivities (i.e., functional derivatives) of a computational model’s responses (i.e., quantities of interest) to the respective model’s imprecisely known (i.e., uncertain) parameters have been considered when assessing the uncertainties induced in the respective responses by the parameter uncertainties. The second- and higher-order sensitivities could not be computed, except for very simple models comprising a handful of parameters, so these sensitivities were ignored. The Second-Order Adjoint Sensitivity Analysis Methodology (2^nd-ASAM) recently conceived by Cacuci [1] is the only practical method that enables the exact computation of the large number of 2^nd-order sensitivities arising in large-scale problems comprising many parameters. The application of the 2^nd-ASAM to a multiplying nuclear system with source [2] [3] [4] has opened the way for the large-scale application presented in [5] - [10] to a polyethylene-reflected plutonium (acronym: PERP) OECD/NEA reactor physics benchmark [11]. The numerical model of the PERP benchmark comprises 21,976 uncertain parameters, of which 7477 parameters have nonzero values, as follows: 180 group-averaged total microscopic cross sections, 7101 non-zero group-averaged scattering microscopic cross sections (the other scattering cross sections, of which there are 21,600 in total, are zero); 120 fission process parameters; 60 fission spectrum parameters; 10 parameters describing the experiment’s nuclear sources; and 6 isotopic number densities.

All of the non-zero first-order sensitivities and second-order sensitivities of the PERP leakage response with respect to the benchmark’s parameters were computed, ranked, and analyzed in [5] - [10]. The results obtained in [5] - [10] showed that the 2^nd-order sensitivities of the leakage response with respect to the group-averaged microscopic total cross sections are the largest (by comparison to the 2^nd-order sensitivities of the leakage response with respect to the other uncertain model parameters) and have therefore the largest impact on the uncertainties induced in the leakage response. For example, neglecting the 2^nd-order sensitivities of the leakage response with respect to the total cross sections could cause an error as large as 2000% in the expected value of the leakage response and up to 6000% in the variance of the leakage response [5]. Given that the effects of these 2^nd-order sensitivities are much larger than the effects of the 1^st-order sensitivities [5] [10], it is logical to quantify the magnitudes and contributions that would stem from the 3^rd-order sensitivities of the PERP benchmark’s total leakage response with respect to the microscopic total cross sections. To enable the computation of such 3^rd-order sensitivities, Cacuci [12] has recently conceived the “third order adjoint sensitivity analysis methodology for reaction rate responses in a multiplying nuclear system with source” and then applied this general theory to the PERP benchmark in order to derive the exact analytical expressions of the 3^rd-order sensitivities of the PERP benchmark’s leakage response with respect to this benchmark’s microscopic total cross sections [13] [14] [15]. These works [13] [14] [15] have shown that the largest 3^rd-order sensitivities involve the microscopic total cross section for the lowest (30^th) energy group of isotope ¹H (i.e., ${\sigma }_{t,6}^{30}$ ). In particular, the largest overall 3^rd-order sensitivity is the mixed 3^rd-order sensitivity $S^{\left(3\right)}\left({\sigma }_{t,1}^{g=30},{\sigma }_{t,6}^{g^{\prime }=30},{\sigma }_{t,6}^{g^{″}=30}\right)=-1.88\times {10}^5$, which also involves the microscopic total cross section for the 30^th energy group of isotope ²³⁹Pu (i.e., ${\sigma }_{t,1}^{g=30}$ ). The largest unmixed 3^rd-order sensitivity is also with respect to ${\sigma }_{t,6}^{30}$, namely $S^{\left(3\right)}\left({\sigma }_{t,6}^{g=30},{\sigma }_{t,6}^{g^{\prime }=30},{\sigma }_{t,6}^{g^{″}=30}\right)=-2.966\times {10}^4$. These sensitivities significantly larger than the largest lower-order sensitivities; by comparison, the largest 1^st-, 2^nd- and 3^rd-order sensitivities are $S^{\left(1\right)}\left({\sigma }_{t,6}^{30}\right)=-9.366$, $S^{\left(2\right)}\left({\sigma }_{t,6}^{30},{\sigma }_{t,6}^{30}\right)=429.6$ and $S^{\left(3\right)}\left({\sigma }_{t,1}^{g=30},{\sigma }_{t,6}^{g^{\prime }=30},{\sigma }_{t,6}^{g^{″}=30}\right)=-1.88\times {10}^5$, respectively. The results obtained in [5] - [10] and [13] [14] [15] indicated that the total microscopic cross section of isotopes ¹H and ²³⁹Pu are the most important parameters affecting the PERP benchmark’s leakage response, since they are responsible for the largest 1^st-, 2^nd- and 3^rd-order sensitivities.

This work is organized as follows: Section 2 presents the computational modeling of the PERP Benchmark. Section 3 presents the main features and the CPU-times required for applying the three fundamental deterministic (as opposed to “statistical”) methods for computing the 1^st-order response sensitivities, namely 1) finite differences; 2) the forward sensitivity analysis method (FSAM); and 3) the 1^st-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1^st-CASAM) for linear systems. Section 4 presents the main features and the CPU-times required for applying the three fundamental deterministic methods (i.e., finite differences, FSAM the 2^nd-CASAM, respectively) for computing the 2^nd-order response sensitivities. The results presented in Section 4 clearly highlight the fact that the 2^nd-CASAM is the only method that can be used in practice to compute 1^st- and 2^nd-order sensitivities of a response of a large-scale system, involving many parameters. Section 5 presents comparisons of the CPU times required for computing the 3^rd-order sensitivities of the PERP benchmark’s leakage response solely to the group-averaged total cross sections. These results underscore the fact that the 3^rd-CASAM is exact, introducing no intrinsic methodological errors in the computation of sensitivities, and is by far more efficient than any other method. Section 6 concludes this work, noting that the expressions and results for the 4^th-order sensitivities will be presented in companion works [16] [17]. Notably, the original mathematical concepts underlying the 1^st-order FSAM and 1^st-CASAM were presented in [18] [19] [20] while the general theory underlying the 4^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology (4^th-CASAM) for linear systems is presented in [21].

2. Mathematical/Computational Model of the PERP Benchmark

The spherical polyethylene-reflected plutonium (acronym: PERP) benchmark comprises a metallic inner sphere (“core”) containing the following 4 isotopes: Isotope 1 (²³⁹Pu), Isotope 2 (²⁴⁰Pu), Isotope 3 (⁶⁹Ga) and Isotope 4 (⁷¹Ga). This core (which is designated as “material 1”) is surrounded by a spherical shell of polyethylene (designated as “material 2”), containing two isotopes, designated as Isotope 5 (C) and Isotope 6 (¹H), respectively. The dimensions and material composition of the polyethylene-reflected plutonium (PERP) metal sphere considered in this work are presented in TableA1 in the Appendix. The quantity of interest in this work, which will be called the “response,” is the leakage of neutrons out of the PERP sphere, which has been measured experimentally [11].

The numerical modeling of the neutron flux distribution within the PERP benchmark, as well as the computation of the PERP leakage response has been performed using the multigroup discrete ordinates particle transport code PARTISN [22] together with neutron sources computed using the code SOURCES4C [23]. The neutron flux distribution within the PERP benchmark, as well as the leakage of neutrons out of the PERP sphere has been modeled using the standard multigroup form of the Boltzmann neutron transport equation subject to the boundary condition of no incoming flux, with an internal spontaneous fission source, which can be written in the following form:

$\begin{array}{l}\Omega \cdot \nabla {\phi }^g\left(r,\Omega \right)+{\Sigma }_t^g\left(\alpha ;r\right){\phi }^g\left(r,\Omega \right)-{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }{\phi }^h\left(r,{\Omega }^{\prime }\right)[{\Sigma }_s^{h\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)}}\\ +{\chi }^g\left(\alpha ;r\right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha ;r\right)]=Q^g\left(\alpha \right),g=1,\cdots ,G,\end{array}$ (1)

${\phi }^g\left(r,\Omega \right)=0,r=r_d,\Omega \cdot n<0,g=1,\cdots ,G$, (2)

where $r_d$ is the radius of the PERP sphere. Note that many works, use the “ $g^{\prime }$ ” as the “dummy” summation index in Equation (1); the index h (instead of $g^{\prime }$ ) is preferred in this work in order to avoid, as much as possible, the appearance of “superscripted superscripts.” In this work, therefore, the superscript-indices “g” and “h” will refer exclusively to “energy groups.”

The source, $Q^g\left(\alpha \right)$, which appears on the right-side of Equation (1) is defined as follows:

$Q^g\left(\alpha \right)\triangleq {\displaystyle \underset{k=1}{\overset{N_f}{\sum }}{\lambda }_k{N}_k{F}_k^{SF}{\nu }_k^{SF}\left(\frac{2}{\sqrt{\pi a_k^3{b}_k}}{\text{e}}^{-\frac{a_k{b}_k}{4}}\right){\displaystyle {\int }_{E^{g+1}}^{E^g}\text{d}E{\text{e}}^{-E/a_k}\sinh \sqrt{b_{k}E}}}$, (3)

where $N_f$ denotes the total number of spontaneous-fission isotopes. The PERP benchmark has $N_f=2$ spontaneous-fission isotopes, namely “isotope 1” (²³⁹Pu) and “isotope 2” (²⁴⁰Pu). The spontaneous fission neutron spectrum of ²³⁹Pu and ²⁴⁰Pu, respectively, is each approximated by a Watts fission spectrum, characterized by: 1) the evaluated parameters $a_k$ and $b_k$; 2) the spontaneous emission, ${\nu }_k^{SF}$, of an average neutron of an isotope k; 3) the decay constant, ${\lambda }_k$, for actinide nuclide k; 4) the fraction, $F_k^{SF}$, of decays that are spontaneous fission (the “spontaneous-fission branching fraction”); 5) the index $k=1$ refers to “isotope 1” (²³⁹Pu) while $k=2$ refers to “isotope 2” (²⁴⁰Pu). The symbols used in Equations (1) through (3) have their usual meanings; for convenient referencing, they are described in the Appendix and are summarized in the Nomenclature.

The mathematical expression of the PERP benchmark’s leakage response, denoted as $L\left(\alpha \right)$, is provided below:

$L\left(\alpha \right)\triangleq {\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n{\phi }^g\left(r,\Omega \right)$. (4)

The PARTISN [22] computations of the neutron flux used the MENDF71X [24] 618-group cross section data collapsed to $G=30$ energy groups, as well as a P₃ Legendre expansion of the scattering cross section and a fine-mesh spacing of 0.005 cm (comprising 759 meshes for the plutonium sphere of radius of 3.794 cm, and 762 meshes for the polyethylene shell of thickness of 3.81 cm). The first- and second-order response sensitivities were computed using an angular quadrature of S₂₅₆. The 3^rd- and 4^th-order sensitivities of the leakage response with respect to the total cross sections were computed using an angular quadrature of S₃₂. The group boundaries of the $G=30$ energy groups are provided in TableA2 in the Appendix. The scattering and fission terms in Equation (1) contain implicitly a factor $1/4\pi$, to conform to the convention used in PARTISN [22].

The vector $\alpha$, which appears in the arguments of the various quantities in Equations (1) and (4), is defined as follows:

$\begin{array}{l}\alpha \triangleq {\left[{\alpha }_1,\cdots ,{\alpha }_{TP}\right]}^{†}\triangleq {\left[q;N;{\sigma }_s;{\sigma }_t;{\sigma }_f;\nu ;p\right]}^{†},\\ \text{where}TP\triangleq JQ+I+JSX+JTX+JFX+JNU+J\chi .\end{array}$ (5)

The components of the TP-dimensional vector $\alpha$ are the imprecisely known (i.e., uncertain) model and response parameters, and are described in the Appendix. In Equation (5) and throughout this work, the dagger will be used to denote “transposition.” The nominal values of the parameters in Equation (5), as well as of all other quantities in this work, will be denoted using the superscript “zero,” e.g., ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\cdots ,{\alpha }_{TP}^0\right]}^{†}$.

For the mathematical derivations to follow in this work, it will be convenient to use the matrix-form of Equations (1) and (2), which is as follows:

$B\left(\alpha \right)\phi \left(r,\Omega \right)=Q\left(\alpha \right),$ (6)

$\phi \left(r,\Omega \right)=0,r=r_d,\Omega \cdot n<0,$ (7)

where

$B\left(\alpha \right)\triangleq \left(\begin{array}{ccccc}{B}_{11}\left(\alpha \right)& ·& B_{1g}\left(\alpha \right)& ·& B_{1G}\left(\alpha \right)\\ ·& ·& ·& ·& ·\\ B_{g1}\left(\alpha \right)& ·& B_{gg}\left(\alpha \right)& ·& B_{gG}\left(\alpha \right)\\ ·& ·& ·& ·& ·\\ B_{G1}\left(\alpha \right)& ·& B_{Gg}\left(\alpha \right)& ·& B_{GG}\left(\alpha \right)\end{array}\right);\phi \left(r,\Omega \right)\triangleq \left(\begin{array}{c}{\phi }^1\\ ·\\ {\phi }^g\\ ·\\ {\phi }^G\end{array}\right);Q\left(\alpha \right)\triangleq \left(\begin{array}{c}{Q}^1\\ ·\\ Q^g\\ ·\\ Q^G\end{array}\right);$ (8)

with components defined below:

$\begin{array}{l}{B}_{gh}\left(\alpha \right)\triangleq {\delta }_{g,h}{B}_{gg}^0\left(\alpha \right)-B_{gh}^1\left(\alpha \right);g,h=1,\cdots ,G;\\ B_{gg}^0\left(\alpha \right){\phi }^g\left(r,\Omega \right)\triangleq \left[\Omega \cdot \nabla +{\Sigma }_t^g\left(\alpha ;r\right)\right]{\phi }^g\left(r,\Omega \right);\\ B_{gh}^1\left(\alpha \right){\phi }^h\left(r,{\Omega }^{\prime }\right)\triangleq {\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }\left[{\Sigma }_s^{h\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)+{\chi }^g\left(\alpha ;r\right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha ;r\right)\right].}\end{array}$ (9)

The notation ${\delta }_{g,h}$, which appears in Equation (9), denotes the Kronecker delta-functional, which is defined as usual, i.e., ${\delta }_{g,h}=1$, if $g=h$ and ${\delta }_{g,h}=0$, if $g\ne h$.

3. Computation of First-Order Sensitivities $\partial L\left(\alpha \right)/\partial {\alpha }_j,j=1,\cdots ,TP$

Basically, there are three methods for computing deterministically (as opposed to computing “statistically”) the response sensitivities to model and response parameters, as follows:

1) The so-called “brute-force re-computations” method, which uses finite-difference formulas to approximate the derivative that expresses the response sensitivity to the parameter under consideration; this method will be presented in Subsection 3.1;

2) The Forward Sensitivity Analysis Methodology (FSAM), which will be applied in Subsection 3.2;

3) The 1^st-Order Comprehensive Adjoint Sensitivity Analysis Methodology (CASAM) conceived by Cacuci [21], which will be applied in Subsection 3.3.

Subsection 3.4 will present a comparative discussion of the computational resources and times (CPU) which are required by the three above-mentioned deterministic methods.

3.1. Re-Computations with Finite-Difference Approximation

The first-order sensitivities of the leakage response, $L\left(\alpha \right)$, can be approximately computed by re-computations using the well-known finite-difference formula presented below:

$\frac{\partial L\left(\alpha \right)}{\partial {\alpha }_j}\approx \frac{1}{2h_j}\left(L_{j+1}-L_{j-1}\right)+{\rm O}\left(h_j^2\right),j=1,\cdots ,TP$, (10)

where $L_{j+1}\triangleq L\left({\alpha }_j+h_j\right)$, $L_{j-1}\triangleq L\left({\alpha }_j-h_j\right)$ and $h_j$ denotes a “judiciously-chosen” variation in the parameter ${\alpha }_j$ around its nominal value ${\alpha }_j^0$. The values $L_{j+1}\triangleq L\left({\alpha }_j+h_j\right)$ and $L_{j-1}\triangleq L\left({\alpha }_j-h_j\right)$ are obtained by re-solving Equations (1) and (2) repeatedly, using the changed parameter values $\left({\alpha }_j\pm h_j\right)$. As is well known, the value of the variation $h_j$ must be chosen by “trial and error” for each parameter ${\alpha }_j$, and varies from parameter to parameter; when $h_j$ is too large or too small, the result produced by Equation (10) will be far off from the exact value of the derivative $\partial L\left(\alpha \right)/\partial {\alpha }_j$. It is important to note that finite difference formulas introduce their intrinsic “methodological errors,” such as the error ${\rm O}\left(h_j^2\right)$ indicated in Equation (10), which are in addition to and independent of the errors that might be incurred in the computation of $L\left(\alpha \right)$. In other words, even if the computation of $L\left(\alpha \right)$ were perfect (error-free), the finite-difference formulas nevertheless introduce their own, intrinsic, numerical errors into the computation of the sensitivity $\text{d}L\left(\alpha \right)/\text{d}{\alpha }_j$.

3.2. Forward Sensitivity Analysis Methodology (FSAM)

If one could solve the neutron transport Equations (1) and (2) in closed form, the respective closed form of ${\phi }^g\left(r,\Omega \right)$ would be an explicit function of the parameters $\alpha$. However, since the neutron transport Equations (1) and (2) cannot be solved in closed form, it follows that ${\phi }^g\left(r,\Omega \right)$ is therefore an implicit function of $\alpha$. Consequently, as Equation (4) indicates, the PERP benchmark’s leakage response, $L\left(\alpha \right)$, is a function of the vector $\alpha$ of model parameters implicitly through the group-flux ${\phi }^g\left(r,\Omega \right)$.

The Forward Sensitivity Analysis Methodology (FSAM) is developed [18] [19] in the vector-space of the model parameters $\alpha$. In this vector-space, the first-order partial sensitivity of $L\left(\alpha \right)$ to a generic model parameter, ${\alpha }_j$, can in principle be obtained by taking the first-order partial G-derivative of Equation (4) with respect to a generic model/respond parameter ${\alpha }_j$, evaluated at the nominal parameter values ${\alpha }^0$, which yields:

$\begin{array}{c}{\left\{\frac{\partial L\left(\alpha \right)}{\partial {\alpha }_j}\right\}}_{{\alpha }^0}\triangleq \frac{\text{d}}{\text{d}\epsilon }{\left\{{\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n\frac{\partial {\phi }^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r,\Omega \right)}{\partial {\alpha }_j}\right\}}_{\epsilon =0}\\ ={\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n{\left\{\frac{\partial {\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_j}\right\}}_{{\alpha }^0},j=1,\cdots ,TP.\end{array}$ (11)

In turn, the derivative $\partial {\phi }^g\left(r,\Omega \right)/\partial {\alpha }_j$ is the solution of the equations obtained by taking the first-order G-derivatives of Equations (1) and (2) with respect to ${\alpha }_j$, evaluated at the nominal parameter values ${\alpha }^0$, which yields, by definition, the following equations:

$\begin{array}{l}\frac{\text{d}}{\text{d}\epsilon }\{\Omega \cdot \nabla {\phi }^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r,\Omega \right)+{\Sigma }_t^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r\right){\phi }^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r,\Omega \right)\\ -{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }{\phi }^h\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r,{\Omega }^{\prime }\right)[{\Sigma }_s^{h\to g}\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r,{\Omega }^{\prime }\to \Omega \right)}}\\ {+{\chi }^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r\right){\left(\nu {\Sigma }_f\right)}^h\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r\right)]-Q^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0\right)\}}_{\epsilon =0}=0,\\ j=1,\cdots ,TP;\end{array}$ (12)

$\frac{\text{d}}{\text{d}\epsilon }{\left\{{\phi }^g\left({\alpha }_j^0+\epsilon h_j;{\alpha }^0;r,\Omega \right)\right\}}_{\epsilon =0}=0,r=r_d,\Omega \cdot n<0,g=1,\cdots ,G;j=1,\cdots ,TP.$ (13)

Carrying out the differentiations with respect to $\epsilon$ in Equations (12) and (13), setting $\epsilon =0$ in the resulting equations and subsequently cancelling everywhere the scalar quantity $h_j$ yields the following equations for the function $\partial {\phi }^g\left(r,\Omega \right)/\partial {\alpha }_j$ :

$\begin{array}{l}{\left\{B^g\left(\alpha \right)\left[\frac{\partial {\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_j}\right]\right\}}_{{\alpha }^0}\\ ={\left\{\frac{\partial S^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_j}\right\}}_{{\alpha }^0},g=1,\cdots ,G;j=1,\cdots ,TP,\end{array}$ (14)

${\left\{\frac{\partial {\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_j}\right\}}_{{\alpha }^0}=0,r=r_d;\Omega \cdot n<0,g=1,\cdots ,G;j=1,\cdots ,TP$, (15)

where

$\frac{\partial S^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_j}\triangleq \frac{\partial Q^g\left(\alpha \right)}{\partial {\alpha }_j}-\frac{\partial B^g\left(\alpha \right)}{\partial {\alpha }_j}\left[{\phi }^g\left(r,\Omega \right)\right],$ (16)

$\begin{array}{l}\frac{\partial B^g\left(\alpha \right)}{\partial {\alpha }_j}\left[{\phi }^g\left(r,\Omega \right)\right]\triangleq \frac{\partial {\Sigma }_t^g\left(\alpha ;r\right)}{\partial {\alpha }_j}{\phi }^g\left(r,\Omega \right)\\ -{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }{\phi }^h\left(r,{\Omega }^{\prime }\right)\left\{\frac{\partial {\Sigma }_s^{h\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)}{\partial {\alpha }_j}+\frac{\partial \left[{\chi }^g\left(\alpha ;r\right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha ;r\right)\right]}{\partial {\alpha }_j}\right\}}}.\end{array}$ (17)

The system of Equations (14) and (15) is called the 1^st-Order Forward Sensitivity System (1^st-OFSS), the solutions of which are required in Equation (11) to compute, in the end, the response sensitivities using the forward sensitivity analysis method. Note that a “place-holder,” denoted as “[]”, has been explicitly used on the left-side of Equation (14) to indicate that the operator $B^g\left(\alpha \right)$ acts (linearly) on the function that appears in this place-holder. Place-holders “[]” have also been used on the right-side of Equation (16) and on the left-side of Equation (17) to indicate that the operator $\partial B^g\left(\alpha \right)/\partial {\alpha }_j$ acts on the function that appears in the respective place-holder. Since there are a total of TP parameters ${\alpha }_j$, it is evident that the computations of all of the sensitivities $\partial L\left(\alpha \right)/\partial {\alpha }_j$ require TP large-scale computations to solve the 1^st-OFSS, cf. Equations (14) and (15), followed by TP small-scale computations for performing the integration represented by Equation (11).

3.3. First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1st-CASAM)

The aim of the First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1^st-CASAM) is to find an alternative way for expressing the contribution of the flux variation (the so-called “indirect-effect term” contribution) to the total response sensitivity, so as to avoid the need for having to compute the derivatives $\partial {\phi }^g\left(r,\Omega \right)/\partial {\alpha }_j$ of the state functions with respect to the model’s parameters, as is the case when using the FSAM. In contradistinction to the FSAM, which is set in the vector-space of the parameters $\alpha$, the 1^st-CASAM is set in the combined space of the parameters and the state-functions, denoted as $u\triangleq \left(\phi ;\alpha \right)={\left[{\phi }^1,\cdots ,{\phi }^G;{\alpha }_1,\cdots ,{\alpha }_{TP}\right]}^{†}$, and does not a priori consider that the state functions $\phi \triangleq {\left[{\phi }^1,\cdots ,{\phi }^G\right]}^{†}$ are implicit functions of the parameters $\alpha$. In its full generality, the 1^st-CASAM also considers that the phase-space domain’s boundaries and internal interfaces are also subject to uncertainties (i.e., are imperfectly well known). For the purpose of this work, however, the boundaries of the phase-space domain (which consist of the boundaries of the spatial domain, i.e., the internal and external radii of the PERP sphere) and the end-points of the energy groups are considered to be free of errors (i.e., are perfectly well known) so that the sensitivities of the leakage response to these phase-space domain boundaries will not be considered. In this case, the 1^st-CASAM reduces to the adjoint sensitivity analysis methodology originally conceived by Cacuci [18] [19]. When the response sensitivities with respect to the boundaries of the phase-space domain are not considered, the G-differential, denoted as $\delta L\left(\alpha ;\delta {\phi }^g\right)$, of the leakage response defined in Equation (4) for a variation $\delta u\triangleq \left(\delta \phi ;\delta \alpha \right)$ around the nominal values $u^0\triangleq \left({\phi }^0;{\alpha }^0\right)$ is obtained, by definition, as follows:

$\begin{array}{c}\delta L\left(\alpha ;\delta \phi \right)\triangleq \frac{\text{d}}{\text{d}\epsilon }{\left\{{\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n\left[{\phi }^{g,0}\left(r,\Omega \right)+\epsilon \delta {\phi }^g\left(r,\Omega \right)\right]\right\}}_{\epsilon =0}\\ ={\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n\left[\delta {\phi }^g\left(r,\Omega \right)\right].\end{array}$ (18)

In general, the total G-variation of the response depends comprises two contributions, one arising directly from the parameter variations $\delta \alpha$ (this contribution is called the “direct-effect term”) and another contribution, stemming from the variations $\delta \phi$ in the state functions (this contribution is called the “indirect-effect term”). In the particular case of the leakage response considered in this work, the variation $\delta L\left(\alpha ;\delta {\phi }^g\right)$ does not depend directly on the parameter variations $\delta \alpha$ but depends only on variations $\delta \phi$ in the state functions, as indicated in Equation (18). Therefore, $\delta L\left(\alpha ;\delta {\phi }^g\right)$ coincides with the “indirect-effect term.”

The variations $\delta {\phi }^g\left(r,\Omega \right)$ are the solutions of the first-order G-differentials of the model Equations (1) and (2), which are determined by applying the definition of the G-differential to these equations, to obtain:

$\begin{array}{l}\frac{\text{d}}{\text{d}\epsilon }\{\Omega \cdot \nabla \left[{\phi }^{g,0}\left(r,\Omega \right)+\epsilon \delta {\phi }^g\left(r,\Omega \right)\right]\\ +{\Sigma }_t^g\left({\alpha }^0+\epsilon \delta \alpha ;r\right)\left[{\phi }^{g,0}\left(r,\Omega \right)+\epsilon \delta {\phi }^g\left(r,\Omega \right)\right]\\ -{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }\left[{\phi }^{h,0}\left(r,\Omega \right)+\epsilon \delta {\phi }^h\left(r,\Omega \right)\right][{\Sigma }_s^{h\to g}\left({\alpha }^0+\epsilon \delta \alpha ;r,{\Omega }^{\prime }\to \Omega \right)}}\\ +{{\chi }^g\left({\alpha }^0+\epsilon \delta \alpha ;r\right){\left(\nu {\Sigma }_f\right)}^h\left({\alpha }^0+\epsilon \delta \alpha ;r\right)]-Q^g\left({\alpha }^0+\epsilon \delta \alpha \right)\}}_{\epsilon =0}=0,\end{array}$ (19)

$\frac{\text{d}}{\text{d}\epsilon }{\left\{\left[{\phi }^{g,0}\left(r,\Omega \right)+\epsilon \delta {\phi }^g\left(r,\Omega \right)\right]\right\}}_{\epsilon =0}=0,r=r_d,\Omega \cdot n<0,g=1,\cdots ,G.$ (20)

Carrying out the differentiations with respect to $\epsilon$ in Equations (19) and (20), and setting $\epsilon =0$ in the resulting equations yields the following system of equations for the variations $\delta {\phi }^g\left(r,\Omega \right)$ in terms of the parameter variations $\delta \alpha$ :

$B^g\left({\alpha }^0\right)\left[\delta {\phi }^g\left(r,\Omega \right)\right]=Q^{\left(1\right)}\left(g;{\alpha }^0,{\phi }^0;\delta \alpha \right),g=1,\cdots ,G,$ (21)

$\delta {\phi }^g\left(r,\Omega \right)=0,r=r_d,\Omega \cdot n<0,g=1,\cdots ,G,$ (22)

where

$Q^{\left(1\right)}\left(g;{\alpha }^0,{\phi }^0;\delta \alpha \right)\triangleq {\left\{\frac{\partial Q^g\left(\alpha \right)}{\partial \alpha }\delta \alpha -\frac{\partial \left\{B^g\left(\alpha \right)\left[{\phi }^g\left(r,\Omega \right)\right]\right\}}{\partial \alpha }\delta \alpha \right\}}_{{\alpha }^0},$ (23)

$\begin{array}{l}\frac{\partial \left\{B^g\left(\alpha \right)\left[{\phi }^g\left(r,\Omega \right)\right]\right\}}{\partial \alpha }\delta \alpha \triangleq {\phi }^g\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial \alpha }\delta \alpha \\ -{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }{\phi }^h\left(r,{\Omega }^{\prime }\right)\left\{\frac{\partial {\Sigma }_s^{h\to g}\left(\alpha ;{\Omega }^{\prime }\to \Omega \right)}{\partial \alpha }\delta \alpha +\frac{\partial \left[{\chi }^g\left(\alpha \right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha \right)\right]}{\partial \alpha }\delta \alpha \right\}}}.\end{array}$ (24)

By analogy to Equations (6) and (7), Equations (21) and (22) can be written in following matrix-vector form:

$B\left(\alpha \right)\left[\delta \phi \left(r,\Omega \right)\right]=Q^{\left(1\right)}\left(\alpha ;\delta \alpha \right),$ (25)

$\delta \phi \left(r,\Omega \right)=0,r=r_d,\Omega \cdot n<0,$ (26)

where

$\delta \phi \left(r,\Omega \right)\triangleq \left(\begin{array}{c}\delta {\phi }^1\left(r,\Omega \right)\\ ·\\ \delta {\phi }^g\left(r,\Omega \right)\\ ·\\ \delta {\phi }^G\left(r,\Omega \right)\end{array}\right);Q^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\triangleq \left(\begin{array}{c}{Q}^{\left(1\right)}\left(1;\alpha ,\phi ;\delta \alpha \right)\\ ·\\ Q^{\left(1\right)}\left(g;\alpha ,\phi ;\delta \alpha \right)\\ ·\\ Q^{\left(1\right)}\left(G;\alpha ,\phi ;\delta \alpha \right)\end{array}\right).$ (27)

The aim of the 1^st-CASAM—to construct an alternative expression for the “indirect-effect term” $\delta L\left(\alpha ;\delta \phi \right)$ which would not depend on $\delta \phi$ —can be achieved by implementing the following sequence of steps, which are all to be carried out at the nominal parameter values ${\alpha }^0$ (although this will not be explicitly indicated in the sequence in order to keep the notation as simple as possible):

1) Consider that the square-integrable functions $\phi \in H_{\left(1\right)}$ and $\delta \phi \in H_{\left(1\right)}$ are elements of a Hilbert space denoted as $H_{\left(1\right)}$, in which the inner product between two elements $u\left(r,\Omega \right)\in H_{\left(1\right)}$ and $v\left(r,\Omega \right)\in H_{\left(1\right)}$ is denoted as ${\langle u\left(r,\Omega \right),v\left(r,\Omega \right)\rangle }_{\left(1\right)}$ and is defined as follows:

${\langle u\left(r,\Omega \right),v\left(r,\Omega \right)\rangle }_{\left(1\right)}\triangleq {\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{u}^g\left(r,\Omega \right)v^g\left(r,\Omega \right)$. (28)

2) Form the inner product of Equation (25) with a yet undefined function ${\psi }^{\left(1\right)}\left(r,\Omega \right)\triangleq {\left[{\psi }^{\left(1\right),g}\left(r,\Omega \right),\cdots ,{\psi }^{\left(1\right),G}\left(r,\Omega \right)\right]}^{†}\in H_{\left(1\right)}$ to obtain the following relation:

$\begin{array}{l}{\langle {\psi }^{\left(1\right)}\left(r,\Omega \right),B\left(\alpha \right)\left[\delta \phi \left(r,\Omega \right)\right]\rangle }_{\left(1\right)}={\langle {\psi }^{\left(1\right)}\left(r,\Omega \right),Q^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\rangle }_{\left(1\right)}\\ ={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{\psi }^{\left(1\right),g}\left(r,\Omega \right)Q^{\left(1\right)}\left(g;{\alpha }^0,{\phi }^0;\delta \alpha \right).\end{array}$ (29)

3) Using the inner product defined in Equation (28) and the definition of the adjoint operator in $H_{\left(1\right)}$, recast the left-side of Equation (29) into the following equivalent form:

$\begin{array}{l}{\langle {\psi }^{\left(1\right)}\left(r,\Omega \right),B\left(\alpha \right)\left[\delta \phi \left(r,\Omega \right)\right]\rangle }_{\left(1\right)}\\ ={\langle \delta \phi \left(r,\Omega \right),A\left(\alpha \right){\psi }^{\left(1\right)}\left(r,\Omega \right)\rangle }_{\left(1\right)}+P\left[\delta \phi ,{\psi }^{\left(1\right)}\right],\end{array}$ (30)

where

$\begin{array}{l}P\left[\delta \phi ,{\psi }^{\left(1\right)}\right]\triangleq {\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n<0}{\int }\text{d}\Omega {\displaystyle \underset{\partial V}{\int }\left|\Omega \cdot n\right|\left[\delta {\phi }^g\left(r,\Omega \right){\psi }^{\left(1\right),g}\left(r,\Omega \right)\right]}}}\\ -{\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega {\displaystyle \underset{\partial V}{\int }\Omega \cdot n\left[\delta {\phi }^g\left(r,\Omega \right){\psi }^{\left(1\right),g}\left(r,\Omega \right)\right]}}},\end{array}$ (31)

denotes the bilinear concomitant evaluated on the PERP sphere’s surface $S_b$ and where $A\left(\alpha \right)$ denotes the operator adjoint to $B\left(\alpha \right)$, having components $A_{gh}\left(\alpha \right)\triangleq {\left[B_{hg}\left(\alpha \right)\right]}^{*}$, where the symbol ${\left[\right]}^{*}$ indicates “formal adjoint operator.”

4) Require the first term on the right-side of Equation (30) to represent the functional defined in Equation (18) to obtain the following relation:

$A\left(\alpha \right){\psi }^{\left(1\right)}\left(r,\Omega \right)=I\left[\Omega \cdot n\delta \left(r-r_d\right)\right];\text{}\text{}\text{}I\triangleq {\left[1,\cdots ,1,\cdots ,1\right]}^{†}$. (32)

5) Use the boundary condition given in Equation (22) to eliminate the corresponding boundary terms in the bilinear concomitant $P\left[\delta \phi ,{\psi }^{\left(1\right)}\right]$ in Equation (31) and eliminate the appearance of the remaining boundary terms in Equation (31) by imposing the following (adjoint) boundary condition:

${\psi }^{\left(1\right)}\left(r,\Omega \right)=0,r=r_d,\Omega \cdot n>0.$ (33)

In component form, Equations (32) and (33) are as follows:

$\begin{array}{l}-\Omega \cdot \nabla {\psi }^{\left(1\right),g}\left(r,\Omega \right)+{\Sigma }_t^g\left(\alpha \right){\psi }^{\left(1\right),g}\left(r,\Omega \right)\\ -{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }}{\psi }^{\left(1\right),h}\left(r,{\Omega }^{\prime }\right)}\left[{\Sigma }_s^{g\to h}\left(\alpha ;\Omega \to {\Omega }^{\prime }\right)+{\left(\nu {\Sigma }_f\right)}^g\left(\alpha \right){\chi }^h\left(\alpha \right)\right]\\ =\Omega \cdot n\delta \left(r-r_d\right),g=1,\cdots ,G;\end{array}$ (34)

${\psi }^{\left(1\right),g}\left(r,\Omega \right)=0,r=r_d,\Omega \cdot n>0,g=1,\cdots ,G.$ (35)

It is important to note that the relations presented in Equations (28) through (35) hold generally, including at the nominal values $u^0\triangleq \left({\phi }^0;{\alpha }^0\right)$. In order to simplify the notation, the superscript “zero” denoting nominal values has been omitted in the foregoing derivations and will occasionally be omitted henceforth, as well. This simplification should not cause any loss of clarity since it will become clear from the context which quantities are to be evaluated/computed using the nominal values for the model parameters.

6) Use Equations (29) through (33) in Equation (18) to obtain the following expression for the “indirect-effect term”:

$\begin{array}{c}\delta L\left(\alpha ;\delta \phi \right)={\langle {\psi }^{\left(1\right)}\left(r,\Omega \right),Q^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\rangle }_{\left(1\right)}\\ ={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{\psi }^{\left(1\right),g}\left(r,\Omega \right)Q^{\left(1\right)}\left(g;{\alpha }^0,{\phi }^0;\delta \alpha \right)\\ \triangleq {\displaystyle \underset{j=1}{\overset{TP}{\sum }}\frac{\partial L\left(\alpha \right)}{\partial {\alpha }_j}\delta {\alpha }_j}.\end{array}$ (36)

7) Replacing in Equation (36) the expression of $Q^{\left(1\right)}\left(g;{\alpha }^0,{\phi }^0;\delta \alpha \right)$ from Equation (23) and identifying the expressions that multiply the arbitrary parameter variations ${\alpha }_j$ yields the expression of each partial sensitivity $\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)/\partial {\alpha }_j$, for all ${\alpha }_j$, $j=1,\cdots ,TP$, as follows:

$\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial q_j}={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right)\frac{\partial Q^g\left(\alpha \right)}{\partial q_j}}\text{}\text{}\text{},j=1,\cdots ,JQ=10;$ (37)

$\begin{array}{l}\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial N_j}={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}\Omega }{\psi }^{(1),g}\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial N_j}{\phi }^g\left(r,\Omega \right)}\text{}\\ -{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }{\phi }^h\left(r,{\Omega }^{\prime }\right)\left\{\frac{\partial {\Sigma }_s^{h\to g}\left(\alpha ;{\Omega }^{\prime }\to \Omega \right)}{\partial N_j}+\frac{\partial \left[{\chi }^g\left(\alpha \right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha \right)\right]}{\partial N_j}\right\}}}\text{},\\ \text{for}j=1,\cdots ,I=6;\end{array}$ (38)

$\begin{array}{l}\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial s_j}={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right)}\\ \times {\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }}\frac{\partial {\Sigma }_s^{h\to g}\left(\alpha ;{\Omega }^{\prime }\to \Omega \right)}{\partial s_j}{\phi }^h\left(r,{\Omega }^{\prime }\right)},\\ \text{for}j=1,\cdots ,JSX=\left(G\times G\right)\times I\times \left(ISCT+1\right)=21600;\end{array}$ (39)

$\begin{array}{l}\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial t_j}=-{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right){\phi }^g\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_j}},\\ \text{for}j=1,\cdots ,JTX=180;\end{array}$ (40)

$\begin{array}{l}\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial f_j}={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right){\chi }^g}\\ \times {\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}{\Omega }^{\prime }}\frac{\partial \left[{\left(\nu {\Sigma }_f\right)}^h\right]}{\partial f_j}{\phi }^h\left(r,{\Omega }^{\prime }\right)},\text{for}j=1,\cdots ,JFX=60;\end{array}$ (41)

$\begin{array}{l}\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial {\nu }_j}={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right){\chi }^g}\\ \times {\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}{\Omega }^{\prime }}\frac{\partial \left[{\left(\nu {\Sigma }_f\right)}^h\right]}{\partial {\nu }_j}{\phi }^h\left(r,{\Omega }^{\prime }\right),}\text{for}j=1,\cdots ,JNU=60;\end{array}$ (42)

$\begin{array}{l}\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial p_j}={\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right)\frac{\partial {\chi }^g\left(\alpha \right)}{\partial p_j}}\\ \times {\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle {\int }_{4\pi }\text{d}{\Omega }^{\prime }}{\left(\nu {\Sigma }_f\right)}^h{\phi }^h\left(r,{\Omega }^{\prime }\right)},\text{for}j=1,\cdots ,J\chi =60.\end{array}$ (43)

The expressions of the sensitivities provided in Equations (37) through (43) are to be evaluated at the nominal parameter values ${\alpha }^0$ (which has not been indicted explicitly, in order to simplify the respective notations). All of these sensitivities have been computed numerically in [5] - [10]. Evidently, the computations of these sensitivities are inexpensive, involving only integrations using quadrature formulas, once the (1^st-level adjoint) function ${\psi }^{\left(1\right)}\left(r,\Omega \right)$ is available.

8) Equations (32) and (33) are called the First-Level Adjoint Sensitivity System (1^st-LASS); the solution ${\psi }^{\left(1\right)}\left(r,\Omega \right)$ of the 1^st-LASS is called the 1^st-level adjoint function. The system of equations is called “first-level” as opposed to “first-order” because it does not contain any derivatives (of first- or higher-order) but its solution is used for computing the first-order sensitivities. This terminology will be also used in the sequel, when deriving the expressions for the 2^nd- and 3^rd-order sensitivities. The 1^st-LASS is independent of parameter variations, so it needs to be solved just once to obtain ${\psi }^{\left(1\right)}\left(r,\Omega \right)$. Thus, using the 1^st-CASAM, the computation of the TP partial sensitivities $\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)/\partial {\alpha }_j$, $j=1,\cdots ,TP$, using Equations (37) through requires just a single large-scale computation in order to determine ${\psi }^{\left(1\right)}\left(r,\Omega \right)$, followed by TP inexpensive computations to perform each integration (quadrature) involving the 1^st-level adjoint function ${\psi }^{\left(1\right)}\left(r,\Omega \right)$.

3.4. Comparison of Computational Requirements

The computer used for the PERP computations is a DELL desktop computer (AMD FX-8350) with an 8-core processor. A forward computation of the flux distribution $\phi \left(r,\Omega \right)$ within the PERP benchmark and of the leakage response takes ca. 2 minutes CPU-time when the forward transport equation is solved using an angular expansion of order ISN = 256 (S₂₅₆).

The 1^st-level adjoint function ${\psi }^{\left(1\right)}\left(r,\Omega \right)$ used to compute the first-order sensitivities given by Equations (37) through (43) has been computed in [5] by solving the 1^st-LASS using an angular expansion of order ISN = 256 (S₂₅₆). The CPU-time for a typical adjoint computation of ${\psi }^{\left(1\right)}\left(r,\Omega \right)$ using S₂₅₆ is ca. 85 seconds, which is comparable to the CPU-time needed for solving the forward transport equation to compute the flux distribution $\phi \left(r,\Omega \right)$. Furthermore, the typical CPU-time needed to perform the integration (quadrature) over the forward and adjoint functions, to compute the respective sensitivity (using S₂₅₆) is ca. 0.0012 seconds. The CPU requirements needed for computing the 7477 non-zero sensitivities using forward computations in conjunction with the 2-point finite difference scheme presented in Section 3.1 and using the FSAM presented in Section 3.2, respectively, are compared in Table 1 with the CPU-requirements needed to compute the same sensitivities using the formulas obtained

Table 1. CPU-time needed for computing $\partial L\left(\alpha \right)/\partial {\alpha }_j,j=1,\cdots ,TP=7477$.

in Equations (37) through (43) in Section 3.3, which were obtained by having applied the 1^st-CASAM. The “95 seconds CPU-time” needed by the 1^st-CASAM comprises 85 seconds for the one adjoint computation plus 10 seconds for computing all of the 7477 1^st-order sensitivities.

The results presented in Table 1 highlights the evident superiority of the 1^st-CASAM over the finite-difference method and/or the Forward Sensitivity Analysis Methodology (FSAM) for computing the sensitivities of one response (i.e., the PERP leakage) to the 7477 PERP parameters that have non-zero nominal values.

4. Computation of ${\partial }^{2}L\left(\alpha \right)/\partial t_{j1}\partial t_{j2},j1=1,\cdots ,JTX;j2=1,\cdots ,j1$

The finite-difference formulas that can be used to compute approximately the second-order response sensitivities with respect to the model parameters will be discussed in Subsection 4.1. The Forward Sensitivity Analysis Methodology will be presented in Subsection 4.2. The application of the 2^nd-CASAM will be presented in Subsection 4.3, while Subsection 4.4 will compare the CPU-times which are required by these three deterministic methods.

4.1. Re-computations with Finite-Difference Approximation

The 2^nd-order unmixed sensitivities, ${{\partial }^{2}L\left(\alpha \right)/\left(\partial {\alpha }_j\right)}^2$, of the leakage response, $L\left(\alpha \right)$, can be approximately computed by re-computations using the well-known finite-difference formula presented below:

$\frac{{\partial }^{2}L\left(\alpha \right)}{{\left(\partial {\alpha }_j\right)}^2}\approx \frac{1}{h_j^2}\left(L_{j+1}-L_j+L_{j-1}\right)+{\rm O}\left(h_j^2\right),j=1,\cdots ,TP$, (44)

where $L_{j+1}\triangleq L\left({\alpha }_j+h_j\right)$, $L_j\triangleq L\left({\alpha }_j\right)$, $L_{j-1}\triangleq L\left({\alpha }_j-h_j\right)$ and $h_j$ denotes a “judiciously-chosen” variation in the parameter ${\alpha }_j$ around its nominal value ${\alpha }_j^0$. The 2^nd-order mixed sensitivities, ${\partial }^{2}L\left(\alpha \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}$, can be calculated by using the following finite-difference formula:

$\begin{array}{l}\frac{{\partial }^{2}L\left(\alpha \right)}{\partial {\alpha }_{j1}\partial {\alpha }_{j2}}\approx \frac{1}{4h_{j1}{h}_{j2}}\left(L_{j1+1,j2+1}-L_{j1-1,j2+1}-L_{j1+1,j2-1}+L_{j1-1,j2-1}\right)+{\rm O}\left(h_{j1}^2,h_{j2}^2\right),\\ j1=1,\cdots ,TP;j2=1,\cdots ,j1,\end{array}$ (45)

where $h_{j1}$ and $h_{j2}$ denote the variations in the parameters ${\alpha }_{j1}$ and ${\alpha }_{j2}$, respectively. The values of the quantities $L_{j1+1,j2+1}\triangleq L\left({\alpha }_{j1}+h_{j1},{\alpha }_{j2}+h_{j2}\right)$, $L_{j1-1,j2+1}\triangleq L\left({\alpha }_{j1}-h_{j1},{\alpha }_{j2}+h_{j2}\right)$, $L_{j1+1,j2-1}\triangleq L\left({\alpha }_{j1}+h_{j1},{\alpha }_{j2}-h_{j2}\right)$, and $L_{j1-1,j2-1}\triangleq L\left({\alpha }_{j1}-h_{j1},{\alpha }_{j2}-h_{j2}\right)$ are obtained by re-solving Equations (1) and (2) repeatedly, using the changed parameter values $\left({\alpha }_{j1}\pm h_{j1}\right)$ and $\left({\alpha }_{j2}\pm h_{j2}\right)$. As also discussed in Subsection 3.1, the values of $h_{j1}$ and $h_{j2}$, respectively, must be chosen “judiciously” by “trial and error” for each of the parameters ${\alpha }_{j1}$ and ${\alpha }_{j2}$; otherwise, the result produced by Equation (45) will be far off from the exact value of the derivative ${\partial }^{2}L\left(\alpha \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}$. These simple finite difference formulas introduce their intrinsic “methodological errors” of order ${\rm O}\left(h_{j1}^2\right)$ and/or ${\rm O}\left(h_{j1}^2,h_{j2}^2\right)$, as indicated in Equations (44) and (45), which are in addition to,and independent of, the errors that might be incurred in the computation of $L\left(\alpha \right)$. In other words, even if the computation of $L\left(\alpha \right)$ were perfect (error-free), the finite-difference formulas nevertheless introduce their own, intrinsic, numerical errors into the computation of ${\partial }^{2}L\left(\alpha \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}$.

4.2. Forward Sensitivity Analysis Methodology

The second-order partial sensitivity of $L\left(\alpha \right)$ to two generic model parameters, ${\alpha }_{j1}$ and ${\alpha }_{j2}$ is obtained by differentiating the expression provided in Equation (11) with respect to (the second parameter) ${\alpha }_{j2}$, to obtain the following expression:

$\begin{array}{l}{\left\{\frac{{\partial }^{2}L\left(\alpha \right)}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0}={\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n{\left\{\frac{{\partial }^2{\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j1}\partial {\alpha }_{j2}}\right\}}_{{\alpha }^0},\\ j1=1,\cdots ,TP;j2=1,\cdots ,j1.\end{array}$ (46)

The second-order derivative ${\partial }^2{\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}$ is the solution of the second-order G-derivative of Equations (1) and (2) with respect to ${\alpha }_{j1}$ and ${\alpha }_{j2}$, namely:

$\begin{array}{l}{\left\{B^g\left(\alpha \right)\left[\frac{{\partial }^2{\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right]\right\}}_{{\alpha }^0}={\left\{\frac{{\partial }^2{S}^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0},\\ \text{for}g=1,\cdots ,G;j1=1,\cdots ,TP;j2=1,\cdots ,j1;\end{array}$ (47)

${\left\{\frac{{\partial }^2{\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0}=0,r=r_d;\Omega \cdot n<0,g=1,\cdots ,G;j1=1,\cdots ,TP;j2=1,\cdots ,j1,$ (48)

where

$\begin{array}{l}\frac{{\partial }^2{S}^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}=\frac{\partial }{\partial {\alpha }_{j2}}\left\{\frac{\partial S^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_{j1}}\right\}\\ =\frac{{\partial }^2{Q}^g\left(r\right)}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}-\frac{{\partial }^2{\Sigma }_t^g\left(\alpha ;r\right)}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}{\phi }^g\left(r,\Omega \right)-\frac{\partial {\Sigma }_t^g\left(\alpha ;r\right)}{\partial {\alpha }_{j1}}\frac{\partial {\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j2}}\\ +{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }\frac{\partial {\phi }^h\left(r,{\Omega }^{\prime }\right)}{\partial {\alpha }_{j2}}\left\{\frac{\partial {\Sigma }_s^{h\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)}{\partial {\alpha }_{j1}}+\frac{\partial \left[{\chi }^g\left(\alpha ;r\right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha ;r\right)\right]}{\partial {\alpha }_{j1}}\right\}}}\\ +{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }{\phi }^{g^{\prime }}\left(r,{\Omega }^{\prime }\right)\left\{\frac{{\partial }^2{\Sigma }_s^{h\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}+\frac{{\partial }^2\left[{\chi }^g\left(\alpha ;r\right){\left(\nu {\Sigma }_f\right)}^h\left(\alpha ;r\right)\right]}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}}.\end{array}$ (49)

It is important to note that, as indicated in Equation (49), the computation of the 2^nd-order flux-derivatives ${\partial }^2{\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j2}\partial {\alpha }_{j1}$ requires knowledge of the 1^st-order flux derivatives $\partial {\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j1}$, which appear in the definition of the right-side (“source”) for Equation (47) and which must therefore by computed prior to commencing the computations for obtaining ${\partial }^2{\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j2}\partial {\alpha }_{j1}$. Thus, the computation of the $TP\left(TP+1\right)/2$ distinct 2^nd-order sensitivities ${\partial }^{2}L\left(\alpha \right)/\partial {\alpha }_{j2}\partial {\alpha }_{j1}$ requires TP large-scale computations for solving Equations (14) and (15) to determine the fluxes $\partial {\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j1}$ which are needed to compute the right-side of Equation (47), followed by at least $TP\left(TP+1\right)/2$ large-scale computations to solve Equations (47) and (48) to determine the quantities ${\partial }^2{\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j2}\partial {\alpha }_{j1}$. Subsequently, $TP\left(TP+1\right)/2$ small-scale computations are needed for performing the integrations represented by Equation (46). To summarize: obtaining the numerical values of the 2^nd-order sensitivities ${\partial }^{2}L\left(\alpha \right)/\partial {\alpha }_i\partial {\alpha }_j$ using the FSAM requires a total of $TP+TP\left(TP+1\right)/2$ large-scale computations followed by $TP\left(TP+1\right)/2$ small-scale computations.

4.3. Second-Order Adjoint Sensitivity Analysis Methodology (2^nd-CASAM)

The numerical results for the complete 2^nd-order sensitivity analysis of the PERP benchmark’s leakage response with respect to all $TP\left(TP+1\right)/2$ nuclear parameters ${\alpha }_j,j=1,\cdots ,TP$, have been reported in [5] - [10]. Since this work ultimately aims at applying the 4^th-Order Comprehensive Sensitivity Analysis Methodology (4^th-CASAM) to obtain the closed form expressions of the 4^th-order sensitivities of the PERP leakage response with respect to the group-averaged total microscopic cross sections, the proliferation of indices, superscripts and subscripts is unavoidable. Nevertheless, “subscripted-subscripts” can be avoided by using subscripts of the form $j1=1,\cdots ,JTX;j2=1,\cdots ,j1$, where the index j1 will replace the index j (which was used in the Section 3) to index the 1^st-order sensitivities, and where the index j2 will be used (in addition to j1) to index the 2^nd-order sensitivities. Furthermore, anticipating the notation to be used in subsequent Subsections, the index j3 will be used (in addition to j1 and j2) to index the 3^rd-order sensitivities; the index j4 will be used (in addition to j1, j2 and j3) to index the 4^th-order sensitivities. Recall that even when only the components ${\sigma }_t\triangleq {\left[t_1,\cdots ,t_{JTX}\right]}^{†}$, $i=1,\cdots ,I=6$, $g=1,\cdots ,G=30$, $JTX=I\times G=180$, of the group-averaged total microscopic cross sections are considered, the numbers of sensitivities of the PERP leakage response with respect to these parameters are as follows: 1) 180 first-order sensitivities; 2) 32,400 second-order sensitivities, of which 16,290 are distinct; and 3) 5,832,000 third-order sensitivities, of which 988,260 are distinct; 4) 1,049,760,000 fourth-order sensitivities, of which 45,212,895 are distinct.

In contradistinction to the FSAM, the 2^nd-CASAM does not use Equation (11) as the starting point for computing 2^nd-order response sensitivities but instead uses the expressions of the sensitivities in terms of the 1^st-level adjoint function ${\psi }^{\left(1\right)}$, as produced by the 1^st-CASAM, namely Equations (37) through (43), as individual starting points for deriving the expressions of the corresponding 2^nd-order sensitivities. This procedure has been illustrated in [5] - [10] which has derived and computed all of the 2^nd-order sensitivities of the PERP leakage response to the nuclear data. In turn, the 3^rd-CASAM uses the expressions of the 2^nd-order response sensitivities derived within the 2^nd-CASAM. The 4^th-CASAM will use the expressions of the 3^rd-order response sensitivities derived within the 3^rd-CASAM, and so on.

In particular, the 2^nd-order sensitivities involving only the total microscopic cross sections, are obtained by G-differentiating the expression provided in Equation (40) for the 1^st-order sensitivities of the leakage response with respect to the group-averaged total microscopic cross sections, which yields the following result:

$\begin{array}{l}\delta {\left\{\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial t_{j1}}\right\}}_{{\alpha }^0}\\ =\frac{\text{d}}{\text{d}\epsilon }\{-{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }{\left\{\left[{\psi }^{\left(1\right),g}\left(r,\Omega \right)+\epsilon \delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)\right]{\phi }^g\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j1}}\right\}}_{{\alpha }^0}}\\ -{\left\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right)\left[{\phi }^g\left(r,\Omega \right)+\epsilon \delta {\phi }^g\left(r,\Omega \right)\right]\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j1}}}\right\}}_{{\alpha }^0}\\ {-{\left\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }{\psi }^{\left(1\right),g}\left(r,\Omega \right){\phi }^g\left(r,\Omega \right)\frac{\partial }{\partial t_{j1}}\left[{\Sigma }_t^g\left(\alpha \right)+\epsilon \delta {\Sigma }_t^g\left(\alpha \right)\right]}\right\}}_{{\alpha }^0}\}}_{\epsilon =0},\\ \text{for}j1=1,\cdots ,JTX=180.\end{array}$ (50)

Carrying out the differentiations with respect to $\epsilon$ in Equation (50) and subsequently setting $\epsilon =0$ in the resulting equations yields the following relation:

$\begin{array}{l}\delta {\left\{\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)}{\partial t_{j1}}\right\}}_{{\alpha }^0}\\ =-{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right){\left\{{\phi }^g\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j1}}\right\}}_{{\alpha }^0}}\\ -{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }\delta {\phi }^g\left(r,\Omega \right)}{\left\{{\psi }^{\left(1\right),g}\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j1}}\right\}}_{{\alpha }^0}\\ -{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }\frac{\partial }{\partial t_{j1}}\delta {\Sigma }_t^g\left(\alpha \right)}{\left\{{\psi }^{\left(1\right),g}\left(r,\Omega \right){\phi }^g\left(r,\Omega \right)\right\}}_{{\alpha }^0},\\ j1=1,\cdots ,JTX=180.\end{array}$ (51)

Taking into account the facts that: 1) the parameter $t_{j1}$ represents the group-averaged total microscopic cross sections ${\sigma }_{t,i}^g$, with $j1=\left(i-1\right)G+g$,

$i=1,\cdots ,I=6$, $g=1,\cdots ,G=30$; and 2) ${\Sigma }_t^g\left(t\right)={\displaystyle \underset{i=1}{\overset{I=6}{\sum }}{N}_i{\sigma }_{t,i}^g}$, it follows that the last term on the left-side of Equation (51) is identically zero.

In Equation (50), the variations $\delta {\phi }^g\left(r,\Omega \right)$ are the solutions of Equations (21) and (22) while the variations $\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)$ are the solutions of the G-differentiated 1^st-LASS, namely:

${\left\{A^g\left(\alpha \right)\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)\right\}}_{{\alpha }^0}={\left\{Q^{\left(2\right)}\left(g;\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\right\}}_{{\alpha }^0},g=1,\cdots ,G,$ (52)

${\left\{\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)\right\}}_{{\alpha }^0}=0,r=r_d,\Omega \cdot n>0,g=1,\cdots ,G,$ (53)

where

$\begin{array}{l}{Q}^{\left(2\right)}\left(g;\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\triangleq -\frac{\partial \left\{A^g\left(\alpha \right)\left[{\psi }^{\left(1\right),g}\left(r,\Omega \right)\right]\right\}}{\partial \alpha }\delta \alpha \\ =-{\psi }^{\left(1\right),g}\left(r,\Omega \right)\frac{\partial {\Sigma }_t^g\left(\alpha ;r\right)}{\partial \alpha }\delta \alpha \\ +{\displaystyle \underset{h=1}{\overset{G}{\sum }}{\displaystyle \underset{4\pi }{\int }\text{d}{\Omega }^{\prime }}{\psi }^{\left(1\right),h}\left(r,{\Omega }^{\prime }\right)\left\{\frac{\partial {\Sigma }_s^{g\to h}\left(\alpha ;\Omega \to {\Omega }^{\prime }\right)}{\partial \alpha }\delta \alpha +\frac{\partial \left[{\left(\nu {\Sigma }_f\right)}^g\left(\alpha \right){\chi }^h\left(\alpha \right)\right]}{\partial \alpha }\delta \alpha \right\},}\\ \text{for}g=1,\cdots ,G.\end{array}$ (54)

In matrix-vector form, which will be used in the sequel, Equations (52) and (53) are expressed as follows:

${\left\{A\left(\alpha \right)\left[\delta {\psi }^{\left(1\right)}\left(r,\Omega \right)\right]\right\}}_{{\alpha }^0}={\left\{Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\right\}}_{{\alpha }^0},$ (55)

$\delta {\psi }^{\left(1\right)}\left(r,\Omega \right)=0,r=r_d,\Omega \cdot n>0.$ (56)

where

$\delta {\psi }^{\left(1\right)}\left(r,\Omega \right)\triangleq \left(\begin{array}{c}\delta {\psi }^{\left(1\right),1}\left(r,\Omega \right)\\ ·\\ \delta {\psi }^{\left(1\right),G}\left(r,\Omega \right)\end{array}\right),Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\triangleq \left(\begin{array}{c}{Q}^{\left(2\right)}\left(1;\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\\ ·\\ Q^{\left(2\right)}\left(G;\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\end{array}\right).$ (57)

The 2^nd-order sensitivities expressed by Equation (51) could be computed, in principle, if the functions $\delta {\phi }^g\left(r,\Omega \right)$ and $\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)$ were available. However, determining the functions $\delta {\phi }^g\left(r,\Omega \right)$ would require solving Equations (21) and (22), while determining the functions $\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)$ would require solving Equations (52) and (53), which would altogether amount to twice as many large-scale computations as there are parameter variations. To avoid the need for performing so many large-scale computations, the right-side of Equation (50) will be recast in terms of the solutions of a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) so as to eliminate the appearance of the functions $\delta {\phi }^g\left(r,\Omega \right)$ and $\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)$. The requisite 2^nd-LASS will be constructed by following the general principles introduced by Cacuci [1], which comprise the following sequence of steps, all of which are to be carried out at the nominal parameter values ${\alpha }^0$ :

1) Consider a Hilbert space, denoted as $H_{\left(2\right)}$, comprising block-vector elements $u^{\left(2\right)}\left(r,\Omega \right)\in H_{\left(2\right)}$ having the structure: $u^{\left(2\right)}\left(r,\Omega \right)\triangleq {\left[u_1^{\left(2\right)}\left(r,\Omega \right),u_2^{\left(2\right)}\left(r,\Omega \right)\right]}^{†}$, where $u_i^{\left(2\right)}\left(r,\Omega \right)\triangleq {\left[u_i^{\left(2\right),1}\left(r,\Omega \right),\cdots ,u_i^{\left(2\right),g}\left(r,\Omega \right),\cdots ,u_i^{\left(2\right),G}\left(r,\Omega \right)\right]}^{†},i=1,2$;

2) Define the inner product in $H_{\left(2\right)}$, denoted as ${\langle u^{\left(2\right)}\left(r,\Omega \right),v^{\left(2\right)}\left(r,\Omega \right)\rangle }_{\left(2\right)}$, of two arbitrary elements in $H_{\left(2\right)}$, denoted as $u^{\left(2\right)}\left(r,\Omega \right)\in H_{\left(2\right)}$ and $v^{\left(2\right)}\left(r,\Omega \right)\in H_{\left(2\right)}$, with $v^{\left(2\right)}\left(r,\Omega \right)\triangleq {\left[v_1^{\left(2\right)}\left(r,\Omega \right),v_2^{\left(2\right)}\left(r,\Omega \right)\right]}^{†}$, and $v_i^{\left(2\right)}\left(r,\Omega \right)\triangleq {\left[v_i^{\left(2\right),1}\left(r,\Omega \right),\cdots ,v_i^{\left(2\right),g}\left(r,\Omega \right),\cdots ,v_i^{\left(2\right),G}\left(r,\Omega \right)\right]}^{†},i=1,2$, as follows:

$\begin{array}{l}{\langle u^{\left(2\right)}\left(r,\Omega \right),v^{\left(2\right)}\left(r,\Omega \right)\rangle }_{\left(2\right)}\triangleq {\displaystyle \underset{i=1}{\overset{2}{\sum }}{\langle u_i^{\left(2\right)}\left(r,\Omega \right),v_i^{\left(2\right)}\left(r,\Omega \right)\rangle }_{\left(1\right)}}\\ ={\displaystyle \underset{i=1}{\overset{2}{\sum }}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{u}_i^{\left(2\right),g}\left(r,\Omega \right)}{v}_i^{\left(2\right),g}\left(r,\Omega \right).\end{array}$ (58)

3) For a 2 × 2 block-matrix linear operator $L^{\left(2\right)}\triangleq \left(\begin{array}{cc}{L}_{11}^{\left(2\right)}& L_{12}^{\left(2\right)}\\ L_{21}^{\left(2\right)}& L_{22}^{\left(2\right)}\end{array}\right)$ in the Hilbert space $H_{\left(2\right)}$, use the inner product defined in Equation (58) to define the formal adjoint operator, denoted as $A^{\left(2\right)}\triangleq \left(\begin{array}{cc}{A}_{11}^{\left(2\right)}& A_{12}^{\left(2\right)}\\ A_{21}^{\left(2\right)}& A_{22}^{\left(2\right)}\end{array}\right)$, of $L^{\left(2\right)}$, through the following relationship:

${\langle v^{\left(2\right)},\text{}{L}^{\left(2\right)}{u}^{\left(2\right)}\rangle }_{\left(2\right)}={\langle u^{\left(2\right)},A^{\left(2\right)}{v}^{\left(2\right)}\rangle }_{\left(2\right)}+P^{\left(2\right)}\left[u^{\left(2\right)},v^{\left(2\right)}\right]$, (59)

where $P^{\left(2\right)}\left[u^{\left(2\right)},v^{\left(2\right)}\right]$ denotes the corresponding bilinear concomitant evaluated on the domain’s boundary.

4) Apply the definition provided in Equation (58) to form the inner product of Equations (25) and (55) with a yet undefined function ${\psi }^{\left(2\right)}\left(j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]}^{†}\in H_{\left(2\right)}$, to obtain the following relation, evaluated at ${\alpha }^0$ :

$\begin{array}{l}{\left\{{\langle {\psi }^{\left(2\right)}\left(j1;r,\Omega \right),\left(\begin{array}{c}B\left(\alpha \right)\left[\delta \phi \left(r,\Omega \right)\right]\\ A\left(\alpha \right)\left[\delta {\psi }^{\left(1\right)}\left(r,\Omega \right)\right]\end{array}\right)\rangle }_{\left(2\right)}\right\}}_{{\alpha }^0}\\ ={\left\{{\langle {\psi }^{\left(2\right)}\left(j1;r,\Omega \right),\left(\begin{array}{c}{Q}^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\\ Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\end{array}\right)\rangle }_{\left(2\right)}\right\}}_{{\alpha }^0}.\end{array}$ (60)

5) Use the relation shown in Equation (59) to recast the left side of Equation (60) into the following equivalent form:

$\begin{array}{l}{\left\{{\langle {\psi }^{\left(2\right)}\left(j1;r,\Omega \right),\left(\begin{array}{c}B\left(\alpha \right)\left[\delta \phi \left(r,\Omega \right)\right]\\ A\left(\alpha \right)\left[\delta {\psi }^{\left(1\right)}\left(r,\Omega \right)\right]\end{array}\right)\rangle }_{\left(2\right)}\right\}}_{{\alpha }^0}\\ ={\left\{{\langle \left[\delta \phi ,\delta {\psi }^{\left(1\right)}\right],\left(\begin{array}{c}A\left(\alpha \right){\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\\ B\left(\alpha \right){\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\end{array}\right)\rangle }_{\left(2\right)}\right\}}_{{\alpha }^0}\\ +P^{\left(2\right)}\left[\delta \phi \left(r,\Omega \right),\delta {\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right],\end{array}$ (61)

where $P^{\left(2\right)}\left[\delta \phi \left(r,\Omega \right),\delta {\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]$ denotes the corresponding bilinear concomitant evaluated on the PERP sphere’s outer boundary $S_b$, at $r=r_d$.

6) Require the first term on the right-side of Equation (61) to represent the same functional as the right-side of Equation (51). This requirement will be satisfied by requiring that the following relations be satisfied by the components of the function ${\psi }^{\left(2\right)}\left(j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]}^{†}\in H_{\left(2\right)}$ for $j1=1,\mathrm{...},JTX$ :

${\left\{A\left(\alpha \right){\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\right\}}_{{\alpha }^0}=-{\left\{S\left(j1;\alpha \right){\psi }^{\left(1\right)}\left(r,\Omega \right)\right\}}_{{\alpha }^0},\text{}\text{}\text{}j1=1,\cdots ,JTX;$ (62)

${\left\{B\left(\alpha \right){\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right\}}_{{\alpha }^0}=-{\left\{S\left(j1;\alpha \right)\phi \left(r,\Omega \right)\right\}}_{{\alpha }^0};j1=1,\cdots ,JTX,$ (63)

where, for each $j1=1,\cdots ,JTX$, $S\left(j1;\alpha \right)$ is a $G\times G$ diagonal matrix having non-zero elements of the form $\partial {\Sigma }_t^g\left(\alpha \right)/\partial t_{j1},g=1,\cdots ,G$ on its diagonal, i.e.,

$S\left(j1;\alpha \right)\triangleq \left(\begin{array}{ccc}\partial {\Sigma }_t^1\left(\alpha \right)/\partial t_{j1}& ·& 0\\ ·& ·& ·\\ 0& ·& \partial {\Sigma }_t^G\left(\alpha \right)/\partial t_{j1}\end{array}\right)\triangleq Diag\left[\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j1}}\right],g=1,\cdots ,G.$ (64)

7) Use in Equation (61) the boundary conditions shown in Equations (22) and (53). Furthermore, set to zero the remaining terms in the bilinear concomitant $P^{\left(2\right)}\left[\delta \phi \left(r,\Omega \right),\delta {\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]$ in Equation (61) by requiring that the components of the function ${\psi }^{\left(2\right)}\left(j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]}^{†}\in H_{\left(2\right)}$ satisfy the following boundary conditions:

${\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)=0,r=r_d,\Omega \cdot n>0;j1=1,\cdots ,JTX;$ (65)

${\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)=0,r=r_d,\Omega \cdot n<0;j1=1,\cdots ,JTX.$ (66)

8) The above boundary conditions complete the well-posed definition of the 2^nd-level adjoint function ${\psi }^{\left(2\right)}\left(j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]}^{†}\in H_{\left(2\right)}$ as the solution of Equations (62) through (66), which are called “the Second-Level Adjoint Sensitivity System (2^nd-LASS)”. The reason for calling this system “Second-Level” (as opposed to “Second-Order”) stems from the fact that this system does not involve any “second-order” differential or derivatives of the dependent variables (i.e., state functions) even though the solution of the 2^nd-LASS is used for computing, efficiently and exactly, the second-order sensitivities of the response with respect to the model parameters.

9) Use Equations (60) through (66) in Equation (51) to obtain the following alternative expression for the differential $\delta {\left\{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)/\partial t_{j1}\right\}}_{{\alpha }^0}$ in terms of the

2^nd-level adjoint function ${\psi }^{\left(2\right)}\left(j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]}^{†}\in H_{\left(2\right)}$ :

$\delta {\left\{\frac{\partial L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)}\right)}{\partial t_{j1}}\right\}}_{{\alpha }^0}$

$\begin{array}{l}={\left\{{\langle {\psi }^{\left(2\right)}\left(j1;r,\Omega \right),\left(\begin{array}{c}{Q}^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\\ Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\end{array}\right)\rangle }_{\left(2\right)}\right\}}_{{\alpha }^0}\\ ={\left\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{\psi }_1^{\left(2\right),g}\left(j1;r,\Omega \right)Q^{\left(1\right)}\left(g;\alpha ,\phi ;\delta \alpha \right)\right\}}_{{\alpha }^0}\\ +{\left\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{\psi }_2^{\left(2\right),g}\left(j1;r,\Omega \right)Q^{\left(2\right)}\left(g;\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\right\}}_{{\alpha }^0}\\ ={\displaystyle \underset{j2=1}{\overset{TP}{\sum }}{\left\{\frac{{\partial }^{2}L\left(\alpha \right)}{\partial {\alpha }_{j2}\partial t_{j1}}\delta {\alpha }_{j2}\right\}}_{{\alpha }^0}},j1=1,\cdots ,JTX.\end{array}$ (67)

10) It is important to note that the differential expression provided in Equation (67) comprises all of the 2^nd-order partial sensitivities that involve the group-averaged microscopic total cross sections. The 2^nd-order partial sensitivities ${\left\{{\partial }^{2}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)/\partial t_{j1}\partial {\alpha }_{j2}\right\}}_{{\alpha }^0}$, $j1=1,\cdots ,JTX$, $j2=1,\cdots ,TP$, which involve the 180 total cross sections ( $t_{j1}$ ) and the other parameters ( ${\alpha }_{j2}$ ), are obtained by replacing in Equation (67) the expressions of $Q^{\left(1\right)}\left(g;\alpha ,\phi ;\delta \alpha \right)$ and $Q^{\left(2\right)}\left(g;\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)$ by their definitions provided in Equations (23) and (54), respectively, and subsequently identifying the expressions that multiply the corresponding parameter variations.

11) In particular, identifying in Equation (67) the quantities that multiply only the variations in the total microscopic cross sections provides the following expression for the 2^nd-order response sensitivities that involve only the group-averaged total microscopic cross sections:

$\begin{array}{l}{\left\{\frac{{\partial }^{2}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)}\right)}{\partial t_{j2}\partial t_{j1}}\right\}}_{{\alpha }^0}\\ =-\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}[{\psi }_1^{\left(2\right),g}\left(j1;r,\Omega \right){\phi }^g\left(r,\Omega \right)\\ +{{\psi }_2^{\left(2\right),g}\left(j1;r,\Omega \right){\psi }^{\left(1\right),g}\left(r,\Omega \right)]\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j2}}\}}_{{\alpha }^0},j1=1,\cdots ,JTX;j2=1,\cdots ,j1.\end{array}$ (68)

4.4. Comparison of Computational Requirements

The numerical results for the 2^nd-order sensitivities of the PERP benchmark’s leakage response to all 7477 nuclear data parameters that characterize the computational model of the PERP benchmark have been presented in [5] - [10]. The 2^nd-order sensitivities were all computed using the PARTISN code with a S₂₅₆ angular quadrature (ISN = 256), using a DELL desktop computer (AMD FX-8350) with an 8-core processor. The CPU-time for a typical adjoint computation was ca. 160 seconds; 0.046 seconds CPU-time was needed for computing a single 2^nd-order sensitivity, including reading the adjoint functions and the integration over the adjoint functions, for ISN = 256. Table 2 presents the comparison of the CPU-times required by the finite-difference, FSAM, and 2^nd-CASAM, respectively,

Table 2. CPU-time needed for computing ${\partial }^{2}L\left(\alpha \right)/\partial t_{j1}\partial t_{j2},j1=1,\cdots ,JTX=180;j2=1,\cdots ,j1$.

for computing the 2^nd-order sensitivities of the leakage response only with respect to the 180 total group-averaged microscopic cross sections.

The number of 16,470 computations required by the FSAM comprises 180 large-scale computations for obtaining the first-order derivatives of the flux with respect to the total microscopic cross sections plus 180(180 + 1)/2 large-scale computations for computing the mixed second-order derivatives of the flux with respect to the total microscopic cross sections. The total CPU time needed for computing all of the 2^nd-order sensitivities using the 2^nd-CASAM was ca. 929 hours, comprising 735 hours used for the 14843 PARTISN computations (S₂₅₆ angular quadrature; ISN = 256) and 194 hours used for performing the integrations needed to compute the respective unmixed and mixed second-order sensitivities. It is evident from the comparison of the CPU-times presented in Table 2 that the 2^nd-CASAM is the only practical methodology for computing 2^nd-order sensitivities for large-scale systems involving many parameters (as illustrated by using the PERP benchmark as a paradigm).

5. Computation of Third-Order Sensitivities ${\partial }^{3}L\left(\alpha \right)/\partial t_{j1}\partial t_{j2}\partial t_{j3},j1=1,\cdots ,JTX;j2=1,\cdots ,j1;j3=1,\cdots ,j2$

Subsection 5.1 presents the finite-difference formulas that can be used to compute approximately the third-order response sensitivities with respect to the model parameters. The Forward Sensitivity Analysis Methodology will be discussed in Subsection 5.2. The 3^rd-CASAM will be discussed in Subsection 5.3, while Subsection 5.4 will present a comparative discussion of the computational resources and times (CPU) which are required by these three deterministic methods.

5.1. Re-Computations with Finite-Difference Approximation

The 3^rd-order unmixed sensitivities, ${\partial }^{3}L\left(\alpha \right)/\partial {\alpha }_j^3$, of the leakage response, $L\left(\alpha \right)$, can be approximately computed by re-computations using the well-known finite-difference formula presented below:

$\frac{{\partial }^{3}L\left(\alpha \right)}{\partial {\alpha }_j^3}\approx \frac{1}{2h_j^3}\left(L_{j+2}-2L_{j+1}+2L_{j-1}-L_{j-2}\right)+{\rm O}\left(h_j^2\right),j=1,\cdots ,TP$, (69)

where $L_{j+2}\triangleq L\left({\alpha }_j+2h_j\right)$, $L_{j+1}\triangleq L\left({\alpha }_j+h_j\right)$, $L_{j-1}\triangleq L\left({\alpha }_j-h_j\right)$, $L_{j-2}\triangleq L\left({\alpha }_j-2h_j\right)$ and $h_j$ denotes a “judiciously-chosen” variation in the parameter ${\alpha }_j$ around its nominal value ${\alpha }_j^0$. The 3^rd-order mixed sensitivities, ${\partial }^{3}L\left(\alpha \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}\partial {\alpha }_{j3}$, can be calculated by using the following finite-difference formula:

$\begin{array}{c}\frac{{\partial }^{3}L\left(\alpha \right)}{\partial {\alpha }_{j1}\partial {\alpha }_{j2}\partial {\alpha }_{j3}}\approx \frac{1}{8h_{j1}{h}_{j2}{h}_{j3}}(L_{j1+1,j2+1,j3+1}-L_{j1+1,j2+1,j3-1}-L_{j1+1,j2-1,j3+1}\\ +L_{j1+1,j2-1,j3-1}-L_{j1-1,j2+1,j3+1}+L_{j1-1,j2+1,j3-1}\\ +L_{j1-1,j2-1,j3+1}-L_{j1-1,j2-1,j3-1})+{\rm O}\left(h_{j1}^2,h_{j2}^2,h_{j3}^2\right),\end{array}$ (70)

where $L_{j1+1,j2+1,j3+1}=L\left({\alpha }_{j1}+h_{j1},{\alpha }_{j2}+h_{j2},{\alpha }_{j3}+h_{j3}\right)$, $L_{j1-1,j2+1,j3+1}=L\left({\alpha }_{j1}-h_{j1},{\alpha }_{j2}+h_{j2},{\alpha }_{j3}+h_{j3}\right)$, etc. The values of the quantities on the rights-sides of the expressions shown in Equations (69) and (70) are obtained by re-solving Equations (1) and (2) repeatedly, using the changed parameter values $\left({\alpha }_{j1}\pm h_{j1}\right)$, $\left({\alpha }_{j1}\pm h_{j2}\right)$ and $\left({\alpha }_{j3}\pm h_{j3}\right)$. As has also been previously discussed, the values of $h_{j1}$, $h_{j2}$ and $h_{j3}$, respectively, must be chosen “judiciously” by “trial and error” for each of the parameters ${\alpha }_{j1}$, ${\alpha }_{j2}$ and ${\alpha }_{j3}$. The finite difference formulas introduce their intrinsic “methodological errors” of order ${\rm O}\left(h_{j1}^2,h_{j2}^2,h_{j3}^2\right)$ which are in addition to, and independent of, the errors that might be incurred in the computation of ${\partial }^{3}L\left(\alpha \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}\partial {\alpha }_{j3}$.

5.2. Forward Sensitivity Analysis Methodology

The third-order partial sensitivity of $L\left(\alpha \right)$ to three generic model parameters, ${\alpha }_{j3},{\alpha }_{j2},{\alpha }_{j1}$, has the following expression:

$\begin{array}{l}{\left\{\frac{{\partial }^{3}L\left(\alpha \right)}{\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0}\triangleq {\left\{{\displaystyle \underset{S_b}{\int }\text{d}S}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{\Omega \cdot n>0}{\int }\text{d}\Omega }}\Omega \cdot n\frac{{\partial }^3{\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0},\\ j1=1,\cdots ,TP;j2=1,\cdots ,j1;j3=1,\cdots ,j2.\end{array}$ (71)

In turn, the third-order derivative ${\partial }^3{\phi }^g\left(r,\Omega \right)/\partial {\alpha }_{j1}\partial {\alpha }_{j2}\partial {\alpha }_{j3}$ is the solution of the third-order G-derivative of Equations (1) and (2) with respect to ${\alpha }_{j3},{\alpha }_{j2},{\alpha }_{j1}$ or, equivalently, the solution of the first G-derivative of the 2^nd-OFSS comprising Equations (47) and (48), which have the following form:

$\begin{array}{l}{\left\{B^g\left(\alpha \right)\left[\frac{{\partial }^3{\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right]\right\}}_{{\alpha }^0}={\left\{\frac{{\partial }^3{S}^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0},\\ g=1,\cdots ,G;j1=1,\cdots ,TP;j2=1,\cdots ,j1;j3=1,\cdots ,j2;\end{array}$ (72)

$\begin{array}{l}{\left\{\frac{{\partial }^3{\phi }^g\left(r,\Omega \right)}{\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}}_{{\alpha }^0}=0,r=r_d;\Omega \cdot n<0,\\ g=1,\cdots ,G;j1=1,\cdots ,TP;j2=1,\cdots ,j1;j3=1,\cdots ,j2;\end{array}$ (73)

where

$\frac{{\partial }^3{S}^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}}=\frac{\partial }{\partial {\alpha }_{j3}}\left\{\frac{{\partial }^2{S}^g\left[{\phi }^g\left(r,\Omega \right),r,\Omega \right]}{\partial {\alpha }_{j2}\partial {\alpha }_{j1}}\right\}.$ (74)

Evidently, the computations of all of the distinctive third-order sensitivities ${\partial }^{3}L\left(\alpha \right)/\partial {\alpha }_{j3}\partial {\alpha }_{j2}\partial {\alpha }_{j1}$ require $TP\left(TP+1\right)\left(TP+2\right)/3!$ large-scale computations to solve Equations (72) and (73), followed by $TP\left(TP+1\right)\left(TP+2\right)/3!$ small-scale computations for performing the integration represented by Equation (71).

5.3. Third-Order Adjoint Sensitivity Analysis Methodology (3rd-CASAM)

In contradistinction to the FSAM, the 3^rd-CASAM does not use Equation (46) as the starting point for computing 3^rd-order response sensitivities but instead uses the expression provided by the 2^nd-CASAM of the sensitivities in terms of the 1^st-and 2^nd-level adjoint functions ${\psi }^{\left(1\right)}$ and ${\psi }^{\left(2\right)}$. Thus, the total G-differential of the expression given in Equation (68) provides the 3^rd-order sensitivities of the leakage response involving the group-averaged total microscopic cross sections, as follows:

$\begin{array}{l}\delta {\left\{\frac{{\partial }^{2}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)}\right)}{\partial t_{j2}\partial t_{j1}}\right\}}_{{\alpha }^0}\\ =-\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}[\delta {\psi }_1^{\left(2\right),g}\left(j1;r,\Omega \right){\phi }^g\left(r,\Omega \right)\\ +{\psi }_1^{\left(2\right),g}\left(j1;r,\Omega \right)\delta {\phi }^g\left(r,\Omega \right)+\delta {\psi }_2^{\left(2\right),g}\left(j1;r,\Omega \right){\psi }^{\left(1\right),g}\left(r,\Omega \right)\\ {+{\psi }_2^{\left(2\right),g}\left(j1;r,\Omega \right)\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)]\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j2}}\}}_{{\alpha }^0},j1=1,\cdots ,JTX;j2=1,\cdots ,j1,\end{array}$ (75)

where: 1) the functions $\delta {\phi }^g\left(r,\Omega \right)$ are the solutions of Equations (25) and (26); 2) the functions $\delta {\psi }^{\left(1\right),g}\left(r,\Omega \right)$ are the solutions of Equations (55) and (56); 3) the functions $\delta {\psi }_{1,j}^{\left(2\right),g}$ and $\delta {\psi }_{2,j}^{\left(2\right),g}$ are the solutions of the G-differentiated 2^nd-LASS defined by Equations (62) through (66), namely:

$\begin{array}{l}{\left\{A\left(\alpha \right)\delta {\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)+S\left(j1;\alpha \right)\delta {\psi }^{\left(1\right)}\left(r,\Omega \right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right);\delta \alpha \right]\right\}}_{{\alpha }^0},\end{array}$ (76)

$\begin{array}{l}{\left\{B\left(\alpha \right)\delta {\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)+S\left(j1;\alpha \right)\delta \phi \left(r,\Omega \right)\right\}}_{{\alpha }^0}\\ ={\left\{Q_2^{\left(3\right)}\left[\alpha ;\phi \left(r,\Omega \right);{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right);\delta \alpha \right]\right\}}_{{\alpha }^0},\end{array}$ (77)

$\delta {\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)=0,r=r_d,\Omega \cdot n>0;j1=1,\cdots ,JTX;$ (78)

$\delta {\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)=0,r=r_d,\Omega \cdot n<0;j1=1,\cdots ,JTX,$ (79)

where

$\begin{array}{l}{Q}_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right);\delta \alpha \right]\\ \triangleq -{\left\{\left[\frac{\partial S\left(j1;\alpha \right)}{\partial \alpha }\delta \alpha \right]{\psi }^{\left(1\right)}\left(r,\Omega \right)\right\}}_{{\alpha }^0}-{\left\{\left[\frac{\partial A\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\right\}}_{{\alpha }^0}\\ \triangleq {\displaystyle \underset{j_3=1}{\overset{TP}{\sum }}\frac{\partial Q_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\right]}{\partial {\alpha }_{j_3}}\delta {\alpha }_{j_3}},\end{array}$ (80)

$\begin{array}{l}{Q}_2^{\left(3\right)}\left[\alpha ;\phi \left(r,\Omega \right);{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right);\delta \alpha \right]\\ \triangleq -{\left\{\left[\frac{\partial S\left(j1;\alpha \right)}{\partial \alpha }\delta \alpha \right]\phi \left(r,\Omega \right)\right\}}_{{\alpha }^0}-{\left\{\left[\frac{\partial B\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right\}}_{{\alpha }^0}\\ ={\displaystyle \underset{j_3=1}{\overset{TP}{\sum }}\frac{\partial Q_2^{\left(3\right)}}{\partial {\alpha }_{j_3}}\delta {\alpha }_{j_3}}.\end{array}$ (81)

Thus, in order to compute the value of the differential shown in Equation (75), it would be first necessary to solve the following system of matrix-equations, subject to the boundary conditions shown in Equations (26), (56), (78) and (79):

$\begin{array}{l}{\left\{\left(\begin{array}{cccc}B\left(\alpha \right)& 0& 0& 0\\ 0& A\left(\alpha \right)& 0& 0\\ 0& S\left(j1;\alpha \right)& A\left(\alpha \right)& 0\\ S\left(j1;\alpha \right)& 0& 0& B\left(\alpha \right)\end{array}\right)\right\}}_{{\alpha }^0}\left(\begin{array}{c}\delta \phi \\ \delta {\psi }^{\left(1\right)}\\ \delta {\psi }_1^{\left(2\right)}\left(j1\right)\\ \delta {\psi }_2^{\left(2\right)}\left(j1\right)\end{array}\right)\\ ={\left\{\left(\begin{array}{c}{Q}^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\\ Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\\ Q_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)};{\psi }_1^{\left(2\right)}\left(j1\right);\delta \alpha \right]\\ Q_2^{\left(3\right)}\left[\alpha ;\phi ;{\psi }_2^{\left(2\right)}\left(j1\right);\delta \alpha \right]\end{array}\right)\right\}}_{{\alpha }^0}.\end{array}$ (82)

It is evident that solving Equation (82) amounts to impractically many large-scale computations, which can be avoided by recasting the right-side of Equation (75) in terms of the solutions of a 3^rd-Level Adjoint Sensitivity System (3^rd-LASS). Following the general principles introduced by Cacuci [12], the construction of the 3^rd-LASS is similar to the construction of the 2^nd-LASS and comprises the following sequence of operations, all of which are to be carried out at the nominal parameter values ${\alpha }^0$ :

1) Consider a Hilbert space, denoted as $H_{\left(3\right)}$, comprising elements denoted generically as $u^{\left(3\right)}\left(r,\Omega \right)\in H_{\left(3\right)}$, which have a “four-component block-vector” structure of the form $u^{\left(3\right)}\left(r,\Omega \right)\triangleq {\left[u_1^{\left(3\right)}\left(r,\Omega \right),u_2^{\left(3\right)}\left(r,\Omega \right),u_3^{\left(3\right)}\left(r,\Omega \right),u_4^{\left(3\right)}\left(r,\Omega \right)\right]}^{†}$, with each of the four vector-components having the following structure: $u_i^{\left(3\right)}\left(r,\Omega \right)\triangleq {\left[u_i^{\left(3\right),1}\left(r,\Omega \right),\cdots ,u_i^{\left(3\right),g}\left(r,\Omega \right),\cdots ,u_i^{\left(3\right),G}\left(r,\Omega \right)\right]}^{†},i=1,2,3,4$.

2) Define the inner product in $H_{\left(3\right)}$, denoted as ${\langle u^{\left(3\right)}\left(r,\Omega \right),v^{\left(3\right)}\left(r,\Omega \right)\rangle }_{\left(3\right)}$, of two arbitrary elements in $H_{\left(3\right)}$, denoted as $u^{\left(3\right)}\left(r,\Omega \right)\in H_{\left(3\right)}$ and $v^{\left(3\right)}\left(r,\Omega \right)\in H_{\left(3\right)}$, with $v^{\left(3\right)}\left(r,\Omega \right)\triangleq {\left[v_1^{\left(3\right)}\left(r,\Omega \right),v_2^{\left(3\right)}\left(r,\Omega \right),v_3^{\left(3\right)}\left(r,\Omega \right),v_4^{\left(3\right)}\left(r,\Omega \right)\right]}^{†}$, and $v_i^{\left(3\right)}\left(r,\Omega \right)\triangleq {\left[v_i^{\left(3\right),1}\left(r,\Omega \right),\cdots ,v_i^{\left(3\right),g}\left(r,\Omega \right),\cdots ,v_i^{\left(3\right),G}\left(r,\Omega \right)\right]}^{†},i=1,2,3,4$, as follows:

$\begin{array}{l}{\langle u^{\left(3\right)}\left(r,\Omega \right),v^{\left(3\right)}\left(r,\Omega \right)\rangle }_{\left(3\right)}\triangleq {\displaystyle \underset{i=1}{\overset{4}{\sum }}{\langle u_i^{\left(3\right)}\left(r,\Omega \right),v_i^{\left(3\right)}\left(r,\Omega \right)\rangle }_{\left(1\right)}}\\ ={\displaystyle \underset{i=1}{\overset{4}{\sum }}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}{u}_i^{\left(3\right),g}\left(r,\Omega \right)}{v}_i^{\left(3\right),g}\left(r,\Omega \right).\end{array}$ (83)

3) Apply the definition provided in Equation (83) to form the inner product of Equation (82) with a yet undefined function ${\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right),\cdots ,{\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)\right]}^{†}\in H_{\left(3\right)}$, to obtain the following relation, evaluated at ${\alpha }^0$ :

$\begin{array}{l}{\left\{{\langle {\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right),\left(\begin{array}{cccc}B\left(\alpha \right)& 0& 0& 0\\ 0& A\left(\alpha \right)& 0& 0\\ 0& S\left(j1;\alpha \right)& A\left(\alpha \right)& 0\\ S\left(j1;\alpha \right)& 0& 0& B\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}\delta \phi \\ \delta {\psi }^{\left(1\right)}\\ \delta {\psi }_1^{\left(2\right)}\left(j1\right)\\ \delta {\psi }_2^{\left(2\right)}\left(j1\right)\end{array}\right)\rangle }_{\left(3\right)}\right\}}_{{\alpha }^0}\\ ={\left\{{\langle {\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right),\left(\begin{array}{c}{Q}^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\\ Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\\ Q_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)};{\psi }_1^{\left(2\right)}\left(j1\right);\delta \alpha \right]\\ Q_2^{\left(3\right)}\left[\alpha ;\phi ;{\psi }_2^{\left(2\right)}\left(j1\right);\delta \alpha \right]\end{array}\right)\rangle }_{\left(3\right)}\right\}}_{{\alpha }^0}.\end{array}$ (84)

4) Use the customary definition of the adjoint operator in $H_{\left(3\right)}$, endowed with the inner product defined in Equation (83) to recast the left side of Equation (84) into the following equivalent form:

$\begin{array}{l}{\left\{{\langle {\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right),\left(\begin{array}{cccc}B\left(\alpha \right)& 0& 0& 0\\ 0& A\left(\alpha \right)& 0& 0\\ 0& S\left(j1;\alpha \right)& A\left(\alpha \right)& 0\\ S\left(j1;\alpha \right)& 0& 0& B\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}\delta \phi \\ \delta {\psi }^{\left(1\right)}\\ \delta {\psi }_1^{\left(2\right)}\left(j1\right)\\ \delta {\psi }_2^{\left(2\right)}\left(j1\right)\end{array}\right)\rangle }_{\left(3\right)}\right\}}_{{\alpha }^0}\\ ={\left\{{\langle {\left(\begin{array}{c}\delta \phi \\ \delta {\psi }^{\left(1\right)}\\ \delta {\psi }_1^{\left(2\right)}\left(j1\right)\\ \delta {\psi }_2^{\left(2\right)}\left(j1\right)\end{array}\right)}^{†},\left(\begin{array}{cccc}A\left(\alpha \right)& 0& 0& S\left(j1;\alpha \right)\\ 0& B\left(\alpha \right)& S\left(j1;\alpha \right)& 0\\ 0& 0& B\left(\alpha \right)& 0\\ 0& 0& 0& A\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}{\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right)\\ {\psi }_2^{\left(3\right)}\left(j2;j1;r,\Omega \right)\\ {\psi }_3^{\left(3\right)}\left(j2;j1;r,\Omega \right)\\ {\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)\end{array}\right)\rangle }_{\left(3\right)}\right\}}_{{\alpha }^0}\\ +P^{\left(3\right)}\left[\delta \phi ,\delta {\psi }^{\left(1\right)};\delta {\psi }^{\left(2\right)}\left(j1\right);{\psi }^{\left(3\right)}\left(j2;j1\right)\right],\end{array}$ (85)

where $P^{\left(3\right)}\left[\delta \phi ,\delta {\psi }^{\left(1\right)};\delta {\psi }^{\left(2\right)}\left(j1\right);{\psi }^{\left(3\right)}\left(j2;j1\right)\right]$ denotes the corresponding bilinear concomitant evaluated on the PERP sphere’s outer boundary $S_b$, at $r=r_d$.

5) Require the first term on the right-side of Equation (85) to represent the same functional as the right-side of Equation (75). This requirement will be satisfied by requiring that the following relations be satisfied by the components of the function ${\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right),\cdots ,{\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)\right]}^{†}\in H_{\left(3\right)}$ :

${\left\{\left(\begin{array}{cccc}A\left(\alpha \right)& 0& 0& S\left(j1;\alpha \right)\\ 0& B\left(\alpha \right)& S\left(j1;\alpha \right)& 0\\ 0& 0& B\left(\alpha \right)& 0\\ 0& 0& 0& A\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}{\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right)\\ {\psi }_2^{\left(3\right)}\left(j2;j1;r,\Omega \right)\\ {\psi }_3^{\left(3\right)}\left(j2;j1;r,\Omega \right)\\ {\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)\end{array}\right)\right\}}_{{\alpha }^0}$

$=-{\left\{\left(\begin{array}{c}S\left(j2;\alpha \right){\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\\ S\left(j2;\alpha \right){\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\\ S\left(j2;\alpha \right)\phi \left(r,\Omega \right)\\ S\left(j2;\alpha \right){\psi }^{\left(1\right)}\left(r,\Omega \right)\end{array}\right)\right\}}_{{\alpha }^0},\text{for}j1=1,\cdots ,JTX;j2=1,\cdots ,j1.$ (86)

6) Use in Equation (86) the boundary conditions shown in Equations (26), (56), (78) and (79). Furthermore, set to zero the remaining terms in the bilinear concomitant $P^{\left(3\right)}\left[\delta \phi ,\delta {\psi }^{\left(1\right)};\delta {\psi }^{\left(2\right)}\left(j1\right);{\psi }^{\left(3\right)}\left(j2;j1\right)\right]$ in Equation (85) by requiring the components of the function ${\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right)\triangleq {\left[{\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right),\cdots ,{\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)\right]}^{†}\in H_{\left(3\right)}$ to satisfy the following boundary conditions:

$\begin{array}{l}{\psi }_i^{\left(3\right)}\left(j2;j1;r,\Omega \right)=0,r=r_d;\Omega \cdot n<0;i=2,3;j1=1,\cdots ,JTX;\text{}\text{}j2=1,\cdots ,j1;\\ {\psi }_i^{\left(3\right)}\left(j2;j1;r,\Omega \right)=0,r=r_d;\Omega \cdot n>0;i=1,4;j1=1,\cdots ,JTX;\text{}\text{}j2=1,\cdots ,j1.\end{array}$ (87)

7) The boundary conditions shown in Equation (87) complete the well-posed definition of the 3^rd-level adjoint function ${\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right)\in H_{\left(3\right)}$ as the solution of Equations (86) and (87), which are called “the Third-Level Adjoint Sensitivity System (3^rd-LASS)”. The reason for calling this system “Third-Level” (as opposed to “Third-Order”) stems from the fact that this system does not involve any 2^nd- and/or 3^rd-order differentials or derivatives of the dependent variables (i.e., state functions) even though the solution of the 3^rd-LASS is used for computing, efficiently and exactly, the 3^rd-order sensitivities of the response with respect to the model parameters.

8) Use Equations (84) through (87) in Equation (75) to obtain the following alternative expression for the differential $\delta {\left\{{\partial }^{2}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)}\right)/\partial t_{j1}\partial t_{j2}\right\}}_{{\alpha }^0}$ in terms of the 3rd-level adjoint function ${\psi }^{\left(3\right)}\left(j2;j1;r,\Omega \right)\in H_{\left(3\right)}$ :

$\begin{array}{l}\delta {\left\{\frac{{\partial }^{2}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)}\right)}{\partial t_{j2}\partial t_{j1}}\right\}}_{{\alpha }^0}\\ ={\left\{{\langle {\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right),Q^{\left(1\right)}\left(\alpha ;\delta \alpha \right)\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle {\psi }_2^{\left(3\right)}\left(j2;j1;r,\Omega \right),Q^{\left(2\right)}\left(\alpha ,{\psi }^{\left(1\right)};\delta \alpha \right)\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle {\psi }_3^{\left(3\right)}\left(j2;j1;r,\Omega \right),Q_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)};{\psi }_1^{\left(2\right)}\left(j1\right);\delta \alpha \right]\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ +{\left\{{\langle {\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right),Q_2^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)};{\psi }_2^{\left(2\right)}\left(j1\right);\delta \alpha \right]\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ ={\displaystyle \underset{j3=1}{\overset{TP}{\sum }}{\left\{\frac{{\partial }^{3}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)};{\psi }^{\left(3\right)}\right)}{\partial {\alpha }_{j3}\partial t_{j2}\partial t_{j1}}\delta {\alpha }_{j3}\right\}}_{{\alpha }^0}},j1=1,\cdots ,JTX;j2=1,\cdots ,j1.\end{array}$ (88)

9) It is important to note that the differential expression provided in Equation (88) comprises all of the 3^rd-order partial sensitivities that involve the group-averaged microscopic total cross sections, i.e., the 3^rd-order partial sensitivities ${\partial }^{3}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)};{\psi }^{\left(3\right)}\right)/\partial {\alpha }_{j3}\partial t_{j2}\partial t_{j1}$, $j1=1,\cdots ,JTX$; $j2=1,\cdots ,j1$, $j3=1,\cdots ,TP$, describes a total of $\left[180\left(180+1\right)/2\right]\times 7477$ 3^rd-order response sensitivities with respect to the microscopic total cross sections (indices $t_{j1}$ and $t_{j2}$ ) and with respect to the other parameters ( ${\alpha }_{j3}$ ).

10) In particular, identifying in Equation (88) the quantities that multiply only the variations in the total microscopic cross sections will reduce the expressions provided in Equations (80) and (81) as follows:

$\begin{array}{l}{Q}_1^{\left(3\right)}\left[\alpha ;{\psi }^{\left(1\right)}\left(r,\Omega \right);{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right);\delta \alpha \right]={\displaystyle \underset{j_3=1}{\overset{JTX}{\sum }}\frac{\partial S_1^{\left(3\right)}\left[\alpha ;{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\right]}{\partial t_{j_3}}\delta t_{j_3}},\\ S_1^{\left(3\right)}\left[\alpha ;{\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\right]\triangleq -{\Sigma }_t^g\left(\alpha \right){\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right),\end{array}$ (89)

$\begin{array}{l}{Q}_2^{\left(3\right)}\left[\alpha ;\phi \left(r,\Omega \right);{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right);\delta \alpha \right]\triangleq {\displaystyle \underset{j_3=1}{\overset{JTX}{\sum }}\frac{\partial S_2^{\left(3\right)}\left[\alpha ;{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]}{\partial t_{j_3}}\delta t_{j_3}},\\ S_2^{\left(3\right)}\left[\alpha ;{\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\right]\triangleq -{\Sigma }_t^g\left(\alpha \right){\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right).\end{array}$ (90)

11) Inserting the results from Equations (89) and (90) into Equation (88) and identifying the specific quantity that multiplies the specific variation $\delta t_{j_3}$ provides the following expression for the 3^rd-order response sensitivities which involve only the group-averaged total microscopic cross sections:

$\begin{array}{l}{\left\{\frac{{\partial }^{3}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)};{\psi }^{\left(3\right)}\right)}{\partial t_{j3}\partial t_{j2}\partial t_{j1}}\right\}}_{{\alpha }^0}\\ =-{\left\{{\langle {\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right),S\left(j3;\alpha \right)\phi \left(r,\Omega \right)\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ -{\left\{{\langle {\psi }_2^{\left(3\right)}\left(j2;j1;r,\Omega \right),S\left(j3;\alpha \right){\psi }^{\left(1\right)}\left(r,\Omega \right)\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ -{\left\{{\langle {\psi }_3^{\left(3\right)}\left(j2;j1;r,\Omega \right),S\left(j3;\alpha \right){\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0}\\ -{\left\{{\langle {\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right),S\left(j3;\alpha \right){\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)\rangle }_{\left(1\right)}\right\}}_{{\alpha }^0},\\ \text{for}j1=1,\cdots ,JTX;j2=1,\cdots ,j1;j3=1,\cdots ,j2.\end{array}$ (91)

In component form, Equation (91) can be written as follows:

$\begin{array}{l}{\left\{\frac{{\partial }^{3}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)};{\psi }^{\left(3\right)}\right)}{\partial t_{j3}\partial t_{j2}\partial t_{j1}}\right\}}_{{\alpha }^0}\\ =-\{{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{0}{\overset{r_d}{\int }}4\pi r^2\text{d}r}{\displaystyle \underset{4\pi }{\int }\text{d}\Omega }}[{\psi }_1^{\left(3\right),g}\left(j2;j1;r,\Omega \right){\phi }^g\left(r,\Omega \right)\\ +{\psi }_2^{\left(3\right),g}\left(j2;j1;r,\Omega \right){\psi }^{\left(1\right),g}\left(r,\Omega \right)+{\psi }_3^{\left(3\right),g}\left(j2;j1;r,\Omega \right){\psi }_1^{\left(2\right),g}\left(j1;r,\Omega \right)\\ {+{\psi }_4^{\left(3\right),g}\left(j2;j1;r,\Omega \right){\psi }_2^{\left(2\right),g}\left(j1;r,\Omega \right)]\frac{\partial {\Sigma }_t^g\left(\alpha \right)}{\partial t_{j3}}\}}_{{\alpha }^0},\\ \text{for}j1=1,\cdots ,JTX;j2=1,\cdots ,j1;j3=1,\cdots ,j2.\end{array}$ (92)

5.4. Comparison of Computational Requirements

As has already been discussed, the number of sensitivities of the PERP leakage response which involve solely the group-averaged total microscopic cross sections are as follows: 1) 180 first-order sensitivities; 2) 32,400 second-order sensitivities, of which 16,290 are distinct; and 3) 5,832,000 third-order sensitivities, of which 988,260 are distinct; 4) 1,049,760,000 fourth-order sensitivities, of which 45,212,895 are distinct. Numerical results for these sensitivities have been obtained by using the code PARTISN with a S₃₂ angular quadrature (ISN = 32). The computer used for these computations is a DELL computer with an 8-core processor (AMD FX-8350).

The CPU-time for a typical adjoint computation using ISN = 32 is ca. 26 seconds, while the CPU-time computing the integrals over the various adjoint functions which appear in the definition of the respective sensitivity, cf. Equation (92), is ca. 0.002 seconds.

The CPU-time for computing the 1^st-, 2^nd- and 3^rd-order sensitivities of the leakage response with respect to the 180 microscopic total cross sections are as follows:

1) The total CPU-time required for computing the 1^st-order sensitivities is ca. 30 Seconds (1 adjoint PARTISN computation, plus 180 integrations over the adjoint functions);

2) The total CPU-time required for computing twice (for verification purposes) the 2^nd-order sensitivities is ca. 3 hours (2 h 55 min for 180 forward-like + 180 adjoint-like PARTISN computations for determining the 2^nd-level adjoint functions, plus 5 minutes for computing the 32,400 integrations over the adjoint functions);

3) The total CPU-time required for computing twice (for verification purposes) the 3^rd-order sensitivities is around 338 hours (2 h 55 min for computing 2^nd-level adjoint functions, 332 hours for computing the 3^rd-level adjoint functions + 3 h 24 min for computing the 180 × 180 × 180 integrals over these functions to obtain the 3^rd-order sensitivities).

Table 3 presents the comparison of the CPU-times for all of the distinct third-order sensitivities ${\partial }^{3}L\left(\alpha \right)/\partial t_{j3}\partial t_{j2}\partial t_{j1},j1=1,\cdots ,JTX;j2=1,\cdots ,j1;j3=1,\cdots ,j2$ required by applying the 3^rd-CASAM, versus using the FSAM procedure using Equation (71) and/or the FD-approximations shown in Equations (69) and (70).

Table 3. Comparison of CPU-times for computing ${\partial }^{3}L\left(\alpha ,\phi ;{\psi }^{\left(1\right)};{\psi }^{\left(2\right)};{\psi }^{\left(3\right)}\right)/\partial t_{j3}\partial t_{j2}\partial t_{j1},j1=1,\cdots ,JTX;j2=1,\cdots ,j1;j3=1,\cdots ,j2$.

Based on Equations (69) and (70), the number of forward computations needed to obtain all of the distinct third-order sensitivities by the FD-approximation method is 7,905,360 [= 180 × 4 forward PARTISN computations for the 180 unmixed 3^rd-order sensitivities + (988,260 − 180) × 8 forward PARTISN computations for the 988,080 mixed 3^rd-order sensitivities].

Based on Equations (72) and (73), the number of forward computations needed to obtain all of the distinct third-order sensitivities ${\partial }^{3}L\left(\alpha \right)/\partial t_{j3}\partial t_{j2}\partial t_{j1}$ by the FSAM method is $JTX+JTX\left(JTX+1\right)/2+JTX\left(JTX+1\right)\left(JTX+2\right)/3!=1004730$.

Based on Equation (86), a total number of 33,301 adjoint computations is needed to obtain all of the distinct third-order sensitivities ${\partial }^{3}L\left(\alpha \right)/\partial t_{j3}\partial t_{j2}\partial t_{j1}$ by the 3^rd-CASAM method, comprising the following large-scale computations:

1) 1 adjoint PARTISN computations to obtain ${\psi }^{\left(1\right)}\left(r,\Omega \right)$;

2) 180 adjoint-like PARTISN computations to obtain ${\psi }_1^{\left(2\right)}\left(j1;r,\Omega \right)$;

3) 180 forward-like PARTISN computations to obtain ${\psi }_2^{\left(2\right)}\left(j1;r,\Omega \right)$;

4) 16,290 adjoint-like PARTISN computations to obtain ${\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right)$;

5) 16,290 forward-like PARTISN computations to obtain ${\psi }_2^{\left(3\right)}\left(j2;j1;r,\Omega \right)$;

6) 180 forward-like PARTISN computations to obtain ${\psi }_3^{\left(3\right)}\left(j2;j1;r,\Omega \right)$;

7) 180 adjoint-like PARTISN computations to obtain ${\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)$.

The computation of each of the adjoint functions ${\psi }_1^{\left(3\right)}\left(j2;j1;r,\Omega \right)$ and ${\psi }_2^{\left(3\right)}\left(j2;j1;r,\Omega \right)$, respectively, requires $180\left(180+1\right)/2=16290$ PARTISN computations, respectively, since the source terms for the equations satisfied by these functions depend on both indices j1 and j2. However, the computation of each of the adjoint functions ${\psi }_3^{\left(3\right)}\left(j2;j1;r,\Omega \right)$ and ${\psi }_4^{\left(3\right)}\left(j2;j1;r,\Omega \right)$, respectively, requires only 180 PARTISN computations (per function), because the source terms for the equations satisfied by these functions only depend on the index j2, as shown in Equation .

The results presented in Table 3 evidently indicate that the 3^rd-CASAM is the only methodology capable of computing 3^rd-order sensitivities, without introducing methodological errors.

6. Concluding Remarks

This work has presented the derivation of the exact mathematical expressions of the following sensitivities of the PERP leakage response with respect to the group-averaged total microscopic cross sections: 1) 180 first-order sensitivities; 2) 32,400 second-order sensitivities, of which 16,290 are distinct; and 3) 5,832,000 third-order sensitivities, of which 988,260 are distinct.

Using a DELL computer with an 8-core processor (AMD FX-8350), the total CPU-time required for computing the 1^st-order sensitivities of the PERP leakage response with respect to all of the PERP benchmark’s parameters by applying the 1^st-CASAM is 95 seconds. By comparison, the forward sensitivity analysis methodology (FSAM) would require 694 hours CPU-time, while a 2-point finite-difference (FD) approximate scheme would require 1,388 hours CPU-time.

Using the same DELL computer (AMD FX-8350), the computations of the 2^nd-order sensitivities of the PERP leakage response with respect to the PERP benchmark’s total microscopic cross sections by applying the 2^nd-CASAM required 19 hours CPU-time. By comparison, the FSAM would require 764 hours CPU-time, while the FD-scheme would require 3006 hours CPU-time.

Finally, using the AMD FX-8350 DELL computer, the computations of the 3^rd-order sensitivities of the PERP leakage response with respect to the PERP benchmark’s total microscopic cross sections by applying the 3^rd-CASAM required 175 hours CPU-time. By comparison, the FSAM would require 12,559 hours CPU-time, while the FD-scheme would require 98,817 hours CPU-time. These results make it abundantly evident that the CASAM is the only practical method for computing sensitivities exactly for large-scale problems involving many parameters.

The formulas derived in this work are valid not only for the PERP benchmark but can also be used for computing the 3^rd-order sensitivities of the leakage response of any nuclear system involving fissionable material and internal or external neutron sources. Subsequent works [16] [17] will use the adjoint-based mathematical expressions obtained in this work to further derive the expressions for the 4^th-order sensitivities, and to compute exactly and efficiently the numerical values of the 1,049,760,000 fourth-order sensitivities, of which 45,212,895 are distinct.

Appendix

The dimensions and material composition of the polyethylene-reflected plutonium (PERP) metal sphere considered in this work are presented in TableA1.

The quantities appearing in Equations (1) through (3) are defined as follows:

1) The quantity ${\phi }^g\left(r,\Omega \right)$ is the customary “group-flux” for group $g,g=1,\cdots ,G$, and is the unknown state-function which is obtained by solving Equations (1) and (2). The group boundaries of the G=30 energy groups are provided in TableA2.

2) The source $Q^g\left(q\right)$ depends on the vector of model parameters $q$, defined as follows:

Table A1. Dimensions and material composition of the PERP benchmark.

Table A2. Group boundaries, $E^g\left[\text{MeV}\right]$, of the $G=30$ energy groups used in the PARTISN forward and adjoint neutron transport computations.

$q\triangleq {\left[q_1,\cdots ,q_{JQ}\right]}^{†}\triangleq {\left[{\lambda }_1,{\lambda }_2;F_1^{SF},F_2^{SF};a_1,a_2;b_1,b_2;{\nu }_1^{SF},{\nu }_2^{SF}\right]}^{†},JQ=10.$ (93)

3) As indicated in TableA1, the PERP benchmark comprises 2 materials: “material 1” comprises 4 isotopes, numbered 1 through 4, while “material 2” comprises 2 isotopes, numbered 5 and 6. In principle, PARTISN allows the same isotope to appear in different materials, in which case the atomic number density $N_{i,m}$ of an isotope i in a material m would be computed by using the formula $N_{i,m}={\rho }_m{w}_{i,m}{N}_A/A_i$, where ${\rho }_m$ denotes the mass density of material m, $m=1,2$, $w_{i,m}$ denotes the weight fraction of isotope i in material m; $A_i$ denotes the atomic weight of isotope i, and $N_A$ denotes Avogadro’s number. However, the two materials in the PERP benchmark contain only isotopes that are distinct from each other, so the subscript m will not be needed if the formula for the atomic number density $N_i$ of an isotope i, $i=1,\cdots ,I=6$, is interpreted as follows:

$N_i\triangleq \frac{{\rho }_1{w}_{i,1}{N}_A}{A_i};\text{for}i=1,2,3,4;\text{}\text{}{N}_i\triangleq \frac{{\rho }_2{w}_{i,2}{N}_A}{A_i};\text{for}i=5,6.\text{}$ (94)

The atomic number densities $N_i,i=1,\cdots ,I=6$ will be considered to be components of the vector $N$, defined below:

$N\triangleq {\left[N_1,N_2,N_3,N_4,N_5,N_6\right]}^{†}.$ (95)

4) The scattering transfer cross section from energy group $g^{\prime },g^{\prime }=1,\cdots ,G$ into energy group $g,g=1,\cdots ,G$ is denoted as ${\Sigma }_s^{g^{\prime }\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)$ and is computed in terms of the l-th order Legendre coefficient ${\sigma }_{s,l,i}^{g^{\prime }\to g}$ using the following 3^rd-order expansion in Legendre functions:

${\Sigma }_s^{g^{\prime }\to g}\left(\alpha ;r,{\Omega }^{\prime }\to \Omega \right)={\displaystyle \underset{i=1}{\overset{I=6}{\sum }}{N}_i}{\displaystyle \underset{l=0}{\overset{ISCT=3}{\sum }}\left(2l+1\right){\sigma }_{s,l,i}^{g^{\prime }\to g}(r)P_l\left({\Omega }^{\prime }\cdot \Omega \right),}\text{}\text{}\text{}g,g^{\prime }=1,\cdots ,G$, (96)

where ISCT = 3 denotes the order of the expansion in Legendre polynomials. The microscopic scattering cross sections ${\sigma }_{s,l,i}^{g^{\prime }\to g}$ for isotope i, and from energy group $g^{\prime }$ into energy group g, are tabulated parameters. The zeroth-order (i.e., $l=0$ ) scattering cross sections must be considered separately from the higher order (i.e., $l\ge 1$ ) scattering cross sections, since the former contribute to the total cross sections (as noted below), while the latter do not. Aiming at reducing the proliferation of superscripts and subscripts when defining response sensitivities with respect to the microscopic scattering cross sections ${\sigma }_{s,l,i}^{g^{\prime }\to g}$, these cross sections will be considered to be components of a vector ${\sigma }_s$ defined below:

$\begin{array}{l}{\sigma }_s\triangleq {\left[s_1,\cdots ,s_{JSX}\right]}^{†}\\ \triangleq {\left[{\sigma }_{s,l=0,i=1}^{g^{\prime }=1\to g=1},{\sigma }_{s,l=0,i=1}^{g^{\prime }=2\to g=1},\cdots ,{\sigma }_{s,l=0,i=1}^{g^{\prime }=G\to g=1},{\sigma }_{s,l=0,i=1}^{g^{\prime }=1\to g=2},{\sigma }_{s,l=0,i=1}^{g^{\prime }=2\to g=2},\cdots ,{\sigma }_{s,l,i}^{g^{\prime }\to g},\cdots ,{\sigma }_{s,ISCT,i=I}^{G\to G}\right]}^{†},\\ l=0,\cdots ,ISCT;i=1,\cdots ,I;g,g^{\prime }=1,\cdots ,G;JSX=\left(G\times G\right)\times I\times \left(ISCT+1\right).\end{array}$ (97)

5) The total cross section for energy group $g,g=1,\cdots ,G$ is denoted as ${\Sigma }_t^g\left(\alpha ;r\right)$ and is computed using the following expression:

${\Sigma }_t^g\left(\alpha ;r\right)={\displaystyle \underset{i=1}{\overset{I}{\sum }}{N}_i{\sigma }_{t,i}^g};{\sigma }_{t,i}^g=[{\sigma }_{f,i}^g\left(r\right)+{\sigma }_{c,i}^g\left(r\right)+{\displaystyle \underset{g^{\prime }=1}{\overset{G}{\sum }}{\sigma }_{s,l=0,i}^{}}$