On the use of the adjoint operator for source reconstruction in particle transport problems

ABSTRACT

In this work, the adjoint to the transport operator is used to estimate the spatial distribution of an isotropic neutral particles source in a homogeneous one-dimensional medium, from readings of internal detectors. An analytical discrete ordinates formulation, the ADO method, is applied to derive a spatially explicit solution for the adjoint flux. Simulations are performed for reconstructing Gaussian and piecewise localized sources. Numerical results indicate that the source estimates may be considered satisfactory although scalar fluxes are reconstructed more precisely.

KEYWORDS:

• Inverse problems
• source reconstruction
• analytical methods
• discrete ordinates
• adjoint transport equation

AMS SUBJECT CLASSIFICATIONS:

• 35Q20
• 65L09
• 65M32
• 65N35

1. Introduction

Several problems of practical interest in science and engineering require the solution of inverse particle transport problems. Among them, we may cite non-destructive material identification [1–4] and nuclear safety [5–7]. In this sense, a series of different approaches can be found in the literature on this subject [8]. For instance, in source reconstruction problems, typical in safety analysis, either readings on the boundaries [1,3,6,9,10] or in the interior of the domain [5,11] is considered in order to estimate the source distribution, usually approximated in given basis functions [1,3,6,9,10]. The adjoint to the transport operator is an additional tool for the solution of such problems, either in iterative schemes, when computing gradients [1,6,12], or to estimate detector readings in a source–detector configuration [5,13].

Our goal, in this paper, is the source term reconstruction in the particle transport equation. In this way, we estimate the projection coefficients of the source function expansion in a given basis, using measurements of particle detectors placed within the domain. The detectors may be distributed over the entire domain as well as concentrated in specific regions. Effects of noise on the measurements are also taken into account in our analysis. We consider that the reconstructed source may be approximated by sectionally constant functions defined over the physical domain. Such a choice leads to low-order systems for the minimization procedure and to simpler implementation schemes. Using the same basis functions, we have successfully reconstructed polynomials and piecewise sources, assuming the basis functions support was defined over intervals where the original source had no discontinuities [14]. Here, when defining the basis function support, no conditions are imposed with respect to possible discontinuities of the source. In addition, thinking of some of the previous applications mentioned, we extend our analysis to a wide class of source functions, trying to represent the occurrence of peaks of localized intensity.

We assume all physical properties of the medium and incoming fluxes at the boundaries are given, then we follow Hykes and Azmy [5] to derive, using the adjoint of the transport equation, a linear model that relates the detector readings with the coefficients of our expansion. The coefficients are estimated by the minimization of the 2-norm of the discrepancy between the numerical readings and the observed ones. As an additional contribution of this paper, we derive analytical expressions with regard to the spatial variable for the matrix entries of the minimization problem, a consequence of the use of the Analytical Discrete Ordinates (ADO) method for solving the adjoint transport equation [15]. Such feature is expected to contribute with computational time reduction and solution accuracy. The ADO method is a fast, accurate and concise variation of the classical Chandrasekhar’s discrete ordinates method ($S_n$ method) [16], which provides spatially closed-form solutions to the discrete ordinates approximation of the transport equation [17,18] and its adjoint [15]. More recently the ADO method has been successfully extended, with aid of nodal techniques, to solve bi-dimensional transport problems [19,20], to which the present derivations may be associated later.

In this work, we consider monoenergetic, slab-geometry transport problems with arbitrarily high anisotropic scattering. As boundary conditions, we assume either known incoming fluxes or reflective boundary conditions. The medium is composed of materials with homogeneous physical properties, although simple modifications would allow us to treat heterogeneous media [15] with sectionally constant physical properties.

The remainder of this paper is organized such that: in Section 2, we briefly discuss the transport operator and relations with its adjoint operator; in Section 3, we present the source reconstruction strategy; in Section 4, we summarize the ADO methodology for obtaining a solution for the one-dimensional discrete ordinates adjoint transport equation and deduce analytical expressions to the matrix related to the inversion problem derived in Section 3 and in Section 5, we apply the source reconstruction strategy in two test problems aiming the reconstruction of a Gaussian and a localized piecewise constant source. To conclude, a few remarks and suggestions for future work are discussed in Section 6.

2. Mathematical formulation

Let us suppose that an isotropic source of neutral particles S lies within a plane-parallel medium, a slab of optical thickness $z_0$. Thus, the monoenergetic distribution of particles within the slab may be modelled by the transport equation(1) $$\begin{array}{c}\hfill \mathcal{L}\psi =S(z)\end{array}$$(1)

with $\mathcal{L}$ being the transport operator, here written as [21,22](2) $$\begin{array}{c}\hfill \mathcal{L}\psi (z,\mu )=\mu \frac{\partial }{\partial z}\psi (z,\mu )+\psi (z,\mu )-\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l(\mu ){\int }_{-1}^1{P}_l({\mu }^{\prime })\psi (z,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime },\end{array}$$(2)

where $\psi$ is the angular flux of neutral particles, which in this work, we assume to be symmetrical with respect to the azimuthal variable and $\mu \in [-1,1]$ is the cosine of the polar angle measured from the positive z-axis, with $z\in (0,z_0)$. Besides that, we represent the mean number of neutral particles emerging from collisions by c and we set ${\beta }_l$ as the coefficient of the lth term in the Lth order expansion of the scattering law in terms of Legendre polynomials $P_l$. Later on, in this text, we mention specific forms of defining those coefficients, depending on the applications. In fact, thinking of general applications of the particle transport equation, even the unknown $\psi$ and the parameter c may named as alternative quantities (angular intensity, albedo, respectively, for instance). Additionally, we prescribe boundary conditions related to more general problems in transport theory [21,23], written as (3a) $$\begin{array}{cc}\hfill \psi (0,\mu )& =g_1(\mu )+{\alpha }_1\psi (0,-\mu ),\hfill \end{array}$$(3a) (3b) $$\begin{array}{cc}\hfill \psi (z_0,-\mu )& =g_2(\mu )+{\alpha }_2\psi (z_0,\mu ),\hfill \end{array}$$(3b) for $\mu \in (0,1]$, with known incoming fluxes $g_1$ and $g_2$ at the boundaries, and reflection coefficients ${\alpha }_1$, ${\alpha }_2\in [0,1]$.

Figure 1. Neutral particle detector placed within the slab.

Before we proceed, it is useful to define an inner product on the phase space $[-1,1]\times [0,z_0]$ as(4) $$\begin{array}{c}\hfill ⟨f,g⟩={\int }_{-1}^1{\int }_0^{z_0}f(z,\mu )g(z,\mu )\mathrm{d}z\mathrm{d}\mu ,\end{array}$$(4)

for $f,g:[-1,1]\times [0,z_0]\to {\mathbb{R}}_{+}$ integrable functions, since, according to the literature [21,24], if we consider a detector of neutral particles, placed within a slab (see Figure 1), then(5) $$\begin{array}{c}\hfill r=⟨{\sigma }_d,\psi ⟩\end{array}$$(5)

is the detector response due to particles migrating in all directions within segment $[z_a,z_b]$. In our (dimensionless) formulation, ${\sigma }_d:[0,z_0]\to {\mathbb{R}}_{+}$ is the ratio between the macroscopic absorption cross-section of the detector and the total cross-section of the medium. For a fixed contiguous region $[z_a,z_b]\subseteq (0,z_0)$, we define(6) $$\begin{array}{c}\hfill {\sigma }_d(z)=\left\{\begin{array}{cc}{\sigma }_d,\hfill & z\in [z_a,z_b],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(6)

For every $\mu \in [-1,1]$, with ${\sigma }_d>0$ a constant value, we refer here as absorption parameter.

Alternatively, an efficient approach to determine the detector response r is given by the solution of the adjoint transport problem, which may be derived by assuming the existence of a smooth function ${\psi }^{†}:[-1,1]\times [0,z_0]\to {\mathbb{R}}_{+}$, which we refer to as the adjoint angular flux. We multiply $\mathcal{L}\psi$ by ${\psi }^{†}$, then we use the inner-product defined in Equation (4(4) $$\begin{array}{c}\hfill ⟨f,g⟩={\int }_{-1}^1{\int }_0^{z_0}f(z,\mu )g(z,\mu )\mathrm{d}z\mathrm{d}\mu ,\end{array}$$(4) ), which yields, after some mathematical manipulation, the result [21,24](7) $$\begin{array}{c}\hfill ⟨{\psi }^{†},\mathcal{L}\psi ⟩=P(\psi ,{\psi }^{†})+⟨{\mathcal{L}}^{†}{\psi }^{†},\psi ⟩\end{array}$$(7)

with $P(\psi ,{\psi }^{†})$ defined as(8) $$\begin{array}{c}\hfill P(\psi ,{\psi }^{†})={\int }_{-1}^1\mu [\psi (z_0,\mu ){\psi }^{†}(z_0,\mu )-\psi (0,\mu ){\psi }^{†}(0,\mu )]\mathrm{d}\mu \end{array}$$(8)

and ${\mathcal{L}}^{†}$, the adjoint transport operator, as [21,24](9) $$\begin{array}{c}\hfill {\mathcal{L}}^{†}{\psi }^{†}(z,\mu )=-\mu \frac{\partial }{\partial z}{\psi }^{†}(z,\mu )+{\psi }^{†}(z,\mu )-\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l(\mu ){\int }_{-1}^1{P}_l({\mu }^{\prime }){\psi }^{†}(z,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime }.\end{array}$$(9)

Here ${\psi }^{†}$ is the adjoint angular flux of particles, and the remaining quantities are defined as in Equation (2(2) $$\begin{array}{c}\hfill \mathcal{L}\psi (z,\mu )=\mu \frac{\partial }{\partial z}\psi (z,\mu )+\psi (z,\mu )-\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l(\mu ){\int }_{-1}^1{P}_l({\mu }^{\prime })\psi (z,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime },\end{array}$$(2) ).

At this point, we write the adjoint transport equation as(10) $$\begin{array}{c}\hfill {\mathcal{L}}^{†}{\psi }^{†}=S^{†}\end{array}$$(10)

for a certain adjoint source $S^{†}$ yet to be chosen. If the adjoint source is chosen to be equal to ${\sigma }_d(z)$, then we may rewrite Equation (5(5) $$\begin{array}{c}\hfill r=⟨{\sigma }_d,\psi ⟩\end{array}$$(5) ) as(11) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},S⟩-P(\psi ,{\psi }^{†})\end{array}$$(11)

using Equation (7(7) $$\begin{array}{c}\hfill ⟨{\psi }^{†},\mathcal{L}\psi ⟩=P(\psi ,{\psi }^{†})+⟨{\mathcal{L}}^{†}{\psi }^{†},\psi ⟩\end{array}$$(7) ), thus, allowing us to compute the detector response using the adjoint angular flux. Furthermore, since Equation (8(8) $$\begin{array}{c}\hfill P(\psi ,{\psi }^{†})={\int }_{-1}^1\mu [\psi (z_0,\mu ){\psi }^{†}(z_0,\mu )-\psi (0,\mu ){\psi }^{†}(0,\mu )]\mathrm{d}\mu \end{array}$$(8) ) depends on $\psi$ and ${\psi }^{†}$ only at the boundaries of the slab, the practical choice of boundary conditions to the adjoint angular flux problem (12a) $$\begin{array}{cc}\hfill {\psi }^{†}(0,-\mu )& ={\alpha }_1{\psi }^{†}(0,\mu ),\hfill \end{array}$$(12a) (12b) $$\begin{array}{cc}\hfill {\psi }^{†}(z_0,\mu )& ={\alpha }_2{\psi }^{†}(z_0,-\mu ),\hfill \end{array}$$(12b) results in the simplification of Equation (8(8) $$\begin{array}{c}\hfill P(\psi ,{\psi }^{†})={\int }_{-1}^1\mu [\psi (z_0,\mu ){\psi }^{†}(z_0,\mu )-\psi (0,\mu ){\psi }^{†}(0,\mu )]\mathrm{d}\mu \end{array}$$(8) ) to [15](13) $$\begin{array}{c}\hfill P(\psi ,{\psi }^{†})=-{\int }_0^1\mu [g_1(\mu ){\psi }^{†}(0,\mu )+g_2(\mu ){\psi }^{†}(z_0,-\mu )]\mathrm{d}\mu .\end{array}$$(13)

We remark that, since the adjoint angular flux ${\psi }^{†}$ does not depend on S, r can be computed by Equation (11(11) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},S⟩-P(\psi ,{\psi }^{†})\end{array}$$(11) ) for any S, given $g_1$ and $g_2$. Thus, if one changes source distribution S and the prescribed boundary conditions $g_1$ and $g_2$, Equation (11(11) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},S⟩-P(\psi ,{\psi }^{†})\end{array}$$(11) ) can still be used to evaluate r, as long as the detector is not changed or displaced. On the other hand, the use of Equation (5(5) $$\begin{array}{c}\hfill r=⟨{\sigma }_d,\psi ⟩\end{array}$$(5) ) requires the transport equation to be solved every time the source distribution and/or the incoming fluxes at the boundaries are changed. This feature along with the bilinearity of Equation (4(4) $$\begin{array}{c}\hfill ⟨f,g⟩={\int }_{-1}^1{\int }_0^{z_0}f(z,\mu )g(z,\mu )\mathrm{d}z\mathrm{d}\mu ,\end{array}$$(4) ) allow us to deal with the linear inverse problem for the source reconstruction, as we will discuss in the next section.

3. Source reconstruction strategy

First of all, we point out some aspects about the configuration of our reconstruction problem. Following Hykes and Azmy [5],

(i)	a set of D particle detectors is placed within the slab $(0,z_0)$;
(ii)	${\sigma }_{d,i}$ is the absorption parameter of the ith detector;
(iii)	for each detector, the adjoint angular flux that solves ${\mathcal{L}}^{†}{\psi }^{†}={\sigma }_{d,i}$ is available;
(iv)	the original source of neutral particles S may be approximated by the projection of S onto a linear space with known basis function $f_j$, $j=1,\cdots ,B$.

According to (iv), S is approximated by(14) $$\begin{array}{c}\hfill \widehat{S}(z)=\sum _{j=1}^B{s}_j{f}_j(z),\end{array}$$(14)

with constants $s_j$, targets of our source reconstruction process, yet to be found. Under these assumptions, we can express the detector response due to neutral particles migrating in all directions within the ith detector region as(15) $$\begin{array}{c}\hfill r_i=⟨{\psi }_i^{†},\widehat{S}⟩-P(\psi ,{\psi }_i^{†})=\sum _{j=1}^B{s}_j⟨{\psi }_i^{†},f_j⟩-P(\psi ,{\psi }_i^{†})\end{array}$$(15)

for $i=1,\cdots ,D$. We consider the vectors $\mathbf{r}=[r_i]\in {\mathbb{R}}^D$, $\mathbf{p}=[P(\psi ,{\psi }_i^{†})]\in {\mathbb{R}}^D$, and the matrix $\mathbf{A}=[A_{i,j}]\in {\mathbb{R}}^{D\times B}$, such that(16) $$\begin{array}{c}\hfill A_{i,j}=⟨{\psi }_i^{†},f_j⟩.\end{array}$$(16)

Then, we are able to rewrite Equation (15(15) $$\begin{array}{c}\hfill r_i=⟨{\psi }_i^{†},\widehat{S}⟩-P(\psi ,{\psi }_i^{†})=\sum _{j=1}^B{s}_j⟨{\psi }_i^{†},f_j⟩-P(\psi ,{\psi }_i^{†})\end{array}$$(15) ) in a matrix form as(17) $$\begin{array}{c}\hfill \mathbf{r}=\mathbf{A}s-\mathbf{p}\end{array}$$(17)

with $\mathbf{s}=[s_j]\in {\mathbb{R}}^B$. At this point, we use Equation (17(17) $$\begin{array}{c}\hfill \mathbf{r}=\mathbf{A}s-\mathbf{p}\end{array}$$(17) ) to write our inverse problem as a minimization problem of the discrepancy between the numerical readings $\mathit{r}$ and observed data. For this, we first suppose the existence of noisy measurements $r_{m,i}$, computed by numerical simulation or through physical experimentation. In this way, the coefficients, components of the vector $\mathbf{s}$ in Equation (17(17) $$\begin{array}{c}\hfill \mathbf{r}=\mathbf{A}s-\mathbf{p}\end{array}$$(17) ), may be estimated by the minimization of the objective function [5](18) $$\begin{array}{c}\hfill h(\mathbf{s})=||{\mathbf{r}}^{\prime }-\mathbf{A}s{||}_2^2\end{array}$$(18)

with ${\mathbf{r}}^{\prime }={\mathbf{r}}_m-\mathbf{p}$, ${\mathbf{r}}_m=[r_{m,i}]\in {\mathbb{R}}^D$, where $||·{||}_2$ is the classical Euclidean in ${\mathbb{R}}^B$.

At this point, basis functions for the reconstruction space must be fixed. In this work, we only consider sectionally constant approximations to the neutral particles source S. Thus, given(19) $$\begin{array}{c}\hfill [0,z_0]=\bigcup _{j=1}^B[z_{j-1},z_j],\end{array}$$(19)

a partition of the domain, basis functions are defined as [5](20) $$\begin{array}{c}\hfill f_j(z)=\left\{\begin{array}{cc}1,\hfill & z\in [z_{j-1},z_j],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(20)

for $j=1,\cdots ,B$.

We remark that the function space generated by the basis functions in Equation (20(20) $$\begin{array}{c}\hfill f_j(z)=\left\{\begin{array}{cc}1,\hfill & z\in [z_{j-1},z_j],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(20) ) may be too poor for representing the true space of the particle source S, in such way that the well-known ill-posedness of inverse problems might negatively affect the quality of the reconstruction. Moreover, we found in our numerical tests that the identifiability [25] of our inverse model is $|\mathit{A}{\mathit{A}}^T|\approx 0$, leading to difficulties to uniquely determine the coefficients. We treat these issues here by searching for Tikhonov regularized solutions for the minimization problem defined by the objective function in Equation (18(18) $$\begin{array}{c}\hfill h(\mathbf{s})=||{\mathbf{r}}^{\prime }-\mathbf{A}s{||}_2^2\end{array}$$(18) ). In fact, we search for solutions that minimize the objective function [26](21) $$\begin{array}{c}\hfill h_{\lambda }(\mathbf{s})=||{\mathbf{r}}^{\prime }-{\mathbf{A}s||}_2^2+{\lambda }^2{||\mathbf{T}s||}_2^2,\end{array}$$(21)

where $\lambda$ is the Tikhonov’s regularization parameter, here chosen by the Morozov’s discrepancy principle [26]. We consider, in this work, three choices for the Tikhonov’s regularization matrix $\mathbf{T}$: $\mathbf{T}=T_0=I$, the identity matrix in ${\mathbb{R}}^{B\times B}$, which selects solutions of minimal norm; $\mathbf{T}=T_1\in {\mathbb{R}}^{(B-1)\times B}$, (22a) $$\begin{array}{c}\hfill {\mathbf{T}}_1=\left(\begin{array}{cccc}-1& 1& & \\ & \ddots & \ddots & \\ & & -1& 1\end{array}\right)\end{array}$$(22a)

which selects solutions that minimizes $\sum _{j=1}^{B-1}{(s_j-s_{j+1})}^2$ and $\mathbf{T}=T_2\in {\mathbb{R}}^{(B-2)\times B}$ (22b) $$\begin{array}{c}\hfill {\mathbf{T}}_2=\left(\begin{array}{ccccc}1& -2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& -2& 1\end{array}\right).\end{array}$$(22b)

which minimizes $\sum _{j=1}^{B-2}{(s_j-2s_{j+1}+s_{j+2})}^2$. In fact, ${\mathit{T}}_1$ and ${\mathit{T}}_2$ are associated with finite difference approximations of the first and the second derivatives of the source term, respectively.

In the next section, before proceeding to the numerical simulations, we describe the Analytical Discrete Ordinates method (ADO), and we use to solve the adjoint to the transport problem and to derive explicit expressions to Equation (16(16) $$\begin{array}{c}\hfill A_{i,j}=⟨{\psi }_i^{†},f_j⟩.\end{array}$$(16) ).

4. An analytical solution of the discrete ordinates adjoint transport equation

The inversion strategy we previously described requires the adjoint transport equation ${\mathcal{L}}^{†}{\psi }^{†}={\sigma }_{d,i}$ to be solved a number of times equal to the number of detectors. A fast, concise and accurate approach is the Analytical Discrete Ordinates (ADO) method, originally derived by Barichello and Siewert [17,23] for solving the discrete ordinates approximation of the transport equation. The extension of such approach for solving the adjoint one-dimensional equation is straightforward, and it was done by Pazinatto et al. [15,27].

Since Equation (10(10) $$\begin{array}{c}\hfill {\mathcal{L}}^{†}{\psi }^{†}=S^{†}\end{array}$$(10) ) is a linear equation, we write the adjoint angular flux ${\psi }^{†}$ as the solution ${\psi }_h^{†}$ of the homogeneous equation ${\mathcal{L}}^{†}{\psi }_h^{†}=0$, plus any particular solution ${\psi }_p^{†}$ of the non-homogeneous equation. For the sake of simplicity, we suppose that all physical parameters in the equation are constant. Otherwise, the domain is partitioned in regions of constant parameters, and we impose continuity conditions on the adjoint angular flux at the interfaces of contiguous regions.

4.1. Homogeneous solution

We start by rewriting the homogeneous version of Equation (10(10) $$\begin{array}{c}\hfill {\mathcal{L}}^{†}{\psi }^{†}=S^{†}\end{array}$$(10) ) in the form [15](23) $$\begin{array}{cc}& -\mu \frac{\partial }{\partial z}{\psi }_h^{†}(z,\mu )+{\psi }_h^{†}(z,\mu )\hfill \\ \hfill & =\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l(\mu ){\int }_0^1{P}_l({\mu }^{\prime })[{\psi }_h^{†}(z,{\mu }^{\prime })+{(-1)}^l{\psi }_h^{†}(z,-{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime },\hfill \end{array}$$(23)

by changing the integration variable from $[-1,1]$ to [0, 1], then applying the Legendre polynomials rule $P_l(-\mu )={(-1)}^l{P}_l(\mu )$ [28]. Furthermore, we consider a set of N nodes and weights $\{{\mu }_i,w_i\}$ of a quadrature scheme for the interval [0, 1], then the angular dependence on Equation (23(23) $$\begin{array}{cc}& -\mu \frac{\partial }{\partial z}{\psi }_h^{†}(z,\mu )+{\psi }_h^{†}(z,\mu )\hfill \\ \hfill & =\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l(\mu ){\int }_0^1{P}_l({\mu }^{\prime })[{\psi }_h^{†}(z,{\mu }^{\prime })+{(-1)}^l{\psi }_h^{†}(z,-{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime },\hfill \end{array}$$(23) ) is discretized in the directions $\mu =\pm {\mu }_i$ as(24) $$\begin{array}{cc}& \mp {\mu }_i\frac{\mathrm{d}}{\mathrm{d}z}{\psi }_h^{†}(z,\pm {\mu }_i)+{\psi }_h^{†}(z,\pm {\mu }_i)\hfill \\ \hfill & =\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l({\mu }_i)\sum _{k=1}^N{w}_k{P}_l({\mu }_k)\left[{\psi }_h^{†}(z,\pm {\mu }_k)+{(-1)}^l{\psi }_h^{†}(z,\mp {\mu }_k)\right]\hfill \end{array}$$(24)

for $i=1,\cdots ,N$. Moving on, we seek for spectral solutions for Equation (24(24) $$\begin{array}{cc}& \mp {\mu }_i\frac{\mathrm{d}}{\mathrm{d}z}{\psi }_h^{†}(z,\pm {\mu }_i)+{\psi }_h^{†}(z,\pm {\mu }_i)\hfill \\ \hfill & =\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l({\mu }_i)\sum _{k=1}^N{w}_k{P}_l({\mu }_k)\left[{\psi }_h^{†}(z,\pm {\mu }_k)+{(-1)}^l{\psi }_h^{†}(z,\mp {\mu }_k)\right]\hfill \end{array}$$(24) ) in the form(25) $$\begin{array}{c}\hfill {\psi }_h^{†}(z,\mu )=\varphi (\nu ,\mu )e^{-z/\nu }\end{array}$$(25)

where $\nu$ is an arbitrary constant. Then, we substitute Equation (25(25) $$\begin{array}{c}\hfill {\psi }_h^{†}(z,\mu )=\varphi (\nu ,\mu )e^{-z/\nu }\end{array}$$(25) ) into Equation (24(24) $$\begin{array}{cc}& \mp {\mu }_i\frac{\mathrm{d}}{\mathrm{d}z}{\psi }_h^{†}(z,\pm {\mu }_i)+{\psi }_h^{†}(z,\pm {\mu }_i)\hfill \\ \hfill & =\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l({\mu }_i)\sum _{k=1}^N{w}_k{P}_l({\mu }_k)\left[{\psi }_h^{†}(z,\pm {\mu }_k)+{(-1)}^l{\psi }_h^{†}(z,\mp {\mu }_k)\right]\hfill \end{array}$$(24) ), to obtain(26) $$\begin{array}{c}\hfill \left(1\pm \frac{{\mu }_i}{\nu }\right)\varphi (\nu ,\pm {\mu }_i)=\frac{c}{2}\sum _{l=0}^L{\beta }_l{P}_l({\mu }_i)\sum _{k=1}^N{w}_k{P}_l({\mu }_k)\left[\varphi (\nu ,\pm {\mu }_k)+{(-1)}^l\varphi (\nu ,\mp {\mu }_k)\right]\end{array}$$(26)

for $i=1,\cdots ,N$. After a series of algebraic operations carefully described in Pazinatto et al. [15], we can write the homogeneous solution to the adjoint transport equation as(27) $$\begin{array}{c}\hfill {\mathit{\psi }}_{\pm ,h}^{†}(z)=\sum _{j=1}^N[a_j{\mathit{\varphi }}_{\pm }({\nu }_j)e^{-z/{\nu }_j}+b_j{\mathit{\varphi }}_{\mp }({\nu }_j)e^{-(z_0-z)/{\nu }_j}],\end{array}$$(27)

where ${\mathit{\psi }}_{\pm ,h}^{†}(z)=\left[{\psi }_h^{†}(z,\pm {\mu }_i)\right]\in {\mathbb{R}}^N$ and ${\mathit{\varphi }}_{\pm }(\nu )=\left[\varphi (\nu ,\pm {\mu }_i)\right]\in {\mathbb{R}}^N$ with eigenfunctions defined by (28a) $$\begin{array}{c}\hfill {\mathit{\varphi }}_{\pm }(\nu )=\frac{1}{2}{\mathit{M}}^{-1}\left(\mathit{I}\mp \nu {\mathit{B}}_{\mathbf{+}}\right)\mathit{x},\end{array}$$(28a)

with $\mathit{M}=\mathrm{d}iag\left({\mu }_i\right)\in {\mathbb{R}}^{N\times N}$, where $\mathit{x}\in {\mathbb{R}}^N$ and $\nu >0$ are such that(28b) $$\begin{array}{c}\hfill {\mathit{B}}_{-}{\mathit{B}}_{+}\mathit{x}=\frac{1}{{\nu }^2}\mathit{x},\end{array}$$(28b)

with ${\mathit{B}}_{\pm }\in {\mathbb{R}}^{N\times N}$ written as(28c) $$\begin{array}{c}\hfill {\mathit{B}}_{+}=\left({\mathit{I}}_N-\frac{c}{2}\sum _{l=0}^L{\beta }_l{\mathbf{\Pi }}_l{\mathbf{\Pi }}_l^T\mathit{W}[1+{(-1)}^l]\right){\mathit{M}}^{-1}\end{array}$$(28c)

and(28d) $$\begin{array}{c}\hfill {\mathit{B}}_{-}=\left({\mathit{I}}_N-\frac{c}{2}\sum _{l=0}^L{\beta }_l{\mathbf{\Pi }}_l{\mathbf{\Pi }}_l^T\mathit{W}[1-{(-1)}^l]\right){\mathit{M}}^{-1}.\end{array}$$(28d) Here ${\mathbf{\Pi }}_l=\left[P_l({\mu }_i)\right]\in {\mathbb{R}}^N$, $\mathit{W}=\mathrm{d}iag(w_i)\in {\mathbb{R}}^{N\times N}$ and ${\mathit{I}}_N$ is the identity matrix in ${\mathbb{R}}^{N\times N}$.

4.2. Particular solution

For the special case where $S^{†}$ is a constant source, we may choose a constant function as particular solution to Equation (10(10) $$\begin{array}{c}\hfill {\mathcal{L}}^{†}{\psi }^{†}=S^{†}\end{array}$$(10) )(29) $$\begin{array}{c}\hfill {\psi }_p^{†}(z,\mu )=\frac{S^{†}}{1-c{\beta }_0},\end{array}$$(29)

which can easily be verified by direct substitution into the equation, as long as the denominator is different from zero. To deal with more general cases, when $S^{†}$ is not a constant, we followed the approach developed in Ref. [29] based on the use of Green’s functions and extended that to the adjoint transport equation [14]. In this case, we express the particular solution in terms of the eigenfunctions ${\mathit{\varphi }}_{\pm }$ defined in Equation (28a(28a) $$\begin{array}{c}\hfill {\mathit{\varphi }}_{\pm }(\nu )=\frac{1}{2}{\mathit{M}}^{-1}\left(\mathit{I}\mp \nu {\mathit{B}}_{\mathbf{+}}\right)\mathit{x},\end{array}$$(28a) )(30) $$\begin{array}{c}\hfill {\mathit{\psi }}_{\pm ,p}^{†}(z)=\sum _{j=1}^N[a_j(z){\mathit{\varphi }}_{\pm }({\nu }_j)+b_j(z){\mathit{\varphi }}_{\mp }({\nu }_j)]\end{array}$$(30)

where, for an isotropically defined source $S^{†}$, (31a) $$\begin{array}{c}\hfill a_j(z)=c_j{\int }_0^z{S}^{†}(z^{\prime })e^{-(z-z^{\prime })/{\nu }_j}\mathrm{d}{z}^{\prime }\end{array}$$(31a)

and(31b) $$\begin{array}{c}\hfill b_j(z)=c_j{\int }_z^{z_0}{S}^{†}(z^{\prime })e^{-(z^{\prime }-z)/{\nu }_j}\mathrm{d}{z}^{\prime }\end{array}$$(31b)

with(31c) $$\begin{array}{c}\hfill c_j=-\frac{\sum _{i=1}^N{w}_i\left[\varphi ({\nu }_j,{\mu }_i)+\varphi ({\nu }_j,-{\mu }_i)\right]}{\sum _{i=1}^N{w}_i{\mu }_i\left[\varphi {({\nu }_j,{\mu }_i)}^2-\varphi {({\nu }_j,-{\mu }_i)}^2\right]}.\end{array}$$(31c)

We remark that this particular solution formulation allows us to deal more easily with source terms defined as piecewise functions, without imposing any restrictions on the domain partitions. Besides, it reduces the order of the linear system to be solved, in comparison with the use of the constant particular solution, previously mentioned. However, the evaluation of the spatial integrals in Equations () may be costly depending on the quadrature scheme. Such problem may be suppressed when analytical expressions for the integrals are available.

4.3. The detector response: an analytical approach

Once the solution of the discrete ordinates approximation of the adjoint transport equation is established, we are able to use Equations (13(13) $$\begin{array}{c}\hfill P(\psi ,{\psi }^{†})=-{\int }_0^1\mu [g_1(\mu ){\psi }^{†}(0,\mu )+g_2(\mu ){\psi }^{†}(z_0,-\mu )]\mathrm{d}\mu .\end{array}$$(13) ) and (16(16) $$\begin{array}{c}\hfill A_{i,j}=⟨{\psi }_i^{†},f_j⟩.\end{array}$$(16) ) to compute the vector ${\mathit{r}}^{\prime }$ and the matrix $\mathit{A}$ in Equation (21(21) $$\begin{array}{c}\hfill h_{\lambda }(\mathbf{s})=||{\mathbf{r}}^{\prime }-{\mathbf{A}s||}_2^2+{\lambda }^2{||\mathbf{T}s||}_2^2,\end{array}$$(21) ). The evaluation of the relevant integrals may be performed using our angular quadrature scheme and defining additional spatial quadrature scheme. However, with aid of the ADO framework, we are able to derive analytical expressions for those equations.

In order to simplify the notation, we consider a single detector where ${\sigma }_d(z)={\sigma }_d>0$ for $z\in [z_a,z_b]$ and ${\sigma }_d(z)=0$ outside the interval, as in Equation (6(6) $$\begin{array}{c}\hfill {\sigma }_d(z)=\left\{\begin{array}{cc}{\sigma }_d,\hfill & z\in [z_a,z_b],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(6) ), and in Equation (20(20) $$\begin{array}{c}\hfill f_j(z)=\left\{\begin{array}{cc}1,\hfill & z\in [z_{j-1},z_j],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(20) ), without loss of generality, a single basis function, $f(z)=1$, for $z\in [z_c,z_d]$, and $f(z)=0$, otherwise. We note that, differently of previous work [14], no conditions are imposed to relate the intervals where the detectors are placed with the ones where the basis functions are defined. Continuing, by the bilinearity of Equation (4(4) $$\begin{array}{c}\hfill ⟨f,g⟩={\int }_{-1}^1{\int }_0^{z_0}f(z,\mu )g(z,\mu )\mathrm{d}z\mathrm{d}\mu ,\end{array}$$(4) ) and the linearity of Equation (13(13) $$\begin{array}{c}\hfill P(\psi ,{\psi }^{†})=-{\int }_0^1\mu [g_1(\mu ){\psi }^{†}(0,\mu )+g_2(\mu ){\psi }^{†}(z_0,-\mu )]\mathrm{d}\mu .\end{array}$$(13) ), we may rewrite Equation (11(11) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},S⟩-P(\psi ,{\psi }^{†})\end{array}$$(11) ) as(32) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},f⟩-P(\psi ,{\psi }^{†})=⟨{\psi }_h^{†},f⟩+⟨{\psi }_p^{†},f⟩-P(\psi ,{\psi }_h^{†})-P(\psi ,{\psi }_p^{†}).\end{array}$$(32)

We then express the first term on the right-hand side of Equation (32(32) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},f⟩-P(\psi ,{\psi }^{†})=⟨{\psi }_h^{†},f⟩+⟨{\psi }_p^{†},f⟩-P(\psi ,{\psi }_h^{†})-P(\psi ,{\psi }_p^{†}).\end{array}$$(32) ) as (33a) $$\begin{array}{cc}\hfill ⟨{\psi }_h^{†},f⟩& ={\int }_0^{z_0}{\int }_{-1}^1{\psi }_h^{†}(z,\mu )f(z)\mathrm{d}\mu \mathrm{d}z\hfill \end{array}$$(33a) (33b) $$\begin{array}{cc}& ={\int }_0^{z_0}f(z){\int }_0^1\left({\psi }_h^{†}(z,\mu )+{\psi }_h^{†}(z,-\mu )\right)\mathrm{d}\mu \mathrm{d}z\hfill \end{array}$$(33b) (33c) $$\begin{array}{ccc}& \approx \hfill & \hfill {\int }_0^{z_0}f(z)\sum _{k=1}^N{w}_k\left({\psi }_h^{†}(z,{\mu }_k)+{\psi }_h^{†}(z,-{\mu }_k)\right)\mathrm{d}z,\end{array}$$(33c) where, looking back to Equation (27(27) $$\begin{array}{c}\hfill {\mathit{\psi }}_{\pm ,h}^{†}(z)=\sum _{j=1}^N[a_j{\mathit{\varphi }}_{\pm }({\nu }_j)e^{-z/{\nu }_j}+b_j{\mathit{\varphi }}_{\mp }({\nu }_j)e^{-(z_0-z)/{\nu }_j}],\end{array}$$(27) ), we get(34) $$\begin{array}{c}\hfill \sum _{k=1}^N{w}_k\left({\psi }_h^{†}(z,{\mu }_k)+{\psi }_h^{†}(z,-{\mu }_k)\right)=\sum _{j=1}^N\left(a_j{e}^{-z/{\nu }_j}+b_j{e}^{-(z_0-z)/{\nu }_j}\right){\varphi }_0({\nu }_j),\end{array}$$(34)

with(35) $$\begin{array}{c}\hfill {\varphi }_0({\nu }_j)=\sum _{k=1}^N{w}_k\left(\varphi ({\nu }_j,{\mu }_k)+\varphi ({\nu }_j,-{\mu }_k)\right).\end{array}$$(35)

After replacing Equation (34(34) $$\begin{array}{c}\hfill \sum _{k=1}^N{w}_k\left({\psi }_h^{†}(z,{\mu }_k)+{\psi }_h^{†}(z,-{\mu }_k)\right)=\sum _{j=1}^N\left(a_j{e}^{-z/{\nu }_j}+b_j{e}^{-(z_0-z)/{\nu }_j}\right){\varphi }_0({\nu }_j),\end{array}$$(34) ) into Equation (33c(33c) $$\begin{array}{ccc}& \approx \hfill & \hfill {\int }_0^{z_0}f(z)\sum _{k=1}^N{w}_k\left({\psi }_h^{†}(z,{\mu }_k)+{\psi }_h^{†}(z,-{\mu }_k)\right)\mathrm{d}z,\end{array}$$(33c) ) and using the definition of f to restrict the integration limits, we obtain (36a) $$\begin{array}{cc}\hfill ⟨{\psi }_h^{†},f⟩& ={\int }_{z_c}^{z_d}\sum _{j=1}^N\left(a_j{e}^{-z/{\nu }_j}+b_j{e}^{-(z_0-z)/{\nu }_j}\right){\varphi }_0({\nu }_j)\mathrm{d}z\hfill \end{array}$$(36a) (36b) $$\begin{array}{cc}& =\sum _{j=1}^N{\nu }_j\left[a_j\left(e^{-z_c/{\nu }_j}-e^{-z_d/{\nu }_j}\right)\right.\hfill \\ \hfill & +\left.b_j\left(e^{-(z_0-z_d)/{\nu }_j}-e^{-(z_0-z_c)/{\nu }_j}\right)\right]{\varphi }_0({\nu }_j).\hfill \end{array}$$(36b) Similarly, for the second term on the right-hand side of Equation (32(32) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},f⟩-P(\psi ,{\psi }^{†})=⟨{\psi }_h^{†},f⟩+⟨{\psi }_p^{†},f⟩-P(\psi ,{\psi }_h^{†})-P(\psi ,{\psi }_p^{†}).\end{array}$$(32) ), we use Equation (30(30) $$\begin{array}{c}\hfill {\mathit{\psi }}_{\pm ,p}^{†}(z)=\sum _{j=1}^N[a_j(z){\mathit{\varphi }}_{\pm }({\nu }_j)+b_j(z){\mathit{\varphi }}_{\mp }({\nu }_j)]\end{array}$$(30) ) to evaluate (37a) $$\begin{array}{ccc}\hfill ⟨{\psi }_p^{†},f⟩& \approx \hfill & \hfill {\int }_0^{z_0}f(z)\sum _{k=1}^N{w}_k\left({\psi }_p^{†}(z,{\mu }_k)+{\psi }_p^{†}(z,-{\mu }_k)\right)\mathrm{d}z\end{array}$$(37a) (37b) $$\begin{array}{cc}& =\sum _{j=1}^N{\varphi }_0({\nu }_j){\int }_{z_c}^{z_d}\left(a_j(z)+b_j(z)\right)\mathrm{d}z.\hfill \end{array}$$(37b) Now, special attention has to be given to the evaluation of the integrals in the right-hand side. In this way, we first substitute Equation (31a(31a) $$\begin{array}{c}\hfill a_j(z)=c_j{\int }_0^z{S}^{†}(z^{\prime })e^{-(z-z^{\prime })/{\nu }_j}\mathrm{d}{z}^{\prime }\end{array}$$(31a) ) along with the definition of ${\sigma }_d$ into the first term of Equation (37b(37b) $$\begin{array}{cc}& =\sum _{j=1}^N{\varphi }_0({\nu }_j){\int }_{z_c}^{z_d}\left(a_j(z)+b_j(z)\right)\mathrm{d}z.\hfill \end{array}$$(37b) ) (38a) $$\begin{array}{cc}\hfill {\int }_{z_c}^{z_d}{a}_j(z)\mathrm{d}z& =c_j{\int }_{z_c}^{z_d}{e}^{-z/{\nu }_j}{\int }_0^z{\sigma }_d(z^{\prime })e^{z^{\prime }/{\nu }_j}\mathrm{d}{z}^{\prime }\mathrm{d}z\hfill \end{array}$$(38a) (38b) $$\begin{array}{cc}& =c_j{\sigma }_d{\int }_{z_c}^{z_d}{e}^{-z/{\nu }_j}{\int }_{min\{z,z_a\}}^{min\{z,z_b\}}{e}^{z^{\prime }/{\nu }_j}\mathrm{d}{z}^{\prime }\mathrm{d}z\hfill \end{array}$$(38b) (38c) $$\begin{array}{cc}& =c_j{\nu }_j{\sigma }_d{\int }_{z_c}^{z_d}\left(e^{-(z-min\{z,z_b\})/{\nu }_j}-e^{-(z-min\{z,z_a\})/{\nu }_j}\right)\mathrm{d}z,\hfill \end{array}$$(38c) taking into account that $0<z^{\prime }<z$ and ${\sigma }_d(z)=0$ when z is outside $[z_a,z_b]$. Here the notation $min\{z,z_a\}$ means the minimum value between z and $z_a$. Similarly, (39a) $$\begin{array}{cc}\hfill {\int }_{z_c}^{z_d}{b}_j(z)\mathrm{d}z& =c_j{\int }_{z_c}^{z_d}{e}^{z/{\nu }_j}{\int }_z^{z_0}{\sigma }_d(z^{\prime })e^{-z^{\prime }/{\nu }_j}\mathrm{d}{z}^{\prime }\mathrm{d}z\hfill \end{array}$$(39a) (39b) $$\begin{array}{cc}& =c_j{\sigma }_d{\int }_{z_c}^{z_d}{e}^{z/{\nu }_j}{\int }_{max\{z,z_a\}}^{max\{z,z_b\}}{e}^{-z^{\prime }/{\nu }_j}\mathrm{d}{z}^{\prime }\mathrm{d}z\hfill \end{array}$$(39b) (39c) $$\begin{array}{cc}& =c_j{\nu }_j{\sigma }_d{\int }_{z_c}^{z_d}\left(e^{-(max\{z,z_a\}-z)/{\nu }_j}-e^{-(max\{z,z_b\}-z)/{\nu }_j}\right)\mathrm{d}z,\hfill \end{array}$$(39c) where, in this case, $z<z^{\prime }<z_0$ and the maximum value between z and $z_a$ is denoted by $max\{z,z_a\}$.

Proceeding, we evaluate the integral in Equation (38c(38c) $$\begin{array}{cc}& =c_j{\nu }_j{\sigma }_d{\int }_{z_c}^{z_d}\left(e^{-(z-min\{z,z_b\})/{\nu }_j}-e^{-(z-min\{z,z_a\})/{\nu }_j}\right)\mathrm{d}z,\hfill \end{array}$$(38c) ), to finally obtain(40) $$\begin{array}{cc}\hfill {\int }_{z_c}^{z_d}{a}_j(z)\mathrm{d}z& =c_j{\nu }_j{\sigma }_d\left(max\{z_b-z_c,0\}-max\{z_a-z_c,0\}\right.\hfill \\ \hfill & -\left.max\{z_b-z_d,0\}+max\{z_a-z_d,0\}\right)\hfill \\ \hfill & +c_j{\nu }_j^2{\sigma }_d\left(e^{-(max\{z_b,z_c\}-z_b)/{\nu }_j}-e^{-(max\{z_b,z_d\}-z_b)/{\nu }_j}\right)\hfill \\ \hfill & -c_j{\nu }_j^2{\sigma }_d\left(e^{-(max\{z_a,z_c\}-z_a)/{\nu }_j}-e^{-(max\{z_a,z_d\}-z_a)/{\nu }_j}\right),\hfill \end{array}$$(40)

where we have tried to properly take into account all possible relations among the values of $z_a$, $z_b$, $z_c$ and $z_d$, with $z_a\le z_b$ and $z_c\le z_d$, when defining the integral limits. Analogously, for Equation (39c(39c) $$\begin{array}{cc}& =c_j{\nu }_j{\sigma }_d{\int }_{z_c}^{z_d}\left(e^{-(max\{z,z_a\}-z)/{\nu }_j}-e^{-(max\{z,z_b\}-z)/{\nu }_j}\right)\mathrm{d}z,\hfill \end{array}$$(39c) ),(41) $$\begin{array}{cc}\hfill {\int }_{z_c}^{z_d}{b}_j(z)\mathrm{d}z& =c_j{\nu }_j{\sigma }_d\left(max\{z_d-z_a,0\}-max\{z_d-z_b,0\}\right.\hfill \\ \hfill & -\left.max\{z_c-z_a,0\}+max\{z_c-z_b,0\}\right)\hfill \\ \hfill & +c_j{\nu }_j^2{\sigma }_d\left(e^{-(z_a-min\{z_a,z_d\})/{\nu }_j}-e^{-(z_a-min\{z_a,z_c\})/{\nu }_j}\right)\hfill \\ \hfill & -c_j{\nu }_j^2{\sigma }_d\left(e^{-(z_b-min\{z_b,z_d\})/{\nu }_j}-e^{-(z_b-min\{z_b,z_c\})/{\nu }_j}\right).\hfill \end{array}$$(41)

Going back now to Equation (32(32) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},f⟩-P(\psi ,{\psi }^{†})=⟨{\psi }_h^{†},f⟩+⟨{\psi }_p^{†},f⟩-P(\psi ,{\psi }_h^{†})-P(\psi ,{\psi }_p^{†}).\end{array}$$(32) ), we use the previously defined angular quadrature scheme to express the third term in the right-hand side as (42a) $$\begin{array}{cc}\hfill P(\psi ,{\psi }_h^{†})& =-{\int }_0^1\mu \left[g_1(\mu ){\psi }_h^{†}(0,\mu )+g_2(\mu ){\psi }_h^{†}(z_0,-\mu )\right]\mathrm{d}\mu \hfill \end{array}$$(42a) (42b) $$\begin{array}{ccc}& \approx \hfill & \hfill -\sum _{k=1}^N{w}_k{\mu }_k\left[g_1({\mu }_k){\psi }_h^{†}(0,{\mu }_k)+g_2({\mu }_k){\psi }_h^{†}(z_0,-{\mu }_k)\right].\end{array}$$(42b) Replacing the expressions for ${\psi }_h^{†}(0,{\mu }_k)$ and ${\psi }_h^{†}(z_0,-{\mu }_k)$ calculated in Equation (27(27) $$\begin{array}{c}\hfill {\mathit{\psi }}_{\pm ,h}^{†}(z)=\sum _{j=1}^N[a_j{\mathit{\varphi }}_{\pm }({\nu }_j)e^{-z/{\nu }_j}+b_j{\mathit{\varphi }}_{\mp }({\nu }_j)e^{-(z_0-z)/{\nu }_j}],\end{array}$$(27) ) into Equation (42b(42b) $$\begin{array}{ccc}& \approx \hfill & \hfill -\sum _{k=1}^N{w}_k{\mu }_k\left[g_1({\mu }_k){\psi }_h^{†}(0,{\mu }_k)+g_2({\mu }_k){\psi }_h^{†}(z_0,-{\mu }_k)\right].\end{array}$$(42b) ), we get, after some algebraic manipulations, (43a) $$\begin{array}{c}\hfill P(\psi ,{\psi }_h^{†})=\sum _{j=1}^N\left[{\varphi }_1({\nu }_j)+{\varphi }_2({\nu }_j)\right]\end{array}$$(43a)

with(43b) $$\begin{array}{c}\hfill {\varphi }_1({\nu }_j)=-a_j\sum _{k=1}^N{w}_k{\mu }_k\left[g_1({\mu }_k)\varphi ({\nu }_j,{\mu }_k)+g_2({\mu }_k)\varphi ({\nu }_j,-{\mu }_k)e^{-z_0/{\nu }_j}\right]\end{array}$$(43b)

and(43c) $$\begin{array}{c}\hfill {\varphi }_2({\nu }_j)=-b_j\sum _{k=1}^N{w}_k{\mu }_k\left[g_1({\mu }_k)\varphi ({\nu }_j,-{\mu }_k)e^{-z_0/{\nu }_j}+g_2({\mu }_k)\varphi ({\nu }_j,{\mu }_k)\right].\end{array}$$(43c) Finally, the last term in the right-hand side of Equation (32(32) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},f⟩-P(\psi ,{\psi }^{†})=⟨{\psi }_h^{†},f⟩+⟨{\psi }_p^{†},f⟩-P(\psi ,{\psi }_h^{†})-P(\psi ,{\psi }_p^{†}).\end{array}$$(32) ), we express, after similar calculations as (44a) $$\begin{array}{ccc}\hfill P(\psi ,{\psi }_p^{†})& \approx \hfill & \hfill -\sum _{k=1}^N{w}_k{\mu }_k\left[g_1({\mu }_k){\psi }_p^{†}(0,{\mu }_k)+g_2({\mu }_k){\psi }_p^{†}(z_0,-{\mu }_k)\right]\end{array}$$(44a) (44b) $$\begin{array}{cc}& ={\sigma }_d\sum _{j=1}^N{c}_j{\nu }_j\left[{\varphi }_3({\nu }_j)+{\varphi }_4({\nu }_j)\right]\hfill \end{array}$$(44b)

with(44c) $$\begin{array}{c}\hfill {\varphi }_3({\nu }_j)=-\left(e^{-z_a/{\nu }_j}-e^{-z_b/{\nu }_j}\right)\sum _{k=1}^N{w}_k{\mu }_k{g}_1({\mu }_k)\varphi ({\nu }_j,-{\mu }_k)\end{array}$$(44c)

and(44d) $$\begin{array}{c}\hfill {\varphi }_4({\nu }_j)=-\left(e^{-(z_0-z_b)/{\nu }_j}-e^{-(z_0-z_a)/{\nu }_j}\right)\sum _{k=1}^N{w}_k{\mu }_k{g}_2({\mu }_k)\varphi ({\nu }_j,-{\mu }_k).\end{array}$$(44d)

To summarize, to use the adjoint angular flux to evaluate the detector response, through Equation (11(11) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},S⟩-P(\psi ,{\psi }^{†})\end{array}$$(11) ), one can either directly apply numerical quadrature schemes to evaluate the requested integrals, a procedure which in general demands more computational time, or to compute r by Equation (32(32) $$\begin{array}{c}\hfill r=⟨{\psi }^{†},f⟩-P(\psi ,{\psi }^{†})=⟨{\psi }_h^{†},f⟩+⟨{\psi }_p^{†},f⟩-P(\psi ,{\psi }_h^{†})-P(\psi ,{\psi }_p^{†}).\end{array}$$(32) ) along with the explicit expressions derived in this section. In the later case, we start by computing Equation (35(35) $$\begin{array}{c}\hfill {\varphi }_0({\nu }_j)=\sum _{k=1}^N{w}_k\left(\varphi ({\nu }_j,{\mu }_k)+\varphi ({\nu }_j,-{\mu }_k)\right).\end{array}$$(35) ) to substitute into Equation (36b(36b) $$\begin{array}{cc}& =\sum _{j=1}^N{\nu }_j\left[a_j\left(e^{-z_c/{\nu }_j}-e^{-z_d/{\nu }_j}\right)\right.\hfill \\ \hfill & +\left.b_j\left(e^{-(z_0-z_d)/{\nu }_j}-e^{-(z_0-z_c)/{\nu }_j}\right)\right]{\varphi }_0({\nu }_j).\hfill \end{array}$$(36b) ), to obtain $⟨{\psi }_h^{†},f⟩$. Then, we evaluate Equations (38c(38c) $$\begin{array}{cc}& =c_j{\nu }_j{\sigma }_d{\int }_{z_c}^{z_d}\left(e^{-(z-min\{z,z_b\})/{\nu }_j}-e^{-(z-min\{z,z_a\})/{\nu }_j}\right)\mathrm{d}z,\hfill \end{array}$$(38c) ) and (39c(39c) $$\begin{array}{cc}& =c_j{\nu }_j{\sigma }_d{\int }_{z_c}^{z_d}\left(e^{-(max\{z,z_a\}-z)/{\nu }_j}-e^{-(max\{z,z_b\}-z)/{\nu }_j}\right)\mathrm{d}z,\hfill \end{array}$$(39c) ) to substitute into (37b(37b) $$\begin{array}{cc}& =\sum _{j=1}^N{\varphi }_0({\nu }_j){\int }_{z_c}^{z_d}\left(a_j(z)+b_j(z)\right)\mathrm{d}z.\hfill \end{array}$$(37b) ), to calculate $⟨{\psi }_p^{†},f⟩$. Next, we compute $P(\psi ,{\psi }_h^{†})$, Equation (43a(43a) $$\begin{array}{c}\hfill P(\psi ,{\psi }_h^{†})=\sum _{j=1}^N\left[{\varphi }_1({\nu }_j)+{\varphi }_2({\nu }_j)\right]\end{array}$$(43a) ), and $P(\psi ,{\psi }_p^{†})$, Equation (44b(44b) $$\begin{array}{cc}& ={\sigma }_d\sum _{j=1}^N{c}_j{\nu }_j\left[{\varphi }_3({\nu }_j)+{\varphi }_4({\nu }_j)\right]\hfill \end{array}$$(44b) ), to have the final value of r established.

We discuss aspects of our computational implementation and present some numerical results in the next section.

5. Computational aspects and numerical results

To begin, we fix a set of D neutral particle detectors, with ${\sigma }_{d,i}$(z), for $i=1,\cdots ,D$, as in Equation (6(6) $$\begin{array}{c}\hfill {\sigma }_d(z)=\left\{\begin{array}{cc}{\sigma }_d,\hfill & z\in [z_a,z_b],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(6) ). To generate the measurements, we first solve Equation (1(1) $$\begin{array}{c}\hfill \mathcal{L}\psi =S(z)\end{array}$$(1) ), given a source S, and compute the numerical response $r_{0,i}$ of each detector, using Equation (5(5) $$\begin{array}{c}\hfill r=⟨{\sigma }_d,\psi ⟩\end{array}$$(5) ). In this work, we use the ADO [23] method for solving the forward transport equation, to compute $r_{0,i}$. After that, we generate noisy readings by computing(45) $$\begin{array}{c}\hfill r_{m,i}=r_{0,i}+{ϵ}_i\end{array}$$(45)

for $i=1,\cdots ,D$, where ${ϵ}_i\sim \mathcal{N}(0,{\sigma }_i)$ is a random number generated from a zero-mean Gaussian distribution with standard deviation ${\sigma }_i$, where, for the numerical tests, we discuss later on in this text, we assumed two values: (A) ${\sigma }_i=0.01r_{0,i}$; (B) ${\sigma }_i=0.05r_{0,i}$.

Continuing, with respect to the minimization procedure, we rewrite the objective function defined in Equation (21(21) $$\begin{array}{c}\hfill h_{\lambda }(\mathbf{s})=||{\mathbf{r}}^{\prime }-{\mathbf{A}s||}_2^2+{\lambda }^2{||\mathbf{T}s||}_2^2,\end{array}$$(21) ) as(46) $$\begin{array}{c}\hfill h_{\lambda }(\mathbf{s})=||{\mathbf{r}}^{\prime }-{\mathbf{A}s||}_2^2+{\lambda }^2{||\mathbf{T}s||}_2^2=||\left[\begin{array}{c}{\mathit{r}}^{\prime }\\ \mathbf{0}\end{array}\right]-\left[\begin{array}{c}\mathit{A}\\ \lambda \mathit{T}\end{array}\right]\mathit{s}|{|}_2^2\end{array}$$(46)

for a given regularization parameter $\lambda >0$, with $\left[\begin{array}{c}{\mathit{r}}^{\prime }\\ \mathbf{0}\end{array}\right]\in {\mathbb{R}}^{D+B-K}$ and $\left[\begin{array}{c}\mathit{A}\\ \lambda \mathit{T}\end{array}\right]\in {\mathbb{R}}^{D+B-K\times B}$, where $K=0,1,2$ depends on the regularization matrix ${\mathbf{T}}_K$. Provided that only non-negative neutral particle sources are physically acceptable, we minimize Equation (46(46) $$\begin{array}{c}\hfill h_{\lambda }(\mathbf{s})=||{\mathbf{r}}^{\prime }-{\mathbf{A}s||}_2^2+{\lambda }^2{||\mathbf{T}s||}_2^2=||\left[\begin{array}{c}{\mathit{r}}^{\prime }\\ \mathbf{0}\end{array}\right]-\left[\begin{array}{c}\mathit{A}\\ \lambda \mathit{T}\end{array}\right]\mathit{s}|{|}_2^2\end{array}$$(46) ) using the non-negative least squares method [30].

To generate matrix $\mathit{A}$, we apply the ADO method for obtaining the adjoint flux, as described in Section 4, solving Equation (10(10) $$\begin{array}{c}\hfill {\mathcal{L}}^{†}{\psi }^{†}=S^{†}\end{array}$$(10) ) D times, considering as source term the absorption function ${\sigma }_d,i(z)$ of each detector. To compute the entries $A_{i,j}$ of $\mathit{A}$, either numerical integration of the inner-product integrals, Equation (16(16) $$\begin{array}{c}\hfill A_{i,j}=⟨{\psi }_i^{†},f_j⟩.\end{array}$$(16) ), can be used or, our choice here, analytical expressions discussed in Section 4.3; in particular Equations (36b(36b) $$\begin{array}{cc}& =\sum _{j=1}^N{\nu }_j\left[a_j\left(e^{-z_c/{\nu }_j}-e^{-z_d/{\nu }_j}\right)\right.\hfill \\ \hfill & +\left.b_j\left(e^{-(z_0-z_d)/{\nu }_j}-e^{-(z_0-z_c)/{\nu }_j}\right)\right]{\varphi }_0({\nu }_j).\hfill \end{array}$$(36b) ) and (37b(37b) $$\begin{array}{cc}& =\sum _{j=1}^N{\varphi }_0({\nu }_j){\int }_{z_c}^{z_d}\left(a_j(z)+b_j(z)\right)\mathrm{d}z.\hfill \end{array}$$(37b) ). Finally, we have to choose the regularization parameter $\lambda$, in Equation (46(46) $$\begin{array}{c}\hfill h_{\lambda }(\mathbf{s})=||{\mathbf{r}}^{\prime }-{\mathbf{A}s||}_2^2+{\lambda }^2{||\mathbf{T}s||}_2^2=||\left[\begin{array}{c}{\mathit{r}}^{\prime }\\ \mathbf{0}\end{array}\right]-\left[\begin{array}{c}\mathit{A}\\ \lambda \mathit{T}\end{array}\right]\mathit{s}|{|}_2^2\end{array}$$(46) ). First choice was such that $r_i$ in Equation (15(15) $$\begin{array}{c}\hfill r_i=⟨{\psi }_i^{†},\widehat{S}⟩-P(\psi ,{\psi }_i^{†})=\sum _{j=1}^B{s}_j⟨{\psi }_i^{†},f_j⟩-P(\psi ,{\psi }_i^{†})\end{array}$$(15) ) would satisfy the Morozov’s principle [26](47) $$\begin{array}{c}\hfill |r_i-r_{0,i}|<{\sigma }_i,\end{array}$$(47)

for $i=1,\cdots ,D$. However, we found in our tests that Equation (47(47) $$\begin{array}{c}\hfill |r_i-r_{0,i}|<{\sigma }_i,\end{array}$$(47) ) was too restrictive, only being satisfied on a few test cases. Instead, we selected $\lambda$ such that a weaker condition(48) $$\begin{array}{c}\hfill ||\mathit{r}-{\mathit{r}}_m{||}_2{<||\mathit{\sigma }||}_2,\end{array}$$(48)

with $\mathit{\sigma }=[{\sigma }_i]\in {\mathbb{R}}^D$, is satisfied.

In order to analyse the quality of our reconstructions, we consider a few measurements related to the reconstruction error:

(1)	the absolute and relative errors between the original font S and the reconstruction $\widehat{S}$ (49a) $$\begin{array}{c}\hfill {\epsilon }_{S,Abs}=||S-\widehat{S}{||}_2,{\epsilon }_{S,Rel}=\frac{||S-\widehat{S}{||}_2}{{||S||}_2},\end{array}$$(49a)
(2)	the absolute and relative errors between the scalar flux $\varphi$ computed with the original source S and the scalar flux $\widehat{\varphi }$ computed with the reconstruction $\widehat{S}$ (49b) $$\begin{array}{c}\hfill {\epsilon }_{\varphi ,Abs}=||\varphi -\widehat{\varphi }{||}_2,{\epsilon }_{\varphi ,Rel}=\frac{||\varphi -\widehat{\varphi }{||}_2}{{||\varphi ||}_2},\end{array}$$(49b)
(3)	the absolute and relative error between the noisy-free readings ${\mathit{r}}_m$ and the reconstruction $\mathit{r}$ (49c) $$\begin{array}{c}\hfill {\epsilon }_{r,Abs}=||{\mathit{r}}_m-{\mathit{r}||}_2=||{\mathit{r}}^{\mathbf{\prime }}-{\mathit{As}||}_2,{\epsilon }_{r,Rel}=\frac{||{\mathit{r}}^{\mathbf{\prime }}-{\mathit{As}||}_2}{||{\mathit{r}}_m{||}_2},\end{array}$$(49c)

where $||·{||}_2$ is the functional 2-norm for the spatial variable in Equations (49a(49a) $$\begin{array}{c}\hfill {\epsilon }_{S,Abs}=||S-\widehat{S}{||}_2,{\epsilon }_{S,Rel}=\frac{||S-\widehat{S}{||}_2}{{||S||}_2},\end{array}$$(49a) ) and (49b(49b) $$\begin{array}{c}\hfill {\epsilon }_{\varphi ,Abs}=||\varphi -\widehat{\varphi }{||}_2,{\epsilon }_{\varphi ,Rel}=\frac{||\varphi -\widehat{\varphi }{||}_2}{{||\varphi ||}_2},\end{array}$$(49b) ), and the vector 2-norm in Equation (49c(49c) $$\begin{array}{c}\hfill {\epsilon }_{r,Abs}=||{\mathit{r}}_m-{\mathit{r}||}_2=||{\mathit{r}}^{\mathbf{\prime }}-{\mathit{As}||}_2,{\epsilon }_{r,Rel}=\frac{||{\mathit{r}}^{\mathbf{\prime }}-{\mathit{As}||}_2}{||{\mathit{r}}_m{||}_2},\end{array}$$(49c) ). Moreover, we define the scalar flux $\varphi :[0,z_0]\to {\mathbb{R}}_{+}$ as [21](50) $$\begin{array}{c}\hfill \varphi (z)={\int }_{-1}^1\psi (z,\mu )\mathrm{d}\mu .\end{array}$$(50)

All tests were performed on a machine equipped with an Intel Core i5-4670 processor with 16 GiB of RAM.

5.1. Problem I

As a first test, we consider an homogeneous medium defined for $z\in [0,50]$ and we analyse two cases of the parameter c, $c=c_1=0.3$ and $c=c_2=0.9$. To be able to explore higher order anisotropy degrees, we chose a known scattering law in radiative transfer applications [22,31] to which a recurrence relation is available to generate the coefficients(51) $$\begin{array}{c}\hfill {\beta }_l=\left(\frac{2l+1}{2l-1}\right)\left(\frac{L+1-l}{L+1+l}\right){\beta }_{l-1},\end{array}$$(51)

for $l=0,\cdots ,200$, with ${\beta }_0=1$. At the boundaries $z=0$ and $z=50$ the incoming fluxes are zero.

We aim to reconstruct a Gaussian source(52) $$\begin{array}{c}\hfill S(z)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(z-25\right)}^2},\end{array}$$(52)

which may be taken as an approximation to the Dirac’s Delta function centred at $z=25$ [32]. We define an uniform partition of the interval [0, 50] to set the basis functions $\{f_j\}$, Equation (20(20) $$\begin{array}{c}\hfill f_j(z)=\left\{\begin{array}{cc}1,\hfill & z\in [z_{j-1},z_j],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(20) ),(53) $$\begin{array}{c}\hfill [0,50]=\bigcup _{j=1}^{50}[j-1,j].\end{array}$$(53)

This basis is clearly too poor to represent the source S. This aspect, however, may be considered as an additional noise to our reconstruction. We found that 0.1058 is the absolute error in the $L^2$ norm between S and its projection onto the space generated by the proposed basis, which leads to a relative error of $61.25\%$.

In our tests, we initially placed a given number ($D>0$) of detectors, with fixed width, equally distributed in the domain. When different from zero, we set the absorption parameter ${\sigma }_{d,i}=0.1$. With respect to the relation with the basis function support, we consider two configurations in [0,50]: (i) a single particle detector within each basis function support and (ii) two detectors within each support. Still, we consider the support of each detector is defined in an arbitrary sub-interval of just one basis function support. It cannot intercept the support of two basis functions. For both configurations, measurements $r_{0,i}$ were computed with at least 10 decimal places of accuracy, by the ADO method, using 64 discrete directions in the interval (0, 1].

In order to compute the matrix $\mathit{A}$, we solve the adjoint transport equation using each one of the detectors as adjoint source term. The solutions here were computed using the ADO method considering 32 discrete directions with 9 to 10 decimal places of accuracy. With respect to the identifiability [25] of our inverse model, we have $|\mathit{A}{\mathit{A}}^T|\approx 0$ for both configurations, indicating the necessity of regularization techniques.

Since we hope at least to be able to approximate the original source S when using noise-free readings, we plot in Figures 2 and 3 the true source S (dashed line), the projection of S onto the span of $\{f_j\}$ (circles) and the reconstruction without regularization using noise-free readings (asterisks) for configurations (i) and (ii) and each choice of c. Table 1 shows the errors we found in this reconstruction.

Figure 2. Exact source S (dashed line), projection of S onto the span of $\{f_j\}$ (circles) and unregularized reconstruction using noise-free readings (asterisks) using 50 particle detectors.

Figure 3. Exact source S (dashed line), projection of S onto the span of $\{f_j\}$ (circles) and unregularized reconstruction using noise-free readings (asterisks) using 100 particle detectors.

Table 1. Absolute error and relative to S error for the noise-free reconstructions, angular flux and detector reading.

	Absolute error		Relative error $(\%)$
$c\backslash$ # detectors	50	100	50	100
	Source error
0.3	0.1087	0.1093	62.07	62.25
0.9	0.1092	0.1099	62.23	62.41
	Scalar flux error
0.3	0.0494	0.0484	1.02	3.74
0.9	0.0539	0.0525	1.11	1.08
	Detector reading error
0.3	$1.6674(-6)$	$3.8726(-4)$	0.01	0.87
0.9	$9.2555(-6)$	$4.2272(-4)$	0.02	0.92

As indicated by Figures 2 and 3, and Table 1, using error-free measurements, the derived method was able to approximate S successfully, in the sense of reproducing a relative error of similar magnitude of its projection onto the basis space.

As a next step, we use Equation (45(45) $$\begin{array}{c}\hfill r_{m,i}=r_{0,i}+{ϵ}_i\end{array}$$(45) ) to add noise to each detector reading. As mentioned previously, we work with cases (A) and (B), respectively, where we set the standard deviation for the normal distribution as $1\%$ and $5\%$ of each reading. Table 2 summarize our test cases.

Table 2. Test cases for problem I.

Perturbation	A		B
c$\backslash$# detectors	50	100	50	100
0.3	I	III	V	VII
0.9	II	IV	VI	VIII

Table 3. Comparisons between exact and perturbed readings.

		Absolute error		Relative error (%)
c$\backslash$# detectors		50	100	50	100
	A
4*0.3	Min	2.0016($-$5)	3.2058($-$5)	1.5482($-$1)	2.6299($-$1)
	Mea	1.1896($-$4)	1.1765($-$4)	9.2018($-$1)	9.6513($-$1)
	Max	3.3041($-$4)	2.5404($-$4)	2.5557($+$0)	2.0840($+$0)
	Std	4.7092($-$5)	3.2667($-$5)	3.6426($-$1)	2.6798($-$1)
4*0.9	Min	1.3413($-$4)	2.2129($-$4)	2.7623($-$1)	4.8336($-$1)
	Mea	4.6931($-$4)	4.4909($-$4)	9.6650($-$1)	9.8095($-$1)
	Max	1.1512($-$3)	8.1245($-$4)	2.3709($+$0)	1.7746($+$0)
	Std	1.2570($-$4)	8.7235($-$5)	2.5886($-$1)	1.9055($-$1)
	B
4*0.3	Min	7.2294($-$5)	1.0652($-$4)	5.5919($-$1)	8.7385($-$1)
	Mea	6.0541($-$4)	5.8953($-$4)	4.6828(+0)	4.8362(+0)
	Max	1.8575($-$3)	1.2368($-$3)	1.4368(+1)	1.0146(+1)
	Std	2.3667($-$4)	1.6388($-$4)	1.8306(+0)	1.3444(+0)
4*0.9	Min	8.9836($-$4)	9.8951($-$4)	1.8501($+$0)	2.1614($+$0)
	Mea	2.3436($-$3)	2.2495($-$3)	4.8264($+$0)	4.9135($+$0)
	Max	5.0127($-$3)	4.0539($-$3)	1.0323(+1)	8.8548($+$0)
	Std	6.2387($-$4)	4.2961($-$4)	1.2848($+$0)	9.3840($-$1)

We list in Table 3 the maximum value, mean, minimum and standard deviation of the perturbed readings for both configurations of the detectors and different values of c. In Table 4, we present the absolute and relative errors between the exact source S and the reconstructed source $\widehat{S}$, calculated by Equation (49a(49a) $$\begin{array}{c}\hfill {\epsilon }_{S,Abs}=||S-\widehat{S}{||}_2,{\epsilon }_{S,Rel}=\frac{||S-\widehat{S}{||}_2}{{||S||}_2},\end{array}$$(49a) ); the absolute and relative errors between the scalar flux calculated using the exact source S and the scalar flux calculated using the reconstructed source $\widehat{S}$, through Equation (49b(49b) $$\begin{array}{c}\hfill {\epsilon }_{\varphi ,Abs}=||\varphi -\widehat{\varphi }{||}_2,{\epsilon }_{\varphi ,Rel}=\frac{||\varphi -\widehat{\varphi }{||}_2}{{||\varphi ||}_2},\end{array}$$(49b) ), and, the absolute and relative errors between the detector readings calculated using the exact source S and the readings computed using the reconstructed source $\widehat{S}$, as in Equation (49c(49c) $$\begin{array}{c}\hfill {\epsilon }_{r,Abs}=||{\mathit{r}}_m-{\mathit{r}||}_2=||{\mathit{r}}^{\mathbf{\prime }}-{\mathit{As}||}_2,{\epsilon }_{r,Rel}=\frac{||{\mathit{r}}^{\mathbf{\prime }}-{\mathit{As}||}_2}{||{\mathit{r}}_m{||}_2},\end{array}$$(49c) ).

Table 4. Problem I: minimum (Min), mean (Mea), maximum (Max) and standard deviation (Std) values for (i) the absolute and relative to S errors of the reconstructed source $\widehat{S}$; (ii) absolute and relative to $\varphi$ error of the scalar flux calculated using the reconstructed source $\widehat{S}$; (iii) absolute and relative to r errors of the detector readings. Results obtained utilizing ${\mathit{T}}_{\mathbf{0}}$ as Tikhonov’s regularization matrix. Test cases described in Table 2.

Test		${ϵ}_{S,Abs}$	${ϵ}_{S,Rel}$ (%)	${ϵ}_{\varphi ,Abs}$	${ϵ}_{\varphi ,Rel}$ (%)	${ϵ}_{r,Abs}$	${ϵ}_{r,Rel}$ (%)
3*I	Min	1.0630($-$1)	6.1387(+1)	4.7640($-$2)	3.6845($+$0)	1.3609($-$5)	1.0526($-$1)
	Mea	1.0822($-$1)	6.1938(+1)	5.3132($-$2)	4.1092($+$0)	1.4355($-$4)	1.1104($+$0)
	Max	1.1229($-$1)	6.3091(+1)	6.6858($-$2)	5.1707($+$0)	3.8869($-$4)	3.0065($+$0)
	Std	8.7269($-$4)	2.4940($-$1)	2.9597($-$3)	2.2890($-$1)	5.7173($-$5)	4.4223($-$1)
3*II	Min	1.0597($-$1)	6.1291(+1)	5.2735($-$2)	1.0860($+$0)	9.7835($-$5)	2.0148($-$1)
	Mea	1.1016($-$1)	6.2487(+1)	7.3209($-$2)	1.5076($+$0)	4.5396($-$4)	9.3489($-$1)
	Max	1.2374($-$1)	6.6231(+1)	1.3696($-$1)	2.8205($+$0)	1.2257($-$3)	2.5243($+$0)
	Std	2.5694($-$3)	7.2360($-$1)	1.1411($-$2)	2.3499($-$1)	1.4515($-$4)	2.9891($-$1)
3*III	Min	1.0732($-$1)	6.1680(+1)	4.7474($-$2)	3.6716($+$0)	3.8742($-$4)	3.1782($+$0)
	Mea	1.0952($-$1)	6.2307(+1)	4.9181($-$2)	3.8036($+$0)	3.9675($-$4)	3.2548($+$0)
	Max	1.1280($-$1)	6.3235(+1)	5.5454($-$2)	4.2888($+$0)	4.5470($-$4)	3.7301($+$0)
	Std	8.3447($-$4)	2.3718($-$1)	1.1216($-$3)	8.6740($-$2)	7.3909($-$6)	6.0632($-$2)
3*IV	Min	1.0611($-$1)	6.1331(+1)	5.1774($-$2)	1.0662($+$0)	4.2641($-$4)	9.3140($-$1)
	Mea	1.1121($-$1)	6.2784(+1)	6.1162($-$2)	1.2595($+$0)	5.1995($-$4)	1.1357($+$0)
	Max	1.2510($-$1)	6.6595(+1)	1.1110($-$1)	2.2879($+$0)	9.5003($-$4)	2.0751($+$0)
	Std	2.3678($-$3)	6.6553($-$1)	6.3606($-$3)	1.3099($-$1)	5.7187($-$5)	1.2491($-$1)
3*V	Min	1.0607($-$1)	6.1320(+1)	4.8383($-$2)	3.7419($+$0)	8.8744($-$5)	6.8643($-$1)
	Mea	1.1452($-$1)	6.3689(+1)	8.4007($-$2)	6.4970($+$0)	6.5251($-$4)	5.0471($+$0)
	Max	1.8061($-$1)	8.0016(+1)	1.8789($-$1)	1.4531(+1)	1.7818($-$3)	1.3782(+1)
	Std	6.7598($-$3)	1.8363($+$0)	2.2124($-$2)	1.7110($+$0)	2.6062($-$4)	2.0159($+$0)
3*VI	Min	1.0700($-$1)	6.1589(+1)	6.1972($-$2)	1.2762($+$0)	3.5917($-$4)	7.3966($-$1)
	Mea	1.4053($-$1)	7.0367(+1)	2.0632($-$1)	4.2487($+$0)	2.0245($-$3)	4.1692($+$0)
	Max	2.6146($-$1)	9.6273(+1)	5.9550($-$1)	1.2263(+1)	5.8991($-$3)	1.2149(+1)
	Std	2.2543($-$2)	5.4788($+$0)	8.3547($-$2)	1.7205($+$0)	8.4428($-$4)	1.7387($+$0)
3*VII	Min	1.0606($-$1)	6.1316(+1)	4.7581($-$2)	3.6799($+$0)	3.9119($-$4)	3.2091($+$0)
	Mea	1.1289($-$1)	6.3244(+1)	6.9039($-$2)	5.3394($+$0)	6.0223($-$4)	4.9404($+$0)
	Max	1.4842($-$1)	7.2536(+1)	1.5115($-$1)	1.1690(+1)	1.3858($-$3)	1.1369(+1)
	Std	5.2983($-$3)	1.4598($+$0)	1.4813($-$2)	1.1456($+$0)	1.3858($-$4)	1.1368($+$0)
3*VIII	Min	1.0633($-$1)	6.1395(+1)	5.7059($-$2)	1.1750($+$0)	4.7954($-$4)	1.0474($+$0)
	Mea	1.2654($-$1)	6.6857(+1)	1.5776($-$1)	3.2487($+$0)	1.4673($-$3)	3.2050($+$0)
	Max	2.5823($-$1)	9.5677(+1)	3.7573($-$1)	7.7375($+$0)	3.5173($-$3)	7.6827($+$0)
	Std	1.5740($-$2)	3.9759($+$0)	4.7458($-$2)	9.7732($-$1)	4.4757($-$4)	9.7761($-$1)

As Table 4 indicates, for the source term reconstruction, we were able to obtain for cases I to IV (less noise) errors of similar order of magnitude as the ones obtained with no perturbation, showed in Table 1.

The errors associated with the particle detector readings may also be considered satisfactory. On the other hand, for cases V to VIII, when we increase the noise, the errors in the source reconstruction are a lot higher, even reaching $96.27\%$ in test VI. However, we were still able to obtain satisfactory results for the error in the scalar flux and detector readings.

Proceeding, we repeat our tests using ${\mathit{T}}_1$, defined in Equation (22a(22a) $$\begin{array}{c}\hfill {\mathbf{T}}_1=\left(\begin{array}{cccc}-1& 1& & \\ & \ddots & \ddots & \\ & & -1& 1\end{array}\right)\end{array}$$(22a) ), as regularization matrix. The errors are very similar to the ones obtained in the previous tests using ${\mathit{T}}_{\mathbf{0}}$; however, we remark a drop between 50% and 60% in computational time, since the regularization parameter is more quickly determined.

Using ${\mathit{T}}_2$ as regularization matrix, Equation (22b(22b) $$\begin{array}{c}\hfill {\mathbf{T}}_2=\left(\begin{array}{ccccc}1& -2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& -2& 1\end{array}\right).\end{array}$$(22b) ), we also found similar results to the errors when comparing to the previous tests using ${\mathit{T}}_{\mathbf{0}}$ and ${\mathit{T}}_{\mathbf{1}}$ as regularization matrices. Similarly with the test using ${\mathit{T}}_1$, the computational time was once again faster than the test using ${\mathit{T}}_0$ as regularization matrix.

We also note, that based on the analysis of the results, the additional computational effort required to run the 100 detectors case do not seem to enhance the reconstruction process (cases III, IV, VII and VIII).

We show in Figures 4 the best and worst reconstructions obtained in our tests with higher noise, $c=0.9$, and ${\mathit{T}}_{\mathbf{0}}$ as regularization matrix. We named the best reconstruction the one which absolute error between the exact and the estimated source is minimum, and similarly, the worst reconstruction is the one with maximum absolute error.

Figure 4. (a) Exact source S (dashed line), best reconstruction (circles), and worst reconstruction (asterisks); (b) Scalar flux calculated using the exact source S (dashed line), using the best reconstruction (dotted line) and using the worst reconstruction (solid line).

In Figure 4(a), both the best and worst reconstruction were able to approximate the exact source shape and the location of the region where it is positive. The scalar flux computed with the best reconstruction almost overlapped, as Figure 4(b) indicates, the exact one. On the other hand, for the worst reconstruction, the scalar flux presented higher values in the same spot where the reconstructed source was far from the exact one.

For our next tests, we refine our partition. We set our basis functions $\{f_j\}$ in Equation (20(20) $$\begin{array}{c}\hfill f_j(z)=\left\{\begin{array}{cc}1,\hfill & z\in [z_{j-1},z_j],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(20) ) using an uniform partition(54) $$\begin{array}{c}\hfill [0,50]=\bigcup _{j=1}^{100}\left[\frac{1}{2}(j-1),\frac{1}{2}j\right]\end{array}$$(54)

of [0, 50]. Using these basis functions, we found an absolute error of 0.0539 between S and its projection into the basis space, and a relative error of $43.70\%$. For the readings, we consider 100 particle detectors, one for each basis function, all with the absorption parameter ${\sigma }_{d,i}=0.1$ within its support. We note, however, that the support of the detector is defined as a subinterval of the basis function support (in general we use about $10\%$ of the interval).

In the same way as before, we compute the real measurements $r_{0,i}$ through solving the transport equation by the ADO method considering 64 discrete directions, obtaining in the process at least 10 decimal places. The inversion matrix $\mathit{A}$ is once again calculated using the ADO method considering 32 discrete directions when solving the adjoint transport equation, which was also able to attain at least 10 decimal places of accuracy.

We plot in Figure 5 the true source S (dashed line), the projection of S onto the span of $\{f_j\}$ (circles) and the reconstruction without regularization using noise-free readings (asterisks) for the two values of c.

Figure 5. Exact source S (dashed line), projection of S onto the span of $\{f_j\}$ (circles) and unregularized reconstruction using noise-free readings (asterisks) using 100 basis functions and 100 particle detectors.

As indicated by Figure 5, the derived method was able to successfully approximate the desired true source S using error-free measurements, once again, in the sense that relative errors between the original source and the reconstructed one are similar to the ones between the original source and its projection.

As we have already done in the previous tests, we use Equation (45(45) $$\begin{array}{c}\hfill r_{m,i}=r_{0,i}+{ϵ}_i\end{array}$$(45) ) to add noisy information to each exact detector reading, then proceed with the reconstructions for different values of c and each regularization matrix.

Similarly as the previous configuration, using noiseless data, we were able to obtain good approximations to the true source. The errors in the readings were also satisfactory. As the noise level increases, we noted again higher errors when estimating the source. But, one more time, we were still able to obtain relatively good results to the scalar flux and detector readings estimates. In Figure 6, we show the best and worst reconstructions using ${\mathit{T}}_{\mathbf{0}}$ as regularization matrix and $c=0.9$, under the conditions just described.

Figure 6. (a) Exact source S (dashed line), best reconstruction (circles), and worst reconstruction (asterisks); (b) Scalar flux calculated using the exact source S (dashed line), using the best reconstruction (dotted line) and using the worst reconstruction (solid line).

Although there is a clear lack in the reconstruction quality for the worst case, Figure 6(a), the scalar flux calculated using this reconstructed source, showed in Figure 6(b), almost overlaps the scalar flux computed using the exact source. Moreover, for almost the totality of the samples (4816 out of 5000), the error ${ϵ}_{S,Abs}\in [\mu -2\sigma ,\mu +2\sigma ]$, with $\mu$ and $\sigma$ being, respectively, the mean and standard deviation of the sample distribution. Restricting the samples to those within this interval, we find a maximum absolute error of $1.3751(-1)$, which is equivalent to a $69.82\%$ relative error to the exact source S error, a value close to the ones obtained using 50 basis functions. We plot in Figure 7 the graphs of the source and scalar flux reconstructions for this case.

Figure 7. (a) Exact source S (dashed line), best (circles) and worst (asterisks) reconstructions within two standard deviations interval to the absolute error mean; (b) Scalar flux calculated using the exact source S (dashed line), using the best (dotted line) and worst reconstructions (solid line) within two standard deviations interval to the absolute error mean.

5.2. Problem II

As a second test problem, we consider the reconstruction of a localized source piecewisely defined for $z\in [0,10]$ by(55) $$\begin{array}{c}\hfill S(z)=\left\{\begin{array}{cc}0.75,\hfill & z\in [4.0,4.7),\hfill \\ 1,\hfill & z\in [4.7,6.2),\hfill \\ 0.25,\hfill & z\in [6.2,7.0],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(55)

We assume that the particles are isotropically migrating within the slab, a medium with physical parameters set as $c=0.99$ and ${\beta }_0=1$. As before, we also assume that there is no incoming flux at the boundaries $z=0$ and $z=10$.

For this test, we assume three different partitions of [0, 10] when defining the basis functions in Equation (20(20) $$\begin{array}{c}\hfill f_j(z)=\left\{\begin{array}{cc}1,\hfill & z\in [z_{j-1},z_j],\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.\end{array}$$(20) ): (56a) $$\begin{array}{cc}\hfill P_1& =\bigcup _{j=1}^{20}[0.5(j-1),0.5j],\hfill \end{array}$$(56a) (56b) $$\begin{array}{cc}\hfill P_2& =\bigcup _{j=1}^{40}[0.25(j-1),0.25j],\hfill \end{array}$$(56b)

and(56c) $$\begin{array}{cc}\hfill P_3& =\left(\bigcup _{j=1}^8[0.5(j-1),0.5j]\right)\bigcup \left(\bigcup _{j=1}^{12}[0.25j+3.75,0.25j+4]\right)\hfill \\ \hfill & \bigcup \left(\bigcup _{j=1}^6[0.5j+6.5,0.5j+7]\right)=\bigcup _{j=1}^{26}[z_{j-1},z_j].\hfill \end{array}$$(56c)

At this time, a detector is placed at the centre of each basis function support, in general occupying about $10\%$ of the support interval, with the absorption parameter defined by ${\sigma }_{d,i}=0.35$ in the detector support. Just as before, for each detector, we compute a reading $r_{0,i}$ using Equation (5(5) $$\begin{array}{c}\hfill r=⟨{\sigma }_d,\psi ⟩\end{array}$$(5) ) and, thereafter, we add normally distributed noise, with zero mean and $5\%$ of each reading as standard deviation. For this test problem, we used only ${\mathit{T}}_1$, defined in Equation (22a(22a) $$\begin{array}{c}\hfill {\mathbf{T}}_1=\left(\begin{array}{cccc}-1& 1& & \\ & \ddots & \ddots & \\ & & -1& 1\end{array}\right)\end{array}$$(22a) ), as regularization matrix and consider 5000 noisy samples. Table 5 displays all errors obtained in the simulation.

Table 5. Problem II: minimum (Min), mean (Mea), maximum (Max) and standard deviation (Std) values for (i) the absolute and relative to S errors of the reconstructed source $\widehat{S}$; (ii) absolute and relative to $\varphi$ error of the scalar flux calculated using the reconstructed source $\widehat{S}$; (iii) absolute and relative to r errors of the detector readings. Results obtained utilizing ${\mathit{T}}_{\mathbf{1}}$ as Tikhonov’s regularization matrix. Number of detectors equal to the number of basis functions for each partition.

Partition		${ϵ}_{S,Abs}$	${ϵ}_{S,Rel}$ (%)	${ϵ}_{\varphi ,Abs}$	${ϵ}_{\varphi ,Rel}$ (%)	${ϵ}_{r,Abs}$	${ϵ}_{r,Rel}$ (%)
4*$P_1$	Min	$2.3377(-1)$	$1.7216(+1)$	$6.0061(+0)$	$1.1841(+1)$	$1.4872(-1)$	$1.1845(+1)$
	Mea	$4.2186(-1)$	$3.1068(+1)$	$8.5285(+0)$	$1.6813(+1)$	$2.1114(-1)$	$1.6816(+1)$
	Max	$8.5100(-1)$	$6.2673(+1)$	$1.1905(+1)$	$2.3471(+1)$	$2.9470(-1)$	$2.3472(+1)$
	Std	$8.2019(-2)$	$6.0404(+0)$	$8.3088(-1)$	$1.6380(+0)$	$2.0559(-2)$	$1.6375(+0)$
4*$P_2$	Min	$3.0231(-1)$	$2.2264(+1)$	$6.5285(+0)$	$1.2870(+1)$	$1.1426(-1)$	$1.2872(+1)$
	Mea	$4.0231(-1)$	$2.9629(+1)$	$8.4779(+0)$	$1.6713(+1)$	$1.4837(-1)$	$1.6714(+1)$
	Max	$6.3641(-1)$	$4.6869(+1)$	$1.0549(+1)$	$2.0797(+1)$	$1.8462(-1)$	$2.0797(+1)$
	Std	$4.3632(-2)$	$3.2133(+0)$	$5.8572(-1)$	$1.1547(+0)$	$1.0250(-2)$	$1.1546(+0)$
4*$P_3$	Min	$3.0819(-1)$	$2.2697(+1)$	$5.5612(+0)$	$1.0963(+1)$	$1.0737(-1)$	$1.0312(+1)$
	Mea	$3.8339(-1)$	$3.0445(+1)$	$8.3075(+0)$	$1.6378(+1)$	$1.2421(-1)$	$1.1929(+1)$
	Max	$6.2880(-1)$	$4.6309(+1)$	$1.0288(+1)$	$2.0282(+1)$	$1.8170(-1)$	$1.7450(+1)$
	Std	$4.2910(-2)$	$3.1601(+0)$	$6.6038(-1)$	$1.3019(+0)$	$0.9301(-2)$	$8.9327(-1)$

Similarly to the previous test problem, we were able to obtain, even when adding more noise, satisfactory approximations to the original source S. Upon using the finer partition $P_2$ instead of $P_1$, the method presented a slightly improvement of the error values. The reconstructions based on the $P_3$ partition provided better results than the ones using $P_2$, probably as consequence of the combination of a finer discretization when defining our basis function and smaller number of detectors adding noise information into the system. In Figure 8, we plot the best reconstruction result using the partition $P_3$.

Figure 8. (a) Exact source S (dashed line) and best reconstruction using the partition $P_3$ (circles); (b) Scalar flux calculated using the exact source S (solid line) and using the reconstruction (dashed line).

As both Figures, 8(a) and (b), indicate, our method is able to provide a qualitative idea of the true particle source S shape and magnitude, although, as numerical simulations indicate, it was not able to reproduce the original source as accurately as the Gaussian source in the previous test problem.

As a variation of the previous test, we now reduce the number of particle detectors, still uniformly distributed along the physical domain, to only 10. Partitions $P_1$ and $P_2$ are considered for defining the basis functions. We show in Table 6 the results.

Table 6. Problem II: minimum (Min), mean (Mea), maximum (Max) and standard deviation (Std) values for (i) the absolute and relative to S errors of the reconstructed source $\widehat{S}$; (ii) absolute and relative to $\varphi$ error of the scalar flux calculated using the reconstructed source $\widehat{S}$; (iii) absolute and relative to r errors of the detector readings. Results obtained utilizing ${\mathit{T}}_{\mathbf{1}}$ as Tikhonov’s regularization matrix. Uniform distribution of ten detectors along the domain.

Partition		${ϵ}_{S,Abs}$	${ϵ}_{S,Rel}$ (%)	${ϵ}_{\varphi ,Abs}$	${ϵ}_{\varphi ,Rel}$ (%)	${ϵ}_{r,Abs}$	${ϵ}_{r,Rel}$ (%)
4*$P_1$	Min	$3.1922(-1)$	$2.3509(+1)$	$1.5326(+0)$	$3.0214(+0)$	$7.4221(-3)$	$2.0983(-1)$
	Mea	$1.0589(+0)$	$7.7987(+1)$	$2.1770(+0)$	$4.2911(+0)$	$3.0374(-2)$	$8.5870(-1)$
	Max	$1.6399(+0)$	$1.2077(+2)$	$2.8483(+0)$	$5.6152(+0)$	$7.3452(-2)$	$2.0765(+0)$
	Std	$2.2969(-1)$	$1.6916(+1)$	$1.9316(-1)$	$3.8082(-1)$	$9.8640(-3)$	$2.7886(-1)$
4*$P_2$	Min	$3.6717(-1)$	$2.7040(+1)$	$1.6109(+0)$	$3.1758(+0)$	$6.5031(-3)$	$1.8384(-1)$
	Mea	$1.3753(+0)$	$1.0128(+2)$	$2.1961(+0)$	$4.3295(+0)$	$3.1419(-2)$	$8.8821(-1)$
	Max	$2.2420(+1)$	$1.6512(+2)$	$2.9022(+0)$	$5.7215(+0)$	$7.9257(-2)$	$2.2406(+0)$
	Std	$4.0663(-1)$	$2.9947(+1)$	$1.8317(-1)$	$3.6110(-1)$	$9.6237(-3)$	$2.7207(-1)$

As the analysis of Table 6 indicates, the reduction of the number of detectors to half and one quarter, respectively, of the number of basis functions, led to higher errors in the source reconstruction.

We also tested the reconstruction with a nonuniform distribution of detectors throughout the domain. We first considered 20 particle detectors placed near to the boundaries – half of them within the range [0, 1], and the remaining ones in [9, 10]. We still assume noisy measurements normally distributed, with zero mean and $5\%$ of each reading as standard deviation.

Table 7. Problem II: minimum (Min), mean (Mea), maximum (Max) and standard deviation (Std) values for (i) the absolute and relative to S errors of the reconstructed source $\widehat{S}$; (ii) absolute and relative to $\varphi$ error of the scalar flux calculated using the reconstructed source $\widehat{S}$; (iii) absolute and relative to r errors of the detector readings. Results obtained utilizing ${\mathit{T}}_{\mathbf{1}}$ as Tikhonov’s regularization matrix. Twenty detectors placed (half of them) close to each one of the boundaries of the domain.

Partition		${ϵ}_{S,Abs}$	${ϵ}_{S,Rel}$ (%)	${ϵ}_{\varphi ,Abs}$	${ϵ}_{\varphi ,Rel}$ (%)	${ϵ}_{r,Abs}$	${ϵ}_{r,Rel}$ (%)
4*$P_1$	Min	$2.9801(-1)$	$2.1947(+1)$	$1.9221(+0)$	$3.8681(+0)$	$1.4783(-4)$	$9.2035(-2)$
	Mea	$1.8879(+0)$	$1.3903(+2)$	$2.6329(+1)$	$5.1906(+1)$	$9.2840(-4)$	$5.7799(-1)$
	Max	$2.4158(+0)$	$1.7792(+2)$	$3.4605(+1)$	$6.8220(+1)$	$2.1780(-3)$	$1.3559(+0)$
	Std	$3.2429(-1)$	$2.3883(+1)$	$7.0962(+0)$	$1.3990(+1)$	$2.8343(-4)$	$1.7646(-1)$
4*$P_2$	Min	$7.4106(-1)$	$5.4576(+1)$	$3.4435(+0)$	$6.7887(+0)$	$2.5367(-4)$	$1.5792(-1)$
	Mea	$2.2179(+0)$	$1.6334(+2)$	$3.0290(+1)$	$5.9716(+1)$	$1.1573(-3)$	$7.2049(-1)$
	Max	$3.2934(+0)$	$2.4254(+2)$	$3.6712(+1)$	$7.2376(+1)$	$2.2418(-3)$	$1.3956(+0)$
	Std	$2.3199(-1)$	$1.7085(+1)$	$6.5066(+0)$	$1.2827(+1)$	$2.8321(-4)$	$1.7631(-1)$

As we see in Table 7, using partition $P_1$ to generate the reconstruction basis, situation where we have as many detectors as basis functions, the errors were much higher than the ones we found using uniformly distributed detectors. Proceeding, using $P_2$, where we have twice as many basis functions as detectors, the error values we obtained were even higher.

Proceeding, we then consider 20 particle detectors now placed all of them within the range [4, 7], exactly the region where the source S is positive.

Table 8. Problem II: minimum (Min), mean (Mea), maximum (Max) and standard deviation (Std) values for (i) the absolute and relative to S errors of the reconstructed source $\widehat{S}$; (ii) absolute and relative to $\varphi$ error of the scalar flux calculated using the reconstructed source $\widehat{S}$; (iii) absolute and relative to r errors of the detector readings. Results obtained utilizing ${\mathit{T}}_{\mathbf{1}}$ as Tikhonov’s regularization matrix. Twenty detectors distributed in a subinterval of the domain.

Partition		${ϵ}_{S,Abs}$	${ϵ}_{S,Rel}$ (%)	${ϵ}_{\varphi ,Abs}$	${ϵ}_{\varphi ,Rel}$ (%)	${ϵ}_{r,Abs}$	${ϵ}_{r,Rel}$ (%)
4*$P_1$	Min	$2.4867(-1)$	$1.8314(+1)$	$1.3703(+0)$	$2.7015(+0)$	$3.9344(-3)$	$3.6219(-1)$
	Mea	$4.9092(-1)$	$3.6154(+1)$	$2.1978(+0)$	$4.3329(+0)$	$6.8760(-3)$	$6.3298(-1)$
	Max	$2.4726(+0)$	$1.8210(+2)$	$9.3313(+0)$	$1.8396(+1)$	$1.2980(-2)$	$1.1949(+0)$
	Std	$1.6400(-1)$	$1.2078(+1)$	$5.4522(-1)$	$1.0749(+0)$	$1.4082(-3)$	$1.2963(-1)$
4*$P_2$	Min	$3.6541(-1)$	$2.6911(+1)$	$1.3900(+0)$	$2.7404(+0)$	$2.7180(-3)$	$2.5020(-1)$
	Mea	$7.6422(-1)$	$5.6282(+1)$	$2.0860(+0)$	$4.1124(+0)$	$7.6468(-3)$	$7.0394(-1)$
	Max	$1.3314(+0)$	$9.8054(+1)$	$3.5097(+0)$	$6.9190(+0)$	$1.3599(-2)$	$1.2519(+0)$
	Std	$1.4551(-1)$	$1.0716(+1)$	$2.8490(-1)$	$5.6166(-1)$	$1.5916(-3)$	$1.4652(-1)$

Table 8 indicates that the errors, in the source reconstruction, are slightly higher than the ones we found using uniformly distributed detectors. However, results obtained through this configuration seems better than the previous one, where the detectors were concentrated closer to the boundaries. In fact, we note the shape of the source still may be recovered.

6. Discussion

The ADO method, an analytical approach to solve the discrete ordinates approximation of the adjoint transport equation, was applied in the inverse problem of estimating the location and magnitude of sources of neutral particles. The use of analytical expressions when calculating the matrix $\mathit{A}$, basic for the minimization problem, Equation (18(18) $$\begin{array}{c}\hfill h(\mathbf{s})=||{\mathbf{r}}^{\prime }-\mathbf{A}s{||}_2^2\end{array}$$(18) ), allowed us to avoid all spatial discretizations error inherent of quadrature schemes. We remark, however, that those analytical expressions are only valid for the particular choice of basis functions and absorption function used in this work, requiring to be re-evaluated for distinct configurations. In both test problems, each reconstruction took less than a second in our machine. In the first test problem, when computing the matrix $\mathit{A}$, our code took an average of 9.445 s to perform all the calculations using quadrature schemes. In contrast, using the analytical expressions proposed in Section 4.3, we observed a reduction in the computational time to 0.007 s, a gain of three orders of magnitude. We should stress here that although we have carefully implemented the computational code for both versions, the codes were not particularly optimized, still having room for improvements.

For the two source reconstruction model problems considered, with localized sources and detectors scattered through the domain, we were able to obtain good results for the absolute error between the reconstructed source and the original one: values of order ${10}^{-2}$ for the Gaussian source and order ${10}^{-1}$ for the localized piecewise constant source. Relative errors were in general of same order of magnitude than the ones between the original source and its projection. Still, scalar fluxes and the detector readings recovered with the reconstructions were very satisfactory. In the reconstruction of the Gaussian source, we found similar errors using ${\mathit{T}}_0$, ${\mathit{T}}_1$ and ${\mathit{T}}_2$ as regularization matrices; however, we were able to get an appropriate regularization parameter faster when using ${\mathit{T}}_1$ and ${\mathit{T}}_2$. For the second test problem, when using a more refined partition only in the interval where the source is positive, we were able to easily achieve same order of magnitude error than the other partitions, using smaller number of basis functions and detectors. In regard to the detectors location, although the cases where they were uniformly distributed through the domain result in smaller errors, it may not be realistic in some applications. When they were concentrated in some intervals of the domain, the configuration with detectors closer to the boundaries shown to be worse than the one where they were located in coincidence with the positive source region.

Perhaps an obvious improvement is the use of others basis functions for the reconstruction space, like polynomials and splines, which may be considered in the future, still requiring our analytical expressions to be re-evaluated for such cases. However, it is our main idea to move, firstly, in the direction of extending our analysis to more complex model problems.

Currently, we are working on the use of Bayesian framework in the inverse process as well as the inclusion of energy dependence in the model problem. Furthermore, studies developed by our research group found out that the ADO method along with nodal techniques, when applied in bi-dimensional problems, are able to reproduce results as accurate as other techniques from coarser meshes [20]. It might provide a speed up in computational time and suggests a natural extension of this work to two-dimensional problems.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by CAPES and CNPq of Brazil [grant number 308233/2014-0].