Stability of Stationary Inverse Transport Equation in Diffusion Scaling

TL;DR

This paper investigates how the stability of reconstructing optical parameters from the stationary radiative transfer equation deteriorates as the Knudsen number approaches zero, revealing increased ill-posedness in the diffusive limit.

Contribution

It provides a quantitative analysis of the stability loss in inverse transport problems as the system transitions to the diffusive regime, highlighting the challenges in parameter reconstruction.

Findings

Discrepancy in measurements is amplified by a factor of $K_n^p$ in the reconstructed parameters.

Inverse problem becomes increasingly ill-posed as $K_n o 0$.

Theoretical results are validated with numerical tests.

Abstract

We consider the inverse problem of reconstructing the optical parameters for stationary radiative transfer equation (RTE) from velocity-averaged measurement. The RTE often contains multiple scales characterized by the magnitude of a dimensionless parameter---the Knudsen number ($K_n$). In the diffusive scaling ($K_n \ll 1$), the stationary RTE is well approximated by an elliptic equation in the forward setting. However, the inverse problem for the elliptic equation is acknowledged to be severely ill-posed as compared to the well-posedness of inverse transport equation, which raises the question of how uniqueness being lost as $K_n \rightarrow 0$. We tackle this problem by examining the stability of inverse problem with varying $K_n$. We show that, the discrepancy in two measurements is amplified in the reconstructed parameters at the order of $K_n^p~ (p = 1\text{ or} ~2)$, and as a result lead to ill-posedness in the zero limit of $K_n$. Our results apply to both continuous and discrete settings. Some numerical tests are performed in the end to validate these theoretical findings.

Equations

$$v\cdot {\nabla }_{x}f=\int k(x,v,v^{\prime })f(x,v^{\prime })\mathrm{d}{v}^{\prime }-\sigma (x,v)f.$$

$${\Gamma }_{\pm }=\{(x,v):x\in \partial \Omega ,\pm v\cdot n_x>0\},$$

$$f{|}_{{\Gamma }_{-}}=\varphi (x,v).$$

$${\sigma }_a(x,v)=\sigma (x,v)-{\int }_{{\mathbb{S}}^{n-1}}k(x,v,v^{\prime })\mathrm{d}{v}^{\prime },$$

$$-C{\nabla }_x\cdot \left(\frac{1}{{\sigma }_s}{\nabla }_{x}u\right)+{\sigma }_{a}u=0,x\in \Omega$$

$$\{\begin{array}{l}v\cdot {\nabla }_{x}f=\frac{1}{\mathsf{Kn}}{\sigma }_s\mathcal{L}f-\mathsf{Kn}{\sigma }_{a}f,\\ f{|}_{{\Gamma }_{-}}=\varphi (x,v),\end{array}$$

$$\mathcal{L}f=\int f(x,v^{\prime })\mathrm{d}{v}^{\prime }-f=⟨f⟩-f.$$

$$\{\begin{array}{l}C{\nabla }_x\cdot \left(\frac{1}{{\sigma }_s}{\nabla }_x\rho \right)-{\sigma }_a\rho =0,\\ \rho {|}_{\partial \Omega }={\xi }_f,\end{array}$$

$${\xi }_f(x_0)=f_{z\to \infty }^l.$$

$$\{\begin{array}{l}{v}_z{\partial }_z{f}^l={\sigma }_s\mathcal{L}[f^l],z\in [0,\infty )\\ f^l{|}_{z=0}=\varphi (x_0,v).\end{array}$$

$$f(x,v)=f_{\text{bd}}+f_{\text{in}},$$

$$f_{\text{in}}=f_0+\mathsf{Kn}{f}_1+{\mathsf{Kn}}^2{f}_2+\cdots .$$

$$\int v\cdot {\nabla }_x\left(-\frac{v}{{\sigma }_s}\cdot {\nabla }_x\rho \right)\mathrm{d}v=-{\sigma }_a\rho \Rightarrow C{\nabla }_x\cdot \left(\frac{1}{{\sigma }_s}{\nabla }_x\rho \right)={\sigma }_a\rho .$$

$$\{\begin{array}{l}{v}_z{\partial }_{z}f={\sigma }_s\mathcal{L}[f],z\in [0,\infty )\\ f{|}_{z=0}=\varphi (x_0,v).\end{array}$$

$$\rho (x_0)=f_{z\to \infty }.$$

$$\Vert f-\rho {\Vert }_{L_{\infty }(\mathrm{d}{x}_i\mathrm{d}v)}=\mathcal{O}({\mathsf{Kn}}^{2/3})\to 0,$$

$$v\cdot {\nabla }_{x}f(x,v)=\frac{{\sigma }_s(x_3)}{\mathsf{Kn}}\mathcal{L}f(x,v)-\mathsf{Kn}{\sigma }_a(x_3)f(x,v)$$

$$\mathcal{L}f(x,v)=\frac{1}{4\pi }{\int }_0^{\pi }{\int }_0^{2\pi }[f(x,v^{\prime })-f(x,v)]\sin {\theta }^{\prime }\mathrm{d}{\varphi }^{\prime }\mathrm{d}{\theta }^{\prime },$$

$$\{\begin{array}{l}\cos \theta {\partial }_{x_3}f(x_3,\cos \theta )=\frac{{\sigma }_s(x_3)}{\mathsf{Kn}}(\frac{1}{4\pi }{\int }_0^{\pi }{\int }_0^{2\pi }(f(x_3,\cos {\theta }^{\prime })-f(x_3,\cos \theta ))\sin {\theta }^{\prime }\mathrm{d}{\varphi }^{\prime }\mathrm{d}{\theta }^{\prime })-\mathsf{Kn}{\sigma }_a(x_3)f(x_3,\cos \theta ),\\ f{|}_{{\Gamma }_{-}}=\varphi (x_3,\cos \theta ).\end{array}$$

$$\{\begin{array}{l}v{\partial }_{x}f(x,v)=\frac{{\sigma }_s(x)}{\mathsf{Kn}}{\int }_{-1}^{1}f(x,v^{\prime })-f(x,v)\frac{\mathrm{d}{v}^{\prime }}{2}-\mathsf{Kn}{\sigma }_a(x)f(x,v),(x,v)\in (0,1)\times [-1,1],\\ f{|}_{{\Gamma }^{-}}=\varphi (x,v),\end{array}$$

$$\mathcal{M}f(x)={\int }_{{\Gamma }_{+}(x)}v\cdot n(x)f(x,v)\mathrm{d}v,$$

$${\mathcal{A}}_{\mathsf{Kn}}({\sigma }_a,{\sigma }_s):\varphi {|}_{{\Gamma }_{-}}\to \mathcal{M}f,$$

$$\{\begin{array}{ll}v\cdot {\nabla }_{x}f& =\frac{1}{\mathsf{Kn}}\mathcal{L}f-\mathsf{Kn}{\sigma }_{a}f,(x,v)\in \Omega \times \mathbb{S},\\ f{|}_{{\Gamma }_{-}}& =\varphi .\end{array}$$

$$\mathcal{A}({\sigma }_a):\varphi \to \mathcal{M}f(\cdot ).$$

$${\tilde{\sigma }}_a(x)={\sigma }_a(x)-{\sigma }_{a0}(x)$$

$$\{\begin{array}{l}v\cdot {\nabla }_x{f}_0=\frac{1}{\mathsf{Kn}}\mathcal{L}{f}_0-\mathsf{Kn}{\sigma }_{a0}{f}_0,\\ f_0{|}_{{\Gamma }_{-}}=\varphi ,\end{array}$$

$$\tilde{f}(x,v)=f(x,v)-f_0(x,v)$$

$$\{\begin{array}{l}v\cdot {\nabla }_x\tilde{f}=\frac{1}{\mathsf{Kn}}\mathcal{L}\tilde{f}-\mathsf{Kn}{\sigma }_{a0}\tilde{f}-\mathsf{Kn}{\tilde{\sigma }}_a{f}_0,\\ \tilde{f}{|}_{{\Gamma }_{-}}=0,\end{array}$$

$$\{\begin{array}{l}-v\cdot {\nabla }_{x}g=\frac{1}{\mathsf{Kn}}\mathcal{L}g-\mathsf{Kn}{\sigma }_{a0}g,\\ g{|}_{{\Gamma }_{+}}={\delta }_y(x).\end{array}$$

$$\int v\cdot n\tilde{f}g{|}_{{\Gamma }_{+}}\mathrm{d}v\mathrm{d}x=\mathsf{Kn}{\int }_{\Omega }{\tilde{\sigma }}_a\int f_{0}g\mathrm{d}v\mathrm{d}x,$$

Full text

All three authors thank the two anonymous referees for the very careful reading of the paper. It leads to significant improvement. The work of K.C. and Q. L. is supported in part by a start-up fund of Q.L. from UW-Madison and National Science Foundation under the grant DMS-1619778. The work of L.W. is supported in part by a start-up fund from SUNY Buffalo and the National Science Foundation under the grant DMS-1620135.

1. Introduction

In this paper, we study the stability of inverse stationary radiative transfer equation (RTE) in different regimes. RTE is a stereotype kinetic equation that describes the dynamics of photon particles in materials with various optical properties [12]. The optical properties are characterized by two parameters—the scattering coefficient and absorption coefficient. Generally speaking, we denote $f(x,v)$ the distribution of particles at location $x$ moving with velocity $v$, and it obeys:

[TABLE]

Here $x\in\Omega\subset\mathbb{R}^{d}$ with $d=2,3$ depending on the dimension of the problem, and $v\in\mathbb{S}^{d-1}$, the unit sphere in $\mathbb{R}^{d}$. It indicates that the particles move with fixed speed (unified to $1$) and therefore has one fewer dimension than $x$. $k(x,v,v^{\prime})$ is termed the scattering coefficient, representing the probability of particles that move in direction $v^{\prime}$ changing to direction $v$. $\mathrm{d}{v}$ has the normalized unit measure. $\sigma(x,v)$ is the total absorption coefficient that represents certain amount of photon particles being absorbed and scattered by the material. The boundary condition is typically imposed as Dirichlet type. We separate the “out-going” and “in-coming” part of boundary by defining:

[TABLE]

where $n_{x}$ is the normal direction pointing out of $\Omega$ at point $x\in\partial\Omega$. In this way, $\Gamma_{-}$ collects all boundary coordinates that represent particles coming into the domain where $\Gamma_{+}$ collects the opposite. It is on the first set we impose Dirichlet boundary condition:

[TABLE]

The well-posedness is summarized from [16] and the solution is proved to be unique with full Dirichlet data.

Although the investigation of RTE is already enormous and generally acknowledged to be well-understood a long time ago, the inverse stationary RTE still attracted lots of attention in the past decade. The new life of RTE lies in the booming of medical techniques and a vast of medical images that require mathematical interpretation. In diffuse optical tomography, for example, near infra-red light (NIR) are sent into biological tissues, and by measuring the outgoing photon current at the surfaces of the tissues, scientists expect to “invert” the problem for the optical properties of the tissue. Mathematically that means one adjusts the incoming data $f|_{\Gamma_{-}}$, and measure a certain form of the out-going data $f|_{\Gamma_{+}}$ seeking for $k(x,v,v^{\prime})$ and $\sigma(x,v)$. Technically it is equivalent to seek for scattering coefficient $k(x,v,v^{\prime})$ and absorption coefficient defined as

[TABLE]

which represents the probability of particles being absorbed by the material only. A vast literature has addressed the problem from various perspectives. The well-posedness of the inverse problem in a generic setting with $\sigma_{a}(x)$ independent of $v$ was addressed in a pioneering paper [15], and the uniqueness based on gauge-invarience when $\sigma_{a}(x,v)$ presents $v$ dependence was shown in [31]. The idea was to decompose the albedo operator according to the singularities. Another approach is to linerize the equation before applying inverse Born series, and show the convergence of the series [26]. The results on the stability of the “inverse” dates back to [33] and was made systematic in [6, 7, 9]. Many papers concern the time-dependent case and the associated stability analysis has also been conducted [21, 30, 14], and also [5] for a review.

Aside from the analytical studies, various numerical techniques are explored accordingly. Numerical treatments could be separated into two categories depending on whether linearization is conducted, both on the original physical domain, or via Green’s function representation [25, 26]. Either way, the resulting numerical problem is typically not well-posed: it is either under-determined or over-determined. The ill-posedness may be inherited from the continuous problem or due, in part, to the lack of data in experiments. The latter reason induces a purely numerical problem that can often be handled via optimization along with some regularization technique. Examples include the standard $L_{2}$ regularization [1] for the smallness, TV regularization [32] for the least variance, $H_{1}$ norm for some regularity, and $L_{1}$ regularization [28] for sparsity, or Tikhonov type on each element in the inverse Born series [26]. The optimization techniques are borrowed accordingly, and efficiency and memory cost for both the Jacobian-type method and the Gradient based method have been compared. The way to set the regularization coefficient, on the other hand, is usually guided by the tolerance on the error [18, 17] and convergence speed.

One very interesting phenomenon associated with inverse stationary transport equation is its connection to the diffusion limit. It has been longly known in the area [29, 2, 3] that the diffusion approximation could serve as a substitute under certain scaling, and instead of inverting the transport equation, one studies the Calderón-like problems. That is, recovering the diffusion and attenuation coefficients $\frac{1}{\sigma_{s}}$ and $\sigma_{a}$ in the following equation

[TABLE]

using the Dirichlet boundary condition $u|_{\partial\Omega}$ and the measurement of flux at boundaries $\frac{\partial u}{\partial n}|_{\partial\Omega}$. Here $\sigma_{s}$ and $\sigma_{a}$ are parameters derived from $\sigma$ and $k$ from the RTE equation and $C$ is a constant that only depends on the dimension. The ill-posedness of the Calderón problem has been shown [19, 20] and its failure in capturing RTE based phenomenon in some scenario has been demonstrated [3], or in a similar problem on diffusion approximation in recovering the doping profile in the Boltzmann-Poisson system [13].

Despite the popularity of both problems, to the best knowledge of the authors, there has not been much study on exploring the connections between the two [4]. Specifically, some questions need to be addressed, such as: when and to what extent can diffusion approximation be used for RTE-based inverse problem? Will such approximation affect the stability in the inverse problem? Considering one type is ill-posed while the other one is well-posed, what is lost when such approximation is performed? In this paper, we make a first attempt to tackle these questions. In particular, we adopt a linearized framework (detailed in the next section), and study the well-posedness and the stability issues when passing to the diffusion limit for three different scenarios: recovering absorption coefficient, and recovering scattering coefficient in both critical and subcritical cases. We would like to mention that there are works on the change of stability with respect to certain parameters in the equation. In [8, 10] the authors particularly studied the stability of the inversion with respect to the modulation frequency in time-harmonic setting, and found that the increasing of the frequency brings more details in the recovery. In [27] the authors studied the stability of recovering acoustic equation.

We emphasize that the current paper concentrates on the illposedness as the transport equation approaches the diffusion regime (losing the stability). Except one special example in 1D (Theorem 5), we assume injectivity in our results. This is not a bizarre assumption as injectivity is shown for the associated nonlinear version of the problem [15]. Indeed, injectivity and stability are two separate issues in inverse problem: the former one concerns the uniqueness in recovery whereas the latter measures the accuracy in recovery when small perturbation in measurement is allowed. In terms of the spectral theory, the injectivity requires that the spectrum is away from zero, and stability studies the whole span of the spectrum. While stability naturally being the next step after injectivity, it is not uncommon that stability can be studied by assuming injectivity. As pointed out in [5], stability is studied in the setting of isotropic source and angularly averaged measurements while the uniqueness is not available [9, 8, 10].

The rest of paper is organized as follows. We present some preliminaries in the next section, including the derivation of the diffusion equation from the RTE and the set-up of the inverse problem in full generality. Section 3 and 4 are devoted to the three scenarios described above respectively. In all three cases, we utilize the linearization approach, study the well-posedness of the problem in both regimes, and examine the change of stability while passing to the diffusion limit. We also introduce a distinguishability parameter to indicate the stability and we justify that the inverse problem becomes more and more indistinguishable in the diffusion limit. Numerical tests are exploited to demonstrate the statements on the properties.

2. Preliminaries

Some preliminaries are collected in this section. The first subsection demonstrates the derivation of the diffusion equation from the RTE in the forward problem and the second subsection sets up the inverse problem we study in a general framework. For the conciseness of the paper we assume $\sigma_{s}(x)$ and $\sigma_{a}(x)$ do not have $v$ dependence.

2.1. Diffusion limit

The diffusion limit sets in when scattering is strong and absorption is weak. The equation in the dimensionless form reads:

[TABLE]

where $\mathcal{L}$ is the collision operator and in the velocity independent case it writes as:

[TABLE]

Here we have re-grouped the gain term and loss term in the collision for the ease of later presentation. A more general collision takes the form $\mathcal{L}f=\int k(v^{\prime},v)(f(x,v^{\prime})-f(x,v))\mathrm{d}{v}^{\prime}$, but our analysis in the rest of the paper can be easily adapted to this case. We therefore keep it in the simplest form, and assume that $\sigma_{s}$ has no $v$ dependence. There are two key features of the collision operator:

• •

Mass conservation: $\int\mathcal{L}[f]\mathrm{d}{v}=0$. If we apply this property to the original equation, immediately we see that $\int vf\mathrm{d}{v}$ is a divergence free field if $\sigma_{a}=0$;

• •

One dimensional Null space: By setting $\mathcal{L}[f]=0$, one gets $f=\langle f\rangle$, meaning that $f$ is a constant in velocity domain. We denote it as ${\rm Null}\,\mathcal{L}=\{\rho(x)\}$, the collection of functions that depend on $x$ only. This property is unique for RTE compared to other linear kinetic equation and it is the main reason that the asymptotic limit only consist of a scalar equation instead of a system.

As $\mathsf{Kn}\to 0$, the equation falls into the diffusion limit and we have the following theorem:

Theorem 1**.**

Suppose $f$ solves (5). As $\mathsf{Kn}\to 0$, $f(x,v)$ converges to $\rho(x)$, which solves the diffusion equation:

[TABLE]

Here $C$ is a constant depending on the dimension of the problem. The boundary condition is determined by:

[TABLE]

with $f^{l}$ solving

[TABLE]

Proof.

The proof follows the standard asymptotic expansion with boundary layer analysis. In the zero limit of $\mathsf{Kn}$, the distribution in the interior will stabilize whereas the boundary condition $\phi$ being away from the equilibrium function will prevent the solution converging near the boundary. To separate the two, we first expand the solution:

[TABLE]

where $f_{\text{bd}}$ is the solution adjacent to boundary accounting for the boundary layer, while $f_{\text{in}}$ characterizes the interior away from the layer. We study $f_{\text{in}}$ first. As $\mathsf{Kn}\rightarrow 0$, we apply the standard asymptotic expansion technique and write:

[TABLE]

Here we only consider the expansion away from the boundary layer so that $f_{\text{bd}}$ is negligible. Inserting the expansion in the equation (5) and equate like powers of $\mathsf{Kn}$:

• $\mathcal{O}(1)$

$\mathcal{L}f_{0}=0$. This immediately indicates that $f_{0}\in{\rm Null}\,\mathcal{L}$. With the form given in (6), ${\rm Null}\,\mathcal{L}$ consists functions that are constants in $v$ domain, and thus $f_{0}(x,v)=\rho(x)$.

• $\mathcal{O}(\mathsf{Kn})$

$v\cdot\nabla_{x}f_{0}=\sigma_{s}\mathcal{L}[f_{1}]$. This indicates that $f_{1}=\mathcal{L}^{-1}\left(v\cdot\nabla_{x}f_{0}\right)$. $\mathcal{L}$ is not a one-to-one map unless the domain is confined in ${\rm Null}\,\mathcal{L}^{\perp}$, and the inverse on $\mathcal{L}$ is pseudo-inverse. Considering the form of $\mathcal{L}$ in (6), then ${\rm Null}\,\mathcal{L}^{\perp}=\{f:\int f\mathrm{d}{v}=0\}$, and thus $f_{1}=-\frac{v}{\sigma_{s}}\cdot\nabla_{x}\rho$.

• $\mathcal{O}(\mathsf{Kn}^{2})$

$v\cdot\nabla_{x}f_{1}=\sigma_{s}\mathcal{L}[f_{2}]-\sigma_{a}f_{0}$. Here we integrate the equation with respect to $v$. The second term will vanish and the left hand side becomes:

[TABLE]

Here the constant $C$ depends on the dimension of the velocity space.

Summarizing up the analysis we obtain $f_{\text{in}}(x,v)\to f_{0}(x,v)=\rho(x)$ that solves the diffusion equation (7).

We then need to provide $\rho$ a correct boundary condition and this comes from the treatment of $f_{\text{bd}}$. For this we follow  [11]. At each point $x\in\partial\Omega$, we perform tangential approximation and by stretching coordinates we can locally change the problem into a half space problem. More specifically, let $x_{0}\in\partial\Omega$ and $n_{x}$ the normal direction pointing out of the domain, we denote $z=-\frac{n_{x}\cdot(x-x_{0})}{\mathsf{Kn}}$. For every fixed point $x=x_{0}-n_{x}y(y\in[0,\infty))$ away from the layer along the ray pointing into the domain, $z\to\infty$ as $\mathsf{Kn}\to 0$. After moving the coordinate frame to $x_{0}$ with $n_{x}$ direction, $z=0$ stands for boundary point and $z=\infty$ is mapped to the interior. Then along $z$ direction the equation reads, in the leading order of $\mathsf{Kn}$:

[TABLE]

It is a half space problem in $z$ with boundary condition given only at $z=0$. The problem is proved to have a unique solution and is computed in [23] and the infinite data on $z$ will be a constant in $v$ direction, which will be used to serve as the Dirichlet boundary condition for $\rho$, meaning:

[TABLE]

This is done at each grid point along the boundary and we end up with the boundary condition for $x\in\partial\Omega$, denoted by $\xi_{f}[\phi](x)$ (or $\xi_{f}(x)$ for short). ∎

Remark 1**.**

The proof here is formal and is not specific on certain norm. In fact with general geometry and boundary condition it is believed to be correct but not proved yet. In 2-D physical domain and 1-D velocity sphere, due to the joint force of [34, 22], it can be made rigorous in $L_{\infty}(\mathrm{d}{x}\mathrm{d}{v})$ norm when boundary layer is excluded, meaning that as $\mathsf{Kn}\to 0$:

[TABLE]

and the decay from $\mathcal{O}(\mathsf{Kn})$ to $\mathcal{O}(\mathsf{Kn}^{2/3})$ is mainly due to the curvature correction, which is controllably small but nontrivial. Here $L_{\infty}(\mathrm{d}{x_{i}}\mathrm{d}{v})$ stands for $L_{\infty}$ norm in the interior only. We exclude a fixed small boundary layer of $\mathcal{O}(\mathsf{Kn})$ width.

Remark 2**.**

Very frequently, RTE in 3D is simplified under some symmetry assumption. Specifically, we assume it with slab plan geometry. Denote $x=(x_{1},x_{2},x_{3})\in\Omega=\mathbb{R}^{2}\times(0,1)$ and $v=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)\in\mathcal{S}^{2}$, then we assume that all parameters and conditions are homogenized along $x_{3}$ direction, i.e. $\sigma_{a}(x)=\sigma_{a}(x_{3})$, $\sigma_{s}(x)=\sigma_{s}(x_{3})$, and that boundary condition is $\phi(x,v)=\phi(x_{3},\cos\theta)$. The stationary RTE becomes:

[TABLE]

with:

[TABLE]

and the boundary condition $f|_{\Gamma_{-}}=\phi(x_{3},\cos\theta)$. With $v\cdot\nabla_{x}=\cos\theta\partial_{x_{3}}$, the stationary RTE becomes:

[TABLE]

For simplicity, we denote $x=x_{3}\in(0,1)$ and make change of variable $v=\cos\theta\in(-1,1)$ to obtain:

[TABLE]

where $v$ is often termed as direction of flight. Sending $\mathsf{Kn}\to 0$ and follow the same asymptotic derivation in the theorem above, one obtains $C=\frac{1}{3}$ in the zero Knudsen number limit in (10). From here on, we always refer to equation (15) as the 1-D RTE. Notice here the velocity domain is interval $[-1,1]$.

2.2. Inverse problem

The inverse problem can be set up associated with a map from the input data on one portion of the boundary to the measurement on the other portion. For RTE specifically, this map is often termed the albedo operator. Depending on the data-acquisition method in the experiments, the measurement can take various forms. Here we assume that only velocity-averaged measurement is available and define the measurement operator:

[TABLE]

where $\Gamma_{+}(x)$ is the “outgoing” semisphere at $x$ defined in (1). And the albedo operator reads:

[TABLE]

where $\phi$ is the Dirichlet boundary condition (2) and $\mathcal{M}f(x)$ is the intensity of light propagating out of the domain at boundary point $x\in\partial\Omega$. Then the inverse problem is to recover $\sigma_{a}$ and $\sigma_{s}$ given the information of the map $\mathcal{A}_{\mathsf{Kn}}$.

On the theoretical level, one concerns about the well-posedness and stability. The well-posedness problem states the following: given the full information on the map $\mathcal{A}$, can one uniquely recover $\sigma_{a}$ and $\sigma_{s}$? The answer is positive if a velocity-resolved measurement (i.e., $f|_{\Gamma_{+}}$ for any $v$ instead of (17)) or time-dependent measurement is available [15, 31]. The analysis is based on singularity separation. Since the albedo operator is a forward map, an explicit form can be obtained and it consists of three parts, separated according to their singular level: the most singular part is a delta function that could be used to recover $\sigma_{a}$ through the inverse X-ray transform, and the secondly singular term is used to recover $\sigma_{s}$, leaving the third term in $L_{\infty}$. The stability problem, on the other hand, asks: if the entire map is off from the accurate one by a small amount of error ($\|\mathcal{A}-\tilde{\mathcal{A}}\|<\epsilon$), how accurate the recovering could be? That is, will $\|\sigma_{a}-\tilde{\sigma}_{a}\|$ remain small? The problem is examined in [7, 6, 9] using the same kind of singular decomposition.

On the numerical level, the well-posedness problem takes a slightly different form: let $\sigma_{s,a}$ be discretized at $N_{x}$ grid points, then if for each incoming data $\phi_{d}$, one could take $N_{p}$ measurements at the boundary, how many incoming data is needed to fully recover $\sigma_{s,a}$? This is intrinsically the same as in the theoretical level. The stability problem, however, has two sides. One is, assume that the data obtained is exact, how much error one obtains in recovering $\sigma_{s,a}$ if $\mathsf{A}$ gets perturbed a little. Here $\mathsf{A}$ is the matrix representation of $\mathcal{A}$. This amounts to analyzing the condition number of the problem. Another is that if the measurement in (17) is off by some error, how accurate will $\sigma_{s,a}$ be? The answer to that lies in analyzing the norm of $\mathsf{A}^{-1}$. The first question is aligned with theoretical stability stated above, while the second question is a pure numerical issue. Concisely, since solving the inverse problem numerically often involves inverting a matrix (or a series of matrices in the nonlinear setting), even if the problem is well-posed, the numerical error is still hard to guarantee. All these questions are indirectly or partially resolved in many papers [2, 28].

We ask a different question in this paper. Our aim is to investigate the dependence of stability and condition number for the inverse problem on the Knudsen number. It is generally acknowledged that the inverse transport problem is well-posed (for most kinds of measurements and under mild assumptions on the scattering/absorption coefficients) whereas the inverse diffusion equation as a limit is ill-posed. Considering the two sets of equations are connected by simply passing to the $\mathsf{Kn}\to 0$ limit, is there a more explicit explanation of the loss of uniqueness?

To answer this question, we adopt a linearization framework [28]. By assuming that the to-be-recovered coefficients $\sigma_{a}$ and $\sigma_{s}$ are only slightly deviated from some given functions from a priori knowledge, we linearize the transport equation around them and write down the relationship between the unknown parameters and the map $\mathcal{A}$. As a result and as will be more clear later, we only need to invert a Fredholm operator of the first kind to recover those parameters. More importantly, this Fredholm operator reveals an explicit dependence on the Knudsen number, which allow us to analyze the stability of the inversion when varying the magnitude of $\mathsf{Kn}$. We would also like to point out that the measurement we took in (17) contains minimum information that no uniqueness results on the reconstruction based on the singular decomposition are available. Nevertheless, the linearization approach we take allows us to quantify the stability explicitly in $\mathsf{Kn}$, and is amenable for numerical schemes. Details will be provided along the paper.

As discussed before, depending on the assumptions on the media, different scenarios could take place. We present here three examples. The following sections are respectively devoted to recovering $\sigma_{a}$ when $\sigma_{s}$ is known, and recovering $\sigma_{s}$ in the critical ($\sigma_{a}=0$) and subcritical ($\sigma_{a}>0$) cases, respectively. It is not our aim in this paper to give a complete analysis for all scenarios but rather to investigate the problem for the first time and nail down the techniques that could make it possible.

3. Recover Absorption Coefficient $\sigma_{a}$

3.1. Inverse problem set-up

In this section, we assume that the scattering coefficient is known (and we set it as $\sigma_{s}=1$ for simplicity), and recover the absorption coefficient $\sigma_{a}$. The equation (5) reads:

[TABLE]

Here $\phi$ is the inflow boundary condition and the solution is denoted by $f(x,v;\phi)$. Experimentally suppose one could measure data at one grid point $y\in\partial\Omega$, then the mapping becomes:

[TABLE]

We follow a linearization framework [28] and set a background absorbing coefficient $\sigma_{a0}(x)$ by assuming that the residue

[TABLE]

is much smaller than $\sigma_{a}$: $\|\tilde{\sigma}_{a}\|\ll\|\sigma_{a}\|$. Then the linearized problem with the same inflow boundary condition reads as

[TABLE]

where $f_{0}(x,v;\phi)$ is the function we linearize upon. Let

[TABLE]

be the fluctuation, then it solves the following equation by subtracting (19) from (18)

[TABLE]

where we have omitted the higher order terms. The incoming boundary information $\phi$ is implicitly contained in $f_{0}$. To make use of the boundary condition and measurement, we write the adjoint problem of (19) and assign it a Delta function boundary condition:

[TABLE]

Multiply the above equation by $\tilde{f}$ and subtract it from the product of (20) and $g$, and integrate in both $x$ and $v$, we get

[TABLE]

which defines a linear mapping from $\phi$ and $\delta_{y}$ with the solutions of (19) (21) to the measured data. Note that the LHS of (22) could be easily obtained by subtracting the computed flux of $f_{0}$ from the measurement, i.e.,

[TABLE]

and the RHS is a Fredholm operator of the first kind with known $\int f_{0}g\mathrm{d}{v}$ and unknown $\tilde{\sigma}_{a}$. Denote

[TABLE]

then the map $\mathcal{A}(\tilde{\sigma}_{a}):~{}\phi\rightarrow b(\delta_{y},\phi)$ rewrites as

[TABLE]

Now it amounts to invert the First type Fredholm integral to recover $\tilde{\sigma}_{a}$.

3.2. Ill-conditioning in the diffusion limit (continuous level)

We study the stability of the inverse problem in terms of the Knudsen number in this subsection. More specifically, consider the linear mapping

[TABLE]

where $\gamma_{\mathsf{Kn}}$ is defined in (23) with $f_{0}$ and $g$ solving (19) and (21), respectively. We aim to understand the influence in recovering $\tilde{\sigma}_{a}$ if a small purturbation in $b(\delta_{y},\phi)$ is introduced. In particular, we would like to check such sensitivity’s dependence on the Knudsen number $\mathsf{Kn}$. To see this, we first define a distinguishability coefficient to measure the “condition number” for a given error $\delta$ on data.

Definition 1**.**

Consider linear mapping in (25), we define the distinguishability coefficient as

[TABLE]

where

[TABLE]

and $\tilde{\sigma}_{a}$ is the solution to (25).

As written, $\Gamma_{\delta}$ is the collection of all possible solutions to the map given that the measurement is contaminated by $\delta$ error. Noticed that both $\langle\gamma_{\mathsf{Kn}},\sigma_{a}\rangle_{L^{2}(\mathrm{d}{x})}$ and $b(\delta_{y},\phi)$ are linearly dependent on $\phi$, we take the sup-norm to normalize in the definition of $\Gamma_{\delta}$. Then $\kappa_{a}$ measures the relative error that could be seen in the recovery — smaller $\kappa_{a}$ leads to better distinguishability. Recall here that the stability defined in [5] says if the two measurement are distinguished by $\delta$, i.e., $\|\mathcal{A}\sigma_{a}-\mathcal{A}\tilde{\sigma}_{a}\|\leq\delta$, then the discrepancy in parameters can be bounded as $\|\sigma_{a}-\tilde{\sigma}_{a}\|\leq\|\mathcal{A}^{-1}\|\|\mathcal{A}(\sigma_{a}-\tilde{\sigma}_{a})\|=\|\mathcal{A}^{-1}\|\delta$, then $\kappa_{a}$ defined in (26) is seen as an estimate of $\|\mathcal{A}^{-1}\|$.

We expect to show two things: 1) smaller error tolerance $\delta$ results in better distinguishability; 2) smaller $\mathsf{Kn}$ drives the problem into the diffusion limit, leading to worse distinguishability. The following theorem groups the two things together.

Theorem 2**.**

We study the inverse problem of recovering $\tilde{\sigma}_{a}$ in  (25). Assume the map from $\tilde{\sigma}_{a}$ to $b$ is injective, given an error tolerance $\delta$ on the measurement, then the distinguishability coefficient grows as $\delta$ grows and $\mathsf{Kn}$ shrinks in the sense that there exists a constant $C$ such that

[TABLE]

Proof.

Let $c(x)$ be an arbitrary function such that vanishes in the boundary layer and

[TABLE]

where $\gamma_{\mathsf{Kn}}$ is defined in (23). Then for such a $c(x)$, one certainly has

[TABLE]

According to Theorem 1, when $\mathsf{Kn}$ is small, $f_{0}$ and $g$ approach the diffusion limit. Taking boundary layers into account, we write

[TABLE]

where $f_{L}$ and $g_{L}$ stand for the layers of the two functions, and $f_{I}$ and $g_{I}$ are the interior solutions. The boundary layers are supported in a thin layer (denoted as $\Omega_{L}$) in the vicinity of the boundary $\partial\Omega$ with $\mathcal{O}(\mathsf{Kn})$ width. The interior $f_{I}$ and $g_{I}$ supported on $\Omega_{I}=\Omega\backslash\Omega_{L}$ are well-approximated by the diffusion limit:

[TABLE]

where $\rho_{f},\rho_{g}$ solve the following equations:

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

Plugging it back into (28) and using the fact $c(x)=0,\forall x\in\Omega_{L}$, we obtain:

[TABLE]

Let $G(x,y)$ be the Green’s function for the operator $C\Delta-\sigma_{a0}$, i.e.,

[TABLE]

Then

[TABLE]

where $\mathrm{d}\mu(y)$ is the surface measure on $\partial\Omega$. Note that one can find $c(x)$ such that

[TABLE]

Indeed, using quadrature rule for (33), one has

[TABLE]

Thus

[TABLE]

then one just need to choose $c(x)$ such that it is perpendicular to the space spanned by $\{\frac{\partial G(x,y_{i})}{\partial n}\frac{\partial G(x,y_{j})}{\partial n},~{}\forall i,j\}$. Plugging (34) into (31), we derive that

[TABLE]

It is seen that to ensure $|\langle\gamma,c\rangle_{L^{2}(\mathrm{d}{x})}|<\delta$, one simply needs

[TABLE]

Considering $\tilde{\sigma}_{a}$ is of $\mathcal{O}(1)$, there exists a constant $C$ such that $\kappa_{a}\geq C\frac{\delta}{\mathsf{Kn}^{2}}$, we finish the proof.

∎

Remark 3**.**

Two immediate take-away information from the theorem:

• •

When $\delta$ is small, meaning that the measurement is relatively accurate, then one gets better recovery of the absorption coefficient as expected.

• •

When $\mathsf{Kn}$ shrinks, the system approaches to the diffusion limit, and the distinguishability coefficient grows dramatically, indicating that the linearized inverse problem is highly ill-conditioned. This phenomena is aligned with the ill-posedness of the Calderón problem.

3.3. Ill-conditioning in the diffusion limit (discrete level)

In this subsection, we revisit the above observation in the discrete setting when solving the inverse problem numerically, in which case the matrix to be inverted becomes highly ill-conditioned as $\mathsf{Kn}$ shrinks. Let us sample $N_{x}$ quadrature points in $\Omega$: $x_{i},i=1,\ldots,N_{x}$, and their corresponding weights are denoted by $w_{i}$. Suppose the measurements are collected at discrete boundary points $\{y_{k}\}_{k=1}^{N_{y}}\subset\partial\Omega$, and there are $\phi_{d},d=1,2,\cdots,N_{\phi}$ different sets of incoming data where $\sup_{d}\|\phi_{d}\|_{L_{\infty}(\Gamma^{-})}\leq 1$. Then the linear system (24) can be rewritten into the form of

[TABLE]

with

[TABLE]

and $f_{0}$ and $g$ solve (19) and (21) respectively with $\phi_{d}$ and $\delta_{y_{k}}$ as boundary condition, here $\delta_{y_{k}}$ is hat function centered at $x=y_{k}$. For simplicity we use index $p$ to denote sub-index $k$ and $d$. Then the linear system (37) can be further written into a compact form:

[TABLE]

where $\mathsf{A}\in\mathbb{R}^{N_{p}\times N_{x}}$ with entries $\mathsf{A}_{pi}=\gamma_{\mathsf{Kn}}(x_{i};y_{k},\phi_{d})w_{i}$ such that $(k,d)=p$. $\mathbf{\tilde{\sigma}_{a}}\in\mathbb{R}^{N_{x}}$ is the discretization of $\tilde{\sigma}_{a}$ on $N_{x}$ quadrature nodes, and data $\mathbf{b}\in\mathbb{R}^{N_{p}}$. Here $N_{p}=N_{y}\times N_{\phi}$ denotes the total number of data points we have, and it is a product of $N_{\phi}$, the number of inflow data and $N_{y}$, the number of measurement positions. We show below that as $\mathsf{Kn}\to 0$ the matrix $\mathsf{A}$ becomes more and more singular, making the inversion impossible.

Theorem 3**.**

Assume that $\mathsf{A}$ is nonsingular. When $\mathsf{Kn}$ is small, the condition number of matrix $\mathsf{A}^{T}\mathsf{A}$ scales as

[TABLE]

Moreover, in 1D, $\mathsf{A}$ is approximately low rank, in the sense that it only has no more than 3 singular values of size $\mathcal{O}(\mathsf{Kn})$, and all the rest are of size $\mathcal{O}(\mathsf{Kn}^{3/2})$.

Proof.

We first the prove the theorem in any dimension and refine the result in 1D. According to the diffusion theory, in the zero limit of $\mathsf{Kn}$, $\mathsf{A}$ can be decomposed into two parts:

[TABLE]

where $\mathsf{A}_{L}\in\mathbb{R}^{N_{p}\times N_{L}},~{}(N_{L}\ll N_{x})$ are the sampled points in the layer, and $\mathsf{A}_{I}\in\mathbb{R}^{N_{p}\times(N_{x}-N_{L})}$ represents the sampled points in the interior. To analyze the rank of the matrix, we rewrite $\mathsf{A}^{T}\mathsf{A}$ as

[TABLE]

When $\mathsf{Kn}$ is small, $\mathsf{A}_{I}=\mathsf{A}^{0}_{I}+\mathcal{O}(\mathsf{Kn})$, where each row in $\mathsf{A}_{I}^{0}$ are:

[TABLE]

Therefore, $\mathsf{A}^{T}\mathsf{A}=\mathsf{Kn}^{2}\mathsf{P}+\mathcal{O}(\mathsf{Kn}^{3})$, where

[TABLE]

Since $\mathsf{A}_{I}^{0T}\mathsf{A}_{I}^{0}$ is symmetric, it is diagonizable. Denote:

[TABLE]

with $\mathsf{Q}\in\mathbb{R}^{(N_{x}-N_{L})\times(N_{x}-N_{L})}$ the collection of eigenvectors and $\mathsf{D}$ the diagonal matrix of eigenvalues. Then we multiply both sides by matrix $\mathsf{X}$ and $\mathsf{X}^{-1}$ defined as below

[TABLE]

we derive that

[TABLE]

Note that the number of elements in $A_{L}$ is of order $\mathcal{O}(\mathsf{Kn})$ due to the fact that the layer length is order of $\mathcal{O}(\mathsf{Kn})$, and its elements are order $\mathcal{O}(1)$ thanks to the maximal principle: the integrand function $|f_{L}g_{L}|\leq\|\phi\|_{L^{\infty}(\Gamma_{-})}\|\delta_{y_{k}}\|_{L^{\infty}(\mathrm{d}{x})}=\mathcal{O}(1)$. Therefore, we could rewrite the equation above:

[TABLE]

By Gershgorin circle theorem, all eigenvalues lie in Gershgorin disc, meaning that $|\lambda(\mathsf{A}^{T}\mathsf{A})-\lambda(\mathsf{Kn}^{2}\mathsf{D}_{1})|<\mathcal{O}(\mathsf{Kn}^{3})$. Since the eigenvalues in $\mathsf{D}_{1}$ can be $\mathcal{O}(1)$ or 0, the largest eigenvalue in $\mathsf{A}^{T}\mathsf{A}$ is $\mathcal{O}(\mathsf{Kn}^{2})$ and the smallest is $\mathcal{O}(\mathsf{Kn}^{3})$, the condition number of $\mathsf{A}^{T}\mathsf{A}$ is larger than $\mathcal{O}\left(\frac{1}{\mathsf{Kn}}\right)$.

Moreover, in 1D, we show in the appendix that $\mathsf{D}$ is indeed a low rank matrix itself with rank not exceeding 3. Therefore, there are at most 3 nonzero eigenvalues in $\mathsf{D}$ and all the rest are zeros. Putting this information back to $\mathsf{A}$, the result directly follows.

∎

Remark 4**.**

We emphasize the difference between Theorem 2 and Theorem 3. Theorem 2 is the study of $\kappa_{a}$, which represents $\|\mathcal{A}^{-1}\|$: it tells that suppose the measurement is different, how different could $\sigma_{a}$ be. However, numerically it is the condition number of $\mathsf{A}$ that is playing the role. If $\mathsf{A}$ is low rank, for example, rank $k$ out of an $n$ dimensional space, then there are infinite many $\sigma_{a}$ that could lead to the same measurement and the space they span is $n-k$ dimensional.

Similar to the analysis on the continuous level, we immediately conclude that the ill-conditioned matrix $\mathsf{A}$ in $\mathsf{Kn}\to 0$ limit leads to the fact that $\sigma_{a}$ is hard to be recovered accurately, which is consistent with the ill-posedness of the Calderón problem. More precisely, we have the following estimate theorem.

Theorem 4**.**

Define the distinguishability coefficient in the discrete setting as

[TABLE]

where $\|\cdot\|$ denote vector $l^{2}$-norm and $\mathsf{A}\tilde{\sigma}_{a}=b$, and $\Gamma_{\delta}=\{\sigma_{a}:\|\mathsf{A}\sigma_{a}-b\|\leq\delta\}$. Assume that $\mathsf{A}$ is nonsingular, then there exists a constant $C$ such that

[TABLE]

Proof.

We again decompose $\mathsf{A}$ as in (40). Note that for each row of $\mathsf{A}_{I}$, we have

[TABLE]

where

[TABLE]

Here $G(x,y)$ is the Green’s function defined in (32) and $\Delta y=\max_{i}\Delta y_{i}$. Denote $\mathbf{c}=\mathbf{\sigma}_{a}-\mathbf{\tilde{\sigma}}_{a}$, then

[TABLE]

Since

[TABLE]

and from (3.3), one can always choose $\mathbf{c}_{I}$ as long as $N_{\phi}\times N_{y}<N_{x}$ (just pick $c_{I}$ from (35) ) such that

[TABLE]

we have

[TABLE]

Therefore, the requirement (44) implies (43). ∎

4. Recover Scattering Coefficient $\sigma_{s}$

In this section we discuss how to recover the scattering coefficient given $\sigma_{a}$. In subsection 4.1 we set up the inverse problem, and the following two subsections are devoted to the non-injectivity in 1D and ill-conditioning in multi-dimension in the zero limit of the Knudsen number.

4.1. Inverse problem set-up

We recall the equation again:

[TABLE]

Here $\sigma_{a}$ as known, and we make a guess for $\sigma_{s}$, and linearize the equation around $\sigma_{s0}$, assuming the deviation $|\tilde{\sigma}_{s}|:=|\sigma_{s}-\sigma_{s0}|$ is much smaller than $\sigma_{s}$. The background solution $f_{0}$ solves

[TABLE]

As done in the last section, we drop the higher order terms, and the fluctuation $\tilde{f}:=f-f_{0}$ satisfies the following equation:

[TABLE]

To recover $\tilde{\sigma}_{s}$, we look for the linear mapping from the incoming information $\phi$ to some computable quantity $b(\delta_{y},\phi)$ (to be determined below), i.e.,

[TABLE]

To this end, we consider an auxiliary function $g(x,v)$ that satisfies the adjoint problem:

[TABLE]

Multiply Equation (47) with $g$ and (48) with $\tilde{f}$ and compare these two equations, with the Green’s identity, one gets:

[TABLE]

Then we define:

[TABLE]

then (49) becomes:

[TABLE]

where

[TABLE]

which is again the difference between the measured data and computed data, given by the boundary condition $\phi$, evaluated at $y$. The inverse problem then is equivalent to the Fredholm first type problem: for all $y$ and $\phi$, we prepare $\gamma(x;\delta_{y},\phi)$ and $b(\delta_{y},\phi)$, and use them to invert for $\sigma_{s}$. For the ease of notation we write $b(y,\phi)$ as $b(\delta_{y},\phi)$ with $\delta_{y}$ representing the boundary condition for $g$.

4.2. Non-injectivity in 1D

Similar to the case of recovering $\sigma_{a}$, recovering $\sigma_{s}$ becomes harder as $\mathsf{Kn}$ shrinks to zero. In 1D, formulae can be made explicitly.

We first restrict our attention to the critical case by setting $\sigma_{a}\equiv 0$. When there is no absorption, the only interaction between particles is scattering, and thus mass is preserved. We show below that in this case, the problem is non-injective in the sense that $\gamma$ cannot provide enough variations to distinguish $\sigma_{s}$ at different $x$, and that different $\sigma_{s}$ could lead to the same measurement provided the same data.

Proposition 1**.**

Given arbitrary $\phi_{d}$ and $\delta_{y}$, let $f_{0}$ solve (46) and $g$ solve (48). Let $\gamma_{\mathsf{Kn}}$ be defined in (50). Then in 1D, if $\sigma_{a}\equiv 0$, $\gamma_{\mathsf{Kn}}$ is a constant independent of $x$.

Proof.

Here we drop the $\mathsf{Kn}$ dependence in the proof as it will change the result. Denote $\langle f\rangle=\frac{1}{2}\int_{-1}^{1}f(x,v)\mathrm{d}{v}$, then

[TABLE]

Notice that

[TABLE]

thus

[TABLE]

which readily implies that $\frac{\mathrm{d}}{\mathrm{d}x}\gamma_{\mathsf{Kn}}=0$ .

∎

The non-injectivity is immediate:

Theorem 5**.**

In 1D critical case, the inverse problem for (51) is non-injective.

Proof.

According to proposition 1, $\gamma_{\mathsf{Kn}}$ is a constant, and thus one can only get $\int\tilde{\sigma}_{s}\mathrm{d}{x}$. Therefore, the variation in $\tilde{\sigma}_{s}$ is not recoverable, and the problem is non-injective. ∎

In the subcritical case with $\sigma_{a}>0$, the recovery of the scattering coefficient is very similar.

Proposition 2**.**

For arbitrary inflow data $\phi$ and Dirac delta function $\delta_{y}$, as the Knudsen number $\mathsf{Kn}\rightarrow 0$, $\frac{\mathrm{d}}{\mathrm{d}x}\gamma_{\mathsf{Kn}}\rightarrow 0$.

Proof.

Notice that

[TABLE]

we can derive that

[TABLE]

this is equivalent to

[TABLE]

therefore we have

[TABLE]

∎

This proposition indicates that as $\mathsf{Kn}\to 0$, whatever boundary condition we provide for $f_{0}$ and $g$, $\gamma_{\mathsf{Kn}}$ is going to be flat, and thus not able to distinguish the variation in $\sigma_{s}$, making the problem more and more non-injective.

Theorem 6**.**

When Kundsen number $\mathsf{Kn}\rightarrow 0$, solving linear system $\langle\gamma_{\mathsf{Kn}},\tilde{\sigma}_{s}\rangle_{L^{2}(\mathrm{d}{x})}=b$ is non-injective. More specifically, in the zero limit of $\mathsf{Kn}$, the space spaned by the kernel $\gamma_{\mathsf{Kn}}$ has finite rank, impossible to reflect full information of $\tilde{\sigma}_{s}$.

Proof.

The non-injectivity of the linear system directly follows from Proposition 2. Moreover, we show in below that the space that $\gamma$ resides in is of low rank in the zero limit of $\mathsf{Kn}$. Let the domain of $x$ be $x_{l}\leq x\leq x_{r}$ and denote

[TABLE]

As Knudsen number $\mathsf{Kn}$ goes to zero, $\rho_{f}$ and $\rho_{g}$ solves the diffusion equation in the leading order:

[TABLE]

with velocity averaged boundary data $\xi_{f}$ and $\xi_{g}$ at only two points: left and right end points. That is,

[TABLE]

In 1D, there are two Green’s functions:

[TABLE]

and

[TABLE]

Then $\rho_{f}$ and $\rho_{g}$ can be written as

[TABLE]

As a result, $\rho_{f}(x)\rho_{g}(x)\in\text{span}\{G_{1}(x)^{2},G_{1}(x)G_{2}(x),G_{2}(x)^{2}\}$, which means that the function space of $\rho_{f}\rho_{g}$ is of low rank (rank 3). Combining with Proposition 2, we see that:

[TABLE]

meaning $\frac{\mathrm{d}}{\mathrm{d}x}\gamma_{\mathsf{Kn}}$ is low rank as well. Therefore it is impossible to recover $\tilde{\sigma}_{s}$ from the linear mapping $\langle\gamma_{\mathsf{Kn}},\cdot\rangle_{L^{2}(\mathrm{d}{x})}$. ∎

Remark 5**.**

Theorem 6 coincides with our intuition in Section 1 since when Knudsen number is small, the scattering will dominate absorption as neutron travels through the medium, and the case described here converges to the critical case.

4.3. Ill-conditioning in higher dimensions

In higher dimension, we show that the inverse problem become more and more singular as the Knudsen number approaches zero.

Recall the linear mapping in this scenario:

[TABLE]

where $f_{0}$ and $g$ solve (46) (48), respectively. We then investigate the sensitivity of reconstructing $\tilde{\sigma}_{s}$. The main theorem states as follows.

Theorem 7**.**

Assume that the map from $\tilde{\sigma}_{s}$ to $b$ is injective. Given an error tolerance $\delta$ on the measurement, the distinguishability coefficient in reconstructing $\tilde{\sigma}_{s}$ grows as $\delta$ grows and $\mathsf{Kn}$ shrinks. Namely, there exists a constant $C$ such that

[TABLE]

where

[TABLE]

Proof.

When $\mathsf{Kn}$ is small, we decompose $f_{0}$ and $g$ into a layer part that accounts for the boundary layer supported in the vicinity of the boundary $\partial\Omega$ with $\mathcal{O}(\mathsf{Kn})$ width and an interior part:

[TABLE]

Decompose also $\gamma_{\mathsf{Kn}}$ as $\gamma_{\mathsf{Kn}}=\left(\gamma_{\mathsf{Kn}}\right)_{L}+\left(\gamma_{\mathsf{Kn}}\right)_{I}$, then we examine the interior part $(\gamma_{\mathsf{Kn}})_{I}$ by means of $f_{I}$ and $g_{I}$. Upon asymptotic expansion, the interior solutions $f_{I}$ and $g_{I}$ are approximated as follows

[TABLE]

hence,

[TABLE]

where $C$ is a constant depending on the dimension of the problem. Now consider $c(x)$ such that it vanishes in the boundary layer and that

[TABLE]

then one certainly has

[TABLE]

Since $\langle\gamma_{\mathsf{Kn}},c\rangle_{L^{2}(\mathrm{d}{x})}=\langle(\gamma_{\mathsf{Kn}})_{I},c_{I}\rangle_{L^{2}(\mathrm{d}{x})}=\langle-\frac{\mathsf{Kn}}{\sigma_{s0}^{2}}C\nabla_{x}\rho_{f}\cdot\nabla_{x}\rho_{g},c_{I}\rangle_{L^{2}(\mathrm{d}{x})}+\mathcal{O}(\mathsf{Kn}^{2})$, we immediately have

[TABLE]

which implies that

[TABLE]

∎

Remark 6**.**

Notice that in the expression (53), the leading term in $\gamma$ has the structure of $\nabla_{x}\rho_{f}\cdot\nabla_{x}\rho_{g}$, where $\rho_{f}$ and $\rho_{g}$ solve the diffusion equation with corresponding boundary condition, which has a Green’s function formation. Therefore, one can show that it is asymptotically “low rank” in the sense that one can find a nonzero $c(x)$ such that

[TABLE]

Details follow the argument in the proof of Theorem 2 .

5. Numerical test

In this section, we conduct a few numerical experiments in 1D to check the conditioning of the inverse problem, and show that it degrades when $\mathsf{Kn}$ goes to zero, as indicated by the theory above. More specifically, we check the variation of $\gamma(x;\delta_{y},\phi)$ in $x$ and examine the singular values of the matrix $\mathsf{A}$, whose element takes values $\mathsf{A}_{pi}=\gamma(x_{i};\delta_{y_{k}},\phi_{d})$ with $(k,d)=p$.

Recall the definition of $\gamma(x;\delta_{y_{k}},\phi_{d})$ in (23) and (50), one needs to solve the forward problems for $f_{0}$ and $g$ with boundary conditions $\phi_{d}$, and $\delta_{y}$, respectively. The forward solver we adopt is the fully implicit solver [24] with Generalized minimal residual method (GMRES) under a tolerance of $10^{-10}$. For all the examples below, the spatial domain is chosen as $\Omega=[0,1]$ and discretized with $N_{x}=200$ uniformly distributed nodes. We also sample $N_{v}=80$ grid points in the velocity domain $\mathcal{S}=[-1,1]$. Note that for the $\mathsf{Kn}$ we considered here, the spatial mesh is fine enough to resolve the boundary layers.

Now we need to decide the boundary condition $f_{0}|_{\Gamma_{-}}=\phi_{d}(x,v)$ and $g|_{\Gamma_{-}}=\delta_{y}(x)$ such that we extract the most information from $f_{0}$ and $g$. Since $f_{0}$ solves a linear equation, we can set delta functions as its inflow data, i.e.,

[TABLE]

The boundary condition for $g$ is easier to set up. Note that in 1D, there are only two boundary points, $x=0$ and $x=1$, therefore, the boundary for $g$ reads

[TABLE]

Then the associated matrix is of size $\mathsf{A}\in\mathbb{R}^{160\times 200}$.

5.1. Recover Absorption Coefficient

The first test addresses the problem of recovering $\sigma_{a}$. Here the background absorption is set to be $\sigma_{s}=1+\frac{1}{1.5+\sin(2\pi x)}$ and the scattering coefficient takes the form $\sigma_{a0}=4+\frac{1}{2}\sin(4\pi x)$. Here

[TABLE]

As predicted in Theorem 3, in 1D, the matrix $\mathsf{A}$ is approximately low rank with rank 3, which indicates that its singular value decays to zero quickly. This is demonstrated in Fig. 1, wherein we plot the singular values for a variety of $\mathsf{Kn}$. It is easy to see that only three singular values stand out, coincide with the rank 3 argument. And the difference between the first three singular values with the rest are more pronounced with smaller $\mathsf{Kn}$.

Moreover, we see in the proof of Theorem 2 (see equation (31)) that $\gamma_{\mathsf{Kn}}\approx\mathsf{Kn}\rho_{f}\rho_{g}$ for small $\mathsf{Kn}$, which lives in a space spanned by $G^{2}_{1}$, $G^{2}_{2}$ and $G_{1}G_{2}$ with $G_{1,2}$ standing for the Green’s functions. The plots in 2 show the first three eigenvectors of $\mathsf{A}$ as $\mathsf{Kn}\to 0$.

5.2. Recover Scattering Coefficient

In the second test, we aim to recover $\sigma_{s}$. Here

[TABLE]

We first compute $\gamma_{\mathsf{Kn}}$ (stored in $\mathsf{A}$) with $\sigma_{s0}=1+\frac{1}{1.5+\sin(2\pi x)}$ and $\sigma_{a}=0$. It is seen from Proposition 1 that $\gamma_{\mathsf{Kn}}$ stays unchanged in $x$ direction, i.e., $\frac{\mathrm{d}\gamma_{\mathsf{Kn}}}{\mathrm{d}x}=0$. To show this, we plot the matrix $\mathsf{A}$ in Fig. 3 and see that the entire matrix is roughly flat in $x$ direction.

Next we consider the case when $\sigma_{a}$ is nontrivial: $\sigma_{a}=2^{-4}(4+\frac{1}{2}\sin(4\pi x))$, and $\sigma_{s0}=\frac{1}{2^{-4}}(1+\frac{1}{1.5+\sin(2\pi x)})$. As predicted in Proposition 2, $\frac{\mathrm{d}\gamma_{\mathsf{Kn}}}{\mathrm{d}x}\sim\mathsf{Kn}\frac{\mathrm{d}\rho_{f}\rho_{g}}{\mathrm{d}x}\to 0$ as $\mathsf{Kn}\to 0$, and this is demonstrated in Fig. 4, where we plot $\frac{\mathrm{d}\gamma_{\mathsf{Kn}}/\mathrm{d}x}{\mathrm{d}\rho_{f}\rho_{g}/\mathrm{d}x}$. As $\mathsf{Kn}$ decreases by $1/2$, the ratio decrease by $1/2$ as well.

In the end we test the singular value decay of $\mathsf{A}_{I}$, the interior part of $\mathsf{A}$. As $\mathsf{Kn}\rightarrow 0$, we expect that its element $\gamma_{\mathsf{Kn}}\approx\mathsf{Kn}\nabla_{x}\rho_{f}\cdot\nabla_{x}\rho_{g}$ (see equation (53)), and thus $\mathsf{A}_{I}$ is a low rank matrix. We plot the singular values Fig. 5, and as before, only the first three singular values dominates, and the rest vanishes at the order of $\mathsf{Kn}$.

We also show that for small $\mathsf{Kn}$, $\gamma_{\mathsf{Kn}}\approx\mathsf{Kn}\nabla_{x}\rho_{f}\cdot\nabla_{x}\rho_{g}$, which lives in a space spanned by $(\partial_{x}G_{1})^{2}$, $(\partial_{x}G_{2})^{2}$ and $\partial_{x}G_{1}\partial_{x}G_{2}$ with $G_{1,2}$ standing for the Green’s functions. The plots in 6 show the first three eigenvectors of $\mathsf{A}$ as $\mathsf{Kn}\to 0$.

6. Appendix

Appendix I: Rank of matrix $\mathsf{A}_{I}^{0}$ Here we examine the rank of matrix $\mathsf{A}_{I}^{0}$. We show below that it has rank less than 4 in 1D, but not necessarily low rank in higher dimensions. Denote

[TABLE]

as the solutions to the discrete form of equations (29) (30):

[TABLE]

where $\mathsf{B}_{Nx\times Nx}$ is the discrete version of operator $C\nabla_{x}\cdot\nabla_{x}-\sigma_{a0}$. Write $\mathsf{B}^{-1}=(b_{ij})_{N_{x}\times N_{x}}$, and recall that the entries of $\mathsf{A}_{I}^{0}$ is (41), we have

[TABLE]

Here without abuse of notation, we still the denote number of interior points $N_{x}$.

Then in 1D, since we only have two boundary points, $N_{y}=2$, and $\xi_{j}^{(l)}=0$ for all $j\neq 1$ or $N_{x}$, and $l=1,2,\cdots,N_{\phi}$. Therefore, $\mathsf{A}_{I}^{0}$ can be simplified to

[TABLE]

Then it is easy to see that each component in the summation is rank one, therefore in total the rank of $\mathsf{A}_{I}^{0}$ does not exceed $3$.

In higher dimensions, on the contrary, does not guarantee a low rank property of $\mathsf{A}_{I}^{0}$. To see this, let $N_{J}=\{j_{0},j_{1},\cdots,j_{J}\}\subset\{1,2,\cdots,N_{x}\}$ be the collection of indices such that $\xi_{j}^{l}\neq 0$, $j\in N_{J}$, $\forall l$, then $J<N_{x}$. Without loss of generality, we let $\mathbf{\delta_{y}^{(1)}}=\cdots=\mathbf{\delta_{y}^{(N_{\phi})}}=\mathbf{\delta_{y_{j_{1}}}}$, $\mathbf{\delta_{y}^{(N_{\phi}+1)}}=\cdots=\mathbf{\delta_{y_{j_{2}}}}$, $\cdots$, $\mathbf{\delta_{y}}^{(Ny-1)N_{\phi}+1}=\cdots=\mathbf{\delta_{y_{j_{N_{y}}}}}$, then

[TABLE]

where

[TABLE]

and

[TABLE]

Therefore, although $\text{rank}(\mathsf{AA})\leq N_{J}$, adding them together may still give a full rank matrix $\mathsf{A}_{I}^{0}$.