The nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (nth-CASAM-L): II. Illustrative Application

This work illustrates the application of the nth-order comprehensive adjoint sensitivity analysis methodology for response-coupled forward/adjoint linear systems (abbreviated as “nth-CASAM-L”) to a paradigm model that describes the transmission of particles (neutrons and/or photons) through homogenized materials, as encountered in radiation protection and shielding. The first-, second-, and third-order sensitivities of responses that depend on both the forward and adjoint particle fluxes are obtained exactly, in closed-form, underscoring the principles and methodology underlying the nth-CASAM-L. The results presented in this work underscore the fundamentally important role of the nth-CASAM-L in the quest to overcome the “curse of dimensionality” in sensitivity analysis, uncertainty quantification and predictive modeling.


Introduction
The n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (abbreviated as "n th -CASAM-L"), which is presented in the accompanying work [1], enables the most efficient computation of exactly obtained 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 −hence the qualifier "comprehensive") underlying the respective forward/adjoint systems. The application of the n th -CASAM-L is illustrated in this work by considering paradigm model which describes the transmission of particles produced by a distributed source through a shield which surrounds the source. Such models are of interest in radiation and/or particle protection and shielding [2][3][4][5]. The particular model chosen in this work is sufficiently simple to admit closed-form exact expressions for the response sensitivities of any order yet contains sufficiently many imprecisely known parameters to pose extreme, if not insurmountable, challenges to conventional statistical or finite-difference methods. The mathematical equations underlying this illustrative model are presented in Section 2, including the most important types of responses for such models, namely point-detector responses, particle leakage responses, reaction rate responses and "contributon-response fluxes" [2][3][4][5]. Sections 3-5 illustrate the application of the n th -CASAM-L (for n = 1,2,3) to the paradigm model, highlighting the efficient computation of exact expressions for the 1 st , 2 nd -and 3 rd -order response sensitivities to model parameters, including imprecisely known internal (interfaces) and external boundaries of physical domains/materials. The discussion presented in Section 6 concludes this work, emphasizing the salient points underlying the computation of arbitrarily-high order of sensitivities to imprecisely known model parameters, interfaces and domain boundaries.

Paradigm Model for Illustrating the Application of the n th -CASAM-L: Particle Transmission through a Medium with Imprecisely Known Properties, Internal and External Boundaries
The paradigm model chosen to illustrate the application of the n th -CASAM-L describes the transmission of particles produced by a distributed source through a shield which surrounds it. For simplicity, the geometry is chosen to be one-dimensional. A source of monoenergetic particles (e.g., gamma ray photons), denoted as q, is distributed within a one-dimensional homogenized inner slab of thickness 2b 1 [cm]. The various quantities that characterize this inner slab will be denoted using the subscript "s" to indicate "source" material properties. Thus, the inner slab is considered to comprise M s distinct nuclides, each nuclide being characterized by its total microscopic cross interaction cross section with the source particles, denoted as σ s,i , and by its atomic number density, denoted as N s,i , for i = 1, . . . , M s . The inner slab is shielded on each of its outer surfaces by an adjacent outer slab of thickness b 2 [cm], neither of which contain sources. Each of the outer slabs is made of a homogenized material comprising M distinct nuclides, each nuclide being characterized by its total microscopic cross section, denoted as σ i , and its atomic number density, denoted as N i , for i = 1, . . . , M. All of the number densities and cross sections are considered to be imprecisely known. The thicknesses of the slabs are also considered to be subject to manufacturing tolerances; hence, the internal boundaries (i.e., internal interface between the slabs) as well as the outer boundaries are also imprecisely known. The result of interest (i.e., model response) is the dose of uncollided particles at the outer boundaries. This paradigm model captures the essential physical processes involved in the transmission (or evolution) of the uncollided flux of particle through materials yet will be shown to admit closed-form exact solutions for all sensitivities, in order to illustrate exactly the impact of the response sensitivities to the model parameters (including the physical material boundaries), enabling their importance ranking in affecting changes in the response.
Since the chosen physical model is symmetrical with respect to the midplane of the inner slab, it is convenient to choose the origin of the one-dimensional coordinate system, which will de denoted as the "z-direction", to be at the inner slab's midplane. Consequently, the transport equation governing the uncollided angular flux of monoenergetic photons, denoted as u s (z, ω), through the inner slab in the positive z-direction has the following form: where ω denotes the cosine of the angle between the gamma ray's direction and the z-axis internal sources while µ s denotes the interaction coefficient of photons with the inner slab's homogenized material and is defined as follows: Due to the model's symmetry when choosing the origin of the coordinate system at the inner slab's midplane, the appropriate boundary condition for Equation (1) is: The solution of Equations (1) and (3) is the following constant function: Due to the model's symmetry, the expression obtained in Equation (4) remains valid for all particle directions, i.e., for −1 ≤ ω ≤ 1.
The photon transport equation governing the uncollided angular flux of monoenergetic photons, denoted as ϕ(z, ω), through the outer slab in particle directions 0 ≤ ω ≤ 1 has the following form: where the interaction coefficient µ is defined as follows: The uncollided flux is continuous across the interface at z = b 1 : For subsequent verification of the expressions to be obtained for the response sensitivities to parameters, the closed-form expression of the uncollided flux ϕ(z, ω), which is the solution of Equations (5) and (7), is provided below: The linear particle transport model described by Equations (5) and (7) is naturally set in a Hilbert space which is endowed with the following inner-product, denoted as ϕ(z, ω), ψ(z, ω) 0 , between two square integrable functions ϕ(z, ω) and ψ(z, ω): ϕ(z, ω), ψ(z, ω) 0 The Hibert endowed with the inner product defined in Equation (9) will be denoted as H 0 , where the subscript "zero" denotes "zeroth-level" or "original". Higher-level Hilbert spaces, which will be denoted as H 1 , H 2 , etc., will also be introduced and used in this work.
In the Hibert space H 0 , the adjoint model corresponding to the original forward model represented by Equations (5) and (7) is readily constructed by forming using the inner product defined in Equation (9) to obtain the following relation for the adjoint function ψ(z, ω): It follows from the relation in Equation (10) that the adjoint model for the adjoint function ψ(z, ω) is as follows: where q * (z, ω) is a source-term that is usually determined by the model response under consideration, while the boundary condition provided in Equation (12) has been chosen in order to eliminate the appearance of the unknown function ϕ(b 2 , ω) from the bilinear concomitant represented by the second term on the right-side of Equation (10).

Point-Detector Response
Doses of uncollided particles, including from point-doses, particle leakage and reaction rates, are of interest in many shielding applications [2][3][4][5]. In particular, the result of a measurement of the particle flux at any phase-space point, denoted as (z d , ω d ), within or on the surface of the medium, b 1 ≤ z d ≤ b 2 , 0 < ω d < 1, will be denoted as R(ϕ; α) and is represented the following mathematical expression: where the phase-space coordinates (z d , ω d ) of the detector's location are also imprecisely known, and where the detector's effective macroscopic cross section is considered to be an imprecisely known parameter denoted as µ d and defined as follows: where M d denotes the number of distinct nuclides comprising the detector's response function, and where σ d,i denotes the total microscopic interaction cross section of the respective nuclide with the source particles, while N d,i denotes the respective nuclide's atomic number density, for i = 1, . . . , M d . The uncertain parameters that describe the properties (microscopic cross sections, atomic number densities, source, slab boundaries.) characterizing the two materials and the detector, will be generically denoted as α i , i = 1, . . . , TP, where TP denotes the "total number of model parameters". These parameters will be considered to be the components of the (column) "vector of imprecisely known model parameters" defined as follows In this work, the dagger " †" denotes "transposition" and all vectors are column vectors except if the contrary is explicitly stated.
All of the components of the vector of parameters α are considered to be imprecisely known (i.e., uncertain). It is assumed that only the parameters' nominal values, which will be denoted using the superscript "0" (i.e., α 0 i , q 0 , z 0 d , a 0 , etc.), are known/available. The nominal values of the model parameters are considered to be components of the column vector of "nominal parameter values", denoted as α 0 and defined as follows: The parameters which appear as components of the vector defined in Equation (15) are the model's fundamental parameters. A variation in any of these fundamental parameters is independent of any variation in any of the other fundamental parameters. Not all of the fundamental parameters necessarily appear in explicitly in the model's equations or in definition of the model's response. A model and/or response may depend explicitly on derived parameters or correlations, which are functions of the fundamental parameters. For example, in the particle transport model under consideration, the interaction coefficients appear explicitly in the definition of the model's response and in the particle conservation equations that underly the paradigm model. Conversely, the nuclide number densities and/or microscopic cross sections (which are independent/fundamental parameters) do not appear explicitly in either the definition of the response or in the equations that underly the model. Parameters such as the interaction coefficients could be called derived parameters and they are important because they may act like independent parameters in the definition of the model's underlying equations and/or response, in which case it suffices to determine the sensitivities of the response to these derived parameters using the model. Subsequently, the sensitivities of the response to the underlying fundamental parameters can be computed separately, by simple differentiations. For example, in the paradigm model under consideration, after having determined the sensitivities of the form ∂R/∂µ of the detector response with respect to the interaction coefficient using the model's equation, the sensitivities to the underlying fundamental parameters (i.e., corresponding nuclide number densities and microscopic cross sections) can be determined by using the chain rule without involving the model's equations, since ∂R/∂σ i = [∂R/∂µ][∂µ/∂σ i ] = N i [∂R/∂µ], and so on. This observation may enable a significant reduction in the number of computations needed to determine the sensitivities of the response to the fundamental parameters. Notably, such simplifications occur only in linear systems since the parameters that characterize such systems cannot be functions of the dependent variable (e.g., the interaction coefficients are independent of the particle flux). In nonlinear systems, this simplification is not possible when the derived parameters depend on the state function (i.e., dependent variable). Obvious examples are heat conduction models: when the conductivity is independent of the medium's temperature (dependent variable), the respective model is linear, and the conductivity can be considered either a primary parameter (if it does not stem from some correlation) or a derived parameters (if it stems from a temperature-independent correlation). In a nonlinear conduction model, the heat conduction coefficient would depend on the temperature to be determined and could not be used as either a fundamental or a derived parameter. On the other hand, for linear systems, such as the paradigm particle transport model under consideration, the advantage of using the derived parameters will become evident in the process of determining the response sensitivities, which will be performed below.
In view of the above discussion, it is convenient to refer to the model parameters which appear explicitly in the model's equations, boundary/initial conditions, and in the model's response as primary parameters, since the response sensitivities will be computed initially with respect to these parameters. The primary parameters could be either fundamental parameters or derived parameters, as is also the case for the illustrative paradigm particle transport model considered in this work. After obtaining the response sensitivities to the model's primary parameters, if any of these parameters are "derived parameters", the sensitivities of the response to the primary parameters which are comprised in the respective "derived primary parameter" can be computed efficiently and exactly without needing to use the model's complete set of equations, thus gaining considerable accuracy and savings of computational resources. These issues will be illustrated by the derivations which will be performed in the remainder of this work.
In view of the above discussion, the interaction coefficients will be denoted as µ(α), µ s (α), and µ d (α), even though they only depend on some, but not all, of the components of α. This short-hand notation conveniently simplifies the list of the actual arguments (i.e., the corresponding nuclide number densities and microscopic cross sections) of the respective interaction coefficients.
Inserting the expression obtained in Equation (8) into Equation (13) yields the following exact, closed-form, expression for the model's response: In practice, the closed-form solution provided in Equation (8) is unavailable. Instead, the nominal value ϕ 0 (z, ω) is obtained by solving Equations (5) and (7) using the nominal parameters values represented by Equation (16). The exact, closed-form expression of the model response provided in Equation (17) is also unavailable in practice, so the value of the response must be computed by evaluating numerically the expression provided in Equation (13) using the nominal value ϕ 0 (z, ω) of the model's state function together with the nominal parameter values α 0 .
Notably, the response R(z d , ω d ; α) can also be represented in terms of an adjoint function ψ(z, ω) by using the relation provided in Equation (10); the resulting expression is: where the adjoint function ψ(z, ω) is the solution of the following adjoint system: The adjoint system comprising Equations (19) and (20) indicates that the adjoint function ψ(z, ω) plays the role of a Green's function for the response R(z d , ω d ; α). As highlighted by the expression on the right-side of Equation (18), the adjoint function ψ(z, ω) also plays the role of an "importance function" in weighing the particle transmitted from the source q to the response R(z d , ω d ; α). Notably, the model response R(z d , ω d ; α) itself plays the role of a Green's function in slab geometry, in that it can be integrated over its position coordinates the phase-space coordinates (z d , ω d ) to obtain other responses of interest involving uncollided doses from distributed sources. This property will be shown in Sections 2.2 and 2.3, below.

Particle Leakage Response
Another response of significant interest is the leakage of particles through the outer slab's surface at z = b 2 . Such a response will be denoted as R L (ϕ; α) and is represented mathematically by the following expression: Replacing the expression of ϕ(z, ω) obtained in Equation (8) into Equation (21) yields the following expression for the response R L (ϕ; α): where the exponential-integral function is defined as follows: The expression of the leakage response obtained in Equation (22) can be obtained by setting z d = b 2 in the expression of the point-detector response R(z d , ω d ; α) and by considering the point ω d in Equation (17) to be a phase-space variable (rather than a fixedpoint). Thus, integrating Equation (17) over ω d yields the expression in Equation (22), i.e.,

Reaction Rate Response
Interaction/reaction rates of particles within the medium they travel through are also of significant interest. In particular, the interaction/reaction rate of uncollided particles within the entire slab will be denoted as R r (ϕ; α) and is represented mathematically by the following expression: Replacing the expression of ϕ(z, ω) obtained in Equation (8) into Equation (25) yields the following expression for the response R r (ϕ; α): Just as in the case of the leakage response, cf. Equation (24), the reaction-rate response R r (ϕ; α) can also be obtained by considering the phase-space points (z d , ω d ) in the expression of the point-detector response R(z d , ω d ; α) as being variable, and integrating over (z d , ω d ) to obtain the following result: The results obtained in Equations (24) and (27) underscore the role as "Green's Function" (or "kernel") played by the point-detector response R(z d , ω d ; α).

Contributon-Response Flux
Bilinear functionals of the forward and adjoint particle fluxes of the form O(x)ϕ(x)ψ(x)dx, where x denotes the vector of independent variable, ϕ(x) and ψ(x) denote the forward and adjoint particle fluxes, respectively, and O(x) denotes a function that does not depend on the fluxes, occur in most Lagrangian, Raleigh, Roussopoulos and Schwinger functionals (to mention just the most prominent such functionals), which are used in many practical applications [2][3][4][5]. Perhaps the simplest example of a response that depends simultaneously on both the forward and the adjoint particle flux is the socalled "contributon-response flux", which has the form δ(x − x 0 )ϕ(x)ψ(x)dx and which is encountered in "channel theory" or "contributon theory", with applications in particle shielding analysis [4,5]. The contributon-response flux retains the essential features which characterize the sensitivity analysis of such responses.
For the illustrative paradigm model of uncollided particle transmission through the slab b 1 ≤ z ≤ b 2 governed by Equations (5) and (7), the "contributon-response flux" at a location z is defined as follows: where the function ψ(z, ω) is the solution of the following adjoint particle transport model: The adjoint source parameter q * characterizes particles that reach a detector within or on the outside surface of the slab. The adjoint particle transport model comprising Equations (29) and (30) has the following closed-form solution: The model response considered in this section is the "contributon-response flux" at a location z 0 within the slab b 1 ≤ z 0 ≤ b 2 , which is defined as follows: The closed-form expression of ρ(z 0 ) is readily obtained by inserting the expressions of ϕ(z, ω) and ψ(z, ω) from Equations (8) and (31), respectively, to obtain the following result:

Application of the 1 st -CASAM-L to Compute First-Order Response Sensitivities to Imprecisely Known Parameters
Consider arbitrary parameter variations δα around the parameters' nominal values, α 0 , where the vector δα is defined as follows: Such parameter variations will cause variations in the model's response, both directly and also indirectly, through variations they cause, via the model's underlying equations, in the model's state functions (i.e., dependent variables). The 1 st -order sensitivities of the model responses considered in Section 2 will be determined in this section by applying the general principled underlying the First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1 st -CASAM-L).

Point-Detector Response
The total first-order sensitivity of the response defined in Equation (13) to the model's primary parameters (including interface and boundary locations) is proved by the Gateaux-(G-) differential {δR(ϕ; α; δϕ; δα)} α 0 of the response R(ϕ; α) for arbitrary variations (δϕ, δα) around the nominal parameter and state functions values, which is defined as follows: where the "direct-effect" term depends only on variations δα in the parameters, is defined as follows: and where the "indirect-effect" term depends only on variations δϕ(z, ω) in the dependent variable (state function) and is defined as follows: The direct-effect term in Equation (36) can be computed immediately at this stage, since the parameter variations are known, i.e.,: . (38) On the contrary, the "indirect-effect" terms in Equation (37) cannot be computed at this stage since the 1 st -order variation, δϕ(z, ω), is not available at this stage, but could be determined by solving the 1 st -Level Variational Sensitivity System (1 st -LVSS) obtained by G-differentiating Equations (5) and (7) at the nominal parameter values, to obtain: Carrying out the differentiations with respect to ε in Equations (39) and (40) yields the following 1 st -LVSS to be satisfied by the 1 st -order variation δϕ(z, ω): The quantity {dϕ/dz} z=b 0 1 , which appears in Equation (42) can be evaluated by using either Equation (5) or Equation (8), to obtain the following boundary condition for the variation δϕ(z, ω): In principle, it is possible to solve the 1 st -LVSS for each parameter variation, δα i , i = 1, . . . , TP, as follows: (i) For each parameter variation δα i , i = 1, . . . , TP, the 1 st -LVSS is equivalent to the following system of equations: (ii) If the parameter variations δα i , i = 1, . . . , TP, are independent of each other, Equations (44) and (45) can be solved by setting δα i = δ ij , i, j = 1, . . . , TP, where δ ij represents the Kronecker delta-functional, which is defined as follows: It is evident that solving Equations (41) and (42) would require at least TP large-scale computations (i.e., at least one large-scale computation for each parameter variation) to obtain the 1 st -order variation δϕ(z, ω) for every parameter variation. Performing so many large-scale computations is impractical for large-scale systems involving many parameters. The alternative to solving repeatedly the 1 st -LVSS for every possible variation in the system's imprecisely known parameters is to use the First-Level Adjoint Sensitivity System (1 st -LASS), which is constructed by applying the following sequence of steps:

1.
Consider that the function δϕ(z, ω) belongs to a Hilbert space, which will be denoted as This Hilbert space is endowed with an inner product of two functions ϕ(z, ω) ∈ H 1 and ψ(z, ω) ∈ H 1 , which is denoted as ϕ(z, ω), ψ(z, ω) 1 and is defined as follows: 2.
Using the definition provided in Equation (46), construct the inner product of a function ψ 1 (z, ω) ∈ H u with Equation (41) to obtain 3. Integrate by parts, over the independent variable z, the left side of Equation (47) to obtain

4.
Require the last term on the right side of Equation (48) to represent the indirect-effect term defined in Equation (37) by imposing the following relationship: The relation in Equation (49) implies that the following equation holds in the weak sense at the nominal parameter values:

5.
Complete the definition of the function ψ 1 (z, ω) by requiring that the unknown value of the function δϕ(z, ω) be eliminated from appearing on the right side of Equation (48). This requirement is met by imposing the following condition be satisfied by the function ψ 1 (z, ω) on the model's outer boundary: Together, Equations (50) and (51) constitute a well-posed system, called the First-Level Adjoint Sensitivity System (1 st -LASS), for determining the function ψ 1 (z, ω) ∈ H 1 . The function ψ 1 (z, ω) ∈ H u is called the first-level adjoint sensitivity function. Since the 1 st -LASS does not depend on the parameter variations, it needs to be solved only once in order to obtain the 1 st -level adjoint sensitivity function ψ 1 (z, ω). For the paradigm model under consideration, the 1 st -LASS given by Equations (50) and (51) can be solved exactly to obtain the following closed-form expression for the 1 st -level adjoint function ψ 1 (z, ω): where H(z − z d ) denotes the Heaviside functional defined as Collecting the results obtained in Equations (47) through (51) leads to the following expression for the indirect-effect term: The appearance of the 1 st -level adjoint function ψ 1 (z, ω) in the list of arguments of the indirect-effect term {δR(ϕ; ψ 1 ; α; δα)} ind in Equation (54) emphasizes the fact the appearance of the function δϕ(z, ω) which depends implicitly on parameter variations has been eliminated. Instead, the indirect-effect term is expressed in terms of the 1 st -level adjoint function ψ 1 (z, ω), which does not depend on the model parameter variations. Hence, solving the 1 st -LASS to obtain the function ψ 1 (z, ω), which requires a single largescale computation comparable to solving the original equations underlying the model, suffices to determine subsequently all of the partial response sensitivities included in the indirect-effect term.
Adding the results obtained in Equations (54) and (36) yields the complete expression for the total 1 st -order sensitivity {δR(ϕ; α; δϕ; δα)} α 0 . The specific expression of each 1 st -order partial sensitivity of the response R(ϕ; α) to each uncertain parameter is obtained by identifying the expression that multiplies the respective parameter variation, which yields the following results: If the detector is placed inside the medium, the sensitivity of the response R(ϕ; α) to the boundary parameter b 2 is zero, as would be expected, since the response would not depend on b 2 in such a situation. On the other hand, if the detector lies on the outer surface, at z d = b 2 , then the response sensitivity to b 2 is provided by the expression in Equation (60), which would be evaluated at z d = b 2 , i.e.,: The following conclusions can be drawn based on the results obtained in Equations (55)-(62): (i) The indirect-effect term provides the complete 1 st -order partial sensitivities of the response R(ϕ; α) with respect to the parameters µ(α), q and µ s (α), which are parameters that occur solely in the equations underlying the paradigm model. (ii) The direct-effect term provides the complete 1 st -order partial sensitivities of the response R(ϕ; α) with respect to the parameters µ d (α), z d , and ω d , which are parameters that occur solely in the definition of the model's response. (iii) Both the direct and the indirect effect terms may, in general, contribute to the partial sensitivities of the response R(ϕ; α) with respect to the interface/boundary parameters b 1 and b 2 . In this particular case, the indirect-effect provides the entire contribution for the response sensitivity with respect to the boundary parameter b 1 , while the direct-effect term provides the entire contribution for the response sensitivity with respect to the boundary parameter b 2 , when the detector is placed on the outer boundary.
Notably, obtaining all of the 1 st -order partial response sensitivities to the primary model parameters has necessitated a single "large-scale" computation for determining the 1 st -level adjoint function ψ 1 (z, ω). This is in contradistinction with the use of the 1 st -LVSS, which would require at least TP large-scale computations to obtain the corresponding TP first-order sensitivities. Furthermore, the expressions obtained for these sensitivities were exact, as opposed to approximate, as would have been the case if these sensitivities would have been computed by, e.g., finite-difference or statistical procedures. Also, many more large-scale computations would have been necessary to compute these 1 st -order partial response sensitivities by any other (finite-difference or statistical) procedure. Thus, the 1 st -CASAM-L is the most efficient computational method for obtaining the exact expressions of the 1 st -order partial response sensitivities to the primary model parameters.
The primary model parameters q, z d , ω d , b 1 , b 2 are fundamental parameters. On the other hand, the primary model parameters µ s (α), µ(α) and µ d (α) are derived parameters, being themselves functions of the respective nuclide number densities and cross sections as fundamental parameters. Hence, the response sensitivities to the corresponding fundamental parameters can now be determined by taking into account the respective functional dependencies of the primary model parameters µ s (α), µ(α) and µ d (α) on the correspond-ing nuclide number densities and cross sections. Thus, the response sensitivities to the nuclide number densities and cross sections included in µ s (α) are obtained as follows: The response sensitivities to the remaining nuclide number densities and cross sections are determined similarly, to obtain the following expressions: In practice, the integrals involving the 1 st -level adjoint function ψ 1 (z, ω) in Equations (55)-(68) would be performed numerically, after having obtained ψ 1 (z, ω) by solving numerically the 1 st -LASS (which would be the only large-scale computation required to obtain all of the TP first-order sensitivities.
It is noteworthy that the relative sensitivities of the response with respect to the nuclide densities have the same values as the relative sensitivities of the response with respect to the nuclide densities, i.e., The fact that so many relative sensitivities have the same values, as indicated in Equations (69)-(71), respectively, render quasi-useless the statistical methods that attempt to rank the sensitivities by their relative values, because such statistical methods break down in situations when many relative sensitivities have equal values. The only reliable method for computing response sensitivities accurately (and also most efficiently, from a computational standpoint) in such situations is the 1 st -CASAM-L. For the paradigm model considered here, the closed form expressions of ϕ(z, ω) and ψ 1 (z, ω) are available (by design, of course), so the closed form expressions of the 1 st -order sensitivities are obtained by inserting the expression of ϕ(z, ω) from Equation (8) and the expression of the 1 st -level adjoint function ψ 1 (z, ω) from Equation (52) into Equations (55)-(62), and performing the respective integrations. The results of these operations are the following closed-form expressions for the partial first-order sensitivities of the response R(ϕ; α) with respect to the primary parameters: The response sensitivities with respect to the nuclide number densities and cross sections included in µ s (α), µ(α) and µ d (α), respectively, have the following closed-form expressions: The closed-form expressions obtained in Equations (72)-(85) can be verified directly by differentiating the closed-form expression of the response provided in Equation (17) with respect to each of the model parameters.

Particle Leakage Response
The 1 st -order total sensitivity of R L (ϕ; α) is provided by the first-order G-differential δR L (ϕ; α; δϕ; δα) α 0 of R L (ϕ; α) at the nominal parameter values, which is determined by applying the definition of the G-differential to Equation (21). This operation yields the following expression: where the direct-effect term {δR L (ϕ; α; δα) α 0 } dir depends only on the parameter variations δα and is defined as follows: and where the indirect-effect term {δR L (ϕ; α; δϕ) α 0 } ind depends only on the variations δϕ and is defined as follows: Following the principles outlined in Section 3.1, the need for computing the function δϕ(z, ω) is circumvented by expressing the indirect-effect term {δR L (ϕ; α; δϕ) α 0 } ind defined in Equation (88) in terms of a 1 st -level adjoint function, which will be denoted as ψ L (z, ω), and which is the solution of a 1 st -Level Adjoint Sensitivity System (1 st -LASS) constructed by performing the same sequence of operations as indicated in Equations (46)-(48) but replacing the right-side of Equation (49) by the right-side of Equation (88). Performing this sequence of operations leads to the following 1 st -LASS for the 1 st -level adjoint function ψ L (z, ω): In terms of the 1 st -level adjoint function ψ L (z, ω), the indirect-effect term {δR L (ϕ; α; δϕ) α 0 } ind will have the same formal expression as in Equation (54) except that the adjoint function ψ 1 (z, ω) is replaced by the adjoint function ψ L (z, ω), i.e., The total 1 st -order G-differential {δR L (ϕ; ψ L ; α; δα)} is obtained by adding the expression for the indirect-effect term obtained in Equation (91) with the expression of the direct-effect term obtained in Equation (87). Identifying in the resulting expression for {δR L (ϕ; ψ L ; α; δα)} the quantities that multiply the individual parameter variations yields the following expressions for the partial sensitivities of the response R L (ϕ; α) with respect to the various parameters: The sensitivities provided in Equations (92)-(94) stem solely from the indirect-effect term given by Equation (91). On the other hand, the direct-effect term defined in Equation (87) gives rise to the following sensitivities: Solving Equations (89) and (90) yields the following expression for the 1 st -level adjoint function ψ L (z, ω): where H(z − b 2 ) denotes the Heaviside functional. Replacing the expressions obtained in Equations (98) and (8) into Equations (92)-(97) and performing the respective integrations yields the following closed-form expressions for the 1 st -order partial sensitivities of the response R L (ϕ; α) with respect to the model parameters: The expressions obtained in Equations (99)-(104) may be verified by differentiating directly (with respect to the various parameters) the expression of R L (ϕ; α) provided in Equation (22), while also considering the following relations: As expected, the expressions of the sensitivities of R L (ϕ; α) with respect to the model parameters can be obtained directly from the corresponding sensitivities of the pointdetector response R(ϕ; α) provided in Equations (72)-(79), as follows: (i) set z d = b 2 in Equations (72)-(79); (ii) consider that the parameter ω d plays the role of a "phase-space variable" and integrate over ω d the resulting expressions. The specific results of these operations are as follows: The explicit indication that the expressions in Equations (106)-(111) are to be evaluated at the nominal parameter and state function values has been omitted, in order to simplify the respective notation.

Reaction Rate Response
The 1 st -order total sensitivity of R r (ϕ; α) is obtained by determining the first-order Gdifferential δR r (ϕ; α; δϕ; δα) α 0 of R r (ϕ; α) at the nominal parameter values. By definition, the G-differential of Equation (25) is obtained as follows: where (114) The need for computing the function δϕ(z, ω) is circumvented by expressing the indirect-effect term {δR L (ϕ; α; δϕ) α 0 } ind defined in Equation (114) in terms of a 1 st -level adjoint function, which will be denoted as ψ r (z, ω), and which is the solution of a 1 st -Level Adjoint Sensitivity System (1 st -LASS) obtained by applying the same principles and sequence of steps as in Sections 3.1 and 3.2. Performing a sequence of operations similar to the sequence followed to obtain Equations (46)-(48) but replacing the right-side of Equation (49) by the right-side of Equation (114) leads to the following 1 st -LASS for the 1 st -level adjoint function ψ r (z, ω): In terms of the 1 st -level adjoint function ψ r (z, ω), the indirect-effect term {δR L (ϕ; α; δϕ) α 0 } ind will have the same formal expression as in Equation (54) except that the adjoint function ψ 1 (z, ω) is replaced by the adjoint function ψ r (z, ω), i.e., The total 1 st -order G-differential {δR r (ϕ; ψ L ; α; δα)} is obtained by adding the expression for the indirect-effect term obtained in Equation (117) with the expression of the direct-effect term obtained in Equation (113). Identifying in the resulting expression for {δR r (ϕ; ψ L ; α; δα)} the quantities that multiply the individual parameter variations yields the following expressions for the partial sensitivities of the response R r (ϕ; α) with respect to the various parameters: The sensitivities provided in Equations (118)-(123) stem solely from the indirect-effect term given by Equation (117) and have expressions that are formally identical to the corresponding expressions of the respective sensitivities of the point-detector response R(ϕ; α) and leakage response R L (ϕ; α), except that the corresponding 1 st -level adjoint functions differ from one another.
On the other hand, the sensitivities ∂R r (ϕ; α)/∂µ d (α) and ∂R r (ϕ; α)/∂b 2 stems solely from the direct-effect term defined in Equation (113), namely: Both the direct-effect and indirect effect terms contribute to the sensitivity ∂R r (ϕ; α)/∂b 1 , which has the following expression: Solving Equations (115) and (116) yields the following expression for the 1 st -level adjoint function ψ r (z, ω): Replacing the expressions obtained in Equations (124) and (8) into Equations (118)-(123) and performing the respective integrations yields the following closed-form expressions for the 1 st -order partial sensitivities of the response R r (ϕ; α) with respect to the model parameters: The explicit indication that the expressions in Equations (125)-(130) are to be evaluated at the nominal parameter and state function values has been omitted, in order to simplify the respective notation.
As expected, the expressions of the sensitivities of R r (ϕ; α) with respect to the model parameters can be obtained directly by integrating the corresponding sensitivities of the point-detector response R(ϕ; α) provided in Equations (72)-(79) over z d and ω d , which are considered for this purpose to be independent phase-space variables (as done when using Green's functions).

Contributon-Response
The model response considered in this section is the "contributon-response flux" at a location z 0 within the slab b 1 ≤ z 0 ≤ b 2 , which was defined (28) and depends on a "vector of fundamental model parameters" α, which, in this case, has the following components to be considered for determining the sensitivities of ρ(z 0 ): The contributon-response ρ(z 0 ) depends on the following primary model parameters: q, q * , z 0 , b 1 , b 2 , µ s (α), µ(α). The interaction coefficients µ s (α) and µ(α) are "derived primary parameters" that depend, as before, on the respective number densities and microscopic cross sections, which remain "fundamental parameters".
The 1 st -order total sensitivity of ρ(ϕ, ψ; α) is obtained by determining the first-order G-differential δρ(ϕ, ψ; α; δϕ; δψ; δα) α 0 of ρ(ϕ, ψ; α) at the nominal parameter values. By definition, the G-differential of Equation (32) is obtained as follows: and The function δϕ(z, ω), which appears in the first term on the right side of Equation (134), is the solution of Equations (41) and (42). Moreover, the function δψ(z, ω), which appears in the second term on the right side of Equation (134), is the solution of the equations obtained by G-differentiating Equations (29) and (30). Applying the definition of the G-differential to Equations (29) and (30) yields the following system of equations: The 1 st -Level Variational Sensitivity System (1 st -LVSS) to be solved for obtaining the functions δϕ(z, ω) and δψ(z, ω) comprises the boundary conditions provided in Equations (43) and (136), together with Equations (41) and (135), which are written in the following matrix-form: where The need for solving the 1 st -LVSS can be circumvented by replacing the appearance of the function δu (1) (z, ω) in the indirect effect term δρ ϕ, ψ; α; δu (1) α 0 ind with a which is independent of parameter variations. This 1 st -level adjoint function is the solution of a 1 st -Level Adjoint Sensitivity System (1 st -LASS) which is constructing by introducing a Hilbert space, denoted as H 1 , which comprises square-integrable functions vector-valued elements of the form η (1) (z, ω) η 2 (z, ω) † , and which is endowed with an inner product between two elements, η (1) (z, ω) ∈ H 1 , ξ (1) (z, ω) ∈ H 1 , denoted as η (1) (x), ξ (1) (x) 1 and defined as follows: In the Hilbert H 1 , form the inner product of Equation (137) with a yet undefined vector-valued function a (1) 2 (z, ω) † ∈ H 1 to obtain the following relation: Integrating by parts the terms containing derivatives with respect to z in Equation (141) and using the boundary conditions provided in Equations (43) and (136) yields the following relation: The first two terms on the left-side of Equation (142) are required to represent the indirect-effect term {δρ(ϕ, ψ; α; δϕ; δψ) α 0 } ind defined in Equation (134) by imposing the following relations: ω da The unknown quantities δϕ(b 2 ,ω) and δψ(b 1 ,ω), which appear in the last two terms on the right-side of Equation (142), are eliminating by choosing the following boundary conditions for the 1 st -level adjoint functions a (1) 1 (z, ω) and a (1) 2 (z, ω): The system of equations comprising Equations (143)-(146) constitutes the 1 st -Level Adjoint Sensitivity System (1 st -LASS). The solution a (1) 2 (z, ω) † 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 functionalderivatives, but its solution, a (1) (z, ω), 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.
Using the relations underlying the 1 st -LASS together with boundary conditions provided in Equations (43) and (136) in Equation (142) yields the following expression for the indirect-effect term {δρ(ϕ, ψ; α; δϕ; δψ) α 0 } ind in terms of the 1 st -level adjoint function The last equality in Equation (147) highlights the fact that the no longer depends on the function δu (1) (z, ω), but depends on the adjoint 1 st -level adjoint function a (1) 2 (z, ω) † , which is independent of parameter variations. By identifying the expressions multiplying the individual parameter variations in Equation (147), it follows that the indirect-effect term {δρ(ϕ, ψ; α; δϕ; δψ) α 0 } ind contributes the expressions of the following 1 st -order sensitivities: It follows from Equation (133) that the direct-effect term contributes the expressions of the following 1 st -order sensitivity: The 1 st -LASS needs to be solved only once since it is independent of parameter variations (as opposed to the 1 st -LVSS, which would need to be solved anew for each parameter variation). Subsequently, the sensitivities of the contributon-response flux ρ(ϕ,ψ;α) to the primary model parameters can be computed efficiently and exactly by simply performing the integrations over the adjoint function a (1) 2 (z, ω) † , as indicated by the expressions provided in Equations (148)-(153). The sensitivities of ρ(ϕ,ψ;α) to the respective nuclide number densities and microscopic cross sections can be obtained directly from the sensitivities ∂ρ(ϕ,ψ;α)/∂µ(α) and ∂ρ(ϕ,ψ;α)/∂µ s (α), respectively, and have the following expressions: The closed-form expressions of the solutions a 1 (z, ω) and a 2 (z, ω) of the 1 st -LASS are obtained as follows: The sensitivities expressed by Equations (148)-(154) can now be obtained by inserting the expressions of ϕ(z,ω), ψ(z,ω), a 1 (z, ω) and a 2 (z, ω) into these equations. Using the expressions of provided in Equations (8), (31), (159) and (160), respectively, into Equations (148)-(154) and performing the respective integrations over the independent variables leads to the following expressions for the 1 st -order sensitivities of the "contributon-response flux" to the model's parameters: It is noteworthy that, even though the contributon-response flux ρ(z 0 ) vanishes at the external boundary z 0 = b 2 , i.e., ρ(z 0 = b 2 ) = 0, not all of the first-order sensitivities of ρ(z 0 ) with respect to the model parameters and boundaries vanish. The following results are obtained by setting z 0 = b 2 in Equations (161)-(167): All of the first-order sensitivities of ρ(ϕ,ψ;α) have finite values for finite non-zero parameter values, except for the sensitivity ∂ρ(ϕ,ψ;α)/∂z 0 , which becomes unbounded at the inner interface, as z 0 → b 1 . The sensitivity ∂ρ(ϕ,ψ;α)/∂b 2 is nonzero for all finite non-zero parameter values.

Application of the 2 nd -CASAM-L to Compute Second-Order Response Sensitivities to Imprecisely Known Parameters
The 2 nd -order response sensitivities to the model parameters will be obtained by applying the general principles presented in [1]. In essence, each of the first-order response sensitivities will be considered a new "model response" and the principles underlying the 1 st -CASAM-L will be applied to determine the "1 st -order sensitivities of the 1 st -order sensitivities", which by definition are the sought after 2 nd -order sensitivities. The expressions obtained in Equations (148)-(154) indicate the following characteristics: (i) The sensitivity ∂ρ(ϕ,ψ;α)/∂z 0 depends on the product ϕ(z,ω) ψ(z,ω) of the original forward and adjoint state functions, just as the contributon-response flux ρ(ϕ,ψ;α) does. Hence, the determination of the second-order sensitivities derived from ∂ρ(ϕ,ψ;α)/∂z 0 will involve the exact same steps as those that were involved in determining the first-order sensitivities of ρ(ϕ,ψ;α).   2 (z, ω) of the 1 st -level adjoint function and also on the original forward and adjoint functions ϕ(z,ω) and ψ(z,ω). Hence, the determination of the 2 nd -order sensitivities stemming from ∂ρ(ϕ,ψ;α)/∂µ(α) will require the complete sequence of steps underlying the 2 nd -CASAM-L, requiring the determination of a four-component vector-valued 2 nd -level adjoint function, as will be shown in Section 4.1.

Determination of the Second-Order Sensitivities of the Form
The explicit determination of the 2 nd -order sensitivities of the form ∂ 2 ρ(ϕ, ψ; α)/∂α i ∂b 2 , which involve the imprecisely known domain boundary at z = b 2 , will be presented in this section. These 2 nd -order sensitivities will be determined by considering that the expression provided in Equation (153) for the 1 st -order sensitivity ∂ρ(ϕ, ψ; α)/∂b 2 plays the role of a "model response". The expressions of the 2 nd -order sensitivities ∂ 2 ρ(ϕ, ψ; α)/∂α i ∂b 2 , i = 1, . . . , TP, will thus be determined by taking the first-order G-differential of the 1 st -order sensitivity ∂ρ(ϕ, ψ; α)/∂b 2 . Applying the definition of the first-order total G-differential to Equation (153) yields the following relation: and The direct-effect term defined in Equation (209) (176) and (41) yields the following 2 nd -LVSS:

But Equation (176) involves the function δϕ(z, ω), which is the solution of Equation (41) subject to the boundary condition provided in Equation (43). Concatenating Equations
(211) The 2 nd -LVSS represented by Equation (211), together with the corresponding boundary conditions provided in Equations (178) and (43), depends on the parameter variations and would need to be solved as many times as there are parameter variations. The need for solving the 2 nd -LVSS can be circumvented by applying the principles of the 2 nd -CASAM-L to construct a 2 nd -LASS. In view of the structure of Equation (211), the Hilbert space appropriate for constructing the corresponding adjoint system will be endowed with the inner product defined in Equation (140). The 2 nd -LASS corresponding to Equation (211) is constructed as follows: with Equation (211) to obtain the following relation: (ii) Integrate by parts the left-side of Equation (212) to obtain the following relation: −ω dc −ω dc (iv) Determine the boundary conditions for the functions c (2) 1 (z, ω) and c (2) 2 (z, ω) by eliminating the unknown quantities δϕ(b 2 , ω) and δa (1) 2 (b 2 , ω) from Equation (213). This is accomplished by imposing the following boundary conditions: The relations provided in Equations (214)-(217) constitute the 2 nd -LASS for the 2 ndlevel adjoint function c (2) Notably, both components 1 (z, ω) and c and δa (1) 2 (z, ω) to obtain the following expression for the indirect-effect term: Inserting the expression obtained in Equation (218) together with the expression of the direct-effect term obtained in Equation (209) into Equation (208) and identifying in the resulting expression the quantities that multiply the various parameter variations yields the following expressions for the corresponding 2 nd -order sensitivities: Solving the 2 nd -LASS yields the following expressions for the components of the 2 nd -level adjoint function c (2) Replacing the expressions obtained in Equations (226) and (227) together with the expressions provided for ϕ(z, ω) and a 2 (z, ω) in Equations (8) and (160) respectively, into Equation (219) yields the following expression: The expression obtained in Equation (228) is identical to the expression that would be obtained by taking the derivative ∂/∂b 2 of Equation (161), which provides a mutual verification of the correctness of the derivations (in particular: the derivation/computation of the adjoint functions involved) underlying the determination of the mixed 2 nd -order sensitivity ∂ 2 ρ/∂µ(α)∂b 2 .
Replacing the expression obtained in Equation (227) into Equation (220) and performing the respective integration yields the following expression: The expression obtained in Equation (229) is identical to the expression that would be obtained by taking the derivative ∂/∂b 2 of Equation (162), which provides a mutual verification of the correctness of the derivations (in particular: the derivation/computation of the adjoint functions involved) underlying the determination of the mixed 2 nd -order sensitivity ∂ 2 ρ/∂µ s (α)∂b 2 .
Replacing the expression obtained in Equation (227) into Equation (221) and performing the respective integration yields the following expression: The expression obtained in Equation (230) is identical to the expression that would be obtained by taking the derivative ∂/∂b 2 of Equation (163), which provides a mutual verification of the correctness of the derivations (in particular: the derivation/computation of the adjoint functions involved) underlying the determination of the mixed 2 nd -order sensitivity ∂ 2 ρ/∂q∂b 2 .
Replacing the expression obtained in Equation (160) into Equation (222) and performing the respective integration yields the following expression: The expression obtained in Equation (231) is identical to the expression that would be obtained by taking the derivative ∂/∂b 2 of Equation (164), which provides a mutual verification of the correctness of the derivations (in particular: the derivation/computation of the adjoint functions involved) underlying the determination of the mixed 2 nd -order sensitivity ∂ 2 ρ/∂q * ∂b 2 .
Replacing the expression obtained in Equation (227) into Equation (223) and performing the respective integration yields the following expression: The expression obtained in Equation (232) is identical to the expression that would be obtained by taking the derivative ∂/∂b 2 of Equation (165), which provides a mutual verification of the correctness of the derivations (in particular: the derivation/computation of the adjoint functions involved) underlying the determination of the mixed 2 nd -order sensitivity ∂ 2 ρ/∂b 1 ∂b 2 .
Using the expressions obtained in Equations (8) and (226) in Equation (224) and performing the respective integrations yields the following expression: The expression obtained in Equation (233) is identical to the expression that would be obtained by taking the derivative ∂/∂b 2 of Equation (167), which provides a mutual verification of the correctness of the derivations (in particular: the derivation/computation of the adjoint functions involved) underlying the determination of the mixed 2 nd -order sensitivity ∂ 2 ρ/∂z 0 ∂b 2 .
Using the expression obtained in Equation (160) in Equation (225) and performing the respective integration yields the following expression: Since ∂ 2 ρ/∂b 2 ∂b 2 is an unmixed 2 nd -order sensitivity, its expression cannot be verified independently in terms of an equivalent expression using different adjoint functions, as is the case for the mixed 2 nd -order sensitivities, the expressions of which can always be independently verified in terms of alternative expressions involving different adjoint and/or forward functions.
It has been shown in this section that the exact computation of all of the partial secondorder sensitivities requires 7 large-scale (adjoint) computations to determine the 7 s-level adjoint functions a (2) 2 (z, ω) † , g (2) (z, ω) g 1 (z, ω), g 2 (z, ω) † and h (2) 2 (z, ω) † , by solving the corresponding 2 nd -Level Adjoint Sensitivity Systems. As has been mentioned under item (iv), above, the 2 nd -order mixed sensitivities 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 nd -order 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 1 st -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.

Application of the 3 rd -CASAM-L to Compute Third-Order Response Sensitivities to Imprecisely Known Parameters
The 3 rd -order response sensitivities to the model parameters will be obtained by applying the general principles presented in [1]. In essence, each of the second-order response sensitivities will be considered a new "model response" and the principles underlying the 1 st -CASAM-L will be applied to determine the "1 st -order sensitivities of the 2 nd -order sensitivities", which by definition are the sought-after 3 rd -order sensitivities. The application of these principles will be illustrated in this section by considering as "model response" the 2 nd -order unmixed sensitivity ∂ 2 ρ(ϕ, ψ; α)/∂µ(α)∂µ(α) of the "contributon-response flux" with respect to the interaction coefficient µ(α). This sensitivity is chosen to demonstrate the principles underlying the n th -CASAM-L for computing high-order sensitivities because its expression, as provided in Equation (197), involves all of the original (forward and adjoint) and 1 st -level, and 2 nd -level adjoint state functions, and consequently embodies all of the conceptual aspects underlying the computation of 3 rd -order sensitivities. The determination of the other 3 rd -order sensitivities involves only partially the mathematical complexities that are needed to determine the 3 rd -order sensitivities which will arise from ∂ 2 ρ(ϕ, ψ; α)/∂µ(α)∂µ(α).
Applying the definition of the first-order total G-differential to Equation (235) yields the following relation: where δρ (2) and δρ (2) The direct-effect term δρ (2) u (3) (z, ω); α; δα α 0 dir can be computed immediately, since all of the quantities appearing in its definition on the right-side of Equation (238) are known.

4.
Use the function a 5 (z, ω) and determine the 3 rd -level adjoint function a 8 (z, ω) by solving the following equation:
The various 3 rd -order sensitivities with respect to the nuclide number densities and microscopic cross sections, respectively, can be obtained by using in Equations (280)-(286) the relations provided in Equations (155)-(158) together with the relations provided in Equations (204)-(206). The relative 3 rd -order sensitivities corresponding to these parameters have equal values.

Determination of the Third-Order
Sensitivities of the Form ∂ 3 ρ(ϕ, ψ;α)/∂α i ∂q∂b 2 , i = 1, . . . , TP This section will present the explicit determination of 3 rd -order sensitivities which involve both the inner (interface) and the outer boundary parameters, b 1 and b 2 , in order to illustrate the general applicability of the n th -CASAM-L (for the particular case n = 3) to determine high-order (in this case: 3 rd -order) sensitivities that involve boundary and interface parameters. These 3 rd -order sensitivities are obtained by considering the sensitivity ∂ 2 ρ/∂b 1 ∂b 2 as the "model response" and by taking its first G-differential. The expression of ∂ 2 ρ/∂b 1 ∂b 2 is given in Equation (223). By definition, the G-differential of ∂ 2 ρ/∂b 1 ∂b 2 is obtained as follows: where and The direct-effect term defined in Equation (288) can be computed immediately, but the computation of the indirect-effect term defined in Equation (289) requires the prior determination of the function δc (2) 1 (z, ω). The function δc (2) 1 (z, ω) is the solution of the system obtained by taking the first-order G-differential of Equations (214)-(217). Applying the definition of the G-differential to Equations (214)-(217) yields the following system of equations: −ω d dz δc −ω d dz δc The 3 rd -LVSS represented by Equations (290)-(293) depends on the parameter variations and would need to be solved as many times as there are parameter variations. The need for solving this 3 rd -LVSS can be circumvented by applying the principles of the 3 rd -CASAM-L to construct a corresponding 3 rd -LASS. In view of the structure of Equations (290)-(293), the Hilbert space appropriate for constructing the corresponding adjoint system is the Hilbert space H 1 , which comprises two-component vector-valued elements and which is endowed with the inner product defined in Equation (140). The 3 rd -LASS corresponding to Equations (290)-(293) is constructed as follows: (i) Form the inner product of a 2-component vector function 2 (z, ω) † with Equations (290) and (291) to obtain the following relation: 1 (z, ω) + µδc (2) 1 (z, ω) − δ(z − z 0 )δc (2) 2 (z, ω) dz α 0 2 (z, ω) −ω d dz δc (2) 2 (z, ω) + µδc 1 (z, ω) + µδc 2 (z, ω) + µδc

Discussion and Conclusions
The n th -Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems (abbreviated as "n th -CASAM-L"), which is presented in the accompanying work [1], enables the most efficient computation of exactly obtained 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 −hence the qualifier "comprehensive") underlying the respective forward/adjoint systems. The application of the n th -CASAM-L has been illustrated in this work by considering a paradigm model which describes the transmission of particles produced by a distributed source through a shield which surrounds the source. The sensitivities of the most important types of responses for particle transport (namely: point-detector responses, particle leakage responses, reaction rate responses and "contributon-response fluxes") have been analyzed, highlighting differences and similarities among them. This paradigm model can be used as a benchmark for comparing the n th -CASAM-L to conventional statistical and/or finite-difference methods, since the model is sufficiently simple to admit closed-form exact expressions for the response sensitivities of any order yet contains many parameters (7 to 10 primary parameters and tens of additional fundamental parameters describing imprecisely known cross sections and number densities of the nuclides corresponding to the interaction coefficients of the particles within the materials in which the particles propagate). It has been also shown that the relative first-order sensitivities of the responses to the nuclide number densities have the same values as the corresponding first-order sensitivities of the responses to the microscopic cross sections. The precise and efficient computation of many relative sensitivities that have identical values does not pose any difficulties to the n th -CASAM-L methodology but is known to pose extreme (even insurmountable) challenges to conventional statistical methods which rely on sensitivities being readily "sortable" by having relative values that clearly differ from each other.
Application of the 1 st -CASAM-L yields exact expressions of the 1 st -order sensitivities in terms of a two-component first-level adjoint function of the form a (1) 2 (z, ω) † , which is the solution of a First-Level Adjoint Sensitivity System (1 st -LASS) that is independent of parameter variations. Hence, the 1 st -LASS needs to be solved only once, regardless of the number of model parameters. Subsequently, the 1 st -order response sensitivities are computed efficiently and exactly in terms of the first-level adjoint function by simply performing inexpensive quadratures rather than solving operator (differential) equations, which represent "large-scale computations".
The 2 nd -order sensitivities are expressed in terms of a 2 nd -level adjoint function, which is the solution of a 2 nd -Level Adjoint Sensitivity System (2 nd -LASS) and which is also independent of parameter variations. The 2 nd -LASS comprises at most 4 differential equations, having as solution a 2 nd -level adjoint function of the form a (2) 4 (z, ω) † . For the computation of most 2 nd -order sensitivities, however, the 2 nd -LASS comprises only 2 differential equations. In all cases, though, the equations underlying this 2 nd -LASS need not be solved simultaneously. Instead, they are solved successively, decoupled from one another. The 2 nd -order mixed sensitivities can 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 1 st -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.
The 3 rd -order sensitivities are determined in terms of 3 rd -level adjoint functions, which are solutions of 3 rd -Level Adjoint Sensitivity Systems (3 rd -LASS), which can comprise at most 8 operator-equations, in which case the corresponding 3 rd -level adjoint function has the form a (3) 1 (z, ω), . . . , a 8 (z, ω) † . As in the case of the 2 nd -LASS, the equations underlying the 3 rd -LASS need not be solved simultaneously since they can be readily decoupled when written in component form and can therefore be solved successively. The computation of the 3 rd -order sensitivities can be prioritized by using the magnitudes/importance of the 1 st -order and 2 nd -order sensitivities as guiding indicators. The sensitivities of order higher than third were not explicitly determined for the illustrative paradigm model presented in work, but the principles for obtaining them have been amply illustrated. The illustrative model analyzed in this work has confirmed the fundamental advantage of applying the n th -CASAM-L (as opposed to statistical and/or finite-difference methods) for sensitivity analysis of large-scale models with many parameters, which can be summarized as follows: for a model having TP-parameters (where "TP" denotes "total number of parameters") for every response and/or sensitivity of interest, the n th -CASAM-L computes the "TP next-higher-order" sensitivities needing just 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 various-order sensitivities.
The n th -CASAM-L 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 various 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 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-L-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-L-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.