The accompanying work (“part I”) has presented the newly developed mathematical framework, the “n^th-order feature adjoint sensitivity analysis methodology for response-coupled forward/adjoint linear systems” (abbreviated as “n^th-FASAM-L”), conceived by Cacuci (2024c). This work illustrates the application of the n^th-FASAM-L to a representative energy-dependent neutron-slowing down model of fundamental importance to reactor physics. The physical considerations underlying this model are presented in Section 2, which briefly reviews the concept of “contributon-flux density response” and particularizes this concept within the modeling of neutron slowing down in a mixture of materials. This physical model is of fundamental importance in nuclear reactor physics and enables the derivation of exact closed-form results for the application of the n^th-FASAM-L. Section 2 also defines the “features” inherent to this model, which enable the advantageous application of the n^th-FASAM-L. By definition, there are considerably fewer “feature functions” of the primary model parameters than there are primary model parameters.

Section 3 presents the first-order adjoint sensitivity analysis of the contributon flux with respect to the features and primary model parameters of the slowing-down model, comparing the application of the 1^st-FASAM-L versus the first-order comprehensive adjoint sensitivity analysis methodology for response-coupled forward/adjoint linear systems (1^st-CASAM-L). Using either the 1^st-FASAM-L or 1^st-CASAM-L involves solving the same operator equations and boundary conditions within the respective 1^st-LASS but with differing source terms. For the computation of the first-order sensitivities, the 1^st-FASAM-L enjoys only a slight computational advantage since it requires only one quadrature per component of the feature function, whereas the 1^st-CASAM-L requires one quadrature per primary model parameter.

Section 4 presents the second-order adjoint sensitivity analysis of the contributon flux with respect to the features and primary model parameters of the slowing-down model, comparing the application of the 2^nd-FASAM-L versus the 2^nd-CASAM-L. It is shown that the 2^nd-FASAM-L requires as many large-scale “adjoint” computations as there are non-vanishing first-order response sensitivities with respect to the components of the feature functions, whereas the 2^nd-CASAM-L requires as many large-scale computations as there are non-vanishing first-order response sensitivities with respect to the primary model parameters. Hence, the 2^nd-FASAM-L is inherently more efficient than the 2^nd-CASAM-L. In particular, one of the three distinct second-order sensitivities with respect to the model’s features vanishes identically within the 2^nd-FASAM-L but none of the ca. 100 second-order sensitivities with respect to the primary model parameters vanish within the 2^nd-CASAM-L.

Section 5 presents the third-order adjoint sensitivity analysis of the contributon flux with respect to the features and primary model parameters of the slowing-down model, comparing the application of the 3^rd-FASAM-L versus the 3^rd-CASAM-L. For computing the exact expressions of the third-order contributon-response sensitivities, the 3^rd-FASAM-L requires only two large-scale computations, whereas the 3^rd-CASAM-L would require hundreds of large-scale computations.

The concluding discussion presented in Section 6 emphasizes the fact that the unparalleled efficiency of the n^th-FASAM-N increases as the order of computed sensitivities increases, and the probability of encountering vanishing sensitivities is much higher when using the n^th-FASAM-L rather than any other methodology. Both the n^th-FASAM-L and n^th-CASAM-L overcome the limitation of dimensionality in the sensitivity analysis of linear systems, being incomparably more efficient and more accurate than any other method (statistical, finite differences, etc.) for computing exact expressions of response sensitivities (of any order) with respect to the uncertain parameters, boundaries, and internal interfaces of the model.

Fundamentally important responses of linear models depend simultaneously on both the forward and adjoint state functions governing the respective linear model, which makes it necessary to treat linear models/systems in their own right since such responses cannot be treated as particular cases of responses of nonlinear models. Typical examples of such responses arise in the modeling of self-diffusion processes in which the interaction mean free path is independent of the phase-space density. Such processes are modeled by linear equations of the Lorentz–Boltzmann type, and they occur in neutron, electron, and photon transport through media, as well as in certain types of transport processes in gas or plasma dynamics. Numerically solving such time-dependent integro-differential equations, albeit linear, is representative of “large-scale” computations and will be used in the sequel for illustrating the application of the n^th-FASAM-L. In particular, the distribution of neutrons in a medium is modeled by the following standard form of the linear Boltzmann equation: L r , E , Ω , t φ r , E , Ω , t = Q r , E , Ω , t , (1) where the linear integro-differential operator L r , E , Ω , t is defined below: L r , E , Ω , t φ r , E , Ω , t ≜ 1 v ∂ φ r , E , Ω , t ∂ t + Ω · ∇ φ r , E , Ω , t + Σ t r , E   φ r , E , Ω , t − ∫ 0 E f d E ′ ∫ 4 π d Ω ′ Σ s r , E ′ → E , Ω ′ → Ω   φ r , E ′ , Ω ′ , t − ∫ 0 E f d E ′ ∫ 4 π χ r , E ′ → E ν Σ r , E ′   φ r , E ′ , Ω ′ , t d Ω ′ . (2)

The quantities that appear in the standard notation used in Equation 2 are defined as follows:

(i) r denotes the three-dimensional position vector in space; E denotes the energy-independent variable; the directional vector Ω denotes the scattering solid angle; t denotes the time-independent variable; and v denotes the neutron particle speed.

The adjoint Boltzmann transport equation is formulated in the Hilbert space denoted as H B and is endowed with the following inner product, denoted as φ r , E , Ω , t , ψ r , E , Ω , t B , between two elements φ r , E , Ω , t ∈ H B and ψ r , E , Ω , t ∈ H B : φ , ψ B ≜ ∫ 0 t f d t ∫ 0 ∞ d E ∫ 4 π d Ω ∫ V d V φ r , E , Ω , t ψ r , E , Ω , t   . (3)

In the Hilbert space H B , the generic adjoint Boltzmann transport equation is as follows: L * r , E , Ω , t ψ r , E , Ω , t = Q * r , E , Ω , t , (4) where the (adjoint) linear integro-differential operator L * r , E , Ω , t is defined below: L * r , E , Ω , t ψ r , E , Ω , t ≜ − 1 v ∂ ψ r , E , Ω , t ∂ t − Ω · ∇ ψ r , E , Ω , t + Σ t r , E   ψ r , E , Ω , t − ∫ 0 E f d E ′ ∫ 4 π d Ω ′ Σ s r , E → E ′ , Ω → Ω ′   ψ r , E ′ , Ω ′ , t − ν Σ r , E ∫ 0 E f d E ′ ∫ 4 π χ r , E → E ′   ψ r , E ′ , Ω ′ , t d Ω ′ . (5)

By construction, the forward and adjoint transport equations satisfy the following relation: φ , L * ψ B − ψ , L φ B = P φ , ψ = φ , Q * B − ψ , Q B , (6) where P φ , ψ denotes the bilinear concomitant evaluated on the boundary of the phase-space domain under consideration. The “generalized reciprocity relation” expressed by Equation 6 relates the bilinear concomitant, which is a functional of the forward and adjoint fluxes at the initial and final times along the incoming and outgoing directions at the surface of the medium, to the fluxes in the interior of the medium comprising fixed sources. This reciprocity relation provides a physical interpretation of the adjoint flux as an “importance function,” which quantifies the contribution of a source to a detector and enables transport problems to be posed either in the forward or adjoint descriptions. These reciprocity relations also restrict the combination of forward and adjoint boundary conditions to those that ensure both the forward and adjoint formulations are mathematically “well posed.” The reciprocity relation expressed by Equation 6 is extensively used in the so-called “source-detector” problems in steady-state subcritical systems, where Q * r , E , Ω models the detector properties (cross section) in the sub-region occupied by the respective detector.

When the boundary conditions for Equation 1 are homogeneous and there is no external source, i.e., when Q r , E , Ω , t = 0 , the stationary neutron transport problem becomes an eigenvalue problem. The largest (i.e., fundamental) eigenvalue in such a case is called the “effective multiplication factor” and, depending on its value, corresponds to a critical, subcritical, or supercritical physical system (e.g., nuclear reactor). This eigenvalue (multiplication factor) is an important system (model) response, and its mathematical expression is a functional (“Raleigh quotient”) of the forward and the adjoint fluxes. Additional important model responses that are functionals of both the forward and adjoint fluxes include the reactivity, generation time, and lifetime of the system, along with several other Lagrangian functionals used in variational principles for developing efficient Raleigh–Ritz type numerical methods (see, e.g., Lewins, 1965; Stacey, 1974; Stacey, 2001). Perhaps the simplest quantity that depends on both the forward and adjoint fluxes—and has important applications in particle transport (particularly in particle shielding)—is the so-called “contributon flux” (Williams and Engle, 1977), which arises as follows:

(i) Multiplying the stationary form of Equation 1 by ψ r , E , Ω , multiplying the stationary form of Equation 4 by φ r , E , Ω , subtracting the resulting equations from each other, and integrating the resulting equation over only the energy- and solid angle-independent variables yield the following relation:

The application of the n^th-FASAM-L is illustrated in this section by considering the simplified model of the distribution in the asymptotic energy range of neutrons produced by a source of neutrons placed in an isotropic medium comprising a homogeneous mixture of “M” non-fissionable materials having constant (i.e., energy-independent) properties. For simplicity, but without diminishing the applicability of the n^th-FASAM-L, this medium is considered to be infinitely large. The simplified form of the Boltzmann neutron transport equation, as shown in Equation 1, that models the energy distributions of neutrons within a mixture of materials is called the “neutron slowing-down equation.” This equation is written using neutron lethargy (rather than the neutron energy) as the independent variable. Neutron lethargy is customarily denoted using the variable/letter “u” and is defined as u ≜ ⁡ ln E 0 / E , where E denotes the energy variable and E 0 denotes the highest energy in the system. Thus, the neutron slowing-down model (see, e.g., Meghreblian and Holmes, 1960; Lamarsh, 1966) for the energy distribution of the neutron flux in a homogeneous mixture of non-fissionable materials of infinite extent takes the following simplified form of Equation 1: d φ u d u + Σ a   ξ ¯   Σ t φ u = S u ξ ¯   Σ t ;   0 < u ≤ u t h ; (12) φ 0 = 0 ;   a t   u = 0 . (13)

The quantities that appear in Equation 12 are defined as follows.

(i) The lethargy-dependent neutron flux is denoted as φ u ; u t h denotes a cut-off lethargy, usually corresponding to the thermal neutron energy (ca. 0.0024 electron volts).

Mirroring the considerations for the Boltzmann transport equation presented in Equations 1‒6, the “adjoint slowing-down model” is constructed in the Hilbert space H B of square-integrable functions φ u ∈ H B and ψ u ∈ H B endowed with the following inner product, denoted as φ u , ψ u B : φ u , ψ u B ≜ ∫ 0 u t h φ u ψ u d u . (19)

Using the inner product φ u , ψ u B defined in Equation 19, the adjoint slowing-down model is constructed by the usual procedure, i.e., by (i) constructing the inner product of Equation 12 with a function ψ u ∈ H B ; (ii) integrating by parts the resulting relation so as to transfer the differential operation from the forward function φ u onto the adjoint function ψ u ; (iii) using the initial condition provided in Equation 13 and eliminating the unknown function φ u t h by choosing the final-value condition ψ u t h = 0 ; and (iv) choosing the source for the resulting adjoint slowing-down model so as to satisfy the generalized reciprocity relation shown in Equation 6. The result of these operations is the following adjoint slowing-down model for the adjoint slowing-down function ψ u : − d ψ u d u + f 1 α ψ u = δ u − u d , (20) ψ u t h = 0 , a t   u = u t h . (21)

The “contributon-flux response density” R c φ , ψ , as generally defined in Equation 8, specialized for the neutron slowing-down model, coincides with the inner product used in this context, i.e., R c φ , ψ ≜ ∫ 0 u t h φ u ψ u   d u ≡ φ u , ψ u B . (22)

It is important to note that R c φ , ψ does not depend explicitly on either the feature function f α or any primary model parameter. Therefore, the G-differential of R c φ , ψ will not comprise a direct-effect term but will consist entirely of the indirect-effect term.

For this “contributon-flux response density” model, the following primary model parameters are subject to experimental uncertainties.

(i) For each material “i,” i = 1 , . . . , M , included in the homogeneous mixture, the following are primary model parameters: the atomic number densities N m i ; the microscopic radiative-capture cross section σ γ i ; and the scattering cross section σ s i ;

The above primary parameters are considered to constitute the components of a “vector of primary model parameters” defined as follows: α ≜ N m 1 , σ γ 1 , σ s 1 , . . . , N m M , σ γ M , σ s M , λ 1 S , λ 2 S , N 1 S , N 2 S , F 1 S , F 2 S , ν 1 S , ν 2 S , W 1 S , W 2 S † ≜ α 1 , … , α T P † ;   T P ≜ 3 M + 10 . (23)

The first-level forward/adjoint system (1^st-LFAS) for the “first-level forward/adjoint function” u 1 2 ; u ≜ φ u , ψ u † comprises Equations 12, 13, 20, and 21. The structure of the 1^st-LFAS suggests that the components f i α of the feature function f α can be defined as follows: f α ≜ f 1 α , f 2 α † ;   f 1 α ≜ Σ a α   ξ ¯   α Σ t α ;   f 2 α ≜ S 0 α ξ ¯   α Σ t α   . (24)

Solving Equations 12, 13 while using the definitions introduced in Equation 24 yields the following expression for the flux φ u in terms of the components f i α of the feature function f α : φ u = H u f 2 α exp − u f 1 α ;   H 0 = 0 ;   H u = 1 , i f   u > 0 . (25)

Solving the above adjoint slowing-down model yields the following closed-form expression for the adjoint slowing-down function ψ u : ψ u = H u d − u exp u − u d f 1 α   . (26)

In terms of the components f i α of the feature function f α , the closed-form expression of the “contributon response density” is obtained by substituting the expressions provided in Equations 25, 26 into Equation 22 and performing the integration over lethargy, which yields R c φ , ψ = ∫ 0 u t h H u f 2 α exp − u f 1 α H u d − u exp u − u d f 1 α   d u = u d f 2 α exp − u d f 1 α   . (27)

In terms of the primary model parameters, the closed-form expression of the “contributon response density” is R c φ , ψ = u d S 0 α ξ ¯   α Σ t α exp − u d Σ a α   ξ ¯   α Σ t α   . (28)

As Equation 28 indicates, the model response can be considered to depend directly on T P ≜ 3 M + 10 primary model parameters. In view of Equation 27, however, the model response can alternatively be considered to depend directly on two feature functions and only indirectly (through the two feature functions) on the primary model parameters. In the former consideration/interpretation, the response sensitivities to the primary model parameters will be obtained by applying the n^th-CASAM-L. In the later consideration/interpretation, the response sensitivities to the primary model parameters will be obtained by applying the n^th-FASAM-L, which will involve two stages: (a) the response sensitivities with respect to the feature functions will be obtained in the first stage; (b) the subsequent computation of the response sensitivities to the primary model parameters will be performed in the second stage by using the response sensitivities with respect to the feature functions obtained in the first stage. The computational distinctions that stem from these differing considerations/interpretations underlying the n^th-CASAM-L versus the n^th-FASAM-L will become evident in the next section by using a paradigm neutron slowing-down model, which is representative of the general situation for any linear system.

The first-order sensitivities of the response R c u 1 2 ; u , where u 1 2 ; u ≜ φ u , ψ u † , are obtained by determining the first-order Gateaux (G-)- differential, denoted as δ R c u 1 2 ; u , v 1 2 ; u α 0 , of this response for variations v 1 2 ; u ≜ δ φ u , δ ψ u † around the phase-space point φ 0 , ψ 0 . By definition, the first-order G-differential δ R c u 1 2 ; u , v 1 2 ; u α 0 is obtained as follows: δ R c u 1 2 ; u , v 1 2 ; u α 0 ≜ d d ε ∫ 0 u t h φ 0 u + ε v 1 u   ψ 0 u + ε δ ψ u d u ε = 0 = ∫ 0 u t h v 1 u ψ u + φ u δ ψ u   d u α 0 . (29)

The sensitivities of R c u 1 2 ; u with respect to the feature functions (and subsequently to the primary model parameters) will be determined in Section 3.1 by applying the 1^st-FASAM-L. Alternatively, the sensitivities of R c u 1 2 ; u directly with respect to the primary model parameters will be determined in Section 3.2 by applying the 1^st-CASAM-L.

In practice, closed-form expressions such as those shown in Equations 53, 54 are unavailable. The 1^st-FASAM-L yields the expressions provided in Equations 48, 49, while the 1^st-CASAM-L yields the expressions provided in Equation 57. Hence, these expressions will provide the starting points for obtaining the second-order sensitivities that stem from the respective first-order sensitivities. As outlined within the general frameworks of both the n^th-FASAM-L and n^th-CASAM-L methodologies, the second-order sensitivities are obtained by conceptually considering them to arise as the “first-order sensitivities of the first-order sensitivities.”