A Fourier approach to the inverse source problem in an absorbing and anisotropic scattering medium

TL;DR

This paper introduces a non-iterative Fourier-based method for reconstructing an unknown isotropic source in a 2D absorbing, scattering medium, extending Bukhgeim's approach to complex media with practical optical tomography applications.

Contribution

It develops a novel Fourier approach for inverse source problems in scattering media, generalizing previous non-scattering techniques to more realistic optical environments.

Findings

Successful numerical demonstration with Henyey-Greenstein scattering kernel

Extension of Bukhgeim's method to scattering media

Feasibility shown for optical tomography scenarios

Abstract

We revisit the inverse source problem in a two dimensional absorbing and scattering medium and present a non-iterative reconstruction method using measurements of the radiating flux at the boundary. The attenuation and scattering coefficients are known and the unknown source is isotropic. The approach is based on the Cauchy problem for a Beltrami-like equation for the sequence valued maps, and extends the original ideas of A. Bukhgeim from the non-scattering to scattering media. We demonstrate the feasibility of the method in a numerical experiment in which the scattering is modeled by the two dimensional Henyey-Greenstein kernel with parameters meaningful in Optical Tomography.

Equations

$$\mathbf{\theta }\cdot \nabla u(z,\mathbf{\theta })+a(z)u(z,\mathbf{\theta })$$

$$u\in H^1(\Omega \times {\mathbf{S}}^1).$$

$$\Omega \ni z\mapsto \mathbf{u}(z)$$

$$\overline{\partial }{u}_1(z)+\partial u_{-1}(z)+a(z)u_0(z)=k_0(z)u_0(z)+f(z),$$

$$\overline{\partial }{u}_n(z)+\partial u_{n-2}(z)+a(z)u_{n-1}(z)=k_{n-1}(z)u_{n-1}(z),n\ne 1.$$

$$\overline{\partial }\mathbf{u}(z)+L^2\partial \mathbf{u}(z)+a(z)L\mathbf{u}(z)=LK\mathbf{u}(z),z\in \Omega ,$$

$$K\mathbf{u}(z):=(k_0(z)u_0(z),k_{-1}(z)u_{-1}(z),k_{-2}(z)u_{-2}(z),...)$$

$$\mathbf{g}=\mathbf{u}{|}_{\Gamma }=⟨g_0,g_{-1},g_{-2},...⟩.$$

$$\overline{\partial }\mathbf{u}+L^2\partial \mathbf{u}=0,$$

$$u_{-n}(\zeta )=\frac{1}{2\pi i}{\int }_{\Gamma }\frac{u_{-n}(z)}{z-\zeta }dz+\frac{1}{2\pi i}{\int }_{\Gamma }\left\{\frac{dz}{z-\zeta }-\frac{d\overline{z}}{\overline{z}-\overline{\zeta }}\right\}\sum _{j=1}^{\infty }{u}_{-n-2j}(z){\left(\frac{\overline{z-\zeta }}{z-\zeta }\right)}^j,\zeta \in \Omega .$$

$$h(z,\mathbf{\theta }):=Da(z,\mathbf{\theta })-\frac{1}{2}\left(I-iH\right)Ra(z\cdot {\mathbf{\theta }}^{\perp },{\mathbf{\theta }}^{\perp }),$$

$$\mathbf{\theta }\cdot \nabla h(z,\mathbf{\theta })=-a(z),(z,\mathbf{\theta })\in \Omega \times {\mathbf{S}}^1.$$

$$e^{-h(z,\theta )}:=\sum _{k=0}^{\infty }{\alpha }_k(z)e^{ik\phi },e^{h(z,\theta )}:=\sum _{k=0}^{\infty }{\beta }_k(z)e^{ik\phi },(z,\theta )\in \overline{\Omega }\times {\mathbf{S}}^1,$$

$$z\mapsto \mathbf{\alpha }(z):=⟨{\alpha }_0(z),{\alpha }_1(z),{\alpha }_2(z),{\alpha }_3(z),...,⟩\in C^{p,s}(\Omega ;l_1)\cap C(\overline{\Omega };l_1),$$

$$z\mapsto \mathbf{\beta }(z):=⟨{\beta }_0(z),{\beta }_1(z),{\beta }_2(z),{\beta }_3(z),...,⟩\in C^{p,s}(\Omega ;l_1)\cap C(\overline{\Omega };l_1).$$

$$\overline{\partial }{\beta }_0(z)=0,\overline{\partial }{\beta }_1(z)=-a(z){\beta }_0(z),$$

$$\overline{\partial }{\beta }_{k+2}(z)+\partial {\beta }_k(z)+a(z){\beta }_{k+1}(z)=0,k\ge 0.$$

$$\overline{\partial }{\alpha }_0(z)=0,\overline{\partial }{\alpha }_1(z)=a(z){\alpha }_0(z),$$

$$\overline{\partial }{\alpha }_{k+2}(z)+\partial {\alpha }_k(z)-a(z){\alpha }_{k+1}(z)=0,k\ge 0.$$

$${\alpha }_0{\beta }_0=1,\sum _{m=0}^k{\alpha }_m{\beta }_{k-m}=0,k\ge 1.$$

$$(e^{-G}\mathbf{u}{)}_m=(\mathbf{\alpha }\ast \mathbf{u}{)}_m=\sum _{k=0}^{\infty }{\alpha }_k{u}_{m-k},\text{and}(e^G\mathbf{u}{)}_m=(\mathbf{\beta }\ast \mathbf{u}{)}_m=\sum _{k=0}^{\infty }{\beta }_k{u}_{m-k},$$

$$e^{-G}\mathbf{u}=\sum _{k=0}^{\infty }{\alpha }_k{L}^k\mathbf{u},\text{and}{e}^G\mathbf{u}=\sum _{k=0}^{\infty }{\beta }_k{L}^k\mathbf{u},$$

$$l^{2,p}(\mathbb{N};H^q(\Omega ))$$

$${\Vert \mathbf{u}\Vert }^2:={\Vert \mathbf{u}\Vert }_{0,0}^2=\sum _{j=0}^{\infty }{\Vert u_{-j}\Vert }_{L^2(\Omega )}^2.$$

$${\Vert \mathbf{g}\Vert }_{p,\frac{1}{2}}^2:=\sum _{j=0}^{\infty }(1+j^2{)}^p{\Vert g_{-j}\Vert }_{H^{1/2}(\Gamma )}^2.$$

$${u_{-j}|}_{\Gamma }(e^{i\beta })=\sum _{k=-\infty }^{\infty }{u}_{-j,k}{e}^{ik\beta },\text{for}{e}^{i\beta }\in \Gamma ,$$

$${\Vert u_{-j}\Vert }_{H^{1/2}(\Gamma )}^2=\sum _{k=-\infty }^{\infty }(1+k^2{)}^{\frac{1}{2}}|u_{-j,k}{|}^2.$$

$${\Vert \mathbf{g}\Vert }_{\frac{1}{2},\frac{1}{2}}^2=\sum _{j=0}^{\infty }\sum _{n=-\infty }^{\infty }(1+j^2{)}^{\frac{1}{2}}(1+n^2{)}^{\frac{1}{2}}|g_{-j,n}{|}^2.$$

$${\Vert \mathbf{u}\Vert }^2\le \mu \left({\Vert \partial \mathbf{u}\Vert }^2+{\Vert \mathbf{u}{|}_{\Gamma }\Vert }_{0,\frac{1}{2}}^2\right),$$

$$l_{\infty }^{1,1}(\overline{\Omega }):=\left\{\mathbf{\alpha }:=⟨{\alpha }_0,{\alpha }_1,{\alpha }_2,...,⟩:{\Vert \mathbf{\alpha }\Vert }_{l_{\infty }^{1,1}(\overline{\Omega })}:=\underset{z\in \overline{\Omega }}{\sup }\sum _{j=1}^{\infty }j|{\alpha }_j(z)|<\infty \right\}.$$

$$e^{\pm G}:l^{2,p}(\mathbb{N};H^q(\Omega ))\to l^{2,p}(\mathbb{N};H^q(\Omega ))$$

$$\Vert e^{-G}\mathbf{u}\Vert$$

Full text

1. Introduction

This work concerns a Fourier approach to the inverse source problem for radiative transport in a strictly convex domain $\Omega$ in the Euclidean plane. The attenuation and scattering coefficients are known real valued functions. Generated by an unknown source $f$, in the steady state case, the density of particles $u(z,\boldsymbol{\theta})$ at $z$ traveling in the direction $\boldsymbol{\theta}$ solve the stationary transport equation

[TABLE]

where ${{\mathbf{S}}^{1}}$ denotes the unit sphere.

Let $\varGamma_{\pm}:=\{(z,\boldsymbol{\theta})\in\varGamma\times{{\mathbf{S}}^{1}}:\,\pm\nu(z)\cdot\boldsymbol{\theta}>0\}$ be the incoming (-), respectively outgoing (+), unit tangent sub-bundles of the boundary; where $\nu(z)$ is the outer unit normal at $z\in\varGamma$. The (forward) boundary value problem for (1) assumes a given incoming flux $u$ on ${\varGamma_{-}}$, In here we assume that there is no incoming radiation from outside the domain, $\displaystyle u|_{\Gamma_{-}}=0$. The boundary value problem is know to be well-posed under various admissibility and subcritical assumptions, e.g, in [12, 10, 11, 1, 26], with the most general result for a generic pair of coefficients obtained by Stefanov and Uhlmann [46]. In here we assume that the forward problem is well-posed, and that the outgoing radiation $u|_{\Gamma_{+}}$ is measured, and thus the trace $u\lvert_{\Gamma\times{{\mathbf{S}}^{1}}}$ is known. Without loss of generality $\Omega$ is the unit disc.

In here we show how to recover $f$ from knowledge of $u$ on the torus $\Gamma\times{{\mathbf{S}}^{1}}$ and provide an error, and stability estimates.

When $a=k=0$, this is the classical $X$-ray tomography problem of Radon [37], where $f$ is to be recovered from its integrals along lines, see also [32, 18, 25]. For $a\neq 0$ but $k=0$, this is the problem of inversion of the Attenuated Radon transform in two dimensions, solved successfully by Arbuzov, Bukhgeim and Kazantsev  [2], and Novikov [34]; see [33, 8, 5] for later approaches.

The inverse source problem in an absorbing and scattering media, $a,k\neq 0$, has also been considered (e.g., [24, 44]) in the Euclidean setting, and in [42] in the Riemannian setting. The most general result ($k$ may vary with two independent directions) on the stable determination of the source was obtained by Stefanov and Uhlmann [46]. The reconstruction of the source based on [46] is yet to be realized. When the anisotropic part of scattering is sufficiently small, a convergent iterative method for source reconstruction was proposed in [7]. Based on a perturbation argument to the non-scattering case in [34], the method does not extend to strongly anisotropic scattering. In addition, it requires solving one forward problem (a computationally extensive effort) at each iteration.

The main motivation of this work is to provide a source reconstruction method that applies to the anisotropic scattering media, with non-small anisotropy. In here we propose such a non-iterative method. Our approach extends the original ideas in [2] from the non-scattering to the scattering media.

Throughout we assume that $a\in C^{2,s}(\overline{\Omega})$, and $k$ and its angular derivative are periodic in the angular variable, $k\in C^{1,s}_{per}([-1,1];Lip(\Omega))$, $s>1/2$, and that the forward problem is well posed. It is known from [46] that for pairs of coefficients $(a,k)$ in an open and dense sets of $C^{2}\times C^{2}$, and for any $f\in L^{2}(\Omega)$, there is a unique solution $u\in L^{2}(\Omega\times{{\mathbf{S}}^{1}})$ to the forward boundary value problem. However, our approach requires a smooth solution $u\in H^{1}(\Omega\times{{\mathbf{S}}^{1}})$. As a direct consequence of [46, Proposition 3.4] the regularity of $u$ is dictated by its ballistic term. In particular, if $f\in H^{1}(\Omega)$, then $u\in H^{1}(\Omega\times{{\mathbf{S}}^{1}})$. With the exception of the numerical examples in Section 7, we assume that the unknown source $f\in H^{1}(\Omega)$, and thus the unknown solution

[TABLE]

In the numerical experiment we use a discontinuous source, whose successful quantitative reconstruction indicates robustness of the method.

Let $u(z,\boldsymbol{\theta})=\sum_{-\infty}^{\infty}u_{n}(z)e^{in\theta}$ be the formal Fourier series representation of $u$ in the angular variable $\boldsymbol{\theta}=(\cos\theta,\sin\theta)$. Since $u$ is real valued, $u_{-n}=\overline{u_{n}}$ and the angular dependence is completely determined by the sequence of its nonpositive Fourier modes

[TABLE]

Let $k_{n}(z)=\frac{1}{2\pi}\int_{-\pi}^{\pi}k(z,\cos\theta)e^{-in\theta}d\theta$, $n\in\mathbb{Z}$, be the Fourier coefficients of the scattering kernel. Since $k(z,\cos\theta)$ is both real valued and even in $\theta$, $k_{n}(z)$ are real valued and $k_{n}(z)=k_{-n}(z)$.

Throughout this paper the Cauchy-Riemann operators $\overline{\partial}=(\partial_{x}+i\partial_{y})/2$ and $\partial=(\partial_{x}-i\partial_{y})/2$ refer to derivatives in the spatial domain. By using the advection operator $\boldsymbol{\theta}\cdot\nabla=e^{-i\theta}\overline{\partial}+e^{i\theta}\partial$, and identifying the Fourier coefficients of the same order, the equation (1) reduces to the system:

[TABLE]

and

[TABLE]

In particular, the sequence valued map (3) solves the Beltrami-like equation

[TABLE]

where $L{\bf u}(z)=L(u_{0}(z),u_{-1}(z),u_{-2}(z),...):=(u_{-1}(z),u_{-2}(z),...)$ denotes the left translation, and

[TABLE]

is a Fourier multiplier operator determined by the scattering kernel.

Our data $u\lvert_{\varGamma\times{{\mathbf{S}}^{1}}}$ yields the trace of a solution of (6) on the boundary,

[TABLE]

Bukhgeim’s original theory in [9] concerns solutions of (6) for $a=0$ and $K=0$. Solutions of

[TABLE]

(called $L^{2}$-analytic) satisfy a Cauchy-like integral formula, which recovers ${\bf u}$ in $\Omega$ from its trace ${\bf u}|_{\varGamma}$. In the explicit form in [15], for each $n\geq 0$,

[TABLE]

In Section 2 we review the absorbing, non-scattering case. While we follow the treatment in [38], it is in this section that the new analytical framework and notation is introduced. Section 3 describe the reconstruction method for scattering kernels of polynomial dependence in the angular variable. Except for the numerical section in the end, the remaining of the paper analyzes the error made by the polynomial approximation of the scattering kernel. In Section 4 we exhibit the gain in smoothness due to the scattering, in particular, the $1/2$-gain in (61) below has been known (with different proofs) see [26], and in a more general case than considered here in [46].

The key ingredient in our analysis is an a priori gradient estimate for solutions of the inhomogeneous Bukhgeim-Beltrami equations, see Theorem 4.2 in Section 4. Our starting point is an energy identity, an idea originated in the work of Mukhometov [30], and an equivalent of Pestov’s identity [42, 49] for the Bukhgeim-Beltrami equation. The proof of the gradient estimate in Theorem 4.2 uses essentially $1$ derivative gain in smoothness due to scattering. As a consequence, in Theorem 6.1 we establish an error estimate, which yields a stability result for scattering with polynomial angular dependence (see Corollary 6.1). Furthermore, in a weakly anisotropic scattering medium the method is convergent (see Theorem 6.2), thus recovering the result in [7].

The feasibility of the proposed method is implemented in two numerical experiments in Section 7. Among the several models for the scattering kernel used in Optical Tomography [3], we work with the two dimensional version of the Henyey-Greenstein kernel for its simplicity. In this kernel we chose the anisotropy parameter to be $1/2$ ( half way between the ballistic and isotropic regime), and a mean free path of $1/5$ units of length e.g. (meaningful value for fluorescent light scattering in the fog).

2. A brief review of the absorbing non-scattering medium

In the case for $a\neq 0$ and $K=0$, the Beltrami equation (6) can be reduced to (9) via an integrating factor. While this idea originates in [2], in here (as in [38]) we use the special integrating factor proposed by Finch in [15], which enjoys the crucial property of having vanishing negative Fourier modes. This special integrating factor is $e^{-h}$, where

[TABLE]

and $\boldsymbol{\theta}^{\perp}$ is orthogonal to $\boldsymbol{\theta}$, $Da(z,\boldsymbol{\theta})=\displaystyle\int_{0}^{{\infty}}a(z+t\boldsymbol{\theta})dt$ is the divergent beam transform of the attenuation $a$, $Ra(s,\boldsymbol{\theta}^{\perp})=\displaystyle\int_{-{\infty}}^{{\infty}}a\left(s\boldsymbol{\theta}^{\perp}+t\boldsymbol{\theta}\right)dt$ is the Radon transform of the attenuation $a$, and the classical Hilbert transform $Hh(s)=\displaystyle\frac{1}{\pi}\int_{-{\infty}}^{{\infty}}\frac{h(t)}{s-t}dt$ is taken in the first variable and evaluated at $s=z\cdotp\boldsymbol{\theta}^{\perp}$. The function $h$ appeared first in the work of Natterer [32]; see also [8] for elegant arguments that show how $h$ extends analytically (in the angular variable on the unit circle ${{\mathbf{S}}^{1}}$) inside the unit disc. We recall some properties of $h$ from [40, Lemma 4.1], while establishing notations.

Lemma 2.1**.**

[40*, Lemma 4.1]** Assume $\Omega\subset\mathbb{R}^{2}$ is $C^{2,s},s>1/2$, convex domain. For $p=1,2$, let $a\in C^{p,s}(\overline{\Omega})$, $s>1/2$, and $h$ defined in (11). Then $h\in C^{p,s}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$ and the following hold

(i) $h$ satisfies*

[TABLE]

(ii) $h$ has vanishing negative Fourier modes yielding the expansions

[TABLE]

with (iii)

[TABLE]

(iv) For any $z\in\Omega$

[TABLE]

(v) For any $z\in\Omega$

[TABLE]

(vi) The Fourier modes $\alpha_{k},\beta_{k},k\geq 0$ satisfy

[TABLE]

The Fourier coefficients of $e^{\pm h}$ define the integrating operators $e^{\pm G}{\bf u}$ component-wise for each $m\leq 0$ by

[TABLE]

where $\alpha_{k}$ and $\beta_{k}$ are the Fourier modes of $e^{-h}$ and $e^{h}$ in (13), and $\boldsymbol{\alpha},\boldsymbol{\beta}$ as in (14), respectively, (15). Note that $e^{\pm G}$ can also be written in terms of left translation operator as

[TABLE]

where $L^{k}=\underbrace{L\circ\cdots\circ L}_{k}$ is the $k$-th composition of left translation operator. It is important to note that the operators $e^{\pm G}$ commute with the left translation, $[e^{\pm G},L]=0$.

Different from [40], in this work we carry out the analysis in the Sobolev spaces $l^{2,p}(\mathbb{N};H^{q}(\Omega))$ with the respective norm $\left\|\cdot\right\|_{p,q}$ given by

[TABLE]

In Proposition 2.1 below we revisit the mapping properties of $e^{\pm G}$ relative to these new spaces.

Throughout, in the notation of the norms, the first index $p\in\{0,\frac{1}{2},1\}$ refers to the smoothness in the angular variable (expressed as decay in the Fourier coefficient), while the second index $q\in\{0,1\}$ shows the smoothness in the spatial variable. The most often occurring is the space $l^{2}(\mathbb{N};L^{2}(\Omega))$, when $p=q=0$. To simplify notation, in this case we drop the double zero subindexes,

[TABLE]

The traces on the boundary $\varGamma$ of functions in $l^{2,p}(\mathbb{N};H^{1}(\Omega))$ are in $l^{2,p}(\mathbb{N};H^{\frac{1}{2}}(\varGamma))$, endowed with the norm

[TABLE]

Since $\Gamma$ is the unit circle, the $H^{1/2}(\varGamma)$-norm can be defined in the Fourier domain as follows. For each integer $j\geq 0$, if we consider the Fourier expansion of the trace $\displaystyle u_{-j}\lvert_{\varGamma}$,

[TABLE]

then

[TABLE]

In view of (25), if ${\bf g}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{\frac{1}{2}}(\varGamma))$, then

[TABLE]

In the estimates we need the following variant of the Poincaré inequality obtained by component-wise summation: If ${\bf u}\in l^{2}(\mathbb{N};H^{1}(\Omega))$, then

[TABLE]

where $\mu$ is a constant depending only on $\Omega$; for the unit disc $\mu=2$.

Consider the Banach space

[TABLE]

Proposition 2.1**.**

Let $a\in C^{2,s}(\overline{\Omega})$, $s>1/2$. Then $\boldsymbol{\alpha},\partial\boldsymbol{\alpha},\boldsymbol{\beta},\partial\boldsymbol{\beta}\in l^{1,1}_{{\infty}}(\overline{\Omega})$, and for any $p\in\{0,\frac{1}{2},1\}$, $q\in\{0,1\}$, the operators

[TABLE]

are bounded, and satisfy the following estimates

[TABLE]

The same estimate works for $e^{G}{\bf u}$ with $\boldsymbol{\alpha}$ replaced by $\boldsymbol{\beta}$.

The proof of the Proposition 2.1 can be found in the Appendix.

We remark that cases $(p=0=q)$, $(p=1,q=0)$, and $(p=0,q=1)$ in Proposition 2.1 hold for $a\in C^{1,s}(\overline{\Omega})$, $s>1/2$, and these are the only properties needed for the lemma below.

Lemma 2.2**.**

Let $a\in C^{1,s}(\overline{\Omega})$, $s>1/2$, and $e^{\pm G}$ as defined in (21).

(i) If ${\bf u}\in l^{2}(\mathbb{N};H^{1}(\Omega))$ solves $\displaystyle\overline{\partial}{\bf u}+L^{2}\partial{\bf u}+aL{\bf u}=0$, then $\displaystyle{\bf v}=e^{-G}{\bf u}\in l^{2}(\mathbb{N};H^{1}(\Omega))$ solves $\overline{\partial}{\bf v}+L^{2}\partial{\bf v}=0$.

(ii) Conversely, if ${\bf v}\in l^{2}(\mathbb{N};H^{1}(\Omega))$ solves $\overline{\partial}{\bf v}+L^{2}\partial{\bf v}=0$, then $\displaystyle{\bf u}=e^{G}{\bf v}\in l^{2}(\mathbb{N};H^{1}(\Omega))$ solves $\displaystyle\overline{\partial}{\bf u}+L^{2}\partial{\bf u}+aL{\bf u}=0$.

Proof.

(i) Let $\displaystyle{\bf v}=e^{-G}{\bf u}=\sum_{k=0}^{\infty}\alpha_{k}L^{k}{\bf u}$. Since ${\bf u}\in l^{2}(\mathbb{N};H^{1}(\Omega))$, then from Proposition 2.1, $\displaystyle{\bf v}\in l^{2}(\mathbb{N};H^{1}(\Omega))$. Then ${\bf v}$ solves

[TABLE]

where in the last equality we have used (18) and (19).

An analogue calculation using the properties in Lemma 2.1 (iv) shows the converse.

∎

3. Source reconstruction for scattering of polynomial type

This section contains the basic idea of reconstruction in the special case of scattering kernel of polynomial type,

[TABLE]

for some fixed integer $M\geq 1$. Recall that since $k(z,\cos\theta)$ is both real valued and even in $\theta$, $k_{n}(z)$ are real valued and $k_{n}(z)=k_{-n}(z)$, $0\leq n\leq M$. We stress here that no smallness assumption on $k_{0},k_{-1},k_{-2},\cdots,k_{-M}$ is assumed. Let $u^{(M)}$ be the solution of (1) with $k$ as in (34) and ${\bf u}^{(M)}$ denote the sequence valued map

[TABLE]

Let also $K^{(M)}$ denote the corresponding Fourier multiplier operator

[TABLE]

The transport equation (1) reduces to the system

[TABLE]

In sequence valued notation, the system (38) and (39) rewrites:

[TABLE]

where $K^{(M)}$ as in (36).

Since $f\in H^{1}(\Omega)$, the solution $u\in H^{1}(\Omega\times{{\mathbf{S}}^{1}})$, and consequently ${\bf u}^{(M)}\in l^{2,1}(\mathbb{N};H^{1}(\Omega))$. We note that in our method we only use ${\bf u}^{(M)}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$, indicating that it may apply to rougher sources.

Let the transformation ${\bf v}^{(M)}=e^{-G}L^{M}{\bf u}^{(M)}$, then by Proposition 2.1, ${\bf v}^{(M)}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$, and ${\bf v}^{(M)}$ is $L^{2}$-analytic:

[TABLE]

The trace of the boundary ${\bf v}^{(M)}\lvert_{\varGamma}$ is determined by the trace of ${\bf u}^{(M)}\lvert_{\varGamma}={\bf g}=(g_{0},g_{-1},g_{-2},...)\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1/2}(\varGamma))$, by

[TABLE]

By Proposition 2.1, ${\bf v}^{(M)}\lvert_{\varGamma}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1/2}(\varGamma))$.

The Bukhgeim-Cauchy integral formula (10) extends ${\bf v}^{(M)}$ from $\varGamma$ to $\Omega$ as $L^{2}$-analytic map. From the uniqueness of an $L^{2}$-analytic map with a given trace, we recovered for $n\geq 0,$

[TABLE]

Thus ${\bf v}^{(M)}=\langle v^{(M)}_{0},v^{(M)}_{-1},v^{(M)}_{-2},...\rangle$ is recovered in $l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$.

We recover $L^{M}{\bf u}^{(M)}=\langle u^{(M)}_{-M},u^{(M)}_{-M-1},u^{(M)}_{-M-2},...\rangle$ in $\Omega$ by using the convolution formula (21)

[TABLE]

where $\beta_{k}$’s as in (13). In particular we recovered $u^{(M)}_{-M-1},u^{(M)}_{-M}\in H^{1}(\Omega)$.

By applying $4\partial$ to (38), the mode $u^{(M)}_{-M+1}$ is then the solution to the Dirichlet problem for the Poisson equation

[TABLE]

Since by construction $u^{(M)}_{-M},u^{(M)}_{-M-1}\in H^{1}(\Omega)$, we have

[TABLE]

Since $g_{-M+1}\in H^{1/2}(\varGamma)$, the solution $u^{(M)}_{-M+1}\in H^{1}(\Omega)$ and

[TABLE]

where the constant $C$ depends only on $\Omega$ and $\displaystyle\max\left\{1,\max_{0\leq j\leq M}\left\|a-k_{-j}\right\|_{L^{\infty}(\Omega)}^{2}\right\}$. Successively all the other modes $u^{(M)}_{-M+j}$ for $j=2,\cdots,M$ are computed by solving the corresponding Dirichlet problem for the Poisson equation. To account for the successive accumulation of error we note the following result which can be proven by induction.

Lemma 3.1**.**

Let $\{a_{n}\}$ and $\{b_{n}\}$ be sequences of nonnegative numbers, such that

[TABLE]

where $c>0$ is a constant, then

[TABLE]

We applying Lemma 3.1 to (46) and estimate

[TABLE]

and

[TABLE]

The source $f^{(M)}$ is computed by

[TABLE]

and we estimate

[TABLE]

This method is implemented in the numerical experiments in Section 7. Next we analyze the error introduced by truncation.

4. Gradient estimates of solutions to nonhomogeneous Bukhgeim-Beltrami equation

When applying the reconstruction method to the data arising from a general scattering kernel $\displaystyle k(z,\cos\theta)=\sum_{n=0}^{{\infty}}k_{-n}(z)\cos(n\theta)$ an error is made due to the truncation in the Fourier series of $k$. This error is controlled by the gradient of the solution to the Cauchy problem for the inhomogeneous Bukhgeim-Beltrami equation

[TABLE]

for some specific ${\bf f}$ and operator coefficient $B$. Estimates of the gradient for solutions of (51) may be of separate interest, reason for which we treat them here independently of the transport problem.

We start with an energy identity (see [49] for ${\bf f}=0$), à la Mukhometov [30] or Pestov [42].

Theorem 4.1**.**

*(Energy identity) Let ${\bf f}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$ and let $B$ be a bounded operator such that $B:l^{2}(\mathbb{N};L^{2}(\Omega))\to l^{2,1}(\mathbb{N};L^{2}(\Omega))$ and $B:l^{2}(\mathbb{N};H^{1}(\Omega))\to l^{2,1}(\mathbb{N};H^{1}(\Omega))$.

If ${\bf v}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$ is a solution to (51), then*

[TABLE]

Proof.

Using the Green’s identity $\displaystyle 2\int_{\Omega}\left\|\partial{\bf v}\right\|_{l_{2}}^{2}dx=2\int_{\Omega}\left\|\overline{\partial}{\bf v}\right\|_{l_{2}}^{2}dx+i\int_{\varGamma}\langle{\bf v},\partial_{s}{\bf v}\rangle ds,$ where $\partial_{s}$ is the tangential derivative at the boundary, it follows that

[TABLE]

For each $n\in\mathbb{N}$,

[TABLE]

By summing in $n$, and using $\displaystyle\lim_{n\to\infty}\int_{\Omega}\left\|L^{2n}\partial{\bf v}\right\|_{l_{2}}^{2}dx=0$, we conclude the theorem. ∎

We note the general identity [38, Lemma 2.1], for a sequence of nonnegative numbers:

Lemma 4.1**.**

Let $\{c_{n}\}$ be a sequence of nonnegative numbers. Then

[TABLE]

whenever one of the sides in (i) and (ii) is finite.

Proof.

(i) By changing the index $j=m+n$, for $m\geq 0$, ($j-n\geq 0,$ and $n\leq j$), we get

[TABLE]

(ii) Similarly, by changing the index $j=m+2n$, for $m\geq 0$, $\left(j-2n\geq 0,\;\text{and}\;\displaystyle n\leq\left\lfloor\frac{j}{2}\right\rfloor\right)$, we get

[TABLE]

∎

Theorem 4.2** (Gradient estimate).**

Let ${\bf f}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$, and $B$ be a smoothing operator such that $B:l^{2}(\mathbb{N};L^{2}(\Omega))\to l^{2,1}(\mathbb{N};L^{2}(\Omega))$ and $B:l^{2}(\mathbb{N};H^{1}(\Omega))\to l^{2,1}(\mathbb{N};H^{1}(\Omega))$ are bounded, and let $C_{B}>0$ be such that

[TABLE]

Assume that $B$ is such that

[TABLE]

where $\mu$ is the factor in Poincaré inequality (27).

If ${\bf v}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$ is a solution to the inhomogeneous Bukgheim-Beltrami equation (51), then

[TABLE]

where

[TABLE]

Proof.

We estimate each term on the right hand side of the energy identity (52). For brevity if we denote $B{\bf v}=\langle b_{0},b_{-1},b_{-2},\cdots\rangle$, then $\displaystyle\left\|B{\bf v}\right\|_{1,0}^{2}=\sum_{j=0}^{\infty}(1+j^{2})\int_{\Omega}|b_{-j}|^{2}$.

1. I.

We estimate the first term in (52):

[TABLE]

where in the first inequality we use Lemma 4.1 part (i), in the second inequality we use Cauchy-Schwarz inequality, in the third inequality we use $(1+x)^{2}\leq 2(1+x^{2})$, in the fifth inequality we use $\left\|B{\bf v}\right\|_{1,0}\leq C_{B}\left\|{\bf v}\right\|$, and in the next to last inequality we use the Poincaré inequality (27). 2. II.

We estimate the second term in (52):

[TABLE]

where in the first inequality we use Lemma 4.1 part (i), and then Cauchy-Schwarz. 3. III.

We estimate the third term in (52):

[TABLE]

where in the first inequality we use Lemma 4.1 part (i), in the next to the last inequality we use $\left\|B{\bf v}\right\|_{1,0}\leq C_{B}\left\|{\bf v}\right\|$, while in the last we use the Poincaré inequality (27). 4. IV.

We estimate the fourth term in (52):

[TABLE]

where in the first inequality we have used Lemma 4.1 part (i). 5. V.

We estimate the next to the last term in (52):

[TABLE]

where in the first inequality we have used Cauchy-Schwarz, in the second inequality we have used the fact $xy\leq\frac{1}{2}(x^{2}+y^{2})$, and in the last inequality we have used estimates from (III) and (IV). 6. VI.

To estimate the last term in (52), we use the fact that $\varGamma$ is the unit circle and consider the Fourier expansion of the traces of the modes $v_{-j}\lvert_{\varGamma}$. For $e^{i\beta}\in\varGamma$,

[TABLE]

Using the parametrization, the last term in (52) becomes

[TABLE]

where in the second equality we have used Lemma 4.1 part (ii), in the first inequality we have used the fact $\left(1+\left\lfloor\frac{j}{2}\right\rfloor\right)\leq(1+j^{2})^{1/2}$, and in the last equality we have used the definition of the norm (26).

Using the above estimates (I)-(VI) for the expressions in (52), we have proved for $\displaystyle\tau=\left(\int_{\Omega}\left\|\partial{\bf v}\right\|_{l_{2}}^{2}\right)^{1/2}$, that $a\tau^{2}-b\tau-c\leq 0$. Assumption on $\epsilon$ as in (54) yield $a>0$ and we have the estimate (55). ∎

For the case when $B=0$ we obtain the immediate corollary.

Corollary 4.1**.**

Let ${\bf f}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$. If ${\bf v}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$ solves

[TABLE]

then

[TABLE]

Proof.

This is the case $\epsilon=0$ and $a=1$ in (55). We also use $\displaystyle\left(\frac{b+\sqrt{b^{2}+4c}}{2}\right)^{2}\leq b^{2}+2c.$ ∎

5. Smoothing due to scattering

In this section we explicit the smoothing properties of the Fourier multiplier operator $K$ in (7) as determined by the appropriate decay of the Fourier coefficients of the scattering kernel $k_{n}(z)=\frac{1}{2\pi}\int_{-\pi}^{\pi}k(z,\cos\theta)e^{-in\theta}d\theta$, $n\in\mathbb{Z}$. The gain of $1/2$ smoothness in the angular variable (see (61) below) has been shown before by different methods in [26], and in a more general case than considered here in [46].

Lemma 5.1**.**

(Smoothing due to scattering) Let $M\geq 1$ be a positive integer and $K$ be the Fourier multiplier in (7). Assume that $k$ is such that its Fourier coefficients starting from index $M$ onward satisfy

[TABLE]

(i) If ${\bf v}\in l^{2}(\mathbb{N};L^{2}(\Omega))$, then

[TABLE]

(ii) $K:l^{2}(\mathbb{N};H^{1}(\Omega))\to l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$ is bounded. More precisely,

[TABLE]

*(iii) Moreover, if (59) holds for $p\geq 1$, then $\displaystyle K:l^{2}(\mathbb{N};L^{2}(\Omega))\to l^{2,1}(\mathbb{N};L^{2}(\Omega))$ and

$\displaystyle K:l^{2}(\mathbb{N};H^{1}(\Omega))\to l^{2,1}(\mathbb{N};H^{1}(\Omega))$ are bounded, and*

[TABLE]

Proof.

(i) Let ${\bf v}\in l^{2}(\mathbb{N};L^{2}(\Omega))$. Then

[TABLE]

where in the second equality we have used Lemma 4.1 part (i), and in the first inequality we have used (59).

To prove (ii), let ${\bf v}\in l^{2}(\mathbb{N};H^{1}(\Omega))$. Then

[TABLE]

We estimate the first term in (64),

[TABLE]

where in the first inequality we have used decay of $k_{j}$’s.

Similarly, following the same proof as of the first term (65), the term

[TABLE]

and the last term

[TABLE]

Thus the expression in (64) becomes

[TABLE]

To prove (iii), assume that $k$ is such that (59) holds for $p\geq 1$. Let ${\bf v}\in l^{2}(\mathbb{N};L^{2}(\Omega))$, then

[TABLE]

where in the first inequality we have used decay of $k_{j}$’s.

To prove the last part, let ${\bf v}\in l^{2}(\mathbb{N};H^{1}(\Omega))$. Then

[TABLE]

Following the similar proof of (ii), $\left(\text{with}\;(1+j^{2})\;\text{instead of}\;(1+j^{2})^{1/2}\right)$, we have

[TABLE]

∎

6. Error estimates

Recall the sequence valued maps $\displaystyle{\bf u}=\langle u_{0},u_{-1},u_{-2},...\rangle$ and $\displaystyle{\bf u}^{(M)}=\langle u^{(M)}_{0},u^{(M)}_{-1},u^{(M)}_{-2},...\rangle$ solutions of $\displaystyle\overline{\partial}{\bf u}+L^{2}\partial{\bf u}+aL{\bf u}=LK{\bf u}$, respectively of $\overline{\partial}{\bf u}^{(M)}+L^{2}\partial{\bf u}^{(M)}+aL{\bf u}^{(M)}=LK^{(M)}{\bf u}^{(M)}$, where $K$ is the multiplier operator in (7), and $K^{(M)}$ is its truncated version in (36) corresponding to the $M$-th order polynomial approximation of $k$.

The error we make in the source reconstruction is controlled by the sequence valued map

[TABLE]

which solves

[TABLE]

Recall that the integrating operators $e^{\pm G}$ commute with the left translation, $[e^{\pm G},L]=0$. From the equation (36) is easy to see that the translated sequence

[TABLE]

solves

[TABLE]

For $a\in C^{2,s}(\overline{\Omega})$, $s>1/2$, let $\boldsymbol{\alpha},\boldsymbol{\beta}$ be as in Proposition 2.1. For brevity, let us define

[TABLE]

Under sufficient regularity on $k$, its Fourier coefficients satisfy (59) with $p=1$.

Lemma 6.1**.**

If $k\in C^{1,s}_{per}([-1,1];Lip(\Omega))$, $s>1/2$, then

[TABLE]

Proof.

The proof follows directly by Bernstein’s lemma [19, Chap. I, Theorem 6.3]. ∎

As a consequence we obtain the following a priori error estimate for the source reconstructed by the method in Section 3.

Theorem 6.1**.**

(Error estimate) Let $M\geq 1$ be a positive integer, $a\in C^{2,s}(\overline{\Omega})$, $k\in C^{1,s}_{per}([-1,1];Lip(\Omega))$, for $s>1/2$, be known, and $f\in H^{1}(\Omega)$ and $u\in H^{1}(\Omega\times{{\mathbf{S}}^{1}})$ be unknown functions satisfying (1). Let ${\bf u}$ be the unknown sequence of nonpositive Fourier modes of $u$ as in (3), and $f^{(M)}$ in (49) be the reconstructed source using the noisy boundary data $\widetilde{{\bf g}}\lvert_{\varGamma}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1/2}(\varGamma))$. Let $\delta{\bf g}={\bf u}\lvert_{\varGamma}-\widetilde{{\bf g}}$ be the error in the data. The error $\delta f=f-f^{(M)}$ in the reconstructed source is estimated by

[TABLE]

where $c_{1}=96(1+\mu)\sigma^{4}\gamma_{M}^{2}$, $c_{2}=4\pi(1+\mu)\sigma^{4}$, $c_{3}=\max\left\{2\mu\sigma^{4},2\right\}$, with $\sigma$ in (70), $\gamma_{M}$ in (71), and $C$ depends only on $\Omega$ and $\displaystyle\max\left\{1,\max_{0\leq j\leq M}\|a-k_{-j}\|_{L^{\infty}(\Omega)}^{2}\right\}$.

Proof.

Since $u\in H^{1}(\Omega\times{{\mathbf{S}}^{1}})$, ${\bf u}\in l^{2,1}(\mathbb{N};H^{1}(\Omega))$, in particular ${\bf u}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$.

From linearity of the problem, the error $\delta f$ also satisfy (50):

[TABLE]

where $\delta g_{m}$ for $-M+1\leq m\leq 0$, is the $m$-th Fourier coefficient of $\delta{\bf g}$.

Recall that the $\delta u_{-M}$ and $\delta u_{-M-1}$ are the first two components of $e^{G}L^{M}{\bf q}^{(M)}$. Using the gradient estimate in (58) we have

[TABLE]

where the first inequality uses (31) in Proposition 2.1, and the second inequality uses (62) with $p=1$. By Poincare and using (74),

[TABLE]

From (75) and (73), the expression in (73) becomes

[TABLE]

where constant $c_{1}=96(1+\mu)\sigma^{4}\gamma_{M}^{2}$, $c_{2}=4\pi(1+\mu)\sigma^{4}$, $c_{3}=\max\left\{2\mu\sigma^{4},2\right\}$, and $C$ depends only on $\Omega$ and $\displaystyle\max\left\{1,\max_{0\leq j\leq M}\left\|a-k_{-j}\right\|_{L^{\infty}(\Omega)}^{2}\right\}$. Thus (72) follows. ∎

For the case when $k$ is a polynomial in the angular variable, $k_{-j}=0$ for $j\geq M+1$, we obtain the following stability result.

Corollary 6.1**.**

(Stability) Assume the hypotheses in Theorem 6.1 with $k$ being a polynomial of degree $M$, i.e., $k_{-j}=0$ for $j\geq M+1$. Then we have the following stability estimate

[TABLE]

where $c_{2}=4\pi(1+\mu)\sigma^{4}$, $c_{3}=\max\left\{2\mu\sigma^{4},2\right\}$, with $\sigma$ in (70), $\gamma_{M}$ in (71), and $C$ depends only on $\Omega$ and $\displaystyle\max\left\{1,\max_{0\leq j\leq M}\left\|a-k_{-j}\right\|_{L^{\infty}(\Omega)}^{2}\right\}$.

Proof.

This is the case $\gamma_{M}=0$ in (72). So $c_{1}=0$ and result (76) follows. ∎

We stress that if $k$ is not of polynomial type, then $(1+C)^{M+1}\left\|L^{M+1}{\bf u}\right\|^{2}$ can not be made arbitrarily small. However, if the anisotropic part of scattering is small enough (same case as in [7]), then we can prove a convergence result.

We have the following convergence result for the weakly anistropic scattering media.

Theorem 6.2**.**

Let $a\in C^{2,s}(\overline{\Omega})$, $k\in C^{1,s}_{per}([-1,1];Lip(\Omega))$, for $s>1/2$, and $f\in H^{1}(\Omega)$. Assume that

[TABLE]

*where $\mu$ as in Poincaré inequality (27) and $\sigma$ as in (70).

For each arbitrarily fixed $M$, let $f^{(M)}$ in (49) be the reconstructed source using the data*

[TABLE]

where $g_{-j}$, for $0\leq j\leq M$, are assumed exact data and $\widetilde{g_{-j}}$, for $j>M$, are allowed to be noisy. Then

[TABLE]

Proof.

Since $f\in H^{1}(\Omega)$, the solution $u\in H^{1}(\Omega\times{{\mathbf{S}}^{1}})$, and consequently ${\bf u}\in l^{2,1}(\mathbb{N};H^{1}(\Omega))$, in particular ${\bf u}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$.

The error we make in the source reconstruction is controlled by the sequence valued map ${\bf q}^{M}=e^{-G}({\bf u}-{\bf u}^{M})\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$, which solves

[TABLE]

where $B^{M}$ and $A^{M}$ are operators given by

[TABLE]

with $e^{\pm G}$ as in (21), $K$ as in (7) and $K^{(M)}$ is the truncated Fourier multiplier as in (36).

The trace

[TABLE]

is determined by the error in the data

[TABLE]

For $c_{1}=\max\left\{\left\|(a-k_{0})\right\|^{2}_{{\infty}},2\right\}$, the norm

[TABLE]

where in the fourth and last inequality we use Proposition 2.1, in the next to last inequality we use the Poincaré inequality (27) and in the last inequality we have used (82).

To estimate the gradient term $\left\|\partial{\bf q}^{(M)}\right\|^{2}$ in (84), we employ Theorem 4.1 with ${\bf v}={\bf q}^{M}$, $B=B^{M}$, and ${\bf f}=A^{M}{\bf u}$ in there as follows. For ${\bf u}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$, we have

[TABLE]

where the first inequality uses Proposition 2.1, and the last inequality uses Lemma 5.1 (iii) for $p=1$. Since ${\bf q}^{M}\in l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))$, we have

[TABLE]

where the first inequality uses Proposition 2.1, and the last inequality uses Lemma 5.1 (iii) for $p=1$. We can apply Theorem 4.1 with $C_{B}=\sqrt{2}\gamma\sigma^{2}$ and $\epsilon=\sqrt{2\mu}C_{B}<\sqrt{2}-1$ by (77), to obtain

[TABLE]

where

[TABLE]

Since the data is assumed exact for the first $M$ modes, $\displaystyle({\bf g}-\widetilde{{\bf g}})\lvert_{\varGamma}$ in (83) has the first $M$ modes equal to zero. Thus

[TABLE]

as $M\to{\infty}$. Since ${\bf u}$ is fixed, the expression $\left\|L^{M+1}{\bf u}\right\|\to 0$, as $M\to{\infty}$. Therefore, both $b_{M}$ and $c_{M}$ tends to zero and $\displaystyle\lim_{M\to{\infty}}\left\|\partial{\bf q}^{(M)}\right\|\to 0.$ Using (84), the norm $\displaystyle\lim_{M\to{\infty}}\left\|f-f^{(M)}\right\|^{2}_{L^{2}(\Omega)}\to 0$. ∎

7. Numerical Experiments

In this section we demonstrate the numerical feasibility of the method on two numerical examples. We focus on the numerical results, and leave their realization details for a separate discussion  [17]. Although the theoretical results assume a source of square integrable gradient, we shall see below that the numerical reconstruction works even in the case of a discontinuous source.

The numerical experiments consider

[TABLE]

with the two dimensional Henyey-Greenstein (Poisson) model of scattering kernel

[TABLE]

and the attenuation coefficient

[TABLE]

where $\mu_{\text{s}}$ is the scattering coefficient, and $\mu_{\text{a}}$ is the absorption coefficient. The parameter $g$ in (85) models a degree of anisotropy, with $g=1$ for the ballistic regime and $g=0$ for the isotropic scattering regime. In our numerical experiments we work with $g=1/2$ and $\mu_{\text{s}}(z)=5$, the latter value shows that, on average, the particle scatters within $1/5$ units of length along the path.

Let $R=(-0.25,0.5)\times(-0.15,0.15)$ be the rectangular region, and

[TABLE]

be the circular region inside the unit disc $\Omega$ as shown in Figure 1.

We consider two examples for the attenuation coefficient $a(z)$ for two different functions $\mu_{\text{a}}(z)$, while $\mu_{s}=5$ is the same. In the first example, we work with a $C^{2}$ smooth absorption, whereas in the second example we consider attenuation to be discontinuous. In the first example, the attenuation $a$ is $C^{2}$-smooth via some scaled and translated quartic polynomial $(|z|+1)^{2}(|z|-1)^{2}$, $0\leq|z|\leq 1$, in the $\epsilon$-neighborhoods of $\partial B_{1}$ and $\partial B_{2}$:

[TABLE]

with $\epsilon=0.025$.

In both examples the data is generated with the same source

[TABLE]

This is the function we reconstruct by implementing the method in Section 3.

We generate the boundary measurement by the numerical computation of the forward problem by piecewise constant approximation following the method in [16]. In the forward problem the triangular mesh has $1,966,690$ triangles, while the velocity direction is split into $360$ equi-spaced intervals.

To obtain the data, we disregard the value of the solution inside and only keep the boundary values. The boundary data, $u(z,\boldsymbol{\theta})$ on $\partial\Omega\times S^{1}$, is represented in Figure 2:

For each $z\in\partial\Omega$ (indicated by a cross), we represent the graph \{\bigl{(}u(z,\boldsymbol{\theta}),\boldsymbol{\theta}\bigr{)}:\;\boldsymbol{\theta}\in S^{1}\} as the (red) closed curve in the polar coordinates \big{\{}z+2u(z,\boldsymbol{\theta})\boldsymbol{\theta}\>:\>\boldsymbol{\theta}\in S^{1}\bigr{\}}. Note that, since $u|_{\Gamma_{-}}=0$ (there is no incoming radiation), the red curve lies outside the domain $\Omega$. The nearly tangential behavior at the boundary is the numerical evidence that the regime is far from ballistic, while the irregular shapes show that it is far from isotropic.

The reconstruction starts with solving (39) with $M=8$ via the Bukhgeim-Cauchy integral formula (10), where the infinite series is replaced by a finite sums of up to $64$ terms.

In the reconstruction procedure, all the steps, including the integral transforms $D[a]$, $R[a]$ and H\bigl{[}R[a]\bigr{]}, have been processed numerically [17]. In solving the Poisson equation (45a) and (45b), the standard $P_{1}$ finite element method (constant and piecewise linear elements) is employed.

The triangulation used in the reconstruction is different from that in the forward problem. In particular the reconstruction mesh consists of $6,998$ triangles (much less than the $1,966,690$ triangles used in the forward problem), and is generated without any information of the location of the subsets $R$, $B_{1}$, and $B_{2}$.

In the first numerical experiment the attenuation is given by (86). The reconstructed source $f(z)$ is shown in Figure 3 on the left. On the right we show the cross section of $f(z)$ along the dotted diameter $y=-\sqrt{3}x$ passing through the origin and the center of $B_{2}$.

The reconstructed $f(z)$ shows a quantitative agreement with the exact source in (87). Similar to the attenuated X-ray tomography, the artifacts appear due to the co-normal singularities in the source, but also in the attenuation. This is because the most singular is the ballistic part (attenuated X-ray transform of the source).

In the second example, we apply the proposed algorithm to the case of a discontinuous attenuation $a=\mu_{\text{s}}+\mu_{\text{a}}$ with

[TABLE]

While different, the boundary data in the second experiment is graphically indistinguishable from the data in Figure 2 in the smooth case.

The numerically reconstructed source, shown in Figure 4, agrees well with the exact source (87), implying that the proposed algorithm has a potential in reconstruction even if the attenuation admits the discontinuity.

Figure 5 shows the projection on the $x$-axis of the reconstructed $f(z)$ for $z\in[-0.4,0.6]\times[-0.1,0.1]$. The arrows show the location of $\partial B_{1}$ in $R$, where the attenuation has a jump.

In both reconstructions two types of artifacts appear, one type due to the singular support of $f$ and the other due to the singular support of $a$. Observe how the effect of the latter type of singularities diminishes with an increase in the smoothness of $a$.

At the level of singular support, the image is similar to that in the non-scattering case ($k=0$). This is expected due to the smoothing effect of scattering. A quantitative understanding of the effect of singularities in the attenuation, in the presence of scattering is subject to future work.

Acknowledgment

The authors thank D. Kawagoe for useful discussions of regularity of solution of the forward problem. The work of H.  Fujiwara was supported by JSPS KAKENHI Grant Numbers 16H02155, 18K18719, and 18K03436. The work of K.  Sadiq was supported by the Austrian Science Fund (FWF), Project P31053-N32.

Appendix: New mapping properties of the integrating operator $e^{\pm G}$

In this section we prove Proposition 2.1.

Proof.

Since $a\in C^{2,s}(\overline{\Omega})$, $s>1/2$, it follows from Lemma 2.1 that $e^{\pm h}\in C^{2,s}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$. For $z\in\overline{\Omega}$, let $\boldsymbol{\alpha}(z)=\langle\alpha_{0}(z),\alpha_{1}(z),...\rangle$, and $\partial\boldsymbol{\alpha}(z)=\langle\partial\alpha_{0}(z),\partial\alpha_{1}(z),...\rangle$, where $\partial$ is the derivative in the spatial variable. By the proof in Bernstein’s lemma [19, Chap. I, Theorem 6.3],

[TABLE]

Similarly, $\boldsymbol{\beta},\partial\boldsymbol{\beta}\in l^{1,1}_{{\infty}}(\overline{\Omega})$.

To prove the case $p=0=q$, let ${\bf u}\in l^{2}(\mathbb{N};L^{2}(\Omega))$. By using Young’s inequality

[TABLE]

where the discrete convolution (with respect to the index) is defined in (21).

To prove the case $p=1$ and $q=0$, let ${\bf u}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$. We estimate the norm

[TABLE]

where in the first inequality we use estimate (88) from the case $p=0=q$.

To estimate the last term in (89), we need to account for both positive and negative Fourier modes. Let $\widehat{u}$, and $\widehat{e^{-h}}$ denote the sequence valued of the Fourier modes of $u$, and $e^{-h}$, respectively, i.e.,

[TABLE]

By using Plancherel we estimate

[TABLE]

To estimate the last two terms in (91), let $\widehat{\partial_{\theta}u}$ denote the (double) sequence of the Fourier modes of $\partial_{\theta}u$, and let $\widehat{\partial_{\theta}e^{-h}}$ denote the sequence valued of the Fourier modes of $\partial_{\theta}e^{-h}$, i.e.,

[TABLE]

By using Plancherel the expression in (91) becomes

[TABLE]

where in the second inequality we use Young’s inequality, in the first equality we use $\displaystyle\sup_{\overline{\Omega}}\left\|\widehat{e^{-h}}\right\|^{2}_{l_{1}(\mathbb{Z})}=\sup_{\overline{\Omega}}\left\|\boldsymbol{\alpha}\right\|^{2}_{l_{1}}$, and $\displaystyle\sup_{\overline{\Omega}}\left\|\widehat{\partial_{\theta}e^{-h}}\right\|^{2}_{l_{1}(\mathbb{Z})}=\sup_{z\in\overline{\Omega}}\left(\sum_{k=1}^{{\infty}}k\lvert\alpha_{k}(z)\rvert\right)^{2}$, and in the last equality we use the norm in (28) and $u_{-n}=\overline{u_{n}}$ as $u$ is real valued.

From the above estimate of (91), the expression in (89) becomes

[TABLE]

To prove the case $p=0$ and $q=1$, let ${\bf u}\in l^{2}(\mathbb{N};H^{1}(\Omega))$, then ${\bf u},\partial{\bf u}\in l^{2}(\mathbb{N};L^{2}(\Omega))$. For ${\bf u}\in l^{2}(\mathbb{N};L^{2}(\Omega))$, $e^{\pm G}{\bf u}\in l^{2}(\mathbb{N};L^{2}(\Omega))$ from the case $p=0=q$, so it suffices to show that $\partial(e^{-G}{\bf u})\in l^{2}(\mathbb{N};L^{2}(\Omega))$.

The discrete convolution below is with respect to the index, while $\partial$ acts on the spatial variable.

[TABLE]

where in the second inequality we use discrete Young’s inequality.

From (94) and estimate (88) from the case $p=0=q$, we estimate the norm

[TABLE]

To prove the case $p=1=q$, let ${\bf u}\in l^{2,1}(\mathbb{N};H^{1}(\Omega))$, then ${\bf u},\partial{\bf u}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$. Since ${\bf u}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$, $e^{-G}{\bf u}\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$ by case $p=1$ and $q=0$. Thus suffices to show that $\partial(e^{-G}{\bf u})\in l^{2,1}(\mathbb{N};L^{2}(\Omega))$. We estimate the norm

[TABLE]

where in the first inequality we have used estimate (94) from the case $p=0$ and $q=1$.

To estimate the last term in (96), we need to account for both positive and negative Fourier modes. Let $\widehat{u}$, and $\widehat{e^{-h}}$ denote the sequence valued of the Fourier modes of $u$, and $e^{-h}$, respectively, as in (90), and let $\widehat{\partial_{\theta}u}$ and $\widehat{\partial_{\theta}e^{-h}}$ denote the sequence valued of the Fourier modes of $\partial_{\theta}u$, and $\partial_{\theta}e^{-h}$, respectively, as in (92). Moreover, let $\widehat{\partial u}$, and $\widehat{\partial_{\theta}\partial u}$ denote the sequence valued of the Fourier modes of $\partial u$, respectively, $\partial_{\theta}\partial u$, i.e.

[TABLE]

Similarly, let $\widehat{\partial e^{-h}}$, and $\widehat{\partial_{\theta}\partial e^{-h}}$ denote the sequence valued of the Fourier modes of $\partial e^{-h}$, respectively, $\partial_{\theta}\partial e^{-h}$, and recall that their negative Fourier modes vanish:

[TABLE]

Using Plancherel and the notations above, we estimate the last term in (96) as follows.

[TABLE]

Next we estimate each term on the right hand side of (97) separately.

1. I.

Estimate the first term in (97):

[TABLE]

where in the first equality we use Plancherel, in the first inequality we use Young’s inequality, in the second equality we have used the fact that $\sup_{\overline{\Omega}}\left\|\widehat{\partial_{\theta}\partial e^{-h}}\right\|^{2}_{l_{1}(\mathbb{Z})}=\sup_{z\in\overline{\Omega}}\left(\sum_{k=1}^{{\infty}}k\lvert\partial\alpha_{k}(z)\rvert\right)^{2}$, and in the last equality we use the norm in (28) and the fact that $u_{-n}=\overline{u_{n}}$ as $u$ is real valued. 2. II.

Estimate the second term in (97):

[TABLE]

where in the first equality we have use Plancherel, in the first inequality we use Young’s inequality, in the second equality we have used the fact that $\sup_{\overline{\Omega}}\left\|\widehat{\partial e^{-h}}\right\|^{2}_{l_{1}(\mathbb{Z})}=\sup_{z\in\overline{\Omega}}\left(\sum_{k=0}^{{\infty}}\lvert\partial\alpha_{k}(z)\rvert\right)^{2}$, and in the last equality we use the fact that $u_{-n}=\overline{u_{n}}$ as $u$ is real valued. 3. III.

Estimate the third term in (97):

[TABLE]

where in the first equality we have use Plancherel, in the first inequality we use Young’s inequality, in the second equality we have used the fact that $\sup_{\overline{\Omega}}\left\|\widehat{\partial_{\theta}e^{-h}}\right\|^{2}_{l_{1}(\mathbb{Z})}=\sup_{z\in\overline{\Omega}}\left(\sum_{k=1}^{{\infty}}k\lvert\alpha_{k}(z)\rvert\right)^{2}$, and in the last equality we use the norm in (28) and the fact that $u_{-n}=\overline{u_{n}}$ as $u$ is real valued. 4. IV.

Estimate the last term in (97):

[TABLE]

where in the first equality we use Plancherel, in the first inequality we use Young’s inequality, in the second equality we have used the fact that $\sup_{\overline{\Omega}}\left\|\widehat{e^{-h}}\right\|^{2}_{l_{1}(\mathbb{Z})}=\sup_{z\in\overline{\Omega}}\left(\sum_{k=0}^{{\infty}}\lvert\alpha_{k}(z)\rvert\right)^{2}$, and in the second to last equality we use $u_{-n}=\overline{u_{n}}$ as $u$ is real valued.

The last term of (96) is thus estimated by

[TABLE]

Therefore the estimate (96) yields,

[TABLE]

Using the above estimate of $||\partial(e^{-G}{\bf u})||^{2}_{1,0}$ and estimate (93), the norm

[TABLE]

The mapping property for the case $p=\frac{1}{2}$ and $q=1$, follows from interpolation and the continuous embeddings $\displaystyle l^{2,1}(\mathbb{N};H^{1}(\Omega))\xhookrightarrow{}l^{2,\frac{1}{2}}(\mathbb{N};H^{1}(\Omega))\xhookrightarrow{}l^{2}(\mathbb{N};H^{1}(\Omega)).$ ∎