The n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (n th -CASAM-L): I. Mathematical Framework

: This work presents the mathematical framework of the n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (abbreviated as “n th -CASAM-L”), which is conceived for obtaining the exact expressions of arbitrarily-high-order (n th -order) sensitivities of a generic system response with respect to all of the parameters (including boundary and initial conditions) underlying the respective forward/adjoint systems. Since many of the most important responses for linear as used for solving the original forward and adjoint systems, thus requiring relatively minor additional software development for computing the various-order sensitivities.


Introduction
The functional derivatives of results (customarily called "responses") produced by computational models of physical systems with respect to the model's parameters are referencing, the pattern underlying the mathematical framework of the 4 th -CASAM-L is succinctly reviewed in Sections 4 and 5. The main features of the structures of the mathematical frameworks for n = 1, 2, 3, 4, are summarized in Tables 1-4. This information will be used in the next Subsection to surmise the general mathematical framework underlying the n th -CASAM-L, for any order n, to enable the computations of response sensitivities of any order, n. Table 1. First-Order Sensitivities (n = 1).

Particular Case: The 5 th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Coupled Forward/Adjoint Linear Systems (5 th -CASAM-L)
As an additional application of the general framework of the n th -CASAM, this Subsection will present the "5 th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Coupled Forward/Adjoint Linear Systems" (5 th -CASAM-L). The derivation of the 5 th -CASAM-L will be performed explicitly, by starting with the results produced by the 4 th -CASAM-L, which are presented in Section 5 and summarized in Table 4. It will be shown that the results thus obtained for the 5 th -CASAM-L coincide with the results that would be predicted by the n th -CASAM-L (Table 5) for the particular value n = 5.

Mathematical Modeling of Response-Coupled Linear Forward and Adjoint Systems
The mathematical model of a physical system comprises independent variables (e.g., space, time, etc.), dependent variables (aka "state functions"; e.g., temperature, mass, momentum, etc.) and various parameters (appearing in correlations, coordinates of physical boundaries, etc.), which are all interrelated by equations that usually represent conservation laws. The model parameters usually stem from processes that are external to the system under consideration and are seldom, if ever, known precisely. Without loss of generality, the model parameters can be considered to be real scalar quantities, having known nominal (or mean) values and, possibly, known higher-order moments or cumulants (i.e., variance/covariances, skewness, kurtosis), which are usually determined from experimental data and/or processes external to the physical system under consideration. These imprecisely known model parameters will be denoted as α 1 , . . . ,α TP , where the subscript "TP" indicates "Total (number of) Parameters". It will be convenient to consider that these parameters are components of a "vector of parameters" denoted as α (α 1 , . . . , α TP ) † ∈ R TP , where R TP denotes the TP-dimensional subset of the set of real scalars. These model parameters are considered to include imprecisely known geometrical parameters that characterize the physical system's boundaries in the phase-space of the model's independent variables. Thus, α ∈ R TP is a TP-dimensional column vector and its components are considered to include all of the uncertain model parameters, including those that may enter into defining the system's boundary in the phase-space of independent variables. For subsequent developments, it is convenient to denote matrices and vectors using capital and, respectively, lower case bold letters. The symbol " " will be used to denote "is defined as" or "is by definition equal to". Transposition will be indicated by a dagger ( †) superscript.
The model is considered to comprise TI independent variables which will be denoted as x i , i = 1, . . . , TI, and which will be considered as the components of a TI-dimensional column vector denoted as x (x 1 , . . . , x TI ) † ∈ R TI , where the sub/superscript "TI" denotes the "Total (number of) Independent variables". The vector x ∈ R TI of independent variables is considered to be defined on a phase-space domain, denoted as Ω(α), The lower boundary-point of an independent variable is denoted as λ i (α) (e.g., the inner radius of a sphere or cylinder, the lower range of an energy-variable, etc.), while the corresponding upper boundarypoint is denoted as ω i (α) (e.g., the outer radius of a sphere or cylinder, the upper range of an energy-variable, etc.). A typical example of boundary conditions that depend on imprecisely known model and geometrical parameters is provided by models based on diffusion theory where the boundaries of the physical domain are facing vacuum. For such models, the boundary conditions for the respective state variables (i.e., particle flux and/or current) are imposed not on the physical boundary but on the "extrapolated boundary" of the respective spatial domain. The "extrapolated boundary" depends both on the imprecisely known physical dimensions of the problem's domain and also on the medium's properties, i.e., atomic number densities and microscopic transport cross sections. The boundary of Ω(α), which will be denoted as ∂Ω(α), comprises the set of all of the endpoints λ i (α), ω i (α), i = 1, . . . , TI, of the respective intervals on which the components of A linear physical system is generally modeled by a system of coupled linear operatorequations which can be generally represented as follows: where ϕ(x) [ϕ 1 (x), . . . , ϕ TD (x)] † is a TD-dimensional column vector of dependent variables, where the sub/superscript "TD" denotes the "Total (number of) Dependent variables". The functions ϕ i (x), i = 1, . . . , TD, denote the system's "dependent variables" (also called "state functions"). The matrix L(x; α) L ij (x; α) , i, j = 1, . . . , TD, has dimensions TD × TD. The components L ij (x; α) are operators that act linearly on the dependent variables ϕ j (x) and also depend (in general, nonlinearly) on the uncertain parameters α. The components of the TD-dimensional column vector q ϕ (x; α) , where the subscript "ϕ" indicates sources associated with the "forward" system of equations, are also nonlinear functions of α. Since the right-side of Equation (63) may contain distributions, the equality in this equation is considered to hold in the weak (i.e., "distributional") sense. Similarly, all of the equalities that involve differential equations in this work will be considered to hold in the distributional sense. When L(x; α) contains differential operators, a set of boundary and/or initial conditions which define the domain of L(x; α) must also be given. Since the complete mathematical model is considered to be linear in ϕ(x), the boundary and/or initial conditions needed to define the domain of L(x; α) must also be linear in ϕ(x). Such linear boundary and/or initial conditions are represented in the following operator form: In Equation (64), the components B ij (x; α); i = 1, . . . , N B ; j = 1, . . . , TD, of the N B × TD-dimensional matrix-operator B ϕ (x; α) are operators that act linearly on ϕ(x) and nonlinearly on α; the quantity N B denotes the total number of boundary and initial conditions. The operator c ϕ (x; α) vector comprising components that are operators which in general act nonlinearly on α. The subscript "ϕ" in Equation (64) indicates boundary conditions associated with the forward state function ϕ(x). In this work, capital bold letters will be used to denote matrices while lower case bold letter will be used to denote vectors; the components of these matrices and vectors may be operators rather than just functions. Physical problems modeled by linear systems and/or operators are naturally defined in Hilbert spaces. The dependent variables ϕ i (x), i = 1, . . . , TD, for the physical system represented by Equations (63) and (64) are considered to be square-integrable functions of the independent variables and are considered to belong to a Hilbert space which will be denoted as H 0 , where the subscript "zero" denotes "zeroth-level" or "original". Higherlevel Hilbert spaces, which will be denoted as H 1 , H 2 , etc., will also be introduced and used in this work. The Hilbert space H 0 is considered to be endowed with the following inner product, denoted as ϕ(x), ψ(x) 0 , between two elements ϕ(x) ∈ H 0 and ψ(x) ∈ H 0 : The dot in Equation (65) indicates the "scalar product of two vectors", which is defined as follows: while the product-notation [ ] dx i in Equation (65) denotes the respective multiple integrals. The linear operator L(x; α) admits an adjoint operator, denoted as L * (x; α), which is defined through the following relation for a vector ψ(x) ∈ H 0 : In Equation (67), the formal adjoint operator L * (x; α) is the TD × TD matrix comprising elements L * ji (x; α) obtained by transposing the formal adjoints of the forward operators L ij (x; α), i.e., Hence, the system adjoint to the linear system represented by Equations (63) and (64) can generally be represented as follows: When L(x; α) comprises differential operators, the operations (e.g., integration by parts) that implement the transition from the left-side to the right side of Equation (67) give rise to boundary terms which are collectively called the "bilinear concomitant". The domain of L * (x; α) is determined by selecting adjoint boundary and/or initial conditions so as to ensure that the bilinear concomitant vanishes when the selected adjoint boundary conditions are implemented together with the forward boundary conditions given in Equation (64). The adjoint boundary conditions thus selected are represented in operator form by Equation (70), where the subscript "ψ" indicates adjoint boundary and/or initial conditions associated with the adjoint state function ψ(x). In most practical situations, the Hilbert space H 0 is self-dual.
The system of equations comprising Equations (63), (64), (69) and (70) will be called the 1 st -Level Forward System (1 st -LFS), since it is the starting point for the computation of the model's response, and its solution will be needed for the subsequent computation of the first-order sensitivities of the model's response with respect to the model parameters. The 1 st -LFS can be written compactly in block-matrix form as follows: where the following definitions were used: The nominal solution of Equations (63) and (64), which will be denoted as ϕ 0 (x), and the nominal solution of Equations (69) and (70), which will be denoted as ψ 0 (x), are obtained by solving these equations at the nominal (or mean) values of the model parameter α 0 . The superscript "zero" will henceforth be used to denote "nominal" (or, equivalently, "expected" or "mean" values). Thus, the vectors ϕ 0 (x), ψ 0 (x) and α 0 satisfy the following equations: or, equivalently, In Equations (75) and (76), the notation { } α 0 has been used to indicate that the quantity within the brackets is to be evaluated at the nominal values of the parameters and state functions. This simplified notation is justified by the fact that when the parameters take on their nominal values, it implicitly means that the corresponding state functions also take on their corresponding nominal values. This simplified notation will be used throughout this work.
The relationship shown in Equation (67), which is the basis for defining the adjoint operator, also provides the following fundamental "reciprocity-like" relation between the sources of the forward and the adjoint equations, cf. Equations (63) and (69), respectively: The functional on the right-side of Equation (77) represents a "detector response", i.e., reaction-rate between the particles and the medium represented by q ψ (x; α) which is equivalent to the "number of counts" of particles incident on a detector of particles that "measures" the particle flux ϕ(x). Thus, the source term q ψ (x; α) q ψ,1 (x; α), . . . ., q ψ,TD (x; α) † in Equation (77) is usually associated with the "result of interest" to be measured and/or computed, which is customarily called the system's "response". In particular, if , which means that, in such a case, the right-side of Equation (77) provides the value of the dependent variable (particle flux, temperature, velocity, etc.) at the point in phase-space where the respective measurement is performed. The relation given in Equation (77) indicates that for linear physical systems, the system's response may depend not only on the model's state-functions and on the system parameters, but may also depend on the adjoint function The system's response will be denoted as R[ϕ(x), ψ(x); α] and, in the most general case, is a nonlinear operator acting on the model's forward and adjoint state functions, as well as on imprecisely known parameters, both directly and indirectly through the state functions. The nominal value of the response, R ϕ 0 (x), ψ 0 (x); α 0 , is determined by using the nominal parameter values α 0 , the nominal value ϕ 0 (x) of the forward state function and the nominal value ψ 0 (x) of the adjoint function.
In general, a function-valued (operator) response R op [ϕ(x), ψ(x); α] can be represented by a spectral expansion in multidimensional orthogonal polynomials or Fourier series of the form: where the quantities P m i (x i ), i = 1, . . . , TI, denote the corresponding spectral functions (e.g., orthogonal polynomials or Fourier exponential/trigonometric functions) and where the spectral Fourier coefficients c m 1 ...m TI are defined as follows: The coefficients c m 1 ...m TI can themselves be considered as system responses since the spectral polynomials P m i (x i ) are perfectly well known while the expansion coefficients will contain all of the dependencies of the respective response on the imprecisely known model and response parameters. This way, the sensitivity analysis of operator-valued responses can be reduced to the sensitivity analysis of scalar-valued responses. A measurement of a physical quantity can be represented as a response R p ϕ x p , ψ x p ; α located at a specific point, x p , in phase-space, having the following form: where δ x − x p denotes the multidimensional Dirac-delta functional. Responses that occur in many fields (e.g., optimization, control, model verification, data assimilation, model validation, model calibration, predictive modeling) are Lagrangians having the following functional form: A particular form of the response defined in Equation (81) is the Schwinger "normalization-free Lagrangian", which takes on the following form [1,2]: The system responses defined in Equation (78) through (82) indicate that their sensitivities (i.e., functional derivatives) with respect to the model parameters can be studied by considering a generic scalar-valued response, which will be denoted as R[ϕ(x), ψ(x); α], of the following form: where S [ϕ(x), ψ(x); α] denotes a suitably Gateaux-(G-) differentiable function of the indicated arguments. In general, S [ϕ(x), ψ(x); α] is nonlinear in ϕ(x), ψ(x), and α, and the components of α are considered to also include parameters that may specifically appear only in the definition of the response under consideration (but which might not appear in the definition of the model). Thus, the (physical) "system" is defined in this work to comprise both the system's computational model and the system's response.
The generic response defined in Equation (83) provides the basis for constructing any other responses of specific interest and will therefore be used for the generic "n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems", which will be abbreviated as "n th -CASAM-L" (where "n" is a finite number that indicates any desired, arbitrarily-high, order) and which will be developed in the remainder of this work. Note that the generic response defined in Equation (83) is, in general, a nonlinear function of all of its arguments, i.e., S [ϕ(x), ψ(x); α] is nonlinear in ϕ(x), ψ(x), and α. . . , TD in the forward and, respectively, adjoint state functions. In turn, the variations δα, δϕ, and δψ will cause a response variation R ϕ 0 + δϕ; ψ 0 + δψ; α 0 + δα around the nominal response value The first-order G-variation δ f (η; h) is an operator defined on the same domain as F(η), and has the same range as F(η). The G-variation δ f (η; h) satisfies the relation [∆(εh)]/ε = 0. The existence of the Gvariation δ f (η; h) does not guarantee its numerical computability. Numerical methods (e.g., Newton's method and variants thereof) most often require that δ f (η; h) be linear in the variations h in a neighborhood η 0 + εh around η 0 . The necessary and sufficient conditions for the G-differential δ f (η; h) of a nonlinear operator f (η) to be linear in the variations h in a neighborhood (η + εh) around η are as follows: (i). f (η) satisfies a weak Lipschitz condition at η 0 ; i.e., (ii). For two arbitrary vectors of variations h 1 and h 2 , the operator f (η) satisfies the relation In practice, it is not necessary to investigate if f (η) satisfies the conditions presented in Equations (85) and (86) since it is usually evident if the right-side of Equation (84) is linear (or not) in the variations h. A total variation δ f (η; h) which is linear in h is called the "total differential of f (u)" and is usually denoted as D f (η; h). In this case, the partial variation δ f (η 1 , . . . , η K ; h i ) is called the "partial differential δ f (η 1 , . . . , η K ; h i ) of F(η) with respect to η i ". Thus, the total differential D f (η; h) is linear in the variations h i and the following representation holds: where the quantities ∂ f (η 1 , . . . , η K )/∂η i denote the partial derivatives of f (η) with respect to its arguments η i . It will henceforth be assumed that all of the operators considered in this work satisfy the conditions presented in Equations (85) and (86), and therefore admit partial G-derivatives such that the representation shown in Equation (87) exists.
Since this work ultimately aims at deriving the explicit expressions of the n th -order (i.e., arbitrarily-high order) sensitivities of the response R[ϕ(x), ψ(x); α] with respect to the model parameters, the proliferation of indices, superscripts and subscripts is unavoidable. To avoid the use of "subscripted-subscripts" as much as possible, the subscripts which will be used for denoting the sensitivities (i.e., partial G-derivatives) of the model response with respect to the model parameters will have the following form: (i) j 1 = 1, . . . , TP, for firstorder sensitivities; (ii) j 1 ; j 2 = 1, . . . , TP, for 2 nd -order sensitivities; (iii) j 1 ; j 2 ; j 3 = 1, . . . , TP, for 3 rd -order sensitivities, and so on. Thus, the indices j 1 ; j 2 ; . . . ; j n = 1, . . . , TP will be used for the n th -order sensitivities.

∂[L(α)ϕ]/∂α, respectively.
Considering the definitions provided in Equations (91), (94) and (95), it follows from Equation (92) that the state functions ϕ(x) and ψ(x) are implicit functions of the model parameters α; therefore, the variations of the variations δϕ(x) and δψ(x) can be considered to represent the following total differentials: Using the results provided in Equations (91), (94)−(97) into Equations (92) and (93), and identifying the expressions that multiply each parameter variation δα i yields the following system of equations: with the functions ∂ϕ/∂α j 1 and ∂ψ/∂α j 1 subject to the corresponding boundary conditions derived from Equation (93). It is evident from Equation (98) that the determinations of the functions ∂ϕ/∂α j 1 and ∂ψ/∂α j 1 by solving the 1 st -LVSS would require TP large-scale computations, which would be prohibitively expensive in terms of computational resources.
The alternative path for determining the indirect-effect term {[δR(ϕ, ψ, α; δϕ, δψ)] α 0 } ind defined in Equation (89), which avoids the need to compute the functions ∂ϕ/∂α j 1 and ∂ψ/∂α j 1 , is provided by the first-order adjoint sensitivity analysis methodology which was conceived by Cacuci [7,8], which will be presented in this subsection in its extended framework, including the computation of response sensitivities to parameters underlying boundary and interface conditions. This extended framework is called the "Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward and Adjoint Linear Systems," and is abbreviated as 1 st -CASAM-L. The starting point for the 1 st -CASAM-L is the 1 st -LVSS, cf. Equations (92) and (93), which will be written in the following block-matrix form: where The solution δu (1) (x) of the "1 st -Level Variational Sensitivity System" (1 st -LVSS) will be called the "1 st -level variational sensitivity function, as highlighted by the superscript "(1)".

The system of equations comprising Equation (105) together with the boundary conditions represented by Equation (107) constitutes the 1 st -Level Adjoint Sensitivity System
(1 st -LASS). The solution a (1) 2 (x) † of the 1 st -LASS is called the 1 st -level adjoint function. The 1 st -LASS is called "first-level" (as opposed to "first-order") because it does not contain any differential or functional-derivatives, but its solution, a (1) (x), will be used below to compute the first-order sensitivities of the response with respect to the model parameters. This terminology will be also used in the sequel, when deriving the expressions for the higher-order sensitivities. It follows from Equations (103) and (105) that the indirect-effect term {[δR(ϕ, ψ, α; δϕ, δψ)] α 0 } ind can be expressed in terms of the 1 st -level adjoint function 2 (x) † as follows: As indicated by the last equality in Equation (108), the variations δϕ and δψ have been eliminated from the original expression of the indirect-effect term, which now instead depends on the 1 st -level adjoint function a (1)  2 (x). It is very important to note that the 1st-LASS is independent of parameter variations δα. Hence, the 1 st -LASS needs to be solved only once (as opposed to the 1st-LVSS, which would need to be solved anew for each parameter variation) to determine the 1 stlevel adjoint function a (1) 2 (x) † . Subsequently, the "indirect-effect term" {[δR(ϕ, ψ, α; δϕ, δψ)] α 0 } ind is computed efficiently and exactly by simply performing the integrations over the adjoint function a (1) 2 (x) † , as indicated on the right-side of Equation (108).
The total 1 st -order sensitivity of the response R[ϕ(x), ψ(x); α] to the model parameters is obtained by adding the expressions of the direct-effect term defined in Equation (90) and indirect-effect term as obtained in Equation (108), which yields the following expression: As indicated in the last equality in Equation (109), the expression of the first-order Gvariation δR u (1) (x); a (1) (x); α; δα no longer depends on the variations δϕ and δψ, which are expensive to compute, but instead depends on the 1 st -level adjoint function a (1) , which is as inexpensive to compute as the original state function(s). In particular, the expression in Equation (109) also reveals that the sensitivities of the response R[ϕ(x), ψ(x); α] to parameters that characterize the system's boundary and/or internal interfaces can arise both from the direct-effect and indirect-effect terms. It also follows from Equation (109) that the total 1 st -order G-variation δR u (1) (x); a (1) (x); α; δα can be expressed in terms of the 1 st -level adjoint functions as follows: where the quantities ∂R u (1) (x); a (1) (x); α /∂α j 1 α 0 represent the 1 st -order response sensitivities with respect to the model parameters, evaluated at the nominal values of the state functions and model parameters. For the subsequent computation of higher-order responses sensitivities, it is convenient to represent the 1 st -order response sensitivities as follows: In particular, if the residual bilinear concomitant is non-zero, the functions S (1) j1; u (1) (x); a (1) (x); α would contain suitably defined Dirac delta-functionals for expressing the respective non-zero boundary terms as volume-integrals over the phase-space of the independent variables. Dirac-delta functionals would also be needed to represent, within S (1) j1; u (1) (x); a (1) (x); α , the terms containing the derivatives of the boundary end-points with respect to the model and/or response parameters. As indicated by Equation (111), the computation of the 1 st -order response sensitivities to all model parameters requires two large-scale computations (which is just as many as solving the original system) namely: one large-scale computation for computing the first-level adjoint state function a (1) 1 (x) and one large-scale computation for computing the first-level adjoint state function a (1) 2 (x).

The 2 nd -CASAM-L: Summary
For each value of the index j 1 = 1, . . . , TP, the first-order sensitivity R (1) j 1 ; u (1) (x); a (1) (x); α will be considered to play the role of a "model response." The starting point for the development of the 2 nd -CASAM-L is provided by the expressions of the sensitivities provided in terms of the 1 st -level adjoint functions, namely Equation (111), which indicates that the computation of the first-order sensitivities R (1) j 1 ; u (1) (x); a (1) (x); α using the 1 st -CASAM-L requires the determination of the functions u (1) (x) and a (1) (x). In view of the developments underlying the 1 st -CASAM-L, the functions u (1) (x) and a (1) (x) are the solutions of the following system of equations, written in block-matrix form: where (114) The system of equations represented by Equations (112) and (113) will be called the 2 nd -Level Forward System, (2 nd -LFS) since it plays the same role as the First-Level (original) Forward System (1 st -LFS) has played for the computation of the first-order response sensitivities. In terms of the function u (2) (x), each 1 st -order sensitivity can be compactly written, for each j 1 = 1, . . . , TP, as follows: There are TP(TP + 1)/2! distinct 2 nd -order sensitivities of the response with respect to the model and response parameters. The 2 nd -order total sensitivity of the model response is obtained by applying the definition of the 1 st -order G-differential to Equation (111), which yields the following result, for each j 1 = 1, . . . , TP: where the direct-effect term δR (1) j 1 ; u (2) (x); α; δα α 0 dir is defined as follows: while the indirect-effect term δR (1) j 1 ; u (2) (x); α; δu (2) is defined as follows: (1) ; a (1) ; α ∂u (2) δu (2) where As indicated by the subscript "zero," the expression in Equation (118) is to be computed at nominal parameter and state function values.
The results provided in Equation (130), (131) and (128) are employed in Equation (118) to obtain the following expression for the indirect-effect term in terms of the 2 nd -level adjoint functions a (2) (j 1 ; x): Inserting the expressions that define the vector q (2) u (2) ; α; δα from Equation (122) into Equation (133) and adding the resulting expression for the indirect-effect term with the expression of the direct-effect term given in Equation (117) yields the following expression for the total second-order G-differential of the response R[ϕ(x), ψ(x); α]: where R (2) j 2 , j 1 ; u (2) (x); a (2) (j 1 ; x); α denotes the 2 nd -order partial sensitivity of the response with respect to the model parameters, evaluated at the nominal parameter values α 0 , and has the following expression for j 1 = 1, . . . , TP; j 2 = 1, . . . , j 1 : Since the 2 nd -LASS is independent of parameter variations δα, the exact computation of all of the partial second-order sensitivities R (2) j 2 , j 1 ; u (2) (x); a (2) (j 1 ; x); α requires at most TP large-scale (adjoint) computations using the 2 nd -LASS, rather than O TP 2 largescale computations as would be required by forward methods. It is important to note that by solving the 2 nd -LASS TP-times, the "off-diagonal" 2 nd -order mixed sensitivities ∂ 2 R/∂α j 1 ∂α j 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 1 st -and 2 ndorder response sensitivities are computed accurately. The information provided by the 1 st -order sensitivities usually indicates which 2 nd -order sensitivities are important and which could be neglected. Therefore, it is useful to prioritize the computation of the 2 nd -order sensitivities by using the rankings of the relative magnitudes of the 1st-order sensitivities as a "priority indicator": the larger the magnitude of the relative 1 st -order sensitivity, the higher the priority for computing the corresponding 2 nd -order sensitivities. Also, since vanishing 1 st -order sensitivities may indicate critical points of the response in the phase-space of model parameters, it is also of interest to compute the 2 nd -order sensitivities that correspond to vanishing 1 st -order sensitivities. In practice, only those 2 nd -order partial sensitivities which are deemed important would need to be computed.
Dirac delta-functionals in Equation (135) may need to be used for expressing the non-zero residual terms in the respective residual bilinear concomitant and/or the terms containing derivatives with respect to the lower-and upper-boundary points, so that the expression of the partial second-order sensitivities R (2) j 2 , j 1 ; u (2) (x); a (2) (j 1 ; x); α can be written in the following form, in preparation for computing the 3 rd -order sensitivities: (2) (j 1 ; x); α . (136)

Discussion and Conclusions
This work has presented the "n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Linear Systems" (abbreviation: n th -CASAM-L), where "n" is a finite number that indicates any desired, arbitrarily-high, order. The n th -CASAM-L enables the efficient computation of the exact expressions of the n th -order functional derivatives ("sensitivities") of a general system response, which depends on both the forward and adjoint state functions, with respect to all of the parameters (including boundary and initial conditions) underlying the respective forward and adjoint systems. Since nonlinear operators do not admit adjoint operators (only linearized versions of nonlinear operators can admit a bona-fide adjoint operator), responses that simultaneously depend on forward and adjoint functions can arise only in conjunction with linear systems. Particularly important model responses that involve both the forward and adjoint functions are various forms of Lagrangian functionals, the Roussopoulos functional for computing reaction rates, the Raleigh quotient for computing eigenvalues and/or separation constants when solving partial differential equations, the Schwinger functional for "normalization-free" solutions, and many others, (e.g., [1,2]). These functionals play a fundamental role in various optimization and control procedures, derivation of numerical methods for solving (differential, integral, integro-differential) equations, etc. The sensitivity analysis of responses that simultaneously involve both forward and adjoint state functions makes it necessary to treat linear models/systems in their own right, rather than treating them as particular cases of nonlinear systems (in which the responses can depend only on the forward functions).
This work has shown that the mathematical framework underlying the n th -CASAM-L is set in linearly increasing higher-dimensional Hilbert spaces, as opposed to exponentially increasing parameter-dimensional spaces. In particular, for a scalar-valued valued response associated with a nonlinear model comprising TP parameters, the 1 st -CASAM-L requires 1 additional large-scale adjoint computation (as opposed to TP large-scale computations, as required by other methods) for computing exactly all of the 1 st -order response sensitivities. All of the (mixed) 2 nd -order sensitivities are computed exactly by the 2 nd -CASAM-L in at most TP computations, as opposed to TP(TP + 1)/2 computations required by all other methods, and so on. For every lower-order sensitivity of interest, the n th -CASAM-L computes the "TP next-higher-order" sensitivities in one adjoint computation performed in a linearly increasing higher-dimensional Hilbert space. Very importantly, the n th -CASAM-L computes the higher-level adjoint functions using the same forward and adjoint solvers (i.e., computer codes) as used for solving the original forward and adjoint systems, thus requiring relatively minor additional software development for computing the variousorder sensitivities.
The n th -CASAM presented in this work is the only practically implementable methodology for obtaining the exact expressions (i.e., free of methodologically-introduced approximations) of arbitrarily-high order sensitivities (functional derivatives) of model responses to model parameters, for coupled forward/adjoint linear systems. By enabling the practical computation of any arbitrarily-order response sensitivities to model parameters, the n th -CASAM-L makes it possible to compare the relative values of the sensitivities of var-ious order, in order to assess which sensitivities are important and which may actually be neglected, thus enabling future investigations of the convergence of the (multivariate) Taylor series expansion of the response in terms of parameter variations, as well as investigating the actual validity of expressions that are derived from Taylor-expansion of the response (e.g., response variances/covariance, skewness, kurtosis, etc.) as a function of the model's parameters. The larger the number of model parameters, the more efficient the C-ASAM-L becomes for computing arbitrarily high-order response sensitivities. The n th -CASAM-L also enables the direct derivation, on paper, of the expression of a specific high-order sensitivity of interest, which can subsequently be computed directly; such a direct computation is not possible with any statistical method.
The n th -CASAM-L presented in this work provides a fundamentally important step in the quest to overcome the "curse of dimensionality" in sensitivity analysis, uncertainty quantification and predictive modeling. Ongoing work aims at developing the "n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems" (abbreviation: n th -CASAM-N) for the practical, efficient, and exact computation of arbitrarily-high order sensitivities of responses to model parameters for nonlinear systems. The n th -CASAM-L, together with the n th -CASAM-N, are expected to revolutionize all of the fields of activities which require response sensitivities, including the fields of uncertainty quantification, model validation, optimization, data assimilation, model calibration, sensor fusion, reduced-order modeling, inverse problems, and predictive modeling.