The Second-Order Features Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (2nd-FASAM-L): Mathematical Framework and Illustrative Application to an Energy System

1. Introduction

The analysis of the accuracy of results (usually called “responses”) computed by models relies fundamentally on the functional derivatives (usually called “sensitivities”) of the respective model responses with respect to the parameters in the respective computational model. Such sensitivities are needed for many purposes, including: (i) understanding the model by ranking the importance of the various parameters; (ii) performing “reduced-order modeling” by eliminating unimportant parameters and/or processes; (iii) quantifying the uncertainties induced in a model response due to model parameter uncertainties; (iv) performing “model validation,” by comparing computations to experiments to address the question “does the model represent reality?” (v) prioritizing improvements in the model; (vi) performing data assimilation and model calibration as part of forward “predictive modeling” to obtain best-estimate predicted results with reduced predicted uncertainties; (vii) performing inverse “predictive modeling”; (viii) designing and optimizing the system.

Response sensitivities are computed by using either deterministic or statistical methods. The simplest deterministic method for computing response sensitivities is to use finite-difference schemes in conjunction with re-computations using the model with “judiciously chosen” altered parameter values. Evidently, such methods can at best compute approximate values of a very limited number of sensitivities. Deterministic methods that can compute more exactly the values of first-order sensitivities include the “Green’s function method” [1], the “forward sensitivity analysis methodology” [2], and the “direct method” [3], which rely on analytical or numerical differentiation of the computational model under investigation to compute local response sensitivities exactly. However, for a computational model comprising many parameters, the conventional deterministic methods become impractical for computing sensitivities higher than first-order because they are subject to the “curse of dimensionality,” a term coined by Belmann [4] to describe phenomena in which the number of computations increases exponentially in the respective phase-space. In the particular case of sensitivity analysis using conventional deterministic methods, the number of large-scale computations increases exponentially in the phase-space of model parameter as the order of sensitivities increases.

The alternatives to “deterministic methods” are the “statistical methods”, which construct an approximate response distribution (often called “response surface”) in the parameters space, and subsequently use scatter plots, regression, rank transformation, correlations, and so-called “partial correlation analysis,” in order to identify approximate expectation values, variances and covariances for the responses. These statistical quantities are subsequently used to construct quantities that play the role of (approximate) first-order response sensitivities. Thus, statistical methods commence with “uncertainty analysis” and subsequently attempt an approximate “sensitivity analysis” of the approximately computed model response (called a “response surface”) in the phase-space of the parameters under consideration. The currently popular statistical methods for uncertainty and sensitivity analysis are broadly categorized as sampling-based methods [5,6], variance-based methods [7,8], and Bayesian methods [9]. Various variants of the statistical methods for uncertainty and sensitivity analysis are reviewed in the book edited by Saltarelli et al. [10]. The main advantage of using statistical methods for uncertainty and sensitivity analysis is that they are conceptually easy to implement. On the other hand, statistical methods for uncertainty and sensitivity analysis have three major inherent practical drawbacks, as follows:

(i)

Even first-order sensitivities cannot be computed exactly.

(ii)

Statistical methods are (also) subject to the curse of dimensionality and have not been developed for producing second- and higher-order sensitivities.

(iii)

Since the response sensitivities and parameter uncertainties are amalgamated, inherently and inseparably, within the results produced by statistical methods, improvements in parameter uncertainties cannot be directly propagated to improve response uncertainties; rather, the entire set of simulations and statistical post-processing must be repeated anew.

(iv)

A “fool-proof” statistical method for analyzing correctly models involving highly correlated parameters does not seem to exist currently, so that particular care must be used when interpreting regression results obtained using such models.

It is known that the “adjoint method of sensitivity analysis” has been the most efficient method for computing exactly first-order sensitivities, since it requires a single large-scale (adjoint) computation for computing all of the first-order sensitivities, regardless of the number of model parameters. The idea underlying the computation of response sensitivities with respect to model parameters using adjoint operators was first used by Wigner [11] to analyze first-order perturbations in nuclear reactor physics and shielding models based on the linear neutron transport (or diffusion) equation, as subsequently described in textbooks on these subjects [12,13,14,15,16]. Cacuci [2] is credited (see, e.g., [17,18]) for having conceived the rigorous “1^st-order adjoint sensitivity analysis methodology” for generic large-scale nonlinear (as opposed to linearized) systems involving generic operator responses and having introduced these principles to the earth, atmospheric and other sciences.

Cacuci [19,20] has extended his 1^st-order adjoint sensitivity analysis methodology to enable the comprehensive computation of 2^nd-order sensitivities of model responses to model parameters (including imprecisely known domain boundaries and interfaces) for large-scale linear and nonlinear systems. The 2^nd-order adjoint sensitivity analysis methodology for linear systems was applied [21] to compute exactly the 21,976 first-order sensitivities and 482,944,576 second-order sensitivities (of which 241,483,276 are distinct from each other) for an OECD/NEA reactor physics benchmark [22]. This benchmark is modeled by the neutron transport equation involving 21,976 uncertain parameters, the solving of which is representative of “large-scale computations.” The neutron transport equation was solved using the software package PARTISN [23] in conjunction with the MENDF71X cross section library [24], which comprises 618-group cross sections based on ENDF/B-VII.1 nuclear data [25]. The spontaneous fission source was computed using the code SOURCES4C [26]. This work has demonstrated that, contrary to the widely held belief that second- and higher-order sensitivities are negligeable for reactor physics systems, many 2^nd-order sensitivities of the OECD benchmark’s response to the benchmark’s uncertain parameters were much larger than the largest 1^st-order ones. This finding has motivated the investigation of the largest 3^rd-order sensitivities, many of which were found to be even larger than the 2^nd-order ones. Subsequently, the mathematical framework for determining and computing the 4^th-order sensitivities was developed, and many of these were found to be larger than the 3^rd-order ones. This sequence of findings has motivated the development by Cacuci [27] of the “n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (abbreviated as “n^th-CASAM-L”), which was developed specifically for linear systems because important model responses produced by such systems are various Lagrangian functionals which depend simultaneously on both the forward and adjoint state functions governing the respective linear system. Among the most important such responses are the Raleigh quotient for computing eigenvalues and/or separation constants when solving partial differential equations, and the Schwinger functional for first-order “normalization-free” solutions [28,29]. These functionals play a fundamental role in optimization and control procedures, derivation of numerical methods for solving equations (differential, integral, integro-differential), etc.

In parallel with developing the n^th-CASAM-L, Cacuci [30] has also developed the n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (n^th-CASAM-N). Just like the n^th-CASAM-L, the n^th-CASAM-N is also formulated in linearly increasing higher-dimensional Hilbert spaces (as opposed to exponentially increasing parameter-dimensional spaces), thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems, enabling the most efficient computation of exactly-determined expressions of arbitrarily high-order sensitivities of generic nonlinear system responses with respect to model parameters, uncertain boundaries and internal interfaces in the model’s phase-space.

Recently, Cacuci [31] has introduced the “Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems” (2^nd-FASAM-N), which enables a considerable reduction (by comparison to the 2^nd-CASAM-N) of the number of large-scale computations needed to compute the second-order sensitivities of a model response with respect to the model parameters, thereby becoming the most efficient methodology known for computing second-order sensitivities exactly. Paralleling the construction of the 2^nd-FASAM-N, this work introduces the “First- and Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Adjoint/Forward Linear Systems” (1^st and 2^nd-FASAM-L). The mathematical methodology of the 1^st-FASAM-L is presented in Section 3, while the mathematical methodology of the 2^nd-FASAM-L is presented in Section 4. The application of the 1^st-FASAM-L and the 2^nd-FASAM-L is illustrated in Section 5 by means of a simplified yet representative energy-dependent neutron-slowing down model which is fundamental importance to reactor physics and design [32,33,34]. The concluding discussion presented in Section 6 prepares the ground for the subsequent generalization of the present work to enable the most efficient possible computation of exact sensitivities of any (arbitrarily-high) order with respect to “feature functions” of model parameters and, hence, to the model’s parameters.

2. Mathematical Modeling of Response-Coupled Linear Forward and Adjoint Systems

The generic mathematical model considered in this work is fundamentally the same as was considered in [27], but with the major difference that functions (‘features”) of the primary model parameters will be generically identified within the model. The primary model parameters will be denoted as ${\alpha }_1$,…,${\alpha }_{TP}$, where the subscript “TP” indicates “Total number of Primary Parameters;” the qualifier “primary” indicates that these parameters do not depend on any other parameters within the model. These model parameters are considered to include imprecisely known geometrical parameters that characterize the physical system’s boundaries in the phase-space of the model’s independent variables. These boundaries depend on the physical system’s geometrical dimensions, which may be imprecisely known because of manufacturing tolerances. In practice, these primary model parameters are subject to uncertainties. It will be convenient to consider that these parameters are components of a “vector of primary parameters” denoted as $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}\in {ℝ}^{TP}$, where ${ℝ}^{TP}$ denotes the TP-dimensional subset of the set of real scalars. For subsequent developments, matrices and vectors will be denoted using capital and lower-case bold letters, respectively. The symbol “$\triangleq$” will be used to denote “is defined as” or “is by definition equal to.” Transposition will be indicated by a dagger $(†)$ superscript. The nominal parameter values will be denoted as ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\mathrm{...},{\alpha }_i^0,\mathrm{..},{\alpha }_{TP}^0\right]}^{†}$; the superscript “0” will be used throughout this work to denote “nominal” or “mean” values.

The model is considered to comprise $TI$ independent variables which will be denoted as $x_i,\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI$, and are considered to be the components of a $TI$-dimensional column vector denoted as $x\triangleq {\left(x_1,\dots ,x_{TI}\right)}^{†}\in {ℝ}^{TI}$, where the sub/superscript “$TI$” denotes the “Total number of Independent variables.” The vector $x\in {ℝ}^{TI}$ of independent variables is considered to be defined on a phase-space domain, denoted as $\Omega \left(\alpha \right)$, $\Omega \left(\alpha \right)\triangleq \left\{-\infty \le {\lambda }_i\left(\alpha \right)\le x_i\le {\omega }_i\left(\alpha \right)\le \infty ;\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI\right\}$, the boundaries of which may depend on some of the model parameters $\alpha$. The lower boundary-point of an independent variable is denoted as ${\lambda }_i\left(\alpha \right)$ (e.g., the inner radius of a sphere or cylinder, the lower range of an energy-variable, the initial time-value, etc.), while the corresponding upper boundary-point is denoted as ${\omega }_i\left(\alpha \right)$ (e.g., the outer radius of a sphere or cylinder, the upper range of an energy-variable, the final time-value, etc.). A typical example of boundary conditions that depend on imprecisely-known parameters that pertain to the geometry of the model and also to parameters that pertain to the material properties of the respective model occur when modeling particle diffusion within a medium, the boundaries of which are facing vacuum. For such models, the boundary conditions for the respective state (dependent) variables (i.e., particle flux and/or current) are imposed not on the physical boundary but on the “extrapolated boundary” of the respective spatial domain. The “extrapolated boundary” depends both on the imprecisely known physical dimensions of the medium’s domain/extent and also on the medium’s properties, i.e., atomic number densities and microscopic transport cross sections. The boundary of $\Omega \left(\alpha \right)$, which will be denoted as $\partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$, comprises the set of all of the endpoints ${\lambda }_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}{\omega }_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI,$of the respective intervals on which the components of $x$ are defined, i.e., $\partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]\triangleq \left\{{\lambda }_i\left(\alpha \right)\cup \hspace{0.17em}{\omega }_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI\right\}$.

The mathematical model that underlies the numerical evaluation of a process and/or state of a physical system comprises equations that relate the system's independent variables and parameters to the system's state/dependent variables. A linear physical system can generally be modeled by a system of coupled operator-equations as follows:

$$L\left[x;g\left(\alpha \right)\right]\phi \left(x\right)=Q\left[x;g\left(\alpha \right)\right]\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left(\alpha \right).$$

(1)

In Equation (1), the vector $\phi \left(x\right)\triangleq {\left[{\phi }_1\left(x\right),\dots ,{\phi }_{TD}\left(x\right)\right]}^{†}$ is a $TD$-dimensional column vector of dependent variables and where the sub/superscript “$TD$” denotes the “Total (number of) Dependent variables.” The functions ${\phi }_i\left(x\right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD$, denote the system’s “dependent variables” (also called “state functions”). The matrix $L\left(x;\alpha \right)\hspace{0.17em}\triangleq \left[L_{ij}\left(x;\alpha \right)\hspace{0.17em}\right],$ $\hspace{0.17em}i,j=1,\mathrm{...},TD$, has dimensions $TD\times TD$. The components $L_{ij}\left(x;\alpha \right)\hspace{0.17em}$ are operators that act linearly on the dependent variables ${\phi }_j\left(x\right)$ and also depend (in general, nonlinearly) on the uncertain model parameters $\alpha$. Furthermore, the vector $g\left(\alpha \right)\triangleq \left[g_1\left(\alpha \right),\dots ,g_{TG}\left(\alpha \right)\right]$ is a $TG$-dimensional vector having components $g_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TG$, which are real-valued functions of (some of) the primary model parameters $\alpha \in {ℝ}^{TP}$. The quantity $TG$ denotes the total number of such functions which appear exclusively in the definition of the model’s underlying equations. Such functions customarily appear in models in the form of correlations that describe “features” of the system under consideration, such as material properties, flow regimes. etc. Usually, the number of functions $g_i\left(\alpha \right)$ is considerably smaller than the total number of model parameters, i.e., $TG\ll TP$. For example, the numerical model (Cacuci and Fang, 2023) of the OECD/NEA reactor physics benchmark (Valentine 2006) comprises 21,976 uncertain primary model parameters (including microscopic cross sections and isotopic number densities) but the neutron transport equation, which is solved to determine the neutron flux distribution within the benchmark, does not use these primary parameters directly but instead uses just several hundreds of “group-averaged macroscopic cross sections” which are functions/features of the microscopic cross sections and isotopic number densities (which in turn are uncertain quantities that would be components of the vector of primary model parameters). In particular, a component $g_j\left(\alpha \right)$ may simply be one of the primary model parameters ${\alpha }_j$, i.e., $g_j\left(\alpha \right)\equiv {\alpha }_j$.

The $TD$-dimensional column vector $Q\left[x;g\left(\alpha \right)\right]\triangleq {\left(q_1,\mathrm{....},q_{TD}\right)}^{†}$, having components $q_i\left[x;g\left(\alpha \right)\right],\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD$, denotes inhomogeneous source terms, which usually depend nonlinearly on the uncertain parameters $\alpha$. Since the right-side of Equation (1) may contain distributions, the equality in this equation is considered to hold in the weak (i.e., “distributional”) sense. Similarly, all of the equalities that involve differential equations in this work will be considered to hold in the distributional sense.

When $L\left[x;g\left(\alpha \right)\right]$ contains differential operators, a set of boundary and initial conditions which define the domain of $L\left[x;g\left(\alpha \right)\right]$ must also be given. Since the complete mathematical model is considered to be linear in $\phi \left(x\right)$, the boundary and/or initial conditions needed to define the domain of $L\left[x;g\left(\alpha \right)\right]$ must also be linear in $\phi \left(x\right)$. Such linear boundary and initial conditions are represented in the following operator form:

$$B\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]\hspace{0.17em}\phi \left(x\right)=C\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$$

(2)

In Equation (2), the quantity $B\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$ denotes a matrix of dimensions $N_B\times TD$ having components denoted as $B_{ij}\left(x;\alpha \right);\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},N_B;\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TD$, which are operators that act linearly on $\phi \left(x\right)$ and nonlinearly on the components of $g\left(\alpha \right)$; the quantity $N_B$ denotes the total number of boundary and initial conditions. The $N_B$-dimensional column vector $C\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$ comprises components that are operators which, in general, act nonlinearly on the components of $g\left(\alpha \right)$.

Physical problems modeled by linear systems and/or operators are naturally defined in Hilbert spaces. The dependent variables ${\phi }_i\left(x\right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD$, for the physical system represented by Eqs. (1) and (2) are considered to be square-integrable functions of the independent variables and are considered to belong to a Hilbert space which will be denoted as $H_0\left(\Omega \right)$, where the subscript “zero” denotes “zeroth-level“ or “original.” Higher-level Hilbert spaces, which will be denoted as $H_1\left(\Omega \right)\hspace{0.17em}$ and $H_2\left(\Omega \right)$, will also be used in this work. The Hilbert space $H_0\left(\Omega \right)$ is considered to be endowed with the following inner product, denoted as ${\langle \phi \left(x\right),\psi \left(x\right)\rangle }_0$, between two elements $\phi \left(x\right)\in H_0\left(\Omega \right)$ and $\psi \left(x\right)\in H_0\left(\Omega \right)$:

$${\langle \phi \left(x\right),\psi \left(x\right)\rangle }_0\triangleq {\displaystyle \prod _{i=1}^{TI}{\displaystyle \underset{{\lambda }_i\left(\alpha \right)}{\overset{{\omega }_i\left(\alpha \right)}{\int }}\phi \left(x\right)·\psi \left(x\right)\hspace{0.17em}dx}}\hspace{0.17em}={\displaystyle \sum _{j=1}^{TD}{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_i\left(\alpha \right)}{\overset{{\omega }_i\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}{\phi }_j\left(x\right){\psi }_j\left(x\right)\hspace{0.17em}d{x}_1\hspace{0.17em}\mathrm{...}d{x}_i\mathrm{...}d{x}_{TI}}}}}.$$

(3)

The “dot” in Equation (3) indicates the “scalar product of two vectors,” which is defined as follows:

$$\phi \left(x\right)·\psi \left(x\right)\triangleq {\displaystyle \sum _{i=1}^{TD}{\phi }_i\left(x\right){\psi }_i\left(x\right)}.$$

(4)

The product-notation ${\displaystyle \prod _{i=1}^{TI}{\displaystyle \underset{{\lambda }_i\left(\alpha \right)}{\overset{{\omega }_i\left(\alpha \right)}{\int }}\left[\right]\hspace{0.17em}d{x}_i}}$ in Equation (3) denotes the respective multiple integrals.

The linear operator $L\left[x;g\left(\alpha \right)\right]$ is considered to admit an adjoint operator, which will be denoted as $L^{*}\left[x;g\left(\alpha \right)\right]$ and which is defined through the following relation for a vector $\psi \left(x\right)\in H_0$

$${\langle \psi \left(x\right),\hspace{0.17em}L\left[x;g\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em}\phi \left(x\right)\rangle }_0={\langle L^{*}\left[x;g\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em}\psi \left(x\right),\hspace{0.17em}\phi \left(x\right)\rangle }_0$$

(5)

In Equation (5), the formal adjoint operator $L^{*}\left[x;g\left(\alpha \right)\right]$ is the $TD\times TD$ matrix comprising elements $L_{ji}^{*}\left[x;g\left(\alpha \right)\right]$ which are obtained by transposing the formal adjoints of the forward operators $L_{ij}\left[x;g\left(\alpha \right)\right]$. Hence, the system adjoint to the linear system represented by (1) and (2) can generally be represented as follows:

$$L^{*}\left[x;g\left(\alpha \right)\right]\hspace{0.17em}\psi \left(x\right)=Q^{*}\left[x;g\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left(\alpha \right)\hspace{0.17em},$$

(6)

$$B^{*}\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]\hspace{0.17em}\psi \left(x\right)=C^{*}\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right].$$

(7)

When the forward operator $L\left[x;g\left(\alpha \right)\right]$ comprises differential operators, the operations (e.g., integration by parts) that implement the transition from the left-side to the right side of Equation (5) give rise to boundary terms which are collectively called the “bilinear concomitant.” The domain of $L^{*}\left[x;g\left(\alpha \right)\right]$ is determined by selecting adjoint boundary and/or initial conditions so as to ensure that the adjoint system is well-posed mathematically. It is also desirable that the selected adjoint boundary conditions should cause the bilinear concomitant to vanish when implemented in Equation (5) together with the forward boundary conditions given in Equation (2). The adjoint boundary conditions thus selected are represented in operator form by Equation (7).

The relationship shown in Equation (5), which is the basis for defining the adjoint operator, also provides the following fundamental “reciprocity-like” relation between the sources of the forward and the adjoint equations, i.e. Eqs. (1) and (6), respectively:

$${\langle \psi \left(x\right),\hspace{0.17em}Q\left[x;g\left(\alpha \right)\right]\rangle }_0={\langle Q^{*}\left[x;g\left(\alpha \right)\right],\hspace{0.17em}\phi \left(x\right)\rangle }_0$$

(8)

The functional on the right-side of Equation (8) represents a “detector response”, i.e., a reaction-rate between the particles and the medium represented by $Q^{*}\left[x;g\left(\alpha \right)\right]$ which is equivalent to the “number of counts” of particles incident on a detector of particles that “measures” the particle flux $\phi \left(x\right)$. In view of the relation provided in (8), the vector-valued source term $Q^{*}\left[x;g\left(\alpha \right)\right]\hspace{0.17em}\triangleq {\left\{q_1^{*}\left[x;g\left(\alpha \right)\right],\mathrm{....},q_{TD}^{*}\left[x;g\left(\alpha \right)\right]\right\}}^{†}$ in the adjoint equation Equation (6) is usually associated with the “result of interest” to be measured and/or computed, which is customarily called the system’s “response.” In particular, if $q_i^{*}\left[x;g\left(\alpha \right)\right]=\delta \left(x-x_d\right)$ and $q_{j\ne i}^{*}\left[x;g\left(\alpha \right)\right]=0$, then ${\langle Q^{*}\left[x;g\left(\alpha \right)\right],\hspace{0.17em}\phi \left(x\right)\rangle }_0={\phi }_i\left(x_d\right)$, which means that, in such a case, the right-side of Equation (8) provides the value of the i^th-dependent variable (particle flux, temperature, velocity, etc.) at the point in phase-space where the respective measurement is performed.

The results computed using a mathematical model are customarily called “model responses” (or “system responses” or “objective functions” or “indices of performance”). For linear physical systems, the system’s response may depend not only on the model’s state-functions and on the system parameters but may simultaneously also depend on the adjoint state function. As has been discussed by Cacuci [27,30], any response of a linear system can be formally represented (using expansions or interpolation, if necessary) and fundamentally analyzed in terms of the following generic integral representation:

$$R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]\triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]d{x}_1\mathrm{...}d{x}_{TI}}},$$

(9)

where $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ is a suitably differentiable nonlinear function of $\phi \left(x\right),\psi \left(x\right)$, and $\alpha$. The integral representation of the response provided in Equation (9) can represent “averaged” and/or “point-valued” quantities in the phase-space of independent variables. For example, if $R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$ represents the computation or the measurement (which would be a “detector-response”) of a quantity of interest at a point $x_d$ in the phase-space of independent variables, then $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ would contain a Dirac-delta functional of the form $\delta \left(x-x_d\right)$. Responses that represent “differentials/derivatives of quantities” would contain derivatives of Dirac-delta functionals in the definition of $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$. The vector $h\left(\alpha \right)\triangleq \left[h_1\left(\alpha \right),\dots ,h_{TH}\left(\alpha \right)\right]$, having components $h_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TH$, which appears among the arguments of the function $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$, represents functions of primary parameters that often appear solely in the definition of the response but do not appear in the mathematical definition of the model, i.e., in Eqs. (1), (2), (6) and (7). The quantity $TH$ denotes the total number of such functions which appear exclusively in the definition of the model’s response. Evidently, the response will depend directly and/or indirectly (through the “feature”-functions) on all of the primary model parameters. This fact has been indicated in Equation (9) by using the vector-valued function $f\left(\alpha \right)$ as an argument in the definition of the response $R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$ to represent the concatenation of all of the “features” of the model and response under consideration. The vector $f\left(\alpha \right)$ of “model features” is thus defined as follows:

$$f\left(\alpha \right)\triangleq {\left[g\left(\alpha \right);h\left(\alpha \right);\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]}^{†}\triangleq {\left[f_1\left(\alpha \right),\mathrm{...},f_{TF}\left(\alpha \right)\right]}^{†};\hspace{0.17em}\hspace{0.17em}TF\triangleq TG+TH+2TI.$$

(10)

As defined in Equation (10), the quantity $TF$ denotes the total number of “feature functions of the model’s parameters” which appear in the definition of the nonlinear model’s underlying equations and response.

Solving Eqs. (1) and (2), at the nominal (or mean) values, denoted as ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\mathrm{...},{\alpha }_i^0,\mathrm{..},{\alpha }_{TP}^0\right]}^{†}$, of the model parameters, yields the nominal forward solution, which will be denoted as ${\phi }^0\left(x\right)$. Solving Eqs. (6) and (7) at the nominal values, ${\alpha }^0$, of the model parameters yields the nominal adjoint solution, which will be denoted as ${\psi }^0\left(x\right)$. The nominal value of the response, $R\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right);f\left({\alpha }^0\right)\right]$, is determined by using the nominal parameter values ${\alpha }^0$, the nominal value ${\phi }^0\left(x\right)$ of the forward state function, and the nominal value ${\psi }^0\left(x\right)$ of the adjoint state function.

The definition provided by Equation (9) implies that the model response $R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$ depends on the components of the feature function $f\left(\alpha \right)$, and would therefore admit a Taylor-series expansion around the nominal value $f^0\triangleq f\left({\alpha }^0\right)$, having the following form:

$$R\left[f\left(\alpha \right)\right]=R\left(f^0\right)+{\displaystyle \sum _{j_1=1}^{TF}{\left\{\frac{\partial R\left(f\right)}{\partial f_{j_1}}\right\}}_{f^0}\delta f_{j_1}}+\frac{1}{2}{\displaystyle \sum _{j_1=1}^{TF}{\displaystyle \sum _{j_2=1}^{TF}{\left\{\frac{{\partial }^{2}R\left(f\right)}{\partial f_{j_1}\partial f_{j_2}}\right\}}_{f^0}}}\delta f_{j_1}\delta f_{j_2}+\mathrm{...}$$

(11)

where $\delta f_j\triangleq \left[f_j\left(\alpha \right)-f_j^0\right];\hspace{0.17em}\hspace{0.17em}{f}_j^0\triangleq f_j\left({\alpha }^0\right)\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TF$. The “sensitivities of the model response with respect to the (feature) functions” are naturally defined as being the functional derivatives of $R\left[f\left(\alpha \right)\right]$ with respect to the components (“features”) $f_j\left(\alpha \right)$ of $f\left(\alpha \right)$. The notation ${\{·\}}_{f^0}$ indicates that the quantity enclosed within the braces is to be evaluated at the nominal values $f^0\triangleq f\left({\alpha }^0\right)$. Since $TF\ll TP$, the computations of the functional derivatives of $R_k\left[f\left(\alpha \right)\right]$ with respect to the functions $f_j\left(\alpha \right)$, which appear in Equation (11), will be considerably less expensive computationally than the computation of the functional derivatives involved in the Taylor-series of the response with respect to the model parameters. The functional derivatives of the response with respect to the parameters can be obtained from the functional derivatives of the response with respect to the “feature” functions $f_j\left(\alpha \right)$ by simply using the chain rule, i.e.:

$${\left\{\frac{\partial R\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}={\displaystyle \sum _{i_1=1}^{TF}{\left\{\frac{\partial R\left(f\right)}{\partial f_{i_1}}\frac{\partial f_{i_1}\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}};\hspace{0.17em}\hspace{0.17em}{\left\{\frac{{\partial }^{2}R\left(\alpha \right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}}\right\}}_{{\alpha }^0}=\frac{\partial }{\partial {\alpha }_{j_2}}{\displaystyle \sum _{i_1=1}^{TF}{\left\{\frac{\partial R\left(f\right)}{\partial f_{i_1}}\frac{\partial f_{i_1}\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}};$$

(12)

and so on. The evaluation/computation of the functional derivatives $\partial f_{i_1}\left(\alpha \right)/\partial {\alpha }_{j_1}$, ${\partial }^2{f}_{i_1}\left(\alpha \right)/\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}$, etc., does not require computations involving the model, and is therefore computationally trivial by comparison to the evaluation of the functional derivatives (“sensitivities”) of the response with respect to either the functions (“features”) $f_j\left(\alpha \right)$ or the model parameters ${\alpha }_i,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TP$.

The range of validity of the Taylor-series shown in Equation (11) is defined by its radius of convergence. The accuracy −as opposed to the “validity”− of the Taylor-series in predicting the value of the response at an arbitrary point in the phase-space of model parameters depends on the order of sensitivities retained in the Taylor-expansion: the higher the respective order, the more accurate the respective response value predicted by the Taylor-series. In the particular cases when the response happens to be a polynomial function of the “feature” functions $f_j\left(\alpha \right)$, the Taylor series represented by Equation (11) is finite and exactly represents the respective model response.

In turn, the functions $f_i\left(\alpha \right)$ can also be formally expanded in a multivariate Taylor-series around the nominal (mean) parameter values ${\alpha }^0$, namely:

$$\begin{array}{l}{f}_i\left(\alpha \right)=f_i\left({\alpha }^0\right)+{\displaystyle \sum _{j_1=1}^{TP}{\left\{\frac{\partial f_i\left(\alpha \right)}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}\delta {\alpha }_{j_1}}+\frac{1}{2}{\displaystyle \sum _{j_1=1}^{TP}{\displaystyle \sum _{j_2=1}^{TP}{\left\{\frac{{\partial }^2{f}_i\left(\alpha \right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}}\right\}}_{{\alpha }^0}}}\delta {\alpha }_{j_1}\delta {\alpha }_{j_2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{3!}{\displaystyle \sum _{j_1=1}^{TP}{\displaystyle \sum _{j_2=1}^{TP}{\displaystyle \sum _{j_3=1}^{TP}{\left\{\frac{{\partial }^3{f}_i\left(\alpha \right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}\partial {\alpha }_{j_3}}\right\}}_{{\alpha }^0}}}}\delta {\alpha }_{j_1}\delta {\alpha }_{j_2}\delta {\alpha }_{j_3}+\mathrm{...},\end{array}$$

(13)

The choice of feature functions $f_i\left(\alpha \right)$ is not unique but can be tailored by the user to the problem at hand. The two most important guiding principles for constructing the feature functions $f_i\left(\alpha \right)$ based on the primary parameters are as follows:

(i)

As will be shown below in Section 4 while establishing the mathematical framework underlying the 2^nd-FASAM-L, the number of large-scale computations needed to determine the numerical value of the second-order sensitivities is proportional to the number of first-order sensitivities of the model’s response with respect to the feature functions f_i(α). Consequently, it is important to minimize the number of feature functions f_i(α), while ensuring that all of the primary model parameters are considered within the expressions constructed for the feature functions f_i(α). In the extreme case when some primary parameters, α_j, cannot be grouped into the expressions of the feature functions f_i(α), each of the respective primary model parameters α_j becomes a feature function f_j(α).

(ii)

The expressions of the features functions f_i(α) must be independent of the model’s state functions; they must be exact, closed-form, scalar-valued functions of the primary model parameters α_j, so the exact expressions of the derivatives of f_i(α) with respect to the primary model parameters α_j can be obtained analytically (with “pencil and paper”) and, hence, inexpensively from a computational standpoint. The motivation for this requirement is to ensure that the numerical determination of the subsequent derivatives of the features functions f_i(α) with respect to the primary model parameters α_j becomes trivial computationally.

The domain of validity of the Taylor-series in Equation (13) is defined by its own radius of convergence. Of course, in the extreme case when no feature function can be constructed, the feature functions will be the primary parameters themselves, in which case the n^th-FASAM-L methodology becomes identical to the previously established n^th-CASAM-L methodology [27].

3. The First-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward and Adjoint Linear Systems (1^st-FASAM-L)

The “First-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (1^st-FASAM-L) aims at enabling the most efficient computation of the first-order sensitivities of a generic model response of the form $R\left[\phi \left(x\right),\psi \left(x\right);\alpha \right]$ with respect to the components of the “features” function $f\left(\alpha \right)$. In preparation for subsequent generalizations towards establishing the generic pattern for computing sensitivities of arbitrarily high-order, the function $u^{\left(1\right)}\left(2;x\right)\triangleq {\left[\phi \left(x\right),\psi \left(x\right)\right]}^{†}$ will be called the “1^st-level forward/adjoint function” and the system of equations satisfied by this function (which is obtained by concatenating the original forward and adjoint equations together with their respective boundary/initial conditions) will be called “the 1^st-Level Forward/Adjoint System (1^st-LFAS)” and will be re-written in the following concatenated matrix-form:

$$F^{\left(1\right)}\left[2\times 2;x;f\right]u^{\left(1\right)}\left(2;x\right)=\hspace{0.17em}{q}_F^{\left(1\right)}\left(2;x;f\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left(\alpha \right)\hspace{0.17em};$$

(14)

$$\hspace{0.17em}{b}_F^{\left(1\right)}\left[u^{\left(1\right)}\left(2;x\right);\hspace{0.17em}f\right]=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right];$$

(15)

where the following definitions were used:

$$F^{\left(1\right)}\left[2\times 2;x;f\right]\triangleq \left(\begin{array}{cc}L\left(x;f\right)& 0\\ 0& L^{*}\left(x;f\right)\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right)\triangleq {\left[\phi \left(x\right),\psi \left(x\right)\right]}^{†};$$

(16)

$$q_F^{\left(1\right)}\left(2;x;f\right)\triangleq \hspace{0.17em}\left(\begin{array}{c}Q\left(x;g\right)\\ Q^{*}\left(x;g\right)\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{b}_F^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);\hspace{0.17em}f\right]\triangleq \hspace{0.17em}\left(\begin{array}{c}B\left(x;f\right)\phi \left(x\right)-C\left(f\right)\\ B^{*}\left(x;f\right)\psi \left(x\right)-C^{*}\left(f\right)\end{array}\right).$$

(17)

In the list of arguments of the matrix $F^{\left(1\right)}\left[2\times 2;x;f\right]$, the argument “$2\times 2$” indicates that this square matrix comprises four component sub-matrices, as indicated in Equation (16). Similarly, the argument “2” that appears in the block-vectors $\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right)$, $q_F^{\left(1\right)}\left(2;x;f\right)$, and $\hspace{0.17em}{b}_F^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);\hspace{0.17em}f\right]$ indicates that each of these column block-vectors comprises two sub-vectors as components. Also, throughout this work, the quantity “$0$” will be used to denote either a vector or a matrix with zero-valued components, depending on the context. For example, the vector “$0$” in Equation (15) is considered to have as many components as the vector $\hspace{0.17em}{b}_F^{\left(1\right)}\left[u^{\left(1\right)}\left(2;x\right);\hspace{0.17em}f\right]$. On the other hand, the quantity “$0$” which appears in Equation (16) may represent either a (sub) matrix or a vector of the requisite dimensions.

The nominal (or mean) parameter values, ${\alpha }^0$, are considered to be known, but these values will differ from the true values $\alpha$, which are unknown, by variations $\delta \alpha \triangleq {\left(\delta {\alpha }_1,\dots ,\delta {\alpha }_{TP}\right)}^{†}$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$. The parameter variations $\delta \alpha$ will induce variations $\delta f\left(\alpha \right)\triangleq {\left[\delta f_1\left(\alpha \right),\dots ,\delta f_{TF}\left(\alpha \right)\right]}^{†}$ in the vector-valued function $f\left(\alpha \right)$, around the nominal value $f^0\triangleq f\left({\alpha }^0\right)$, and will also induce variations $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$, respectively, around the nominal solution $\left({\phi }^0,{\psi }^0\right)$ through the equations underlying the model. All of these variations will induce variations in the model response $R\left[u^{\left(1\right)}\left(2;x\right);f\right]\equiv R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$, in a neighborhood $\left[{\phi }^0\left(x\right)+\epsilon \delta \phi \left(x\right),{\psi }^0\left(x\right)+\epsilon \delta \psi \left(x\right);f^0+\epsilon \delta f\right]$ around $\left({\phi }^0,{\psi }^0;f^0\right)$, where $\epsilon$ is a real-valued scalar.

Formally, the first-order sensitivities of the response $R\left[u^{\left(1\right)}\left(2;x\right);f\right]$ with respect to the components of the feature function $f\left(\alpha \right)$ are provided by the first-order Gateaux (G-)variation of $R\left(\phi ,\psi ,f\right)$ at the phase-space point $\left({\phi }^0,{\psi }^0,f^0\right)$, which is defined as follows:

$$\begin{array}{l}\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)\triangleq {\left\{\frac{d}{d\epsilon }R\left[{\phi }^0\left(x\right)+\epsilon \delta \phi \left(x\right),{\psi }^0\left(x\right)+\epsilon \delta \psi \left(x\right);f^0+\epsilon \delta f\right]\right\}}_{\epsilon =0}\\ \equiv {\left\{\frac{d}{d\epsilon }R\left[u^{\left(1,0\right)}\left(2;x\right)+\epsilon v^{\left(1\right)}\left(2;x\right);f^0+\epsilon \delta f\right]\right\}}_{\epsilon =0}\equiv \delta R\left[u^{\left(1,0\right)}\left(2;x\right);f^0;v^{\left(1\right)}\left(2;x\right),\delta f\right],\end{array}$$

(18)

where the following definitions were used:

$$u^{\left(1,0\right)}\left(2;x\right)\triangleq {\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right)\right]}^{†};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{v}^{\left(1\right)}\left(2;x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \right]}^{†}.$$

(19)

In general, the G-variation $\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)$ is nonlinear in the variations $\delta f\left(\alpha \right)$, $\delta \phi \left(x\right)$ and/or $\delta \psi \left(x\right)$. In such cases, the partial functional Gateaux (G-)derivatives of the response $R\left(\phi ,\psi ,f\right)$ with respect to the functions $\phi ,\psi ,f$ do not exist, which implies that the response sensitivities to the model parameters do not exist, either. Therefore, it will be henceforth assumed in this work that $\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)$ is linear in the respective variations, so the corresponding partial G-derivatives exist and $\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)$ is actually the first-order G-differential of the response. The usual numerical methods (e.g., Newton’s method and variants thereof) for solving the equations underlying the model also require the existence of the first-order G-derivatives of the original model equations; these will also be assumed to exist. When the 1^st-order G-derivatives exists, the G-differential $\delta R\left[u^{\left(1,0\right)}\left(2;x\right);f^0;v^{\left(1\right)}\left(2;x\right),\delta f\right]$ can be written as follows:

$$\begin{array}{l}\delta R\left[u^{\left(1,0\right)}\left(2;x\right);f^0;v^{\left(1\right)}\left(2;x\right),\delta f\right]={\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{dir}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}.\end{array}$$

(20)

In Equation (20), the “direct-effect” term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{dir}$ comprises only dependencies on $\delta f\left(\alpha \right)$ and is defined as follows:

$${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{dir}\triangleq {\left\{\frac{\partial R\left(u^{\left(1\right)};f\right)}{\partial f}\delta f\hspace{0.17em}\right\}}_{{\alpha }^0}$$

(21)

The following convention/definition was used in Equation (21):

$$\frac{\partial \left[\right]}{\partial f}\delta f\triangleq {\displaystyle \sum _{i=1}^{TF}\frac{\partial \left[\right]}{\partial f_i}\delta f_i}={\displaystyle \sum _{i=1}^{TG}\frac{\partial \left[\right]}{\partial g_i}\delta g_i}+{\displaystyle \sum _{i=1}^{TH}\frac{\partial \left[\right]}{\partial h_i}\delta h_i}+{\displaystyle \sum _{i=1}^{TI}\frac{\partial \left[\right]}{\partial {\omega }_i}\delta {\omega }_i}+{\displaystyle \sum _{i=1}^{TI}\frac{\partial \left[\right]}{\partial {\lambda }_i}\delta {\lambda }_i}$$

(22)

The above convention implies that:

(a) For $j=1,\mathrm{...},TG$:

$$\frac{\partial R\left(u^{\left(1\right)};f\right)}{\partial f_j\hspace{0.17em}}\delta f_j\hspace{0.17em}\triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}\frac{\partial S\left(\phi ,\psi ;g;h\right)}{\partial g_i}d{x}_1\mathrm{...}d{x}_{TI}}}\right\}}_{{\alpha }^0}\delta g_i;\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TG;$$

(23)

(b) For $j=TG+1,\mathrm{...},TG+TH$:

$$\frac{\partial R\left(u^{\left(1\right)};f\right)}{\partial f_j\hspace{0.17em}}\delta f_j\hspace{0.17em}\triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}\frac{\partial S\left(\phi ,\psi ;g;h\right)}{\partial h_i}d{x}_1\mathrm{...}d{x}_{TI}}}\right\}}_{{\alpha }^0}\delta h_i;\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TH\hspace{0.17em};$$

(24)

(c) For $j=TG+TH+1,\mathrm{...},TG+TH+TI$:

$$\frac{\partial R\left(u^{\left(1\right)};f\right)}{\partial f_j\hspace{0.17em}}\delta f_j\hspace{0.17em}\triangleq {\left\{\frac{\partial }{\partial {\omega }_i}{\displaystyle \underset{{\lambda }_1}{\overset{{\omega }_1}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}}{\overset{{\omega }_{TI}}{\int }}d{x}_{TI}S\left(\phi ,\psi ;g;h\right)}}\right\}}_{{\alpha }^0}\delta {\omega }_i\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI;$$

(25)

(d) For $j=TG+TH+TI+1,\mathrm{...},TG+TH+2TI$:

$$\frac{\partial R\left(u^{\left(1\right)};f\right)}{\partial f_j\hspace{0.17em}}\delta f_j\hspace{0.17em}\triangleq {\left\{\frac{\partial }{\partial {\lambda }_i}{\displaystyle \underset{{\lambda }_1}{\overset{{\omega }_1}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}}{\overset{{\omega }_{TI}}{\int }}d{x}_{TI}S\left(\phi ,\psi ;g;h\right)}}\right\}}_{{\alpha }^0}\delta {\lambda }_i\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI.$$

(26)

The notation on the left-side of Equation (22) represents the inner product between two vectors, but the symbol “($†$)” which indicates “transposition” has been omitted in order to keep the notation as simple as possible. “Daggers” indicating transposition will also be omitted in other inner products, whenever possible, while avoiding ambiguities.

In Equation (20), the “indirect-effect” term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}$ depends only on the variations $v^{\left(1\right)}\left(2;x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \right]}^{†}$ in the state functions, and is defined as follows:

$$\begin{array}{l}{\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}\triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}\frac{\partial S\left(\phi ,\psi ;g;h\right)}{\partial u^{\left(1\right)}\left(2;x\right)}{v}^{\left(1\right)}\left(2;x\right)}}\right\}}_{{\alpha }^0}\\ \triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}\frac{\partial S\left(\phi ,\psi ;g;h\right)}{\partial \phi }\delta \phi }}\right\}}_{{\alpha }^0}+{\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}\frac{\partial S\left(\phi ,\psi ;g;h\right)}{\partial \psi }\delta \psi }}\right\}}_{{\alpha }^0}.\end{array}$$

(27)

In Eqs. (21) and (27), the notation ${\left\{\right\}}_{{\alpha }^0}$ has been used to indicate that the quantity within the brackets is to be evaluated at the nominal values of the parameters and state functions. This simplified notation is justified by the fact that when the parameters take on their nominal values, it implicitly means that the corresponding state functions also take on their corresponding nominal values. This simplified notation will be used throughout this work.

The direct-effect term can be computed after having solved the forward system modeled by Eqs. (1) and (2), as well as the adjoint system modeled by Eqs. (6) and (7), to obtain the nominal values ${\phi }^0,{\psi }^0$ of the forward and adjoint dependent variables.

On the other hand, the indirect-effect term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}$ defined in Equation (27) can be quantified only after having determined the variations $v^{\left(1\right)}\left(2;x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \right]}^{†}$ in the state functions of the 1^st-Level Forward/Adjoint System (1^st-LFAS). The variations $v^{\left(1\right)}\left(2;x\right)$ are obtained as the solutions of the system of equations obtained by taking the first-order G-differentials of the 1^st-LFAS defined by Eqs. (14) and (15), which are obtained by definition as follows:

$${\left\{\frac{d}{d\epsilon }{F}^{\left(1\right)}\left[2\times 2;x;f^0+\epsilon \delta f\right]\left[u^{\left(1,0\right)}\left(2;x\right)+\epsilon v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{\epsilon =0}={\left\{\frac{d}{d\epsilon }{q}_F^{\left(1\right)}\left[2;x;f^0+\epsilon \delta f\right]\right\}}_{\epsilon =0}\hspace{0.17em},$$

(28)

$${\left\{\frac{d}{d\epsilon }\hspace{0.17em}{b}_F^{\left(1\right)}\left[2;u^{\left(1,0\right)}\left(2;x\right)+\epsilon v^{\left(1\right)}\left(2;x\right);\hspace{0.17em}{f}^0+\epsilon \delta f\right]\right\}}_{\epsilon =0}=0\left[2\right]\hspace{0.17em}.$$

(29)

Carrying out the differentiations with respect to $\epsilon$ in the above equations and setting $\epsilon =0$ in the resulting expressions yields the following matrix-vector equations:

$${\left\{V^{\left(1\right)}\left[2\times 2;x;f\right]v^{\left(1\right)}\left(2;x\right)\right\}}_{{\alpha }^0}={\left\{q_V^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left({\alpha }^0\right);$$

(30)

$${\left\{b_v^{\left(1\right)}\left(u^{\left(1\right)};v^{\left(1\right)};f\hspace{0.17em};\delta f\right)\right\}}_{{\alpha }^0}=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right];$$

(31)

where:

$$V^{\left(1\right)}\left[2\times 2;x;f\right]\triangleq \left(\begin{array}{cc}L\left(x;f\right)& 0\\ 0& L^{*}\left(x;f\right)\end{array}\right)=F^{\left(1\right)}\left[2\times 2;x;f\right];$$

(32)

$$q_V^{\left(1\right)}\left[2;u^{\left(1\right)};f;\delta f\right]\triangleq \hspace{0.17em}\left(\begin{array}{c}{q}_1^{\left(1\right)}\left(\phi ;f;\delta f\right)\\ q_2^{\left(1\right)}\left(\psi ;f;\delta f\right)\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{b}_v^{\left(1\right)}\left(u^{\left(1\right)};\hspace{0.17em}{v}^{\left(1\right)}\hspace{0.17em};f;\delta f\right)\triangleq \hspace{0.17em}\left(\begin{array}{c}{b}_1^{\left(1\right)}\left(\phi ;\delta \phi ;f;\delta f\right)\\ b_2^{\left(1\right)}\left(\psi ;\delta \psi ;f;\delta f\right)\end{array}\right);$$

(33)

$$$$

(34)

$$q_2^{(1)}\left(\psi ,f;\delta f\right)\triangleq \frac{\partial \left[Q^{*}-L^{*}\psi \left(x\right)\right]}{\partial f}\delta f\triangleq \hspace{0.17em}{\displaystyle \sum _{j_1=1}^{TF}{s}_2^{\left(1\right)}\left(j_1;\psi ;f\right)\delta f_{j_1}}$$

(35)

$$b_1^{\left(1\right)}\left(\phi ;\delta \phi ;f;\delta f\right)\triangleq B\delta \phi +\hspace{0.17em}\frac{\partial \left(B\phi -C\right)}{\partial f}\delta f;$$

(36)

$$b_2^{\left(1\right)}\left(\psi ;\delta \psi ;f;\delta f\right)\triangleq B^{*}\delta \psi +\hspace{0.17em}\frac{\partial \left(B^{*}\psi -C^{*}\right)}{\partial f}\delta f.$$

(37)

In order to keep the notation as simple as possible in Eqs. (30)‒(37), the differentials with respect to the various components of the feature function $f\left(\alpha \right)$ have all been written in the form $\left(\partial \left[\right]/\partial f\right)\delta f$, keeping in mind the convention/notation introduced in Equation (22). The system of equations comprising Eqs. (30) and (31) will be called the “1^st-Level Variational Sensitivity System (1^st-LVSS)” and its solution, $v^{\left(1\right)}\left(2;x\right)$, will be called the “1^st-level variational sensitivity function,” which is indicated by the superscript “(1)”. The solution, $v^{\left(1\right)}\left(2;x\right)$, of the 1^st-LVSS will be a function of the components of the vector of variations $\delta f$. In principle, therefore, if the response sensitivities with respect to the components of the feature function $f\left(\alpha \right)$ are of interest, then the 1^st-LVSS would need to be solved as many times as there are components in the variational features-function $\delta f$. On the other hand, if the response sensitivities with respect to the primary parameters are of interest, then the 1^st-LVSS would need to be solved as many times as there are primary parameters. Solving the 1^st-LVSS involves “large-scale computations.”

On the other hand, solving the 1^st-LVSS can be avoided altogether by using the ideas underlying the “adjoint sensitivity analysis methodology” originally conceived by Cacuci [2] and subsequently generalized by Cacuci [27,30] to enable the computation of arbitrarily high-order response sensitivities with respect to primary model parameters for both linear and nonlinear models. Thus, the need for solving repeatedly the 1^st-LVSS for every variation in the components of the feature function (or for every variation in the model’s parameters) is eliminated by expressing the indirect-effect term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}$ defined in Equation (27) in terms of the solutions of the “1^st-Level Adjoint Sensitivity System” (1^st-LASS), which will be constructed by implementing the following sequence of steps:

• Introduce a Hilbert space, denoted as $H_1$, comprising vector-valued elements of the form ${\chi }^{\left(1\right)}\left(2;x\right)\triangleq {\left[{\chi }_1^{\left(1\right)}\left(x\right),{\chi }_2^{\left(1\right)}\left(x\right)\right]}^{†}$, where the components ${\chi }_i^{\left(1\right)}\left(x\right)\triangleq {\left[{\chi }_{i,1}^{\left(1\right)}\left(x\right),\mathrm{...},{\chi }_{i,j}^{\left(1\right)}\left(x\right),\mathrm{...},{\chi }_{i,TD}^{\left(1\right)}\left(x\right)\right]}^{†}$, $i=1,2$, are square-integrable functions. Consider further that this Hilbert space is endowed with an inner product denoted as ${\langle {\chi }^{\left(1\right)}\left(2;x\right),{\theta }^{\left(1\right)}\left(2;x\right)\rangle }_1$ between two elements, ${\chi }^{\left(1\right)}\left(2;x\right)\in H_1$, ${\theta }^{\left(1\right)}\left(2;x\right)\in H_1$, which is defined as follows:

  $${\langle {\chi }^{\left(1\right)}\left(2;x\right),{\theta }^{\left(1\right)}\left(2;x\right)\rangle }_1\triangleq {\displaystyle \sum _{i=1}^2{\langle {\chi }_i^{\left(1\right)}\left(x\right),{\theta }_i^{\left(1\right)}\left(x\right)\rangle }_0}.$$

  (38)
• In the Hilbert $H_1$, form the inner product of Equation (30) with a yet undefined vector-valued function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}\in H_1$ to obtain the following relation:

  $${\left\{{\langle a^{\left(1\right)}\left(2;x\right),V^{\left(1\right)}\left[2\times 2;x;f^0\right]v^{\left(1\right)}\left(2;x\right)\rangle }_1\right\}}_{{\alpha }^0}={\left\{{\langle a^{\left(1\right)}\left(2;x\right),q_V^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}.$$

  (39)
• Using the definition of the adjoint operator in the Hilbert space $H_1$, recast the left-side of Equation (39) as follows:

  $$\begin{array}{l}{\left\{{\langle a^{\left(1\right)}\left(2;x\right),V^{\left(1\right)}\left[2\times 2;x;f\right]v^{\left(1\right)}\left(2;x\right)\rangle }_1\right\}}_{{\alpha }^0}={\left\{{\langle v^{\left(1\right)}\left(2;x\right),A^{\left(1\right)}\left[2\times 2;x;f\right]a^{\left(1\right)}\left(2;x\right)\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}\\ +{\left\{P^{\left(1\right)}\left[v^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0},\end{array}$$

  (40)

  where ${\left\{P^{\left(1\right)}\left[v^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}$ denotes the bilinear concomitant defined on the phase-space boundary $x\in \partial \Omega \left({\alpha }^0\right)$, and where $A^{\left(1\right)}\left[2\times 2;x;f\right]$ is the operator formally adjoint to $V^{\left(1\right)}\left[2\times 2;x;f\right]$, i.e.,

  $$A^{\left(1\right)}\left[2\times 2;x;f\right]\triangleq {\left\{V^{\left(1\right)}\left[2\times 2;x;f\right]\right\}}^{*}=\left(\begin{array}{cc}{L}^{*}\left(x;f\right)& 0\\ 0& L\left(x;f\right)\end{array}\right).$$

  (41)
• Require the first term on right-side of Equation (40) to represent the indirect-effect term defined in Equation (27),by imposing the following relation:

  $$A^{\left(1\right)}\left[2\times 2;x;f\right]a^{\left(1\right)}\left(2;x\right)\hspace{0.17em}=q_A^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right)\hspace{0.17em};f\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left({\alpha }^0\right);$$

  (42)

  where:

  $$q_A^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right)\hspace{0.17em};f\right]\triangleq {\left[\frac{\partial S\left(u^{\left(1\right)};f\right)}{\partial u^{\left(1\right)}\left(2;x\right)}\right]}^{†}\triangleq \left(\begin{array}{c}{\left[\partial S\left(\hspace{0.17em}{u}^{\left(1\right)};f\right)/\partial \phi \right]}^{†}\\ {\left[\partial S\left(\hspace{0.17em}{u}^{\left(1\right)};f\right)/\partial \psi \right]}^{†}\end{array}\right)\hspace{0.17em}\hspace{0.17em}.$$

  (43)
• Implement the boundary conditions represented by Equation (31) into Equation (40) and eliminate the remaining unknown boundary-values of the function $v^{\left(1\right)}\left(2;x\right)$ from the expression of the bilinear concomitant ${\left\{P^{\left(1\right)}\left[v^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}$ by selecting appropriate boundary conditions for the function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$, to ensure that Equation (42) is well-posed while being independent of unknown values of $v^{\left(1\right)}\left(2;x\right)$ and of $\delta f$. The boundary conditions thus chosen for the function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$ can be represented in operator form as follows:

  $${\left\{b_A^{\left(1\right)}\left[\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\right]\right\}}_{{\alpha }^0}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right].$$

  (44)
• The selection of the boundary conditions for $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$ represented by Equation (44) eliminates the appearance of the unknown values of $v^{\left(1\right)}\left(2;x\right)$ in ${\left\{P^{\left(1\right)}\left[v^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}$ and reduces this bilinear concomitant to a residual quantity that contains boundary terms involving only known values of $u^{\left(1\right)}\left(2;x\right)$, $a^{\left(1\right)}\left(2;x\right)$, $f$, and $\delta f$. This residual quantity will be denoted as ${\left\{{\widehat{P}}^{\left(1\right)}\left[\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}$. In general, this residual quantity does not automatically vanish, although it may do so occasionally.
• The system of equations comprising Equation (42) together with the boundary conditions represented Equation (44) will be called the 1^st-Level Adjoint Sensitivity System (1^st-LASS). The solution $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$ of the 1^st-LASS will be called the 1^st-level adjoint sensitivity function. The 1^st-LASS is called “first-level” (as opposed to “first-order”) because it does not contain any differential or functional-derivatives, but its solution, $a^{\left(1\right)}\left(2;x\right)$, will be used below to compute the first-order sensitivities of the response with respect to the components of the feature function $f\left(\alpha \right)$.
• Using Equation (39) together with the forward and adjoint boundary conditions represented by Eqs. (31) and (44) in Equation (40) reduces the latter to the following relation:

  $$\begin{array}{l}{\left\{{\langle a^{\left(1\right)}\left(2;x\right),q_V^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}={\left\{{\langle v^{\left(1\right)}\left(2;x\right),A^{\left(1\right)}\left[2\times 2;x;f\right]a^{\left(1\right)}\left(2;x\right)\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\left\{{\widehat{P}}^{\left(1\right)}\left[\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}.\end{array}$$

  (45)
• In view of Eqs. (27) and (42), the first term on the right-side of Equation (45) represents the indirect-effect term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\right]\right\}}_{ind}$. It therefore follows from Equation (45) that the indirect-effect term can be expressed in terms of the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$ as follows:

  $$\begin{array}{l}{\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}={\left\{{\langle a^{\left(1\right)}\left(2;x\right),q_V^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}\hspace{0.17em}\\ -{\left\{{\widehat{P}}^{\left(1\right)}\left[\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}\equiv {\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{ind}.\end{array}$$

  (46)

As indicated by the identity shown in Equation (46), the variations $\delta \phi$ and $\delta \psi$ have been eliminated from the original expression of the indirect-effect term, which now depends on the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$. Adding the expression obtained in Equation (46) with the expression for the direct-effect term defined in Equation (21) yields, according to Equation (20) the following expression for the total 1^st-order sensitivity ${\left\{\delta R\left(\phi ,\psi ,f;\delta \phi ,\delta \psi ,\delta f\right)\right\}}_{{\alpha }^0}$ of the response $R\left[\phi \left(x\right),\psi \left(x\right);f\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\begin{array}{l}{\left\{\delta R\left(\phi ,\psi ,f;\delta \phi ,\delta \psi ,\delta f\right)\right\}}_{{\alpha }^0}={\left\{\frac{\partial R\left(u^{\left(1\right)};f\right)}{\partial f}\delta f\hspace{0.17em}\right\}}_{{\alpha }^0}+{\left\{{\langle a^{\left(1\right)}\left(2;x\right),q_V^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}\\ \hspace{0.17em}-{\left\{{\widehat{P}}^{\left(1\right)}\left[\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}\equiv {\displaystyle \sum _{j_1=1}^{TF}{\left\{R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]\delta f_{j_1}\right\}}_{{\alpha }^0}.}\end{array}$$

(47)

The identity which appears in Equation (47) emphasizes the fact that the variations $\delta \phi$ and $\delta \psi$, which are expensive to compute, have been eliminated from the final expressions of the 1^st-order sensitivities $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$ of the response with respect to the components $f_{j_1}\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}{j}_1=1,\mathrm{...},TF$, of the features functions. The dependence on the variations $\delta \phi$ and $\delta \psi$ has been replaced in the expression of $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$ by the dependence on the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$. It is very important to note that the 1^st-LASS is independent of variations $\delta f\left(\alpha \right)$ in the components of the feature function and is consequently also independent of any variations $\delta \alpha$ in the primary model parameters. Hence, the 1^st-LASS needs to be solved only once to determine the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$. Subsequently, the “indirect-effect term” is computed efficiently and exactly by simply performing the integrations required to compute the inner product over the adjoint function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$, as indicated on the right-side of Equation (47). Solving the 1^st-Level Adjoint Sensitivity System (1^st-LASS) requires the same computational effort as solving the original coupled linear system, entailing the following operations: (i) inverting (i.e., solving) the left-side of the original adjoint equation with the source ${\left[\partial S\left(u^{\left(1\right)};\alpha \right)/\partial \phi \right]}^{†}$ to obtain the 1^st-level adjoint sensitivity function $a_1^{\left(1\right)}\left(x\right)$; and (ii) inverting the left-side of the original forward equation with the source ${\left[\partial S\left(u^{\left(1\right)};\alpha \right)/\partial \psi \right]}^{†}$ to obtain the 1^st-level adjoint sensitivity function $a_2^{\left(1\right)}\left(x\right)$.

The 1^st-order sensitivities $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$, $j_1=1,\mathrm{...},TF$, can be expressed as an integral over the independent variables as follows:

$$R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]\triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}}}{S}^{\left(1\right)}\hspace{0.17em}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right].$$

(48)

In particular, if the residual bilinear concomitant is non-zero, the functions $S^{\left(1\right)}\hspace{0.17em}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$ would contain suitably defined Dirac delta-functionals for expressing the respective non-zero boundary terms as volume-integrals over the phase-space of the independent variables. Dirac-delta functionals would also be used in the expression of $S^{\left(1\right)}\hspace{0.17em}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$ to represent terms containing the derivatives of the boundary end-points with respect to the model and/or response parameters.

The response sensitivities with respect to the primary model parameters would be obtained by using the expression obtained in Equation (48) in conjunction with the “chain rule” of differentiation provided in Equation (12).

It is important to compare the results produced by the 1^st-FASAM-L (for obtaining the sensitivities of the model response with respect to the model’s features) with the 1^st-CASAM (the 1^st-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems) methodology, which provides the expressions of the responses sensitivities directly with respect to the model’s primary parameters. Recall that the 1^st-CASAM-L [27] yields the following expression for the 1^st-order sensitivities of the response with respect to the primary model parameters:

$$\begin{array}{l}{\left\{\frac{\partial R\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);\alpha \right]}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}\hspace{0.17em}={\left\{{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}}}\frac{\partial \hspace{0.17em}S\left[u^{\left(1\right)}\left(2;x\right);\alpha \right]}{\partial {\alpha }_{j_1}}\right\}}_{{\alpha }^0}\\ +{\displaystyle \sum _{k=1}^{TI}{\displaystyle \prod _{m=1,k\ne j}^{TI}{\left\{{\displaystyle \underset{{\lambda }_m\left(\alpha \right)}{\overset{{\omega }_m\left(\alpha \right)}{\int }}d{x}_{m}S\hspace{0.17em}\left[u^{\left(1\right)}\left(2;\mathrm{...},{\omega }_k,\mathrm{...}\right);\alpha \right]\frac{\partial {\omega }_k\left(\alpha \right)}{\partial {\alpha }_{j_1}}-S\hspace{0.17em}\left[u^{\left(1\right)}\left(2;\mathrm{...},{\lambda }_k,\mathrm{...}\right);\alpha \right]\frac{\partial {\lambda }_k\left(\alpha \right)}{\partial {\alpha }_{j_1}}}\right\}}_{{\alpha }^0}}}\\ +{\left\{{\langle a^{\left(1\right)}\left(2;x\right),\frac{\partial }{\partial {\alpha }_{j_1}}{q}^{\left(1\right)}\left[u^{\left(1\right)}\left(2;x\right);\alpha \right]\hspace{0.17em}\rangle }_1\right\}}_{{\alpha }^0}-{\left\{\frac{\partial }{\partial {\alpha }_{j_1}}{\widehat{P}}^{\left(1\right)}\left[u^{\left(1\right)};a^{\left(1\right)};\alpha \right]\right\}}_{{\alpha }^0};\hspace{0.17em}\hspace{0.17em}{j}_1=1,\mathrm{...},TP.\end{array}$$

(49)

The same 1^st-level adjoint sensitivity function, denoted as $a^{\left(1\right)}\left(2;x\right)$, appears in Equation (49) as well as in Equation (48). Therefore, the same number of “large-scale computations” (which are needed to solve the 1^st-LASS to determine the 1^st-level adjoint sensitivity function) is needed for obtaining either the response sensitivities with respect to the components, $f_j\left(\alpha \right)$, $j=1,\mathrm{...},TF$, of the feature function $f\left(\alpha \right)$ using the c, or for obtaining the response sensitivities directly with respect to the primary model parameters ${\alpha }_j$, $j=1,\mathrm{...},TP$, by using the 1^st-CASAM-L. The use of the 1^st-CASAM-L would also require performing a number of $TP$ integrations to compute all of the response sensitivities with respect to the primary parameters; in contradistinction, the use of the 1^st-FASAM-L would require only $TF$ integrations ($TF<TP$) to compute all of the response sensitivities with respect to the components $f_j\left(\alpha \right)$ of the feature function. Since integrations using quadrature-scheme are significantly less expensive computationally by comparison to solving systems of equations (e.g., the original equations underlying the model and the 1^st-LASS), the computational savings provided by the use of the 1^st-FASAM-L is small by comparison to using the 1^st-CASAM-L. However, this conclusion is valid only for the computation of 1^st-order sensitivities. As will be shown in Section 4, below, the computational savings are significantly larger when computing the second-order sensitivities by using the 2^nd-FASAM-L rather than using the 2^nd-CASAM-L (or any other method).

4. The Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward and Adjoint Linear Systems (2^nd-FASAM-L)

The “Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems“ (2^nd-FASAM-L) determines the 2^nd-order sensitivities ${\partial }^{2}R\left[u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]/\partial f_{j_2}\partial f_{j_1}$ of the response with respect to the components of the “feature” function $f\left(\alpha \right)$ by conceptually considering that the first-order sensitivities $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]\triangleq \partial R\left[u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]/\partial f_{j_1}$, which were obtained in Equation (48), are “model responses.” Consequently, the 2^nd-order sensitivities are obtained as the “1^st-order sensitivities of the 1^st-order sensitivities” by applying the concepts underlying 1^st-FASAM to each 1^st-order sensitivity $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$, $j_1=1,\mathrm{...},TF$, which depends on both the vector $u^{\left(1\right)}\left(2;x\right)$, which comprises the original state variables, as well as on the 1^st-level adjoint function $a^{\left(1\right)}\left(2;x\right)$.

To establish the pattern underlying the computation of sensitivities of arbitrarily high-order, it is useful to introduce a systematic classification of the systems of equations that will underly the computation of the sensitivities of various orders. As has been shown in the Section 2.1, above, the 1^st-order response sensitivities $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$ depend on the original state functions $u^{\left(1\right)}\left(2;x\right)\triangleq {\left[\phi \left(x\right),\psi \left(x\right)\right]}^{†}$ and on the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(2;x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}$. The system of equations satisfied by these functions will be called “the 2^nd-Level Forward/Adjoint System (2^nd-LFAS)” and will be re-written in the following concatenated form:

$$F^{\left(2\right)}\left[2^2\times 2^2;f\left(\alpha \right)\right]u^{\left(2\right)}\left(2^2;x\right)=\hspace{0.17em}{q}_F^{\left(2\right)}\left[2^2;u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left(\alpha \right)\hspace{0.17em};$$

(50)

$$\hspace{0.17em}{b}_F^{\left(2\right)}\left(2^2;u^{\left(2\right)};\hspace{0.17em}f\right)\triangleq {\left(\hspace{0.17em}{b}_F^{\left(1\right)},b_A^{\left(1\right)}\right)}^{†}=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right];$$

(51)

where the following definitions were used:

$$F^{\left(2\right)}\left[2^2\times 2^2;f\left(\alpha \right)\right]\triangleq diag\left(F^{\left(1\right)},A^{\left(1\right)}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{u}^{\left(2\right)}\left(2^2;x\right)\triangleq {\left[u^{\left(1\right)}\left(2;x\right),a^{\left(1\right)}\left(2;x\right)\right]}^{†};$$

(52)

$$\hspace{0.17em}{q}_F^{\left(2\right)}\left[2^2;u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]\triangleq {\left[q_F^{\left(1\right)}\left(2;x;f\right),\hspace{0.17em}{q}_A^{\left(1\right)}\left[2;u^{\left(1\right)};f\right]\right]}^{†}.$$

(53)

The notation used for the matrix $F^{\left(2\right)}\left[2^2\times 2^2;f\left(\alpha \right)\right]$ indicates the following characteristics: (i) the superscript “2” indicates “2^nd-level”; (ii) the argument “$2^2\times 2^2$” indicates that this square matrix comprises 4x4=16 component sub-matrices. Similarly, the argument “2²” that appears in the block-vectors $u^{\left(2\right)}\left(2^2;x\right)$, $q_F^{\left(2\right)}\left[2^2;u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$ and $\hspace{0.17em}{b}_F^{\left(2\right)}\left(2^2;u^{\left(2\right)};\hspace{0.17em}\alpha \right)$ indicates that each of these column block-vectors comprises four sub-vectors as components.

The first-order G-differential of a first-order sensitivity $\hspace{0.17em}{R}^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);f\left(\alpha \right)\right]$, $j_1=\hspace{0.17em}1,\mathrm{...},TF$, is obtained by definition as follows:

$$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}\\ \triangleq {\left\{\frac{d}{d\epsilon }\delta R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right)+\epsilon v^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right)+\epsilon \delta a^{\left(1\right)}\left(2;x\right);f+\epsilon \delta f\right]\right\}}_{\epsilon =0}\\ ={\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right);f\right]\right\}}_{ind}+{\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);\delta f\right]\right\}}_{dir}\end{array}$$

(54)

The direct-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);\delta f\right]\right\}}_{dir}$ in Equation (54) is defined as follows

$$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);\delta f\right]\right\}}_{dir}\\ \triangleq {\displaystyle \sum _{j_2=1}^{TF}{\left\{\frac{\partial }{\partial f_{j_2}}{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}}}{S}^{\left(1\right)}\hspace{0.17em}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}}\delta f_{j_2}.\end{array}$$

(55)

and can be computed immediately. The indirect-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right);f\right]\right\}}_{ind}$ in Equation (54) depends on the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2^2;x\right)\triangleq \left[v^{\left(1\right)}\left(2;x\right),\delta a^{\left(1\right)}\left(2;x\right)\right]$ and is defined as follows:

$${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right);f\right]\right\}}_{ind}\triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}\left[s^{\left(2\right)}\left(2^2;j_1;u^{\left(2\right)};f\right)·v^{\left(2\right)}\left(2^2;x\right)\right]}}\hspace{0.17em},$$

(56)

where:

$$s^{\left(2\right)}\left(2^2;j_1;u^{\left(2\right)};f\right)\triangleq \partial R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right);f\right]/\partial u^{\left(2\right)}\hspace{0.17em}.$$

(57)

Evidently, the functions $v^{\left(1\right)}\left(2;x\right)$ and $\delta a^{\left(1\right)}\left(2;x\right)$ are needed in order to evaluate the above indirect-effect term. These functions are the solutions of the system of equations obtained by taking the first-G-differential of the 2^nd-LFAS defined by Eqs. (52) and (53). Applying the definition of the first G-differential the 2^nd-LFAS yields the following 2^nd-Level Variational Sensitivity System (2^nd-LVSS)” for the 2^nd-Level variational sensitivity function $v^{\left(2\right)}\left(2^2;x\right)\triangleq {\left[v^{\left(1\right)}\left(2;x\right),\delta a^{\left(1\right)}\left(2;x\right)\right]}^{†}$:

$$\begin{array}{l}{\left\{\frac{d}{d\epsilon }{F}^{\left(2\right)}\left[2^2\times 2^2;f^0+\epsilon \delta f\hspace{0.17em}\right]\left[u^{\left(2,0\right)}\left(2^2;x\right)+\epsilon v^{\left(2\right)}\left(2^2;x\right)\right]\right\}}_{\epsilon =0}\\ ={\left\{\frac{d}{d\epsilon }{q}_F^{\left(2\right)}\left[2^2;u^{\left(1,0\right)}\left(2;x\right)+\epsilon v^{\left(1\right)}\left(2;x\right);f^0+\epsilon \delta f\hspace{0.17em}\right]\right\}}_{\epsilon =0}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left({\alpha }^0\right)\hspace{0.17em};\end{array}$$

(58)

$${\left\{\frac{d}{d\epsilon }{b}_F^{\left(2\right)}\left[u^{\left(2,0\right)}\left(2^2;x\right)+\epsilon v^{\left(2\right)}\left(2^2;x\right);f^0+\epsilon \delta f\right]\right\}}_{\epsilon =0}\hspace{0.17em}=0;\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right].$$

(59)

Carrying out the differentiation with respect to $\epsilon$ in Eqs. (58) and (59), and setting $\epsilon =0$ in the resulting expressions yields the following 2^nd-LVSS:

$${\left\{V^{\left(2\right)}\left[2^2\times 2^2;x;f\right]v^{\left(2\right)}\left(2^2;x\right)\right\}}_{{\alpha }^0}={\left\{q_V^{\left(2\right)}\left[2^2;u^{\left(2\right)}\left(2^2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0};\hspace{0.17em}x\in \Omega \left({\alpha }^0\right);$$

(60)

$${\left\{b_v^{\left(2\right)}\left(u^{\left(2\right)};v^{\left(2\right)}\hspace{0.17em};f;\delta f\right)\right\}}_{{\alpha }^0}=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left[\lambda \left({\alpha }^0\right);\omega \left({\alpha }^0\right)\right];$$

(61)

where the following definitions were used:

$$\begin{array}{l}{V}^{\left(2\right)}\left[2^2\times 2^2;x;f\right]=\left(\begin{array}{cc}{V}^{\left(1\right)}\left[2\times 2;x;f\right]& 0\\ V_{21}^{\left(2\right)}\left(2\times 2;u^{\left(1\right)};f\right)& A^{\left(1\right)}\left[2\times 2;x;f\right]\end{array}\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ V_{21}^{\left(2\right)}\left(2\times 2;u^{\left(1\right)};f\right)\triangleq \left(\begin{array}{cc}-\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \phi \partial \phi }& -\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \phi \partial \psi }\\ -\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \psi \partial \phi }& -\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \psi \partial \psi }\end{array}\right)\hspace{0.17em};\end{array}$$

(62)

$$\begin{array}{l}{q}_V^{\left(2\right)}\left[2^2;u^{\left(2\right)}\left(2^2;x\right);f;\delta f\right]\hspace{0.17em}\triangleq \left(\begin{array}{c}{q}_V^{\left(1\right)}\left[2;u^{\left(1\right)};f;\delta f\right]\\ p_A^{\left(1\right)}\left[2;u^{\left(1\right)};a_1^{\left(1\right)};f;\delta f\right]\end{array}\right);\hspace{0.17em}\hspace{0.17em}\\ p_A^{\left(1\right)}\left[2;u^{\left(1\right)};a_1^{\left(1\right)};f;\delta f\right]\triangleq \hspace{0.17em}\left(\begin{array}{c}{p}_1^{\left(1\right)}\left(u^{\left(1\right)};a_1^{\left(1\right)};\delta f\right)\\ p_2^{\left(1\right)}\left(u^{\left(1\right)};a_1^{\left(1\right)};\delta f\right)\end{array}\right);\end{array}$$

(63)

$$p_1^{\left(1\right)}\left(u^{\left(1\right)};a_1^{\left(1\right)};f;\delta f\right)\triangleq \frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial f\partial \phi }\delta f-\frac{\partial \left[L^{*}\left(f\right)a_1^{\left(1\right)}\right]}{\partial f}\delta f;$$

(64)

$$p_2^{\left(1\right)}\left(u^{\left(1\right)};a_2^{\left(1\right)};f;\delta f\right)\triangleq \frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial f\partial \psi }\delta f-\frac{\partial \left[L\left(f\right)a_2^{\left(1\right)}\right]}{\partial f}\delta f;$$

(65)

$$b_v^{\left(2\right)}\left(u^{\left(2\right)};v^{\left(2\right)}\hspace{0.17em};f;\delta f\right)\triangleq \left(\begin{array}{c}\hspace{0.17em}{b}_v^{\left(1\right)}\left(u^{\left(1\right)};v^{\left(1\right)};f\hspace{0.17em};\delta f\right)\\ \delta b_A^{\left(1\right)}\left(u^{\left(2\right)};v^{\left(2\right)};f;\delta f\right)\end{array}\right);$$

(66)

The matrix $V_{21}^{\left(2\right)}\left(2\times 2;u^{\left(1\right)};f\right)$ depends only the system’s response and is responsible for coupling the forward and adjoint systems. Although the forward and adjoint systems are coupled, they could nevertheless be solved successively rather than simultaneously, because the matrix $V^{\left(2\right)}\left[2^2\times 2^2;x;f\right]$ is block-diagonal. All of the components of the matrices and vectors underlying the 2^nd-LVSS are to be computed at nominal parameter and state function values, as indicated in Eqs. (60) and (61).

Computing the indirect-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right)\hspace{0.17em};f\right]\right\}}_{ind}$ by solving the 2^nd-LVSS would require at least $2TF\left(TF+1\right)$ large-scale computations (to solve the 2^nd-LVSS) for every component of the feature function $f\left(\alpha \right)$.

The need for solving the 2^nd-LVSS will be circumvented by deriving an alternative expression for the indirect-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right)\hspace{0.17em};f\right]\right\}}_{ind}$, as defined in Equation (56), in which the second-level variational function $v^{\left(2\right)}\left(2^2;x\right)$ will be replaced by a 2^nd-level adjoint function which is independent of variations in the model parameter and state functions. This 2^nd-level adjoint function will be the solution of a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS), which will be constructed by using the same principles as employed for deriving the 1^st-LASS. The 2^nd-LASS is constructed in a Hilbert space, denoted as $H_2$, which will comprise as elements block-vectors of the same form as $v^{\left(2\right)}\left(2^2;x\right)$, i.e., a vector in $H_2$ has the generic structure${\chi }^{\left(2\right)}\left(2^2;x\right)\triangleq {\left[{\chi }_1^{\left(2\right)}\left(x\right),{\chi }_2^{\left(2\right)}\left(x\right),{\chi }_3^{\left(2\right)}\left(x\right),{\chi }_4^{\left(2\right)}\left(x\right)\right]}^{†}$, comprising four vector-valued components of the form ${\chi }_i^{\left(2\right)}\left(x\right)\triangleq {\left[{\chi }_{i,1}^{\left(2\right)}\left(x\right),\mathrm{...},{\chi }_{i,j}^{\left(2\right)}\left(x\right),\mathrm{...},{\chi }_{i,TD}^{\left(2\right)}\left(x\right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,3,4=2^2$. The inner product between two elements, ${\chi }^{\left(2\right)}\left(2;x\right)\in H_2$ and ${\theta }^{\left(2\right)}\left(2;x\right)\in H_2$, of this Hilbert space, will be denoted as ${\langle {\chi }^{\left(2\right)}\left(2;x\right),{\theta }^{\left(2\right)}\left(2;x\right)\rangle }_{2^2}$ and is defined as follows:

$${\langle {\chi }^{\left(2\right)}\left(2;x\right),{\theta }^{\left(2\right)}\left(2;x\right)\rangle }_{2^2}\triangleq {\displaystyle \sum _{i=1}^{2^2}{\langle {\chi }_i^{\left(2\right)}\left(x\right),{\theta }_i^{\left(2\right)}\left(x\right)\rangle }_0}$$

(67)

Note that there are $j_1=1,\mathrm{...},TF$ distinct indirect-effect terms ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right)\hspace{0.17em};f\right]\right\}}_{ind}$. Each of these indirect-effect terms will serve as a “source” for a “2^nd-Level Adjoint Sensitivity System (2^nd-LASS)” that will be constructed by applying the same sequence of steps that were used in Section 2.1, above, to construct the 1^st-LASS. This implies that a distinct 2^nd-level adjoint sensitivity function, of the form $a^{\left(2\right)}\left(2^2;j_1;x\right)\triangleq {\left[a_1^{\left(2\right)}\left(j_1;x\right),a_2^{\left(2\right)}\left(j_1;x\right),a_3^{\left(2\right)}\left(j_1;x\right),a_4^{\left(2\right)}\left(j_1;x\right)\right]}^{†}\in H_2$, $j_1=1,\mathrm{...},TF$, corresponding to each distinct indirect-effect term, will be needed for constructing each of the corresponding 2^nd-LASS, as follows:

• For each $j_1=1,\mathrm{...},TP$, form the inner product in the Hilbert space $H_2$ of Equation (60) with a yet undefined function $a^{\left(2\right)}\left(2^2;j_1;x\right)$ to obtain the following relation:

  $$\begin{array}{l}{\left\{{\langle a^{\left(2\right)}\left(2^2;j_1;x\right),V^{\left(2\right)}\left[2^2\times 2^2;x;f\right]v^{\left(2\right)}\left(2^2;x\right)\rangle }_{2^2}\right\}}_{{\alpha }^0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ ={\left\{{\langle a^{\left(2\right)}\left(2^2;j_1;x\right),q_V^{\left(2\right)}\left[2^2;u^{\left(2\right)}\left(2^2;x\right);f;\delta f\right]\hspace{0.17em}\rangle }_{2^2}\right\}}_{{\alpha }^0};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \Omega \left({\alpha }^0\right).\end{array}$$

  (68)
• Using the definition of the adjoint operator in the Hilbert space $H_2$, recast the left-side of Equation (68) as follows:

  $$\begin{array}{l}{\left\{{\langle a^{\left(2\right)}\left(2^2;j_1;x\right),V^{\left(2\right)}\left[2^2\times 2^2;x;f\right]v^{\left(2\right)}\left(2^2;x\right)\rangle }_{2^2}\right\}}_{{\alpha }^0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ ={\left\{{\langle v^{\left(2\right)}\left(2^2;x\right),A^{\left(2\right)}\left[2^2\times 2^2;x;f\right]a^{\left(2\right)}\left(2^2;j_1;x\right)\hspace{0.17em}\rangle }_{2^2}\right\}}_{{\alpha }^0}\\ +{\left\{P^{\left(2\right)}\left[v^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f;\delta f\right]\right\}}_{{\alpha }^0},\end{array}$$

  (69)

  where ${\left\{P^{\left(2\right)}\left[v^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}$ denotes the bilinear concomitant defined on the phase-space boundary $x\in \partial {\Omega }_x\left({\alpha }^0\right)$ and where $A^{\left(2\right)}\left[2^2\times 2^2;x;f\right]\triangleq {\left[V^{\left(2\right)}\left[2^2\times 2^2;x;f\right]\right]}^{*}$ is the operator formally adjoint to $V^{\left(2\right)}\left[2^2\times 2^2;x;f\right]$.
• The first term on right-side of Equation (69) is now required to represent the indirect-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right)\hspace{0.17em};f\right]\right\}}_{ind}$ defined in Equation (56). This requirement is satisfied by recalling Equation (57) and imposing the following relation on each function $a^{\left(2\right)}\left(2^2;j_1;x\right)$, $j_1=1,\mathrm{...},TF$:

  $${\left\{A^{\left(2\right)}\left[2^2\times 2^2;x;f\right]a^{\left(2\right)}\left(2^2;j_1;x\right)\right\}}_{{\alpha }^0}\hspace{0.17em}={\left\{s^{\left(2\right)}\left(2^2;j_1;u^{\left(2\right)};f\right)\right\}}_{{\alpha }^0}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{j}_1=1,\mathrm{...},TF,$$

  (70)
• The definition of the vector $a^{\left(2\right)}\left(2^2;j_1;x\right)$ will now be completed by selecting boundary conditions which will be represented in operator form as follows:

  $${\left\{b_A^{\left(2\right)}\left[u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f\right]\right\}}_{{\alpha }^0}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial \Omega \left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}{j}_1=1,\mathrm{...},TF.$$

  (71)
• The boundary conditions represented by Equation (71) are selected so as to satisfy the following requirements: (a) these boundary conditions together with Equation (70) constitute a well posed problem for the functions $a^{\left(2\right)}\left(2^2;j_1;x\right)$; (b) the implementation in Equation (69) of these boundary conditions together with those provided in Equation (61) eliminates all of the unknown values of the functions $v^{\left(2\right)}\left(2^2;x\right)$ and $a^{\left(2\right)}\left(2^2;j_1;x\right)$ in the expression of the bilinear concomitant ${\left\{P^{\left(2\right)}\left[v^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}$. This bilinear concomitant may vanish after these boundary conditions are implemented, but if it does not, it will be reduced to a residual quantity which will be denoted as ${\widehat{P}}^{\left(2\right)}\left[u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f;\delta f\right]$ and which will comprise only known values of $u^{\left(2\right)}\left(2^2;x\right)$, $a^{\left(2\right)}\left(2^2;j_1;x\right)$, $f$ and $\delta f$.

The system of equations represented by Equation (70) together with the boundary conditions represented by Equation (71) constitute the 2^nd-Level Adjoint Sensitivity System (2^nd-LASS). The solution of the 2^nd-LASS, i.e., the four-component vector $a^{\left(2\right)}\left(2^2;j_1;x\right)$, $j_1,\mathrm{...},TP$, will be called the 2^nd-level adjoint sensitivity function. It is important to note that the 2^nd-LASS is independent of any variations, $\delta f$, in the components of the feature function and, hence, is independent of any parameter variations, $\delta \alpha$, as well.

The equations underlying the 2^nd-LASS, represented by Eqs. (70), (71), together with the equations underlying the 2^nd-LVSS, represented by Eqs. (60) and (61), are now employed in Equation (69) in conjunction with Equation (56) to obtain the following expression for the indirect-effect term ${\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right)\hspace{0.17em};f\right]\right\}}_{ind}$ in terms of the 2^nd-level adjoint sensitivity functions $a^{\left(2\right)}\left(2^2;j_1;x\right)$, for $j_1=1,\mathrm{...},TP$:

$$\begin{array}{l}\hspace{0.17em}{\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);v^{\left(2\right)}\left(2^2;x\right)\hspace{0.17em};f\right]\right\}}_{ind}=\hspace{0.17em}{\left\{{\langle a^{\left(2\right)}\left(2^2;j_1;x\right),q_V^{\left(2\right)}\left[2^2;u^{\left(2\right)}\left(2^2;x\right);f;\delta f\right]\hspace{0.17em}\rangle }_2\right\}}_{{\alpha }^0}\\ -{\left\{{\widehat{P}}^{\left(2\right)}\left[u^{\left(2\right)};a^{\left(2\right)};f;\delta f\right]\right\}}_{{\alpha }^0}\equiv {\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right)\hspace{0.17em};f;\delta f\right]\right\}}_{ind}.\end{array}$$

(72)

As the last equality (identity) in Equation (72) indicates, the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2^2;x\right)$ has been eliminated from appearing in the expression of the indirect-effect term, having been replaced by the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2^2;j_1;x\right)$, for each $j_1=1,\mathrm{...},TF$.

Inserting the expressions that define the vector $q_V^{\left(2\right)}\left[2^2;u^{\left(2\right)}\left(2^2;x\right);f;\delta f\right]$ from Equation (63)‒(65) into Equation (72) and adding the resulting expression for the indirect-effect term to the expression of the direct-effect term given in Equation (54) yields the following expression for the total second-order G-differential of the response $R\left[\phi \left(x\right),\psi \left(x\right);f\right]$:

$$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(j_1;2^2;x\right);f;\delta f\right]\right\}}_{{\alpha }^0}\\ ={\displaystyle \sum _{j_2=1}^{TF}{\left\{R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(j_1;2^2;x\right);f\right]\right\}}_{{\alpha }^0}\delta f_{j_2}}\hspace{0.17em},\end{array}$$

(73)

where $R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(x\right);a^{\left(2\right)}\left(j_1;x\right);f\right]\equiv {\partial }^{2}R\left[\phi \left(x\right),\psi \left(x\right);f\right]/\partial f_{j_1}\partial f_{j_2}$ denotes the 2^nd-order partial sensitivity of the response $R\left[\phi \left(x\right),\psi \left(x\right);f\right]$ with respect to the components $f_{j_2}\left(\alpha \right)$ of the feature function $f\left(\alpha \right)$, evaluated at the nominal parameter values ${\alpha }^0$, and has the following expression for $j_1,j_2=1,\mathrm{...},TP$:

$$\begin{array}{l}{R}^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(j_1;2^2;x\right);f\right]\\ ={\left\{\frac{\partial }{\partial f_{j_2}}{\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}}}{S}^{\left(1\right)}\hspace{0.17em}\left[j_1;u^{\left(2\right)}\left(2^2;x\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}\\ +{\left\{{\langle a_1^{\left(2\right)}\left(j_1;x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial \left[Q\left(f\right)-L\left(f\right)\phi \left(x\right)\right]}{\partial f_{j_2}}\hspace{0.17em}\rangle }_0\right\}}_{{\alpha }^0}\\ +\hspace{0.17em}\hspace{0.17em}{\left\{{\langle a_2^{\left(2\right)}\left(j_1;x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial \left[Q^{*}\left(f\right)-L^{*}\left(f\right)\psi \left(x\right)\right]}{\partial f_{j_2}}\rangle }_0\right\}}_{{\alpha }^0}\hspace{0.17em}\\ +\hspace{0.17em}{\left\{{\langle a_3^{\left(2\right)}\left(j_1;x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{{\partial }^{2}S\left(u^{\left(1\right)}\left(x\right);f\right)}{\partial f_{j_2}\partial \phi }-\frac{\partial \left[L^{*}\left(f\right)a_1^{\left(1\right)}\right]}{\partial f_{j_2}}\hspace{0.17em}\rangle }_0\right\}}_{{\alpha }^0}\\ +\hspace{0.17em}\hspace{0.17em}{\left\{{\langle a_4^{\left(2\right)}\left(j_1;x\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{{\partial }^{2}S\left(u^{\left(1\right)}\left(x\right);f\right)}{\partial f_{j_2}\partial \psi }-\frac{\partial \left[L\left(f\right)a_2^{\left(1\right)}\right]}{\partial f_{j_2}}\hspace{0.17em}\rangle }_0\right\}}_{{\alpha }^0}\\ \hspace{0.17em}-\hspace{0.17em}{\left\{\frac{\partial }{\partial f_{j_2}}{\widehat{P}}^{\left(2\right)}\left[u^{\left(2\right)}\left(x\right);a^{\left(2\right)}\left(j_1;x\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}.\end{array}$$

(74)

Since the 2^nd-LASS is independent of variations in the components of the feature-functions (and, hence, variations in the model parameters), the exact computation of all of the partial second-order sensitivities $R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(j_1;2^2;x\right);f\right]\equiv {\partial }^{2}R\left[\phi \left(x\right),\psi \left(x\right);f\right]/\partial f_{j_1}\partial f_{j_2}$ requires at most $TF$ large-scale (adjoint) computations using the 2^nd-LASS. When the 2^nd-LASS is solved $TF$-times, the “off-diagonal” 2^nd-order mixed sensitivities ${\partial }^{2}R/\partial f_{j_1}\partial f_{j_2}$ will be computed twice, in two different ways, using two distinct 2^nd-level adjoint sensitivity functions, thereby providing an independent intrinsic (numerical) verification that the 1^st- and 2^nd-order response sensitivities with respect to the components of the feature functions are computed accurately. In component form, the equations comprising the 2^nd-LASS are solved, for each $\hspace{0.17em}{j}_1=1,\mathrm{...},TF$, in the following order:

$$L\left(f\right)a_3^{\left(2\right)}\left(j_1;x\right)=\frac{\partial S^{\left(1\right)}\left(j_1;u^{\left(2\right)};\hspace{0.17em}f\right)}{\partial a_1^{\left(1\right)}},$$

(75)

$$L^{*}\left(f\right)a_4^{\left(2\right)}\left(j_1;x\right)=\frac{\partial S^{\left(1\right)}\left(j_1;u^{\left(2\right)};f\right)}{\partial a_2^{\left(1\right)}},$$

(76)

$$L^{*}\left(f\right)a_1^{\left(2\right)}\left(j_1;x\right)=\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \phi \partial \phi }{a}_3^{\left(2\right)}\left(j_1;x\right)+\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \psi \partial \phi }{a}_4^{\left(2\right)}\left(j_1;x\right)+\frac{\partial S^{\left(1\right)}\left(j_1;u^{\left(2\right)};f\right)}{\partial \phi },$$

(77)

$$L\left(f\right)a_2^{\left(2\right)}\left(j_1;x\right)=\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \phi \partial \psi }{a}_3^{\left(2\right)}\left(j_1;x\right)+\frac{{\partial }^{2}S\left(u^{\left(1\right)};f\right)}{\partial \psi \partial \psi }{a}_4^{\left(2\right)}\left(j_1;x\right)+\frac{\partial S^{\left(1\right)}\left(j_1;u^{\left(2\right)};f\right)}{\partial \psi }.$$

(78)

Dirac delta-functionals may need to be used in Equation (74) in order to express in integral form the eventual non-zero residual terms in the residual bilinear concomitant and/or the terms containing derivatives with respect to the lower- and upper-boundary points. Ultimately, the expression of the partial second-order sensitivities $R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f\right]$obtained in Equation (74) is written in the following integral form, which mirrors Equation (48):

$$\begin{array}{l}{R}^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f\left(\alpha \right)\right]\\ \triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}d{x}_1\mathrm{...}{\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}d{x}_{TI}}}{S}^{\left(2\right)}\hspace{0.17em}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f\left(\alpha \right)\right].\end{array}$$

(79)

The computation of the partial second-order sensitivities $R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f\right]$ using Equation (74) requires quadratures for performing the integrations over the four components of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2^2;j_1;x\right)$, which are obtained by solving the 2^nd-LASS for $j_1=1,\mathrm{...},TF$. Thus, obtaining all of the second-order sensitivities $R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)}\left(2^2;x\right);a^{\left(2\right)}\left(2^2;j_1;x\right);f\right]\equiv {\partial }^{2}R/\partial f_{j_1}\partial f_{j_2}$ with respect to the components $f_{j_1}$ of the feature function $f\left(\alpha \right)$ requires performing at most $TF$ large-scale computations for solving the 2^nd-LASS.

By comparison, if the 2^nd-CASAM-L [27] had been applied to compute the second-order sensitivities of the response directly with respect to the model parameters, $TP$ (instead of $TF$) large-scale computations for solving the corresponding 2^nd-LASS would have been required, where $TP$ denotes the total number of primary model parameters. Since $TF<TP$, fewer large-scale computations are needed when using the 2^nd-FASAM-L rather than the 2^nd-CASAM-L. Notably, the left-sides of the 2^nd-LASS to be solved within the 2^nd-FASAM-L are the same as those to be solved within the 2^nd-CASAM-L. However, the source terms on the right-sides of these 2^nd-LASS are different from each other: there are as many source-terms on the right-sides as there are components of the feature function within the 2^nd-FASAM-L, and there are as many right-side sources as there are primary model parameters within the 2^nd-CASAM-L.

5. Illustrative High-Order Feature Adjoint Sensitivity Analysis of Energy-Dependent Particle Detector Response

The application of the n^th-FASAM-L methodology will be 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 methodology, this medium is considered to be infinitely large. The simplified neutron transport equation that models the energy-distributions of neutrons in such materials is called the “neutron slowing-down equation” and is written using the neutron lethargy (rather than the neutron energy) as the independent variable, which is denoted as “u” and is defined as follows: $u\triangleq \ln \left(E_0/E\right)$, where $E$ denotes the energy-variable and $E_0$ denotes the highest energy in the system. Thus, the neutron slowing-down model [32,33,34] for the energy-distribution of the neutron flux in a homogeneous mixture of non-fissionable materials of infinite extent takes on the following drastically simplified form of the neutron transport balance equation:

$$\frac{d\phi \left(u\right)}{du}+\frac{{\Sigma }_a}{\hspace{0.17em}\overline{\xi }\hspace{0.17em}{\Sigma }_t}\phi \left(u\right)=\frac{S\left(u\right)}{\overline{\xi }\hspace{0.17em}{\Sigma }_t};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<u\le u_{th};$$

(80)

$$\phi \left(0\right)=0;\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=0.$$

(81)

The quantities which appear in Equation (80) are defined below.

(1)

The lethargy-dependent neutron flux is denoted as $\phi \left(u\right)$; $u_{th}$ denotes a cut-off lethargy, usually taken to be the lethargy that corresponds to the thermal neutron energy (ca. 0.0024 electron-volts).

(2)

The macroscopic elastic scattering cross section for the homogeneous mixture of “M” materials is denoted as ${\Sigma }_s$ and is defined as follows:

$${\Sigma }_s\triangleq {\displaystyle \sum _{i=1}^M{N}_m^{\left(i\right)}}{\sigma }_s^{\left(i\right)}\hspace{0.17em};\hspace{0.17em}$$

(82)

where ${\sigma }_s^{\left(i\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M$ denotes the elastic scattering cross section of material “i”, and where the atomic or molecular number density of material “i” is denoted as $N_m^{\left(i\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M$ and is defined as follows: $N_m^{\left(i\right)}\triangleq {\rho }_i{N}_A/A_i$, where $N_A$ is Avogadro’s number $\left(0.602\times {10}^{24}\hspace{0.17em}nuclei/mole\right)$, while $A_i$ and ${\rho }_i$ denote the respective material’s mass number and density.

(3)

The average gain in lethargy of a neutron per collision is denoted as $\overline{\xi }$ and is defined as follows for the homogeneous mixture:

$$\overline{\xi }\triangleq \frac{1}{{\Sigma }_s}{\displaystyle \sum _{i=1}^M{\xi }_i{N}_m^{\left(i\right)}}{\sigma }_s^{\left(i\right)};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\xi }_i\triangleq 1+\hspace{0.17em}\frac{a_i\ln \hspace{0.17em}{a}_i}{1-\hspace{0.17em}{a}_i};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_i\triangleq {\left(\frac{A_i-1}{A_i+1}\right)}^2.$$

(83)

(4)

The macroscopic absorption cross section is denoted as ${\Sigma }_a$ and is defined as follows for the homogeneous mixture:

$${\Sigma }_a\triangleq {\displaystyle \sum _{i=1}^M{N}_m^{\left(i\right)}}{\sigma }_{\gamma }^{\left(i\right)}\hspace{0.17em},$$

(84)

where ${\sigma }_{\gamma }^{\left(i\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M$, denotes the microscopic radiative-capture cross section of material “i”.

(5)

The macroscopic total cross section is denoted as ${\Sigma }_t$ and is defined as follows for the homogeneous mixture:

$${\Sigma }_t\triangleq {\Sigma }_a+{\Sigma }_s.$$

(85)

(6)

The source $S\left(u\right)$ is considered to be a simplified “spontaneous fission” source stemming from fissionable actinides, such as ²³⁹Pu and ²⁴⁰Pu, emitting monoenergetic neutrons at the highest energy (i.e., zero lethargy). Such a source is comprised within the OECD/NEA polyethylene-reflected plutonium (PERP) OECD/NEA reactor physics benchmark [21,22] which can be modeled by the following simplified expression:

$$S\left(u\right)=S_0\delta \left(u\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{S}_0\triangleq {\displaystyle \sum _{k=1}^2{\lambda }_k^S{N}_k^S{F}_k^S{\nu }_k^S}{W}_k^S,$$

(86)

where the superscript “S” indicates “source;” the subscript index k=1 indicates material properties pertaining to the isotope ²³⁹Pu; the subscript index k=2 indicates material properties pertaining to the isotope ²⁴⁰Pu; ${\lambda }_k^S$ denotes the decay constant;$N_k^S$ denotes the atomic density of the respective actinide; $F_k^S$ denotes the spontaneous fission branching ratio; ${\nu }_k^S$ denotes the average number of neutrons per spontaneous fission; $W_k^S$ denotes a function of parameters used in a Watt’s fission spectrum to approximate the spontaneous fission neutron spectrum of the respective actinide. The detailed forms of the parameters $W_k^S$ are unimportant for illustrating the application of the n^th-FASAM-L methodology. The nominal values for these imprecisely known parameters are available from a library file contained in SOURCES4C [26].

The response considered for the above neutron slowing-down model is the reaction rate, denoted as $R$, of neutrons of energy $u=u_d$ that would be measured by a detector characterized by an interaction cross section ${\Sigma }_d\triangleq N_d{\sigma }_d$, where $N_d$ denotes the atomic or molecular number density of the detector’s material while ${\sigma }_d$ denotes the detector’s microscopic interaction cross section. Mathematically, the detector’s reaction rate can be represented by the following functional of the neutron flux $\phi \left(u\right)$:

$$R={\Sigma }_d\phi \left(u_d\right)={\Sigma }_d{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\Sigma }_d=N_d{\sigma }_d.$$

(87)

For this “source-detector” model, the following primary model parameters are subject to experimental uncertainties:

(i) the atomic number densities $N_m^{\left(i\right)}$; the microscopic radiative-capture cross section ${\sigma }_{\gamma }^{\left(i\right)}$; the scattering cross section ${\sigma }_s^{\left(i\right)}$, for each material “i”, $i=1,\mathrm{...},M$, included in the homogeneous mixture;

(ii) the source parameters ${\lambda }_k^S$,$N_k^S$, $F_k^S$, ${\nu }_k^S$, $W_k^S$, for k=1,2;

(iii) the atomic density $N_d$ and the microscopic interaction cross section ${\sigma }_d$ that characterize the detector’s material.

These above primary parameters are considered to constitute the components of a “vector of primary model parameters” defined as follows:

$$\begin{array}{l}\alpha \triangleq {\left(N_m^{\left(1\right)},{\sigma }_{\gamma }^{\left(1\right)},{\sigma }_s^{\left(1\right)},\mathrm{...},N_m^{\left(M\right)},{\sigma }_{\gamma }^{\left(M\right)},{\sigma }_s^{\left(M\right)},{\lambda }_1^S,{\lambda }_2^S,N_1^S,N_2^S,F_1^S,F_2^S,{\nu }_1^S,{\nu }_2^S,W_1^S,W_2^S,N_d,{\sigma }_d\right)}^{†}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}TP\triangleq 3M+12.\end{array}$$

(88)

On the other hand, the structure of the computational model comprising Eqs. (80), (81) and (87) suggests that the components $f_i\left(\alpha \right)$ of the feature function $f\left(\alpha \right)$ can be defined as follows:

$$\begin{array}{l}f\left(\alpha \right)\triangleq {\left[f_1\left(\alpha \right),f_2\left(\alpha \right),f_3\left(\alpha \right)\right]}^{†};\hspace{0.17em}\\ f_1\left(\alpha \right)\triangleq \frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{f}_2\left(\alpha \right)\triangleq \frac{S_0\left(\alpha \right)}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{f}_3\left(\alpha \right)\triangleq \hspace{0.17em}{\Sigma }_d\left(\alpha \right).\end{array}$$

(89)

Solving Eqs. (80) and (81) while using the definitions introduced in Equation (89) yields the following expression for the flux $\phi \left(u\right)$ in terms of the components $f_i\left(\alpha \right)$ of the feature function $f\left(\alpha \right)$:

$$\phi \left(u\right)=H\left(u\right)f_2\left(\alpha \right)\exp \left[-u{f}_1\left(\alpha \right)\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}H(0)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}H\left(u\right)=1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}u>0.$$

(90)

In terms of the components $f_i\left(\alpha \right)$ of the feature function $f\left(\alpha \right)$, the model’s response takes on the following expression:

$$R\left(\alpha \right)=f_3\left(\alpha \right)f_2\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}.$$

(91)

As Equation (91) indicates, the model response can be considered to depend directly on $TP\triangleq 3M+12$ primary model parameters. Alternatively, the model response can be considered to depend directly on 3 feature functions and only indirectly (through the 3 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 methodology. In the later consideration/interpretation, the response sensitivities to the primary model parameters will be obtained by applying the n^th-FASAM-L methodology, which will involve two stages, as follows: the response sensitivities with respect to the feature functions will be obtained in the first stage, while 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 within the n^th-CASAM-L methodology versus the n^th-FASAM-L methodology will become evident in the remainder of this Section by means of using the illustrative neutron slowing-down model, which is representative of the general situation for any linear system.

According to the “reciprocity relation” for linear systems highlighted in Equation (8), the detector response defined in Equation (87) can be alternatively expressed in terms of the solution of the “adjoint slowing-down model”, i.e., the model that would be adjoint to the forward slowing-down model represented by Eqs. (80) and (81). The “adjoint slowing-down model” is constructed in the Hilbert space $H_B$ of all square-integrable functions $\phi \left(u\right)\in H_B$, $\psi \left(u\right)\in H_B$ endowed with the following inner product, denoted as ${\langle \phi \left(u\right),\psi \left(u\right)\rangle }_B$:

$${\langle \phi \left(u\right),\psi \left(u\right)\rangle }_B\triangleq {\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\psi \left(u\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(92)

Using the inner product ${\langle \phi \left(u\right),\psi \left(u\right)\rangle }_B$ defined in Equation (92), the adjoint slowing-down model is constructed by the usual procedure, namely: (i) construct the inner product of Equation (80) with a function $\psi \left(u\right)\in H_B$; (ii) integrate by parts the resulting relation so as to transfer the differential operation from the forward function $\phi \left(u\right)$ onto the adjoint function $\psi \left(u\right)$; (iii) use the initial condition provided in Equation (81) and eliminate the unknown function $\phi \left(u_{th}\right)$ by choosing the final-value condition $\psi \left(u_{th}\right)=0$; (iv) choose the source for the resulting adjoint slowing-down model so as to satisfy the reciprocity relation shown in Equation (8). The result of these operations is the following adjoint slowing down model for the adjoint slowing-down function $\psi \left(u\right)$:

$$-\frac{d\psi \left(u\right)}{du}+f_1\left(\alpha \right)\psi \left(u\right)=f_3\left(\alpha \right)\delta \left(u-u_d\right)\hspace{0.17em},$$

(93)

$$\psi \left(u_{th}\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=u_{th}.$$

(94)

In terms of the adjoint slowing-down function $\psi \left(u\right)$, the detector response takes on the following alternative expression:

$$R=f_2\left(\alpha \right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\psi \left(u\right)\delta \left(u\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(95)

The correctness of the alternative expression for the detector response provided in Equation (95) can be readily verified by solving the adjoint slowing down equation to obtain the following closed form expression for the adjoint slowing-down function $\psi \left(u\right)$:

$$\psi \left(u\right)=H\left(u_d-u\right)f_3\left(\alpha \right)\exp \left[\left(u-u_d\right)f_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},$$

(96)

and subsequently inserting the above expression into Equation (95) to obtain the same final result as was obtained in Equation (91) in terms of the forward slowing-down flux $\phi \left(u\right)$.

5.1. First-Order Adjoint Sensitivity Analysis: 1^st-FASAM-L Versus 1^st-CASAM-L

In this Subsection, the computation of the first-order sensitivities of the response $R\left(\alpha \right)$ with respect to the primary model parameters will first be demonstrated by using the 1^st-FASAM-L. Subsequently, the same first-order sensitivities will be obtained by using the 1^st-CASAM-L and the two alternative paths will be compared to each other, showing that the same expressions are obtained for the respective sensitivities, as expected. Although the computational efforts are not identical, they are comparable in terms of efficiency, with a slight advantage for the 1^st-FASAM-L methodology.

5.1.1. Application of the 1^st-FASAM-L

The 1^st-FASAM-L will be applied to the neutron slowing-down paradigm illustrative model by following the principles outlined in Section 3. In this case, the model response is written in terms of the feature functions as follows:

$$R\left(\phi ,f\right)=f_3{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em},$$

(97)

where the flux $\phi \left(u\right)$ is the solution of the 1^st-Level Forward/Adjoint System (1^st-LFAS) comprising Eqs. (80) and (81), where Equation (80) is written in terms of the feature functions as follows:

$$\frac{d\phi \left(u\right)}{du}+f_1\phi \left(u\right)=f_2\delta \left(u\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<u\le u_{th}.$$

(98)

The first-order sensitivities of the response $R\left(\phi ,f\right)$ with respect to the components of the feature function $f\left(\alpha \right)$ are provided by the first-order Gateaux (G-)variation $\delta R\left({\phi }^0,f^0;v^{\left(1\right)},\delta f\right)$ of $R\left(\phi ,f\right)$, for variations $v^{\left(1\right)}\left(u\right)\triangleq \delta \phi \left(u\right)$ and $\delta f_3$ around the phase-space point $\left({\phi }^0,f^0\right)$, as shown in Equation (18), to obtain:

$$\begin{array}{l}\delta R\left({\phi }^0,f^0;v^{\left(1\right)},\delta f\right)\triangleq {\left\{\frac{d}{d\epsilon }\left[\left(f_3^0+\epsilon \delta f_3\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left({\phi }^0+\epsilon v^{\left(1\right)}\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\right]\right\}}_{\epsilon =0}\\ \triangleq {\left\{\delta R\left({\phi }^0,f^0;\delta f_3\right)\right\}}_{dir}+{\left\{\delta R\left({\phi }^0,f^0;v^{\left(1\right)}\right)\right\}}_{ind},\end{array}$$

(99)

where the “direct-effect” term ${\left\{\delta R\left({\phi }^0,f^0;\delta f_3\right)\right\}}_{dir}$ and, respectively, the “indirect-effect” term ${\left\{\delta R\left({\phi }^0,f^0;v^{\left(1\right)}\right)\right\}}_{ind}$ are defined as follows:

$${\left\{\delta R\left({\phi }^0,f^0;\delta f_3\right)\right\}}_{dir}\triangleq \left(\delta f_3\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{\phi }^0\left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em},$$

(100)

$${\left\{\delta R\left({\phi }^0,f^0;v^{\left(1\right)}\right)\right\}}_{ind}\triangleq f_3^0{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(101)

The “1^st-level variational sensitivity function” $v^{\left(1\right)}\left(u\right)$ is obtained as the solution of the “1^st-Level Variational Sensitivity System (1^st-LVSS)” obtained by taking the first-order G-differentials of the 1^st-LFAS defined by Eqs. (98) and (81), which are derived as shown in Eqs. (28) and (29), to obtain:

$${\left\{\frac{d}{d\epsilon }\left[\frac{d\left({\phi }^0+\epsilon v^{\left(1\right)}\right)}{du}+\left(f_1^0+\epsilon \delta f_1\right)\left({\phi }^0+\epsilon v^{\left(1\right)}\right)\right]\right\}}_{\epsilon =0}=\delta \left(u\right){\left\{\frac{d}{d\epsilon }\left(f_2^0+\epsilon \delta f_2\right)\right\}}_{\epsilon =0},$$

(102)

$${\left\{\frac{d}{d\epsilon }\left[{\phi }^0\left(u\right)+\epsilon v^{\left(1\right)}\left(u\right)\right]\right\}}_{\epsilon =0}=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=0.$$

(103)

Carrying out the differentiations with respect to $\epsilon$ in the above equations and setting $\epsilon =0$ in the resulting expressions yields the following 1^st-LVSS:

$$\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+f_1\left({\alpha }^0\right)v^{\left(1\right)}\left(u\right)=\left(\delta f_2\right)\delta \left(u\right)-\left(\delta f_1\right){\phi }^0\left(u\right)\hspace{0.17em},$$

(104)

$$v^{\left(1\right)}\left(u\right)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=0.$$

(105)

For further reference, the closed-form solution of the above 1^st-LVSS has the following expression:

$$v_1\left(u\right)={\left\{\left[\left(\delta f_2\right)-\left(\delta f_1\right)u{f}_2\right]H\left(u\right)\exp \left(-u{f}_1\right)\right\}}_{{\alpha }^0}.$$

(106)

In principle, the above expression for $v^{\left(1\right)}\left(u\right)$ could be used in Equation (101) to obtain the value if the indirect-effect term. In practice, however, the 1^st-LVSS cannot be solved analytically so the closed form expression of $v^{\left(1\right)}\left(u\right)$ is not available. Consequently, rather than (numerically) solve repeatedly the 1^st-LVSS for every possible variation induced by the primary parameters in the component feature functions, the alternative route to determining the expression of the indirect-effect term is to develop the 1^st-Level Adjoint sensitivity System (1^st-LASS) by following the procedure described in Section 3. The Hilbert space, denoted as $H_1$, appropriate for this illustrative model is the space of all square-integrable functions endowed with the following inner product, denoted as ${\langle {\chi }^{\left(1\right)}\left(u\right),{\theta }^{\left(1\right)}\left(u\right)\rangle }_1$, between two elements, ${\chi }^{\left(1\right)}\left(u\right)\in H_1$, ${\theta }^{\left(1\right)}\left(2;x\right)\in H_1$, belonging to this Hilbert space:

$${\langle {\chi }^{\left(1\right)}\left(u\right),{\theta }^{\left(1\right)}\left(u\right)\rangle }_1\triangleq {\displaystyle \underset{0}{\overset{u_{th}}{\int }}{\chi }^{\left(1\right)}\left(u\right){\theta }^{\left(1\right)}\left(u\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(107)

In this particular instance, the Hilbert space $H_1$ coincides with the original Hilbert space $H_B$ used for the original forward and adjoint slowing down models. More generally, similar situations occur when the response depends either just on the forward or just on the adjoint state function(s).

Using Equation (107), construct in the Hilbert space $H_1$ the inner product of Equation (104) with a square-integrable function $a^{\left(1\right)}\left(u\right)\in H_1$ to obtain the following relation:

$${\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+f_1{v}^{\left(1\right)}\left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}={\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\left(\delta f_2\right)\delta \left(u\right)-\left(\delta f_1\right)\phi \left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}.$$

(108)

Using the definition of the adjoint operator in the Hilbert space $H_1$, which in this case amounts to integration by parts of the left-side of Equation (108), obtain the following relation:

$$\begin{array}{l}{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+f_1{v}^{\left(1\right)}\left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}={\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)}\left[-\frac{d{a}^{\left(1\right)}\left(u\right)}{du}+f_1{a}^{\left(1\right)}\left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}\\ +{\left\{a^{\left(1\right)}\left(u_{th}\right)v^{\left(1\right)}\left(u_{th}\right)-a^{\left(1\right)}\left(0\right)v^{\left(1\right)}\left(0\right)\right\}}_{{\alpha }^0}.\end{array}$$

(109)

Require the first term on right-side of Equation (109) to represent the indirect-effect term defined in Equation (101), to obtain the following relation:

$${\left\{-\frac{d{a}^{\left(1\right)}\left(u\right)}{du}+f_1{a}^{\left(1\right)}\left(u\right)\hspace{0.17em}\right\}}_{{\alpha }^0}={\left\{f_3\delta \left(u-u_d\right)\right\}}_{{\alpha }^0};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<u\le u_{th}.$$

(110)

Implement the boundary conditions represented by Equation (105) into Equation (109) and eliminate the unknown boundary-value $v^{\left(1\right)}\left(u_{th}\right)$ from this relation by imposing the following boundary condition:

$$a^{\left(1\right)}\left(u_{th}\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=u_{th}.$$

(111)

The system of equations comprising Equation (110) together with the boundary condition represented Equation (111) is the 1^st-Level Adjoint Sensitivity System (1^st-LASS) and its solution $a^{\left(1\right)}\left(u\right)$ is the 1^st-level adjoint sensitivity function.

Using Equation (101), together with the equations underlying the 1^st-LASS and 1^st-LVSS in Equation (108) reduces the latter to the following expression for the indirect-effect term:

$${\left\{\delta R\left({\phi }^0,f^0;a^{\left(1\right)}\right)\right\}}_{ind}={\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\left(\delta f_2\right)\delta \left(u\right)-\left(\delta f_1\right)\phi \left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}.$$

(112)

Adding the expression obtained in Equation (112) with the expression for the direct-effect term defined in Equation (100) yields the following expression for the total 1^st-order variation $\delta R\left({\phi }^0,f^0;v^{\left(1\right)},\delta f\right)$ of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\begin{array}{l}\delta R\left({\phi }^0,f^0;v^{\left(1\right)},\delta f\right)={\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\left(\delta f_2\right)\delta \left(u\right)-\left(\delta f_1\right)\phi \left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}\\ \hspace{0.17em}\hspace{0.17em}+{\left\{\left(\delta f_3\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\right\}}_{{\alpha }^0}\equiv {\left\{\delta R\left(\phi ,f;a^{\left(1\right)},\delta f\right)\right\}}_{{\alpha }^0}.\end{array}$$

(113)

The identity which appears in Equation (113) emphasizes the fact that “1^st-level variational sensitivity function” $v^{\left(1\right)}\left(u\right)$, which is expensive to compute, has been eliminated from the final expression of the 1^st-order total variation ${\left\{\delta R\left(\phi ,f;a^{\left(1\right)},\delta f\right)\right\}}_{{\alpha }^0}$, being replaced by the dependence on the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$, which is independent of variations $\delta f\left(\alpha \right)$ in the components of the feature function and is consequently also independent of any variations $\delta \alpha$ in the primary model parameters. Hence, the 1^st-LASS needs to be solved only once to determine the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$, which requires the same amount of computational effort as solving the original forward system for the function $\phi \left(u\right)$.

The expressions of the sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$ are given by the expressions that multiply the respective components of $f\left(\alpha \right)$ in Equation (113), namely:

$$\frac{\partial R\left(\phi ;f\right)}{\partial f_1}=-{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\hspace{0.17em};$$

(114)

$$\frac{\partial R\left(\phi ;f\right)}{\partial f_2}={\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em},$$

(115)

$$\frac{\partial R\left(\phi ;f\right)}{\partial f_3}={\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(116)

The above expressions are to be evaluated at the nominal parameter values ${\alpha }^0$ but the indication ${\left\{\hspace{0.17em}\right\}}_{{\alpha }^0}$ has been omitted for simplicity.

Solving the 1^st-LASS yields the following closed-form expression for the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$:

$$a^{\left(1\right)}\left(u\right)=H\left(u_d-u\right)f_3\left(\alpha \right)\exp \left[\left(u-u_d\right)f_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},$$

(117)

where the Heaviside functional has the usual meaning, namely: $H\left(u_d-u\right)=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}u>u_d$ and $H\left(u_d-u\right)=1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}u<u_d$.

Inserting the expression obtained in Equation (117) into Eqs. (114)‒(116) yields the following closed-form expressions for the sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\frac{\partial R\left(\phi ;f\right)}{\partial f_1}=-u_d{f}_2\left(\alpha \right)f_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em};$$

(118)

$$\frac{\partial R\left(\phi ;f\right)}{\partial f_2}=f_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},$$

(119)

$$\frac{\partial R\left(\phi ;f\right)}{\partial f_3}=f_2\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right].$$

(120)

The correctness of the expressions obtained in Eqs. (118)‒(120) can be verified by directly differentiating the closed-form expression given in Equation (91).

Alternatively, the 1^st-FASAM-L methodology could have been applied to the alternative expression, in terms of the adjoint slowing-down function, for the detector response provided in Equation (95). It can be verified that the final expressions for the response sensitivities with respect to the feature functions $f_i\left(\alpha \right)$, $i=1,2,3$, obtained by using Equation (95) as the starting point in conjunction with the adjoint slowing-down model are the same as obtained in Eqs. (118)‒(120).

The expressions of the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the primary model parameters ${\alpha }_j$, $\hspace{0.17em}j=TP\triangleq 3M+12$, as defined in Equation (88), are obtained by using the expressions obtained in Eqs. (118)‒(120) in conjunction with the chain rule of differentiation of the compound functions $f_i\left(\alpha \right)$, $i=1,2,3$. Note that the feature function $f_3\left(\alpha \right)\triangleq f_3\left(N_d,{\sigma }_d\right)$ depends only on the parameters that characterize the detector, so the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the primary model parameters $N_d$ and ${\sigma }_d$ can be readily obtained by using Equation (120), as follows:

$$\frac{\partial R\left(\phi ;f\right)}{\partial N_d}=\frac{\partial R\left(\phi ;f\right)}{\partial f_3}\frac{\partial f_3}{\partial N_d}={\sigma }_d{f}_2\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]$$

(121)

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\sigma }_d}=\frac{\partial R\left(\phi ;f\right)}{\partial f_3}\frac{\partial f_3}{\partial {\sigma }_d}=N_d{f}_2\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]$$

(122)

Similarly, the primary model parameters $\left({\lambda }_1^S,{\lambda }_2^S,N_1^S,N_2^S,F_1^S,F_2^S,{\nu }_1^S,{\nu }_2^S,W_1^S,W_2^S\right)$ that characterize the neutron source distribution only appear through the definition of the feature function $f_2\left(\alpha \right)\triangleq S_0\left(\alpha \right)/\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)$. It therefore follows that the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to these primary model parameters are obtained as follows:

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\lambda }_i^S}=\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial {\lambda }_i^S}=\frac{N_k^S{F}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(123)

$$\frac{\partial R\left(\phi ;f\right)}{\partial N_k^S}=\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial N_k^S}=\frac{{\lambda }_i^S{F}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(124)

$$\frac{\partial R\left(\phi ;f\right)}{\partial F_k^S}=\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial F_k^S}=\frac{{\lambda }_i^S{N}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(125)

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\nu }_k^S}=\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial {\nu }_k^S}=\frac{{\lambda }_i^S{N}_k^S{F}_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(126)

$$\frac{\partial R\left(\phi ;f\right)}{\partial W_k^S}=\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial W_k^S}=\frac{{\lambda }_i^S{N}_k^S{F}_k^S{\nu }_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2.$$

(127)

On the other hand, the primary model parameters $\left(N_m^{\left(1\right)},{\sigma }_{\gamma }^{\left(1\right)},{\sigma }_s^{\left(1\right)},\mathrm{...},N_m^{\left(M\right)},{\sigma }_{\gamma }^{\left(M\right)},{\sigma }_s^{\left(M\right)}\right)$ that characterize the composition of the homogenized material in which the neutrons slow-down appear through the definitions of both feature functions $f_1\left(\alpha \right)\triangleq {\Sigma }_a\left(\alpha \right)/\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)$ and $f_2\left(\alpha \right)\triangleq S_0\left(\alpha \right)/\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)$. It therefore follows that the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to these primary model parameters are obtained as follows:

$$\frac{\partial R\left(\phi ;f\right)}{\partial N_m^{\left(i\right)}}=\frac{\partial R\left(\phi ;f\right)}{\partial f_1}\frac{\partial f_1}{\partial N_m^{\left(1\right)}}+\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial N_m^{\left(1\right)}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;$$

(128)

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\sigma }_{\gamma }^{\left(i\right)}}=\frac{\partial R\left(\phi ;f\right)}{\partial f_1}\frac{\partial f_1}{\partial {\sigma }_{\gamma }^{\left(i\right)}}+\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial {\sigma }_{\gamma }^{\left(i\right)}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;$$

(129)

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\sigma }_s^{\left(i\right)}}=\frac{\partial R\left(\phi ;f\right)}{\partial f_1}\frac{\partial f_1}{\partial {\sigma }_s^{\left(i\right)}}+\frac{\partial R\left(\phi ;f\right)}{\partial f_2}\frac{\partial f_2}{\partial {\sigma }_s^{\left(i\right)}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;$$

(130)

The explicit differentiations in Eqs. (128)‒(130) are straightforward to perform, but are too lengthy to be presented here and are not material to applying the principles of the 1^st-FASAM-L methodology.

In summary, the application of the 1^st-FASAM-L to compute the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the primary model parameters ${\alpha }_j$, $\hspace{0.17em}j=1,\mathrm{...},TP\triangleq 3M+12$, requires the following computations:

• One “large-scale” computation to solve the 1^st-LASS to obtain the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$.
• Three “quadratures”, as indicated in Eqs. (114)‒(116), involving the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$ to obtain the three sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components $f_i\left(\alpha \right)$, $i=1,2,3$, of the feature function $f\left(\alpha \right)$. These computations are inexpensive.
• Chain-rule type differentiations using the definitions of the components $f_i\left(\alpha \right)$, $i=1,2,3$ of the feature function $f\left(\alpha \right)$, and the three sensitivities obtained in Eqs. (114)‒(116). These computations are inexpensive.

5.1.2. Application of the 1^st-CASAM-L

The application of the 1^st-CASAM-L methodology will yield the first-order response sensitivities directly with respect to the primary model parameters. These sensitivities will be obtained by determining the first-order Gateaux (G-) variation $\delta R\left({\phi }^0,{\alpha }^0;v^{\left(1\right)},\delta \alpha \right)$ of the response $R\left(\phi ,\alpha \right)$ as for variations $v^{\left(1\right)}\left(u\right)\triangleq \delta \phi \left(u\right)$ and $\delta \alpha$ around the phase-space point $\left({\phi }^0,{\alpha }^0\right)$, using the definition provided in Equation (87), to obtain:

$$\begin{array}{l}\delta R\left({\phi }^0,{\alpha }^0;v^{\left(1\right)},\delta \alpha \right)\triangleq {\left\{\delta R\left({\phi }^0,{\alpha }^0;\delta \alpha \right)\right\}}_{dir}+{\left\{\delta R\left({\phi }^0,{\alpha }^0;v^{\left(1\right)}\right)\right\}}_{ind}\\ \triangleq {\left\{{\left\{\frac{d}{d\epsilon }\left[\left(N_d+\epsilon \delta N_d\right)\left({\sigma }_d+\epsilon \delta {\sigma }_d\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left({\phi }^0+\epsilon v^{\left(1\right)}\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\right]\right\}}_{\epsilon =0}\right\}}_{{\alpha }^0},\end{array}$$

(131)

where the “direct-effect” term ${\left\{\delta R\left({\phi }^0,{\alpha }^0;\delta \alpha \right)\right\}}_{dir}$ and, respectively, the “indirect-effect” term ${\left\{\delta R\left({\phi }^0,{\alpha }^0;v^{\left(1\right)}\right)\right\}}_{ind}$ are defined as follows:

$${\left\{\delta R\left({\phi }^0,{\alpha }^0;\delta f_3\right)\right\}}_{dir}\triangleq {\left\{\left[\left(\delta N_d\right){\sigma }_d+\left(\delta {\sigma }_d\right)N_d\right]{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0},$$

(132)

$${\left\{\delta R\left({\phi }^0,{\alpha }^0;v^{\left(1\right)}\right)\right\}}_{ind}\triangleq {\left\{N_d{\sigma }_d{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}.$$

(133)

The “1^st-level variational sensitivity function” $v^{\left(1\right)}\left(u\right)$ is obtained as the solution of the “1^st-Level Variational Sensitivity System (1^st-LVSS)” obtained by taking the first-order G-differentials of the 1^st-LFAS defined by Eqs. (80) and (81) to obtain:

$$\begin{array}{l}{\left\{\frac{d}{d\epsilon }\left[\frac{d\left({\phi }^0+\epsilon v^{\left(1\right)}\right)}{du}+\frac{{\Sigma }_a\left({\alpha }^0+\epsilon \delta \alpha \right)}{\hspace{0.17em}\overline{\xi }\left({\alpha }^0+\epsilon \delta \alpha \right)\hspace{0.17em}{\Sigma }_t\left({\alpha }^0+\epsilon \delta \alpha \right)}\left({\phi }^0+\epsilon v^{\left(1\right)}\right)\right]\right\}}_{\epsilon =0}\\ =\delta \left(u\right){\left\{\frac{d}{d\epsilon }\frac{S_0\left({\alpha }^0+\epsilon \delta \alpha \right)}{\overline{\xi }\left({\alpha }^0+\epsilon \delta \alpha \right)\hspace{0.17em}{\Sigma }_t\left({\alpha }^0+\epsilon \delta \alpha \right)}\right\}}_{\epsilon =0},\end{array}$$

(134)

$$$$

(135)

Carrying out the differentiations with respect to $\epsilon$ in the above equations and setting $\epsilon =0$ in the resulting expressions yields the following 1^st-LVSS:

$$\begin{array}{l}\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+\frac{{\Sigma }_a\left({\alpha }^0\right)}{\hspace{0.17em}\overline{\xi }\left({\alpha }^0\right)\hspace{0.17em}{\Sigma }_t\left({\alpha }^0\right)}{v}^{\left(1\right)}\left(u\right)=\delta \left(u\right){\left\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i\right\}}_{{\alpha }^0}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\left\{\phi \left(u\right){\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i\right\}}_{{\alpha }^0},\end{array}$$

(136)

$$v^{\left(1\right)}\left(u\right)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=0.$$

(137)

To avoid solving the above the 1^st-LVSS repeatedly, for every possible variation in the primary parameters, the appearance of the function $v^{\left(1\right)}\left(u\right)$ will be eliminated for the expression of the indirect-effect term by replacing it with the solution of the 1^st-Level Adjoint Sensitivity System (1^st-LASS) which will be constructed in the Hilbert space $H_1$, as before: use Equation (107) to form the inner product of Equation (136) with a square-integrable function $a^{\left(1\right)}\left(u\right)\in H_1$ to obtain the following relation:

$$\begin{array}{l}{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}{v}^{\left(1\right)}\left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}=\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i\\ \times {{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\}}_{{\alpha }^0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\left\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}.\end{array}$$

(138)

Integration by parts of the left-side of Equation (138) yields the following relation:

$$\begin{array}{l}{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\left[\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}{v}^{\left(1\right)}\left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}={\left\{a^{\left(1\right)}\left(u_{th}\right)v^{\left(1\right)}\left(u_{th}\right)\right\}}_{{\alpha }^0}\\ -{\left\{a^{\left(1\right)}\left(0\right)v^{\left(1\right)}\left(0\right)\right\}}_{{\alpha }^0}{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)}\left[-\frac{d{a}^{\left(1\right)}\left(u\right)}{du}+\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}{a}^{\left(1\right)}\left(u\right)\right]\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}.\end{array}$$

(139)

Requiring the first term on right-side of Equation (139) to represent the indirect-effect term defined in Equation (133) yields the following relation:

$${\left\{-\frac{d{a}^{\left(1\right)}\left(u\right)}{du}+\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}{a}^{\left(1\right)}\left(u\right)\hspace{0.17em}\right\}}_{{\alpha }^0}={\left\{N_d{\sigma }_d\delta \left(u-u_d\right)\right\}}_{{\alpha }^0}.$$

(140)

Eliminate the unknown boundary-value $v^{\left(1\right)}\left(u_{th}\right)$ from Equation (139) by imposing the following boundary condition:

$$a^{\left(1\right)}\left(u_{th}\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=u_{th}.$$

(141)

The system of equations comprising Eqs. (140) and (141) is the 1^st-Level Adjoint Sensitivity System (1^st-LASS) and its solution $a^{\left(1\right)}\left(u\right)$ is the 1^st-level adjoint sensitivity function. As already shown in the general 1^st-FASAM methodology presented in Section 3, the 1^st-LASS which arises within the framework of the 1^st-CASAM-L is identical to the 1^st-LASS that arises within the 1^st-FASAM methodology, which is the reason underlying the use of the same notation for the 1^st-level adjoint sensitivity function, namely $a^{\left(1\right)}\left(u\right)$, in both cases.

Implementing the equations underlying the 1^st-LASS and 1^st-LVSS into Equation (138) and recalling the expression of the indirect-effect term provided in Equation (133) yields the following expression for the indirect-effect term:

$$\begin{array}{l}{\left\{\delta R\left({\phi }^0,{\alpha }^0;a^{\left(1\right)}\right)\right\}}_{ind}={\left\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}\\ \hspace{0.17em}\hspace{0.17em}-{\left\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}.\end{array}$$

(142)

Adding the expression obtained in Equation (142) with the expression for the direct-effect term defined in Equation (132) yields the following expression for the total 1^st-order variation $\delta R\left({\phi }^0,f^0;v^{\left(1\right)},\delta f\right)$ of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\begin{array}{l}\delta R\left({\phi }^0,{\alpha }^0;v^{\left(1\right)},\delta \alpha \right)={\left\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}\\ \hspace{0.17em}\hspace{0.17em}-{\left\{{\displaystyle \sum _{i=1}^{TP}\frac{\partial }{\partial {\alpha }_i}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\delta {\alpha }_i{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}\\ \hspace{0.17em}\hspace{0.17em}+{\left\{\left[\left(\delta N_d\right){\sigma }_d+\left(\delta {\sigma }_d\right)N_d\right]{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}\right\}}_{{\alpha }^0}\equiv {\left\{\delta R\left(\phi ,\alpha ;a^{\left(1\right)},\delta \alpha \right)\right\}}_{{\alpha }^0}.\end{array}$$

(143)

The identity which appears in Equation (143) emphasizes the fact that “1^st-level variational sensitivity function” $v^{\left(1\right)}\left(u\right)$, which is expensive to compute, has been eliminated from the final expression of the 1^st-order total variation ${\left\{\delta R\left(\phi ,\alpha ;a^{\left(1\right)},\delta \alpha \right)\right\}}_{{\alpha }^0}$, being replaced by the dependence on the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$, which is independent of any variations $\delta \alpha$ in the primary model parameters. Hence, the 1^st-LASS needs to be solved only once to determine the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$, which requires the same amount of computational effort as solving the original forward system for the function $\phi \left(u\right)$.

The expressions of the first-order sensitivities of the response $R\left[\phi \left(u\right);\alpha \right]$ with respect to the primary model parameters are the expressions that multiply the corresponding parameter variations $\delta {\alpha }_i$ in Equation (143). In particular, the (two) first-order sensitivities of the response $R\left[\phi \left(u\right);\alpha \right]$ with respect to the primary model parameters underlying the detector’s interaction cross section arise solely from the direct-effect term and have the following expressions:

$$\frac{\partial R\left(\phi ;\alpha \right)}{\partial N_d}={\sigma }_d{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du={\sigma }_d\phi \left(u_d\right)\hspace{0.17em};$$

(144)

$$\frac{\partial R\left(\phi ;\alpha \right)}{\partial {\sigma }_d}=N_d{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}=N_d\phi \left(u_d\right)\hspace{0.17em}.$$

(145)

The above expressions are to be evaluated at the nominal parameter values ${\alpha }^0$ but the indication ${\left\{\hspace{0.17em}\right\}}_{{\alpha }^0}$ has been omitted for simplicity. As expected, the above expressions are identical to the corresponding expressions obtained using the 1^st-FASAM-L, as provided in Eqs. (121) and (122), respectively.

The first-order sensitivities of the response $R\left[\phi \left(u\right);\alpha \right]$ with respect to the primary model parameters underlying the spontaneous fission source arise solely from the first term on the right-side of Equation (143) and have the following expressions in terms of the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$:

$$\frac{\partial R\left(\phi ;\alpha \right)}{\partial {\lambda }_k^S}=\frac{N_k^S{F}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(146)

$$\frac{\partial R\left(\phi ;f\right)}{\partial N_k^S}=\frac{{\lambda }_i^S{F}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(147)

$$\frac{\partial R\left(\phi ;f\right)}{\partial F_k^S}=\frac{{\lambda }_i^S{N}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(148)

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\nu }_k^S}=\frac{{\lambda }_i^S{N}_k^S{F}_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(149)

$$\frac{\partial R\left(\phi ;f\right)}{\partial W_k^S}=\frac{{\lambda }_i^S{N}_k^S{F}_k^S{\nu }_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}k=1,2.$$

(150)

As expected, the above expressions are identical to the corresponding expressions obtained using the 1^st-FASAM-L, as provided in Eqs. (123)‒(127), respectively.

The first-order sensitivities of the response $R\left[\phi \left(u\right);\alpha \right]$ with respect to the primary model parameters $\left(N_m^{\left(1\right)},{\sigma }_{\gamma }^{\left(1\right)},{\sigma }_s^{\left(1\right)},\mathrm{...},N_m^{\left(M\right)},{\sigma }_{\gamma }^{\left(M\right)},{\sigma }_s^{\left(M\right)}\right)$ that characterize the composition of the homogenized material in which the neutrons slow-down arise from both the first and the second terms on the right-side of Equation (143) and have the following expressions in terms of the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$:

$$\begin{array}{l}\frac{\partial R\left(\phi ;f\right)}{\partial N_m^{\left(i\right)}}=\left[{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right]\frac{\partial }{\partial N_m^{\left(i\right)}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}-\left[{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right]\frac{\partial }{\partial N_m^{\left(i\right)}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;\end{array}$$

(151)

$$\begin{array}{l}\frac{\partial R\left(\phi ;f\right)}{\partial {\sigma }_{\gamma }^{\left(i\right)}}=\left[{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right]\frac{\partial }{\partial {\sigma }_{\gamma }^{\left(i\right)}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}-\left[{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right]\frac{\partial }{\partial {\sigma }_{\gamma }^{\left(i\right)}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;\end{array}$$

(152)

$$\begin{array}{l}\frac{\partial R\left(\phi ;f\right)}{\partial {\sigma }_s^{\left(i\right)}}=\left[{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right]\frac{\partial }{\partial {\sigma }_s^{\left(i\right)}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}-\left[{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right]\frac{\partial }{\partial {\sigma }_s^{\left(i\right)}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M.\end{array}$$

(153)

As expected, the above expressions are identical to the corresponding expressions obtained using the 1^st-FASAM-L, as provided in Eqs. (128)‒(130), respectively.

In summary, the application of the 1^st-CASAM-L to compute the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the primary model parameters ${\alpha }_j$, $\hspace{0.17em}j=1,\mathrm{...},TP\triangleq 3M+12$, requires the following computations:

• One “large-scale” computation to solve the 1^st-LASS to obtain the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$. As has been already remarked, this 1^st-LASS is exactly the same as the 1^st-LASS needed within the 1^st-FASAM-L methodology for computing the first-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components $f_i\left(\alpha \right)$, $i=1,2,3$, of the feature function $f\left(\alpha \right)$.
• A total of $TP\triangleq 3M+12$ “quadratures” involving the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$ to obtain numerically the $TP\triangleq 3M+12$ sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the primary model parameters ${\alpha }_j$, $\hspace{0.17em}j=1,\mathrm{...},TP\triangleq 3M+12$. These numerical computations are inexpensive by comparison to solving the 1^st-LASS but are more expensive than performing “chain-rule”-type differentiation “on paper,” as performed if applying the 1^st-FASAM-L. Hence, the 1^st-FASAM-L methodology enjoys a slight computational advantage over the 1^st-CASAM-L methodology.

5.2. Second-Order Adjoint Sensitivity Analysis: 2^nd-FASAM-L Versus 2^nd-CASAM-L

In this Subsection, the computation of the first-order sensitivities of the response $R\left(\alpha \right)$ with respect to the primary model parameters will first be demonstrated by using the 2^nd-FASAM-L. Subsequently, the same first-order sensitivities will be obtained by using the 2^nd-CASAM-L and the two alternative paths will be compared to each other, showing that the same expressions are obtained for the respective sensitivities, as expected. Both the 2^nd-FASAM-L and the 2^nd-CASAM-L methodologies obtain the second-order sensitivities by considering the first-order G-differential of each of the first-order sensitivities. Therefore, the 2^nd-FASAM-L methodology will provide significant computational advantages by comparison to the 2^nd-CASAM-L methodology the since it will require at most 3 large-scale computations, i.e., the same number of large-scale computations as the number of components $f_i\left(\alpha \right)$, $i=1,2,3$, of the feature function $f\left(\alpha \right)$. In contradistinction, the 2^nd-CASAM-L methodology will require one large-scale (adjoint) computation for each primary model parameter ${\alpha }_j$, $\hspace{0.17em}j=1,\mathrm{...},TP\triangleq 3M+12$, amounting to a total of number of $TP\triangleq 3M+12$ large-scale computations.

5.2.1. Application of the 2^nd-FASAM-L

As has been shown in Section 4, the 2^nd-FASAM-L methodology generically determines the 2^nd-order sensitivities ${\partial }^{2}R\left[u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]/\partial f_{j_2}\partial f_{j_1}$ of the response with respect to the components of the “feature” function $f\left(\alpha \right)$ by conceptually considering that the first-order sensitivities $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]\triangleq \partial R\left[u^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]/\partial f_{j_1}$, are “model responses.” Consequently, the 2^nd-order sensitivities are obtained as the “1^st-order sensitivities of the 1^st-order sensitivities” by applying the concepts underlying 1^st-FASAM to each 1^st-order sensitivity $R^{\left(1\right)}\left[j_1;u^{\left(1\right)}\left(2;x\right);a^{\left(1\right)}\left(2;x\right);f\left(\alpha \right)\right]$, $j_1=1,\mathrm{...},TF$.

5.2.1.1. Second-Order Sensitivities Stemming From

$\partial R\left(\phi ;f\right)/\partial f_1$The above principles will be applied to the first-order sensitivity $R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\triangleq \partial R\left(\phi ;f\right)/\partial f_1$ expressed by Equation (114) to obtain the second-order sensitivities of the form ${\partial }^{2}R\left(\phi ;f\right)/\partial f_j\partial f_1$, $j=1,2,3$. The argument “1” in the notation $R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]$ indicates that this sensitivity is with respect to the first component, namely $f_1\left(\alpha \right)$, of the feature function $f\left(\alpha \right)$, while also depending on the functions $\phi \left(u\right)$ and $a^{\left(1\right)}\left(u\right)$. These functions are the solutions of the “2^nd-Level Forward/Adjoint System (2^nd-LFAS)” which is obtained by concatenating the original 1^st-Level Forward/Adjoint System (1st-LFAS) with the 1^st-Level Adjoint Sensitivity System (1^st-LASS), cf. Eqs. (98), (81), (110) and (111), as reproduced below:

$$\frac{d\phi \left(u\right)}{du}+f_1\phi \left(u\right)=f_2\delta \left(u\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<u\le u_{th}\hspace{0.17em}\hspace{0.17em};$$

(154)

$$-\frac{d{a}^{\left(1\right)}\left(u\right)}{du}+f_1{a}^{\left(1\right)}\left(u\right)=f_3\delta \left(u-u_d\right);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<u\le u_{th}\hspace{0.17em};$$

(155)

$$\phi \left(0\right)=0;\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}^{\left(1\right)}\left(u_{th}\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=u_{th}.$$

(156)

The first-order G-differential of $R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]$ is obtained from Equation (114), by definition, as follows:

$$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f;\delta f\right]\right\}}_{{\alpha }^0}\\ \triangleq -{\left\{\frac{d}{d\epsilon }{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[a^{\left(1\right)}\left(u\right)+\epsilon \delta a^{\left(1\right)}\left(u\right)\right]\left[\phi \left(u\right)+\epsilon v^{\left(1\right)}\left(u\right)\right]du}\right\}}_{{\alpha }^0,\epsilon =0}\\ =-{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[v^{\left(1\right)}\left(u\right)a^{\left(1\right)}\left(u\right)+\delta a^{\left(1\right)}\left(u\right)\phi \left(u\right)\right]}\hspace{0.17em}du\right\}}_{{\alpha }^0}\equiv {\displaystyle \sum _{j=1}^3\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_j\partial f_1}}\left(\delta f_j\hspace{0.17em}\right).\end{array}$$

(157)

Note that the first-order G-differential ${\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$ consists solely of the “indirect-effect term”; there is no “direct-effect” term since there is no explicit dependence on variations $\delta f\left(\alpha \right)$.

The variational functions $v^{\left(1\right)}\left(u\right)$ and $\delta a^{\left(1\right)}\left(u\right)$ are the solutions of the system of equations obtained by taking the first-G-differential of the 2^nd-LFAS. Applying the definition of the first G-differential to the equations underlying the 2^nd-LFAS yields the following 2^nd-Level Variational Sensitivity System (2^nd-LVSS)” for the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2;u\right)\triangleq {\left[v^{\left(1\right)}\left(u\right),\delta a^{\left(1\right)}\left(u\right)\right]}^{†}$:

$${\left\{\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+f_1{v}^{\left(1\right)}\left(u\right)\right\}}_{{\alpha }^0}={\left\{\left(\delta f_2\right)\delta \left(u\right)-\left(\delta f_1\right)\phi \left(u\right)\right\}}_{{\alpha }^0}\hspace{0.17em},$$

(158)

$${\left\{-\frac{d\left[\delta a^{\left(1\right)}\left(u\right)\right]}{du}+f_1\left[\delta a^{\left(1\right)}\left(u\right)\right]\right\}}_{{\alpha }^0}={\left\{\left(\delta f_3\right)\delta \left(u-u_d\right)-\left(\delta f_1\right)a^{\left(1\right)}\left(u\right)\right\}}_{{\alpha }^0};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<u\le u_{th}\hspace{0.17em};$$

(159)

$$v^{\left(1\right)}\left(u\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=0\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\delta a^{\left(1\right)}\left(u_{th}\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}u=u_{th}.$$

(160)

The above 2^nd-LVSS would need to be solved repeatedly, for every possible variation in the feature functions $f_i\left(\alpha \right)$, $i=1,2,3$. This need is circumvented by deriving an alternative expression for ${\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$, in which the 2^nd-level variational function $v^{\left(2\right)}\left(2;u\right)\triangleq {\left[v^{\left(1\right)}\left(u\right),\delta a^{\left(1\right)}\left(u\right)\right]}^{†}$ is replaced by a 2^nd-level adjoint sensitivity function which will be independent of variations in the feature functions $f_i\left(\alpha \right)$, $i=1,2,3$. This 2^nd-level adjoint sensitivity function will be the solution of a 2^nd-Level Adjoint Sensitivity System (2nd-LASS) to be constructed below by following the steps generally outlined in Section 5, in a Hilbert space, denoted as $H_2$, which is endowed with the following inner product, denoted as ${\langle {\chi }^{\left(2\right)}\left(2;u\right),{\theta }^{\left(2\right)}\left(2;u\right)\rangle }_2$, between two elements, ${\chi }^{\left(2\right)}\left(2;u\right)\triangleq {\left[{\chi }_1^{\left(2\right)}\left(u\right),{\chi }_2^{\left(2\right)}\left(u\right)\right]}^{†}\in H_2$ and ${\theta }^{\left(2\right)}\left(2;u\right)\triangleq {\left[{\theta }_1^{\left(2\right)}\left(u\right),{\theta }_2^{\left(2\right)}\left(u\right)\right]}^{†}\in H_2$:

$${\langle {\chi }^{\left(2\right)}\left(2;u\right),{\theta }^{\left(2\right)}\left(2;u\right)\rangle }_2\triangleq {\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[{\chi }_1^{\left(2\right)}\left(u\right){\theta }_1^{\left(2\right)}\left(u\right)+{\chi }_2^{\left(2\right)}\left(u\right){\theta }_2^{\left(2\right)}\left(u\right)\right]}\hspace{0.17em}du\hspace{0.17em}.$$

(161)

Using Equation(161) to form the inner product in the Hilbert space $H_2$ of the 2^nd-LVSS, cf. Eqs. (158) and (159), with a yet undefined function $a^{\left(2\right)}\left(2;1;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;1;u\right),a_2^{\left(2\right)}\left(2;1;u\right)\right]}^{†}\in H_2$ yields the following relation:

$$\begin{array}{l}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\left[\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+f_1{v}^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}du+\hspace{0.17em}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)\left[-\frac{d}{du}\delta a^{\left(1\right)}\left(u\right)+f_1\delta a^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}du\hspace{0.17em}\\ ={\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\left[\left(\delta f_2\right)\delta \left(u\right)-\left(\delta f_1\right)\phi \left(u\right)\right]}\hspace{0.17em}du\hspace{0.17em}\\ +{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)\left[\left(\delta f_3\right)\delta \left(u-u_d\right)-\left(\delta f_1\right)a^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}du\hspace{0.17em}.\end{array}$$

(162)

The above relation holds for the nominal parameter values, but the notation ${\left\{\right\}}_{{\alpha }^0}\hspace{0.17em}$ has been omitted, for simplicity.

Using the definition of the adjoint operator in the Hilbert space $H_2$, which amounts to integration by parts, recast the left-side of Equation (162) into the form below:

$$\begin{array}{l}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\left[\frac{d{v}^{\left(1\right)}\left(u\right)}{du}+f_1{v}^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}du+\hspace{0.17em}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)\left[-\frac{d}{du}\delta a^{\left(1\right)}\left(u\right)+f_1\delta a^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}du\hspace{0.17em}\\ ={\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)\left[-\frac{d{a}_1^{\left(2\right)}\left(2;1;u\right)}{du}+f_1{a}_1^{\left(2\right)}\left(2;1;u\right)\right]}\hspace{0.17em}du+{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\delta a^{\left(1\right)}\left(u\right)\left[\frac{d{a}_2^{\left(2\right)}\left(2;1;u\right)}{du}+f_1{a}_2^{\left(2\right)}\left(2;1;u\right)\right]}\hspace{0.17em}du\hspace{0.17em}\\ +a_1^{\left(2\right)}\left(2;1;u_{th}\right)v^{\left(1\right)}\left(u_{th}\right)-a_1^{\left(2\right)}\left(2;1;0\right)v^{\left(1\right)}\left(0\right)-\hspace{0.17em}{a}_2^{\left(2\right)}\left(2;1;u_{th}\right)\delta a^{\left(1\right)}\left(u_{th}\right)+a_2^{\left(2\right)}\left(2;1;0\right)\delta a^{\left(1\right)}\left(0\right).\end{array}$$

(163)

The first two terms on right-side of Equation (163) are now required to represent the G-differential ${\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$ defined in Equation (157), which yields the following relations:

$$\left(\begin{array}{cc}-d/du+f_1& 0\\ 0& d/du+f_1\end{array}\right)\left(\begin{array}{c}{a}_1^{\left(2\right)}\left(2;1;u\right)\\ a_2^{\left(2\right)}\left(2;1;u\right)\end{array}\right)=\left(\begin{array}{c}-a^{\left(1\right)}\left(u\right)\\ -\phi \left(u\right)\end{array}\right).$$

(164)

The definition of the vector $a^{\left(2\right)}\left(2;1;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;1;u\right),a_2^{\left(2\right)}\left(2;1;u\right)\right]}^{†}$ will now be completed by selecting boundary conditions so as to eliminate the unknown values $v^{\left(1\right)}\left(u_{th}\right)$ and $\delta a^{\left(1\right)}\left(0\right)$ in Equation (163). This is accomplished by imposing the following boundary conditions:

$$a_1^{\left(2\right)}\left(2;1;u_{th}\right)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2^{\left(2\right)}\left(2;1;0\right)=0.$$

(165)

The system of equations represented by Eqs. (164) and (165) constitute the 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;1;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;1;u\right),a_2^{\left(2\right)}\left(2;1;u\right)\right]}^{†}$. It is important to note that the 2^nd-LASS is independent of any variations, $\delta f$, in the components of the feature function and, hence, is independent of any parameter variations, $\delta \alpha$, as well.

The equations underlying the 2^nd-LASS, together with the equations underlying the 2nd-LVSS, are now employed in Equation(162), in conjunction with Equation (163), to obtain the following expression for the G-differential ${\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$ in terms of the 2nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;1;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;1;u\right),a_2^{\left(2\right)}\left(2;1;u\right)\right]}^{†}$:

$$\begin{array}{l}\hspace{0.17em}{\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}=\hspace{0.17em}\left(\delta f_2\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\delta \left(u\right)}\hspace{0.17em}du\\ \hspace{0.17em}-\left(\delta f_1\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\phi \left(u\right)}\hspace{0.17em}du\hspace{0.17em}+\left(\delta f_3\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\\ -\left(\delta f_1\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)a^{\left(1\right)}\left(u\right)}\hspace{0.17em}du\hspace{0.17em}\equiv {\displaystyle \sum _{j=1}^3\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_j\partial f_1}}\left(\delta f_j\hspace{0.17em}\right)\\ \equiv {\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);a^{\left(2\right)}\left(2;1;u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}.\end{array}$$

(166)

As the last equality (identity) in Equation (166) indicates, the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2;u\right)\triangleq {\left[v^{\left(1\right)}\left(u\right),\delta a^{\left(1\right)}\left(u\right)\right]}^{†}$ has been eliminated from the final expression of the G-differential ${\left\{\delta R^{\left(1\right)}\left[1;\phi \left(u\right);a^{\left(1\right)}\left(u\right);v^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$, having been replaced by the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;1;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;1;u\right),a_2^{\left(2\right)}\left(2;1;u\right)\right]}^{†}$. Identifying in Equation (166) the expressions that multiply the variations $\delta f_i$, $i=1,2,3$, yields the following expressions for the second-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_1\partial f_1}=-{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\phi \left(u\right)}\hspace{0.17em}du\hspace{0.17em}-{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)a^{\left(1\right)}\left(u\right)}\hspace{0.17em}du\hspace{0.17em};$$

(167)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_2\partial f_1}=\hspace{0.17em}{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;1;u\right)\delta \left(u\right)}\hspace{0.17em}du\hspace{0.17em};$$

(168)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_3\partial f_1}={\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;1;u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(169)

The 2^nd-LASS can be solved to obtain the following closed-form expressions for the components of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;1;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;1;u\right),a_2^{\left(2\right)}\left(2;1;u\right)\right]}^{†}$;

$$a_1^{\left(2\right)}\left(2;1;u\right)=-f_3\left(\alpha \right)\left(u_d-u\right)H\left(u_d-u\right)\exp \left[\left(u-u_d\right)f_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},$$

(170)

$$a_2^{\left(2\right)}\left(2;1;u\right)=-u{f}_2\left(\alpha \right)\exp \left[-u{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em},$$

(171)

Inserting the expressions obtained in Eqs. (170) and (171) into Eqs. (167)‒(169) and performing the respective integrations yields the following expressions for the respective second-order sensitivities:

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_1\partial f_1}={\left(u_d\right)}^2{f}_2\left(\alpha \right)f_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em};$$

(172)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_2\partial f_1}=-u_d{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em};$$

(173)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_3\partial f_1}=-u_d{f}_2\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em};$$

(174)

The correctness of the expressions obtained in Eqs. (172)‒(174) can be verified by directly differentiating the closed-form expressions provided in Eqs. (118)‒(120).

5.2.1.2. Second-Order Sensitivities Stemming From $\partial R\left(\phi ;f\right)/\partial f_2$

Applying the procedure used in the previous Subsection to the first-order sensitivity $R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\triangleq \partial R\left(\phi ;f\right)/\partial f_2$ expressed by Equation (115) will provide the second-order sensitivities of the form ${\partial }^{2}R\left(\phi ;f\right)/\partial f_j\partial f_2$, $j=1,2,3$. The argument “2” in the notation $R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]$ indicates that this sensitivity is with respect to the second component, namely $f_2\left(\alpha \right)$, of the feature function $f\left(\alpha \right)$. Remarkably, this sensitivity does not depend on the forward function $\phi \left(u\right)$ but only depends on the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$. As previously discussed, these functions are the solutions of the “2^nd-Level Forward/Adjoint System (2^nd-LFAS).”

The first-order G-differential of $R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]$ is obtained by applying its definition to Equation (115), as follows:

$${\left\{\delta R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}={\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\delta a^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}\equiv {\displaystyle \sum _{j=1}^3\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_j\partial f_2}}\left(\delta f_j\right)\hspace{0.17em}..$$

(175)

As indicated in Equation (175), the first-order G-differential ${\left\{\delta R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);\delta a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$ consists solely of the indirect-effect term which depends on the 1^st-level variational function $\delta a^{\left(1\right)}\left(u\right)$. As before, the need for computing this variational function is circumvented by constructing a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for a 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;2;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;2;u\right),a_2^{\left(2\right)}\left(2;2;u\right)\right]}^{†}$, by implementing the same steps as outlined above for obtaining the 2^nd-order sensitivities that stem from the first-order sensitivity $\partial R\left(\phi ;f\right)/\partial f_2$, namely Eqs. (162)‒(165). These steps will not be repeated here in detail; they lead to the following 2^nd-LASS for 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;2;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;2;u\right),a_2^{\left(2\right)}\left(2;2;u\right)\right]}^{†}$:

(176)

$$a_1^{\left(2\right)}\left(2;2;u_{th}\right)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2^{\left(2\right)}\left(2;2;0\right)=0.$$

(177)

Solving the 2^nd-LASS defined by Eqs. (176) and (177) yields the following closed-form expressions for the components of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;2;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;2;u\right),a_2^{\left(2\right)}\left(2;2;u\right)\right]}^{†}$:

$$a_1^{\left(2\right)}\left(2;2;u\right)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2^{\left(2\right)}\left(2;2;u\right)=H\left(u\right)\exp \left[-u{f}_1\left(\alpha \right)\right]\hspace{0.17em}$$

(178)

.

The alternative expression of the G-differential ${\left\{\delta R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);a^{\left(2\right)}\left(2;2;u\right);\delta f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}$ in terms of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;2;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;2;u\right),a_2^{\left(2\right)}\left(2;2;u\right)\right]}^{†}$ has the following form (which is obtained by implementing the same steps as those leading to Equation (166), above):

$$\begin{array}{l}{\left\{\delta R^{\left(1\right)}\left[2;a^{\left(1\right)}\left(u\right);a^{\left(2\right)}\left(2;2;u\right);\delta f\left(\alpha \right)\right]\right\}}_{{\alpha }^0}=\left(\delta f_3\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;2;u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\\ -\left(\delta f_1\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;2;u\right)a^{\left(1\right)}\left(u\right)}\hspace{0.17em}du\hspace{0.17em}\equiv {\displaystyle \sum _{j=1}^3\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_j\partial f_2}}\left(\delta f_j\right)\hspace{0.17em}.\end{array}$$

(179)

Identifying in Equation (179) the expressions that multiply the variations $\delta f_i$, $i=1,2,3$, yields the following second-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_1\partial f_2}=-{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;2;u\right)a^{\left(1\right)}\left(u\right)}\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em};$$

(180)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_2\partial f_2}=0\hspace{0.17em};$$

(181)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_3\partial f_2}={\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_2^{\left(2\right)}\left(2;2;u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em}.$$

(182)

Inserting the expression obtained for $a_2^{\left(2\right)}\left(2;2;u\right)$ in Equation (178) into Eqs. (180) and (182), and performing the respective integrations yields the following expressions for the respective second-order sensitivities:

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_1\partial f_2}=-u_d{f}_3\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}\hspace{0.17em};$$

(183)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_3\partial f_2}=\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em};$$

(184)

The correctness of the expressions obtained in Eqs. (183) and (184) can be verified by directly differentiating accordingly the closed-form expressions given in Eqs. (118)‒(120).

5.2.1.3. Second-Order Sensitivities Stemming From $\partial R\left(\phi ;f\right)/\partial f_3$

Applying the above principles to the first-order sensitivity $R^{\left(1\right)}\left[3;\phi \left(u\right);a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]\triangleq \partial R\left(\phi ;f\right)/\partial f_3$ obtained in Equation (116) will provide the second-order sensitivities of the form ${\partial }^{2}R\left(\phi ;f\right)/\partial f_j\partial f_3$, $j=1,2,3$. The argument “3” in the notation $R^{\left(1\right)}\left[3;\phi \left(u\right);a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]$ indicates that this sensitivity is with respect to the third component, namely $f_3\left(\alpha \right)$, of the feature function $f\left(\alpha \right)$. Notably, this sensitivity depends on the forward function $\phi \left(u\right)$ but does not depend on the 1st-level adjoint sensitivity function $a^{\left(1\right)}\left(u\right)$.

The first-order G-differential of $R^{\left(1\right)}\left[3;\phi \left(u\right);a^{\left(1\right)}\left(u\right);f\left(\alpha \right)\right]$ is obtained, by definition, as follows:

$${\left\{\delta R^{\left(1\right)}\left[3;\phi \left(u\right);v^{\left(1\right)}\left(u\right);f;\delta f\right]\right\}}_{{\alpha }^0}={\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\right\}}_{{\alpha }^0}\equiv {\displaystyle \sum _{j=1}^3\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_j\partial f_3}}\left(\delta f_j\right)\hspace{0.17em}.$$

(185)

Note that the first-order G-differential ${\left\{\delta R^{\left(1\right)}\left[3;\phi \left(u\right);v^{\left(1\right)}\left(u\right);f;\delta f\right]\right\}}_{{\alpha }^0}$ consists solely of the indirect-effect term. As before, the need for computing the variational function $v^{\left(1\right)}\left(u\right)$ is circumvented by constructing a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for a 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;3;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;3;u\right),a_2^{\left(2\right)}\left(2;3;u\right)\right]}^{†}$, by implementing the same steps were used for obtaining the previous 2^nd-order sensitivities. These steps lead to the following 2^nd-LASS for 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;3;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;3;u\right),a_2^{\left(2\right)}\left(2;3;u\right)\right]}^{†}$:

$$\left(\begin{array}{cc}-d/du+f_1& 0\\ 0& d/du+f_1\end{array}\right)\left(\begin{array}{c}{a}_1^{\left(2\right)}\left(2;3;u\right)\\ a_2^{\left(2\right)}\left(2;3;u\right)\end{array}\right)=\left(\begin{array}{c}\delta \left(u-u_d\right)\\ 0\end{array}\right).$$

(186)

$$a_1^{\left(2\right)}\left(2;2;u_{th}\right)=0;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2^{\left(2\right)}\left(2;2;0\right)=0.$$

(187)

Solving the 2^nd-LASS defined by Eqs. (186) and (187) yields the following closed-form expressions for the components of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;3;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;3;u\right),a_2^{\left(2\right)}\left(2;3;u\right)\right]}^{†}$:

$$a_1^{\left(2\right)}\left(2;3;u\right)=H\left(u_d-u\right)\exp \left[\left(u-u_d\right)f_1\left(\alpha \right)\right];\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_2^{\left(2\right)}\left(2;3;u\right)=\hspace{0.17em}0$$

(188)

In terms of the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(2;3;u\right)\triangleq {\left[a_1^{\left(2\right)}\left(2;3;u\right),a_2^{\left(2\right)}\left(2;3;u\right)\right]}^{†}$ the alternative expression of the G-differential ${\left\{\delta R^{\left(1\right)}\left[3;\phi \left(u\right);a^{\left(2\right)}\left(2;3;u\right);f;\delta f\right]\right\}}_{{\alpha }^0}$ has the following form (which is obtained by implementing the same steps as those leading to Equation (166), above):

$$\begin{array}{l}\hspace{0.17em}{\left\{\delta R^{\left(1\right)}\left[3;\phi \left(u\right);a^{\left(2\right)}\left(2;3;u\right);f;\delta f\right]\right\}}_{{\alpha }^0}=\hspace{0.17em}\left(\delta f_2\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;3;u\right)\delta \left(u\right)}\hspace{0.17em}du\\ \hspace{0.17em}-\left(\delta f_1\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;3;u\right)\phi \left(u\right)}\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em}\equiv {\displaystyle \sum _{j=1}^3\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_j\partial f_3}}\left(\delta f_j\hspace{0.17em}\right).\end{array}$$

(189)

Identifying in Equation (189) the expressions that multiply the variations $\delta f_i$, $i=1,2,3$, yields the following second-order sensitivities of the response $R\left[\phi \left(u\right);f\left(\alpha \right)\right]$ with respect to the components of the feature function $f\left(\alpha \right)$:

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_1\partial f_3}=-{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;3;u\right)\phi \left(u\right)}\hspace{0.17em}du\hspace{0.17em}\hspace{0.17em};$$

(190)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_2\partial f_3}={\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}_1^{\left(2\right)}\left(2;3;u\right)\delta \left(u\right)}\hspace{0.17em}du\hspace{0.17em};$$

(191)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_3\partial f_3}=0\hspace{0.17em}.$$

(192)

Inserting the expression obtained for $a_1^{\left(2\right)}\left(2;3;u\right)$ in Equation (188) into Eqs. (190) and (191), and performing the respective integrations yields the following expressions for the respective second-order sensitivities:

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_1\partial f_3}=-u_d{f}_2\left(\alpha \right)\exp \left[-u_d{f}_1\left(\alpha \right)\right];$$

(193)

$$\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial f_2\partial f_3}=\exp \left[-u_d{f}_1\left(\alpha \right)\right]\hspace{0.17em}$$

(194)

The correctness of the expressions obtained in Eqs. (193) and (194) can be verified by directly differentiating the closed-form expressions given in Eqs. (118)‒(120).

Summarizing the results obtained in Subsection 5.2.1 leads to the following conclusions:

• The second-order sensitivities ${\partial }^{2}R\left(\phi ;f\right)/\partial f_i\partial f_j$, $i,j=1,2,3$, of the model response with respect to the three features components $f_i\left(\alpha \right)$, $i=1,2,3$, of the feature function $f\left(\alpha \right)$ are obtained by performing 3 “large-scale” computations to solve the 3 corresponding 2^nd-LASS which all have the same left-side but have differing sources on their right-sides. The source-term for each of these 2^nd-LASS corresponds to one of the 3 first-order sensitivities. Thus, computing the second-order sensitivities ${\partial }^{2}R\left(\phi ;f\right)/\partial f_i\partial f_j$ requires as many “large-scale” computations as there are non-zero first-order sensitivities, i.e., at most as many “large-scale” computations as there are components $f_i\left(\alpha \right)$, $i=1,2,3$, of the feature function $f\left(\alpha \right)$.
• The mixed second-order sensitivities ${\partial }^{2}R\left(\phi ;f\right)/\partial f_i\partial f_j$, $i\ne j=1,2,3$, are computed twice, involving distinct 2^nd-level adjoint sensitivity functions. Therefore, the symmetry property ${\partial }^{2}R\left(\phi ;f\right)/\partial f_i\partial f_j={\partial }^{2}R\left(\phi ;f\right)/\partial f_j\partial f_i$ provides an intrinsic mechanism for verifying the accuracy of the computations of the respective 2^nd-level adjoint sensitivity functions.
• The unmixed second order sensitivities ${\partial }^{2}R\left(\phi ;f\right)/\partial f_i\partial f_i$$i=1,2,3$, are computed just once.

5.2.2. Application of the 2^nd-CASAM-L

The principles underlying the application of the 2^nd-CASAM-L methodology are the same as those underlying the 2^nd-FASAM-L methodology: both methodologies obtain the second-order sensitivities by considering the first-order G-differential of each of the first-order sensitivities. As has been shown in the foregoing, the 2^nd-FASAM-L methodology requires at most 3 large-scale computations (i.e., the same number of large-scale computations as the number of components $f_i$, $i=1,2,3$, of the feature function $f$) for solving the three 2^nd-Level Adjoint Sensitivity Systems that arise by considering the three first-order sensitivities of the detector response with respect to the three components of the feature function $f\left(\alpha \right)$. In contradistinction, the 2^nd-CASAM-L methodology requires one large-scale (adjoint) computation for each primary model parameter ${\alpha }_j$, $\hspace{0.17em}j=1,\mathrm{...},TP\triangleq 3M+12$, amounting to a total of number of $TP\triangleq 3M+12$ large-scale computations. The specific computations are described below.

5.2.2.1. Second-Order Sensitivities Stemming from the First-Order Sensitivities with Respect to the Medium’s Material Properties

The expressions of the first-order sensitivities of the detector response with respect to the material properties (i.e., microscopic cross sections and atomic number densities) of the medium in which the neutrons are slowing-down (i.e., losing energy or, equivalently, gaining lethargy) are provided in Eqs. (151)‒(153). These expressions have the following generic form:

$$\frac{\partial R\left(\phi ;f\right)}{\partial a_j^{\left(i\right)}}=g_j^{\left(i\right)}\left(\alpha \right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du-h_j^{\left(i\right)}\left(\alpha \right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du;\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;\hspace{0.17em}j=1,2,3,$$

(195)

where $a_1^{\left(i\right)}\triangleq N_m^{\left(i\right)}$, $a_2^{\left(i\right)}\triangleq {\sigma }_{\gamma }^{\left(i\right)}$, $a_3^{\left(i\right)}\triangleq {\sigma }_s^{\left(i\right)}$, and

$$g_1^{\left(i\right)}\left(\alpha \right)\triangleq \frac{\partial }{\partial N_m^{\left(i\right)}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{h}_1^{\left(i\right)}\left(\alpha \right)\triangleq \frac{\partial }{\partial N_m^{\left(i\right)}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em};$$

(196)

$$g_2^{\left(i\right)}\left(\alpha \right)\triangleq \frac{\partial }{\partial {\sigma }_{\gamma }^{\left(i\right)}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{h}_2^{\left(i\right)}\left(\alpha \right)\triangleq \frac{\partial }{\partial {\sigma }_{\gamma }^{\left(i\right)}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em};$$

(197)

$$g_3^{\left(i\right)}\left(\alpha \right)\triangleq \frac{\partial }{\partial {\sigma }_s^{\left(i\right)}}\left[\frac{S_0\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{h}_3^{\left(i\right)}\left(\alpha \right)\triangleq \frac{\partial }{\partial {\sigma }_s^{\left(i\right)}}\left[\frac{{\Sigma }_a\left(\alpha \right)}{\hspace{0.17em}\overline{\xi }\left(\alpha \right)\hspace{0.17em}{\Sigma }_t\left(\alpha \right)}\right]\hspace{0.17em}.$$

(198)

The second-order sensitivities stemming from the first-order sensitivities represented by Equation (195) are obtained from the first G-differential of this equation, which has the following expression, by definition, for each $i=1,\mathrm{...},M;\hspace{0.17em}j=1,2,3$:

$$\begin{array}{l}\delta \left[\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}\right]\triangleq {\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\left\{\frac{d}{d\epsilon }\left[g_j^{\left(i\right)}\left({\alpha }^0+\epsilon \delta \alpha \right)\right]\right\}}_{{\alpha }^0,\epsilon =0}\\ +{\left\{g_j^{\left(i\right)}\left({\alpha }^0\right)\frac{d}{d\epsilon }{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[a^{\left(1\right)}\left(u\right)+\epsilon \delta a^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0,\epsilon =0}\\ -{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\left\{\frac{d}{d\epsilon }\left[h_j^{\left(i\right)}\left({\alpha }^0+\epsilon \delta \alpha \right)\right]\right\}}_{{\alpha }^0,\epsilon =0}\\ -{\left\{h_j^{\left(i\right)}\left({\alpha }^0\right)\frac{d}{d\epsilon }{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[a^{\left(1\right)}\left(u\right)+\epsilon \delta a^{\left(1\right)}\left(u\right)\right]\left[\phi \left(u\right)+\epsilon v^{\left(1\right)}\left(u\right)\right]du}\right\}}_{{\alpha }^0,\epsilon =0}\\ ={\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}\right]\right\}}_{dir}+{\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}\right]\right\}}_{ind}={\displaystyle \sum _{n=1}^{TP}{\left\{\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial {\alpha }_n\partial a_j^{\left(i\right)}}\delta {\alpha }_n\right\}}_{{\alpha }^0}},\end{array}$$

(199)

where:

$$\begin{array}{l}{\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}\right]\right\}}_{dir}\triangleq {\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\displaystyle \sum _{n=1}^{TP}{\left\{\frac{\partial g_j^{\left(i\right)}\left(\alpha \right)}{\partial {\alpha }_n}\right\}}_{{\alpha }^0}\delta {\alpha }_n}\\ -{\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\phi \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\displaystyle \sum _{n=1}^{TP}{\left\{\frac{\partial h_j^{\left(i\right)}\left(\alpha \right)}{\partial {\alpha }_n}\right\}}_{{\alpha }^0}\delta {\alpha }_n},\end{array}$$

(200)

$$\begin{array}{l}{\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}\right]\right\}}_{ind}\triangleq {\left\{g_j^{\left(i\right)}\left({\alpha }^0\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\delta a^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}\\ -{\left\{h_j^{\left(i\right)}\left({\alpha }^0\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[a^{\left(1\right)}\left(u\right)v^{\left(1\right)}\left(u\right)+\phi \left(u\right)\delta a^{\left(1\right)}\left(u\right)\right]du}\right\}}_{{\alpha }^0}.\end{array}$$

(201)

The direct-effect term can be computed immediately, since all quantities are known. The indirect-effect term can be computed, in principle, after solving the 2^nd-LVSS to obtain the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2;u\right)\triangleq {\left[v^{\left(1\right)}\left(u\right),\delta a^{\left(1\right)}\left(u\right)\right]}^{†}$. As has been repeatedly discussed in the foregoing, solving the 2^nd-LVSS is expensive computationally, so this variational function is replaced in the expression of the indirect-effect term by a 2^nd-level adjoint sensitivity function, by following the same steps as outlined in Section 5.2.1. Since there are $3M$ first-order sensitivities of the form $\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}$, $\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},M;\hspace{0.17em}j=1,2,3,$there will be $3M$ distinct 2^nd-level adjoint sensitivity functions, one corresponding to each first-order sensitivity. These $3M$ distinct 2^nd-level adjoint sensitivity functions will be the solutions of the corresponding $3M$ distinct 2^nd-Level Adjoint Sensitivity Systems (2^nd-LASS). Each of these $3M$ 2^nd-LASS will have a distinct source-term on the right-side (each distinct source stemming from the corresponding first-order sensitivities of the form $\partial R\left(\phi ;f\right)/\partial a_j^{\left(i\right)}$), but all of these $3M$ 2^nd-LASS will have the same left-sides, which will have the same form as the left-side of the 2^nd-LASS needed previously, in Subsection 5.2.1 for the computations of the 2^nd-order sensitivities of the response with respect to the components of the feature functions, cf. Equation (164), (176), and (186). Since the left-sides of these 2^nd-LASS represent the (differential) operators that need to be inverted, the actual inversion of these operators needs to be performed once only, and the inverted operator should be stored; subsequently, the inverted operator can be used $3M$ times, operating on $3M$ distinct source terms, to compute the respective $3M$ distinct 2^nd-level adjoint sensitivity functions.

5.2.2.2. Second-Order Sensitivities Stemming from the First-Order Sensitivities with Respect to the Source Properties

The expressions of the first-order sensitivities of the detector response with respect to the parameters that characterize the source which emits the neutrons into the medium are provided in Eqs. (146)‒(150). These expressions have the following generic form:

$$\frac{\partial R\left(\phi ;f\right)}{\partial b_k^{\left(i\right)}}={\omega }_k^{\left(i\right)}\left(\alpha \right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},5;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,2,$$

(202)

where:

$$b_k^{\left(1\right)}\triangleq {\lambda }_k^S;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\omega }_k^{\left(1\right)}\left(\alpha \right)\triangleq \frac{N_k^S{F}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(203)

$$b_k^{\left(2\right)}\triangleq N_k^S;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\omega }_k^{\left(2\right)}\left(\alpha \right)\triangleq \frac{{\lambda }_i^S{F}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(204)

$$b_k^{\left(3\right)}\triangleq F_k^S;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\omega }_k^{\left(3\right)}\left(\alpha \right)\triangleq \frac{{\lambda }_i^S{N}_k^S{\nu }_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(205)

$$b_k^{\left(4\right)}\triangleq {\nu }_k^S;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\omega }_k^{\left(4\right)}\left(\alpha \right)\triangleq \frac{{\lambda }_i^S{N}_k^S{F}_k^S{W}_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,2;$$

(206)

$$b_k^{\left(5\right)}\triangleq {\nu }_k^S;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\omega }_k^{\left(5\right)}\left(\alpha \right)\triangleq \frac{{\lambda }_i^S{N}_k^S{F}_k^S{\nu }_k^S}{\overline{\xi }\hspace{0.17em}\left(\alpha \right){\Sigma }_t\left(\alpha \right)}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=1,2.$$

(207)

The second-order sensitivities stemming from the first-order sensitivities represented by Equation (202) are obtained from the first G-differential of this relation, which has the following expression, by definition, for each $i=1,\mathrm{...},5$; $k=1,2$:

$$\begin{array}{l}\delta \left[\partial R\left(\phi ;f\right)/\partial b_k^{\left(i\right)}\right]\triangleq {\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\left\{\frac{d}{d\epsilon }\left[{\omega }_k^{\left(i\right)}\left({\alpha }^0+\epsilon \delta \alpha \right)\right]\right\}}_{{\alpha }^0,\epsilon =0}\\ +{\left\{{\omega }_k^{\left(i\right)}\left({\alpha }^0\right)\frac{d}{d\epsilon }{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[a^{\left(1\right)}\left(u\right)+\epsilon \delta a^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0,\epsilon =0}\\ =={\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial b_k^{\left(i\right)}\right]\right\}}_{dir}+{\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial b_k^{\left(i\right)}\right]\right\}}_{ind}={\displaystyle \sum _{n=1}^{TP}{\left\{\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial {\alpha }_n\partial b_k^{\left(i\right)}}\delta {\alpha }_n\right\}}_{{\alpha }^0}},\end{array}$$

(208)

where:

$${\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial b_k^{\left(i\right)}\right]\right\}}_{dir}\triangleq {\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}{a}^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\displaystyle \sum _{n=1}^{TP}{\left\{\frac{\partial {\omega }_k^{\left(i\right)}\left(\alpha \right)}{\partial {\alpha }_n}\right\}}_{{\alpha }^0}\delta {\alpha }_n}\hspace{0.17em},$$

(209)

$${\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial b_k^{\left(i\right)}\right]\right\}}_{ind}\triangleq {\left\{{\omega }_k^{\left(i\right)}\left({\alpha }^0\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\delta a^{\left(1\right)}\left(u\right)}\hspace{0.17em}\delta \left(u\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}.$$

(210)

The direct-effect term can be computed immediately, since all quantities are known. The indirect-effect term can be computed, in principle, after solving the 2^nd-LVSS to obtain the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2;u\right)\triangleq {\left[v^{\left(1\right)}\left(u\right),\delta a^{\left(1\right)}\left(u\right)\right]}^{†}$, but this path is expensive computationally, so this variational function is replaced in the expression of the indirect-effect term by a 2^nd-level adjoint sensitivity function, by following the same steps as outlined in Section 5.2.1. Since there are 10 first-order sensitivities of the form $\partial R\left(\phi ;f\right)/\partial b_k^{\left(i\right)}$, $i=1,\mathrm{...},5$; $k=1,2$, there will be 10 distinct 2^nd-level adjoint sensitivity functions, one corresponding to each first-order sensitivity. Thus, there will be 10 distinct 2^nd-Level Adjoint Sensitivity Systems (2^nd-LASS) to be solved, each having a distinct source-term on the right-side, but all of them having the same left-sides as the left-side of the 2^nd-LASS needed previously, in Subsection 5.2.1 for the computations of the 2^nd-order sensitivities of the response with respect to the components of the feature functions, namely Eqs. (164), (176), and (186).

5.2.2.3 Second-Order Sensitivities Stemming from the First-Order Sensitivities with Respect to the Detector Properties

The expressions of the first-order sensitivities of the detector response with respect to the detector’s material properties (i.e., microscopic cross section and atomic number density) are provided in Eqs. (144) and (145). These expressions have the following generic form:

$$\frac{\partial R\left(\phi ;f\right)}{\partial {\zeta }^{\left(i\right)}}=c^{\left(i\right)}\left(\alpha \right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)}\hspace{0.17em}du\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}i=1,2,$$

(211)

where:

$${\zeta }^{\left(1\right)}=N_d;\hspace{0.17em}\hspace{0.17em}{\zeta }^{\left(2\right)}={\sigma }_d;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}^{\left(1\right)}\left(\alpha \right)\hspace{0.17em}=\hspace{0.17em}{\sigma }_d;\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}^{\left(2\right)}\left(\alpha \right)\hspace{0.17em}=\hspace{0.17em}{N}_d.$$

(212)

The second-order sensitivities stemming from the first-order sensitivities represented by Equation (211) are obtained by determining the first G-differential of this relation, which has the following expression, by definition, for each $i=1,2$:

$$\begin{array}{l}\delta \left[\partial R\left(\phi ;f\right)/\partial {\zeta }^{\left(i\right)}\right]\triangleq {\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)\delta \left(u-u_d\right)du}\hspace{0.17em}\right\}}_{{\alpha }^0}{\left\{\frac{d}{d\epsilon }\left[c^{\left(i\right)}\left({\alpha }^0+\epsilon \delta \alpha \right)\right]\right\}}_{{\alpha }^0,\epsilon =0}\\ +{\left\{c^{\left(i\right)}\left({\alpha }^0\right)\frac{d}{d\epsilon }{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\left[{\phi }^0\left(u\right)+\epsilon v^{\left(1\right)}\left(u\right)\right]}\hspace{0.17em}\delta \left(u-u_d\right)\hspace{0.17em}du\right\}}_{{\alpha }^0,\epsilon =0}\\ ={\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial {\zeta }^{\left(i\right)}\right]\right\}}_{dir}+{\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial {\zeta }^{\left(i\right)}\right]\right\}}_{ind}={\displaystyle \sum _{n=1}^{TP}{\left\{\frac{{\partial }^{2}R\left(\phi ;f\right)}{\partial {\alpha }_n\partial {\zeta }^{\left(i\right)}}\delta {\alpha }_n\right\}}_{{\alpha }^0}},\end{array}$$

(213)

where:

$${\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial {\zeta }^{\left(i\right)}\right]\right\}}_{dir}\triangleq {\left\{{\displaystyle \underset{0}{\overset{u_{th}}{\int }}\phi \left(u\right)}\hspace{0.17em}\delta \left(u-u_d\right)\hspace{0.17em}du\right\}}_{{\alpha }^0}{\displaystyle \sum _{n=1}^{TP}{\left\{\frac{\partial c^{\left(i\right)}\left({\alpha }^0\right)}{\partial {\alpha }_n}\right\}}_{{\alpha }^0}\delta {\alpha }_n}\hspace{0.17em},$$

(214)

$${\left\{\delta \left[\partial R\left(\phi ;f\right)/\partial {\zeta }^{\left(i\right)}\right]\right\}}_{ind}\triangleq {\left\{c^{\left(i\right)}\left({\alpha }^0\right){\displaystyle \underset{0}{\overset{u_{th}}{\int }}{v}^{\left(1\right)}\left(u\right)\delta \left(u-u_d\right)\hspace{0.17em}du}\hspace{0.17em}\right\}}_{{\alpha }^0}.$$

(215)

The direct-effect term can be computed immediately, since all quantities are known. The indirect-effect term can be computed, in principle, after solving the 2^nd-LVSS to obtain the 2^nd-level variational sensitivity function $v^{\left(2\right)}\left(2;u\right)\triangleq {\left[v^{\left(1\right)}\left(u\right),\delta a^{\left(1\right)}\left(u\right)\right]}^{†}$, but this path is expensive computationally. As before, this variational function is replaced in the expression of the indirect-effect term by a 2^nd-level adjoint sensitivity function, by following the same steps as outlined in Section 5.2.1. Since there are 2 first-order sensitivities of the form $\partial R\left(\phi ;f\right)/\partial {\zeta }^{\left(i\right)}$, $i=1,2$, there will be 2 distinct 2^nd-level adjoint sensitivity functions, one corresponding to each first-order sensitivity. As before, the two 2^nd-LASS to be solved have distinct source-terms on their right-sides, but both have the same left-sides as the left-side of the 2^nd-LASS needed previously, as in Eqs. (164), (176), and (186).

In summary, the results discussed in Subsection 5.2.2 indicate that computing the 2^nd-order sensitivities of the model response directly with respect to the $TP\triangleq 3M+12$ primary model parameters, ${\alpha }_j$, by applying the 2^nd-CASAM-L methodology requires one large-scale (adjoint) computation for each primary model parameter ${\alpha }_j$, amounting to a total of number of $TP\triangleq 3M+12$ large-scale computations for solving the respective 2^nd-LASS. All of these 2^nd-LASS have the same left-side (which is also the same as needed for computing the 2^nd-order sensitivities of the response with respect to the feature functions by applying the 2^nd-FASAM-L) but have differing sources on their right-sides. The unmixed second order sensitivities ${\partial }^{2}R\left(\phi ;f\right)/\partial {\alpha }_i\partial {\alpha }_i$$i=1,\mathrm{...},TP\triangleq 3M+12$, are computed just once. The mixed second order sensitivities ${\partial }^{2}R\left(\phi ;f\right)/\partial {\alpha }_i\partial {\alpha }_j$, $i\ne j$, are computed twice, involving distinct 2^nd-level adjoint sensitivity functions. Therefore, the symmetry property ${\partial }^{2}R\left(\phi ;f\right)/\partial {\alpha }_i\partial {\alpha }_j={\partial }^{2}R\left(\phi ;f\right)/\partial {\alpha }_j\partial {\alpha }_i$ provides an intrinsic mechanism for verifying the accuracy of the computations of the respective 2^nd-level adjoint sensitivity functions.

6. Concluding Discussion

This work has presented the mathematical framework of the “2^nd-Order Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (2^nd-FASAM-L)” along with an illustrative application to a paradigm model of energy slowing-down of neutrons in an infinitely large homogeneous mixture of materials, as found in many energy-related systems. It has been shown that the 2^nd-FASAM-L is the most efficient methodology for computing exactly the first- and second-order sensitivities of model responses with respect to the features (functions) of model parameters and, subsequently, to the primary model parameters themselves. This efficiency stems from the maximal reduction of the number of adjoint computations (which are “large-scale” computations) which are needed for obtaining these sensitivities. In the extreme case when the model presents no features (functions) of the primary model parameters, the 2^nd-FASAM-L reduces to the 2^nd-CASAM-L (“2^nd-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems) developed by Cacuci [27]. Comparing the mathematical framework of the 2^nd-FASAM-L methodology to the framework of the 2^nd-CASAM-L methodology indicates the following commonalities and distinctions:

• The components $f_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TF$, of the “feature function” $f\left(\alpha \right)\triangleq {\left[f_1\left(\alpha \right),\mathrm{...},f_{TF}\left(\alpha \right)\right]}^{†}$ play within the 2^nd-FASAM-L the same role as played by the components ${\alpha }_j,\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TP$, of the “vector of primary model parameters” $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}$ within the framework of the 2^nd-CASAM-L. Notably, the total number of model parameters is always larger (usually by wide margin) than the total number of components of the feature function $f\left(\alpha \right)$, i.e., $TP\gg TF$.
• The 1^st-FASAM-L and the 1^st-CASAM-L methodologies require a single large-scale “adjoint” computations for solving the 1^st-LASS (1^st-Level Adjoint Sensitivity System), so they are similarly efficient for computing the exact expressions of the first-order sensitivities of a model response to the model’s uncertain parameters, boundaries, and internal interfaces, with a slight computational advantage towards the 1^st-FASAM-L, which requires only $TP$ quadratures, as opposd to $TF$ quadratures required by the 1^st-CASAM-L methodology.
• For computing the exact expressions of the second-order response sensitivities with respect to the primary model’s parameters, the 2^nd-FASAM-L methodology requires as many large-scale “adjoint” computations as there are “feature functions of parameters” $f_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TF$, for solving the left-side of the 2^nd-LASS with $TF$ distinct sources on its right-side. By comparison, the 2^nd-CASAM-L methodology requires $TP$ large-scale computations for solving the same left-side of the 2^nd-LASS but with $TP$ distinct sources. Since $TF\ll TP$, the 2^nd-FASAM-L methodology is considerably more efficient than the 2^nd-CASAM-L methodology for computing the exact expressions of the second-order sensitivities of a model response to the model’s uncertain parameters, boundaries, and internal interfaces.
• Both the 2^nd-FASAM-L and the 2^nd-CASAM-L methodologies are formulated in linearly increasing higher-dimensional Hilbert spaces −as opposed to exponentially increasing parameter-dimensional spaces− thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems. Both the 2^nd-FASAM-L and the 2^nd-CASAM-L methodologies are incomparably more efficient and more accurate than any other methods (statistical, finite differences, etc.) for computing exact expressions of response sensitivities (of any order) with respect to the model’s uncertain parameters, boundaries, and internal interfaces.

Ongoing work aims at generalizing the 2^nd-FASAM-L methodology to enable the exact and most efficient computation of response sensitivities of arbitrarily-high (n^th-) order with respect to features (functions) of model parameters, thus becoming the companion for ‒and most efficient alternative to‒ the n^th-CASAM-L methodology [27], whenever the model comprises features (functions) of model parameters.