Second-Order Adjoint Sensitivity Analysis Methodology for Computing Exactly Response Sensitivities to Uncertain Parameters and Boundaries of Linear Systems: Mathematical Framework

Dan Gabriel Cacuci

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

DOI: 10.4236/ajcm.2020.103018                

Abstract

This work presents the “Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2^nd-CASAM)” for the efficient and exact computation of 1^st- and 2^nd-order response sensitivities to uncertain parameters and domain boundaries of linear systems. The model’s response (i.e., model result of interest) is a generic nonlinear function of the model’s forward and adjoint state functions, and also depends on the imprecisely known boundaries and model parameters. In the practically important particular case when the response is a scalar-valued functional of the forward and adjoint state functions characterizing a model comprising N parameters, the 2^nd-CASAM requires a single large-scale computation using the First-Level Adjoint Sensitivity System (1^st-LASS) for obtaining all of the first-order response sensitivities, and at most N large-scale computations using the Second-Level Adjoint Sensitivity System (2^nd-LASS) for obtaining exactly all of the second-order response sensitivities. In contradistinction, forward other methods would require (N2/2 + 3 N/2) large-scale computations for obtaining all of the first- and second-order sensitivities. This work also shows that constructing and solving the 2^nd-LASS requires very little additional effort beyond the construction of the 1^st-LASS needed for computing the first-order sensitivities. Solving the equations underlying the 1^st-LASS and 2^nd-LASS requires the same computational solvers as needed for solving (i.e., “inverting”) either the forward or the adjoint linear operators underlying the initial model. Therefore, the same computer software and “solvers” used for solving the original system of equations can also be used for solving the 1^st-LASS and the 2^nd-LASS. Since neither the 1^st-LASS nor the 2^nd-LASS involves any differentials of the operators underlying the original system, the 1^st-LASS is designated as a “first-level” (as opposed to a “first-order”) adjoint sensitivity system, while the 2^nd-LASS is designated as a “second-level” (rather than a “second-order”) adjoint sensitivity system. Mixed second-order response sensitivities involving boundary parameters may arise from all source terms of the 2^nd-LASS that involve the imprecisely known boundary parameters. Notably, the 2^nd-LASS encompasses an automatic, inherent, and independent “solution verification” mechanism of the correctness and accuracy of the 2nd-level adjoint functions needed for the efficient and exact computation of the second-order sensitivities.

Keywords

Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2^nd-CASAM), First-Level Adjoint Sensitivity System (1st-LASS), Second-Level Adjoint Sensitivity System (2^nd-LASS), Operator-Type Response, Second-Order Sensitivities to Uncertain Model Boundaries, Second-Order Sensitivities to Uncertain Model Parameters

1. Introduction

The earliest use of adjoint operators for computing exactly and efficiently the first-order sensitivities of responses of a large-scale linear system comprising many parameters has appeared in the report by Wigner [2] on his work on the “nuclear pile” (i.e., nuclear reactor) using the linear neutron transport or diffusion equations. Cacuci [2] [3] conceived the rigorous 1^st-order adjoint sensitivity analysis methodology for generic large-scale nonlinear systems involving operator responses, comprising functional-type responses as particular cases. Representative developments of the adjoint sensitivity analysis methodology were reviewed in the books by Cacuci [4] and Cacuci, Ionescu-Bujor and Navon [5]. The first-order adjoint sensitivity analysis methodology was extended by Cacuci [6] [7] [8] to the “Second-Order Adjoint Sensitivity Analysis Methodology (2^nd-ASAM)” for linear and nonlinear systems, to enable the exact and efficient computation of all of the second-order sensitivities (i.e., functional derivatives) of model responses to parameters. However, all of the above works, as well as the vast majority of the application of the first- and second-order adjoint sensitivity analysis methodology considered that the phase-space location of the physical model’s boundary is perfectly well known, and did not investigate the consequences of imprecisely known (i.e., uncertain) model boundaries.

Komata [9] seems to have been the first author to use an adjoint equation (in this case; the adjoint neutron diffusion equation) for investigating the effects of perturbations in boundary conditions for simple idealized geometries in neutron diffusion theory. First-order perturbation theory in conjunction with adjoint operators was subsequently used [10] - [15], either formally or in conjunction with variational formulations, to obtain approximate first-order sensitivities to boundary parameters of responses that were linear functionals (or ratios thereof) of the neutron flux in the context of linear neutron diffusion or neutron transport problems. Furthermore, the works referenced in the foregoing have considered the computation of response sensitivities to model parameters separately (rather than simultaneously) from the computations of response sensitivities to imprecisely known domain boundary parameters. Also, none of the above works considered responses that are simultaneously functions of the forward and adjoint state functions.

Noteworthy, Cacuci [16] [17] has presented methodologies for the exact computation of 1^st- and 2^nd-order sensitivities of responses that are functionals of both the forward and the adjoint fluxes in a multiplying nuclear system with source. However, these works [16] [17] have specifically considered the neutron transport equation (rather than generic linear problems) within precisely known domains, without considering response sensitivities to uncertain domain parameters.

The aim of this work is to present the novel “Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2^nd-CASAM),” which has the following features that generalize and extend all previously published works on this topic:

2) The system response considered within the 2^nd-CASAM framework is an operator-valued response that depends on both the forward and adjoint state-functions. Functional-valued responses are subsumed as particular cases.

2) The 2^nd-CASAM framework enables the efficient and exact computation of the 1^st- and 2^nd-order response sensitivities not only to uncertain model parameters but also to uncertain domain boundaries of generic linear systems.

This work is structured as follows: Section 2 presents the mathematical framework of the generic linear physical system comprising imprecisely known parameters and boundaries. Section 3 presents the mathematical framework of the 2^nd-CASAM. Section 4 offers concluding remarks. An accompanying work [18] presents an illustrative application of the general 2^nd-CASAM methodology to a paradigm model of generic particle/radiation transmission and/or evolution which has fundamental applications in many fields, including nuclear engineering (e.g., radiation detection, chemical reprocessing of spent reactor fuel, radioactive decay, etc.).

2. Mathematical Modeling of Response-Coupled Forward and Adjoint Linear Systems Comprising Uncertain Responses, Parameters and Boundaries

For the mathematical derivations to be presented in this work, the vectors are column vectors unless explicitly stated otherwise. Bold letters will be used to denote matrices and vectors. Transposition will be indicated by a dagger $(†)$ superscript. A linear physical system is generally represented by means of $N_u$ coupled operator equations of the form

${\displaystyle \underset{j=1}{\overset{N_u}{\sum }}{L}_{ij}\left(x,\alpha \right){\phi }_j\left(x\right)}=Q_i\left(x,\alpha \right),i=1,\cdots ,N_u,x\in {\Omega }_x$ in which the operators $L_{ij}\left(\alpha \right)$ act linearly on the state functions $u_j\left(x\right)$. In matrix form, a linear system can be written as follows:

$L\left(x,\alpha \right)\phi \left(x\right)=Q_L\left(x,\alpha \right),x\in {\Omega }_x$. (1)

The quantities appearing in the Equation (1) are defined as follows:

2) $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{N_{\alpha }}\right)}^{†}$ denotes a $N_{\alpha }$ -dimensional column vector whose components are the physical system’s imprecisely known parameters, which are subject to uncertainties; $\alpha \in {ℝ}^{N_{\alpha }}$, where ${ℝ}^{N_{\alpha }}$ denotes a $N_{\alpha }$ -dimensional subset of the set of real scalars. The vector $\alpha \in {ℝ}^{N_{\alpha }}$ is considered to include any imprecisely known model parameters that may enter into defining the system’s boundary in the phase-se space of independent variables. The symbol “ $\triangleq$ ” will be used to denote “is defined as” or “is by definition equal to.”

2) $x\triangleq \left(x_1,\cdots ,x_{N_x}\right)\in {ℝ}^{N_x}$ denotes the $N_x$ -dimensional phase-space position vector, defined on a phase-space domain denoted as ${\Omega }_x$ which is defined as follows: ${\Omega }_x\triangleq \left\{-\infty \le x_i^{\mathcal{l}}\left(\alpha \right)\le x_i\le x_i^u\left(\alpha \right)\le \infty ;i=1,\cdots ,N_x\right\}$. The lower boundary-point of an independent variable is denoted as $x_i^{\mathcal{l}}\left(\alpha \right)$, where the superscript “ $\mathcal{l}$ ” denotes “lower” (e.g., the inner radius of a sphere or cylinder, the lower range of an energy-variable, etc.), while the corresponding upper boundary-point is denoted as $x_i^u\left(\alpha \right)$, where the superscript “u” denotes “upper” (e.g., the outer radius of a sphere or cylinder, the upper range of an energy-variable, etc.). A typical example of “boundaries that depend on imprecisely known parameters” is provided by the boundary conditions needed for models based on diffusion theory, in which the respective “flux and/or current conditions” for the “boundaries facing vacuum” are imposed on the “extrapolated boundary” of the respective spatial domain. As is well known, the “extrapolated boundary” depends not only on the imprecisely known physical dimensions of the problem’s domain, but also on the medium’s microscopic transport cross sections and atomic number densities. For subsequent derivations, it is convenient to define the “vectors of endpoints” $x^{\mathcal{l}}\left(\alpha \right)\triangleq {\left[x_1^{\mathcal{l}}\left(\alpha \right),\cdots ,x_i^{\mathcal{l}}\left(\alpha \right),\cdots ,x_{N_x}^{\mathcal{l}}\left(\alpha \right)\right]}^{†}$ and $x^u\left(\alpha \right)\triangleq {\left[x_1^u\left(\alpha \right),\cdots ,x_i^u\left(\alpha \right),\cdots ,x_{N_x}^u\left(\alpha \right)\right]}^{†}$. The boundary of ${\Omega }_x$, which will be denoted as $\partial {\Omega }_x$, comprises all of the endpoints $x^{\mathcal{l}}\left(\alpha \right)$ and $x^u\left(\alpha \right)$ of the respective intervals on which the components of $x$ are defined, i.e., $\partial {\Omega }_x\triangleq \left\{x^{\mathcal{l}}\left(\alpha \right)\cup x^u\left(\alpha \right)\right\}$.

3) $\phi \left(x\right)\triangleq {\left[{\phi }_1\left(x\right),\cdots ,{\phi }_{N_{\phi }}\left(x\right)\right]}^{†}$ denotes a $N_{\phi }$ -dimensional column vector whose components represent the system’s dependent variables (also called “state functions”);

4) $Q_L\left(x,\alpha \right)\triangleq {\left[Q_{L,1}\left(x,\alpha \right),\cdots ,Q_{L,N_{\phi }}\left(x,\alpha \right)\right]}^{†}\in E_Q$ denotes a $N_{\phi }$ -dimensional column vector defined on a vector space denoted as $E_Q$, whose elements represent inhomogeneous source terms that depend either linearly or nonlinearly on $\alpha$ ; the components of $Q_L$ may be operators, rather than just functions, acting on $\alpha$ and $x$ ;

5) $L\left(x,\alpha \right)\triangleq \left[L_1\left(x,\alpha \right),\cdots ,L_{N_{\phi }}\left(x,\alpha \right)\right]$ denotes a $N_{\phi }$ -component column vector whose components are operators (including differential, difference, integral, distributions, and/or infinite matrices) acting linearly on $\phi$ and nonlinearly on $\alpha$.

6) All of the equalities in this work are considered to hold in the weak (“distributional”) sense, since the right-sides (“sources”) of Equation (1) and of other various equations to be derived in this work may contain distributions (“generalized functions/functionals”), particularly Dirac-distributions and derivatives and/or integrals thereof.

When differential operators appear in Equation (1), a corresponding set of boundary and/or initial conditions, which are essential to define the domain of $L\left(x,\alpha \right)$, must also be given. Since this work considers only systems that are linear in the state function $\phi \left(x\right)$, the accompanying boundary and/or initial conditions must also be linear in $\phi \left(x\right)$, which means that they can be represented in operator form as follows:

$B_L\left[x^{\mathcal{l}}\left(\alpha \right),x^u\left(\alpha \right)\right]\phi \left[x^{\mathcal{l}}\left(\alpha \right),x^u\left(\alpha \right)\right]-C_L\left[x^{\mathcal{l}}\left(\alpha \right),x^u\left(\alpha \right)\right]=0,x\in \partial {\Omega }_x.$ (2)

In Equation (2), the operator $B_L\left(x,\alpha \right)\triangleq \left[B_{ij}\left(x,\alpha \right)\right];i=1,\cdots ,N_B;j=1,\cdots ,N_u$ is a matrix comprising, as components, operators that act linearly on $\phi \left(x\right)$ and nonlinearly on $\alpha$ ; the quantity $N_B$ denotes the total number of boundary and initial conditions. The operator $C_L\left(x,\alpha \right)\triangleq {\left[C_1\left(x,\alpha \right),\cdots ,C_{N_B}\left(x,\alpha \right)\right]}^{†}$ is a $N_B$ -dimensional vector comprising components that are operators acting, in general, nonlinearly on $\alpha$.

The nominal solution of Equations (1) and (2) is denoted as ${\phi }^0\left(x\right)$, and is obtained by solving these equations at the nominal parameter values ${\alpha }^0$. In other words, the vectors ${\phi }^0\left(x\right)$ and ${\alpha }^0$ satisfy the following equations:

$L\left(x,{\alpha }^0\right){\phi }^0\left(x\right)=Q\left(x,{\alpha }^0\right),x\in {\Omega }_x,$ (3)

$\begin{array}{l}{B}_L\left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right]{\phi }^0\left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right]\\ -C_L\left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right]=0,x\in \partial {\Omega }_x.\end{array}$ (4)

Equations (3) and (4) represent the “base-case” or nominal state of the physical system. Throughout this work, the superscript “0” will be used to denote “nominal” or “expected” values.

The vector-valued function $\phi \left(x\right)$ is considered to be the unique nontrivial solution of the physical problem described by Equations (1) and (2). The linear operator $L\left(x,\alpha \right)$ is considered in this work to admit an adjoint operator, which will be denoted as $L^{*}\left(x,\alpha \right)$. The specific form of the adjoint operator $L^{*}\left(x,\alpha \right)$ depends upon the specific inner product and associated Hilbert space chosen for defining it. In view of Equations (1) and (2), the adjoint operator $L^{*}\left(x,\alpha \right)$ is most often defined in a Hilbert space (which will henceforth be denoted as $H_{\phi }$ ) endowed with the following inner product, denoted as ${\langle {\phi }^{\left(a\right)},{\phi }^{\left(b\right)}\rangle }_{\left(\phi \right)}$, between two the vectors ${\phi }^{\left(a\right)}\in H_{\phi }$ and ${\phi }^{\left(b\right)}\in H_{\phi }$, each having square integrable components and the same structure as the $N_u$ -dimensional column vector $\phi \left(x\right)$ :

$\begin{array}{l}{\langle {\phi }^{\left(a\right)},{\phi }^{\left(b\right)}\rangle }_{\left(\phi \right)}\triangleq {\displaystyle \underset{x_1^{\mathcal{l}}\left(\alpha \right)}{\overset{x_1^u\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{x_i^{\mathcal{l}}\left(\alpha \right)}{\overset{x_i^u\left(\alpha \right)}{\int }}\mathrm{...}{\displaystyle \underset{x_{N_x}^{\mathcal{l}}\left(\alpha \right)}{\overset{x_{N_x}^u\left(\alpha \right)}{\int }}\left({\displaystyle \underset{i=1}{\overset{N_u}{\sum }}{\phi }_i^{\left(a\right)}{\phi }_i^{\left(b\right)}}\right)\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_i\cdots \text{d}{x}_{N_x}}}}\\ \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\triangleq {\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left(\alpha \right)}{\overset{x_i^u\left(\alpha \right)}{\int }}{\phi }^{\left(a\right)}\left(x\right)\cdot {\phi }^{\left(b\right)}\left(x\right)\text{d}x}.}\end{array}$ (5)

In Equation (5), the “product” notation ${\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left(\alpha \right)}{\overset{x_i^u\left(\alpha \right)}{\int }}\left[\right]\text{d}{x}_i}}$ compactly denotes the respective multiple integrals, while the dot indicates the “scalar product of two vectors” defined as follows:

${\phi }^{\left(a\right)}\left(x\right)\cdot {\phi }^{\left(b\right)}\left(x\right)\triangleq {\displaystyle \underset{i=1}{\overset{N_u}{\sum }}{\phi }_i^{\left(a\right)}{\phi }_i^{\left(b\right)}}$. (6)

Formally, the inner product introduced in Equation (5) defines a self-dual Hilbert space, which shall be denoted as $H_{\phi }$. It is important to note that the inner product defined in Equation (5) is continuous in $\alpha$, holding at any value particular value of $\alpha$ in the neighborhood of the nominal parameter values ${\alpha }^0$, including at ${\alpha }^0$.

The formal adjoint operator $L^{*}\left(x,\alpha \right)$ is defined through the following relationship, which holds for an arbitrary vector $\psi \left(x\right)\in H_{\phi }$ :

${\langle \psi \left(x\right),L\left(x,\alpha \right)\phi \left(x\right)\rangle }_{\left(\phi \right)}={\langle L^{*}\left(x,\alpha \right)\psi \left(x\right),\phi \left(x\right)\rangle }_{\left(\phi \right)}+{\left\{P^{\left(1\right)}\left(\phi ;\alpha ;\psi \right)\right\}}_{\partial {\Omega }_x}$ (7)

In Equation (7), the formal adjoint operator $L^{*}\left(x,\alpha \right)$ is the $N_u\times N_u$ matrix

$L^{*}\left(x,\alpha \right)\triangleq \left[L_{ji}^{*}\left(x,\alpha \right)\right],\left(i,j=1,\cdots ,N_u\right),$ (8)

comprising elements $L_{ji}^{*}\left(x,\alpha \right)$ obtained by transposing the formal adjoints of the operators $L_{ij}\left(x,\alpha \right)$, while ${\left\{P^{\left(1\right)}\left(\phi ;\alpha ;\psi \right)\right\}}_{\partial {\Omega }_x}$ denotes the associated bilinear form evaluated on $\partial {\Omega }_x$. The adjoint state function $\psi \left(x\right)$ is the solution of an adjoint system which can be written as follows:

$L^{*}\left(x,\alpha \right)\psi \left(x\right)=Q_A\left(x,\alpha \right),x\in {\Omega }_x,$ (9)

where the source $Q_A\left(x,\alpha \right)\triangleq {\left[Q_{A,1}\left(x,\alpha \right),\cdots ,Q_{A,N_{\phi }}\left(x,\alpha \right)\right]}^{†}$ for the adjoint equation is usually related to the system response under consideration. The domain of $L^{*}\left(x,\alpha \right)$ is determined by selecting appropriate adjoint boundary and/or initial conditions, represented here in operator form as follows:

$B_A\left[x^{\mathcal{l}}\left(\alpha \right),x^u\left(\alpha \right)\right]\psi \left[x^{\mathcal{l}}\left(\alpha \right),x^u\left(\alpha \right)\right]-C_A\left[x^{\mathcal{l}}\left(\alpha \right),x^u\left(\alpha \right)\right]=0,x\in \partial {\Omega }_x.$ (10)

In Equation (10), the subscript “A” indicates “adjoint,” and the letter “B” indicates “boundary and/or initial conditions.” The nominal value of the adjoint state function $\psi \left(x\right)$ will be denoted as ${\psi }^0\left(x\right)$, and is obtained by solving the adjoint system at the nominal parameter values ${\alpha }^0$ :

$L^{*}\left(x,{\alpha }^0\right){\psi }^0\left(x\right)=Q_A\left(x,{\alpha }^0\right),x\in {\Omega }_x,$ (11)

$\begin{array}{l}{B}_A\left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right]{\psi }^0\left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right]\\ -C_A\left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right]=0,x\in \partial {\Omega }_x.\end{array}$ (12)

The system response (i.e., result of interest) associated with the problem modeled by Equations (1) and (2) is typically a real-valued nonlinear operator (function) of the system’s forward and adjoint state-functions (i.e., dependent variables) and parameters, denoted as $R\left[\phi \left(x\right),\psi \left(x\right);\alpha ;x\right]$, which can be generally represented in the form

$R\left[\phi \left(x\right),\psi \left(x\right);\alpha ;x\right]=F\left[\phi \left(x\right),\psi \left(x\right);\alpha ;x\right],$ (13)

where $F\left(\phi ,\psi ;\alpha ;x\right)$ is a suitably differentiable function with respect to its arguments.

3. The Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2^nd-CASAM)

As has been already mentioned in the foregoing, the model and boundary parameters are considered to be imprecisely known quantities. Their true values may differ from their nominal values by quantities denoted as $\delta \alpha \triangleq \left(\delta {\alpha }_1,\cdots ,\delta {\alpha }_{N_{\alpha }}\right)$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$. Since the forward and the adjoint state functions $\phi \left(x\right)$ and $\psi \left(x\right)$ are related to the model and boundary parameters $\alpha$ through Equations (1), (2), (9) and (10), it follows that variations $\delta \alpha$ in the model and boundary parameter will cause corresponding variations

$\delta \phi \left(x\right)\triangleq {\left[\delta {\phi }_1\left(x\right),\cdots ,\delta {\phi }_{N_{\phi }}\left(x\right)\right]}^{†}$ around the nominal value ${\phi }^0\left(x\right)$ of the forward state function, as well as variations $\delta \psi \left(x\right)\triangleq {\left[\delta {\psi }_1\left(x\right),\cdots ,\delta {\psi }_{N_{\phi }}\left(x\right)\right]}^{†}$ around the nominal value ${\psi }^0\left(x\right)$ of the adjoint state function. In turn, the variations in the model’s parameters, boundaries and state functions will induce variations in the system’s response. For the derivations to follow, it is convenient to define the following vectors in the combined phase-space of state-functions, model and boundary parameters: 1) the vector $e\left(x\right)\triangleq {\left[\phi \left(x\right),\psi \left(x\right);\alpha \right]}^{†}$, comprising as components the state functions and parameters; 2) the vector $e^0\left(x\right)\triangleq {\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right),{\alpha }^0\right]}^{†}$ of nominal values of the state functions and parameters; 3) the vector $h\left(x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \left(x\right);\delta \alpha \right]}^{†}$ of variations $\delta \phi \left(x\right)\triangleq \phi \left(x\right)-{\phi }^0\left(x\right)$, $\delta \psi \left(x\right)\triangleq \psi \left(x\right)-{\psi }^0\left(x\right)$, and $\delta \alpha \triangleq \alpha -{\alpha }^0$ in the state functions and parameters around their respective nominal values. The 2^nd-CASAM formalism underlying the computation of the 1^st-order sensitivities of the response $R\left(\phi ,\psi ;\alpha \right)$ is presented in Section 3.2, while the mathematical methodology for computing the 2^nd-order sensitivities is presented in Section 3.2.

3.1. Derivation of the 1^st-Level Adjoint Sensitivity System (1^st-LASS) for Computing Exactly and Efficiently the 1^st-Order Response Sensitivities to Model and Boundary Parameters

As shown by Cacuci [2] [3], in front of the sentence, the most general definition of the 1^st-order total sensitivity of the operator-valued model response $R\left(\phi ,\psi ;\alpha ;x\right)$ to variations $h\left(x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \left(x\right);\delta \alpha \right]}^{†}$ around the nominal values $e^0\left(x\right)\triangleq {\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right),{\alpha }^0\right]}^{†}$, respectively, is provided by the first-order Gateaux-variation $\delta R\left(e^0;h;x\right)$ of $R\left(\phi ,\psi ;\alpha ;x\right)$, which is defined as follows:

$\delta R\left(e^0;h;x\right)\triangleq {\left\{\frac{\text{d}}{\text{d}\epsilon }\left[F\left({\phi }^0+\epsilon \delta \phi ;{\psi }^0+\epsilon \delta \psi ;{\alpha }^0+\epsilon \delta \alpha ;x\right)\right]\right\}}_{\epsilon =0},$ (14)

for $\epsilon \in F$, where $F$ denotes the field of real scalars, and all (i.e., arbitrary) vectors $h\left(x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \left(x\right);\delta \alpha \right]}^{†}$ in a neighborhood $\left(e^0+\epsilon h\right)$ around $e^0\left(x\right)\triangleq {\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right),{\alpha }^0\right]}^{†}$. The G-differential $\delta R\left(e^0;h;x\right)$ is an operator defined on the same domain as $R\left(e;x\right)$, has the same range as $R\left(e;x\right)$, and provides the total 1^st-order sensitivity of $R\left(e;x\right)$ with respect to variations in the model’s parameters and state functions. The G-differential $\delta R\left(e^0;h;x\right)$ satisfies the relation $R\left(e^0+\epsilon h;x\right)-R\left(e^0;x\right)=\delta R\left(e^0;h;x\right)+\Delta \left(h;x\right)$, with $\underset{\epsilon \to 0}{\lim }\left[\Delta \left(\epsilon h;x\right)\right]/\epsilon =0$.

The existence of the G-variation $\delta R\left(e^0;h;x\right)$ does not guarantee its numerical computability. Numerical methods most often require that $\delta R\left(e^0;h;x\right)$ be linear in $h$ in a neighborhood $\left(e^0+\epsilon h\right)$ around $e^0$. The necessary and sufficient conditions for the G-variation $\delta R\left(e^0;h;x\right)$ of a nonlinear operator $R\left(e;x\right)$ to be linear in $h$ in a neighborhood $\left(e^0+\epsilon h\right)$ around $e^0$, and thus admit a total first-order G-derivative, are as follows:

1) $F\left(e;x\right)$ must satisfy a weak Lipschitz condition at $e^0$ ; (15)

2) $F\left(e;x\right)$ must satisfy the following condition for $h_1\text{,}{h}_2\in E\text{;}\epsilon \in F$ :

$F\left(e^0+\epsilon h_1+\epsilon h_2;x\right)-F\left(e^0+\epsilon h_1;x\right)-F\left(e^0+\epsilon h_2;x\right)+F\left(e^0;x\right)=o\left(\epsilon \right).$ (16)

Numerical methods (e.g., Newton’s method and variants thereof) for solving Equations (1) and (2) also require the existence of the 1^st-order G-derivatives of original model equations, in which case the components of the operators which appear in these equations must also satisfy the conditions described in Equations (15) and (16). Of course, if the first-order G-derivatives of the system’s response do not exist, the computation of higher-order response sensitivities (G-derivatives) would be moot. Therefore, the conditions described in Equations (15) and (16) will henceforth be considered to be satisfied by the operators underlying the physical system, in which case the partial G-derivatives of $R\left(e\right)$ at $e^0$ with respect to $\phi$, $\psi$ and $\alpha$ exist. These derivatives are row vectors defined as $\frac{\partial R\left(e;x\right)}{\partial \phi }\triangleq \left[\frac{\partial F\left(e;x\right)}{\partial {\phi }_1},\cdots ,\frac{\partial F\left(e;x\right)}{\partial {\phi }_{N_{\phi }}}\right]$, $\frac{\partial R\left(e;x\right)}{\partial \psi }\triangleq \left[\frac{\partial F\left(e;x\right)}{\partial {\psi }_1},\cdots ,\frac{\partial F\left(e;x\right)}{\partial {\psi }_{N_{\phi }}}\right]$, and $\partial R\left(e;x\right)/\partial \alpha \triangleq \left[\partial F\left(e;x\right)/\partial {\alpha }_1,\cdots ,\partial F\left(e;x\right)/\partial {\alpha }_{N_{\alpha }}\right]$, respectively, so that the first-order G-variation $\delta R\left(e^0;h;x\right)$ can be written in the form

$\delta R\left(e^0;h;x\right)={\left\{\delta R\left(e^0;h;x\right)\right\}}_{direct}+{\left\{\delta R\left(e^0;h;x\right)\right\}}_{indirect},$ (17)

where

${\left\{\delta R\left(e^0;h;x\right)\right\}}_{direct}={\left\{\frac{\partial F\left(e;x\right)}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha ,$ (18)

${\left\{\delta R\left(e^0;h;x\right)\right\}}_{indirect}\triangleq {\left\{\frac{\partial F\left(e;x\right)}{\partial \phi }\right\}}_{\left(e^0\right)}\delta \phi \left(x\right)+{\left\{\frac{\partial F\left(e;x\right)}{\partial \psi }\right\}}_{\left(e^0\right)}\delta \psi \left(x\right).$ (19)

The notation ${\left\{\right\}}_{\left(e^0\right)}$ indicates that the quantity within the braces is to be evaluated at $e^0\left(x\right)\triangleq {\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right),{\alpha }^0\right]}^{†}$. The quantity ${\left\{\delta R\left(e^0;h;x\right)\right\}}_{direct}$ in Equation (18) is called the “direct-effect term” because it can be computed after the base-case values $e^0$ become available. On the other hand, the “indirect-effect term” ${\left\{\delta R\left(e^0;h;x\right)\right\}}_{indirect}$ defined in Equation (19) can be quantified only after having determined the variations $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$ in terms of the variations $\delta \alpha$. The first-order relationship between the variations $\delta \phi \left(x\right)$ and $\delta \alpha$ is determined by determining the G-differentials of Equations (1) and (2). Thus, applying the definition of the G-differential to Equations (1) and (2) yields the following equations:

$L\left(x,{\alpha }^0\right)\delta \phi \left(x\right)=Q_L^{\left(1\right)}\left({\phi }^0,{\alpha }^0;\delta \alpha \right),x\in {\Omega }_x,$ (20)

$B_L^{\left(1\right)}\left\{\phi \left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right];{\alpha }^0;\delta \phi ,\delta \alpha \right\}=0,x\in \partial {\Omega }_x.$ (21)

where the superscript “(1)” indicates “1^st-Level” and the letter “B” indicates “boundary and/or initial conditions,” and where the following definitions were used:

$Q_L^{\left(1\right)}\left({\phi }^0,{\alpha }^0;\delta \alpha \right)\triangleq {\left\{\frac{\partial \left[Q_L\left(\alpha \right)-L\left(\alpha \right)\phi \right]}{\partial \alpha }\delta \alpha \right\}}_{\left(e^0\right)},$ (22)

$\begin{array}{l}{B}_L^{\left(1\right)}{\left\{\phi \left(x^{\mathcal{l}},x^u\right);{\alpha }^0;\delta \phi ,\delta \alpha \right\}}_{\left(e^0\right)}\\ \triangleq {\left\{\frac{\partial \left[B_L\left(\alpha \right)\phi -C_L\left(\alpha \right)\right]}{\partial \alpha }\delta \alpha \right\}}_{\left(e^0\right)}+{\left\{B_L\left(\alpha \right)\delta \phi \right\}}_{\left(e^0\right)}\\ +{\left\{{\left[\frac{\partial \left[B_L\left(\alpha \right)\phi \right]}{\partial x^u}\frac{\partial x^u\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]}_{x=x^u}\right\}}_{\left(e^0\right)}\\ +{\left\{{\left[\frac{\partial \left[B_L\left(\alpha \right)\phi \right]}{\partial x^{\mathcal{l}}}\frac{\partial x^{\mathcal{l}}\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]}_{x=x^{\mathcal{l}}}\right\}}_{\left(e^0\right)}.\end{array}$ (23)

The partial G-derivatives appearing in Equation (22) are matrices defined as follows:

$\frac{\partial Q_L\left(\alpha \right)}{\partial \alpha }\triangleq \left(\begin{array}{ccc}\frac{\partial Q_{L,1}}{\partial {\alpha }_1}& \cdots & \frac{\partial Q_{L,1}}{\partial {\alpha }_{N_{\alpha }}}\\ ⋮& \ddots & ⋮\\ \frac{\partial Q_{L,N_u}}{\partial {\alpha }_1}& \cdots & \frac{\partial Q_{L,N_u}}{\partial {\alpha }_{N_{\alpha }}}\end{array}\right),$ (24)

$\frac{\partial \left[L\left(\alpha \right)\phi \right]}{\partial \alpha }\triangleq \left(\begin{array}{ccc}\frac{\partial \left[{\displaystyle \underset{j=1}{\overset{N_u}{\sum }}{L}_{1,j}\left(\alpha \right){\phi }_j}\right]}{\partial {\alpha }_1}& \cdots & \frac{\partial \left[{\displaystyle \underset{j=1}{\overset{N_u}{\sum }}{L}_{1,j}\left(\alpha \right){\phi }_j}\right]}{\partial {\alpha }_{N_{\alpha }}}\\ ⋮& \ddots & ⋮\\ \frac{\partial \left[{\displaystyle \underset{j=1}{\overset{N_u}{\sum }}{L}_{N_u,j}\left(\alpha \right){\phi }_j}\right]}{\partial {\alpha }_1}& \cdots & \frac{\partial \left[{\displaystyle \underset{j=1}{\overset{N_u}{\sum }}{L}_{N_u,j}\left(\alpha \right){\phi }_j}\right]}{\partial {\alpha }_{N_{\alpha }}}\end{array}\right).$ (25)

The partial G-derivatives that appear in Equation (23) are also matrices, having structures and components similar to those defined in Equations (24) and (25).

The first-order relationship between the variations $\delta \psi \left(x\right)$ and $\delta \alpha$ is determined by taking the G-differentials of Equations (9) and (10), which yields the following system of equations:

$L^{*}\left(x,{\alpha }^0\right)\delta \psi \left(x\right)=Q_A^{(1)}\left({\psi }^0,{\alpha }^0;\delta \alpha \right),x\in {\Omega }_x,$ (26)

$B_A^{\left(1\right)}\left\{\psi \left[x^{\mathcal{l}}\left({\alpha }^0\right),x^u\left({\alpha }^0\right)\right];{\alpha }^0;\delta \phi ,\delta \alpha \right\}=0,x\in \partial {\Omega }_x.$ (27)

where the superscript “(1)” indicates “1^st-Level,” the letter “A” indicates “adjoint,” the letter “B” indicates “boundary and/or initial conditions,” and where the following definitions were used:

$Q_A^{\left(1\right)}\left({\psi }^0,{\alpha }^0;\delta \alpha \right)\triangleq {\left\{\frac{\partial \left[Q_A\left(\alpha \right)-L^{*}\left(\alpha \right)\psi \right]}{\partial \alpha }\delta \alpha \right\}}_{\left(e^0\right)},$ (28)

$\begin{array}{l}{B}_A^{\left(1\right)}{\left\{\psi \left(x^{\mathcal{l}},x^u\right);{\alpha }^0;\delta \psi ,\delta \alpha \right\}}_{\left(e^0\right)}\\ \triangleq {\left\{\frac{\partial \left[B_A\left(\alpha \right)\psi -C_A\left(\alpha \right)\right]}{\partial \alpha }\delta \alpha \right\}}_{\left(e^0\right)}+{\left\{B_A\left(\alpha \right)\delta \psi \right\}}_{\left(e^0\right)}\\ +{\left\{{\left[\frac{\partial \left[B_A\left(\alpha \right)\psi \right]}{\partial x^u}\frac{\partial x^u\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]}_{x=x^u}\right\}}_{\left(e^0\right)}\\ +{\left\{{\left[\frac{\partial \left[B_A\left(\alpha \right)\psi \right]}{\partial x^{\mathcal{l}}}\frac{\partial x^{\mathcal{l}}\left(\alpha \right)}{\partial \alpha }\delta \alpha \right]}_{x=x^{\mathcal{l}}}\right\}}_{\left(e^0\right)}.\end{array}$ (29)

The system of equations comprising Equations (20), (21), (26) and (27) are called the “1^st-Level Forward Sensitivity System” (1^st-LFSS). In principle, the 1^st-LFSS could to be solved for each possible vector of parameter variations $\delta \alpha$ in a neighborhood around ${\alpha }^0$ to obtain $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$ as functions of $\delta \alpha$. In turn, the functions $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$ thus obtained could be used in Equation (17) to compute the total response sensitivity $\delta R\left(e^0;h\right)$. Computing the (total) response sensitivity $\delta R\left(e^0;h\right)$ by using the $\delta \alpha$ -dependent solutions $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$ of the 1^st-LFSS is called the Forward Sensitivity Analysis Method/Procedure (FSAM). From the standpoint of computational costs and effort, the FSAM would require $2N_{\alpha }$ large-scale computations to solve the 1^st-LFSS. Therefore, the FSAM is advantageous to employ only if, in the problem under consideration, the number $N_{\alpha }$ of model and boundary parameters is considerably less than the number of responses of interest. This is rarely the case in practice, however, since most problems of practical interest are characterized by many model parameters and comparatively few responses. In such situations, it is not economical to employ the FSAM since it becomes prohibitively expensive to solve repeatedly the $\delta \alpha$ -dependent 1^st-LFSS in order to determine the variations $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$ for all possible variations $\delta \alpha$.

In most practical situations, the number of model parameters exceeds significantly the number of responses of interest. In such cases, the exact computation of the 1^st-order response sensitivities can be performed most efficiently by using the 1^st-order adjoint sensitivity analysis methodology introduced by Cacuci [2] [3] and extended to 2^nd-order by Cacuci [4] [5] [6], which expresses the indirect effect defined in Equation (35) in terms of a suitably defined 1^st-level adjoint function, thus eliminating the need to compute the variations $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$ by solving the 1^st-LFSS. As has been first shown by Cacuci [3], the adjoint sensitivity analysis methodology cannot be applied directly to an operator-valued response, but only to responses that are functionals (i.e., scalar-valued operators) of the state functions. For this purpose, the operator-valued response $R\left(\phi ,\psi ;\alpha \right)$ is expressed through its generalized Fourier (spectral) expansion

$R\left(e\right)={\displaystyle \underset{f=0}{\overset{NF}{\sum }}F\left(f;e\right)p_f}\left(x\right),$ (30)

where the Fourier (spectral) coefficients $F\left(f;e\right)$ are functionals of the parameters, the forward and adjoint state variables, defined as follows:

$F\left(f;e\right)\triangleq {\displaystyle \underset{x_1^{\mathcal{l}}\left(\alpha \right)}{\overset{x_1^u\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{x_i^{\mathcal{l}}\left(\alpha \right)}{\overset{x_i^u\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{x_{N_x}^{\mathcal{l}}\left(\alpha \right)}{\overset{x_{N_x}^u\left(\alpha \right)}{\int }}F\left(e\right)p_f\left(x\right)\text{d}{x}_1\cdots \text{d}{x}_i\cdots \text{d}{x}_{N_x}}}}.$ (31)

and where the quantities $p_f\left(x\right)$ are the basis functions (e.g., orthogonal polynomials, trigonometric functions, etc.) chosen for the domain ${\Omega }_x$. Although the expansion/summation in Equation (30) comprises in principle infinitely many terms, only a finite number, denoted as NF, of terms will be computed in practice. It is therefore important to choose the basis-functions $p_f\left(x\right)$ so as to minimize the total number of terms NF that would be needed to reproduce the response $R\left(e\right)$ within an a priori chosen error criterion. In particular, if the model response is ab initio simply a functional (as opposed to a function or operator) of the forward and/or adjoint functions $\phi$, $\psi$, then the expansion in Equation (31) would contain just the single $f=0$ term. The nominal value of the functional defined in Equation (31) is denoted as follows:

$F\left(f;e^0\right)\triangleq {\displaystyle \underset{x_1^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_1^u\left({\alpha }^0\right)}{\int }}\cdots {\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\cdots {\displaystyle \underset{x_{N_x}^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_{N_x}^u\left({\alpha }^0\right)}{\int }}F\left(e^0\right)p_f\left(x\right)\text{d}{x}_1\cdots \text{d}{x}_i\cdots \text{d}{x}_{N_x}}}}.$ (32)

The 1^st-order G-differential of the response $R\left(\phi ,\psi ;\alpha \right)$ defined in Equation (30) is

$\delta R\left(e^0;h;x\right)={\displaystyle \underset{f=0}{\overset{NF}{\sum }}\delta F\left(f;e^0;h\right)p_f}\left(x\right),$ (33)

where the functional $\delta F\left(f;e^0;h\right)$ has the following expression:

$\begin{array}{l}\delta F\left(f;e^0;h\right)\\ \triangleq {\left\{\frac{\text{d}}{\text{d}\epsilon }{\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0+\epsilon \delta \alpha \right)}{\overset{x_i^u\left({\alpha }^0+\epsilon \delta \alpha \right)}{\int }}F\left({\phi }^0+\epsilon \delta \phi ,{\psi }^0+\epsilon \delta \psi ;{\alpha }^0+\epsilon \delta \alpha \right)p_f\left(x\right)\text{d}{x}_i}}\right\}}_{\epsilon =0}\\ ={\left\{\delta F\left(f;e^0;\delta \alpha \right)\right\}}_{dir}+{\left\{\delta F\left(f;e^0;h\right)\right\}}_{ind},\end{array}$ (34)

with the “indirect-effect term” ${\left\{\delta F\left(f;e^0;\delta \alpha \right)\right\}}_{ind}$ defined as

$\begin{array}{l}{\left\{\delta F\left(f;e^0;h\right)\right\}}_{ind}\\ \triangleq {\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{\left[\frac{\partial F\left(f;e;\alpha \right)}{\partial \phi }{p}_f\left(x\right)\right]\delta \phi +\left[\frac{\partial F\left(f;e;\alpha \right)}{\partial \psi }{p}_f\left(x\right)\right]\delta \psi \right\}}_{\left(e^0\right)}}},\end{array}$ (35)

and the “direct-effect term” ${\left\{\delta F\left(f;e^0;\delta \alpha \right)\right\}}_{dir}$ defined as

$\begin{array}{l}{\left\{\delta F\left(f;e^0;\delta \alpha \right)\right\}}_{dir}\triangleq {\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{\frac{\partial F\left(f;e;\alpha \right)}{\partial \alpha }{p}_f\left(x\right)\right\}}_{\left(e^0\right)}\delta \alpha }}\\ +{\displaystyle \underset{j=1}{\overset{N_x}{\sum }}{\displaystyle \underset{\begin{array}{l}i=1\\ i\ne j\end{array}}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{F\left[f;\phi \left(.\text{},x_j^u\left(\alpha \right),.\right),\psi \left(.,x_j^u\left(\alpha \right),.\right);\alpha \right]p_f\left[x_j^u\left(\alpha \right)\right]\frac{\partial x_j^u\left(\alpha \right)}{\partial \alpha }\delta \alpha \right\}}_{\left(e^0\right)}}}}\\ -{\displaystyle \underset{j=1}{\overset{N_x}{\sum }}{\displaystyle \underset{\begin{array}{l}i=1\\ i\ne j\end{array}}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{F\left[f;\phi \left(.\text{},x_j^{\mathcal{l}}\left(\alpha \right),.\right),\psi \left(.,x_j^{\mathcal{l}}\left(\alpha \right),.\right);\alpha \right]p_f\left[x_j^{\mathcal{l}}\left(\alpha \right)\right]\frac{\partial x_j^{\mathcal{l}}\left(\alpha \right)}{\partial \alpha }\delta \alpha \right\}}_{\left(e^0\right)}}}}.\end{array}$ (36)

The appearance of the variations $\delta \phi$ and $\delta \psi$ in the indirect-effect term ${\left\{\delta F\left(f;e^0;h\right)\right\}}_{ind}$ defined in Equation (35) can be eliminated by expressing the right-side of Equation (35) in terms of adjoint functions that are the solutions of a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) which is constructed following the steps outlined by Cacuci [4] [5] [6] :

1) Introduce a Hilbert space, denoted as $H^{\left(1\right)}$, having as elements two-component vectors of the form $a^{\left(1\right)}\left(x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}\in H^{\left(1\right)}$, with square integrable components having the same structure as the $N_{\phi }$ -dimensional vectors $\phi \left(x\right)$ and/or $\psi \left(x\right)$. The inner product between two vectors $a^{\left(1\right)}\left(x\right)\triangleq {\left[a_1^{\left(1\right)}\left(x\right),a_2^{\left(1\right)}\left(x\right)\right]}^{†}\in H^{\left(1\right)}$ and $b^{\left(1\right)}\left(x\right)\triangleq {\left[b_1^{\left(1\right)}\left(x\right),b_2^{\left(1\right)}\left(x\right)\right]}^{†}\in H^{\left(1\right)}$ will be denoted as ${\langle a^{\left(1\right)}\left(x\right),b^{\left(1\right)}\left(x\right)\rangle }_{\left(1\right)}$ and is defined as follows:

$\begin{array}{c}{\langle a^{\left(1\right)}\left(x\right),b^{\left(1\right)}\left(x\right)\rangle }_{\left(1\right)}\triangleq {\displaystyle \underset{j=1}{\overset{2}{\sum }}{\langle a_j^{\left(1\right)}\left(x\right)\cdot b_j^{\left(1\right)}\left(x\right)\rangle }_{\left(\phi \right)}}\\ ={\displaystyle \underset{j=1}{\overset{2}{\sum }}{\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left(\alpha \right)}{\overset{x_i^u\left(\alpha \right)}{\int }}{a}_j^{\left(1\right)}\left(x\right)b_j^{\left(1\right)}\left(x\right)\text{d}{x}_i}}}.\end{array}$ (37)

In particular, the inner product defined in Equation (37) also holds at the nominal parameter values ${\alpha }^0$.

2) Using Equation (37), construct the inner product of Equations (20) and (26) with a vector ${\psi }^{\left(1\right)}\left(x\right)\triangleq {\left[{\psi }_1^{\left(1\right)}\left(x\right),{\psi }_2^{\left(1\right)}\left(x\right)\right]}^{†}\in H^{\left(1\right)}$ to obtain the following relation:

$\begin{array}{l}{\left\{{\langle {\left(\begin{array}{c}{\psi }_1^{\left(1\right)}\\ {\psi }_2^{\left(1\right)}\end{array}\right)}^{†},\left(\begin{array}{cc}L\left(\alpha \right)& 0\\ 0& L^{*}\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}\delta \phi \\ \delta \psi \end{array}\right)\rangle }_{\left(1\right)}\right\}}_{\left(e^0\right)}\\ ={\left\{{\langle {\left(\begin{array}{c}{\psi }_1^{\left(1\right)}\\ {\psi }_2^{\left(1\right)}\end{array}\right)}^{†},\left(\begin{array}{c}{Q}_L^{\left(1\right)}\left(\phi ;\alpha ;\delta \alpha \right)\\ Q_A^{\left(1\right)}\left(\psi ;\alpha ;\delta \alpha \right)\end{array}\right)\rangle }_{\left(1\right)}\right\}}_{\left(e^0\right)}\end{array}$ (38)

The relation expressed by Equation (38) is continuous in the neighborhood of the nominal values $e^0\left(x\right)\triangleq {\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right),{\alpha }^0\right]}^{†}$.

3) Using the definition of the adjoint operator $L^{*}\left(x,\alpha \right)$ provided in Equation (7) makes it possible to recast the left-side of Equation (38) into the following form:

$\begin{array}{l}{\left\{{\langle {\left(\begin{array}{c}{\psi }_1^{\left(1\right)}\\ {\psi }_2^{\left(1\right)}\end{array}\right)}^{†},\left(\begin{array}{cc}L\left(\alpha \right)& 0\\ 0& L^{*}\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}\delta \phi \\ \delta \psi \end{array}\right)\rangle }_{\left(1\right)}\right\}}_{\left(e^0\right)}\\ ={\left\{{\langle {\left(\begin{array}{c}\delta \phi \\ \delta \psi \end{array}\right)}^{†},\left(\begin{array}{cc}{L}^{*}\left(\alpha \right)& 0\\ 0& L\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}{\psi }_1^{\left(1\right)}\\ {\psi }_2^{\left(1\right)}\end{array}\right)\rangle }_{\left(1\right)}\right\}}_{\left(e^0\right)}\\ +{\left\{P^{\left(1\right)}\left[{\alpha }^0;{\psi }_1^{\left(1\right)}\left(x\right);{\psi }_2^{\left(1\right)}\left(x\right);\delta \phi \left(x\right),\delta \psi \left(x\right)\right]\right\}}_{\partial {\Omega }_x}.\end{array}$ (39)

where ${\left\{P^{\left(1\right)}\left[{\alpha }^0;{\psi }_1^{\left(1\right)}\left(x\right);{\psi }_2^{\left(1\right)}\left(x\right);\delta \phi \left(x\right),\delta \psi \left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ denotes the associated bilinear form evaluated on $\partial {\Omega }_x$.

4) Require the right-side of Equation (35) and the first term on the right-side of Equation (39) to represent the same functional, which is accomplished by requiring that

${\left\{\left(\begin{array}{cc}{L}^{*}\left(\alpha \right)& 0\\ 0& L\left(\alpha \right)\end{array}\right)\left(\begin{array}{c}{\psi }_1^{\left(1\right)}\left(x\right)\\ {\psi }_2^{\left(1\right)}\left(x\right)\end{array}\right)\right\}}_{\left(e^0\right)}={\left\{\left(\begin{array}{c}\frac{\partial F\left(f;e;\alpha \right)}{\partial \phi }{p}_f\left(x\right)\\ \frac{\partial F\left(f;e;\alpha \right)}{\partial \psi }{p}_f\left(x\right)\end{array}\right)\right\}}_{\left(e^0\right)},$ (40)

which implies that the following system of equations is to be satisfied in the weak (distributional) sense:

${\left\{L^{*}\left(x,\alpha \right){\psi }_1^{\left(1\right)}\left(f;x\right)\right\}}_{\left(e^0\right)}={\left\{\frac{\partial F\left(f;e;\alpha \right)}{\partial \phi }{p}_f\left(x\right)\right\}}_{\left(e^0\right)},x\in {\Omega }_x,$ (41)

${\left\{L\left(x,{\alpha }^0\right){\psi }_2^{\left(1\right)}\left(f;x\right)\right\}}_{\left(e^0\right)}={\left\{\frac{\partial F\left(f;e;\alpha \right)}{\partial \psi }{p}_f\left(x\right)\right\}}_{\left(e^0\right)},x\in {\Omega }_x,$ (42)

Since the right-sides of Equations (41) and (42) depend on the index f, it follows that the functions ${\phi }^{\left(1\right)}\left(f;x\right)$ and ${\psi }^{\left(1\right)}\left(f;x\right)$ will also depend on the index f, which means that they will be distinct for each value of $f=1,\cdots ,NF$, a fact that has been emphasized by showing the index f in the arguments of these functions.

5) Determine the adjoint boundary and/or initial conditions for the functions ${\phi }^{\left(1\right)}\left(f;x\right)$ and ${\psi }^{\left(1\right)}\left(f;x\right)$, represented in operator form as

${\left\{B_A^{\left(1\right)}\left[\phi \left(x\right);\psi \left(x\right);{\psi }_1^{\left(1\right)}\left(f,x\right),{\psi }_2^{\left(1\right)}\left(f,x\right);\alpha \right]\right\}}_{\left(e^0\right)}=0,x\in \partial {\Omega }_x,$ (43)

where the superscript “(1)” indicates “1^st-Level.” The boundary and/or initial conditions represented by Equation (43) are obtained by requiring that they satisfy the following two criteria:

a) They must be independent of unknown values of $\delta \phi$, $\delta \psi \left(x\right)$ and $\delta \alpha$.

b) The substitution of the boundary and/or initial conditions represented by Equations (21), (27) and (43) into the expression of the bilinear concomitant ${\left\{P^{\left(1\right)}\left[{\alpha }^0;{\psi }_1^{\left(1\right)}\left(x\right);{\psi }_2^{\left(1\right)}\left(x\right);\delta \phi \left(x\right),\delta \psi \left(x\right)\right]\right\}}_{\partial {\Omega }_x}$, must cause all terms containing unknown values of $\delta \phi$ and $\delta \psi \left(x\right)$ in this bilinear concomitant to vanish. The selection of the boundary conditions for the functions ${\phi }^{\left(1\right)}\left(f;x\right)$ and ${\psi }^{\left(1\right)}\left(f;x\right)$, based on the two criteria mentioned above, eliminates the appearance of the unknown values of the variations $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$, and consequently reduces the bilinear concomitant ${\left\{P^{\left(1\right)}\left[{\alpha }^0;{\psi }_1^{\left(1\right)}\left(x\right);{\psi }_2^{\left(1\right)}\left(x\right);\delta \phi \left(x\right),\delta \psi \left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ to a residual quantity that contains boundary terms involving only known values of $\delta \alpha$, the index f and the nominal parameter values ${\alpha }^0$. This residual quantity will be denoted by ${\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)$, which is linear in $\delta \alpha$ and can be therefore written in the form

${\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)={\displaystyle \underset{i=1}{\overset{N_{\alpha }}{\sum }}\frac{\partial {\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)}{\partial {\alpha }_i}\delta {\alpha }_i}.$ (44)

In general, the residual quantity ${\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)$ does not automatically vanish, although it may do so in particular instances. By considering extensions of $L\left(x,\alpha \right)$ and $L^{*}\left(x,\alpha \right)$, ${\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)$ could be forced to vanish, if necessary.

6) The system comprising Equations (41)-(43) will be called the First-Level Adjoint Sensitivity System (1^st-LASS) and the functions ${\psi }_1^{\left(1\right)}\left(f;x\right)$ and ${\psi }_2^{\left(1\right)}\left(f;x\right)$ will be called the first-level adjoint sensitivity functions. It is important to note that the 1^st-LASS is independent of the variations $\delta \alpha$ in the model and boundary parameters. It is noteworthy that the function ${\psi }_1^{\left(1\right)}\left(f;x\right)$ satisfies an “adjoint-type” operator equation, i.e., Equation (41), while the function ${\psi }_2^{\left(1\right)}\left(f;x\right)$ satisfies “forward-type” operator equation, i.e., Equation (42). The operator equations underlying the 1^st-LASS are to be solved using the nominal values $e^0$ of the model’s parameters and state functions. By continuity, however, Equations (41)-(43) are valid not only at the nominal values $e^0$, but also within a neighborhood $e$ of $e^0$.

7) Using the right-side of Equation (38) to replace the left-side of Equation (39), replacing the right-sides of Equation (40) into the right-side of Equation (39), and recalling that the boundary conditions in Equations (21), (27), and (43) have reduced the bilinear concomitant

${\left\{P^{\left(1\right)}\left[{\alpha }^0;{\psi }_1^{\left(1\right)}\left(x\right);{\psi }_2^{\left(1\right)}\left(x\right);\delta \phi \left(x\right),\delta \psi \left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ to the residual quantity ${\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)$, yields the following form for Equation (39):

$\begin{array}{l}{\left\{{\langle {\psi }_1^{\left(1\right)}\left(f,x\right),Q_L^{\left(1\right)}\left(\phi ,\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}+{\left\{{\langle {\psi }_2^{\left(1\right)}\left(f,x\right),Q_A^{\left(1\right)}\left(\psi ,\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}\\ ={\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)+{\langle \delta \phi \left(x\right),{\left[\frac{\partial F\left(f;e;\alpha \right)}{\partial \phi }{p}_f\left(x\right)\right]}_{\left(e^0\right)}\rangle }_{\phi }\\ +{\langle \delta \psi \left(x\right),{\left[\frac{\partial F\left(f;e;\alpha \right)}{\partial \psi }{p}_f\left(x\right)\right]}_{\left(e^0\right)}\rangle }_{\phi }.\end{array}$ (45)

8) The second term on the right-side of Equation (45) is actually the indirect-effect term defined in Equation (35), so that Equation (45) can be re-written as follows:

$\begin{array}{l}{\left\{\delta F\left(f;e^0;h\right)\right\}}_{ind}\triangleq {\left\{\delta F\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right);\delta \alpha \right]\right\}}_{ind}\\ ={\left\{{\langle {\psi }_1^{\left(1\right)}\left(f,x\right),Q_F^{\left(1\right)}\left(\phi ,\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}+{\left\{{\langle {\psi }_2^{\left(1\right)}\left(f,x\right),Q_A^{\left(1\right)}\left(\psi ,\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}\\ -{\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right),f=1,\cdots ,NF.\end{array}$ (46)

As the right-side of Equation (46) indicates, all unknown values of have been eliminated from the expression of the indirect effect term, ${\left\{\delta F\left(f;e^0;h\right)\right\}}_{ind}$ by having replaced the variations $\delta \phi$ and $\delta \psi \left(x\right)$ with functions that depend on the 1^st-level adjoint functions ${\phi }^{\left(1\right)}\left(f;x\right)$ and ${\psi }^{\left(1\right)}\left(f;x\right)$. The identity ${\left\{\delta F\left(f;e^0;h\right)\right\}}_{ind}\triangleq {\left\{\delta F\left[e^0;{\phi }^{\left(1\right)}\left(f;x\right);{\psi }^{\left(1\right)}\left(f;x\right)\right];\delta \alpha \right\}}_{ind}$ in Equation (46) explicitly indicates that the indirect-effect term now depends on the 1^st-level adjoint functions ${\phi }^{\left(1\right)}\left(f;x\right)$ and ${\psi }^{\left(1\right)}\left(f;x\right)$ rather than on the unknown values of the variations $\delta \phi$ and $\delta \psi \left(x\right)$.

9) Inserting the results provided in Equations (22), (28) and (44) into Equation (46) yields the following expression for the indirect-effect term:

$\begin{array}{l}{\left\{\delta F\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right];\delta \alpha \right\}}_{ind}\\ ={\displaystyle \underset{i=1}{\overset{N_{\alpha }}{\sum }}\{{\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}{\left\{{\psi }_1^{\left(1\right)}\left(f,x\right)\cdot \frac{\partial Q_F^{\left(1\right)}\left(\phi ,\alpha ;\delta \alpha \right)}{\partial {\alpha }_i}+{\psi }_2^{\left(1\right)}\left(f,x\right)\cdot \frac{Q_A^{\left(1\right)}\left(\psi ,\alpha ;\delta \alpha \right)}{\partial {\alpha }_i}\right\}}_{\left(e^0\right)}\text{d}{x}_i}}}\\ -\frac{\partial {\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)}{\partial {\alpha }_i}\}\delta {\alpha }_i,f=1,\cdots ,NF.\end{array}$ (47)

Adding the results obtained in Equation (47) to the direct-effect term expressed by Equation (36) yields the following expression for the total 1^st-order sensitivity $\delta F\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right]$ of the spectral expansion coefficient $F\left(f;\phi ,\psi ;\alpha \right)$ :

$\begin{array}{l}\delta F\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right]\\ ={\displaystyle \underset{i_1=1}{\overset{N_{\alpha }}{\sum }}{F}_{i_1}^{\left(1\right)}\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right]\delta {\alpha }_{i_1}},f=1,\cdots ,NF.\end{array}$ (48)

where the quantity $F_{i_1}^{\left(1\right)}\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right]$ denotes the 1^st-order partial sensitivities (G-derivatives) of the functional $F\left(f;\phi ,\psi ;\alpha \right)$ with respect to the model/boundary parameters ${\alpha }_i$, evaluated at the nominal values of the 1^st-level adjoint functions and model and boundary parameters, and is given by the following expression, for each $f=1,\cdots ,NF$ :

$\begin{array}{l}{F}_{i_1}^{\left(1\right)}\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right]\\ \triangleq {\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{\frac{\partial F\left(f;\alpha \right)}{\partial {\alpha }_i}{p}_f\left(x\right)\right\}}_{\left(e^0\right)}-{\left\{\frac{\partial {\hat{P}}^{\left(1\right)}\left(f;\alpha \right)}{\partial {\alpha }_i}\right\}}_{\left(e^0\right)}}}\\ +{\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}{\left\{{\psi }_1^{\left(1\right)}\left(f,x\right)\cdot \frac{\partial Q_F^{\left(1\right)}\left(\phi ,\alpha ;\delta \alpha \right)}{\partial {\alpha }_i}+{\psi }_2^{\left(1\right)}\left(f,x\right)\cdot \frac{Q_A^{\left(1\right)}\left(\phi ,\alpha ;\delta \alpha \right)}{\partial {\alpha }_i}\right\}}_{\left(e^0\right)}\text{d}{x}_i}}\\ +{\displaystyle \underset{j=1}{\overset{N_x}{\sum }}{\displaystyle \underset{\begin{array}{l}i=1\\ i\ne j\end{array}}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{F\left[f;\phi \left(\mathrm{..}{x}_j^u\left(\alpha \right)\mathrm{..}\right),\psi \left(\mathrm{..}{x}_j^u\left(\alpha \right)\mathrm{..}\right);\alpha \right]p_f\left[x_j^u\left(\alpha \right)\right]\frac{\partial x_j^u\left(\alpha \right)}{\partial {\alpha }_i}\right\}}_{\left(e^0\right)}}}}\\ -{\displaystyle \underset{j=1}{\overset{N_x}{\sum }}{\displaystyle \underset{\begin{array}{l}i=1\\ i\ne j\end{array}}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left({\alpha }^0\right)}{\overset{x_i^u\left({\alpha }^0\right)}{\int }}\text{d}{x}_i{\left\{F\left[f;\phi \left(\mathrm{..}{x}_j^{\mathcal{l}}\left(\alpha \right)\mathrm{..}\right),\psi \left(\mathrm{..}{x}_j^{\mathcal{l}}\left(\alpha \right)\mathrm{..}\right);\alpha \right]p_f\left[x_j^{\mathcal{l}}\left(\alpha \right)\right]\frac{\partial x_j^{\mathcal{l}}\left(\alpha \right)}{\partial {\alpha }_i}\right\}}_{\left(e^0\right)}}}}.\end{array}$ (49)

In particular, contributions to the sensitivities $F_{i_1}^{\left(1\right)}\left[e^0;{\psi }_1^{\left(1\right)}\left(f,x\right);{\psi }_2^{\left(1\right)}\left(f,x\right)\right]$ with respect to the specific parameters that characterize the imprecisely known model’s domain boundary can arise from the second term of the right-side of Equation (49), i.e., the boundary term ${\hat{P}}^{\left(1\right)}\left(f;\alpha ;\delta \alpha \right)$, and from the last two terms on the right-side of Equation (49) involving the boundary end-points $x_j^{\mathcal{l}}\left(\alpha \right)$ and $x_j^u\left(\alpha \right)$.

3.2. Derivation of the 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for Computing Exactly and Efficiently the 2^nd-Order Response Sensitivities to Model and Boundary Parameters

The 2^nd-order total sensitivity, denoted as ${\delta }^{\left(2\right)}F\left(e^0;h_1,h_2\right)$, of the functional $F\left(f;\phi ,\psi ;\alpha \right)$ is provided by the 2^nd-order G-differential of $F\left(f;\phi ,\psi ;\alpha \right)$ to variations $h\triangleq \left(\delta \phi ,\delta \psi ;\delta \alpha \right)$ in the model/boundary parameters and state functions around the respective nominal values $e^0=\left({\phi }^0,{\psi }^0;{\alpha }^0\right)$. The 2^nd-order total sensitivity ${\delta }^{\left(2\right)}F\left(e^0;h_1,h_2\right)$ is defined inductively as “the total first-order G-variation (or 1^st-order total G-differential) of the 1^st-order total response variation (or 1^st-order G-differential),” as follows:

${\delta }^{\left(2\right)}F\left(e^0;h_1,h_2\right)\triangleq {\left\{\frac{\text{d}}{\text{d}{\epsilon }_1}\left[\delta F\left(e^0+{\epsilon }_1{h}_1;h_2\right)\right]\right\}}_{{\epsilon }_1=0},$ (50)

for ${\epsilon }_1,{\epsilon }_2\in F$, and all (i.e., arbitrary) vectors of variations $h_i,i=1,2$ around $e^0$. The total 2^nd-order sensitivity ${\delta }^{\left(2\right)}F\left(e^0;h_1,h_2\right)$ is bilinear in the vectors of variations $h_i,i=1,2$, and will exists if and only if the 1^st-order partial sensitivities $F_{i_1}^{\left(1\right)}\left[e;{\phi }^{\left(1\right)}\left(f;x\right);{\psi }^{\left(1\right)}\left(f;x\right)\right]$ satisfy the conditions expressed in Equations (15) and (16) in a neighborhood $e$ of $e^0$. In order to proceed, it will be assumed that this is indeed the case. For the derivations to follow, it is convenient to introduce the vectors $g\triangleq {\left[e;{\phi }^{\left(1\right)}\left(f;x\right);{\psi }^{\left(1\right)}\left(f;x\right)\right]}^{†}$, comprising as components the vectors of model parameters, forward and adjoint state functions, as well as the 1^st-level adjoint functions. The nominal value of the vector $g$ will be denoted as $g^0\triangleq {\left[e^0;{\phi }^{\left(1,0\right)}\left(f;x\right);{\psi }^{\left(1,0\right)}\left(f;x\right)\right]}^{†}$, and the variations around $g^0$ will be denoted as $\delta g\triangleq g-g^0$. Using the definition provided in Equation (50), the 1^st-order G-differential $\delta F_{i_1}^{\left(1\right)}\left(g^0;\delta g\right)$ of the functional $F_{i_1}^{\left(1\right)}\left(g\right)$ defined in Equation (49) can be computed, for each $i_1=1,\cdots ,N_{\alpha }$, in the usual manner, namely:

$\begin{array}{l}\delta F_{i_1}^{\left(1\right)}\left(f;g^0;\delta g\right)\triangleq {\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{dir}+{\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}\\ \triangleq {\left\{\frac{\text{d}}{\text{d}\epsilon }{F}_{i_1}^{\left(1\right)}\left({\phi }^0+\epsilon \delta \phi ,{\psi }^0+\epsilon \delta \psi ;{\phi }^{\left(1,0\right)}+\epsilon \delta {\phi }^{\left(1\right)},{\psi }^{\left(1,0\right)}+\epsilon \delta {\psi }^{\left(1\right)};{\alpha }^0+\epsilon \delta \alpha \right)\right\}}_{\epsilon =0}\end{array}$ (51)

where ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{dir}$ denotes the “direct-effect term”, while ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}$ denotes the “indirect-effect term,” which are defined as follows:

${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{dir}\triangleq {\left\{\frac{\partial F_{i_1}^{\left(1\right)}\left[e;{\phi }^{\left(1\right)}\left(f;x\right);{\psi }^{\left(1\right)}\left(f;x\right)\right]}{\partial \alpha }\delta \alpha \right\}}_{\left(g^0\right)},$ (52)

${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}\triangleq {\left\{\frac{\partial F_{i_1}^{\left(1\right)}}{\partial \phi }\delta \phi +\frac{\partial F_{i_1}^{\left(1\right)}}{\partial \psi }\delta \psi +\frac{\partial F_{i_1}^{\left(1\right)}}{\partial {\phi }^{\left(1\right)}}\delta {\phi }^{\left(1\right)}+\frac{\partial F_{i_1}^{\left(1\right)}}{\partial {\psi }^{\left(1\right)}}\delta {\psi }^{\left(1\right)}\right\}}_{\left(g^0\right)}.$ (53)

For notational simplicity, the arguments of the functional $F_{i_1}^{\left(1\right)}\left[e;{\phi }^{\left(1\right)}\left(f;x\right);{\psi }^{\left(1\right)}\left(f;x\right)\right]$ have been omitted in Equation (53). The index f, has been kept in the argument of the direct and indirect-effect terms in order to serve as a reminder that the 2^nd-order sensitivities $\delta F_{i_1}^{\left(1\right)}\left(f;g^0;\delta g\right)$ are distinct for each $f=1,\cdots ,NF$ and each $i_1=1,\cdots ,N_{\alpha }$. Since the direct-effect term, ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{dir}$, is independent of: 1) variations $\delta \phi$ and $\delta \psi$ in the forward and adjoint state functions; and 2) variations $\delta {\phi }^{\left(1\right)}$ and $\delta {\psi }^{\left(1\right)}$ in the 1^st-level adjoint functions, itcan be computed immediately at this stage, without needing any additional large-scale computations. On the other hand, the “indirect-effect term”, ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}$, cannot be computed at this stage because the variations $\left(\delta \phi ,\delta \psi ,\delta {\phi }^{\left(1\right)},\delta {\psi }^{\left(1\right)}\right)$ are not available. The variations $\delta \phi$ and $\delta \psi$ could be determined by solving the 1^st-LFSS, which is computationally expensive. Furthermore, the variations $\delta {\phi }^{\left(1\right)}$ and $\delta {\psi }^{\left(1\right)}$ are the solutions of the system of equations obtained by taking 1^st-order G-differentials of the 1^st-LASS, i.e., Equations (41)-(43), which yields the following equations:

$\begin{array}{l}{\left\{-\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \phi }{p}_f\left(x\right)\delta \phi \left(x\right)-\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \phi }{p}_f\left(x\right)\delta \psi \left(x\right)+L^{*}\left(\alpha \right)\delta {\psi }_1^{\left(1\right)}\left(f;x\right)\right\}}_{\left(g^0\right)}\\ ={\left\{\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \alpha }{p}_f\left(x\right)-\frac{\partial \left[L^{*}\left(\alpha \right){\psi }_1^{\left(1\right)}\left(x\right)\right]}{\partial \alpha }\right\}}_{\left(g^0\right)}\delta \alpha \\ \triangleq {\left\{Q_1^{\left(2\right)}\left(f;\phi ;\psi ;{\psi }_1^{\left(1\right)};\alpha ;\delta \alpha \right)\right\}}_{\left(g^0\right)},x\in {\Omega }_x,\end{array}$ (54)

$\begin{array}{l}{\left\{-\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \psi }\delta \phi \left(x\right)-\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \psi }\delta \psi \left(x\right)+L\left(\alpha \right)\delta {\psi }_2^{\left(1\right)}\left(x\right)\right\}}_{\left(g^0\right)}\\ ={\left\{\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \alpha }{p}_f\left(x\right)-\frac{\partial \left[L\left(\alpha \right){\psi }_2^{\left(1\right)}\left(x\right)\right]}{\partial \alpha }\right\}}_{\left(g^0\right)}\delta \alpha \\ \triangleq {\left\{Q_2^{\left(2\right)}\left(f;\phi ;\psi ;{\psi }_2^{\left(1\right)};\alpha ;\delta \alpha \right)\right\}}_{\left(g^0\right)},x\in {\Omega }_x,\end{array}$ (55)

$\begin{array}{l}{\left\{\delta B_A^{\left(1\right)}\left[\phi \left(x\right);\psi \left(x\right);{\psi }_1^{\left(1\right)}\left(f,x\right),{\psi }_2^{\left(1\right)}\left(f,x\right);\alpha \right]\right\}}_{\left(g^0\right)}\\ \triangleq \{\frac{\partial B_A^{\left(1\right)}\left(\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};\alpha \right)}{\partial \phi }\delta \phi +\frac{\partial B_A^{\left(1\right)}\left(\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};\alpha \right)}{\partial \psi }\delta \psi \\ +\frac{\partial B_A^{\left(1\right)}\left(\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};\alpha \right)}{\partial {\psi }_1^{\left(1\right)}\left(f,x\right)}\delta {\psi }_1^{\left(1\right)}\left(f,x\right)\\ {+\frac{\partial B_A^{\left(1\right)}\left(\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};\alpha \right)}{\partial {\psi }_2^{\left(1\right)}\left(f,x\right)}\delta {\psi }_2^{\left(1\right)}\left(f,x\right)+\frac{\partial B_A^{\left(1\right)}\left(\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};\alpha \right)}{\partial \alpha }\delta \alpha \}}_{\left(g^0\right)}\\ =0,x\in \partial {\Omega }_x.\end{array}$ (56)

Equations (54)-(56) together with the 1^st-LFSS constitute the 2^nd-Level Forward Sensitivity System (2^nd-LFSS). For further developments, it is convenient to write Equations (20), (26), (54) and (55) in the following matrix form:

${\left\{L^{\left(2\right)}\left(f;g\right)u^{\left(2\right)}\left(x\right)\right\}}_{\left(g^0\right)}={\left\{S^{\left(2\right)}\left(f;g\right)\right\}}_{\left(g^0\right)},x\in {\Omega }_x,$ (57)

where

$\begin{array}{l}{L}^{\left(2\right)}\left(f;g;x\right)\triangleq {\left[L_{ij}^{\left(2\right)}\right]}_{\left(4\times 4\right)}\\ \triangleq \left[\begin{array}{cccc}L\left(\alpha \right)& 0& 0& 0\\ 0& L^{*}\left(\alpha \right)& 0& 0\\ -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \phi }{p}_f\left(x\right)& -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \phi }{p}_f\left(x\right)& L^{*}\left(\alpha \right)& 0\\ -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \psi }{p}_f\left(x\right)& -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \psi }{p}_f\left(x\right)& 0& L\left(\alpha \right)\end{array}\right],\end{array}$ (58)

$u^{\left(2\right)}\left(x\right)\triangleq \left[\begin{array}{c}\delta \phi \left(x\right)\\ \delta \psi \left(x\right)\\ \delta {\psi }_1^{\left(1\right)}\left(x\right)\\ \delta {\psi }_2^{\left(1\right)}\left(x\right)\end{array}\right],S^{\left(2\right)}\left(f;g\right)\triangleq \left[\begin{array}{c}{Q}_F^{\left(1\right)}\left(\phi ;\alpha ;\delta \alpha \right)\\ Q_A^{\left(1\right)}\left(\psi ;\alpha ;\delta \alpha \right)\\ Q_1^{\left(2\right)}\left(f;\phi ;\psi ;{\psi }_1^{\left(1\right)};\alpha ;\delta \alpha \right)\\ Q_2^{\left(2\right)}\left(f;\phi ;\psi ;{\psi }_2^{\left(1\right)};\alpha ;\delta \alpha \right)\end{array}\right],$ (59)

Since the source-terms of the 2^nd-LFSS depend on the parameter variations $\delta {\alpha }_i$, the 2^nd-LFSS is computationally expensive to solve. To avoid the need for solving the 2^nd-LFSS, the indirect-effect term defined in Equation (53) will be expressed in terms of a 2^nd-Level Adjoint Sensitivity System (2^nd-LASS), which will be constructed by following the general principles introduced by Cacuci [4] [5] [6], comprising the following sequence of steps:

1) Define a Hilbert space, denoted as $H^{\left(2\right)}$, having vector-valued elements of the form the $a^{\left(2\right)}\left(x\right)\triangleq {\left[a_1^{\left(2\right)}\left(x\right),a_2^{\left(2\right)}\left(x\right),a_3^{\left(2\right)}\left(x\right),a_4^{\left(2\right)}\left(x\right)\right]}^{†}\in H^{\left(2\right)}$, with components that are $N_u$ -dimensional vectors of the form $a_i^{\left(2\right)}\left(x\right)\triangleq {\left[a_{i,1}^{\left(2\right)}\left(x\right),\cdots ,a_{i,j}^{\left(2\right)}\left(x\right),\cdots ,a_{i,N_u}^{\left(2\right)}\left(x\right)\right]}^{†},i=1,2,3,4$ with square-integrable components $a_{i,j}^{\left(2\right)}\left(x\right),i=1,2,3,4;j=1,\cdots ,N_u$. The inner-product in $H^{\left(2\right)}$, denoted as ${\langle a^{\left(2\right)}\left(x\right),b^{\left(2\right)}\left(x\right)\rangle }_{\left(2\right)}$, of two vector-valued functions $a^{\left(2\right)}\left(x\right)\in H^{\left(2\right)}$ and $b^{\left(2\right)}\left(x\right)\in H^{\left(2\right)}$, is defined as follows:

$\begin{array}{c}{\langle a^{\left(2\right)}\left(x\right),b^{\left(2\right)}\left(x\right)\rangle }_{\left(2\right)}\triangleq {\displaystyle \underset{j=1}{\overset{4}{\sum }}{\langle a_k^{\left(2\right)}\left(x\right),b_k^{\left(2\right)}\left(x\right)\rangle }_{\left(\phi \right)}}\\ ={\displaystyle \underset{j=1}{\overset{4}{\sum }}{\displaystyle \underset{k=1}{\overset{N_u}{\sum }}{\displaystyle \underset{i=1}{\overset{N_x}{\prod }}{\displaystyle \underset{x_i^{\mathcal{l}}\left(\alpha \right)}{\overset{x_i^u\left(\alpha \right)}{\int }}{a}_{j,k}^{\left(1\right)}\left(x\right)b_{j,k}^{\left(1\right)}\left(x\right)\text{d}{x}_i}}}}.\end{array}$ (60)

In particular, the inner product defined in Equation (60) also holds at the nominal parameter values ${\alpha }^0$. Using the definition provided in Equation (60), construct the inner product of a vector ${\psi }_{i_1}^{\left(2\right)}\left(x\right)\triangleq {\left[{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]}^{†}\in H^{\left(2\right)}$ with Equation (57) to obtain the following relation:

${\left\{{\langle {\psi }_{i_1}^{\left(2\right)}\left(x\right),L^{\left(2\right)}\left(f;g\right)u^{\left(2\right)}\left(x\right)\rangle }_{\left(2\right)}\right\}}_{\left(g^0\right)}={\left\{{\langle {\psi }_{i_1}^{\left(2\right)}\left(x\right),S^{\left(2\right)}\left(f;g\right)\rangle }_{\left(2\right)}\right\}}_{\left(g^0\right)}.$ (61)

The relation expressed by Equation (60) is continuous in the neighborhood of the nominal values $g^0\triangleq {\left[e^0;{\phi }^{\left(1,0\right)}\left(f;x\right);{\psi }^{\left(1,0\right)}\left(f;x\right)\right]}^{†}$.

2) Using the definition of the adjoint operator in the Hilbert space $H_{\left(2\right)}$, the left-side of Equation (61) is recast as follows:

$\begin{array}{l}{\left\{{\langle {\psi }_{i_1}^{\left(2\right)}\left(x\right),L^{\left(2\right)}\left(f;g\right)u^{\left(2\right)}\left(x\right)\rangle }_{\left(2\right)}\right\}}_{\left(g^0\right)}\\ ={\left\{{\langle u^{\left(2\right)}\left(x\right),{\left[L^{\left(2\right)}\right]}^{*}\left(f;g\right){\psi }_{i_1}^{\left(2\right)}\left(x\right)\rangle }_{\left(2\right)}\right\}}_{\left(g^0\right)}\\ +{\left\{P^{\left(2\right)}\left[f;{\alpha }^0;\delta \phi ;\delta \psi ;\delta {\psi }_1^{\left(1\right)},\delta {\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]\right\}}_{\partial {\Omega }_x}\end{array}$ (62)

where ${\left\{P^{\left(2\right)}\left[f;{\alpha }^0;\delta \phi ;\delta \psi ;\delta {\psi }_1^{\left(1\right)},\delta {\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ denotes the bilinear concomitant defined on the phase-space boundary $x\in \partial {\Omega }_x$, and where ${\left[L^{\left(2\right)}\right]}^{*}$ denotes the formal adjoint of the operator $L^{\left(2\right)}$, having components $L_{ji}^{\left(2\right)*}$ defined as the transposed formal adjoints of the components of $L^{\left(2\right)}$, i.e.,

$\begin{array}{l}{\left[L^{\left(2\right)}\right]}^{*}\triangleq {\left[L_{ji}^{\left(2\right)*}\right]}_{\left(4\times 4\right)}\\ \triangleq \left[\begin{array}{cccc}{L}^{*}\left(\alpha \right)& 0& -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \phi }{p}_f\left(x\right)& -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \psi }{p}_f\left(x\right)\\ 0& L\left(\alpha \right)& -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \phi }{p}_f\left(x\right)& -\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \psi }{p}_f\left(x\right)\\ 0& 0& L\left(\alpha \right)& 0\\ 0& 0& 0& L^{*}\left(\alpha \right)\end{array}\right]\end{array}$ (63)

3) The first term on the right-side of Equation (62) is now required to represent the indirect-effect term defined in Equation (53), which is accomplished by requiring that the following relation be satisfied:

$\begin{array}{l}{\left\{{\langle u^{\left(2\right)}\left(x\right),{\left[L^{\left(2\right)}\right]}^{*}\left(f;g\right){\psi }_{i_1}^{\left(2\right)}\left(x\right)\rangle }_{\left(2\right)}\right\}}_{\left(g^0\right)}\\ ={\left\{{\left[\frac{\partial F_{i_1}^{\left(1\right)}}{\partial \phi },\frac{\partial F_{i_1}^{\left(1\right)}}{\partial \psi },\frac{\partial F_{i_1}^{\left(1\right)}}{\partial {\phi }^{\left(1\right)}},\frac{\partial F_{i_1}^{\left(1\right)}}{\partial {\psi }^{\left(1\right)}}\right]}^{†}\right\}}_{\left(g^0\right)},\text{for}f=1,\cdots ,NF;i_1=1,\cdots ,N_{\alpha }.\end{array}$ (64)

4) The boundary conditions given in Equations (56) are now implemented in Equation (62), thereby reducing by half the number of unknown boundary-values in the bilinear concomitant

${\left\{P^{\left(2\right)}\left[f;{\alpha }^0;\delta \phi ;\delta \psi ;\delta {\psi }_1^{\left(1\right)},\delta {\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]\right\}}_{\partial {\Omega }_x}$. The boundary conditions for the functions ${\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)$ are now chosen so as to eliminate the remaining unknown boundary-values of the functions $\delta \phi$, $\delta \psi$, $\delta {\psi }_1^{\left(1\right)}$ and $\delta {\psi }_2^{\left(1\right)}$ in

${\left\{P^{\left(2\right)}\left[f;{\alpha }^0;\delta \phi ;\delta \psi ;\delta {\psi }_1^{\left(1\right)},\delta {\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ while ensuring that Equation (64) is well posed. The boundary conditions thus chosen for the adjoint functions ${\psi }_1^{\left(1\right)}\left(x\right)$ and ${\psi }_2^{\left(1\right)}\left(x\right)$ can be represented in operator form as follows:

$B_A^{\left(2\right)}\left[f;\alpha ;\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)};{\psi }_{2,i_1}^{\left(2\right)};{\psi }_{3,i_1}^{\left(2\right)};{\psi }_{4,i_1}^{\left(2\right)}\right]=0,x\in \partial {\Omega }_x.$ (65)

In most cases, the above choice of boundary conditions for the 1^st-level adjoint function ${\psi }^{\left(1\right)}\left(x\right)$ will cause the bilinear concomitant ${\left\{P^{\left(2\right)}\left[f;{\alpha }^0;\delta \phi \left(x\right);\delta \psi \left(x\right);\delta {\psi }_1^{\left(1\right)}\left(x\right),\delta {\psi }_2^{\left(1\right)}\left(x\right);{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ in Equation (62) to vanish. Even when it does not vanish, however, this bilinear concomitant will be reduced to a residual quantity, which will be denoted as ${\hat{P}}^{\left(2\right)}\left(f;i_1;{\alpha }^0;\delta \alpha \right),f=1,\cdots ,NF;i_1=1,\cdots ,N_{\alpha }$, and which will contain only known values of its arguments.

5) Replacing the right-side of Equation (61) by the left-side of Equation (62), replacing the right-side of Equation (64) in the first term on the right-side of Equation (62), and recalling that the boundary conditions in Equations (56) and (65) have reduced the bilinear concomitant

${\left\{P^{\left(2\right)}\left[f;{\alpha }^0;\delta \phi ;\delta \psi ;\delta {\psi }_1^{\left(1\right)},\delta {\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right)\right]\right\}}_{\partial {\Omega }_x}$ to the residual quantity ${\hat{P}}^{\left(2\right)}\left(f;i_1;{\alpha }^0;\delta \alpha \right)$, yields the following form for Equation (62):

$\begin{array}{l}{\langle {\psi }_{i_1}^{\left(2\right)}\left(x\right),{\left\{S^{\left(2\right)}\left(f;g\right)\right\}}_{\left(g^0\right)}\rangle }_{\left(2\right)}\\ ={\hat{P}}^{\left(2\right)}\left(f;i_1;{\alpha }^0;\delta \alpha \right)+{\langle u^{\left(2\right)}\left(x\right),{\left\{{\left[\frac{\partial F_{i_1}^{\left(1\right)}}{\partial \phi },\frac{\partial F_{i_1}^{\left(1\right)}}{\partial \psi },\frac{\partial F_{i_1}^{\left(1\right)}}{\partial {\phi }^{\left(1\right)}},\frac{\partial F_{i_1}^{\left(1\right)}}{\partial {\psi }^{\left(1\right)}}\right]}^{†}\right\}}_{\left(g^0\right)}\rangle }_{\left(2\right)}.\end{array}$ (66)

6) The second term on the right-side of Equation (66) is actually the indirect-effect term, ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}$, defined in Equation (53), so that Equation (66) can be re-written in the following detailed form:

$\begin{array}{l}{\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}\\ =-{\hat{P}}^{\left(1\right)}\left(f;{\alpha }^0;\delta \alpha \right)+{\left\{{\langle {\psi }_{1,i_1}^{\left(2\right)}\left(f;x\right),Q_F^{\left(1\right)}\left(\psi ,\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}\\ +{\left\{{\langle {\psi }_{2,i_1}^{\left(2\right)}\left(f;x\right),Q_A^{\left(1\right)}\left(\psi ,\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}\\ +{\left\{{\langle {\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right),Q_1^{\left(2\right)}\left(f;\phi ;\psi ;{\psi }_1^{\left(1\right)};\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)}\\ +{\left\{{\langle {\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right),Q_2^{\left(2\right)}\left(f;\phi ;\psi ;{\psi }_2^{\left(1\right)};\alpha ;\delta \alpha \right)\rangle }_{\phi }\right\}}_{\left(e^0\right)},f=1,\cdots ,N_f;i_1=1,\cdots ,N_{\alpha }.\end{array}$ (67)

As the right-side of Equation (67) indicates, all unknown values of the variations $\delta \phi \left(x\right),\delta \psi \left(x\right),\delta {\psi }_1^{\left(1\right)}\left(x\right),\delta {\psi }_2^{\left(1\right)}\left(x\right)$ have been eliminated from the expression of the indirect effect term, ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}$, by having replaced these unknown variations by quantities that depend on the functions ${\psi }_{1,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right)$.

The system of equations comprising Equations (64) and (65) will be called the 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) for the 2^nd-level adjoint functions ${\psi }_{1,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right)$, $f=1,\cdots ,N_f$, $i=1,\cdots ,N_{\alpha }$. The explicit form of Equation (64) is

$\begin{array}{l}{L}^{*}\left(\alpha \right){\psi }_{1,i_1}^{\left(2\right)}\left(f;x\right)-\left[\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \phi }{p}_f\left(x\right)\right]{\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right)\\ -\left[\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial \psi }{p}_f\left(x\right)\right]{\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right)=\frac{\partial F_{i_1}^{\left(1\right)}\left(f;x\right)}{\partial \phi },\end{array}$ (68)

$\begin{array}{l}L\left(\alpha \right){\psi }_{2,i_1}^{\left(2\right)}\left(f;x\right)-\left[\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \phi }{p}_f\left(x\right)\right]{\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right)\\ -\left[\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial \psi }{p}_f\left(x\right)\right]{\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right)=\frac{\partial F_{i_1}^{\left(1\right)}\left(f;x\right)}{\partial \psi },\end{array}$ (69)

$L\left(\alpha \right){\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right)=\frac{\partial F_{i_1}^{\left(1\right)}\left(f;x\right)}{\partial {\phi }^{\left(1\right)}}$ (70)

$L^{*}\left(\alpha \right){\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right)=\frac{\partial F_{i_1}^{\left(1\right)}\left(f;x\right)}{\partial {\psi }^{\left(1\right)}}$ (71)

It is important to note that the 2^nd-LASS is independent of the variations $\delta \alpha$ in the model and boundary parameters. The 2^nd-level adjoint functions ${\psi }_{2,i_1}^{\left(2\right)}\left(f;x\right)$ and ${\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right)$ are the solutions of equations involving the forward operator $L\left(\alpha \right)$, while the 2^nd-level adjoint functions ${\psi }_{1,i_1}^{\left(2\right)}\left(f;x\right)$ and ${\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right)$ are the solutions of equations involving the adjoint operator $L^{*}\left(\alpha \right)$. Thus, the 2^nd-LASS is solved successively by using two “forward” and two “adjoint” computations, for each of the imprecisely known scalar model parameters. The operator equations underlying the 1^st-LASS are to be solved using the nominal values $e^0$ of the model’s parameters and state functions. By continuity, however, Equations (68)-(71) are valid not only at the nominal parameter values ${\alpha }^0$, but also within a neighborhood $\alpha$ of ${\alpha }^0$.

Replacing the expression obtained in Equation (67) for the indirect-effect term ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{ind}$ together with the expression in Equation (52) for the direct-effect term ${\left\{\delta F_{i_1}^{\left(1\right)}\left(f\right)\right\}}_{dir}$ into Equation (51) yields the following expression for the total 2^nd-order response sensitivity $\delta F_{i_1}^{\left(1\right)}\left(f;g^0;\delta g\right)$ :

$\begin{array}{l}\text{for}f=1,\cdots ,NF;i_1=1,\cdots ,N_{\alpha }:\delta F_{i_1}^{\left(1\right)}\left(f;g;\delta g\right)\\ \triangleq {\displaystyle \underset{i_2=1}{\overset{N_{\alpha }}{\sum }}{F}_{i_1{i}_2}^{\left(2\right)}\left[f;\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right);\alpha \right]\delta {\alpha }_j},\end{array}$ (72)

where $F_{i_1{i}_2}^{\left(2\right)}\left[f;\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right);\alpha \right]\equiv {\partial }^{2}F/\partial {\alpha }_{i_1}\partial {\alpha }_{i_2}$, $i_1,i_2=1,\cdots ,N_{\alpha };f=1,\cdots ,N_f$, denotes the 2^nd-order partial sensitivity of the response to the model parameters and is given by the following expression:

$\begin{array}{l}{F}_{i_1{i}_2}^{\left(2\right)}\left[f;\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right);\alpha \right]\equiv {\partial }^{2}F/\partial {\alpha }_{i_1}\partial {\alpha }_{i_2}\\ ={\langle {\psi }_{1,i_1}^{\left(2\right)}\left(x\right),\frac{\partial \left[Q_F\left(\alpha \right)-L\left(\alpha \right)\phi \right]}{\partial {\alpha }_{i_2}}\rangle }_{\phi }+{\langle {\psi }_{2,i_1}^{\left(2\right)}\left(x\right),\frac{\partial \left[Q_A\left(\alpha \right)-L^{*}\left(\alpha \right)\psi \right]}{\partial {\alpha }_{i_2}}\rangle }_{\phi }\\ +{\langle {\psi }_{3,i_1}^{\left(2\right)}\left(x\right),\frac{{\partial }^{2}F\left(f;e\right)}{\partial \phi \partial {\alpha }_{i_2}}{p}_f\left(x\right)-\frac{\partial \left[L^{*}\left(\alpha \right){\psi }_1^{\left(1\right)}\left(x\right)\right]}{\partial {\alpha }_{i_2}}\rangle }_{\phi }\\ +{\langle {\psi }_{4,i_1}^{\left(2\right)}\left(x\right),\frac{{\partial }^{2}F\left(f;e\right)}{\partial \psi \partial {\alpha }_{i_2}}{p}_f\left(x\right)-\frac{\partial \left[L\left(\alpha \right){\psi }_2^{\left(1\right)}\left(x\right)\right]}{\partial {\alpha }_{i_2}}\rangle }_{\phi }+\frac{\partial R_i^{\left(1\right)}\left(\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};\alpha \right)}{\partial {\alpha }_{i_2}}\\ -\frac{\partial {\hat{P}}^{\left(2\right)}\left[\phi ;\psi ;{\psi }_1^{\left(1\right)};{\psi }_2^{\left(1\right)};{\psi }_{1,i_1}^{\left(2\right)}\left(x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(x\right);\alpha \right]}{\partial {\alpha }_{i_2}},f=1,\cdots ,NF;i_1,i_2=1,\cdots ,N_{\alpha }.\end{array}$ (73)

All quantities appearing in Equation (73) are to be computed at the nominal values for the model and boundary parameters, although this requirement has not been indicated explicitly in Equation (73) in order to simplify the notations. In particular, contributions to the 2^nd-order response sensitivities involving boundary parameters will arise not only from the last two terms in Equation (73), i.e., from the residual term on the domain’s boundary and from the direct-effect term, but may also arise from other terms involving mixed second-order derivatives with respect to the imprecisely known boundary parameters.

Since the 2^nd-LASS is independent of parameter variations $\delta \alpha$, each of the G-differentials $\delta F_{i_1}^{\left(1\right)}$ contains one complete row (running on the index $i_2=1,\cdots ,N_{\alpha }$ ) of the mixed second-order partial sensitivities, ${\partial }^{2}F/\partial {\alpha }_{i_1}\partial {\alpha }_{i_2}$, which can be computed exactly and efficiently after solving only once the 2^nd-Level Adjoint Sensitivity System (2^nd-LASS) to obtain the 2^nd-level adjoint functions ${\psi }_{1,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{2,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{3,i_1}^{\left(2\right)}\left(f;x\right),{\psi }_{4,i_1}^{\left(2\right)}\left(f;x\right)$, for each index $i_1=1,\cdots ,N_{\alpha }$ and index $f=1,\cdots ,NF$. Thus, for each index $f=1,\cdots ,NF$, the exact computation of all of the partial second-order sensitivities, ${\partial }^{2}F/\partial {\alpha }_{i_1}\partial {\alpha }_{i_2}$ $i_1,i_2=1,\cdots ,N_{\alpha }$ requires at most $N_{\alpha }$ large-scale (adjoint) computations using the 2^nd-LASS, rather than $\left(N_{\alpha }^2/2+3N_{\alpha }/2\right)$ large-scale computations as would be required by re-computations and/or forward methods. It is also important to note that by solving the 2^nd-LASS $N_{\alpha }$ -times, the “off-diagonal” 2^nd-order mixed sensitivities ${\partial }^{2}F/\partial {\alpha }_{i_1}\partial {\alpha }_{i_2}$, $i_1\ne i_2$, will be computed twice, in two different ways (i.e., using distinct 2^nd-level adjoint functions), thereby providing an independent intrinsic (numerical) verification that the respective sensitivities are computed accurately.

Another important characteristic of using the 2^nd-LASS is the flexibility it provides for prioritizing the computation of the 2^nd-order sensitivities. The 2^nd-order mixed G-derivatives (sensitivities) corresponding to the largest relative 1^st-order response sensitivity should be computed first; the second largest relative 1^st-order response sensitivity should be considered next, and so on. Computing 2^nd-order partial sensitivities that correspond to vanishing 1^st-order sensitivities may also be of interest, since vanishing 1^st-order sensitivities may indicate critical points of the response in the phase-space of model parameters. Thus, only the 2^nd-order partial sensitivities which are deemed important should be computed; the unimportant ones can be deliberately neglected while knowing the error incurred by neglecting them.

It is noteworthy that the solving the equations underlying the 1^st-LASS and 2^nd-LASS require computational solvers for solving (i.e., “inverting”) either the forward linear operator $L\left(\alpha \right)$ or the adjoint linear operator $L^{*}\left(\alpha \right)$. Only the various source on the right-sides of the 1^st-LASS and 2^nd-LASS differ from one another. Therefore, the same computer program and “solvers” can be used for solving the 1^st-LASS and the 2^nd-LASS as were already used for solving the original system of equations. For this reason, the 2^nd-LASS was designated as a “second-level” rather than a “second-order” adjoint sensitivity system, since the 2^nd-LASS does not involve any explicit 2^nd-order G-derivatives of the operators underlying the original system, but involves the inversion of the same operators that are needed to be inverted for solving the 1^st-LASS and the original system of equations.

4. Conclusions

This work has presented the “Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2^nd-CASAM)” for the efficient and exact computation of 1^st- and 2^nd-order response sensitivities to imprecisely known parameters and domain boundary for a generic/general model for linear physical systems. The model’s response (i.e., model result of interest) is considered to be a generic function/operator that depends on the model’s forward and adjoint state functions, as well as on the imprecisely known boundary and model parameters. It has been shown that the novel 2^nd-CASAM has the following features that generalize and extend all previously published works on this topic:

1) The system response considered within the 2^nd-CASAM framework is an operator-valued response that depends on both the forward and adjoint state-functions. Functional-valued responses are subsumed as particular cases.

2) The 2^nd-CASAM framework enables the efficient and exact computation of the 1st- and 2^nd-order response sensitivities not only to uncertain model parameters but also to uncertain domain boundaries of generic linear systems. In particular, contributions to the first-order response sensitivities with respect to the specific parameters that characterize the imprecisely known model’s domain boundary can arise solely from boundary terms or directly from the response if its definition includes the domain’s boundary. On the other hand, mixed second-order response sensitivities involving boundary parameters may arise from all source terms (of the 2^nd-LASS) that involve the imprecisely known boundary parameters.

Companion works [18] [19] present illustrative applications of the 2^nd-CASAM to various responses for a paradigm model of generic particle/radiation transmission and/or evolution which has fundamental applications in many fields, including nuclear engineering (e.g., radiation detection, chemical reprocessing of spent reactor fuel, radioactive decay, etc.).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Wigner, E.P. (1945) Effect of Small Perturbations on Pile Period. Chicago Report CP-G-3048.
[2]	Dan Cacuci, G. (1981) Sensitivity Theory for Nonlinear Systems: I. Nonlinear Functional Analysis Approach. Journal of Mathematical Physics, 22, 2794-2802. https://doi.org/10.1063/1.525186
[3]	Dan Cacuci, G. (1981) Sensitivity Theory for Nonlinear Systems: II. Extensions to Additional Classes of Responses. Journal of Mathematical Physics, 22, 2803-2812. https://doi.org/10.1063/1.524870
[4]	Dan Cacuci, G. (2003) Sensitivity and Uncertainty Analysis: Theory, Volume 1. Chapman & Hall/CRC, Boca Raton.
[5]	Cacuci, D.G., Ionescu-Bujor, M. and Navon, M.I. (2005) Sensitivity and Uncertainty Analysis: Applications to Large Scale Systems, Volume 2. Chapman & Hall/CRC, Boca Raton. https://doi.org/10.1201/9780203483572
[6]	Dan Cacuci, G. (2015) Second Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) for Computing Exactly and Efficiently First- and Second-Order Sensitivities in Large-Scale Linear Systems: I. Computational Methodology. Journal of Computational Physics, 284, 687-699. https://doi.org/10.1016/j.jcp.2014.12.042
[7]	Dan Cacuci, G. (2016) Second-Order Adjoint Sensitivity Analysis Methodology (2nd-ASAM) for Large-Scale Nonlinear Systems: I. Theory. Nuclear Science and Engineering, 184, 16-30. https://doi.org/10.13182/NSE16-16
[8]	Dan Cacuci, G. (2018) The Second-Order Adjoint Sensitivity Analysis Methodology. CRC Press, Taylor & Francis Group, Boca Raton. https://doi.org/10.1201/9781315120270
[9]	Komata, M. (1977) A Generalized Perturbation Theory Applicable to Reactor Boundary Changes. Nuclear Science and Engineering, 64, 811. https://doi.org/10.13182/NSE77-A14496
[10]	Larsen, E.W. and Pomraning, G.C. (1981) Boundary Perturbation Theory. Nuclear Science and Engineering, 77, 415. https://doi.org/10.13182/NSE81-A18954
[11]	Rahnema, F. and Pomraning, G.C. (1983) Boundary Perturbation Theory for Inhomogeneous Transport Equations. Nuclear Science and Engineering, 84, 313-319. https://doi.org/10.13182/NSE83-A15451
[12]	Rahnema, F. (1996) Boundary Condition Perturbations in Transport Theory. Nuclear Science and Engineering, 124, 320. https://doi.org/10.13182/NSE96-A28581
[13]	McKinley, M.S. and Rahnema, F. (2002) High-Order Boundary Condition Perturbation Theory for the Neutron Transport Equation. Nuclear Science and Engineering, 140, 285-294. https://doi.org/10.13182/NSE02-A2261
[14]	Favorite, J.A. (2006) Variational Estimates of Internal Interface Perturbations and a New Variational Functional for Inhomogeneous Transport Problems. Nuclear Science and Engineering, 152, 180. https://doi.org/10.13182/NSE06-A2574
[15]	Favorite, J.A. and Gonzalez, E. (2017) Revisiting Boundary Perturbation Theory for Inhomogeneous Transport Problems. Nuclear Science and Engineering, 185, 445-459. https://doi.org/10.1080/00295639.2016.1277108
[16]	Dan Cacuci, G. (2019) Second-Order Sensitivities of a General Functional of the Forward and Adjoint Fluxes in a Multiplying Nuclear System with Source. Nuclear Engineering and Design, 344, 83-106. https://doi.org/10.1016/j.nucengdes.2019.01.007
[17]	Dan Cacuci, G. (2019) Second Order Adjoint Sensitivity Analysis of a General Ratio of Functionals of the Forward and Adjoint Fluxes in a Multiplying Nuclear System with Source. Annals of Nuclear Energy, 135, Article ID: 106956. https://doi.org/10.1016/j.anucene.2019.106956
[18]	Cacuci, D.G. (2020) Illustrative Application of the 2nd-Order Adjoint Sensitivity Analysis Methodology to a Paradigm Linear Evolution/Transmission Model: Point-Detector Response. AJCM, 10, 355-381.
[19]	Cacuci, D.G. (2020) Illustrative Application of the 2nd-Order Adjoint Sensitivity Analysis Methodology to a Paradigm Linear Evolution/Transmission Model: Reaction-Rate Detector Response. AJCM, 10, 382-397.