On the Need to Determine Accurately the Impact of Higher-Order Sensitivities on Model Sensitivity Analysis, Uncertainty Quantification and Best-Estimate Predictions

Abstract

This work aims at underscoring the need for the accurate quantification of the sensitivities (i.e., functional derivatives) of the results (a.k.a. “responses”) produced by large-scale computational models with respect to the models’ parameters, which are seldom known perfectly in practice. The large impact that can arise from sensitivities of order higher than first has been highlighted by the results of a third-order sensitivity and uncertainty analysis of an OECD/NEA reactor physics benchmark, which will be briefly reviewed in this work to underscore that neglecting the higher-order sensitivities causes substantial errors in predicting the expectation and variance of model responses. The importance of accurately computing the higher-order sensitivities is further highlighted in this work by presenting a text-book analytical example from the field of neutron transport, which impresses the need for the accurate quantification of higher-order response sensitivities by demonstrating that their neglect would lead to substantial errors in predicting the moments (expectation, variance, skewness, kurtosis) of the model response’s distribution in the phase space of model parameters. The incorporation of response sensitivities in methodologies for uncertainty quantification, data adjustment and predictive modeling currently available for nuclear engineering systems is also reviewed. The fundamental conclusion highlighted by this work is that confidence intervals and tolerance limits on results predicted by models that only employ first-order sensitivities are likely to provide a false sense of confidence, unless such models also demonstrate quantitatively that the second- and higher-order sensitivities provide negligibly small contributions to the respective tolerance limits and confidence intervals. The high-order response sensitivities to parameters underlying large-scale models can be computed most accurately and most efficiently by employing the high-order comprehensive adjoint sensitivity analysis methodology, which overcomes the curse of dimensionality that hampers other methods when applied to large-scale models involving many parameters.

1. Introduction

The state-of-the-art computational tool for sensitivity analysis, uncertainty quantification and “data adjustment” of model responses (effective multiplication factors, reaction rates) in reactor physics and reactor criticality safety is embodied in the module TSURFER [1] of the code system SCALE, developed at the Oak Ridge National Laboratory. The methodology underlying TSURFER is the so-called “generalized linear least squares adjustment (GLLSA)” procedure [2], which is limited to incorporating just the first-order sensitivities of the response with respect to the model parameters. The first-order response sensitivities to the many thousands of model parameters (cross-sections, number densities, etc.) involved in reactor physics computations are computed efficiently and accurately using the first-order adjoint sensitivity analysis method.

Since a general-purpose adjoint sensitivity analysis methodology capable of computing second-order sensitivities efficiently, reliably and accurately for large-scale systems (such as encountered in reactor physics) was unavailable in the past, the potential effects of higher-order sensitivities were not considered. However, a general-purpose, efficient and accurate methodology was recently developed by Cacuci [3,4,5] and was applied, in a sequence of pioneering works [6,7,8,9,10,11], to perform a comprehensive second-order sensitivity and uncertainty analysis of an OECD/NEA polyethylene-reflected plutonium (acronym: PERP) reactor physics benchmark [12], which comprised 21,976 uncertain models parameters (including 180 group-averaged total microscopic cross-sections, 21,600 group-averaged scattering microscopic cross-sections, 120 parameters describing the fission process, 60 parameters describing the fission spectrum, 10 parameters describing the system’s sources and six isotopic number densities). Thus, the PERP benchmark comprises 21,976 first-order sensitivities of the leakage response with respect to the model parameters and 482,944,576 second-order sensitivities (of which 241,483,276 are distinct from each other). These fundamental works [6,7,8,9,10,11] have demonstrated that, for this benchmark, many second-order sensitivities were significantly larger (by over an order of magnitude) than the first-order sensitivities. Consequently, the cumulative effects of the second-order sensitivities on the predicted uncertainties of the PERP benchmark’s response far exceeded the effects of the first-order sensitivities. For example, for fully correlated total cross-sections, neglecting the second-order sensitivities would 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.

The results obtained in [6,7,8,9,10,11] have motivated the quest for investigating the third-order sensitivities of a large-scale linear system such as the PERP benchmark. For this purpose, Cacuci [12] developed “The Third-Order Adjoint Method for Sensitivity Analysis of Response-Coupled Linear Forward/Adjoint Systems”, which was applied to compute [13,14,15] the largest third-order sensitivities of the PERP benchmark response, namely the 5,832,000 third-order sensitivities of the PERP response to the 180 total cross-sections. The results obtained in [13,14,15] indicated that, for a uniform standard deviation of 10% of the microscopic total cross-sections, the third-order sensitivities contribute 80% to the response variance, whereas the contribution stemming from the first- and second-order sensitivities amount only to 2% and 18%, respectively. Neglecting the third-order sensitivities could cause a very large non-conservative error by under-reporting the response variance by a factor of 506%.

The results obtained in [14,15,16] have motivated the investigation reported in this work, which demonstrates that significant, non-negligible values of second- and higher-order sensitivities, such as those observed in [6,7,8,9,10,11,14,15,16], are not an isolated phenomenon but are actually generic to the fundamental physical processes of neutron scattering (“slowing down”) in any medium. This work is structured as follows: Section 2 briefly reviews the computational model of the PERP benchmark and the impact of the second- and third-order PERP sensitivities on the expectation and variance of the PERP benchmark’s leakage response. The computational model of the PERP benchmark involves the solution of the neutron slowing-down and transport (Boltzmann) equation, which is generic to modeling any reactor physics or shielding system and, therefore, serves to the introduction of the simplified physical model of neutron slowing down, which was originally developed by Bethe [17] and used ever since [18,19,20] for computing the intra-group fluxes needed to generate multigroup cross-section libraries, such as those contained in SCALE [1], which are in turn needed to solve the multigroup form of the Boltzmann equation. Notably, a sensitivity analysis of neutron slowing down processes does not seem to have been published in the literature thus far. Section 3 presents the exact high-order sensitivity analysis of the neutron slowing down equation/model, showing that all of the relative sensitivities of the neutron flux are unity, so none can be dismissed a priori as “unimportant”. The exact uncertainty analysis if the slowing down model is also presented in Section 3, underscoring the importance of the higher-order sensitivities in contributing to the expectation, variance, skewness and kurtosis of the neutron flux.

Section 4 presents a comparison of the end results for the “adjusted” values produced by the TSURFER-GLLSA procedure [1,2] for the model parameters (in this case: nuclear cross-sections) and responses (in this case: effective multiplication factors) with the end-results produced by the HO-BERRU-PM (high-order best-estimate results with reduced uncertainties predictive modeling) methodology conceived by Cacuci [13,21]. The end-results for the “adjusted” values produced by the TSURFER-GLLSA procedure for the nuclear cross-sections (considered as “model parameters”) and effective multiplication factors (considered as “model responses”) are obtained by minimizing a user-chosen quadratic cost functional, which represents the “discrepancies between measurements and computations”. In contradistinction to the TSURFER-GLLSA procedure [1,2], the HO-BERRU-PM methodology uses the maximum entropy principle [22] to eliminate the need for any “user-defined functional to be minimized”, while incorporating all of the available high-order information, thus yielding higher-order best-estimate predictions for the model parameters and responses. The distinctions between the TSURFER-GSSLA procedure (which can consider only first-order sensitivities) and the HO-BERRU-PM methodology (which can incorporate arbitrarily high order sensitivities) are summarized in Section 4, for both the a priori and posterior (i.e., adjusted vs. predicted) quantities, highlighting the contributions of the higher-order sensitivities.

Section 5 concludes this work by underscoring the need for computing sensitivities of an order higher than one and underscoring the impact of these higher-order sensitivities on predicted values for the moments (expectation, variance, skewness, kurtosis) of the distribution of model responses in the phase space of parameters. The results presented in this work underscore the fundamental conclusion that—at the very least—the second-order sensitivities must be computed quantitatively, in order to show quantitatively if their impact is negligible (or not) on confidence intervals and tolerance limits for the model responses under investigation. In the absence of any information regarding the impact of the second- and higher-order sensitivities, confidence intervals and tolerance limits based on first-order sensitivities only for model responses cannot be trusted.

2. Motivation: Impact of Second- and Third-Order Response Sensitivities on Uncertainty Analysis of the PERP OECD/NEA Reactor Physics Benchmark

The impact of the second- and third-order response sensitivities on the expected value and variance of the model response under consideration is illustrated in this section by recalling results obtained in [6,7,8,9,10,11] for the polyethylene-reflected plutonium (acronym: PERP) OECD/NEA reactor physics benchmark [12], which is a sphere used for benchmark subcritical neutron and gamma measurements. The constitutive materials of the PERP benchmark are specified in

Table 1. Dimensions and material composition of the PoRP benchmark.

Materials	Isotopes	Weight Fraction	Density (g/cm3)	Zones
Material 1 (plutonium metal)	Isotope 1 (239Pu)	9.3804 × 10−1	19.6	Material 1 is assigned to zone 1; inner radius = 3.794 cm.
	Isotope 2 (240Pu)	5.9411 × 10−2
	Isotope 3 (69Ga)	1.5152 × 10−3
	Isotope 4 (71Ga)	1.0346 × 10−3
Material 2 (polyethylene)	Isotope 5 (C)	8.5630 × 10−1	0.95	Material 2 is assigned to zone 2; inner radius = 3.794 cm; outer radius = 7.604 cm.
	Isotope 6 (1H)	1.4370 × 10−1

.

The distribution of neutrons in the PERP system is modeled using the standard form of the time-independent Boltzmann neutron transport equation, namely:

$$\begin{array}{l}\Omega ·\nabla \phi \left(r,\Omega ,E\right)+{\mathsf{\Sigma }}_t\left(r,E\right) \phi \left(r,\Omega ,E\right)={\displaystyle \underset{4\pi }{\int }d{\Omega }^{\prime }{\displaystyle \underset{0}{\overset{E_f}{\int }}d{E}^{\prime } {\mathsf{\Sigma }}_s\left(r,E^{\prime }\to E,{\Omega }^{\prime }\to \Omega \right) \phi \left(r,{\Omega }^{\prime },E^{\prime }\right)}}\\ +{\displaystyle \underset{4\pi }{\int }d{\Omega }^{\prime }{\displaystyle \underset{0}{\overset{E_f}{\int }}d{E}^{\prime } \chi \left(r,E^{\prime }\to E\right)\nu {\mathsf{\Sigma }}_f\left(r,E^{\prime }\right)\phi \left(r,{\Omega }^{\prime },E^{\prime }\right)}} +Q\left(r,E\right),\end{array}$$

(1)

where $Q\left(r,E\right)$ is a spontaneous fission source, having the following expression:

$$Q\left(r,E\right)\triangleq {\displaystyle \sum _{k=1}^{N_f}{\lambda }_k{N}_{k,1}{F}_k^{SF}{\nu }_k^{SF}\left(\frac{2}{\sqrt{\pi a_k{}^3{b}_k}}{e}^{-\frac{a_k{b}_k}{4}}\right)e^{-E/a_k}\sinh \sqrt{b_{k}E}}.$$

(2)

The neutron angular flux $\phi \left(r,\Omega ,E\right)$ depends on the energy variable $E$, the solid angle variable $\Omega$ and the three-dimensional position vector $r$. The angular flux $\phi \left(r,\Omega ,E\right)$ is subject to the customary vacuum boundary condition, which specifies that there is no incoming flux of particles, namely:

$$\phi \left(r_s,\Omega ,E\right)=0,r_s\in \partial V,\Omega \cdot n<0,0<E<\infty ,$$

(3)

where $n$ denotes the unit outward normal vector at any point $r_s\in \partial V$ on the body’s outer surface $\partial V$, and $E_f$ denotes the highest neutron energy considered for the system under consideration. The terms on the left side of Equation (1) represent the neutron streaming and total collision loss processes, while the terms on the right side of Equation (1) represent the neutron scattering and fission production processes within the PERP benchmark.

The neutron flux is computed by solving numerically Equations (1)–(3) using the multigroup discrete ordinates particle and radiation transport code PARTISN [23], which solves the following multi-group approximation of Equation (1):

$$\begin{array}{l}\Omega ·\nabla {\phi }_g\left(r,\Omega \right)+{\mathsf{\Sigma }}_{t,g}\left(\alpha ;r\right) {\phi }_g\left(r,\Omega \right)=\frac{{\chi }_g\left(\alpha ;r\right)}{k_{eff}}\mathsf{\Phi }\left(\alpha ;r\right)+Q_g\left(\alpha ;r\right)\\ +{\displaystyle \sum _{g^{\prime }=1}^G{\displaystyle \underset{4\pi }{\int }d{\Omega }^{\prime }{\phi }_{g^{\prime }}\left(r,{\Omega }^{\prime }\right){\mathsf{\Sigma }}_{s,g^{\prime }\to g}\left(\alpha ;r,{\mu }_0\right) }}, g=1,\dots ,G,\end{array}$$

(4)

$${\phi }_g\left(r_d,\Omega \right)=0,\Omega \cdot n<0, g=1,\dots ,G$$

(5)

where ${\mu }_0\triangleq \Omega ·{\Omega }^{\prime }$, where $r_d$ is the radius of the PERP sphere and where the boundary condition shown in Equation (5) imposes no incoming particle flux. The vector $\alpha$, which appears in the arguments of the various quantities in Equation (4), represents the “vector of imprecisely known model parameters”, including atomic number densities, microscopic group cross-sections, etc. The detailed definition of the components of the vector $\alpha$ for the PERP benchmark is provided in [6,7,8,9,10,11], but this detailed definition is not used in this work. The convention used in PARTISN [23] for obtaining the multigroup transport equation shown in Equation (4) is to partition the range of the energy variable into $G$ intervals of width $\Delta E_g\triangleq E_{g-1/2}-E_{g+1/2}, g=1,\dots ,G$ and, subsequently, integrate Equations (1)–(3) over each interval $\Delta E_g$. By convention, the highest energy is at $E_{1/2}$. Since the normal transport of particles in energy is from high to low energy, as they collide with nuclei in the medium, the index $g$ increases as energy decreases. Integrating Equations (1)–(3) over each energy interval $\Delta E_g$ leads to Equations (4) and (5), where the following definitions for the various group-averaged quantities were used:

$${\phi }_g\left(r,\Omega \right)\triangleq {\displaystyle \underset{\Delta E_g}{\int }dE \phi \left(r,\Omega ,E\right)}$$

(6)

$${\mathsf{\Sigma }}_{t,g}\left(\alpha ;r\right)\triangleq \frac{1}{{\phi }_g\left(r,\Omega \right)}{\displaystyle \underset{\Delta E_g}{\int }dE }{\mathsf{\Sigma }}_t\left(\alpha ;r,E\right) \phi \left(r,\Omega ,E\right)$$

(7)

$${\mathsf{\Sigma }}_{s,g^{\prime }\to g}\left(\alpha ;r,{\mu }_0\right)\triangleq \frac{1}{{\phi }_{g^{\prime }}\left(r,\Omega \right)}{\displaystyle \underset{\Delta E_g}{\int }dE{\displaystyle \underset{\Delta E_{g^{\prime }}}{\int }d{E}^{\prime } }}{\mathsf{\Sigma }}_s\left(\alpha ;r,E^{\prime }\to E,{\mu }_0\right) \phi \left(r,\Omega ,E^{\prime }\right)$$

(8)

$${\chi }_g\left(\alpha ;r\right)\triangleq {\displaystyle \underset{\Delta E_g}{\int }dE }\chi \left(\alpha ;r;E\right)$$

(9)

$$\mathsf{\Phi }\left(\alpha ;r\right)\triangleq \frac{1}{{\displaystyle \underset{4\pi }{\int }{\phi }_g\left(r,\Omega \right)d\Omega }}{\displaystyle \sum _{g=1}^G{\displaystyle \underset{\Delta E_g}{\int }dE }\nu {\mathsf{\Sigma }}_f\left(\alpha ;r,E\right){\displaystyle \underset{4\pi }{\int }\phi \left(r,\Omega ,E\right)d\Omega } }$$

(10)

$$Q_g\left(\alpha ;r\right)\triangleq {\displaystyle \sum _{k=1}^{N_f}{\lambda }_k{N}_{k,1}{F}_k^{SF}{\nu }_k^{SF}\left(\frac{2}{\sqrt{\pi a_k{}^3{b}_k}}{e}^{-\frac{a_k{b}_k}{4}}\right){\displaystyle \underset{\Delta E_g}{\int }dE e^{-E/a_k}\sinh \sqrt{b_{k}E}}}.$$

(11)

The group fluxes ${\phi }^g\left(r,\Omega \right)$ are computed numerically using the MENDF71X [24] 618-group cross-sections collapsed to 30 energy groups, in conjunction with a P₃-Legendre expansion of the scattering cross-section, an angular quadrature of S₂₅₆ 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 other quantities that appear in Equations (4)–(11) are described in detail in [6,7,8,9,10,11].

The response of interest for the PERP benchmark is the leakage of neutrons out of the sphere’s outer surface. The leakage response, denoted as $L\left(\alpha \right)$, is defined mathematically as the following integral over the outer surface of the spherical benchmark, comprising neutrons of all energies leaking through the surface:

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

(12)

As shown in [6,7,8,9,10,11], the distribution of the PERP leakage in each energy interval is depicted in Figure 1. The value of the total leakage computed using Equation (1) for the PERP benchmark is 1.7648 × 10⁶ neutrons/s.

The vector $\alpha$, which appears in the arguments of the various quantities in Equations (4)–(12), represents the “vector of imprecisely known model parameters”. For the PERP benchmark, the vector $\alpha$ comprises 21,976 imprecisely known (uncertain) model parameters, as follows: 180 group-averaged total microscopic cross-sections, 21,600 group-averaged scattering microscopic cross-sections, 120 parameters describing the fission process, 60 parameters describing the fission spectrum, 10 parameters describing the system’s sources and six isotopic number densities. The exact and efficient computation of the 21,976 first-order and 482,944,576 second-order sensitivities of the leakage response to the PERP model parameters was reported in [6,7,8,9,10,11].

The results obtained in [6,7,8,9,10,11,14,15,16] indicated that the largest first-, second- and third-order relative sensitivities of the PERP leakage response were all with respect to the total microscopic group cross-section, ${\sigma }_{t,H}^{30}$, of hydrogen (¹H) in the lowest energy group, $g=30$. Thus, the third-order relative sensitivity, $S^{(3)}\left({\sigma }_{t,H}^{30},{\sigma }_{t,H}^{30},{\sigma }_{t,H}^{30}\right)\triangleq \left[{\partial }^{3}L/{\left(\partial {\sigma }_{t,H}^{30}\right)}^3\right]\left[{\left({\sigma }_{t,H}^{30}\right)}^3/L\right]=-2.966\times {10}^4$, was found to be significantly larger than the second-order sensitivity,

$S^{(2)}\left({\sigma }_{t,H}^{30},{\sigma }_{t,H}^{30}\right)\triangleq \left[{\partial }^{2}L/{\left(\partial {\sigma }_{t,H}^{30}\right)}^2\right]\left[{\left({\sigma }_{t,H}^{30}\right)}^2/L\right]=4.296\times {10}^2$, which in turn was an order of magnitude larger than the first-order relative sensitivity $S^{(1)}\left({\sigma }_{t,H}^{30}\right)\triangleq \left[\partial L/\left(\partial {\sigma }_{t,H}^{30}\right)\right]\left[\left({\sigma }_{t,H}^{30}\right)/L\right]=-9.366\times {10}^0$.

The results presented in

Table 2. Comparison of response moments for various relative standard deviations of uncorrelated and normally distributed microscopic total cross-sections (all numbers denote neutrons/second).

Relative Standard Deviation	5%	10%
$L\left({\alpha }^0\right)$	1.765 × 106	1.765 × 106
${\left[E\left(L\right)\right]}_t^{(2,U,N)}$	1.149 × 106	4.598 × 106
${\left[E\left(L\right)\right]}_t^{(U,N)}=L\left({\alpha }^0\right)+{\left[E\left(L\right)\right]}_t^{(2,U,N)}$	2.914 × 106	6.363 × 106
${\left[\mathrm{var} (L)\right]}_t^{\left(1,U,N\right)}$	8.549 × 1011	3.419 × 1012
${\left[\mathrm{var} (L)\right]}_t^{\left(2,U,N\right)}$	1.799 × 1012	2.879 × 1013
${\left[\mathrm{var} (L)\right]}_t^{\left(3,U,N\right)}$	8.713 × 1012	1.308 × 1014
$\begin{array}{l}{\left[\mathrm{var} (L)\right]}_t^{(U,N)}={\left[\mathrm{var} (L)\right]}_t^{\left(1,U,N\right)}\\ +{\left[\mathrm{var} (L)\right]}_t^{\left(2,U,N\right)}+{\left[\mathrm{var} (L)\right]}_t^{\left(3,U,N\right)}\end{array}$	1.083 × 1013	1.630 × 1014

, reproduced (with permission) from [6], illustrate the impact of uncertainties in the total group cross-sections on the response’s expected value and variance. The group total cross-sections were considered to be uncorrelated and normally distributed, having uniform relative standard deviations of 5% (moderate) and 10% (large), respectively. The impact of the other uncertain cross-sections and nuclear parameter underlying the PERP benchmark are presented in [6,7,8,9,10,11], where they are shown to accentuate the impacts presented in Table 2 and Figure 2, below.

In Table 2, the expected value of the leakage response is denoted as ${\left[E\left(L\right)\right]}_t^{\left(U,N\right)}$ and has the following mathematical expression:

$${\left[E\left(L\right)\right]}_t^{\left(U,N\right)}=L\left({\alpha }^0\right)+{\left[E\left(L\right)\right]}_t^{\left(2,U,N\right)},$$

(13)

where ${\left[E\left(L\right)\right]}_t^{\left(2,U,N\right)}$ denotes the second-order contributions to the expected value of the leakage response and is defined as follows:

$${\left[E\left(L\right)\right]}_t^{\left(2,U,N\right)}\triangleq \frac{1}{2}{\displaystyle \sum _{g=1}^G{\displaystyle \sum _{i=1}^I\frac{{\partial }^{2}L\left(\alpha \right)}{\partial {\sigma }_{t,i}^g\partial {\sigma }_{t,i}^g}{\left(s_{t,i}^g\right)}^2}} , G=30 energy groups, I=6 isotopes.$$

(14)

In Equation (13), the quantity $L\left({\alpha }^0\right)$ represents the leakage response computed using the expected cross-section values; the superscript “U” denotes “uncorrelated” and the superscript “N” denotes “normally distributed”. The quantity $s_{t,i}^g$, which appears in Equation (14), denotes the standard deviation associated with the total microscopic cross-section for energy group $g, g=1,\dots ,G,$ and isotope $i, i=1,\dots ,I$.

The variance of the leakage response of the PERP benchmark is denoted as ${\left[\mathrm{var} (L)\right]}_t^{(U,N)}$ in Table 2 and has the following mathematical expression:

$${\left[\mathrm{var} (L)\right]}_t^{(U,N)}={\left[\mathrm{var} (L)\right]}_t^{\left(1,U,N\right)}+{\left[\mathrm{var} (L)\right]}_t^{\left(2,U,N\right)}+{\left[\mathrm{var} (L)\right]}_t^{\left(3,U,N\right)}.$$

(15)

In Equation (15), the quantity ${\left[\mathrm{var} (L)\right]}_t^{\left(1,U,N\right)}\triangleq {\displaystyle \sum _{g=1}^G{\displaystyle \sum _{i=1}^I{\left[\frac{\partial L\left(\alpha \right)}{\partial {\sigma }_{t,i}^g}\right]}^2{\left(s_{t,i}^g\right)}^2}}$ denotes the contributions stemming from the first-order sensitivities; the quantity ${\left[\mathrm{var} (L)\right]}_t^{\left(2,U,N\right)}\triangleq \frac{1}{2}{\displaystyle \sum _{g=1}^G{\displaystyle \sum _{i=1}^I{\left[\frac{{\partial }^{2}L\left(\alpha \right)}{\partial {\sigma }_{t,i}^g\partial {\sigma }_{t,i}^g}{\left(s_{t,i}^g\right)}^2\right]}^2}}$ denotes the contributions stemming from the second-order sensitivities; the quantity ${\left[\mathrm{var} (L)\right]}_t^{\left(3,U,N\right)}\triangleq {\displaystyle \sum _{g=1}^G{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{g^{\prime }=1}^G{\displaystyle \sum _{j=1}^I\left[\frac{\partial L\left(\alpha \right)}{\partial {\sigma }_{t,i}^g}\frac{{\partial }^{3}L\left(\alpha \right)}{\partial {\sigma }_{t,i}^g\partial {\sigma }_{t,j}^{g^{\prime }}\partial {\sigma }_{t,j}^{g^{\prime }}}\right]}}{\left(s_{t,i}^g\right)}^2}}{\left(s_{t,j}^{g^{\prime }}\right)}^2$ denotes the contributions stemming from the third-order sensitivities, for $G=30$ energy groups and $I=6$ isotopes. The results presented in the second column of Table 2 are illustrated graphically in Figure 2, which is reproduced (with permission) from [6].

As illustrated in Figure 2, the relationship $S{D}^{\left(3\right)}\gg S{D}^{\left(2\right)}>S{D}^{\left(1\right)}$ indicates that the contribution from the second-order sensitivities exceed the contributions from the first-order ones, and even more importantly, the contributions from the third-order sensitivities are dominant, being much larger than the contributions from either the first-order or the second-order sensitivities. Hence, neglecting the third-order sensitivities would cause a significant error in quantifying the standard deviation of the leakage response. The results in Table 2 and Figure 2 also indicate that the contributions from the second-order sensitivities cause a significant shift, by about 40%, of the expected value, ${\left[E\left(L\right)\right]}_t^{\left(U,N\right)}$, of the leakage response by comparison to the leakage computed value $L\left({\alpha }^0\right)$. This means that the customary procedure of neglecting second- (and higher) order sensitivities and considering that the computed value, $L\left({\alpha }^0\right)$, is the same as the expected value ${\left[E\left(L\right)\right]}_t^{\left(U,N\right)}$ of the leakage response would be ca. 40% in error when considering uniform 5% relative standard deviations for uncorrelated total cross-sections. The results presented in Table 2 also underscore the fact that the larger the parameters’ standard deviations, the larger the impact of the second-order and third-order sensitivities. It was shown in [6,7,8,9,10,11] that correlations among the cross-sections enhance the effects of the high-order sensitivities.

The results obtained in [6,7,8,9,10,11] and summarized in this section, which indicated that the PERP leakage response was most sensitive to the hydrogen cross-sections, have motivated the high-order sensitivity and uncertainty analysis to be presented next, in Section 3.

3. High-Order Sensitivity and Uncertainty Analysis of Neutron Slowing Down in Hydrogen

The multigroup fluxes and cross-sections defined in Equations (6)–(10) are used not only in PARTISN [23] but are used in all of the deterministic solvers (e.g., SCALE [2]) of the Boltzmann multi-group neutron transport equation. Although the multigroup Boltzmann equation shown in Equation (4) is formally exact (i.e., no approximations have been made in deriving it), this equation cannot be solved as it stands since the energy-dependent flux $\phi \left(r,\Omega ,E\right)$, which enters into the definitions of the multigroup fluxes and cross-sections defined in Equations (6)–(10), is not known at that stage. The first step in the usual procedure for determining the energy-dependent flux $\phi \left(r,\Omega ,E\right)$ is to partition the energy interval of interest into very many (“fine” or “ultra-fine”) energy groups and solve the “infinite-medium” energy-dependent Boltzmann quasi-analytically to obtain the energy-dependent “infinite-medium” flux [18,19,20]. In view of the results presented in Section 2, it is of interest to investigate the sensitivities of the energy-dependent neutron flux with respect to the total microscopic cross-section of hydrogen in a hydrogen-moderated system, since the energy-dependent neutron flux is directly related, via Equation (12), to the leakage response.

The main features of the energy-dependent neutron flux are determined by considering a medium of large (“infinite”) extent, in which neutrons emitted by a source (which may be an external source or an internal fission source) will lose energy (“slow down”) following scattering collisions and will disappear following absorption collisions with the nuclei in the medium. A widely used model [17,18,19,20], which includes the scattering and absorption processes of neutrons emitted by an energy-dependent source $Q\left(E\right)$ in a hydrogen-moderated system (which includes absorbing nuclide(s) having atomic mass(es) much larger than unity), has the following form:

$${\displaystyle {\int }_E^{\infty }\mathsf{\Phi }\left(E^{\prime }\right) }{\mathsf{\Sigma }}_s\left(E^{\prime }\right)\frac{d{E}^{\prime }}{E^{\prime }}+Q\left(E\right)=\mathsf{\Phi }\left(E\right)\mathsf{\Sigma }\left(E\right); \mathsf{\Sigma }\left(E\right)\triangleq {\mathsf{\Sigma }}_s\left(E\right)+{\mathsf{\Sigma }}_a\left(E\right) .$$

(16)

The integral Equation (16) for the energy-dependent neutron flux, $\mathsf{\Phi }\left(E\right)$, is derived under the assumptions that: (i) the absorbing isotope does not contribute to scattering, so that ${\mathsf{\Sigma }}_s\left(E\right)$ represents the scattering cross-section of hydrogen, and (ii) the absorption cross-section ${\mathsf{\Sigma }}_a\left(E\right)$ stems mostly from the heavy nuclide in the hydrogen-moderated system, since the absorption in hydrogen is negligible by comparison to its scattering cross-section. Without affecting the sensitivity analysis to follow, the macroscopic cross-sections can be considered to be just the corresponding products of the respective atomic number densities and microscopic cross-sections, i.e., ${\mathsf{\Sigma }}_s\left(E\right)=N_H{\sigma }_s\left(E\right)$ and ${\mathsf{\Sigma }}_a\left(E\right)=N_a{\sigma }_a\left(E\right)$, respectively. Additionally, without affecting the sensitivity analysis, the total cross-section characterizing the hydrogen-moderated system can be considered to be simply the sum of the two component macroscopic cross-sections, i.e., $\mathsf{\Sigma }\left(E\right)=N_H{\sigma }_s\left(E\right)+N_a{\sigma }_a\left(E\right)$. Under the additional assumption that the (resonance) absorption sets in below the low-energy limit, $E_s$, of the fission spectrum (an assumption that also does not affect the sensitivity analysis to follow), the solution of Equation (16) has the following form:

$$\mathsf{\Phi }\left(\alpha ;E\right)=\frac{p\left(\alpha ;E\right){\displaystyle {\int }_{E_s\left(\alpha \right)}^{\infty }Q\left(\alpha ;E^{\prime }\right)d{E}^{\prime } }}{E\mathsf{\Sigma }\left(\alpha ;E\right)}; p\left(\alpha ;E\right)\triangleq \exp \left[-{\displaystyle {\int }_E^{E_s\left(\alpha \right)}\frac{{\mathsf{\Sigma }}_a\left(\alpha ;E^{\prime }\right)}{\mathsf{\Sigma }\left(\alpha ;E^{\prime }\right)} }\frac{d{E}^{\prime }}{E^{\prime }}\right].$$

(17)

Evidently, the flux $\mathsf{\Phi }\left(\alpha ;E\right)$ depends on the fission spectrum and the various number densities and microscopic cross-sections of the isotopes contained in the medium under consideration. Therefore, it is convenient to use the generic notation introduced in Section 2 and, thus, to consider all of the uncertain model parameters that influence the flux $\mathsf{\Phi }\left(\alpha ;E\right)$ to be the components of a “vector of parameters” denoted as $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}\in {ℝ}^{TP}$, where $TP$ denotes the “total number of model parameters”. Variations in the model’s parameters (i.e., cross-sections, number densities, etc.) will induce variations in the flux $\mathsf{\Phi }\left(\alpha ;E\right)$, the quantification of which is the scope of sensitivity analysis. Uncertainties in the model parameters will induce uncertainties in the flux, which can be quantified by using the expressions provided in [13].

As indicated by Equation (17), the neutron flux $\mathsf{\Phi }\left(\alpha ;E\right)$ varies as $1/E$ and undergoes sharp dips/depressions at resonance absorption energies, where $\mathsf{\Sigma }\left(\alpha ;E\right)$ becomes very large. In addition, the flux $\mathsf{\Phi }\left(\alpha ;E\right)$ is attenuated by the resonance escape probability, $p\left(\alpha ;E\right)$, as the neutrons undergo moderating scattering collisions. The qualitative form of $\mathsf{\Phi }\left(\alpha ;E\right)$ shown in Equation (17) is also representative of the neutron slowing down flux in different mixtures of materials/isotopes [18,19,20].

The impact of the scattering process can be qualitatively determined by considering a system in which neutron absorption is negligible, which would be the case of a neutron slowing down in pure hydrogen. In such a system, the energy-dependent flux would take on the following particular form of Equation (17):

$${\mathsf{\Phi }}_H\left(\alpha ;E\right)=\frac{Q\left(\alpha \right)}{E{\mathsf{\Sigma }}_s\left(\alpha ;E\right)}; Q\left(\alpha \right)\triangleq {\displaystyle {\int }_{E_s\left(\alpha \right)}^{\infty }Q\left(\alpha ;E^{\prime }\right)d{E}^{\prime } }\left(E\right).$$

(18)

Although Equation (18) is strictly correct only for an infinite hydrogen medium, it is often used to estimate the energy-dependent flux in water-moderated reactors [18,19,20].

The effects of the scattering cross-sections can be quantified (in isolation from the effects of the other model parameters) by considering that the source does not vary from its nominal value denoted as $Q\left({\alpha }^0\right)$ and that the number density of hydrogen, $N_H$, remains unchanged, as well. Consider that the energy-dependent microscopic total cross-section of hydrogen, ${\sigma }_s\left(E\right)$, was measured in a series of measurements at an arbitrary but fixed energy value $E_m$ $0<E_m<\infty$ and was found to have a mean value denoted as ${\sigma }_m$. Under these assumptions, the flux ${\mathsf{\Phi }}_H\left(\alpha ;E\right)$ in Equation (18) takes on the following expression:

$${\mathsf{\Phi }}_H\left(\sigma ;E\right)=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_s\left(E\right)}.$$

(19)

The Taylor-series expansion of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ around the nominal parameter value ${\sigma }_m$ has the following form:

$${\mathsf{\Phi }}_H\left(\sigma ;E\right)=\frac{Q\left({\alpha }^0\right)}{E{N}_H}{\displaystyle \sum _{n=0}^{\infty }\frac{1}{n!}}{\left\{\frac{d^n}{d{\left[{\sigma }_s\left(E\right)\right]}^n}{\left[{\sigma }_s\left(E\right)\right]}^{-1}\right\}}_{\left({\sigma }_m\right)}{\left[{\sigma }_s\left(E\right)-{\sigma }_m\right]}^n.$$

(20)

Recalling that

$$\frac{d^n}{d{x}^n}\left(x^{-1}\right)={\left(-1\right)}^{n}n!x^{-\left(n+1\right)}, n=0,1,\dots .$$

(21)

and using the result from Equation (21) in Equation (20) yields the following Taylor series for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$:

$${\mathsf{\Phi }}_H\left(\sigma ;E\right)=\underset{n\to \infty }{\lim }\left[{\mathsf{\Phi }}_H^{\left(n\right)}\left(\sigma ;E\right)\right],$$

(22)

where the quantity ${\mathsf{\Phi }}_H^{\left(n\right)}\left(\sigma ;E\right)$ denotes the “nth-order approximation of the value of the actual flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$” and is defined as follows:

$${\mathsf{\Phi }}_H^{\left(n\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\displaystyle \sum _{k=0}^n{\left(-1\right)}^n{t}^n;} t\triangleq \frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}.$$

(23)

The Taylor series in Equation (22) converges only for values of ${\sigma }_s\left(E\right)$ confined within the following interval:

$$\left(1-\beta \right){\sigma }_m\le {\sigma }_s\left(E\right)\le \left(1+\beta \right){\sigma }_m; 0<\beta <1.$$

(24)

The results obtained in Equations (22) and (24) are essential for using the probability distribution of ${\sigma }_s\left(E\right)$, which is seldom known completely in practice, to construct statistical indicators (i.e., mean value, standard deviation, etc.) for the distribution of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$. It is especially important to ensure that the values of ${\sigma }_s\left(E\right)$ used for constructing statistical properties for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ do not fall outside the interval of convergence shown in Equation (24).

In view of the relation shown in Equation (24), it follows that if a standard deviation, $SD\left[{\sigma }_s\left(E\right)\right]$, is also provided (in addition to the measured mean value ${\sigma }_m$, then it must be ensured that the Taylor series shown in Equation (22) is used only for the following range of values:

$$\beta {\sigma }_m\le SD\left[{\sigma }_s\left(E\right)\right]\le \beta {\sigma }_m; 0<\beta <1.$$

(25)

If a standard deviation is not provided for the measurement of ${\sigma }_s\left(E_m\right)$, then the most conservative assumption is to consider that ${\sigma }_s\left(E\right)$ is uniformly distributed within the interval provided in Equation (24), in which the Taylor series of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ convergences, since the uniform distribution is the least informative (hence, the least biased) distribution. It will henceforth be assumed that the distribution for ${\sigma }_s\left(E\right)$ is uniform in the interval given in Equation (24). Thus, the uniform probability distribution, $P\left[{\sigma }_s\left(E\right)\right]$, the expectation, $E\left[{\sigma }_s\left(E\right)\right]$, the variance, $Var\left[{\sigma }_s\left(E\right)\right]$, and the standard deviation, $SD\left[{\sigma }_s\left(E\right)\right]$, of ${\sigma }_s\left(E\right)$ have the following expressions:

$$\begin{array}{l}P\left[{\sigma }_s\left(E\right)\right]=\frac{1}{2\beta {\sigma }_m}; E\left[{\sigma }_s\left(E\right)\right]={\sigma }_m; Var\left[{\sigma }_s\left(E\right)\right]=\frac{{\left(\beta {\sigma }_m\right)}^2}{3}; \\ SD\left[{\sigma }_s\left(E\right)\right]=\frac{1}{\sqrt{3}}\beta {\sigma }_m; 0<\beta <1 .\end{array}$$

(26)

The third and fourth moments, ${\mu }_3\left[{\sigma }_s\left(E\right)\right]$ and ${\mu }_4\left[{\sigma }_s\left(E\right)\right]$, of the uniform distribution ${\sigma }_s\left(E\right)$ have the following expressions:

$$\begin{array}{l}{\mu }_3\left[{\sigma }_s\left(E\right)\right]\triangleq {\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left[{\sigma }_s\left(E\right)-{\sigma }_m\right]}^3} d{\sigma }_s\left(E\right)=0,\\ {\mu }_4\left[{\sigma }_s\left(E\right)\right]\triangleq {\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left[{\sigma }_s\left(E\right)-{\sigma }_m\right]}^4} d{\sigma }_s\left(E\right)=\frac{{\left(\beta {\sigma }_m\right)}^4}{5} .\end{array}$$

(27)

The uniform distribution defined in Equation (26) is symmetrical with respect to ${\sigma }_m$ and can be visualized as a rectangle of height $P\left[{\sigma }_s\left(E\right)\right]={\left[2\beta {\sigma }_m\right]}^{-1}$ and width extending from $\left(1-\beta \right){\sigma }_m$ to $\left(1+\beta \right){\sigma }_m$ for $0<\beta <1$. When $\beta$ increases towards unity, the uniform distribution $P\left[{\sigma }_s\left(E\right)\right]={\left[2\beta {\sigma }_m\right]}^{-1}$ widens and its height decreases towards $0.5/{\sigma }_m$. At the other extreme, when $\beta$ decreases towards zero, the distribution $P\left[{\sigma }_s\left(E\right)\right]={\left[2\beta {\sigma }_m\right]}^{-1}$ becomes very narrow and peaked, tending towards a Dirac-delta distribution centered at ${\sigma }_m$.

It follows from the expressions provided in Equations (26) and (24) that, if ${\sigma }_s\left(E\right)$ is distributed uniformly as shown in Equation (26), then the Taylor series shown in Equation (22) for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ converges only in the following interval expressed in terms of “standard deviations of ${\sigma }_s\left(E\right)$”:

$${\sigma }_m-\sqrt{3}\cdot SD\left[{\sigma }_s\left(E\right)\right]\le {\sigma }_s\left(E\right)\le {\sigma }_m+\sqrt{3}\cdot SD\left[{\sigma }_s\left(E\right)\right]; 0<\beta <1.$$

(28)

3.1. Sensitivity Analysis

It follows from the expression given in Equation (20) that the absolute sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, which will be denoted as $S_n\left(E\right), n=1,2,\dots$, with respect to the parameter ${\sigma }_s\left(E\right)$ at energy $E=E_m$, where ${\sigma }_s\left(E_m\right)={\sigma }_m$, have the following expressions:

$$S_n\left(E\right)=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\left(-1\right)}^n{\left[\frac{1}{{\sigma }_m}\right]}^n.$$

(29)

It follows from the expression in Equation (29) that the relative sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to the parameter ${\sigma }_s\left(E\right)$, which will be denoted as $R_n\left(E\right), n=1,2,\dots$, have the following expressions:

$$R_n\left(E\right)=S_n\left(E\right)\frac{{\left({\sigma }_m\right)}^n}{{\mathsf{\Phi }}_H\left({\sigma }_m;E\right)}={\left(-1\right)}^n, n=1,2,\dots .$$

(30)

It is evident from Equation (30) that the relative values of the sensitivities of all orders are all equal to unity; their values do not decrease with increasing order $n=1,2,\dots$.

The relative error, ${\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}$, between the exact value of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ provided in Equation (19) and the Taylor series in Equation (22), truncated to the nth order is given by the following expression:

$${\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}\triangleq \frac{1}{{\mathsf{\Phi }}_H\left(\sigma ;E\right)}\left[{\mathsf{\Phi }}_H\left(\sigma ;E\right)-\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\displaystyle \sum _{k=0}^n{\left(-1\right)}^k{t}^k}\right]={\left(-1\right)}^{n+1}{t}^{n+1}.$$

(31)

It follows from Equation (31) that:

$$\left|{\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}\right|=\left|t^{n+1}\right|=\left|{\left[\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}\right]}^{n+1}\right|<1.$$

(32)

In turn, it follows from Equation (32) that the smallest order, $n$, which is required to ensure that the Taylor-series expansion shown in Equation (22) represents the exact flux given in Equation (19) within an error $\left|{\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}\right|$, must satisfy the following inequality:

$$n>\frac{\left\{-\log \left|{\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}\right|\right\}}{\left\{-\log \left|t\right|\right\}}-1.$$

(33)

In particular, Equations (33) and (24) indicate that when the parameter ${\sigma }_s\left(E\right)$ attains the largest value for which the Taylor-series representation given in Equation (22) still converges, namely, ${\sigma }_s\left(E\right)=\left(1+\beta \right){\sigma }_m$, the smallest order, $n$, which is required to ensure that the Taylor-series expansion shown in Equation (22) represents the exact flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ given in Equation (19) within an error $\left|{\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}\right|$, must satisfy the following inequality:

$$n>\frac{\left\{-\log \left|{\epsilon }_{\mathsf{\Phi }}^{\left(n\right)}\right|\right\}}{\left\{-\log \beta \right\}}-1.$$

(34)

As an example of the use of the Taylor-series expansion in sensitivity analysis, consider that it is desired to predict the value of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ when ${\sigma }_s\left(E\right)$ is “one standard deviation away from the mean”, i.e., when:

$${\sigma }_s^{1SD}\left(E\right)={\sigma }_m+SD\left[{\sigma }_s\left(E\right)\right]={\sigma }_m\left(1+\frac{1}{\sqrt{3}}\beta \right).$$

(35)

Replacing the result obtained in Equation (35) into the definition provided in Equation (23) and taking $\beta =1-{10}^{-4}$(to be extremely close to the outer boundaries of the “one-standard deviation interval”) yields the following result when keeping the first four decimal digits:

$$t=\beta /\sqrt{3}= \left(1-{10}^{-4}\right)/\sqrt{3}\cong 0.5773.$$

(36)

In sensitivity analysis, the value provided in Equation (36) is used in Equation (23) to investigate how accurately the various orders of approximations, ${\mathsf{\Phi }}_H^{\left(n\right)}\left(\sigma ;E\right)$, can predict the actual value of the flux at ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$. The following results, to four significant decimals, are obtained from Equations (19), (23) and (35) for $\beta =1-{10}^{-4}$ and $t=0.5773$:

• The exact value of the flux is:

  $${\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\left(1+\frac{1-{10}^{-4}}{\sqrt{3}}\right)}^{-1}\cong 0.6340\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}.$$

  (37)
• If only the first-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the first-order expansion ($n=1$) in Equation (23) yields the following result:

  $${\mathsf{\Phi }}_H^{\left(1\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t\right)=0.4227\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}.$$

  (38)

It follows from Equations (37) and (38) that the error of the first-order approximation is:

$${\epsilon }_{\mathsf{\Phi }}^{\left(1\right)}\triangleq \frac{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)-{\mathsf{\Phi }}_H^{\left(1\right)}\left(\sigma ;E\right)}{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)}=33.33\%.$$

(39)

Thus, the result obtained by keeping only the first-order sensitivities in the Taylor expansion underestimates the exact flux by 33.34%.

3.

If first-order and second-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the second-order expansion ($n=2$) in Equation (23) yields the following result:

$${\mathsf{\Phi }}_H^{\left(2\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t+t^2\right)=0.7560\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}.$$

(40)

It follows from Equations (40) and (37) that the error of the second-order approximation is as follows:

$${\epsilon }_{\mathsf{\Phi }}^{\left(2\right)}\triangleq \frac{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)-{\mathsf{\Phi }}_H^{\left(2\right)}\left(\sigma ;E\right)}{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)}=-19.24\%.$$

(41)

Thus, the result obtained by keeping the first- and second-order sensitivities in the Taylor expansion overestimates the exact flux by 19.24%.

4.

If all sensitivities up to and including the third-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the third-order expansion ($n=3$) in Equation (23) yields the following result:

$${\mathsf{\Phi }}_H^{\left(3\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t+t^2-t^3\right)=0.5636\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}.$$

(42)

It follows from Equations (42) and (37) that the error of the third-order approximation is as follows:

$${\epsilon }_{\mathsf{\Phi }}^{\left(3\right)}\triangleq \frac{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)-{\mathsf{\Phi }}_H^{\left(3\right)}\left(\sigma ;E\right)}{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)}=11.11\%.$$

(43)

Thus, the result obtained by keeping the first-, second- and third-order sensitivities in the Taylor expansion underestimates the exact flux by 11.11%.

5.

If all sensitivities up to and including the fourth-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the fourth-order expansion ($n=4$) in Equation (23) yields the following result:

$${\mathsf{\Phi }}_H^{\left(4\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t+t^2-t^3+t^4\right)=0.6747\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}.$$

(44)

It follows from Equations (44) and (37) that the error of the fourth-order approximation is as follows:

$${\epsilon }_{\mathsf{\Phi }}^{\left(4\right)}\triangleq \frac{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)-{\mathsf{\Phi }}_H^{\left(4\right)}\left(\sigma ;E\right)}{{\mathsf{\Phi }}_H\left({\sigma }_s^{1SD};E\right)}=-6.42\%.$$

(45)

Thus, the result obtained by keeping the first-, second-, third- and fourth-order sensitivities in the Taylor expansion overestimates the exact flux by 6.42%.

As indicated by the results in Equations (39), (41), (43) and (45) the odd-order approximations underestimate the true value, and the even-order approximations overestimate the true value of the flux, while converging to the actual, true value of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$. This behavior is expected, since the flux is represented by an alternating series, cf. Equation (23). The results obtained in Equations (39), (41), (43) and (45) are consistent with the error estimate provided in Equation (31). The relation provided in Equation (34) can be used to determine which order of approximation would need to be used to obtain an approximate result that would be within a predetermined error by comparison to the exact result. For example, over 46,000 terms in the Taylor-series expansion provided in Equation (23) would be required to represent the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with an accuracy of 1% at a distance of ${\sigma }_s\left(E\right)=\left(1+ 1-{10}^{-4}\right){\sigma }_m$ away from ${\sigma }_m$; the expected extremely slow convergence of the Taylor-series expansion very close to the end of the interval of convergence is evident.

3.2. Uncertainty Quantification: Moments of the Response Distribution

The expectation, $E\left[{\mathsf{\Phi }}_H\right]$, of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, is obtained by integrating either Equation (19) or Equation (22) over the uniform probability distribution for ${\sigma }_s\left(E\right)$ shown in Equation (26). The result of the respective integration is:

$$\begin{array}{l}E\left[{\mathsf{\Phi }}_H\right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H}\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}\frac{d{\sigma }_s\left(E\right)}{{\sigma }_s\left(E\right)}}=\frac{Q\left({\alpha }^0\right)}{E{N}_H}\frac{1}{2\beta {\sigma }_m}\ln \frac{1+\beta }{1-\beta }\\ =\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\displaystyle \sum _{n=0}^{\infty }\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left[\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}\right]}^{n}d{\sigma }_s\left(E\right)}}=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\displaystyle \sum _{k=0}^{\infty }\frac{{\beta }^{2k}}{2k+1}}\\ =\frac{Q\left({\alpha }^0\right)}{E{N}_H\left(\beta {\sigma }_m\right)}\mathrm{Arctanh}\left(\beta \right); 0<\beta <1.\end{array}$$

(46)

Note that $E\left[{\mathsf{\Phi }}_H\right]$ is an even function of $\beta$, i.e.,

$$E\left[{\mathsf{\Phi }}_H,\beta \right]=E\left[{\mathsf{\Phi }}_H,-\beta \right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H}\frac{1}{2\beta {\sigma }_m}\ln \frac{1+\beta }{1-\beta }.$$

(47)

The nth-order approximation, denoted as $E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, of the expectation $E\left[{\mathsf{\Phi }}_H\right]$, is provided by the following truncated Taylor series:

$$E^{\left(2n\right)}\left[{\mathsf{\Phi }}_H\right]= E\left[{\mathsf{\Phi }}_H^{\left(2n\right)}\left(\sigma ;E\right)\right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}{\displaystyle \sum _{k=0}^n\frac{{\beta }^{2k}}{2k+1}}; E^{\left(2n+1\right)}\left[{\mathsf{\Phi }}_H\right]=E^{\left(2n\right)}\left[{\mathsf{\Phi }}_H\right]; n=0,1,\dots ; .$$

(48)

Notably, the convergence of the infinite series representation of $E\left[{\mathsf{\Phi }}_H\right]$ obtained in Equation (46) depends solely on the chosen value for the positive parameter $\beta$: the smaller the chosen value for $\beta$, the faster the convergence of the series representation for $E\left[{\mathsf{\Phi }}_H\right]$.

The variance, denoted as $V\left[{\mathsf{\Phi }}_H\right]$, of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, is obtained by recalling the following relation:

$$V\left[{\mathsf{\Phi }}_H\right]=E\left[{\left({\mathsf{\Phi }}_H\right)}^2\right]-E^2\left[{\mathsf{\Phi }}_H\right].$$

(49)

Using Equations (19) and (46) in Equation (49) yields the following result:

$$\begin{array}{l}V\left[{\mathsf{\Phi }}_H\right]={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H}\right]}^2\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}\frac{d{\sigma }_s\left(E\right)}{{\left[{\sigma }_s\left(E\right)\right]}^2}}-{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H}\frac{1}{2\beta {\sigma }_m}\ln \frac{1+\beta }{1-\beta }\right]}^2\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left[\frac{1}{1-{\beta }^2}-\frac{1}{4{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2\right].\end{array}$$

(50)

Note that $V\left[{\mathsf{\Phi }}_H\right]$ is an even function of $\beta$, i.e.,

$$V\left[{\mathsf{\Phi }}_H,\beta \right]=V\left[{\mathsf{\Phi }}_H,-\beta \right]={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left[\frac{1}{1-{\beta }^2}-\frac{1}{4{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2\right].$$

(51)

As indicated by the result obtained in Equation (50), the convergence of the series representation of $V\left[{\mathsf{\Phi }}_H\right]$ would also depend solely on the chosen value for the positive parameter $\beta$: the smaller the chosen value for $\beta$, the faster the convergence of the series representation for $V\left[{\mathsf{\Phi }}_H\right]$ would be.

The nth-order approximation, denoted as $V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ of the variance $V\left[{\mathsf{\Phi }}_H\right]$, is also obtained by using the relation provided in Equation (49), which takes on the following form:

$$V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]=E^{\left(n\right)}\left[{\left({\mathsf{\Phi }}_H^{\left(n\right)}\right)}^2\right]-{\left(E^{\left(n\right)}\right)}^2\left[{\mathsf{\Phi }}_H^{\left(n\right)}\right]; n=1,2,....; .$$

(52)

The third-order moment, denoted as ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$ of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, is obtained using Equations (19) and (46) in the definition of ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$ and by performing the ensuing integrations to obtain the result shown below:

$$\begin{array}{l}{\mu }_3\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{{\mathsf{\Phi }}_H\left(\sigma ;E\right)-E\left[{\mathsf{\Phi }}_H\right]\right\}}^{3}d{\sigma }_s\left(E\right)}\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^3\frac{1}{2\beta }\left[\frac{2\beta }{{\left(1-{\beta }^2\right)}^2}-\frac{3}{1-{\beta }^2}\ln \frac{1+\beta }{1-\beta }+\frac{1}{2{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^3\right].\end{array}$$

(53)

Note that ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$ is an even function of $\beta$, i.e.,

$$\begin{array}{l}{\mu }_3\left[{\mathsf{\Phi }}_H,\beta \right]={\mu }_3\left[{\mathsf{\Phi }}_H,-\beta \right]\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^3\frac{1}{2\beta }\left[\frac{2\beta }{{\left(1-{\beta }^2\right)}^2}-\frac{3}{1-{\beta }^2}\ln \frac{1+\beta }{1-\beta }+\frac{1}{2{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^3\right].\end{array}$$

(54)

The skewness of the distribution of the response ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ will be denoted as $Skew\left[{\mathsf{\Phi }}_H\right]$ and is defined as follows:

$$Skew\left[{\mathsf{\Phi }}_H\right]=\frac{{\mu }_3\left[{\mathsf{\Phi }}_H\right]}{{\left\{V\left[{\mathsf{\Phi }}_H\right]\right\}}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}=\frac{\frac{1}{2\beta }\left[\frac{2\beta }{{\left(1-{\beta }^2\right)}^2}-\frac{3}{1-{\beta }^2}\ln \frac{1+\beta }{1-\beta }+\frac{1}{2{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^3\right]}{{\left[\frac{1}{1-{\beta }^2}-\frac{1}{4{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2\right]}^{\frac{3}{2}}}.$$

(55)

Skewness indicates the direction and relative magnitude of a distribution’s deviation from the normal distribution. Note that, for this illustrative model, $Skew\left[{\mathsf{\Phi }}_H,\beta \right]=Skew\left[{\mathsf{\Phi }}_H,-\beta \right]$.

The nth-order approximation of the third-order moment, denoted as ${\mu }_3^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, is defined below:

$${\mu }_3^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ {\mathsf{\Phi }}_H^{\left(n\right)}\left(\sigma ;E\right)-E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{3}d{\sigma }_s\left(E\right)} .$$

(56)

It follows from Equations (56), (23) and (48) that the first-order approximation of the third-order moment, denoted as ${\mu }_3^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$, has the following expression:

$$\begin{array}{l}{\mu }_3^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ {\mathsf{\Phi }}_H^{\left(1\right)}\left(\sigma ;E\right)-E^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{3}d{\sigma }_s\left(E\right)}\\ =\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left[1-\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}\right]-\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right\}}^{3}d{\sigma }_s\left(E\right)}\equiv 0.\end{array}$$

(57)

It also follows from Equations (56), (23) and (48) that the second-order approximation of the third-order moment, denoted as ${\mu }_3^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$, has the following expression:

$$\begin{array}{l}{\mu }_3^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ {\mathsf{\Phi }}_H^{\left(2\right)}\left(\sigma ;E\right)-E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{3}d{\sigma }_s\left(E\right)}\\ =\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left[1-\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}+{\left(\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}\right)}^2\right]-\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1+\frac{{\beta }^2}{3}\right)\right\}}^{3}d{\sigma }_s\left(E\right)}\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^3\frac{4}{15}{\beta }^4+O\left({\beta }^6\right).\end{array}$$

(58)

It is evident that high-order approximations for the third-order moment are very laborious to obtain. On the other hand, since the definition of skewness involves not only the response’s third moment but also the response’s variance, it is possible to compute approximate values of the skewness, denoted as $Ske{w}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$, by using approximations of different orders for the response’s variance and third moment, respectively, as follows:

$$Ske{w}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{{\mu }_3^{\left(m\right)}\left[{\mathsf{\Phi }}_H\right]}{{\left\{V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}; m,n=1,2,\dots$$

(59)

The flexible definition for the approximate skewness $Ske{w}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$ provided in Equation (59) enables the use of higher-order approximations for the response variance, which require less computational effort to determine by comparison to the same order of approximation for the third moment. The flexible definition of $Ske{w}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$ could provide higher accuracy for the skewness with less computational effort by using a lower-order approximation for the third moment than the order of approximation used for computing the variance.

The fourth-order moment, denoted as ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$ of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, is defined as follows:

$$\begin{array}{l}{\mu }_4\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{{\mathsf{\Phi }}_H\left(\sigma ;E\right)-E\left[{\mathsf{\Phi }}_H\right]\right\}}^{4}d{\sigma }_s\left(E\right)}\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^4\left[\frac{1+{\beta }^2/3}{{\left(1-{\beta }^2\right)}^3}-\frac{2}{\beta {\left(1-{\beta }^2\right)}^2}\ln \frac{1+\beta }{1-\beta }+\frac{3}{2{\beta }^2\left(1-{\beta }^2\right)}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2-\frac{3}{16{\beta }^4}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^4\right] .\end{array}$$

(60)

Note that ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$ is an even function of $\beta$, i.e.,

$$\begin{array}{l}{\mu }_4\left[{\mathsf{\Phi }}_H,\beta \right]={\mu }_4\left[{\mathsf{\Phi }}_H,-\beta \right]\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^4\left[\frac{1+{\beta }^2/3}{{\left(1-{\beta }^2\right)}^3}-\frac{2}{\beta {\left(1-{\beta }^2\right)}^2}\ln \frac{1+\beta }{1-\beta }+\frac{3}{2{\beta }^2\left(1-{\beta }^2\right)}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2-\frac{3}{16{\beta }^4}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^4\right] .\end{array}$$

(61)

The nth-order approximation of the fourth-order moment, denoted as ${\mu }_4^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, is defined below:

$${\mu }_4^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ {\mathsf{\Phi }}_H^{\left(n\right)}\left(\sigma ;E\right)-E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{4}d{\sigma }_s\left(E\right)} .$$

(62)

It follows from Equations (62), (23) and (48) that the first-order approximation of the fourth-order moment, denoted as ${\mu }_4^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$, has the following expression:

$$\begin{array}{l}{\mu }_4^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ {\mathsf{\Phi }}_H^{\left(1\right)}\left(\sigma ;E\right)-E^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{4}d{\sigma }_s\left(E\right)}\\ =\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left[1-\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}\right]-\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right\}}^{4}d{\sigma }_s\left(E\right)}={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^4\frac{{\beta }^4}{5}.\end{array}$$

(63)

It also follows from Equations (62), (23) and (48) that the second-order approximation of the third-order moment, denoted as ${\mu }_4^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$, has the following expression:

$$\begin{array}{l}{\mu }_4^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ {\mathsf{\Phi }}_H^{\left(2\right)}\left(\sigma ;E\right)-E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^{4}d{\sigma }_s\left(E\right)}\\ =\frac{1}{2\beta {\sigma }_m}{\displaystyle \underset{\left(1-\beta \right){\sigma }_m}{\overset{\left(1+\beta \right){\sigma }_m}{\int }}{\left\{ \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left[1-\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}+{\left(\frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}\right)}^2\right]-\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1+\frac{{\beta }^2}{3}\right)\right\}}^{4}d{\sigma }_s\left(E\right)}\\ ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^3\left[\frac{{\beta }^4}{5}+\frac{6}{7}{\beta }^6+O\left({\beta }^8\right)\right].\end{array}$$

(64)

Evidently, high-order approximations for the fourth-order moment ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$ are very laborious to obtain.

The kurtosis of the distribution of the response ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ will be denoted as $Kurt\left[{\mathsf{\Phi }}_H\right]$ and is defined as follows:

$$\begin{array}{l}Kurt\left[{\mathsf{\Phi }}_H\right]=\frac{{\mu }_4\left[{\mathsf{\Phi }}_H\right]}{{\left\{V\left[{\mathsf{\Phi }}_H\right]\right\}}^2}\\ =\frac{\left[\frac{1+{\beta }^2/3}{{\left(1-{\beta }^2\right)}^3}-\frac{2}{\beta {\left(1-{\beta }^2\right)}^2}\ln \frac{1+\beta }{1-\beta }+\frac{3}{2{\beta }^2\left(1-{\beta }^2\right)}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2-\frac{3}{16{\beta }^4}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^4\right]}{{\left[\frac{1}{1-{\beta }^2}-\frac{1}{4{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2\right]}^2}.\end{array}$$

(65)

The kurtosis of a probability distribution indicates the propensity of the respective distribution to produce outliers or heavy tails.

The approximate kurtosis, denoted as $Kur{t}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$, can be computed by using approximations of different orders for the response’s variance and fourth moment, respectively, as follows:

$$Kur{t}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]=\frac{{\mu }_4^{\left(m\right)}\left[{\mathsf{\Phi }}_H\right]}{{\left\{V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\right\}}^2}; m,n=1,2,\dots .$$

(66)

The flexible definition for the approximate skewness $Kur{t}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$ provided in Equation (66) aims at obtaining higher accuracy for the kurtosis with less computational effort than would be required if the fourth moment were computed at the same high order $m=n$ as the variance. Thus, the variance would be computed relatively inexpensively at an order of approximation $n>m$, while the fourth moment would be computed only up to the order of approximation $m<n$, expecting that $Kur{t}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]\simeq Kur{t}^{\left(n,n\right)}\left[{\mathsf{\Phi }}_H\right]$ for a sufficiently high order $m$.

Section 3.1 has presented the investigation of the order of sensitivities that would be needed to achieve a target accuracy for predicting changes in the flux due to changes in the cross-section ${\sigma }_s\left(E\right)$ by using the Taylor-series expansion given in Equation (23). Analogously, it is of interest to investigate the accuracy provided by the approximations $E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, $V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, $Ske{w}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$ and $Kur{t}^{\left(m,n\right)}\left[{\mathsf{\Phi }}_H\right]$, $n=1,2,\dots ; m=1,2,\dots$, for the expectation $E\left[{\mathsf{\Phi }}_H\right]$, variance $V\left[{\mathsf{\Phi }}_H\right]$, skewness $Skew\left[{\mathsf{\Phi }}_H\right]$ and kurtosis $Kurt\left[{\mathsf{\Phi }}_H\right]$, respectively.

The expressions of the various orders of approximation for $E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ and $V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, $n=1,2,\dots$, will be computed by using the expression of the approximations ${\mathsf{\Phi }}_H^{\left(n\right)}$ in Equation (52), in order to verify independently the partial series that would be obtained by expanding (in truncated Taylor series in $\beta$) the exact results obtained in Equations (46) and (50). The following results are thus obtained:

• If only the first-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the first-order approximation ${\mathsf{\Phi }}_H^{\left(1\right)}\left(\sigma ;E\right)$ of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ takes on the following particular form of Equation (23):

  $${\mathsf{\Phi }}_H^{\left(1\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t\right); t\triangleq \frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}$$

  (67)

The first-order approximate expectation $E^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ has the following expression:

$$E^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}\left(1-t\right)} dt=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}=E^{\left(0\right)}\left[{\mathsf{\Phi }}_H\right]$$

(68)

The expectation $E^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate only to the zeroth-order term, since the contribution from the first-order sensitivities vanishes after the integration of the respective first-order (odd) term over the symmetric-interval around the expectation of the parameter ${\sigma }_s\left(E\right)$.

Replacing the expression given in Equation (67) into Equation (52) yields:

$$V^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left[\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}{\left(1-t\right)}^2}dt\right]-{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\frac{{\beta }^2}{3}$$

(69)

The variance $V^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate up to and including the contributions stemming from the first-order sensitivities.

2.

If first-order and second-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the second-order expansion ${\mathsf{\Phi }}_H^{\left(2\right)}\left(\sigma ;E\right)$ of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ takes on the following particular form of Equation (23):

$${\mathsf{\Phi }}_H^{\left(2\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t+t^2\right); t\triangleq \frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}.$$

(70)

The expectation $E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$ has the following expression:

$$E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}\left(1-t+t^2\right)} dt=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1+\frac{{\beta }^2}{3}\right).$$

(71)

The expectation $E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate up to and including the contributions stemming from the second-order sensitivities.

Replacing the expression given in Equation (70) into Equation (52) and using Equation (71) yields:

$$\begin{array}{ll}{V}^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]& ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left[\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}{\left(1-t+t^2\right)}^2}dt\right]-{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2{\left(1+\frac{{\beta }^2}{3}\right)}^2\\ & ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left(\frac{{\beta }^2}{3}+\frac{4}{45}{\beta }^4\right).\end{array}$$

(72)

The variance $V^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate up to and including the contributions stemming from the first-order sensitivities but the contributions from the second-order sensitivities are incomplete because the term $4{\beta }^4/45$ is incomplete. Additional contributions to the term $4{\beta }^4/45$ will arise from the third-order sensitivities, as will be shown below.

3.

If all sensitivities up to and including the third-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the third-order expansion ${\mathsf{\Phi }}_H^{\left(3\right)}\left(\sigma ;E\right)$ of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ takes on the following particular form of Equation (23):

$${\mathsf{\Phi }}_H^{\left(3\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t+t^2-t^3\right); t\triangleq \frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}.$$

(73)

The expectation $E^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]$ has, therefore, the following expression:

$$E^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}\left(1-t+t^2-t^3\right)} dt=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1+\frac{{\beta }^2}{3}\right)=E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right].$$

(74)

The expectation $E^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]$ is only as accurate as the lower-order approximation $E^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$, i.e., it is accurate up to and including the contributions stemming from the second-order sensitivities, because the contribution from the third-order sensitivities vanishes after the integration of the respective third-order (odd) term over the symmetric interval around the expectation of the parameter ${\sigma }_s\left(E\right)$.

Replacing the expression given in Equation (73) into Equation (52) and using Equation (74) yields:

$$\begin{array}{ll}{V}^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]& ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left[\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}{\left(1-t+t^2-t^3\right)}^2}dt\right]-{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2{\left(1+\frac{{\beta }^2}{3}\right)}^2\\ & ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left(\frac{{\beta }^2}{3}+\frac{22}{45}{\beta }^4+\frac{{\beta }^6}{7}\right).\end{array}$$

(75)

The variance $V^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate up to and including the contributions stemming from the second-order sensitivities but the contributions from the third-order sensitivities are incomplete because the term ${\beta }^6/7$ is incomplete. Additional contributions to the term ${\beta }^6/7$ will arise from the fourth-order sensitivities, as will be shown below.

4.

If all sensitivities up to and including the fourth-order sensitivities of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ with respect to ${\sigma }_s\left(E\right)$ are available, then the fourth-order expansion ${\mathsf{\Phi }}_H^{\left(4\right)}\left(\sigma ;E\right)$ of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ takes on the following particular form of Equation (23):

$${\mathsf{\Phi }}_H^{\left(4\right)}\left(\sigma ;E\right)\triangleq \frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1-t+t^2-t^3+t^4\right); t\triangleq \frac{{\sigma }_s\left(E\right)-{\sigma }_m}{{\sigma }_m}.$$

(76)

The expectation $E^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]$ has the following expression:

$$E^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}\left(1-t+t^2-t^3+t^4\right)} dt=\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\left(1+\frac{{\beta }^2}{3}+\frac{{\beta }^4}{5}\right).$$

(77)

The expectation $E^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate up to and including the contributions stemming from the fourth-order sensitivities.

Replacing the expression given in Equation (77) into Equation (52) and using Equation (77) yields:

$$\begin{array}{ll}{V}^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]& ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left[\frac{1}{2\beta } {\displaystyle \underset{-\beta }{\overset{\beta }{\int }}{\left(1-t+t^2-t^3+t^4\right)}^2}dt\right]-{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2{\left(1+\frac{{\beta }^2}{3}+\frac{{\beta }^4}{5}\right)}^2\\ & ={\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^2\left(\frac{{\beta }^2}{3}+\frac{22}{45}{\beta }^4+\frac{31}{105}{\beta }^6+\frac{16}{225}{\beta }^8\right).\end{array}$$

(78)

The variance $V^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]$ is accurate up to and including the contributions stemming from the third-order sensitivities but the contributions from the fourth-order sensitivities are incomplete because the term $16{\beta }^8/225$ is incomplete. Additional contributions to the term $16{\beta }^8/225$ will arise from the fifth-order sensitivities.

As has been mentioned in the foregoing and as indicated by the results obtained thus far, the convergence of the infinite series representations of the expectation and variance, $E\left[{\mathsf{\Phi }}_H\right]$ and, respectively, $V\left[{\mathsf{\Phi }}_H\right]$ depends solely on the chosen value for the positive parameter $\beta$: the smaller the chosen value for $\beta$, the faster the convergence of the series representation for $E\left[{\mathsf{\Phi }}_H\right]$ and for $V\left[{\mathsf{\Phi }}_H\right]$, respectively. It is, therefore, of interest to investigate numerically the errors associated with the first four orders of approximations $E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ and $V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ for $n=1,2,3,4$, for several representative values for $\beta$, $0<\beta <1$. Three values will be considered for the parameter $\beta$, as follows:

(a) $\beta =0.2$, which corresponds to a uniform distribution of height $P\left[{\sigma }_s\left(E\right)\right]=5/{\sigma }_m$ and (very narrow) width confined to the interval $0.8 {\sigma }_m\le {\sigma }_s\left(E\right)\le 1.2 {\sigma }_m$, with a standard deviation $SD\left[{\sigma }_s\left(E\right)\right]=\frac{1}{\sqrt{3}}\beta {\sigma }_m=\left(11.55\% \right){\sigma }_m$.

(b) $\beta =0.5$, which corresponds to a uniform distribution of height $P\left[{\sigma }_s\left(E\right)\right]=1/{\sigma }_m$ and width confined to the interval $0.5 {\sigma }_m\le {\sigma }_s\left(E\right)\le 1.5 {\sigma }_m$, with a standard deviation $SD\left[{\sigma }_s\left(E\right)\right]=\left(28.87\%\right){\sigma }_m$.

(c) $\beta =0.95$, which corresponds to a uniform distribution of height $P\left[{\sigma }_s\left(E\right)\right]=0.5263/{\sigma }_m$ and width confined to the interval $0.05 {\sigma }_m\le {\sigma }_s\left(E\right)\le 1.95 {\sigma }_m$, with a standard deviation $SD\left[{\sigma }_s\left(E\right)\right]=\left(54.85\%\right){\sigma }_m$.

The exact normalized expectation, $E_N\left[{\mathsf{\Phi }}_H\right]$, and normalized variance, $V_N\left[{\mathsf{\Phi }}_H\right]$, of the flux ${\mathsf{\Phi }}_H(\sigma ;E$ are defined by using Equations (46) and (50), as follows:

$$E_N\left[{\mathsf{\Phi }}_H\right]\triangleq E\left[{\mathsf{\Phi }}_H\right]{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H}\frac{1}{{\sigma }_m}\right]}^{-1}=\frac{1}{2\beta }\ln \frac{1+\beta }{1-\beta },$$

(79)

$$V_N\left[{\mathsf{\Phi }}_H\right]\triangleq V\left[{\mathsf{\Phi }}_H\right]{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^{-2}=\left[\frac{1}{1-{\beta }^2}-\frac{1}{4{\beta }^2}{\left(\ln \frac{1+\beta }{1-\beta }\right)}^2\right].$$

(80)

The normalized exact standard deviation is defined as follows: $S{D}_N\left[{\mathsf{\Phi }}_H\right]=\sqrt{V_N\left[{\mathsf{\Phi }}_H\right]}$. The normalized nth-order approximations $E_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$, $V_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ and $S{D}_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ of the expectation, variance and standard deviation of the flux response ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, respectively, are defined as follows:

$$\begin{array}{l}{E}_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq E^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H}\frac{1}{{\sigma }_m}\right]}^{-1}; V_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq V^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]{\left[\frac{Q\left({\alpha }^0\right)}{E{N}_H{\sigma }_m}\right]}^{-2};\\ S{D}_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]\triangleq \sqrt{V_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]} ; n=1,2,3,4.\end{array}$$

(81)

The results obtained from using the exact expression provided in Equations (79) and (80) are compared to the approximate results obtained by using the expressions provided in Equation (81), for $n=1,2,3,4$ and for $\beta =0.2$, $\beta =0.5$ and $\beta =0.95$, respectively. Since all four moments (i.e., expectation, variance, skewness and kurtosis) of the distribution of the response ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ are even functions of $\beta$, it suffices to consider only positive values for $\beta$. The results obtained for the expectation and variance of the predicted distribution of the response ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ are summarized in

Table 3. Exact and approximate values for the expectation, variance and standard deviation of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ for $\beta =0.2$.

Moment	Exact Value	First-Order Value (Error 1)	Second-Order Value (Error 1)	Third-Order Value (Error 1)	Fourth-Order Value (Error 1)
$E_N\left[{\mathsf{\Phi }}_H\right]$	1.014	1.000 (1.4%)	1.000 (1.4%)	1.013 (0.07%)	1.013 (0.07%)
$V_N\left[{\mathsf{\Phi }}_H\right]$	0.0142	0.0133 (5%)	0.0134 (3.88%)	0.0142 (exact 2)	0.0142 (exact 2)
$S{D}_N\left[{\mathsf{\Phi }}_H\right]$	0.1190	0.1155 (3%)	0.1158 (2.7%)	0.0190 (exact 2)	0.0190 (exact 2)

¹ Error = (exact-approximate)/exact. ² Exact to at least four decimal digits.

,

Table 4. Exact and approximate values for the expectation, variance and standard deviation of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ for $\beta =0.5$.

Moment	Exact Value	First-Order Value (Error 1)	Second-Order Value (Error 1)	Third-Order Value (Error 1)	Fourth-Order Value (Error 1)
$E_N\left[{\mathsf{\Phi }}_H\right]$	1.0990	1.000 (9.01%)	1.083 (1.43%)	1.083 (1.43%)	1.096 (0.27%)
$V_N\left[{\mathsf{\Phi }}_H\right]$	0.1264	0.0833 (34.05%)	0.0889 (29.67%)	0.1161 (8.13%)	0.1188 (6.03%)
$S{D}_N\left[{\mathsf{\Phi }}_H\right]$	0.3555	0.2887 (18.8%)	0.2981 (16.2%)	0.3375 (5.06%)	0.3447 (3.05%)

¹ Error = (exact-approximate)/exact.

and

Table 5. Exact and approximate values for the expectation, variance and standard deviation of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ for $\beta =0.95$.

Moment	Exact Value	First-Order Value (Error 1)	Second-Order Value (Error 1)	Third-Order Value (Error 1)	Fourth-Order Value (Error 1)
$E_N\left[{\mathsf{\Phi }}_H\right]$	1.928	1.000 (9.01%)	1.301 (32.53%)	1.301 (32.53%)	1.464 (24.1%)
$V_N\left[{\mathsf{\Phi }}_H\right]$	6.5385	0.3008 (95.40%)	0.3732 (94.29%)	0.8040 (87.70%)	0.9632 (85.27%)
$S{D}_N\left[{\mathsf{\Phi }}_H\right]$	2.5570	0.5485 (78.55%)	0.6109 (76.11%)	0.8967 (64.93%)	0.9814 (62.62%)

¹ Error = (exact-approximate)/exact.

, below.

The following conclusions can be drawn from the results presented in Table 3, Table 4 and Table 5:

• For very small values $\left(\beta =0.2\right)$ of the parameter $\beta$, the expansions (in powers of $\beta$), which represent the various orders $\left(n=1,2,3,4\right)$ of approximations of the normalized expectation $E_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ and normalized variance $V_N^{\left(n\right)}\left[{\mathsf{\Phi }}_H\right]$ of the unknown distribution of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, converge very quickly to the respective exact values $E_N\left[{\mathsf{\Phi }}_H\right]$ and $V_N\left[{\mathsf{\Phi }}_H\right]$, as indicated by the results presented in Table 3. Already the first-order approximations $E_N^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ and $V_N^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ yield results that are within 5% of the exact values of the expectation and, respectively, variance of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$. The third-order approximation $E_N^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ is within 0.07% of the exact value for the expectation of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, while the third-order approximation $V_N^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]$ of the standard deviation of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ is exact up to at least four decimals. Recall that a very small value of $\beta$ implies that the measurement ${\sigma }_m\triangleq {\sigma }_s\left(E_m\right)$ is extremely accurate, i.e., the uniform distribution of ${\sigma }_s\left(E_m\right)$ around ${\sigma }_m$ is extremely narrow.
• The value $\beta =0.5$ is characteristic of a measurement of moderate precision, so it is representative of the majority of evaluated cross-sections such as the scattering cross-section ${\sigma }_s\left(E_m\right)$. The results in Table 4 indicate that the first-order approximations are not satisfactory: the first-order approximation $E_N^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ of the exact expectation $E_N\left[{\mathsf{\Phi }}_H\right]$ is in error by 9%, while the first-order approximation $S{D}_N^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ of the exact standard deviation $S{D}_N\left[{\mathsf{\Phi }}_H\right]$ is in error by 19%. Only the third-order (or higher-order) approximations, $E_N^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]$ and $S{D}_N^{\left(3\right)}\left[{\mathsf{\Phi }}_H\right]$, would provide values for the expectation and, respectively, standard deviation of the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, which would be within 5% of the respective exact values.
• The value $\beta =0.95$ characterizes either an imprecise measurement or the need to construct “tolerance intervals” that cover a large segment of the unknown distribution of the quantity under investigation, i.e., the flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, in this illustrative example. As the results in Table 5 indicate, the convergence of the various approximations is extremely slow. Even the fourth-order approximations $E_N^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]$ and $S{D}_N^{\left(4\right)}\left[{\mathsf{\Phi }}_H\right]$ are very far off the exact values $E_N\left[{\mathsf{\Phi }}_H\right]$ and $S{D}_N\left[{\mathsf{\Phi }}_H\right]$, respectively, of the expectation and the standard deviation of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$.

The accuracies of the first-order and second-order approximations for the third-order moment, ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$, and the skewness, $Skew\left[{\mathsf{\Phi }}_H\right]$, of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ are illustrated by the results presented in

Table 6. Exact and approximate values for the third-order moment and skewness of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$.

Third-Order Moment	$\mathit{\beta }=0.2$	$\mathit{\beta }=0.5$	$\mathit{\beta }=0.95$
Exact: ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$	0.0005	0.0353	60.20
Exact: $Skew\left[{\mathsf{\Phi }}_H\right]$	0.2955	0.7857	3.6
First-order: ${\mu }_3^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ (Error 1)	0.0000 (100%)	0.0000 (100%)	0.0000 (100%)
First-order: $Ske{w}^{\left(1,1\right)}\left[{\mathsf{\Phi }}_H\right]$ (Error 1)	0.0000 (100%)	0.0000 (100%)	0.0000 (100%)
Second-order: ${\mu }_3^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$ (Error 1)	0.0004 (20%)	0.0267 (32.2%)	0.2172 (99.6%)
Second-order: $Ske{w}^{\left(2,2\right)}\left[{\mathsf{\Phi }}_H\right]$	0.2576	1.0073	0.9527
(Error 1)	(13.2%)	(−28.3%)	(73.6%)
$Ske{w}^{\left(2,4\right)}\left[{\mathsf{\Phi }}_H\right]$	0.2364	0.6521	0.2298
(Error 1)	(20%)	(17.0%)	(93.6%)

¹ Error = (exact-approximate)/exact.

. The exact moment, ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$, of the predicted distribution of the response ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ is positive, indicating a distribution, which would be skewed toward values more positive than the distribution’s expected value ${\sigma }_m$. As the positive values of the parameter $\beta$ increase towards unity, which would represent an increasingly broader distribution for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$, ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$ also increases; hence, the skewness of the distribution for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ also increases. With pronounced skewness, standard statistical inference procedures such as a confidence interval for a mean will be not only incorrect, in the sense that the true coverage level will differ from the nominal (e.g., 95%) level, but they will also result in unequal error probabilities on each side.

Very importantly, the results presented in bold in Table 6 indicate that the first-order approximation for the third moment ${\mu }_3^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ of the predicted response distribution is always zero. Thus, using just the first-order response sensitivities will always be 100% incorrect unless the distribution of the response under consideration happens to be symmetric. Conversely, if only first-order sensitivities are considered, the third-order moment of the response is always zero. Thus, the response distribution is always predicted to be symmetrical if only first-order sensitivities are considered. Hence, a “first-order sensitivity and uncertainty quantification” will always produce an erroneous third moment (and hence skewness) of the predicted response distribution, unless the distribution happens to be symmetrical. At least second-order sensitivities must be used in order to estimate the third-order moment (and, hence, the skewness) of the response distribution.

The overall conclusion that follows from the results presented in Table 6 is that low-order (i.e., first- and second-order) approximations for the third-order moment, ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$, of the response distribution can provide the correct sign (i.e., positive or negative) of the response skewness, but cannot provide reliably correct values for the skewness. Significantly higher-order approximations, which would require the computation of high-order response sensitivities to model parameters, are needed in order to obtain reliable values for the third-order response moment ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$ and, hence, for the skewness of the response distribution. Comparing the results in Table 6 for increasing (positive) values for the parameter $\beta$, the low-order (first- and second-order) approximations become increasingly unreliable as $\beta$ increases toward unity. This trend indicates that, as the response distribution becomes wider, sensitivities of an increasingly higher-order would be needed in order to improve the predicted values for the response skewness. Using a low-order (in Table 6: second-order) approximation for the third moment ${\mu }_3\left[{\mathsf{\Phi }}_H\right]$ but a higher-order (in Table 6: fourth-order) approximation for the variance (so that a comparable amount of computational effort would be needed for computing the respective approximations) does not significantly improve the final results. If accurate values for the response skewness are required, it is not possible to circumvent the need for computing high-order sensitivities, which in turn enable the computation of high-order approximations to the third-order moment and, hence, skewness.

The accuracies of the first-order and second-order approximations for the fourth-order moment ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$ of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ are illustrated by the results presented in

Table 7. Exact and approximate values for the fourth-order moment of ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$.

Fourth-Order Moment	$\mathit{\beta }=0.2$	$\mathit{\beta }=0.5$	$\mathit{\beta }=0.95$
Exact: ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$	0.0004	0.0411	779.47
Exact: $Kurt\left[{\mathsf{\Phi }}_H\right]$	1.9837	2.5725	18.23
First-order: ${\mu }_4^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$ (Error 1)	0.0003 (25%)	0.0125 (69.6%)	0.1629 (100%)
First-order: $Kur{t}^{\left(1,1\right)}\left[{\mathsf{\Phi }}_H\right]$ (Error 1)	1.6960 (14.5%)	1.8014 (29.97%)	1.8004 (90%)
Second-order: ${\mu }_4^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$ (Error 1)	0.0004 (Exact 2)	0.0259 (37%)	0.7930 (99.99%)
Second-order: $Kur{t}^{\left(2,2\right)}\left[{\mathsf{\Phi }}_H\right]$	2.2277	3.2771	5.694
(Error 1)	(−12.3%)	(−27.39)	(68.8%)
$Kur{t}^{\left(2,4\right)}\left[{\mathsf{\Phi }}_H\right]$	1.9837	1.8351	0.8548
(Error 1)	(Exact2)	(28.6%)	(95.3%)

¹ Error = (exact-approximate)/exact. ² Exact to at least four decimal digits.

. The results presented in Table 7 also indicate (similar to that indicated by the results presented in Table 6) that for small values of the parameter $\beta$, the first-order approximation, ${\mu }_4^{\left(1\right)}\left[{\mathsf{\Phi }}_H\right]$, provides an approximate value, which differs by 25% from the exact value ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$, while the second-order approximation ${\mu }_4^{\left(2\right)}\left[{\mathsf{\Phi }}_H\right]$ provides the exact first four significant digits of the exact value of ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$. The main conclusions that follow from the results presented in Table 7 are similar to those that were drawn from the results presented in Table 6 for the response’s third-order moment and skewness, namely: low-order (i.e., first- and second-order) approximations for the fourth-order moment, ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$, of the response distribution can provide the correct sign of the response kurtosis, but cannot provide reliably correct values for the kurtosis. Significantly higher-order approximations, which would require the computation of high-order response sensitivities to model parameters, are needed in order to obtain reliable values for the fourth-order response moment ${\mu }_4\left[{\mathsf{\Phi }}_H\right]$ and, hence, for the kurtosis of the response distribution. The results in Table 7 indicate that the low-order (first- and second-order) approximations become increasingly unreliable as $\beta$ increases (with positive values) toward unity, which means that, as the response distribution becomes wider, sensitivities of an increasingly higher-order would be needed in order to improve the predicted values for the response kurtosis. If accurate values for the response kurtosis are required, it is not possible to circumvent the need for computing high-order sensitivities.

It is known [18,19,20] that the form showed in Equation (18) of the slowing down neutron flux in a hydrogenous medium also holds for the asymptotic (below the cut-off energy of the fission source) slowing down flux in an infinite non-hydrogenous medium, where it takes on the following form:

$${\mathsf{\Phi }}_m\left(\alpha ;E\right)=\frac{Q\left(\alpha \right)}{\overline{\xi }\left(E\right)E{\mathsf{\Sigma }}_s\left(\alpha ;E\right)}.$$

(82)

In Equation (82), the quantity $\overline{\xi }\left(E\right)$ denotes the average increase in lethargy per collision and depends on the atomic/molecular masses and the scattering cross-sections of the nuclides in the medium. The analysis of the dependence of the flux ${\mathsf{\Phi }}_m\left(\alpha ;E\right)$ on the scattering cross-section of the mixture, ${\mathsf{\Sigma }}_s\left(\alpha ;E\right)$, is similar to the analysis of the slowing-down flux ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ in hydrogen [18,19,20].

The expression of the slowing down neutron flux distribution, $\mathsf{\Phi }\left(\alpha ;E\right)$, in a hydrogenous medium in which both scattering and absorption are important was provided in Equation (17). The Taylor-series expansion of $\mathsf{\Phi }\left(\alpha ;E\right)$ at the nominal parameter value ${\sigma }_m$ has the following form:

$$\mathsf{\Phi }\left({\sigma }^H;E\right)={\displaystyle \sum _{n=0}^{\infty }\frac{1}{n!}}{\left\{\frac{d^n\mathsf{\Phi }\left({\sigma }^H;E\right)}{d{\left({\sigma }^H\right)}^n}\right\}}_{\left(E_m;{\sigma }_m;{\alpha }^0\right)}{\left({\sigma }^H-{\sigma }_m\right)}^n,$$

(83)

where

$$\begin{array}{l}\frac{d^n\mathsf{\Phi }\left({\sigma }^H;E\right)}{d{\left({\sigma }^H\right)}^n}={\displaystyle \sum _{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right)}\left\{\frac{d^k}{d{\left({\sigma }^H\right)}^k}\left[\frac{1}{\mathsf{\Sigma }\left({\sigma }^H;E\right)}\right]\right\}\frac{d^{ n-k} G\left(\alpha ;E\right)}{d{\left({\sigma }^H\right)}^{n-k}} =G\left(\alpha ;E\right)\frac{d^{ n}}{d{\left({\sigma }^H\right)}^n}\left[\frac{1}{\mathsf{\Sigma }\left({\sigma }^H;E\right)}\right] \\ +\left(\begin{array}{c}n\\ 1\end{array}\right)\frac{d G\left(\alpha ;E\right)}{d\left({\sigma }^H\right)}\frac{d^{ n-1}}{d{\left({\sigma }^H\right)}^{n-1}}\left[\frac{1}{\mathsf{\Sigma }\left({\sigma }^H;E\right)}\right] +\dots +\frac{1}{\mathsf{\Sigma }\left({\sigma }^H;E\right)}\frac{d^{ n} G\left(\alpha ;E\right)}{d{\left({\sigma }^H\right)}^n};\\ G\left(\alpha ;E\right)\triangleq \frac{p\left(\alpha ;E\right)}{E}{\displaystyle {\int }_{E_s\left(\alpha \right)}^{\infty }Q\left(\alpha ;E^{\prime }\right)d{E}^{\prime } }\mathsf{\Phi }\left(\alpha ;E\right) .\\ \end{array}$$

(84)

It follows from Equations (83) and (84) that the first term, denoted as ${\mathsf{\Phi }}_1\left({\sigma }^H;E\right)$, in the Taylor series of $\mathsf{\Phi }\left({\sigma }^H;E\right)$ is the following infinite series:

$${\mathsf{\Phi }}_1\left({\sigma }^H;E\right)=G\left({\sigma }_m;{\alpha }^0;E_m\right){\displaystyle \sum _{n=0}^{\infty }\frac{1}{n!}}{\left\{\frac{d^{ n}}{d{\left({\sigma }^H\right)}^n}\left[\frac{1}{\mathsf{\Sigma }\left({\sigma }^H;E\right)}\right]\right\}}_{\left(E_m;{\sigma }_m;{\alpha }^0\right)}{\left({\sigma }^H-{\sigma }_m\right)}^n,$$

(85)

which has the same functional form as the series for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ in Equation (20). In fact, each successive term in the series representation of $\mathsf{\Phi }\left(\alpha ;E\right)$ shown in Equation (83) has the same functional form as the series for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ in Equation (20). Therefore, the analysis that has been performed in the foregoing for ${\mathsf{\Phi }}_H\left(\sigma ;E\right)$ is also applicable to perform sensitivity and uncertainty analyses aimed at quantifying the impact of the uncertainties afflicting the other parameters (cross-sections, number densities, source parameters, etc.) on the distribution of the slowing down flux $\mathsf{\Phi }\left(\alpha ;E\right)$ in a scattering and absorbing medium, but these additional analyses are beyond the scope of this work.

4. Use of Model Response Sensitivities in First-Order Data Adjustment and High-Order Predictive Modeling Methodologies

This section reviews the incorporation of sensitivities (of various orders) in data adjustment and predictive modeling methodologies in order to indicate the level of accuracy that could be expected when using these methodologies. Section 4.1 reviews the first-order “generalized linear least squares adjustment (GSSLA)” methodology [2] incorporated in the data adjustment software system TSURFER [1], developed at Oak Ridge National Laboratory (ORNL), which represents the current state of the art in data adjustment used in nuclear engineering for nuclear criticality safety and reactor physics/shielding applications. Section 4.2 compares the TSURFER-GSSLA data adjustment methodology to the high-order predictive modeling (HO-BERRU-PM) methodology recently developed by Cacuci [13,21], highlighting the current capabilities and limitations of these two methodologies.

4.1. First-Order Generalized Least Squares Data Adjustment (GLLSA) Methodology in TSURFER

The current state of the art in “data adjustment” is represented by the so-called “generalized linear least squares adjustment (GLLSA)” procedure [2] implemented in the software module TSURFER [1] of the ORNL-SCALE software system used for reactor criticality safety applications. The GSSLA procedure in TSURFER is based on concepts proposed in the 1970s [25,26,27] and is limited to incorporating just the first-order sensitivities of the response with respect to the model parameters. The end results for the “adjusted” values produced by the TSURFER-GLLSA procedure for the nuclear cross-sections (considered as “model parameters”) and effective multiplication factors (considered as “model responses”) are obtained by minimizing a (user-chosen) quadratic cost functional, which represents the discrepancies (“least-squares”) between measurements and computations. The GSSLA procedure in TSURFER considers the following a priori information, which is described in the paragraphs labeled (A) through (C), below:

(A). A set $\left\{{\alpha }_n|n=1,\dots ,M\right\}$ of $M$ imprecisely known “model parameters”, which in TSURFER are the infinitely dilute multigroup cross-sections for all nuclide-reaction pairs used in the transport calculations of all responses. The total number of model parameters, $M$, is obtained by multiplying the number of unique nuclide-reaction pairs by the number of energy groups. In practice, the values of the cross-sections ${\alpha }_n$ are determined experimentally and are, therefore, imprecisely known. For practical computations, these cross-sections are arranged into an $M$-dimensional column vector, denoted as $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_M\right)}^{†}$, and are considered to be variates that obey an unknown multivariate probability distribution function, denoted as $p_{\alpha }\left(\alpha \right)$, which is formally defined on a domain $D_{\alpha }$. The dagger superscript $(†)$ will be used in this work to denote “transposition”. The various moments (e.g., mean values, covariances and variances, etc.) of $p_{\alpha }\left(\alpha \right)$ are determined from cross-section measurements and subsequent evaluation processes and are formally defined as follows:

(i) The expected (or mean) value of a model parameter ${\alpha }_n$, denoted as $E\left({\alpha }_n\right)$, is formally defined as follows:

$$E\left({\alpha }_n\right)\triangleq {\displaystyle \underset{D_{\alpha }}{\int }{\alpha }_n{p}_{\alpha }\left(\alpha \right)d\alpha }, n=1,\dots ,M.$$

(86)

The vector of expected values of the vector of parameters $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_M\right)}^{†}$ is defined as follows:

$$E\left[\alpha \right]\triangleq {\left[E\left({\alpha }_1\right),\dots ,E\left({\alpha }_M\right)\right]}^{†}.$$

(87)

(ii) The covariance, denoted as $\mathrm{cov}\left({\alpha }_j,{\alpha }_n\right)$, of two parameters ${\alpha }_j$ and ${\alpha }_n$ is formally defined as follows:

$$\mathrm{cov}\left({\alpha }_j,{\alpha }_n\right)\triangleq {\displaystyle \underset{D_{\alpha }}{\int }\left(\delta {\alpha }_j\right)\left(\delta {\alpha }_n\right)p_{\alpha }\left(\alpha \right)d\alpha } ; \delta {\alpha }_n\triangleq {\alpha }_n-E\left({\alpha }_n\right); j,n=1,\dots ,M.$$

(88)

(iii) The variance, denoted as $\mathrm{var}\left({\alpha }_n\right)$, of a parameter ${\alpha }_n$, is defined as follows:

$$\mathrm{var}\left({\alpha }_n\right)\triangleq {\displaystyle \underset{D_{\alpha }}{\int }{\left[{\alpha }_n-E\left({\alpha }_n\right)\right]}^2{p}_{\alpha }\left(\alpha \right)d\alpha }, n=1,\dots ,M.$$

(89)

(iv) The standard deviation ${\sigma }_n$ of ${\alpha }_n$ is defined as follows:

$${\sigma }_n\triangleq \sqrt{\mathrm{var}({\alpha }_n)}, n=1,\dots ,M.$$

(90)

(v) The correlation, denoted as ${\rho }_{jn}$, between two parameters ${\alpha }_m$ and ${\alpha }_n$, is defined as follows:

$${\rho }_{jn}\triangleq \mathrm{cov}\left({\alpha }_j,{\alpha }_n\right)/\left({\sigma }_j{\sigma }_n\right); j,n=1,\dots ,M.$$

(91)

(vi) In TSURFER, the cross-section moments defined in Equations (86)–(89) are used to define “relative” quantities as follows:

$$C_{\alpha \alpha }\left(m,n\right)\triangleq \left[\frac{\mathrm{cov}\left({\alpha }_m,{\alpha }_n\right)}{E\left({\alpha }_m\right)E\left({\alpha }_n\right)}\right] , m,n=1,\dots ,M.$$

(92)

Note that the quantities $C_{\alpha \alpha }\left(m,n\right)$ are not correlation coefficients, since the correlation coefficients (in a correlation matrix) are defined by dividing the absolute covariances by the respective standard deviations, as defined in Equation (91), rather than by the parameter expected values, as done in TSURFER and as (therefore) shown in Equation (92). The elements $C_{\alpha \alpha }\left(m,n\right)$ are used in TSURFER to construct an $M\times M$-dimensional cross-section covariance matrix, denoted as $C_{\alpha \alpha }$, and defined as follows:

$$C_{\alpha \alpha }\triangleq {\left[C_{\alpha \alpha }\left(m,n\right)\right]}_{M\times M}.$$

(93)

(B) TSURFER also considers a set $\left\{m_i|i=1,\dots ,I\right\}$ of $I$ “integral response” (e.g., effective multiplication factors, reaction rates), some of which may have been measured in selected benchmark experiments. The set of measured responses is arranged in an $I$-dimensional column vector denoted as $m\triangleq {\left(m_1,\dots m_i,\dots ,m_I\right)}^{†}$. The measured integral responses are imprecisely known quantities, which are characterized by their expected values and variances/covariances. The expected values $E\left(m_i\right)$ are formally defined as follows:

$$E\left(m_i\right)\triangleq {\displaystyle \underset{D_m}{\int }{m}_i{p}_m\left(m\right)dm}, i=1,\dots ,I,$$

(94)

where $p_m\left(m\right)$ is an unknown multivariate distribution for the measured responses, formally defined on a domain $D_m$. In TSURFER’s “absolute format”, the $I\times I$ covariance matrix for the measured responses is denoted as ${\tilde{C}}_{mm}$ and comprises elements ${\tilde{C}}_{mm}\left(i,j\right)$ defined as follows:

$${\tilde{C}}_{mm}\left(i,j\right)\triangleq {\displaystyle \underset{D_m}{\int }\left[m_i-E\left(m_i\right)\right]\left[m_j-E\left(m_j\right)\right]p_m\left(m\right)dm}, i,j=1,\dots ,I.$$

(95)

(C) Finally, TSURFER considers a set $\left\{k_i\left(\alpha \right)|i=1,\dots ,I\right\}$ of $I$ “calculated responses” (which represent either effective multiplication factors or reaction rates), which correspond to each experimental “integral response”. The notation $k_i\left(\alpha \right)$ indicates that these calculated responses are functions of the model parameters (i.e., nuclear data). The calculated responses $k_i\left(\alpha \right)$ are arranged in an $I$-dimension vector $k\left(\alpha \right)\triangleq {\left[k_1\left(\alpha \right),\dots ,k_I\left(\alpha \right)\right]}^{†}$. Computed responses, which have no measured counterparts, are designated in TSURFER as “applications”. By convention, unknown experimental values corresponding to applications are represented by the corresponding calculated values. Each calculated response, $k_i\left(\alpha \right)$, is considered in TSURFER to be represented by the following first-order (linear) Taylor-series expansion:

$$k_i\left(\alpha \right)=k_i\left({\alpha }^{ref}\right)+{{\displaystyle \sum _{j=1}^M\left[\frac{\partial k_i\left(\alpha \right)}{\partial {\alpha }_j}\right]}}_{\left({\alpha }_j^{ref}\right)}\left({\alpha }_j-{\alpha }_j^{ref}\right), i=1,\dots ,I,$$

(96)

where ${\alpha }^{ref}$ is a vector that comprises “cross-section reference values”. In particular, setting ${\alpha }_n^{ref}\equiv E\left({\alpha }_n\right)$ in Equation (96) and taking the expectation value of the resulting expression yields the following expression for the expected value of $k_i\left(\alpha \right)$:

$$E\left[k_i\left(\alpha \right)\right]=k_i\left({\alpha }^{ref}\right); {\alpha }^{ref}\triangleq \left[E\left(\alpha \right)\right], i=1,\dots ,I.$$

(97)

The relation in Equation (97) indicates that the expected value of the calculated response in TSURFER is the same as the TSURFER-calculated response itself. As has been shown in Section 3 and will also be shown in Section 4.2, however, the assumption made in TSURFER that Equation (97) is correct holds only if the second- and higher-order terms in the Taylor series of $k_i\left(\alpha \right)$ are discarded. Of course, Equation (97) is invalid when second (and/or higher) order sensitivities are not discarded.

The linear first-order expression provided in Equation (96) is also used in TSURFER to obtain expressions for the variances/covariances, $\mathrm{cov}\left[k_i\left(\alpha \right),k_j\left(\alpha \right)\right]$, between two computed responses $k_i\left(\alpha \right)$ and $k_j\left(\alpha \right)$. These covariances are obtained by writing Equation (96) for $k_j\left(\alpha \right)$, multiplying the Taylor series for $k_j\left(\alpha \right)$ by the Taylor series for $k_i\left(\alpha \right)$ and formally determining the expectation value of the expression that results from this multiplication. The expression obtained after this sequence of operations is as follows:

$$\mathrm{cov}\left[k_i\left(\alpha \right),k_j\left(\alpha \right)\right]={\displaystyle \sum _{m=1}^M{{\displaystyle \sum _{n=1}^M\left[\frac{\partial k_i\left(\alpha \right)}{\partial {\alpha }_m}\frac{\partial k_j\left(\alpha \right)}{\partial {\alpha }_n}\right]}}_{E\left(\alpha \right)}\mathrm{cov}\left({\alpha }_m,{\alpha }_n\right)}, i,j=1,\dots ,I,$$

(98)

where the symbol ${[ ]}_{E\left(\alpha \right)}$ indicates “evaluation at the expected values of the model parameters” (cross-sections). The quantity $\partial k_i\left(\alpha \right)/\partial {\alpha }_m$ represents the “absolute” first-order sensitivity (i.e., functional derivative) of the ith response with respect to the mth parameter and has units of “[response]/[parameter]”. The unitless quantities $\left(\partial k_i/\partial {\alpha }_m\right)\times \left({\alpha }_m/k_i\right)$ are called “relative” sensitivities.

Using the quantities defined in Equation (92), TSURFER re-writes Equation (98) in the following form:

$$\begin{array}{l}\mathrm{cov}\left[k_i\left(\alpha \right),k_j\left(\alpha \right)\right]={\displaystyle \sum _{m=1}^M{{\displaystyle \sum _{n=1}^M\left[\left({\alpha }_m\frac{\partial k_i\left(\alpha \right)}{\partial {\alpha }_m}\right)\left({\alpha }_n\frac{\partial k_j\left(\alpha \right)}{\partial {\alpha }_n}\right)\right]}}_{E\left(\alpha \right)}{C}_{\alpha \alpha }\left(m,n\right)}, \\ ={\displaystyle \sum _{m=1}^M{{\displaystyle \sum _{n=1}^M\left[{\tilde{S}}_{k\alpha }\left(i,m\right){\tilde{S}}_{k\alpha }\left(j,n\right)\right]}}_{E\left(\alpha \right)}{C}_{\alpha \alpha }\left(m,n\right)} , i,j=1,\dots ,I,\end{array}$$

(99)

where

$${\tilde{S}}_{k\alpha }\left(i,m\right)\triangleq {\alpha }_m\frac{\partial k_i\left(\alpha \right)}{\partial {\alpha }_m}.$$

(100)

The quantities ${\tilde{S}}_{k\alpha }\left(i,m\right)$ are used in TSURFER to construct the $I\times M$ matrix ${\tilde{S}}_{k\alpha }$, defined as follows:

$${\tilde{S}}_{k\alpha }\triangleq {\left[{\tilde{S}}_{k\alpha }\left(i,m\right)\right]}_{I\times M}.$$

(101)

Although TSURFER calls the matrix ${\tilde{S}}_{k\alpha }$ to be the “I by M absolute sensitivity matrix”, this is a misnomer, since the elements ${\tilde{S}}_{k\alpha }\left(i,m\right)$ of this matrix have the same units as the respective computed response and are, therefore, neither “absolute” sensitivities (which the quantities $\partial k_i/\partial {\alpha }_m$ are) nor “relative” sensitivities (which the unitless quantities $\left(\partial k_i/\partial {\alpha }_m\right)\times \left({\alpha }_m/k_i\right)$ are).

In matrix form, Equation (99) is written in TSURFER in the following form when evaluated at the expected parameter values $E\left[\alpha \right]$:

$${\tilde{C}}_{kk}\triangleq {\left[\mathrm{cov}\left(k_i,k_j\right)\right]}_{I\times I} ={\tilde{S}}_{k\alpha }{C}_{\alpha \alpha }{\tilde{S}}_{k\alpha }^{†}$$

(102)

The matrix ${\tilde{C}}_{kk}$ is called in TSURFER the “I by I absolute covariance matrix for prior calculated responses”, and Equation (102) is called the “sandwich rule”.

The GLLSA procedure [2] implemented in TSURFER “consolidates” the a priori information provided above in items (A) through (C) in order to obtain “adjusted” (i.e., “improved estimates of”) mean values and covariances for the parameters (cross-sections), as well as “adjusted” values for the means and covariances of the responses and the (target) application under consideration. The GLLSA procedure computes the adjusted values for the means and covariances of the parameters and responses by performing a constrained minimization of the user-chosen “chi-square functional” defined below

$${\chi }^2\triangleq {\left[\frac{{\alpha }^{\prime }-\alpha }{\alpha }\right]}^{†}{C}_{\alpha \alpha }^{-1}\left[\frac{{\alpha }^{\prime }-\alpha }{\alpha }\right]+{\left[\frac{m^{\prime }-m}{m}\right]}^{†}{C}_{mm}^{-1}\left[\frac{m^{\prime }-m}{m}\right],$$

(103)

and subject to the constraint that the linear model represented by Equation (96) be satisfied. In Equation (103), the matrix $C_{mm}$ represents the $I\times I$-dimensional relative covariance matrix for the prior measured responses, having the components $C_{mm}\left(i,j\right)\triangleq {\tilde{C}}_{mm}\left(i,j\right)/\left(m_i{m}_j\right),$$i,j=1,\dots ,I$. Furthermore, the quantities ${\alpha }^{\prime }$ and $m^{\prime }$ in Equation (103) represent the “adjusted values” for the model parameters and responses, which are to be determined as the end results of performing the constrained minimization procedure. This minimization procedure is performed by constructing an augmented Lagrangian functional, which is obtained by using Lagrange multipliers to append the “hard constraints” represented by Equation (96) to the user-defined “chi-square functional” shown in Equation (103) and, subsequently, determining the values/expressions for ${\alpha }^{\prime }$ and $m^{\prime }$ at which the first-order derivatives of the resulting augmented Lagrangian vanish.

The TSURFER formulas for the adjusted quantities that result from the GLLSA procedure are reproduced below, written in what TSURFER calls “absolute format”:

• The formula for computing the $M$-dimensional vector of adjusted cross-sections (“model parameters”), denoted in TSURFER as ${\alpha }^{\prime }$, is as follows:

  $${\alpha }^{\prime }=\alpha -\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]{\tilde{C}}_{dd}^{-1}\tilde{d}, \Delta \tilde{\alpha }\equiv {\alpha }^{\prime }-\alpha ,$$

  (104)

  where:

  $${\tilde{C}}_{dd}\triangleq {\tilde{C}}_{kk}+{\tilde{C}}_{mm}={\tilde{S}}_{k\alpha }{C}_{\alpha \alpha }{\tilde{S}}_{k\alpha }^{†}+{\tilde{C}}_{mm}$$

  (105)

  Additionally,

  $$\tilde{d}\triangleq k\left(\alpha \right)-m$$

  (106)
• The formula for computing the $M\times M$-dimensional covariance matrix for the adjusted cross-sections, denoted in TSURFER as $C_{{\alpha }^{\prime }{\alpha }^{\prime }}$, is as follows:

  $$C_{{\alpha }^{\prime }{\alpha }^{\prime }}=C_{\alpha \alpha }-\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]{\tilde{C}}_{dd}^{-1}\left[{\tilde{S}}_{k\alpha }{C}_{\alpha \alpha }\right].$$

  (107)
• The $M$-dimensional vector of adjusted measured responses produced by TSURFER’s GLLSA procedure is denoted as $m^{\prime }$, and the formula for computing it is as follows:

  $$m^{\prime }=m+{\tilde{C}}_{mm}{\tilde{C}}_{dd}^{-1}\tilde{d}; \Delta \tilde{m}\triangleq m^{\prime }-m.$$

  (108)
• The formula for computing the covariance matrix ${\tilde{C}}_{m^{\prime }{m}^{\prime }}$ for the adjusted measured responses is as follows:

  $${\tilde{C}}_{m^{\prime }{m}^{\prime }}={\tilde{C}}_{mm}-{\tilde{C}}_{mm}{\tilde{C}}_{dd}^{-1}{\tilde{C}}_{mm} .$$

  (109)
• The GLLSA procedure in TSURFER forces the following relationships to hold:

  $$k^{\prime }\left({\alpha }^{\prime }\right)=m^{\prime },$$

  (110)

  where $m^{\prime }$ denotes the $I$-dimensional vector of adjusted measured responses produced by the GLLSA procedure and $k^{\prime }\left({\alpha }^{\prime }\right)$ denotes the $I$-dimensional vector of adjusted calculated responses obtained using the adjusted nuclear parameters ${\alpha }^{\prime }$. Consequently, the following relationships are also forced to hold in TSURFER:

  $$\Delta k=k^{\prime }\left({\alpha }^{\prime }\right)-k\left(\alpha \right)=m^{\prime }-m-d=S_{k\alpha }\Delta \alpha .$$

  (111)
• TSURFER also computes the “consistency indicator between the calculations and measurements”, which is given in TSURFER by the following “chi-square formula”:

  $${\chi }_{TSURFER}^2=d^{†}{C}_{dd}^{-}d .$$

  (112)

The application response bias in TSURFER is denoted as ${\beta }_{\alpha }$ and defined as the expected deviation of the original calculated application response, denoted as $k_a\left(\alpha \right)$, from the best estimate of the measured response, denoted as ${m^{\prime }}_a$, which is unknown but is assumed to obey some probability distribution. If the application response did actually have a prior measured value $m_a$, then the best estimate for the experiment value would be the adjusted value ${m^{\prime }}_a$ provided by the GLLSA procedure. Mathematically, TSURFER defines the application bias as the expectation of the difference between $k_a\left(\alpha \right)$ and ${m^{\prime }}_a$, i.e.,

$${\beta }_{\alpha }\triangleq E\left[k_a\left(\alpha \right)-{m^{\prime }}_a\right]\cong k_a-k_a\left({\alpha }^{\prime }\right)\cong -{\tilde{S}}_a^{†}\Delta \alpha .$$

(113)

The first-order linear relation shown in Equation (113) is used in TSURFER as the basis for constructing “confidence intervals” and “tolerance limits” for the predicted multiplication factors (predicted model responses) for which there are no measurements. TSURFER assumes a first-order approximation of the dependence of the effective multiplication factors on the underlying cross-sections and, therefore, assumes a linear relationship between the initial experimental biases $\Delta k_i^{(\exp )}$ for the ith experiment and the TSURFER calculated bias $\Delta k^{(\mathrm{app})}$ for the application. Furthermore, TSURFER assumes that all distributions are Gaussian and that the standard deviations for the experimental measurements and reference calculated values are known.

Of all of the above assumptions made in TSURFER for constructing “confidence intervals” and “tolerance limits” for the predicted multiplication factors, the most questionable (and severe) is the assumption of linearity shown in Equation (113), which holds only if all sensitivities higher than the first order are discarded. However, TSURFER cannot compute second- (or higher) order sensitivities, so the validity of Equation (113) is unproven. Evidently, if second- and/or higher-order sensitivities were not negligible, Equation (113) would be severely in error, and consequently, the “confidence intervals” and “tolerance limits” computed in TSURFER would be severely misleading. In conclusion, the “confidence intervals” and “tolerance limits” computed in TSURFER are not trustworthy until TSURFER implements capabilities for computing higher-order sensitivities and proves quantitatively that these are negligibly small. Otherwise, all of the uncertainty analysis and data adjustment results, including tolerance limits and confidence intervals, produced by TSURFER must be interpreted with due caution.

4.2. High-Order Sensitivities in the HO-BERRU-PM Methodology

The HO-BERRU-PM methodology (high-order best-estimated results with reduced uncertainties predictive modeling) conceived by Cacuci [13,21] combines experimental and computational information in the joint phase space of responses and model parameters, including second-order and third-order response sensitivities to model parameters. The HO-BERRU-PM methodology uses the maximum entropy principle [22] to eliminate the need for introducing and “minimizing” a user-chosen “cost functional quantifying the discrepancies between measurements and computations”, thus yielding results that are free of subjective user-interferences, while generalizing and significantly extending the current dynamic data assimilation procedures [28,29]. Incorporating correlations, including those between the imprecisely known model parameters and computed model responses, the HO-BERRU-PM also provides a quantitative metric, constructed from sensitivity and covariance matrices, for determining the degree of agreement among the various computational and experimental data while eliminating discrepant information.

In contradistinction to TSURFER, which uses the linear (first-order) model shown in Equation (96), Cacuci’s HO-BERRU-PM framework [13,21] uses the third-order multivariate Taylor-series expansion of the model’s computed response with respect to the model’s parameters. Denoting the computed model response as $r_{i_1}^c\left(\alpha \right)$, where the superscript “c” denotes “computed” and where the subscript $i_1=1,\dots ,N_r$ denotes one of a total of $N_r$ responses that would be of interest, the third-order multivariate Taylor-series expansion of the response around the model’s parameters’ expected values ${\alpha }^0\triangleq \left({\alpha }_1^0,\dots ,{\alpha }_{N_{\alpha }}^0\right)$, where $N_{\alpha }$ denotes the total number of model parameters, has the following form:

$$r_{i_1}^c\left(\alpha \right)=r_{i_1}^c\left({\alpha }^0\right)+{{\displaystyle \sum _{i=1}^{N_{\alpha }}\left\{\frac{\partial r_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i}\right\}}}_{{\alpha }^0}\left({\alpha }_i-{\alpha }_i^0\right)+r_{i_1}^{HO}\left(\alpha ;{\alpha }^0\right); i_1=1,\dots ,N_r,$$

(114)

where

$$\begin{array}{l}{r}_{i_1}^{HO}\left(\alpha ;{\alpha }^0\right)\triangleq \frac{1}{2}{\displaystyle \sum _{i,j=1}^{N_{\alpha }}{\left\{\frac{{\partial }^2{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\right\}}_{{\alpha }^0}\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)}\\ +\frac{1}{3!}{\displaystyle \sum _{i,j,k=1}^{N_{\alpha }}{\left\{\frac{{\partial }^3{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j\partial {\alpha }_k}\right\}}_{{\alpha }^0}\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_k-{\alpha }_k^0\right)}\\ + \frac{1}{4!}{\displaystyle \sum _{i,j,k,\ell =1}^{N_{\alpha }}{\left\{\frac{{\partial }^4{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j\partial {\alpha }_k\partial {\alpha }_{\ell }}\right\}}_{{\alpha }^0}\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_k-{\alpha }_k^0\right)\left({\alpha }_{\ell }-{\alpha }_{\ell }^0\right)}; i_1=1,\dots ,N_r.\end{array}$$

(115)

In Equations (114) and (115), the expected parameter values are denoted as ${\alpha }^0\triangleq \left({\alpha }_1^0,\dots ,{\alpha }_{N_{\alpha }}^0\right)$, while the quantity $r_{i_1}^c\left({\alpha }^0\right)$ and the response derivatives (a.k.a. “sensitivities”) with respect to the model parameters are evaluated at the expected parameter values.

The expected value of the third-order Taylor-series expansion provided in Equation (114) is obtained by integrating this expansion formally over the unknown parameter distribution $p_{\alpha }\left(\alpha \right)$. This operation yields the following expression for the expected (mean) value, denoted as $E\left[r_{i_1}^c\left(\alpha \right)\right]$, of a response $r_{i_1}^c\left(\alpha \right)$:

$$E\left[r_{i_1}^c\left(\alpha \right)\right]=r_{i_1}^c\left({\alpha }^0\right)+\frac{1}{2}{\displaystyle \sum _{i,j=1}^{N_{\alpha }}{\left\{\frac{{\partial }^2{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\right\}}_{{\alpha }^0}{\rho }_{ij}}{\sigma }_i{\sigma }_j+\frac{1}{6}{\displaystyle \sum _{i,j,k=1}^{N_{\alpha }}{\left\{\frac{{\partial }^3{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j\partial {\alpha }_k}\right\}}_{{\alpha }^0}{t}_{ijk}{\sigma }_i{\sigma }_j{\sigma }_k}; i_1=1,\dots ,N_r.$$

(116)

The third-order term in Equation (116) involves implicitly the third-order moment,${\mu }_3^{ijk}$, of the multivariate parameter distribution function $p_{\alpha }\left(\alpha \right)$, which is written explicitly in terms of the third-order correlation parameter, $t_{ijk}$; these third-order quantities are defined as follows:

$${\mu }_3^{ijk}\left(\alpha \right)\triangleq {\displaystyle \underset{D_{\alpha }}{\int }\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_k-{\alpha }_k^0\right)p_{\alpha }\left(\alpha \right)d\alpha }\triangleq t_{ijk}{\sigma }_i{\sigma }_j{\sigma }_k; i,j,k=1,\dots ,N_{\alpha };$$

(117)

$$q_{ijkl}{\sigma }_i{\sigma }_j{\sigma }_k{\sigma }_l\triangleq {\langle \left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_k-{\alpha }_k^0\right)\left({\alpha }_l-{\alpha }_l^0\right)\rangle }_{\alpha }\triangleq {\mu }_4^{ijkl}; i,j,k,l=1,\dots ,J_{\alpha }$$

(118)

Comparing the expression in Equation (116) with the expression used in TSURFER for the expected value of the calculated response under consideration, which is given in Equation (97), indicates that for a response $k_i\left(\alpha \right)\equiv r_i^c\left(\alpha \right)$, the HO-BERRU-PM expected value contains higher-order terms in addition to the zeroth-order term $k_i\left({\alpha }^0\right)$ considered by TSURFER, namely:

$$E^{BERRU}\left[r_i^c\left(\alpha \right)\right]=E^{TSURFER}\left[r_i^c\left(\alpha \right)\right]+E^{HO}\left[r_i^c\left({\alpha }^0\right)\right] ; i=1,\dots ,I,$$

(119)

where

$$E^{TSURFER}\left[r_i^c\left(\alpha \right)\right]=r_i^c\left({\alpha }^0\right),$$

(120)

Additionally,

$$E^{HO}\left[r_i^c\left(\alpha \right)\right]\triangleq \frac{1}{2}{\displaystyle \sum _{i,j=1}^{N_{\alpha }}{\left\{\frac{{\partial }^2{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\right\}}_{{\alpha }^0}{\rho }_{ij}}{\sigma }_i{\sigma }_j+\frac{1}{6}{\displaystyle \sum _{i,j,k=1}^{N_{\alpha }}{\left\{\frac{{\partial }^3{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j\partial {\alpha }_k}\right\}}_{{\alpha }^0}{t}_{ijk}{\sigma }_i{\sigma }_j{\sigma }_k}; i=1,\dots ,N_r.$$

(121)

Multiplying the third-order Taylor expansion for the response $r_{i_1}^c\left(\alpha \right)$ provided in Equations (114) by the corresponding third-order Taylor expansion that would correspond to another calculated response, denoted as $r_{i_2}^c\left(\alpha \right)$, and integrating the resulting expression formally over the unknown parameter distribution $p_{\alpha }\left(\alpha \right)$ yields the following expression, which is used in HO-BERRU-PM for the covariance, denoted as $cov\left(r_{i_1}^c, r_{i_2}^c\right)$, of two responses, $r_{i_1}^c\left(\alpha \right)$ and $r_{i_2}^c\left(\alpha \right)$, $for i_1, i_2=1,\dots ,N_r$:

$$\begin{array}{l}\mathrm{cov}\left(r_{i_1}^c, r_{i_2}^c\right)={\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}\left(\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}\frac{\partial r_{i_2}^c}{\partial {\alpha }_j}\right){\rho }_{ij}{\sigma }_i{\sigma }_j}}+\frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}\left(\frac{{\partial }^2{r}_{i_1}^c}{\partial {\alpha }_i\partial {\alpha }_j}\frac{\partial r_{i_2}^c}{\partial {\alpha }_{\mu }}+\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}\frac{{\partial }^2{r}_{i_2}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }}\right)t_{ij\mu }{\sigma }_i{\sigma }_j{\sigma }_{\mu }}}}\\ +\frac{1}{4}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\left(\frac{{\partial }^2{r}_{i_1}^c}{\partial {\alpha }_i\partial {\alpha }_j}\right)}\left(\frac{{\partial }^2{r}_{i_2}^c}{\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}\right)\left(q_{ij\mu \nu }-{\rho }_{ij}{\rho }_{\mu \nu }\right){\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }}}}\\ +\frac{1}{6}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\left(\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}\frac{{\partial }^3{r}_{i_2}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}+\frac{\partial r_{i_2}^c}{\partial {\alpha }_i}\frac{{\partial }^3{r}_{i_1}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}\right)} q_{ij\mu \nu }{\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }}}} .\end{array}$$

(122)

Terms involving response sensitivities higher than the fourth order in the parameters’ standard deviations have been discarded in Equation (122). The expression in Equation (122) involves implicitly the fourth-order moment, ${\mu }_4^{ijkl}$, of the multivariate parameter distribution function $p_{\alpha }\left(\alpha \right)$ and, explicitly, the fourth-order parameter correlation, $q_{ijkl}$; these fourth-order moments are defined as follows:

$${\mu }_4^{ijkl}\left(\alpha \right)\triangleq {\displaystyle \underset{D_{\alpha }}{\int }\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_k-{\alpha }_k^0\right)\left({\alpha }_{\ell }-{\alpha }_{\ell }^0\right)p_{\alpha }\left(\alpha \right)d\alpha } \triangleq q_{ijkl}{\sigma }_i{\sigma }_j{\sigma }_k{\sigma }_{\ell }; i,j,k,\ell =1,\dots ,N_{\alpha }.$$

(123)

In particular, the variance of a response $r_{i_1}^c\left(\alpha \right)$ is obtained by setting $i_1=i_2$ in Equation (122).

When the parameters follow a multivariate Gaussian distribution (as it is often assumed in practical applications), the following relations hold:

$$t_{ij\mu }\equiv 0; q_{ij\mu \nu }={\rho }_{ij}{\rho }_{\mu \nu } +{\rho }_{i\mu }{\rho }_{j\nu }+{\rho }_{i\nu }{\rho }_{j\mu } ; i,j,\mu ,\nu =1,\dots ,N_{\alpha } .$$

(124)

The third-order cumulant for three responses, $r_{i_1}^c\left(\alpha \right)$, $r_{i_2}^c\left(\alpha \right)$ and $r_{i_3}^c\left(\alpha \right)$, which is denoted below as ${\mu }_3\left(r_{i_1}^c,r_{i_2}^c,r_{i_3}^c\right)$, $for i_1,i_2,i_3=1,\dots ,N_r$, is also obtained by using Equation (114) and has the following expression:

$$\begin{array}{l}{\mu }_3\left(r_{i_1}^c,r_{i_2}^c,r_{i_3}^c\right)\triangleq {\langle \left[r_{i_1}^c-E\left(r_{i_1}^c\right)\right]\left[r_{i_2}^c-E\left(r_{i_2}^c\right)\right]\left[r_{i_3}^c-E\left(r_{i_3}^c\right)\right]\rangle }_{\alpha }={\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}\frac{\partial r_{i_2}^c}{\partial {\alpha }_j}}}}\frac{\partial r_{i_3}^c}{\partial {\alpha }_{\mu }}{t}_{ij\mu }{\sigma }_i{\sigma }_j{\sigma }_{\mu }\\ +\frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\frac{\partial r_{i_1}^c{}_k}{\partial {\alpha }_i}\frac{\partial r_{i_2}^c}{\partial {\alpha }_j}}}}\frac{{\partial }^2{r}_{i_3}^c}{\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}\left(q_{ij\mu \nu }-{\rho }_{ij}{\rho }_{\mu \nu }\right){\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }}\\ +\frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}}}}\frac{{\partial }^2{r}_{i_2}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }}\frac{\partial r_{i_3}^c}{\partial {\alpha }_{\nu }}\left(q_{ij\mu \nu }-{\rho }_{i\nu }{\rho }_{j\mu }\right){\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }}\\ +\frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\frac{{\partial }^2{r}_{i_1}^c}{\partial {\alpha }_i\partial {\alpha }_j}\frac{\partial r_{i_2}^c}{\partial {\alpha }_{\mu }}\frac{\partial r_{i_3}^c}{\partial {\alpha }_{\nu }}}}}\left(q_{ij\mu \nu }-{\rho }_{ij}{\rho }_{\mu \nu }\right){\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }} .\end{array}$$

(125)

The covariance of a response, $r_{i_1}^c\left(\alpha \right)$ and a parameter ${\alpha }_{\ell }$, $i_1=1,\dots ,N_r$ and $\ell =1,\dots ,N_{\alpha }$, which is denoted as $\mathrm{cov}\left(r_{i_1}^c, {\alpha }_{\ell }\right)$, is obtained by multiplying Equation (114) by $\left({\alpha }_{\ell }-{\alpha }_{\ell }^0\right)$ and formally integrating the resulting expression over the unknown multivariate parameter distribution function $p_{\alpha }\left(\alpha \right)$ to obtain the following expression:

$$\begin{array}{l}\mathrm{cov}\left(r_{i_1}^c, {\alpha }_{\ell }\right)={{\displaystyle \sum _{i=1}^{N_{\alpha }}\left\{\frac{\partial r_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i}\right\}}}_{{\alpha }^0}\mathrm{cov}\left({\alpha }_i,{\alpha }_{\ell }\right)+\frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\left\{\frac{{\partial }^2{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\right\}}_{{\alpha }^0}{t}_{ij\ell }{\sigma }_i{\sigma }_j{\sigma }_{\ell }}}\\ +\frac{1}{6}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{k=1}^{N_{\alpha }}{\left\{\frac{{\partial }^3{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j\partial {\alpha }_k}\right\}}_{{\alpha }^0}{q}_{ij\mu \nu }{\sigma }_i{\sigma }_j{\sigma }_k{\sigma }_{\ell }}}}; i_1=1,\dots ,N_r; \ell =1,\dots ,N_{\alpha }.\end{array}$$

(126)

The expectation values $E\left[r_k^c\left(\alpha \right)\right]$, $k=1,\dots ,N_r$, given by Equation (116) are considered to be the components of the following vector of “expected values of the computed response”:

$$E\left[r^c\left(\alpha \right)\right]\triangleq {\left[E\left(r_1^c\right),\dots ,E\left(r_{N_r}^c\right)\right]}^{†}.$$

(127)

The response covariances defined in Equation (122) are considered to be the components of an $\left(N_r\times N_r\right)$-dimensional matrix denoted as $C_r$, which is defined as follows:

$$C_r\triangleq {\langle \left[r^c-E\left(r^c\right)\right]{\left[r^c-E\left(r^c\right)\right]}^{†}\rangle }_{\alpha }.$$

(128)

Using Equation (122) indicates that Equation (128) can be written as follows:

$$C_r={\tilde{C}}_{kk}+C_r^{HO}$$

(129)

where the matrix ${\tilde{C}}_{kk}$ is TSURFER’s “sandwich rule” provided in Equation (102) and where the elements $c_r^{HO}(i_1,i_2)$ of the matrix $C_r^{HO}\triangleq {\left[c_r^{HO}(i_1,i_2)\right]}_{N_r\times N_r}$,$i_1 ,i_2=1,..,N_r,$ comprise the contributions stemming from the higher-order response sensitivities and are defined as follows:

$$\begin{array}{l}{c}_r^{HO}(i_1,i_2)\triangleq \frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}\left(\frac{{\partial }^2{r}_{i_1}^c}{\partial {\alpha }_i\partial {\alpha }_j}\frac{\partial r_{i_2}^c}{\partial {\alpha }_{\mu }}+\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}\frac{{\partial }^2{r}_{i_2}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }}\right)t_{ij\mu }{\sigma }_i{\sigma }_j{\sigma }_{\mu }}}}\\ +\frac{1}{4}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\left(\frac{{\partial }^2{r}_{i_1}^c}{\partial {\alpha }_i\partial {\alpha }_j}\right)}\left(\frac{{\partial }^2{r}_{i_2}^c}{\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}\right)\left(q_{ij\mu \nu }-{\rho }_{ij}{\rho }_{\mu \nu }\right){\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }}}}\\ +\frac{1}{6}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{\mu =1}^{N_{\alpha }}{\displaystyle \sum _{\nu =1}^{N_{\alpha }}\left(\frac{\partial r_{i_1}^c}{\partial {\alpha }_i}\frac{{\partial }^3{r}_{i_2}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}+\frac{\partial r_{i_2}^c}{\partial {\alpha }_i}\frac{{\partial }^3{r}_{i_1}^c}{\partial {\alpha }_j\partial {\alpha }_{\mu }\partial {\alpha }_{\nu }}\right)} q_{ij\mu \nu }{\sigma }_i{\sigma }_j{\sigma }_{\mu }{\sigma }_{\nu }}}} .\end{array}$$

(130)

The covariances between the computed responses and the model parameters defined in Equation (126) are considered to be the components of an $\left(N_r\times N_{\alpha }\right)$-dimensional matrix denoted as $C_{r\alpha }$ and defined as follows:

$$C_{r\alpha }\triangleq {\langle \left[r^c-E\left(r^c\right)\right]{\left(\alpha -{\alpha }^0\right)}^{†}\rangle }_{\alpha }=C_{\alpha r}^{†}.$$

(131)

Using Equation (126) indicates that Equation (131) can be written as follows:

$$C_{\alpha r}=\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]+C_{\alpha r}^{HO}.$$

(132)

where the first-order matrix $\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]$ is as it appears in the posterior expressions calculated by TSURFER for the adjusted cross-sections and their adjusted covariances, see Equations (104) and (107), and where the elements $c_{r\alpha }^{HO}(i_1,\ell )$ of the matrix $C_{r\alpha }^{HO}\triangleq {\left[c_{r\alpha }^{HO}(i_1,\ell )\right]}_{N_r\times N_{\alpha }}={\left[C_{\alpha r}^{HO}\right]}^{†}$ comprise the contributions stemming from the higher-order response sensitivities and are defined as follows:

$$\begin{array}{l}{c}_{r\alpha }^{HO}(i_1,\ell )\triangleq \frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\left\{\frac{{\partial }^2{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\right\}}_{{\alpha }^0}{t}_{ij\ell }{\sigma }_i{\sigma }_j{\sigma }_{\ell }}}\\ +\frac{1}{6}{\displaystyle \sum _{i=1}^{N_{\alpha }}{\displaystyle \sum _{j=1}^{N_{\alpha }}{\displaystyle \sum _{k=1}^{N_{\alpha }}{\left\{\frac{{\partial }^3{r}_{i_1}^c\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j\partial {\alpha }_k}\right\}}_{{\alpha }^0}{q}_{ij\mu \nu }{\sigma }_i{\sigma }_j{\sigma }_k{\sigma }_{\ell }}}}; i_1=1,\dots ,N_r; \ell =1,\dots ,N_{\alpha }.\end{array}$$

(133)

The prior information used in the HO-BERRU-PM regarding the measured responses is the same as in TSURFER, namely: (i) the expectation (or mean value) of the experimentally measured responses $r_i^m$, which is denoted as $E\left(r_i^m\right)$; (ii) the covariance between two measured responses, which is denoted as $\mathrm{cov}\left(r_i^m,r_j^m\right)$, $i,j=1,\dots ,N_r$.

In the HO-BERRU-PM methodology, the expected values of the measured responses are considered to constitute the components of a column vector $r^m$ defined as follows:

$$E\left(r^m\right)\triangleq {\left[E\left(r_1^m\right),\dots ,E\left(r_{N_r}^m\right)\right]}^{†},$$

(134)

while the covariance matrix of measured responses is denoted as $C_m$ and is defined as follows:

$$C_m \triangleq {\left[\mathrm{cov}\left(r_i^m,r_j^m\right)\right]}_{N_r\times N_r}$$

(135)

Cacuci’s HO-BERRU-PM [13,21] methodology combines experimental and computational information in the joint phase space of the measured responses, calculated responses and model parameters, including not only the first-order response sensitivities, but also the complete set of second- and third-order sensitivities of the model responses to the model parameters. In contradistinction to the GLLSA procedure in TSURFER, the HO-BERRU-PM formalism uses the maximum entropy principle to eliminate the need for introducing and “minimizing” the user-chosen “cost functional quantifying the discrepancies between measurements and computations”, as done in TSURFER. The use of the maximum entropy principle enables the HO-BERRU-PM to consider the model as an “imprecisely known” entity, as opposed to being appended by Lagrange multipliers as a “hard constraint” devoid of any model errors, as the TSURFER’s GLLSA procedure does when appending the first-order approximation to the model represented by Equation (96) as a hard constraint to the augmented Lagrangian to be minimized. Consequently, the HO-BERRU-PM methodology yields best-estimate results that not only incorporate second- and third-order response sensitivities but are also free of subjective user interferences and can include numerical errors stemming from the numerical solution of the model equations.

After the prior information is assimilated/consolidated within the HO-BERRU-PM methodology, the saddle-point method is used to evaluate the joint posterior distribution obtained in the combined phase spaces of the model parameters, calculated responses and measured responses, to obtain the following best-estimate results:

• The best-estimate posterior expectation values for the vector of predicted model parameters; this vector is denoted as ${\alpha }^{be}$ and is given by the following formula:

  $${\alpha }^{be}={\alpha }^0-C_{\alpha r}{\left(C_m+C_r\right)}^{-1}d,$$

  (136)

  where the vector $d$ denotes the vector of deviations between the expected value of the computed response and the expected value of the measured response, namely:

  $$d\triangleq \left[E\left(r^c\right)-E\left(r^m\right)\right].$$

  (137)
• The best-estimate posterior parameter covariance matrix, denoted as $C_{\alpha }^{be}$, for the best-estimate parameters ${\alpha }^{be}$:

  $$C_{\alpha }^{be}=C_{\alpha }-C_{\alpha r}{\left(C_m+C_r\right)}^{-1}{C}_{r\alpha }.$$

  (138)

In Equation (138), the matrices $C_{\alpha }$ and $C_{\alpha r}{\left(C_m+C_r\right)}^{-1}{C}_{r\alpha }$ are symmetric and positive definite. Therefore, the subtraction indicated in Equation (138) implies that the components of the main diagonal of $C_{\alpha }^{be}$ must have smaller values than the corresponding elements of the main diagonal of $C_{\alpha }$. In this sense, the introduction of new computational and experimental information has produced reduced (by comparison to the a priori variances of the model parameters) values for the posterior best-estimate parameter variances (which are the elements on the diagonal of $C_{\alpha }^{be}$).

3.

The best-estimate posterior expectation values for the vector of predicted responses, which is denoted as $r^{be}$ and which has the following expression:

$$r^{be}=r^m+C_m{\left(C_m+C_r\right)}^{-1}\left[E\left(r^c\right)-E\left(r^m\right)\right],$$

(139)

4.

The best-estimate posterior covariance matrix for the best-estimate responses $r^{be}$, which is denoted as $C_r^{be}$ and which has the following expression:

$$C_r^{be}=C_m-C_m{\left(C_m+C_r\right)}^{-1}{C}_m .$$

(140)

As indicated in Equation (140), the initial covariance matrix $C_m$ is multiplied by the matrix $\left[I-{\left(C_m+C_r\right)}^{-1}{C}_m\right]$, which means that the variances contained on the diagonal of the best-estimate matrix $C_r^{be}$ will be smaller than the experimentally measured variances contained in $C_m$. Hence, the addition of new experimental information will reduce the predicted best-estimate response variances in $C_r^{be}$ by comparison to the measured variances contained a priori in $C_m$.

5.

The posterior covariance matrix comprising the best-estimate correlations between the best-estimate parameters ${\alpha }^{be}$ and the best-estimate responses $r^{be}$, which is denoted as $C_{\alpha r}^{be}$, has the following expression:

$$C_{\alpha r}^{be}=C_{\alpha r}{\left(C_m+C_r\right)}^{-1}{C}_m .$$

(141)

6.

The a “chi-square” indicator, which measures the agreement/disagreement between the experimental responses and the computed responses. This indicator is denoted as ${\chi }_{HO}^2$ (since it includes higher-order sensitivities and correlations) and has the following expression:

$${\chi }_{HO}^2\triangleq \frac{1}{2}{\left[E\left(r^c\right)-E\left(r^m\right)\right]}^{†}{\left(C_m+C_r\right)}^{-1}\left[E\left(r^c\right)-E\left(r^m\right)\right].$$

(142)

As the expression obtained in Equation (142) indicates, the “consistency indicator” ${\chi }^2$ can be evaluated directly from the a priori data, after having inverted the matrix $\left(C_m+C_r\right)$. The indicator ${\chi }^2$ can be interpreted as measuring (in the corresponding metric) the square of the length of the vector $d\triangleq \left[E\left(r^c\right)-E\left(r^m\right)\right]$, thus indicating the degree of consistency among the a priori information available about the model parameters, calculated and measured responses. As the dimension of $\left[E\left(r^c\right)-E\left(r^m\right)\right]$ indicates, the number of degrees of freedom characteristic of the calibration under consideration is equal to the number $N_r$ of experimental responses. In the extreme case of “absence of experimental responses”, no actual calibration occurs. An actual calibration (adjustment) occurs only when including at least one experimental response.

4.3. Comparison: HO-BERRU-PM Methodology versus TSURFER-GLLSA Methodology

The correspondences between the notations used in TSURFER and HO-BERRU-PM, respectively, are as follows:

(i): The number of parameters is denoted as $N_{\alpha }$ in HO-BERRU-PM; the number of parameters is denoted as $M$ in TSURFER; hence, $N_{\alpha }\equiv M$.

(ii) The number of computed responses is denoted as $N_r$ in HO-BERRU-PM; the number of parameters is denoted as $I$ in TSURFER; hence, $N_r\equiv I$.

(iii) The expected parameter values are denoted as ${\alpha }^0\triangleq \left({\alpha }_1^0,\dots ,{\alpha }_{N_{\alpha }}^0\right)$, using the superscript “zero” in HO-BERRU-PM; the expected parameter values are denoted as $E\left[\alpha \right]\triangleq {\left[E\left({\alpha }_1\right),\dots ,E\left({\alpha }_M\right)\right]}^{†}$ in TSURFER; the correspondence is evident.

It is noteworthy that the equivalent of the vector $d\triangleq \left[E\left(r^c\right)-E\left(r^m\right)\right]$ also appears in the posterior quantities produced by the TSURFER-GLLSA procedure, in the form of the vector $\tilde{d}\triangleq k\left[E\left(\alpha \right)\right]-m$, as defined in Equation (106), but $d\ne \tilde{d}$ since $E\left(r^c\right)\ne k\left[E\left(\alpha \right)\right]$, as shown in Equation (119). Recalling the expression of $E\left(r^c\right)$ from Equation (119), it follows that the following relation holds:

$$d\triangleq \left[E\left(r^c\right)-E\left(r^m\right)\right]=\tilde{d}+E^{HO}\left[r_i^c\left({\alpha }^0\right)\right] .$$

(143)

4.3.1. A Priori Information: TSURFER versus HO-BERRU-PM

HO-BERRU-PM methodology considers all of the a priori information incorporated into the TSURFER-GLLSA methodology, as well as the following additional a priori information: (i) second-order and third-order sensitivities of the calculated responses with respect to the model parameters; (ii) third-order correlations between parameters; (iii) higher-order expressions for the expected values of computed responses; (iv) higher-order expressions for the covariances between computed responses; (v) higher-order expressions for the covariances between computed responses and model parameters.

Table 8. A Priori Information: TSURFER-GLLSA vs. HO-BERRU-PM.

A Priori Quantity	TSURFER	HO-BERRU-PM	[HO-BERRU-PM]−[TSURFER]
Second-order sensitivities	N/A	Yes	Equation (115)
Third-order sensitivities	N/A	Yes	Equation (115)
Third-order parameter correlations	N/A	Yes	Equation (117)
Expected value of calculated response	Equation (97)	Equation (119)	Equation (121)
Covariance of two calculated responses	Equation (102)	Equation (129)	Equation (130)
Parameter-calculated response covariance	$\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]$	Equation (132)	Equation (133)
“Vector of deviations”	$\tilde{d}$ [Equation (106) ]	$d$ [Equation (137)]	Equation (143)

summarizes the distinctions between TSURFER-GLLSA and HO-BERRU-PM methodologies regarding the incorporation of a priori information.

4.3.2. Posterior Results: TSURFER Adjustments versus HO-BERRU-PM Predictions

At the most fundamental conceptual level, the TSURFER-GLLSA methodology relies fundamentally on minimizing a user-chosen “chi-square” functional subject to the requirement that the following hard constraint be satisfied: “the first-order Taylor-expansion of the calculated response as a linear function of the model parameters”. In contradistinction, the HO-BERRU-PM methodology eliminates the need for a “user-chosen functional to be minimized” by employing the maximum entropy principle to combine, without imposing any “hard constraints”, the available information regarding the measured and computed responses, including a nonlinear third-order Taylor-series representation of the underlying model, connecting the model parameters to the model’s calculated responses. The distinctions between the best-estimate parameter values, ${\alpha }^{be}$, predicted by the HO-BERRU-PM methodology, and the adjusted parameter values, ${\alpha }^{\prime }$, produced by TSURFER-GLLSA, are obtained by replacing Equations (129), (132) and (143) in Equation (136) and subtracting Equation (104) from the resulting equation. These operations yield the following relation:

$$\begin{array}{l}{\alpha }^{be}-{\alpha }^{\prime }=-C_{\alpha r}{\left(C_m+C_r\right)}^{-1}d+\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]{\tilde{C}}_{dd}^{-1}\tilde{d}\\ =-\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}+C_{\alpha r}^{HO}\right]{\left(C_m+{\tilde{C}}_{kk}+C_r^{HO}\right)}^{-1}\left\{\tilde{d}+E^{HO}\left[r_i^c\left({\alpha }^0\right)\right]\right\} +\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]{\left(C_m+{\tilde{C}}_{kk}\right)}^{-1}\tilde{d}.\end{array}$$

(144)

The difference between the best-estimate covariances $C_{\alpha }^{be}$ for the best-estimate parameter values, ${\alpha }^{be}$, predicted by the HO-BERRU-PM methodology and the adjusted covariances $C_{{\alpha }^{\prime }{\alpha }^{\prime }}$ for the adjusted parameter values, ${\alpha }^{\prime }$, produced by TSURFER, is obtained by replacing Equations (129) and (132) in Equation (138) and subtracting Equation (107) from the resulting equation. These operations yield the following relation:

$$\begin{array}{l}{C}_{\alpha }^{be}-C_{{\alpha }^{\prime }{\alpha }^{\prime }}=-C_{\alpha r}{\left(C_m+C_r\right)}^{-1}{C}_{r\alpha }+\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]{\tilde{C}}_{dd}^{-1}\left[{\tilde{S}}_{k\alpha }{C}_{\alpha \alpha }\right]\\ =-\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}+C_{\alpha r}^{HO}\right]{\left(C_m+{\tilde{C}}_{kk}+C_r^{HO}\right)}^{-1}{\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}+C_{\alpha r}^{HO}\right]}^{†}+\left[C_{\alpha \alpha }{\tilde{S}}_{ka}^{†}\right]{\left(C_m+{\tilde{C}}_{kk}\right)}^{-1}\left[{\tilde{S}}_{k\alpha }{C}_{\alpha \alpha }\right].\end{array}$$

(145)

The difference between the best-estimate response values, $r^{be}$, predicted by the HO-BERRU-PM methodology and the adjusted response values, $m^{\prime }$, produced by TSURFER, is obtained by replacing Equations (129) and (143) in Equation (139) and subtracting Equation (108) from the resulting equation. These operations yield the following relation:

$$\begin{array}{l}{r}^{be}-m^{\prime }=C_m{\left(C_m+C_r\right)}^{-1}\left[E\left(r^c\right)-E\left(r^m\right)\right]-{\tilde{C}}_{mm}{\tilde{C}}_{dd}^{-1}\tilde{d}\\ =C_m{\left(C_m+{\tilde{C}}_{kk}+C_r^{HO}\right)}^{-1}\left\{\tilde{d}+E^{HO}\left[r_i^c\left({\alpha }^0\right)\right]\right\}-C_m{\left(C_m+{\tilde{C}}_{kk}\right)}^{-1}\tilde{d}.\end{array}$$

(146)

The difference between the best-estimate covariances $C_r^{be}$ for the best-estimate response values, $r^{be}$, predicted by the HO-BERRU-PM methodology and the adjusted covariances, ${\tilde{C}}_{m^{\prime }{m}^{\prime }}$, for the adjusted responses $m^{\prime }$ produced by TSURFER-GLLSA is obtained by replacing Equation (129) in Equation (140) and subtracting Equation (109) from the resulting equation. These operations yield the following relation:

$$\begin{array}{l}{C}_r^{be}-{\tilde{C}}_{m^{\prime }{m}^{\prime }}=-C_m{\left(C_m+C_r\right)}^{-1}{C}_m +{\tilde{C}}_{mm}{\tilde{C}}_{dd}^{-1}{\tilde{C}}_{mm}\\ =C_m{\left(C_m+{\tilde{C}}_{kk}+C_r^{HO}\right)}^{-1}{C}_m-C_m{\left(C_m+{\tilde{C}}_{kk}\right)}^{-1}{C}_m.\end{array}$$

(147)

The difference between the high-order indicator, ${\chi }_{HO}^2$, produced by the HO-BERRU-PM methodology and the indicator ${\chi }_{TSURFER}^2$ produced by TSURFER is obtained by replacing Equations (129) and (143) in Equation (142) and subtracting Equation (112) from the resulting equation. These operations yield the following relation:

$$\begin{array}{l}{\chi }_{HO}^2-{\chi }_{TSURFER}^2=\frac{1}{2}{\left[E\left(r^c\right)-E\left(r^m\right)\right]}^{†}{\left(C_m+C_r\right)}^{-1}\left[E\left(r^c\right)-E\left(r^m\right)\right]-d^{†}{C}_{dd}^{-1}d\\ ={\left\{\tilde{d}+E^{HO}\left[r_i^c\left({\alpha }^0\right)\right]\right\}}^{†}{\left(C_m+{\tilde{C}}_{kk}+C_r^{HO}\right)}^{-1}\left\{\tilde{d}+E^{HO}\left[r_i^c\left({\alpha }^0\right)\right]\right\} -{\tilde{d}}^{†}{\left(C_m+{\tilde{C}}_{kk}\right)}^{-1}\tilde{d}.\end{array}$$

(148)

The distinctions between the posterior adjusted quantities produced by TSURFER-GLLSA versus the HO-BERRU-PM framework are summarized in

Table 9. Posterior information: TSURFER-GLLSA vs. HO-BERRU-PM.

Posterior Quantity	HO-BERRU-PM Predicted Quantity	TSURFER Adjusted Quantity	[HO-BERRU-PM]−[TSURFER] Differences
Expected parameter value	Equation (136)	Equation (104)	Equation (144)
Parameter covariance	Equation (138)	Equation (107)	Equation (145)
Expected response value	Equation (139)	Equation (108)	Equation (146)
Response covariance	Equation (140)	Equation (109)	Equation (147)
Parameter-response covariance	Equation (141)	N/A	N/A
Chi-square indicator	Equation (142)	Equation (112)	Equation (148)

.

It is evident from Equations (144)−(148) that the expressions of the quantities predicted by the HO-BERRU-PM methodology comprise, as particular cases, the quantities adjusted by TSURFER-GLLSA. The quantities predicted by the HO-BERRU-PM methodology will become identical to the quantities adjusted by TSURFER-GLLSA if (and only if) all of the sensitivities higher than the first order are ignored. Notably, the computation within the HO-BERRU-PM methodology of the best-estimate parameter and response values, together with their corresponding best-estimate covariance matrices, only requires the computation of ${\left(C_m+C_r\right)}^{-1}$, which entails the inversion of a matrix of size $N_r\times N_r$. The matrix $\left(C_m+C_r\right)$ has the same dimensions (of the number of measured responses) as the corresponding matrix, ${\tilde{C}}_{dd}$, which needs to be inverted in TSURFER, even though the matrices in the HO-BERRU-PM methodology comprise not only the first-order response sensitivities (as in TSURFER-GLLSA) but also comprise all of the second- and third-order response sensitivities to the model parameters. This important HO-BERRU-PM characteristic is essential for applications to large-scale practical problems, since in the overwhelming majority of practical situations, the number of model parameters far exceeds the number of responses of interest.

5. Discussion and Conclusions

This work has underscored the need for computing higher-order (i.e., higher than first-order) sensitivities (functional derivatives) of model responses with respect to the model parameters. The significant potential impact of the higher-order sensitivities on the expected values and variances/covariances for the calculated and predicted model responses have also been highlighted. In particular, if only first-order sensitivities are considered, the third-order moment of the response is always zero. Hence, a “first-order sensitivity and uncertainty quantification” will always produce an erroneous third moment (and, hence, skewness) of the predicted response distribution, unless the unknown response distribution happens to be symmetrical. At least second-order sensitivities must be used in order to estimate the third-order moment (and, hence, the skewness) of the response distribution. Skewness indicates the direction and relative magnitude of a distribution’s deviation from the normal distribution, while kurtosis indicates the propensity of the predicted response distribution to have heavy tails and/or outliers. With pronounced skewness, standard statistical inference procedures such as constructing a confidence interval for the mean (expectation) of a computed/predicted model response will be not only incorrect, in the sense that the true coverage level will differ from the nominal (e.g., 95%) level, but the error probabilities will be unequal on each side of the predicted mean.

The evident conclusion that can be drawn from this work is that the consideration of only the first-order sensitivities is insufficient for making credible predictions regarding the expected values and uncertainties (variances, covariances, skewness) of calculated and predicted/adjusted responses. At the very least, the second-order sensitivities must also be computed in order to enable the quantitative assessment of their impact on the predicted quantities. Since the second-order sensitivities impact decisively the expected values and the skewness of the calculated/predicted responses, they will also impact the computation of confidence intervals and tolerance limits for the predicted expectation of these responses. For large-scale nonlinear systems, the second-order response sensitivities can be computed exactly and efficiently by applying the adjoint methodology developed by Cacuci [5]. Extending the methodology presented in [5] to enable the exact and efficient computation of third-order (and higher-order) sensitivities for large-scale nonlinear systems is currently in progress. For large-scale (response-coupled) linear forward/adjoint systems, all of the response sensitivities up to and including the fourth-order sensitivities can be computed exactly and most efficiently by applying the methodology (fourth CASAM), recently developed by Cacuci [30].