On diffusive scaling in acousto-optic imaging

Francis J. Chung Department of Mathematics, University of Kentucky, Lexington, KY 40506, USA fj.chung@uky.edu , Ru-Yu Lai School of Mathematics, University of Minnesota, Minneapolis, MN 55455, USA rylai@umn.edu and Qin Li Department of Mathematics and Wisconsin Institute for Discoveries, University of Wisconsin-Madison, Madison, WI 53705, USA qinli@math.wisc.edu

Abstract.

Acousto-optic imaging (AOI) is a hybrid imaging process. By perturbing the to-be-reconstructed tissues with acoustic waves, one introduces the interaction between the acoustic and optical waves, leading to a more stable reconstruction of the optical properties. The mathematical model was described in [27], with the radiative transfer equation serving as the forward model for the optical transport. In this paper we investigate the stability of the reconstruction. In particular, we are interested in how the stability depends on the Knudsen number, $\mathsf{Kn}$ , a quantity that measures the intensity of the scattering effect of photon particles in a media. Our analysis shows that as $\mathsf{Kn}$ decreases to zero, photons scatter more frequently, and since information is lost, the reconstruction becomes harder. To counter this effect, devices need to be constructed so that laser beam is highly concentrated. We will give a quantitative error bound, and explicitly show that such concentration has an exponential dependence on $\mathsf{Kn}$ . Numerical evidence will be provided to verify the proof.

Chung is supported in part by Simons Collaboration Grant 582020. Lai is supported in part by the NSF grant DMS-1714490. Li is supported in part by DMS-1619778 and DMS-1750488. The collaboration started during the workshop “Mathematics in Optical Imaging” in IMA, 2019, and was continuously supported by KI-net, DMS-1107291 in the workshop “Forward and Inverse Problems in Kinetic Theory” in UW-Madison, 2019. The authors thank Professor John Schotland for helpful discussions.

Key words and phrases: Acousto-optic imaging, inverse problem, radiative transfer equation.

1. Introduction

High energy light is the classical way to probe optical properties of thick and highly scattering media. In contrast to most destructive experiments, here, high energy photons are injected into biological tissues, and the outgoing light intensities are measured at the surface of the samples by detectors. The map that maps the incoming light intensity to the outgoing light intensity carries optical information of the biological tissues, and is utilized to reconstruct optical parameters.

However, a long standing challenge of the reconstruction process centers around the stability issue. If the injected photons carry low energy, the resulting images are very blurred, and the reconstruction is typically unstable [10, 11, 19, 20, 31, 34]. This phenomenon is the so-called diffuse optical tomography, especially when the injected light is near-infrared. In this situation, the forward model is the classical elliptic type equation, and mathematically, recovering the optical parameter amounts to reconstructing the diffusion coefficient in the elliptic equation using the Dirichlet-to-Neumann map [42], and is proven mathematically to be logarithmically unstable [1]. Multiple strategies are adopted to “stabilize” the problem [36, 40], both by adjusting the experimental modalities upfront, or introducing image deblurring techniques as a post-processing [18, 29, 39, 41].

We follow the first approach in this paper. In particular, we investigate a modality called acousto-optic imaging (AOI), and study the stability of the media reconstruction when biological tissue is optically thick. AOI is one of the state-of-the-art hybrid imaging processes that fuse two or more imaging modalities to form a new imaging technique [8, 30]. In particular, AOI is based on the acousto-optic effect: optical properties of the medium are modified upon interaction with acoustic radiation. Presumably one can obtain more information upon a sequence of such modifications. The procedure combines the contrast advantage inheriting from optical properties and the resolution advantage of ultrasound [15, 32, 44], and is expected to provide a more stable reconstruction.

In general, there are two types of AOI - direct imaging and tomographic imaging. In direct imaging, the sample is simultaneously affected by an ultrasound beam and illuminated with a laser source. The tagged photon, resulting from the interaction between the light and the ultrasound, is detected by external ultrasound transducers and then forms the image [28]. Tomographic imaging, the second type of AOI, makes use of a reconstruction algorithm to track down the optical properties of the medium through various inversion procedures. In the latter case, the inverse problems for incoherent AOI [2, 3, 4, 14, 15, 16] and for coherent AOI with diffuse light [25] were studied and the absorption and the diffusion coefficients of a scattering medium are stably reconstructed. On the other hand, for the problem within the framework of radiative transport theory, relevant results are also investigated in [9, 26, 27]; also see reviews [6, 7, 8, 35].

As a typical hybrid imaging problem, the inverse problem in acousto-optic imaging involves two steps: the first step amounts to reconstructing internal information from triggering the system using ultrasound. Ultrasound waves of different frequencies imposed on a physical system induces a wave-like perturbation to the media, leading to wave-like perturbation in the measurement. From measurements of this perturbation over a big spectrum of frequencies, a Fourier-type calculation provides internal information inside the physical domain. The second step then amounts to reconstructing optical parameters in the modeling equation using the internal data. The first step basically is an inverse Fourier series, and thus is regarded as a stable process, and since the internal data is obtained, it is largely believed that the reconstruction is also stable in the second step.

This is indeed the situation that we will be studying. As a forward model we use the radiative transfer equation (RTE) to characterize the dynamics of photon particles. This equation, besides containing optical parameters, also depends on the Knudsen number $\mathsf{Kn}$ : it is the ratio of mean free path and the typical domain length. In the thick and dense optical environment, the photon particles scatter frequently, and mathematically this amounts to setting this ratio to be small. One can also view this quantity to reflect the energy of the photon particles being used: indeed, for the same media, in the low energy regime, photon particles are scattered more often. One interesting phenomenon associated with low energy, high scattering situation is that, in this regime, mathematically one can asymptotically connects the RTE to a parabolic heat equation. It is a well-accepted fact that reconstruction of parameters are unstable for diffusion type equations, and thus it is a natural expectation that the inverse problem of the RTE is ill-posed accordingly in this regime.

The rest of this paper is organized as follows. In Section 2, we discuss the first step in AOI of determining the internal data, introduce preliminary results for the well-posedness of the transport equation, and present the main result in Theorem 2.4. Section 3 is concerned with the second step in AOI. We derive the reconstruction formula for the absorption coefficient from the reconstructed internal data and focus on the dependence of $\mathsf{Kn}$ in the associated convergence rate. The numerical evidence is presented in Section 4.

2. Problem setup and main results

In this section, we start by introducing some preliminary results and then describe the setup of the problem. Lastly, we present the main theorem.

To begin, consider a smooth bounded domain $\Omega\subset\mathbb{R}^{n}$ . In this paper we use the RTE as the forward model for light propagation on $\Omega$ . We also assume the time scale is significantly larger, and the equation achieves the steady state:

(2.3)

\displaystyle\left\{\begin{array}[]{ll}\theta\cdot\nabla u=\frac{\sigma_{s}(x)}{\mathsf{Kn}}\left(\langle u\rangle-u\right)-\mathsf{Kn}\sigma_{a}(x)u&\hbox{in }\Omega\times\mathbb{S}^{n-1}\,,\\ u|_{\Gamma_{-}}=f\,.\\ \end{array}\right.

The solution $u(x,\theta)$ characterizes the density of photon particles at a physical point $x\in\Omega$ moving in direction $\theta\in\mathbb{S}^{n-1}$ for $n\geq 2$ . The left-hand side of (2.3) is a transport term describing particles at the position $x$ moving in the $\theta$ direction. The terms on the right represent the interaction of the photon particles with the medium, where $\sigma_{a}$ is the absorption coefficient, $\sigma_{s}$ is the scattering coefficient. The scattering operator is $\langle u\rangle-u$ with

\langle u\rangle:=\frac{1}{|\mathbb{S}^{n-1}|}\int_{\mathbb{S}^{n-1}}u(x,\theta)~{}d\theta=:\fint u(x,\theta)~{}d\theta

standing for the angular average of $u$ . Physically it means the particles with different velocities at $x$ are redistributed in an isotropic way at rate $\sigma_{s}(x)$ . It is also common to insert an anisotropic term $k(x,\theta,\theta^{\prime})$ in the integral defined above, but it does not bring extra mathematical difficulty and thus will be neglected here. The parameter $\mathsf{Kn}$ in the equation is called the Knudsen number, and as described above, represents the ratio of mean free path to the typical domain length.

For conciseness of notation, we denote the total absorption coefficient

\sigma=\sigma_{s}+\mathsf{Kn}^{2}\sigma_{a}\,,

and rewrite the right-hand side of (2.3) as a linear operator

(2.4)

\mathcal{L}u:=\frac{\sigma_{s}}{\mathsf{Kn}}\langle u\rangle-{\frac{\sigma}{\mathsf{Kn}}}u\,.

The boundary condition of (2.3) is imposed on $\Gamma_{-}$ : a a set of coordinates at the physical boundary $\partial\Omega$ with inward velocities. More specifically:

(2.5)

\displaystyle\Gamma_{\pm}:=\{(x,\theta)\in\partial\Omega\times\mathbb{R}^{3}:\pm\ n(x)\cdot\theta>0\}\,,

where $n(x)$ is the unit outer normal to $\partial\Omega$ at the point $x\in\partial\Omega$ . Similarly, $\Gamma_{+}$ collects all outgoing coordinates at the boundary.

The well-posedness problem of (2.3) has been extensively studied. We refer the readers to, for instance, [22, 23, 24] for detailed discussion on the existence of solutions to (2.3) under different constraints. For our purpose, we will be using the following boundedness result, with $\mathsf{Kn}$ included and adjusted to fit our setting.

Theorem 2.1 (Theorem 2.1 in [9]).

Suppose $\sigma_{s}$ and $\sigma_{a}$ are bounded and uniformly positive functions in $\Omega$ and the source term $S\in L^{\infty}(\Omega\times\mathbb{S}^{n-1})$ . For any $\mathsf{Kn}>0$ and any boundary value $f\in L^{\infty}(\Gamma_{-})$ , there exists a strong unique solution $u\in L^{\infty}(\Omega\times\mathbb{S}^{n-1})$ to the boundary value problem

(2.8)

\displaystyle\left\{\begin{array}[]{ll}\theta\cdot\nabla u={\frac{\sigma_{s}}{\mathsf{Kn}}}(\langle u\rangle-u)-\mathsf{Kn}\sigma_{a}u+S&\hbox{in }\Omega\times\mathbb{S}^{n-1}\,,\\ u=f&\hbox{on }\Gamma_{-}\,.\\ \end{array}\right.

Furthermore we have the estimate:

(2.9)

\|u\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\leq\frac{C}{\mathsf{Kn}}\left(\|S\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}+\frac{1}{\mathsf{Kn}}\|f\|_{L^{\infty}(\Gamma_{-})}\right),

holds for some positive constant $C$ depending only on $\sigma_{a}$ and $\sigma_{s}$ .

In the low energy regime, $\mathsf{Kn}$ is small, and one can perform the asymptotic expansion in $\mathsf{Kn}$ . In the leading order the scattering term $\frac{\sigma_{s}}{\mathsf{Kn}}(\langle u\rangle-u)$ becomes dominant, and it drives $u\sim\langle u\rangle$ , meaning $u$ is a constant in $\theta$ at every location $x$ , see [19]. This degeneracy of information, as will be explained in details in Section 3, is one of the main reason for the rise of instability.

With the well-posedness result in Theorem 2.1, the trace can be taken and we define the albedo operator $\mathcal{A}$ :

(2.10)

\displaystyle\mathcal{A}:L^{\infty}(\Gamma_{-})\rightarrow L^{\infty}(\Gamma_{+}):\quad u|_{\Gamma_{+}}=\mathcal{A}f\,.

It maps the incoming light intensity to the outgoing light intensity.

2.1. Problem setup

Acousto-optic imaging is an imaging modality that introduces an acoustic wave to the medium, and thus induces a wave-like perturbation to the coefficients $\sigma_{a}$ and $\sigma_{s}$ , that gives rise to new coefficients $\sigma_{a,\varepsilon}$ and $\sigma_{s,\varepsilon}$ . Following the setup presented in [27], the modulated coefficients take the form

(2.11)

\displaystyle\sigma_{a,\varepsilon}:=(1+\varepsilon\cos(q\cdot x))\sigma_{a}\,,\quad\sigma_{s,\varepsilon}:=(1+\varepsilon\cos(q\cdot x))\sigma_{s}\,,

where $0<\varepsilon\ll 1$ is the strength of introduced acoustic effect. The $\cos{(q\cdot x)}$ perturbation comes from the assumption that the introduced acoustic wave is basically sinusoidal, and the wave vector $q$ can be experimentally adjusted.

We denote $u_{\varepsilon}$ to be the solution to the RTE with the modulated coefficients, namely:

(2.12)

\displaystyle\begin{cases}\theta\cdot\nabla u_{\varepsilon}=\frac{\sigma_{s,\varepsilon}}{\mathsf{Kn}}(\langle u_{\varepsilon}\rangle-u_{\varepsilon})-\mathsf{Kn}\sigma_{a,\varepsilon}u_{\varepsilon}\,,\\ u_{\varepsilon}|_{\Gamma_{-}}=f\,.\\ \end{cases}

For small $\varepsilon$ , the well-posedness result for (2.12) still holds true according to Theorem 2.1. Hence, for each $q$ , we can define the corresponding albedo operator:

\mathcal{A}_{\varepsilon,q}:f\mapsto u_{\varepsilon}|_{\Gamma_{+}}\,,

where $u_{\varepsilon}$ is the solution to (2.12). Likewise, we also define the perturbed scattering operator $\mathcal{L}_{\varepsilon}$ :

\mathcal{L}_{\varepsilon}u:=\frac{\sigma_{s,\varepsilon}}{\mathsf{Kn}}\langle u\rangle-\frac{\sigma_{\varepsilon}}{\mathsf{Kn}}u\,,

where the modulated total absorption coefficient is

\sigma_{\varepsilon}:=\sigma_{s,\varepsilon}+\mathsf{Kn}^{2}\sigma_{a,\varepsilon}\,.

To extract internal information, we first introduce the adjoint RTE:

(2.13)

\displaystyle\begin{cases}-\theta\cdot\nabla v=\mathcal{L}v\,,\\ v|_{\Gamma_{+}}=g\,.\\ \end{cases}

It is straightforward to show that if $v$ solves the adjoint RTE, then

\tilde{v}(x,\theta)=v(x,-\theta)

solves the regular RTE. Therefore the well-posedness result applies and we can also define the adjoint albedo operator

\tilde{\mathcal{A}}:L^{\infty}(\Gamma_{+})\rightarrow L^{\infty}(\Gamma_{-}):v|_{\Gamma_{-}}=\tilde{\mathcal{A}}g\,.

Furthermore:

\widetilde{\mathcal{A}(\tilde{g})}=\tilde{\mathcal{A}}(g)\,,

so $\tilde{\mathcal{A}}$ is determined by $\mathcal{A}$ . Here for a given function $\varphi$ , the expression $\tilde{\varphi}$ indicates the reflection of $\varphi$ in the $\theta$ variable, namely, $\tilde{\varphi}(x,\theta):=\varphi(x,-\theta)$ .

Now we can use the adjoint RTE to obtain an internal functional. Using (2.12) and (2.13) we find

	$\displaystyle\int_{\partial\Omega\times\mathbb{S}^{n-1}}\theta\cdot n(x)(u_{\varepsilon}v)~{}dxd\theta$	$\displaystyle=\int_{\Omega}\int_{\mathbb{S}^{n-1}}(v\mathcal{L}_{\varepsilon}u_{\varepsilon}-u_{\varepsilon}\mathcal{L}v)~{}dxd\theta\,$
		$\displaystyle=\int_{\Omega}\int_{\mathbb{S}^{n-1}}v(\mathcal{L}_{\varepsilon}-\mathcal{L})u_{\varepsilon}~{}dxd\theta\,,$
(2.14)			$\displaystyle=\int_{\Omega}\int_{\mathbb{S}^{n-1}}\varepsilon\cos(q\cdot x)(\mathcal{L}u_{\varepsilon})v~{}dxd\theta\,,$

where we used the self-adjoint property of $\mathcal{L}$ and $\mathcal{L}_{\varepsilon}$ as well as $\mathcal{L}_{\varepsilon}-\mathcal{L}=\varepsilon\cos(q\cdot x)\mathcal{L}$ . We can label the left-hand side of this $BT$ (for boundary term) and rewrite it in terms of boundary values of $v$ and $u_{\varepsilon}$ :

BT:=\int_{\Gamma_{-}}\theta\cdot nf\tilde{\mathcal{A}}(g)~{}dS+\int_{\Gamma_{+}}\theta\cdot n\mathcal{A}_{\varepsilon,q}(f)g~{}dS\,.

Since the boundary values are known, this implies that $BT$ is known.

Finally, in (2.1), we can write $\mathcal{L}u_{\varepsilon}=\mathcal{L}u+\mathcal{L}(u_{\varepsilon}-u)$ , so if $u_{\varepsilon}-u$ is small when $\varepsilon\to 0$ , (2.1) becomes

(2.15)

\displaystyle\varepsilon|\mathbb{S}^{n-1}|\int_{\Omega}\cos(q\cdot x)H(x)~{}dx+\mathcal{O}(\varepsilon^{2})=BT\,,

where one defines the internal data $H$ as

(2.16)

H(x):=\fint_{\mathbb{S}^{n-1}}v\mathcal{L}u~{}d\theta\,.

This means that in the leading order, with the measurement at the boundary, we can recover the quantity

\varepsilon\int_{\Omega}\cos(q\cdot x)H(x)~{}dx\,.

As the reconstruction can be done for all $q\in\mathbb{R}^{n}$ , we obtain the entire spectrum of the Fourier transform of $H$ , which is further used to recover $H$ . ¹¹1We assume that $\sigma_{a},\sigma_{s}$ are compactly supported in $\Omega$ .

To justify (2.15), we use the smallness of $u-u_{\varepsilon}$ described in the following Lemma:

Lemma 2.2.

For small $\varepsilon$ , we have

\|u_{\varepsilon}-u\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\leq{\mathcal{O}(\varepsilon)}\|u_{\varepsilon}\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\,.

Proof.

Since $u$ solves the RTE $\theta\cdot\nabla u=\mathcal{L}u$ and $u_{\varepsilon}$ solves the modulated RTE $\theta\cdot\nabla u_{\varepsilon}=\mathcal{L}_{\varepsilon}u_{\varepsilon}$ with the same incoming boundary conditions, we can write

	$\displaystyle\theta\cdot\nabla(u_{\varepsilon}-u)$	$\displaystyle=\mathcal{L}_{\varepsilon}u_{\varepsilon}-\mathcal{L}u$
		$\displaystyle=\mathcal{L}(u_{\varepsilon}-u)+(\mathcal{L}_{\varepsilon}-\mathcal{L})u_{\varepsilon}$
		$\displaystyle=\mathcal{L}(u_{\varepsilon}-u)+\varepsilon\cos(q\cdot x)\mathcal{L}u_{\varepsilon}$

and $u_{\varepsilon}-u=0$ on $\Gamma_{-}$ . It follows from Theorem 2.1 that

\|u_{\varepsilon}-u\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\leq{\mathcal{O}(\varepsilon)}\|\mathcal{L}u_{\varepsilon}\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\,,

and thus from the definition of $\mathcal{L}$ we derive

\|u_{\varepsilon}-u\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\leq{\mathcal{O}(\varepsilon)}\|u_{\varepsilon}\|_{L^{\infty}(\Omega\times\mathbb{S}^{n-1})}\,.

∎

Remark 2.1.

We comment that the proof above may lead to a large constant in the $\mathcal{O}$ notation for small $\mathsf{Kn}$ , but this constant will not affect the later proof. Furthermore, the bounds shown above is not sharp in the diffusion limit. In the diffusion limit of $\mathsf{Kn}\to 0$ , $u$ and $u_{\varepsilon}$ approximately satisfy two elliptic equations. By comparing these two elliptic equations, we can obtain the smallness of $u-u_{\varepsilon}$ in $\varepsilon$ , and such smallness will be independent of $\mathsf{Kn}$ .

This immediately leads to the following proposition:

Proposition 2.3.

Suppose that $f,\ g\in L^{\infty}(\Gamma_{-})$ . Let $u$ be the solution to the RTE with boundary data $u|_{\Gamma_{-}}=f$ , and $v$ be the solution to the adjoint RTE with boundary data $v|_{\Gamma_{+}}=g$ . Then $\tilde{\mathcal{A}}(g)$ and $\mathcal{A}_{\varepsilon,q}(f)$ uniquely determine the internal functional $H\in L^{\infty}(\Omega)$ up to $\mathcal{O}(\varepsilon)$ accuracy:

(2.17)

\displaystyle H(x)=\frac{\sigma_{s}(x)}{\mathsf{Kn}}\langle u\rangle\langle v\rangle-\frac{\sigma(x)}{\mathsf{Kn}}\langle uv\rangle\,.

This directly comes from inserting the definition of $\mathcal{L}$ (2.4) in that of $H$ in (2.16).

2.2. Main results

The main objective of this paper is to analyze the stability of the reconstruction process of AOI. In particular we are interested in the requirement imposed on the probing light when the diffusion effect is strong. More specifically, we consider the incoming light $f$ that is concentrated in angle:

f(x,\theta):=\left\{\begin{array}[]{ll}c_{n}h^{\frac{1}{2}(1-n)}&\hbox{if }|\theta-\theta_{0}|<h\,,\\ 0&\hbox{otherwise}\,,\\ \end{array}\right.

in $\Omega\times\mathbb{S}^{n-1}$ , where $h$ stands for the concentration of the light beam, and $c_{n}$ is chosen to be the constant depending only on $n$ such that

(2.18)

\displaystyle\fint_{\mathbb{S}^{n-1}}f^{2}d\theta=1\,.

With concentrated light beam injected into the tissue, singularities are more likely to be preserved, keeping more information for the reconstruction. As $\mathsf{Kn}$ decreases, media becomes optically thicker, and the light beam needs to be more concentrated for a stable reconstruction. We will carefully trace such dependence of $h$ on $\mathsf{Kn}$ .

The main theorem we will prove is the following, and we leave the proof to Section 3.

Theorem 2.4.

Let $H$ be the internal functional defined in Proposition 2.3. Suppose that $f$ satisfies (3.4). If $\sigma\in C^{\alpha}(\Omega)$ , for some $0<\alpha\leq 1$ , then for all $x\in\Omega$ one has

(2.19)

\frac{\sigma(x)}{\mathsf{Kn}}=-\frac{f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})}{\mathcal{A}(f)(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})}H(x)+\mathcal{O}(h^{\alpha})e^{\frac{3d\sup\sigma}{\mathsf{Kn}}}\,,

where $\tau_{\pm}$ represent the distance from $x$ to $\partial\Omega$ in the direction $\pm\theta$ , and $d$ is the diameter of the domain $\Omega$ . The constant implied in the big $\mathcal{O}$ is independent of $\mathsf{Kn}$ . Moreover, if one chooses

(2.20)

h\ll e^{\frac{-3d\sup\sigma}{\alpha\mathsf{Kn}}}\,,

then

(2.21)

\frac{\sigma(x)}{\mathsf{Kn}}=-\frac{f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})}{\mathcal{A}(f)(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})}H(x)+\mathcal{O}(h^{\alpha^{\prime}})\,,

for some $0<\alpha^{\prime}<\alpha$ .

We make the following five remarks regarding the interpretation of Theorem 2.4.

First, on explicit reconstruction: we note, as in [27], that everything on the right side of (2.21), other than the error term, is explicit: can be directly calculated from boundary measurements. Thus (2.21) can be viewed as a reconstruction formula for $\sigma$ in terms of boundary quantities.

Second, in the case $\mathsf{Kn}\rightarrow 0$ , we have $\mathsf{Kn}^{2}\sigma_{a}\rightarrow 0$ in the total absorption coefficient $\sigma$ . Therefore it makes sense in this case to consider the reconstruction of one coefficient by (2.21), and not attempt to recover $\sigma_{a}$ and $\sigma_{s}$ separately as in [27].

Third, one preliminary assumption of $\sigma$ : the theorem does require Hölder continuity. This is different from what is required in [27] where Hölder continuity is removed. The main difference is we are establishing a rate of convergence for the formula in terms of the angular concentration of the source, and thus some control on the modulus of continuity is needed (see the proof of Proposition 3.1 for details.) Heuristically, it should be reasonable to concede that a sufficiently rough coefficient cannot be well-resolved with any source that is not completely singular.

Fourth, on incoming data concentration: in the case $\mathsf{Kn}\ll 1$ , (2.20) gives us a condition to guarantee that (2.19) can be used to reconstruct $\sigma$ . In other words, in the case $\mathsf{Kn}\ll 1$ , Theorem 2.4 gives a quantitative condition on the degree of angular concentration of the source required for tomography to work at all, with possible implications for the RTE-based optical tomography more generally.

Finally, on stability: this theorem also implies a quantitative stability estimate in terms of $\mathsf{Kn}$ for the acousto-optic reconstruction. In Section 3, we will show that $\frac{f(x-\tau_{-})}{\mathcal{A}f(x+\tau_{+})}\sim e^{\frac{c}{\mathsf{Kn}}}$ ; this means the perturbation in $H$ will be reflected exponentially in the reconstruction of $\sigma$ . Applying (3.10), this means that in the regime where (2.21) applies, we have

\frac{1}{\mathsf{Kn}}|\sigma_{1}(x)-\sigma_{2}(x)|\leq e^{\frac{d\sup\sigma_{j}}{\mathsf{Kn}}}|H_{1}(x)-H_{2}(x)|+\mathcal{O}(h^{\alpha})\,.

This shows that the stability of the reconstruction decays exponentially as $\mathsf{Kn}\rightarrow 0$ , which should be expected from the reliance on singularities in the RTE (see also [31, 45, 33]).

3. Reconstruction of optical properties

We show Theorem 2.4 in this section. The proof relies on a proper decomposition of the transport solution according to the singularity. This technique, termed the singular decomposition, is the classical technique in the study of optics in inverse problems.

This was introduced in [5, 17, 21, 22, 23] and later was applied extensively to treat various settings related to parameter identification of the RTE [12, 13, 24, 35, 37, 38, 43]. For the purpose of this paper, we will decompose the solution to the RTE into two parts - the most singular (ballistic part) and the remainders, and then separate the singular contribution from the rest of $H(x)$ .

To start, we first write the solution $u\in L^{\infty}(\Omega\times\mathbb{S}^{n-1})$ of the equation (2.3) in the form of

u=u_{1}+u_{2}\,,

where $u_{1}$ is the ballistic part that satisfies the X-ray equation and $u_{2}$ takes care of the remaining terms, namely:

(3.1)

\displaystyle\begin{cases}\theta\cdot\nabla u_{1}=-\frac{\sigma}{\mathsf{Kn}}u_{1}\,,\\ u_{1}|_{\Gamma_{-}}=f\,,\\ \end{cases}\quad\text{and}\hskip 28.45274pt\begin{cases}\theta\cdot\nabla u_{2}=\frac{\sigma_{s}}{\mathsf{Kn}}\langle u_{2}\rangle-\frac{\sigma}{\mathsf{Kn}}u_{2}+\frac{\sigma_{s}}{\mathsf{Kn}}\langle u_{1}\rangle\,,\\ u_{2}|_{\Gamma_{-}}=0\,.\\ \end{cases}

The existence of solutions to both equations follows from Theorem 2.1. Note that the entire boundary condition $f$ is applied on $u_{1}$ leaving $u_{2}$ with the trivial boundary condition.

A similar separation can also be performed to the adjoint RTE (2.13). We write $v=v_{1}+v_{2}$ with the components satisfying:

(3.2)

\displaystyle\begin{cases}-\theta\cdot\nabla v_{1}=-\frac{\sigma}{\mathsf{Kn}}v_{1}\,,\\ v_{1}|_{\Gamma_{+}}=g\,,\\ \end{cases}\quad\text{and}\hskip 28.45274pt\begin{cases}-\theta\cdot\nabla v_{2}=\frac{\sigma_{s}}{\mathsf{Kn}}\langle v_{2}\rangle-\frac{\sigma}{\mathsf{Kn}}v_{2}+\frac{\sigma_{s}}{\mathsf{Kn}}\langle v_{1}\rangle\,,\\ v_{2}|_{\Gamma_{+}}=0\,.\end{cases}

With these decompositions, we now rewrite the internal function $H(x)$ as:

	$\displaystyle H(x)$	$\displaystyle=\frac{\sigma_{s}(x)}{\mathsf{Kn}}(\langle u_{1}\rangle\langle v_{1}\rangle+\langle u_{1}\rangle\langle v_{2}\rangle+\langle u_{2}\rangle\langle v_{1}\rangle+\langle u_{2}\rangle\langle v_{2}\rangle)$
(3.3)			$\displaystyle\hskip 28.45274pt-\frac{\sigma(x)}{\mathsf{Kn}}(\langle u_{1}v_{1}\rangle+\langle u_{1}v_{2}\rangle+\langle u_{2}v_{1}\rangle+\langle u_{2}v_{2}\rangle)\,,$

simply by substituting $u=u_{1}+u_{2}$ and $v=v_{1}+v_{2}$ into the definition of $H$ (2.17).

As described previously, boundary conditions are light beams concentrated in the angle variable $\theta_{0}\in\mathbb{S}^{n-1}$ :

(3.4)

g(x,\theta)=f(x,\theta):=\left\{\begin{array}[]{ll}c_{n}h^{\frac{1}{2}(1-n)}&\hbox{if }|\theta-\theta_{0}|<h\,,\\ 0&\hbox{otherwise}\,,\\ \end{array}\right.

in $\Omega\times\mathbb{S}^{n-1}$ , where $c_{n}$ is the normalization constant.

In the rest of the section we will show that with such concentrated incoming data (3.4) and sufficiently small $h$ , the term $\langle u_{1}v_{1}\rangle$ essentially dominates the other terms in $H(x)$ . This can be seen from the estimates of the dominating and remaining terms presented in Proposition 3.1 and Proposition 3.2, respectively. Finally, from the explicit dependence of $\langle u_{1}v_{1}\rangle$ on $\sigma$ , we can derive the reconstruction of $\sigma$ from the internal data.

Proposition 3.1.

Suppose $\sigma\in C^{\alpha}(\Omega)$ , and let $u_{i}$ and $v_{i}$ , $i=1,2$ , satisfy (3.1) and (3.2) respectively with boundary condition defined in (3.4). Then for each $x$ in $\Omega$ , one has

\langle u_{1}v_{1}\rangle(x)=e^{-\int_{\ell(x,\theta_{0})}\frac{\sigma}{\mathsf{Kn}}ds}(1+\mathsf{Kn}^{-1}\mathcal{O}(h^{\alpha}))\,,

where $\ell(x,\theta)$ denotes the section of the straight line through $x$ in the direction of $\theta$ , intersected with $\Omega$ .

Proof.

The ballistic terms of the solution $u$ to the RTE and the solution $v$ to the adjoint RTE are given explicitly by

(3.5)

\displaystyle u_{1}(x,\theta)=e^{-\int_{0}^{\tau_{-}(x,\theta)}{\frac{1}{\mathsf{Kn}}}\sigma(x-s\theta)ds}f(x-\tau_{-}(x,\theta)\theta,\theta)\,,

(3.6)

\displaystyle v_{1}(x,\theta)=e^{-\int_{0}^{\tau_{+}(x,\theta)}{\frac{1}{\mathsf{Kn}}}\sigma(x+s\theta)ds}g(x+\tau_{+}(x,\theta)\theta,\theta)\,,

where $\tau_{\pm}(x,\theta)$ is the distance from the point $x$ to the boundary in the direction $\pm\theta$ .

Substituting these into $\langle u_{1}v_{1}\rangle(x)$ and using (3.4), we find that for any $x\in\Omega$ ,

\langle u_{1}v_{1}\rangle(x)=\fint_{\mathbb{S}^{n-1}}e^{-\int_{-\tau_{-}(x,\theta)}^{\tau_{+}(x,\theta)}{\frac{1}{\mathsf{Kn}}}\sigma(x+s\theta)ds}c_{n}^{2}h^{1-n}\chi^{h}_{\theta_{0}}(\theta)d\theta\,,

where $\chi^{h}_{\theta_{0}}(\theta)$ is the characteristic function that is equal to $1$ in $|\theta-\theta_{0}|<h$ and equal to $0$ otherwise. The integral in the exponent is precisely the integral over $\ell(x,\theta)$ , so we can rewrite this as

\langle u_{1}v_{1}\rangle(x)=c_{n}^{2}h^{1-n}\fint_{\mathbb{S}^{n-1}}e^{-{\frac{1}{\mathsf{Kn}}}\int_{\ell(x,\theta)}\sigma ds}\chi^{h}_{\theta_{0}}(\theta)d\theta\,.

From (2.18) for small $h$ , we have

(3.7)

\langle u_{1}v_{1}\rangle(x)=e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds}+R\,,

where $R$ is the error term

	$\displaystyle R$	$\displaystyle=\fint_{\mathbb{S}^{n-1}}(e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta)}\sigma ds}-e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds})c_{n}^{2}h^{1-n}\chi^{h}_{\theta_{0}}(\theta)d\theta$
		$\displaystyle=e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds}\fint_{\mathbb{S}^{n-1}}(e^{-\frac{1}{\mathsf{Kn}}(\int_{\ell(x,\theta)}\sigma ds-\int_{\ell(x,\theta_{0})}\sigma ds)}-1)c_{n}^{2}h^{1-n}\chi^{h}_{\theta_{0}}(\theta)d\theta\,.$

Since $\sigma\in C^{\alpha}(\Omega)$ , we have for sufficiently small $h$

R=e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds}\fint_{\mathbb{S}^{n-1}}(e^{-\frac{\mathcal{O}(h^{\alpha})}{\mathsf{Kn}}}-1)c_{n}^{2}h^{1-n}\chi^{h}_{\theta_{0}}(\theta)d\theta\,.

Note that the constant in the $\mathcal{O}$ notation merely depends on the regularity of $\sigma$ . We can estimate $|e^{-\frac{\mathcal{O}(h^{\alpha})}{\mathsf{Kn}}}-1|$ by the Taylor remainder formula if $\frac{\mathcal{O}(h^{\alpha})}{\mathsf{Kn}}$ is small, or by $1$ otherwise. In either case, we get

R=e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds}\frac{\mathcal{O}(h^{\alpha})}{\mathsf{Kn}}\fint_{\mathbb{S}^{n-1}}c_{n}^{2}h^{1-n}\chi^{h}_{\theta_{0}}(\theta)d\theta\,.

Then (2.18) implies that

R=\frac{1}{\mathsf{Kn}}e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds}\mathcal{O}(h^{\alpha})\,.

Returning to (3.7), we conclude that

\langle u_{1}v_{1}\rangle(x)=e^{-\frac{1}{\mathsf{Kn}}\int_{\ell(x,\theta_{0})}\sigma ds}(1+\mathsf{Kn}^{-1}\mathcal{O}(h^{\alpha}))

as desired. ∎

The following proposition justifies the smallness of the remainder terms.

Proposition 3.2.

Under the same conditions as in Proposition 3.1, suppose $\mathsf{Kn}\leq 1$ , there exists a constant $C$ independent of $h$ and $\mathsf{Kn}$ such that

\langle u_{j}\rangle\langle v_{k}\rangle(x)\leq C\mathsf{Kn}^{-4}h^{n-1}

for $j,k=1,2$ , and

\langle u_{j}v_{k}\rangle(x)\leq C\mathsf{Kn}^{-4}h^{n-1}

for $j,k=1,2$ and $j+k>2$ for each $x$ in $\Omega$ .

Proof.

First note that it follows directly from (3.5), (3.4), and $\sigma>0$ , the angular average of $u_{1}$ satisfies

(3.8)

\langle u_{1}\rangle\leq\mathcal{O}(h^{\frac{n-1}{2}})\ \hbox{for all }x\in\Omega\,.

Moreover, Theorem 2.1 and (3.1) tell us that

\|u_{2}\|_{L^{\infty}(\Omega\times S^{n-1})}\leq C\mathsf{Kn}^{-2}\left\|\langle u_{1}\rangle\right\|_{L^{\infty}(\Omega)}.

Combining with (3.8), we can deduce that

(3.9)

\|u_{2}\|_{L^{\infty}(\Omega\times S^{n-1})}\leq C\mathsf{Kn}^{-2}\mathcal{O}(h^{\frac{n-1}{2}})\,.

The estimates (3.8) and (3.9) also hold for $v_{1}$ and $v_{2}$ , respectively. Together these estimates yield all of the bounds of Proposition 3.2. For instance, applying (3.8) to $\langle u_{1}\rangle$ and (3.9) to $\|v_{2}\|_{L^{\infty}(\Omega\times S^{n-1})}$ we get

\langle u_{1}v_{2}\rangle\leq\langle u_{1}\rangle\|v_{2}\|_{L^{\infty}(\Omega\times S^{n-1})}\leq\mathsf{Kn}^{-2}\mathcal{O}(h^{n-1})\,.

Similarly, we can derive

\langle u_{2}v_{2}\rangle\leq C\mathsf{Kn}^{-4}\mathcal{O}(h^{n-1})\,,

where $C$ is a constant independent of $\mathsf{Kn}$ and $h$ . ∎

The main result is a natural consequence of the previous two propositions.

Proof of Theorem 2.4.

We first rewrite (3) as

		$\displaystyle H(x)+\frac{\sigma(x)}{\mathsf{Kn}}\langle u_{1}v_{1}\rangle$
		$\displaystyle=\frac{\sigma_{s}(x)}{\mathsf{Kn}}(\langle u_{1}\rangle\langle v_{1}\rangle+\langle u_{1}\rangle\langle v_{2}\rangle+\langle u_{2}\rangle\langle v_{1}\rangle+\langle u_{2}\rangle\langle v_{2}\rangle)$
		$\displaystyle\hskip 28.45274pt-\frac{\sigma(x)}{\mathsf{Kn}}(\langle u_{1}v_{2}\rangle+\langle u_{1}v_{2}\rangle+\langle u_{2}v_{2}\rangle)\,.$

Applying Proposition 3.1 to the $\langle u_{1}v_{1}\rangle$ term and Proposition 3.2 to everything else, for $\mathsf{Kn}\ll 1$ , we get

H(x)+\frac{\sigma(x)}{\mathsf{Kn}}\left(e^{-{\frac{1}{\mathsf{Kn}}}\int_{\ell(x,\theta_{0})}\sigma ds}(1+\mathsf{Kn}^{-1}\mathcal{O}(h^{\alpha}))\right)=\mathsf{Kn}^{-5}\mathcal{O}(h^{n-1})\,.

Multiplying through by $e^{\int_{\ell(x,\theta_{0})}{\frac{\sigma}{\mathsf{Kn}}}ds}\,$ , one has

\frac{\sigma(x)}{\mathsf{Kn}}+e^{\int_{\ell(x,\theta_{0})}{\frac{\sigma}{\mathsf{Kn}}}ds}H(x)=\frac{1}{\mathsf{Kn}^{2}}\mathcal{O}(h^{\alpha})+\frac{1}{\mathsf{Kn}^{5}}\mathcal{O}(h^{n-1})e^{\frac{d\sup\sigma}{\mathsf{Kn}}}\,.

In the zero limit of $\mathsf{Kn}$ , the exponential term dominates, and the equation becomes

(3.10)

\frac{\sigma(x)}{\mathsf{Kn}}+e^{\int_{\ell(x,\theta_{0})}{\frac{\sigma}{\mathsf{Kn}}}ds}H(x)\leq\mathcal{O}(h^{\alpha})e^{\frac{2d\sup\sigma}{\mathsf{Kn}}}\,.

To estimate the exponential factor in front of $H(x)$ , we recall that from (3.5), $u_{1}$ has explicit formula:

u_{1}(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})=e^{-{\frac{1}{\mathsf{Kn}}}\int_{\ell(x,\theta_{0})}\sigma ds}f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})\,,

and from (3.9) we also have:

|u_{2}(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})|\leq\mathsf{Kn}^{-2}\mathcal{O}(h^{\frac{n-1}{2}})\,.

Considering the definition of the albedo operator:

\mathcal{A}(f)(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})=u(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})=u_{1}(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})+u_{2}(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})\,,

we finally have

	$\displaystyle\frac{\mathcal{A}(f)(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})}{f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})}$	$\displaystyle=e^{-{\frac{1}{\mathsf{Kn}}}\int_{\ell(x,\theta_{0})}\sigma ds}+\frac{\mathcal{O}(h^{\frac{n-1}{2}})}{\mathsf{Kn}^{2}f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})}$
(3.11)			$\displaystyle=e^{-{\frac{1}{\mathsf{Kn}}}\int_{\ell(x,\theta_{0})}\sigma ds}+\frac{\mathcal{O}(h^{n-1})}{\mathsf{Kn}^{2}}\,,$

where we also used the fact that $f|_{\Gamma_{-}}$ defined in (3.4). By substituting this back to (3.10), it immediately suggests:

\frac{\sigma(x)}{\mathsf{Kn}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})}{\mathcal{A}(f)(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})}}H(x)=\mathcal{O}(h^{\alpha})e^{\frac{2d\sup\sigma}{\mathsf{Kn}}}+\frac{\mathcal{O}(h^{n-1})}{\mathsf{Kn}^{2}}H(x)\,.

Since (3.10) also implies that $H(x)$ is of order $e^{\frac{2d\sup\sigma}{\mathsf{Kn}}}$ , we conclude

\frac{\sigma(x)}{\mathsf{Kn}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{f(x-\tau_{-}(x,\theta_{0})\theta_{0},\theta_{0})}{\mathcal{A}(f)(x+\tau_{+}(x,\theta_{0})\theta_{0},\theta_{0})}}H(x)=\mathcal{O}(h^{\alpha})e^{\frac{3d\sup\sigma}{\mathsf{Kn}}}\,,

Thus, Theorem 2.4 now follows. ∎

4. Numeric Examples

In this section we provide numerical evidence that verifies the analysis above.

As a numerical setup, we let $\Omega=[0,1]^{2}$ , and mesh grids are sampled to resolve small scales in $\mathsf{Kn}$ . Since we only concern ourselves with the behavior of the different terms in $H(x)$ in terms of $\mathsf{Kn}$ and $h$ , we use the simple media $\sigma=1$ all over the domain. In particular, the incoming light is concentrated at $(x=0,y=1/2)$ with $h$ indicating the concentration in velocity domain. In Figure 1 we plot the velocity average of the ballistic part, $\langle u_{1}\rangle$ , for different $\mathsf{Kn}$ with $h=\frac{2\pi}{25}$ and $h=\frac{14\pi}{25}$ , and it can be clearly seen that smaller $h$ narrows down the spreading of $\langle u_{1}\rangle$ while smaller $\mathsf{Kn}$ gives stronger decay of $\langle u_{1}\rangle$ in space. Moreover, we also plot $\langle u_{2}\rangle$ in Figure 2. Similar to the case of $\langle u_{1}\rangle$ , as $h$ increases the function is more spread out, while smaller $\mathsf{Kn}$ leads to the stronger decay.

We also plot $\langle u_{1}v_{1}\rangle$ as a function of $\mathsf{Kn}$ and $h$ . As seen in Figure 3, the quantity decays exponentially fast with respect to small $\mathsf{Kn}$ and algebraically fast with respect to small $h$ , which confirms with the result in Proposition 3.1. In addition, the remainder term $H(x)-\frac{1}{\mathsf{Kn}}\langle u_{1}v_{1}\rangle$ is plotted in Figure 4.

Lastly, the quantitative dependence of error on $\mathsf{Kn}$ and $h$ is plotted in Figure 5 and the observation confirms with the prediction of our theorem. Specifically, this error is an exponential function of $\frac{1}{\mathsf{Kn}}$ and grows algebraically with respect to $h$ , where the error is defined to be the relative error, namely:

\text{Error}(x)=\frac{\mathsf{Kn}H(x)-\langle u_{1}v_{1}\rangle}{\langle u_{1}v_{1}\rangle}

and the plot shows the $L^{2}$ norm over the space $x\in\Omega$ .

Refer to caption — Figure 1. The plot show $\langle u_{1}\rangle$ with different $\mathsf{Kn}$ and $h$ . From left to right, $\mathsf{Kn}$ decreases as $2^{-k}$ with $k=0,1,2,3,4$ . The two rows are for $h$ being $\frac{2\pi}{25}$ and $\frac{14\pi}{25}$ respectively.

References

[1] G. Alessandrini. Stable determination of conductivity by boundary measurements. Appl. Anal., 27:153–172, 1988.
[2] H. Ammari, E. Bossy, J. Garnier, L. H. Nguyen, and L. Seppecher. A reconstruction algorithm for ultrasound-modulated diffuse optical tomography. Proc. Amer. Math. Soc., 142:3221–3236, 2014.
[3] H. Ammari, J. Garnier, L. H. Nguyen, and L. Seppecher. Reconstruction of a piecewise smooth absorption coefficient by an acousto-optic process. Commun. Part. Diff. Eq., 38:1737–1762, 2013.
[4] H. Ammari, L. H. Nguyen, and L. Seppecher. Reconstruction and stability in acousto-optic imaging for absorption maps with bounded variation. J. Funct. Anal., 267:4361–4398, 2014.
[5] D. Anikonov, I. Prokhorov, and A. Kovtanyuk. Investigation of scattering and absorbing media by the methods of X-ray tomography. Journal of Inverse and Ill-posed Problems, 1:259–281, 1993.
[6] S. Arridge. Optical tomography in medical imaging. Inverse Problems, 15:R41–93, 1999.
[7] S. R. Arridge and J. C. Schotland. Optical tomography: forward and inverse problems. Inverse Problems, 25(12):123010, 2009.
[8] G. Bal. Hybrid inverse problems and internal functionals, Inside Out II, MSRI Publications, G. Uhlmann Editor. Cambridge University Press, 2012.
[9] G. Bal, F. J. Chung, and J. C. Schotland. Ultrasound modulated bioluminescence tomography and controllability of the radiative transport equation. SIAM J. Math. Analysis, 48(2):1332–1347, 2016.
[10] G. Bal and A. Jollivet. Time-dependent angularly averaged inverse transport. Inverse Problems, 25(7):075010, 2009.
[11] G. Bal and A. Jollivet. Stability for time-dependent inverse transport. SIAM J. Math. Anal., 42(2):679–700, 2010.
[12] G. Bal and A. Jollivet. Generalized stability estimates in inverse transport theory. Inverse problems and Imaging, 12(1):59–90, 2018.
[13] G. Bal and F. Monard. Inverse transport with isotropic time-harmonic sources. SIAM Journal on Mathematical Analysis, 44(1):134–161, 2012.
[14] G. Bal and S. Moskow. Local inversions in ultrasound modulated optical tomography. Inverse Problems, 30(2):025005, 2014.
[15] G. Bal and J. C. Schotland. Inverse scattering and acousto-optic imaging. Phys. Rev. Letters, 104:043902, 2010.
[16] G. Bal and J. C. Schotland. Ultrasound-modulated bioluminescence tomography. Phys. Rev. E, 89:031201, 2014.
[17] A. Bondarenko. Singular structure of the fundamental solution of the transport equation, and inverse problems in particle scattering theory. Dokl. Akad. Nauk SSSR, 322:274–276, 1992.
[18] K. Chen, Q. Li, and J.-G. Liu. Online learning in optical tomography: a stochastic approach. Inverse Problems, 34(7):075010, may 2018.
[19] K. Chen, Q. Li, and L. Wang. Stability of stationary inverse transport equation in diffusion scaling. Inverse Problems, 34(2), 2018.
[20] Y. Cheng, I. M. Gamba, and K. Ren. Recovering doping profiles in semiconductor devices with the Boltzmann-Poisson model. J. Comput. Phys., 230(9):3391–3412, May 2011.
[21] M. Choulli and P. Stefanov. Scattering inverse pour l’équation du transport et relations entre les opérateurs de scattering et d’albédo. C. R. Acad. Sci. Paris, 320:947–952, 1995.
[22] M. Choulli and P. Stefanov. Inverse scattering and inverse boundary value problems for the linear Boltzmann equation. Comm. P.D.E., 21:763–785, 1996.
[23] M. Choulli and P. Stefanov. Reconstruction of the coefficients of the stationary transport equation from boundary measurements. Inverse Problems, 12:L19–L23, 1996.
[24] M. Choulli and P. Stefanov. An inverse boundary value problem for the stationary transport equation. Osaka J. Math., 36:87–104, 1998.
[25] F. Chung, J. Hoskins, and J. Schotland. Coherent acousto-optic tomography with diffuse light. preprint, 2018.
[26] F. Chung, J. Hoskins, and J. Schotland. A transport model for multi-frequency acousto-optic tomography. arXiv:1910.04798, 2019.
[27] F. Chung and J. Schotland. Inverse transport and acousto-optic imaging. SIAM J. Math. Analysis, 49:4704–4721, 2017.
[28] D. S. Elson, R. Li, C. Dunsby, R. Eckersley, and M.-X. Tang. Ultrasound-mediated optical tomography: a review of current methods. Interface Focus, 1:632–648, 2011.
[29] B. Jin, Z. Zhou, and J. Zou. On the convergence of stochastic gradient descent for nonlinear ill-posed problems. arXiv:1907.03132, 2019.
[30] P. Kuchment. Mathematics of Hybrid Imaging: A Brief Review. In: Sabadini I., Struppa D. (eds), volume 16. Springer Proceedings in Mathematics, 2012.
[31] R.-Y. Lai, Q. Li, and G. Uhlmann. Inverse problems for the stationary transport equation in the diffusion scaling. SIAM Journal on Applied Mathematics, 79(6):2340–2358, 2019.
[32] M. Lesaffrea, F. Jeana, A. Funkea, P. Santosa, M. Atlana, B. Forgeta, E. Bossya, F. Ramaza, A. Boccaraa, M. Grossb, P. Delayec, and G. Roosenc. Acousto-optic imaging techniques for optical diagnosis. IFAC Proceedings Volumes, 39(18):11–15, 2006.
[33] W. Li, Y. Yang, and Y. Zhong. A hybrid inverse problem in the fluorescence ultrasound modulated optical tomography in the diffusive regime. SIAM Journal on Applied Mathematics, 79(1):356–376, 2019.
[34] L. D. Montejo, J. Jia, H. K. Kim, U. J. Netz, S. Blaschke, G. A. Muller, and A. H. Hielscher. Computer-aided diagnosis of rheumatoid arthritis with optical tomography, part 2: image classification. Journal of Biomedical Optics, 18(7):076002–076002, 2013.
[35] K. Ren. Recent developments in numerical techniques for transport-based medical imaging methods. Commun. Comput. Phys., 8(1):1–50, 2010.
[36] K. Ren, R. Zhang, and Y. Zhong. Inverse transport problems in quantitative pat for molecular imaging. Inverse Problems, 31(12):125012, 2015.
[37] P. Stefanov. Inverse problems in transport theory, volume 47. Inside Out: Inverse Problems; MSRI Publications, edited by G. Uhlmann, 2003.
[38] P. Stefanov and G. Uhlmann. Optical tomography in two dimensions. Methods Appl. Anal., 10:1–9, 2003.
[39] T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. R. Arridge. Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography. Inverse Problems, 28(8):084009, 2012.
[40] T. Tarvainen, V. Kolehmainen, S. R. Arridge, and J. P. Kaipio. Image reconstruction in diffuse optical tomography using the coupled radiative transport–diffusion model. Journal of Quantitative Spectroscopy and Radiative Transfer, 112(16):2600 – 2608, 2011.
[41] T. Tarvainen, V. Kolehmainen, A. Pulkkinen, M. Vauhkonen, M. Schweiger, S. R. Arridge, and J. P. Kaipio. An approximation error approach for compensating for modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography. Inverse Problems, 26(1):015005, 2010.
[42] G. Uhlmann. Electrical impedance tomography and Calderón’s problem. Inverse problems, 25(12):123011, 2009.
[43] J.-N. Wang. Stability estimates of an inverse problem for the stationary transport equation. Ann. Inst. H. Poincaré Phys. Théor., 70(5):473–495, 1999.
[44] L. Wang. Ultrasound-mediated biophotonic imaging: A review of acousto-optical tomography and photo-acoustic tomography. Disease Markers, 19:123–138, 2003, 2004.
[45] H. Zhao and Y. Zhong. Instability of an inverse problem for the stationary radiative transport near the diffusion limit. arXiv:1809.01790, 2018.