# A source reconstruction method in two dimensional radiative transport using
boundary data measured on an arc

Hiroshi Fujiwara  Graduate School of Informatics, Kyoto University, Yoshida
Honmachi, Sakyo-ku, Kyoto 606-8501, Japan<EMAIL_ADDRESS>,
Kamran Sadiq  Johann Radon Institute of Computational and Applied Mathematics
(RICAM), Altenbergerstrasse 69, 4040 Linz, Austria
<EMAIL_ADDRESS>and  Alexandru Tamasan  Department of
Mathematics, University of Central Florida, Orlando, 32816 Florida, USA
<EMAIL_ADDRESS>
###### Abstract.

We consider an inverse source problem in the stationary radiating transport
through a two dimensional absorbing and scattering medium. Of specific
interest, the exiting radiation is measured on an arc. The attenuation and
scattering properties of the medium are assumed known. For scattering kernels
of finite Fourier content in the angular variable, we show how to
quantitatively recover the part of the isotropic sources restricted to the
convex hull of the measurement arc. The approach is based on the Cauchy
problem with partial data for a Beltrami-like equation associated with
$A$-analytic maps in the sense of Bukhgeim, and extends authors’ previous work
to this specific partial data case. The robustness of the method is
demonstrated by the results of several numerical experiments.

###### Key words and phrases:

Radiative transport, source reconstruction, scattering, $A$-analytic maps,
Hilbert transform, Bukhgeim-Beltrami equation, optical tomography, optical
molecular imaging

###### 2010 Mathematics Subject Classification:

Primary 35J56, 30E20; Secondary 45E05

##  1\. Introduction

Let $\Omega$ be a strictly convex planar domain. In the steady state case,
when generated solely by a source $f$ inside $\Omega$, the density
$u(z,\boldsymbol{\theta})$ of particles at $z$ traveling in the direction
$\boldsymbol{\theta}$ through an absorbing and scattering domain $\Omega$
solves the stationary transport boundary value problem

(1) $\displaystyle\boldsymbol{\theta}\cdot\nabla
u(z,\boldsymbol{\theta})+a(z,\boldsymbol{\theta})u(z,\boldsymbol{\theta})-\int_{{{\mathbf{S}}^{1}}}k(z,\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})u(z,\boldsymbol{\theta}^{\prime})d\boldsymbol{\theta}^{\prime}=f(z,\boldsymbol{\theta}),\quad(z,\boldsymbol{\theta})\in\Omega\times{{\mathbf{S}}^{1}},$
$\displaystyle u\lvert_{\Gamma_{-}}=0,$

where
$\varGamma_{-}:=\\{(\zeta,\boldsymbol{\theta})\in\varGamma\times{{\mathbf{S}}^{1}}:\,\nu(\zeta)\cdot\boldsymbol{\theta}<0\\}$
with $\nu$ being the outer unit normal field at the boundary. The boundary
condition indicates that no radiation is coming from outside the domain.
Throughout, the measure on the circle is normalized to
$\int_{{{\mathbf{S}}^{1}}}d\boldsymbol{\theta}=1$.

The boundary value problem (1) is known to have a unique solution under
various “subcritical” assumptions, e.g., [11, 9, 10, 1, 22], with a general
result in [36] showing that, for an open and dense set of coefficients $a\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$ and $k\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}})$, the
problem (1) has a unique solution $u\in L^{2}(\Omega\times{{\mathbf{S}}^{1}})$
for any $f\in L^{2}(\Omega\times{{\mathbf{S}}^{1}})$. Some of our arguments in
the reconstruction method here require solutions $u\in
C^{1,\mu}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$, $\frac{1}{2}<\mu<1$. We
revisit the arguments in [36] and show that such a regularity can be achieved
for sources $f\in W^{2.p}(\Omega\times{{\mathbf{S}}^{1}})$, $p>4$; see Theorem
2.2 (ii) below.

For an arc $\Lambda$ of the boundary of $\Omega$, we consider the inverse
problem of the determining $f$ from measurements $g$ of exiting radiation on
$\Lambda$,

(2) $\displaystyle u|_{\Lambda_{+}}=g,$

where
$\Lambda_{+}:=\\{(z,\boldsymbol{\theta})\in\Lambda\times{{\mathbf{S}}^{1}}:\,\nu(z)\cdot\boldsymbol{\theta}>0\\}$
with $\nu$ being the outer unit normal field at the boundary. This is an
inverse source problem with partial boundary data, which applies to multiple
anisotropic scattering media of non-small anisotropy far from diffusive but
also single-scattering regime with applications to Optical Molecular Imaging
and Optical Tomography [3, 14].

Except for the results in Section 2, which concern the forward problem, in
this work the source and attenuation coefficient are assumed isotropic,
$f=f(z)$ and $a=a(z)$, and that the scattering kernel
$k(z,\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})=k(z,\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime})$
depends polynomially on the angle between the directions,

(3) $\displaystyle
k(z,\cos\theta)=k_{0}+2\sum_{n=1}^{M}k_{-n}(z)\cos(n\theta),$

for some fixed integer $M\geq 1$. Moreover, the functions $a,f,k$ are assumed
real valued. These assumptions occur naturally in radiative transfer models in
Optics, see, e.g., [8].

Our main result, Theorem 4.1, shows that $u|_{\Lambda_{+}}$ determines both
$f$ and $u$ in the convex hull of $\Lambda$, and provides a method of
reconstruction. Specific to two dimensional domains, our approach is based on
the Cauchy problem with partial data for a Beltrami-like equation associated
with $A$-analytic maps in the sense of Bukhgeim [7], and extends the authors’
previous work [16], which used measurements on the entire boundary, to this
specific partial data case. The new insight is that, similar to the non-
scattering case [18], the trace $u|_{\Lambda}$ determines $u|_{L}$ provided
the scattering kernel has finite Fourier content as in (3). The role of the
finite Fourier content has been independently recognized in [23].

We are aware of only one existing work that reconstructs the source in the
presence of scattering using partial boundary data in stationary radiative
transport. Namely, in [35] the domain is a slab and data is available on each
side on sufficiently long intervals. In contrast, our method here assumes only
“one sided” boundary data, and does not require iterative solvability of the
forward problem.

As demonstrated by the numerical experiments in Section 5, the method is
robust in the sense that it reconstructs discontinuous sources, even when
embedded in media with a discontinuous absorption property, and a scattering
kernel of infinite Fourier content in the angular variable.

##  2\. Remarks on the existence and regularity of the forward problem

The well posedness in $L^{p}(\Omega\times{{\mathbf{S}}^{1}})$ of the boundary
value problem (1) relies on the following compactness result, proven in [36,
Lemma 2.4] for the case $p=2$ and $a$ and $k$ twice differentiable. In this
section we revisit the arguments in [36] for any $1<p<\infty$, and show that
they hold if the attenuation is merely _once_ differentiable. We work in two
dimensions but this is not essential. Adopting the notation in [36], let us
consider the operators

(4) $\displaystyle(x,\boldsymbol{\theta})$
$\displaystyle=\int_{-{\infty}}^{0}e^{-\int_{s}^{0}a(x+t\boldsymbol{\theta},\boldsymbol{\theta})dt}g(x+s\boldsymbol{\theta},\boldsymbol{\theta})ds,\quad\text{and}$
$\displaystyle[Kg](x,\boldsymbol{\theta})$
$\displaystyle=\int_{{{\mathbf{S}}^{1}}}k(x,\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})g(x,\boldsymbol{\theta}^{\prime})d\boldsymbol{\theta}^{\prime},$

where the intervening functions are extended by 0 outside $\Omega$.

Using the formal expansion

(5) $\displaystyle u=T_{1}^{-1}f$
$\displaystyle+T_{1}^{-1}KT_{1}^{-1}f+T_{1}^{-1}(KT_{1}^{-1}K)[I-T_{1}^{-1}K]^{-1}T_{1}^{-1}f,$

the well posed-ness in $L^{p}(\Omega\times{{\mathbf{S}}^{1}})$ of the boundary
value problem (1) reduces to the invertibility of $I-T_{1}^{-1}K$ in
$L^{p}(\Omega\times{{\mathbf{S}}^{1}})$.

To further simplify notations, let $\widehat{x}=\frac{x}{|x|}$, so that for
$y\neq x$ we have $y=x-|x-y|(\widehat{x-y})$.

######  Proposition 2.1.

Let $a\in C^{1}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$ and $k\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}})$.
Then the operator

(6) $\displaystyle
KT_{1}^{-1}K:L^{p}(\Omega\times{{\mathbf{S}}^{1}})\longrightarrow
W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})\mbox{ is bounded},\;1<p<\infty.$

###### Proof.

Using the definitions of $T_{1}^{-1}$ and $K$ above, and a change to polar
coordinates, one can write

(7)
$\displaystyle[KT_{1}^{-1}Kg](x,\boldsymbol{\theta})=\int_{{{\mathbf{S}}^{1}}}\int_{\Omega}\frac{\eta\left(x,|x-y|,(\widehat{x-y}),\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)}{|x-y|}g\left(y,\boldsymbol{\theta}^{\prime}\right)dyd\boldsymbol{\theta}^{\prime},$

where, for
$\left(x,r,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)\in\Omega\times[0,\infty)\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}}$,

$\displaystyle\eta\left(x,r,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)=e^{-\int^{r}_{0}a\left(x-t\boldsymbol{\alpha},\boldsymbol{\alpha}\right)dt}k\left(x,\boldsymbol{\theta},\boldsymbol{\alpha}\right)k\left(x-r\boldsymbol{\alpha},\boldsymbol{\alpha},\boldsymbol{\theta}^{\prime}\right).$

An application of the fundamental theorem of calculus to
$r\mapsto\eta\left(x,r,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)$
yields

(8)
$\displaystyle\displaystyle\eta\left(x,r,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)=\eta\left(x,0,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\right)-rk(x,\boldsymbol{\theta},\boldsymbol{\alpha})\eta_{1}\left(x,r,\boldsymbol{\alpha},\boldsymbol{\theta}^{\prime}\right),$

where, for $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2})$, the function
$\eta_{1}(x,r,\boldsymbol{\alpha},\boldsymbol{\theta}^{\prime})$ is defined by

$\displaystyle\int_{0}^{1}e^{-\int_{0}^{r\rho}a(x-t\boldsymbol{\alpha},\boldsymbol{\alpha})dt}\left[a(x-r\rho\boldsymbol{\alpha},\boldsymbol{\alpha})k(x-r\rho\boldsymbol{\alpha},\boldsymbol{\alpha},\boldsymbol{\theta}^{\prime})+\sum_{j=1}^{2}\alpha_{j}\frac{\partial
k}{\partial
x_{j}}(x-r\rho\boldsymbol{\alpha},\boldsymbol{\theta}^{\prime})\right]d\rho.$

Note that in $\eta_{1}$ there are no derivatives taken on $a$, whereas there
are first order derivatives on $k$.

The split of the kernel in (8) induces the split of the operator
$KT_{1}^{-1}K=A-B$ with

$\displaystyle[Ag](x,\boldsymbol{\theta})=\int_{{{\mathbf{S}}^{1}}}\int_{\Omega}\frac{k\left(x,\boldsymbol{\theta},\widehat{x-y}\right)k\left(x,\widehat{x-y},\boldsymbol{\theta}^{\prime}\right)}{|x-y|}g\left(y,\boldsymbol{\theta}^{\prime}\right)dyd\boldsymbol{\theta}^{\prime},\;\mbox{
and }$
$\displaystyle[Bg](x,\boldsymbol{\theta})=\int_{{{\mathbf{S}}^{1}}}\int_{\Omega}k\left(x,\boldsymbol{\theta},\widehat{x-y}\right)\eta_{1}\left(x,|x-y|,\widehat{x-y},\boldsymbol{\theta}^{\prime}\right)g\left(y,\boldsymbol{\theta}^{\prime}\right)dyd\boldsymbol{\theta}^{\prime}.$

Using

$\displaystyle\nabla\frac{1}{|x-y|}=\frac{\widehat{y-x}}{|x-y|^{2}},\quad\frac{\partial}{\partial
x_{j}}\left(\frac{x_{k}-y_{k}}{|x-y|}\right)=\frac{1}{|x-y|}\left(\delta_{jk}-\frac{(x_{j}-y_{j})(x_{k}-y_{k})}{|x-y|^{2}}\right),$

and the regularity $a\in C^{1}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$ and
$k\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}})$, by
a straightforward calculation of the derivatives, one can verify that
$\frac{\partial}{\partial\theta_{j}}B$ is an operator with bounded kernel,
while $\frac{\partial}{\partial x_{j}}B$ and
$\frac{\partial}{\partial\theta_{j}}A$ are operators with weakly singular
kernel, all of which are bounded in $L^{p}(\Omega\times{{\mathbf{S}}^{1}})$,
e.g., [21, Theorem VIII.3.1]. On the other hand, in addition to the terms with
weakly singular kernels, the derivatives $\frac{\partial}{\partial x_{j}}A$
also yield operators of the Calderón-Zygmund type,

(9)
$\displaystyle[\mathcal{C}v](x,\boldsymbol{\theta})=\int_{{\mathbf{S}}^{1}}\int_{\Omega}\frac{\phi(x,\widehat{x-y},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})}{|x-y|^{2}}v(y,\boldsymbol{\theta}^{\prime})dyd\boldsymbol{\theta}^{\prime},$

where the characteristic $\phi$ satisfies

(10)
$\displaystyle\sup_{\Omega\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}}}|\phi(x,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})|<\infty.$

The following lemma concludes the proof of the Proposition 2.1. ∎

######  Lemma 2.1.

Let $\mathcal{C}$ be the operator in (9) with the characteristic $\phi$
satisfying (10). Then $\mathcal{C}$ is bounded in
$L^{p}(\Omega\times{{\mathbf{S}}^{1}})$, $p>1$.

###### Proof.

For brevity let us introduce the following notations

(11)
$\displaystyle\phi_{\infty}(x,\boldsymbol{\alpha}):=\sup_{\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}\in{{\mathbf{S}}^{1}}}|\phi(x,\boldsymbol{\alpha},\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})|,\quad
g(y):=\int_{{\mathbf{S}}^{1}}|v(y,\boldsymbol{\theta}^{\prime})|d\boldsymbol{\theta}^{\prime},$

and note that, via the Hölder inequality,
$\displaystyle\|g\|_{L^{p}(\Omega)}\leq\|v\|_{L^{p}(\Omega\times{{\mathbf{S}}^{1}})}.$
Using the Calderon-Zygmund boundedness theorem [21, Theorem XI.3.1] (the third
inequality below), we estimate

$\displaystyle\|\mathcal{C}v\|^{p}_{L^{p}(\Omega\times{{\mathbf{S}}^{1}})}$
$\displaystyle\leq\int_{\Omega}\sup_{\boldsymbol{\theta}\in{{\mathbf{S}}^{1}}}|[\mathcal{C}v](x,\boldsymbol{\theta})|^{p}dx\leq\int_{\Omega}\left|\int_{{\mathbf{S}}^{1}}\int_{\Omega}\frac{\phi_{\infty}(x,\widehat{x-y})}{|x-y|^{2}}|v(y,\boldsymbol{\theta}^{\prime})|dyd\boldsymbol{\theta}^{\prime}\right|^{p}dx$
$\displaystyle=\int_{\Omega}\left|\int_{\Omega}\frac{\phi_{\infty}(x,\widehat{x-y})}{|x-y|^{2}}g(y)dy\right|^{p}dx\leq
C\|g\|^{p}_{L^{p}(\Omega)}\leq
C\|v\|^{p}_{L^{p}(\Omega\times{{\mathbf{S}}^{1}})}.$

∎

The following simple result is useful.

######  Lemma 2.2.

Let $X$ be a Banach space and $A:X\to X$ be bounded. Then $I-A$ has a bounded
inverse in $X$ if and only if $I-A^{2}$ has a bounded inverse in $X$.

As a consequence, for $\lambda\in\mathbb{C}$, the operator
$I-T_{1}^{-1}(\lambda K)$ is invertible in
$L^{p}(\Omega\times{{\mathbf{S}}^{1}})$ if and only if $I-(T_{1}^{-1}(\lambda
K))^{2}$ is invertible in $L^{p}(\Omega\times{{\mathbf{S}}^{1}})$. By
Proposition 2.1, $\displaystyle(T_{1}^{-1}(\lambda K))^{2}$ is compact for any
$z\in\mathbb{C}$. Since $\displaystyle I-(T_{1}^{-1}(\lambda K))^{2}$ is
invertible for $\lambda$ in a neighborhood of $0$, an application of the
meromorphic Fredholm alternative (in reflexive Banach spaces, e.g., [12])
yields the following result.

######  Theorem 2.1.

Let $p>1$, $a\in C^{1}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$, and $k\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}})$. At
least one of the following statements is true.

(i) $\displaystyle I-T_{1}^{-1}K$ is invertible in
$L^{p}(\Omega\times{{\mathbf{S}}^{1}})$.

(ii) there exists $\epsilon>0$ such that $\displaystyle I-T_{1}^{-1}(\lambda
K)$ is invertible in $L^{p}(\Omega\times{{\mathbf{S}}^{1}})$, for any
$0<|\lambda-1|<\epsilon$.

If $a\in C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$, then the
regularity of the solution $u$ of (1) increases with the regularity of $f$ as
follows.

######  Theorem 2.2.

Consider the boundary value problem (1) with $a\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$. For $p>1$, let $k\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}})$ be
such that $\displaystyle I-T_{1}^{-1}K$ is invertible in
$L^{p}(\Omega\times{{\mathbf{S}}^{1}})$, and let $u\in
L^{p}(\Omega\times{{\mathbf{S}}^{1}})$ in (5) be the solution of (1).

(i) If $f\in W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$, then $u\in
W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$.

(ii) If $f\in W^{2,p}(\Omega\times{{\mathbf{S}}^{1}})$, then $u\in
W^{2,p}(\Omega\times{{\mathbf{S}}^{1}})$.

###### Proof.

(i) Recall the representation (5) of the solution of (1),

$\displaystyle u=T_{1}^{-1}f$
$\displaystyle+T_{1}^{-1}KT_{1}^{-1}f+T_{1}^{-1}[KT_{1}^{-1}K](I-T_{1}^{-1}K)^{-1}T_{1}^{-1}f.$

It is easy to see that $T_{1}^{-1}$ and $K$ preserve the space
$W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$, so that the first two terms belong
to $W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$. Now, by Proposition 2.1, the
third term is also in $W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$.

(ii) For brevity we introduce the operators

(12) $\displaystyle
T_{0}^{-1}u(x,\boldsymbol{\theta})=\int_{-\infty}^{0}u(x+t\boldsymbol{\theta},\boldsymbol{\theta})dt,\quad
K_{j}u(x,\boldsymbol{\theta})=\int_{{\mathbf{S}}^{1}}\frac{\partial
k}{\partial
x_{j}}(x,\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})u(x,\boldsymbol{\theta}^{\prime})d\boldsymbol{\theta}^{\prime},$
$\displaystyle\widetilde{T}_{0}^{-1}u(x,\boldsymbol{\theta})=\int_{-\infty}^{0}u(x+t\boldsymbol{\theta},\boldsymbol{\theta})tdt,\quad\widehat{K}_{j}u(x,\boldsymbol{\theta})=\int_{{\mathbf{S}}^{1}}\frac{\partial
k}{\partial\theta_{j}}(x,\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})u(x,\boldsymbol{\theta}^{\prime})d\boldsymbol{\theta}^{\prime},\;j=1,2.$

It is easy to see that $T_{0}^{-1},\widetilde{T}_{0}^{-1},K_{j}$ and
$\widehat{K}_{j}$ preserve $W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$.

By evaluating (1) at $x+t\boldsymbol{\theta}$ and integrating in $t$ from
$-{\infty}$ to $0$, the problem (1) with zero incoming fluxes is equivalent to
the integral equation.

(13) $\displaystyle u+T_{0}^{-1}(au)-T_{0}^{-1}Ku=T_{0}^{-1}f.$

For $f\in W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$, according to part (i),
$u_{x_{j}}\in L^{p}(\Omega\times{{\mathbf{S}}^{1}})$. In particular
$u_{x_{j}}$ solves the integral equation

(14) $\displaystyle
u_{x_{j}}+T_{0}^{-1}(au_{x_{j}})-T_{0}^{-1}Ku_{x_{j}}=T_{0}^{-1}f_{x_{j}}-T_{0}^{-1}(a_{x_{j}}u)+T_{0}^{-1}K_{j}u$

Moreover, since $a\in C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$, $k\in
C^{2}(\overline{\Omega}\times{{\mathbf{S}}^{1}}\times{{\mathbf{S}}^{1}})$, and
$f\in W^{2,p}(\Omega\times{{\mathbf{S}}^{1}})$, the right-hand-side of (14)
lies in $W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$. By applying part (i) above,
we get that the unique solution to (14)

(15) $\displaystyle u_{x_{j}}\in
W^{1,p}(\Omega\times{{\mathbf{S}}^{1}}),\;j=1,2.$

For $f\in W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$, also according to part
(i), $u_{\theta_{j}}\in L^{p}(\Omega\times{{\mathbf{S}}^{1}})$. In particular
$u_{\theta_{j}}$ is the unique solution of the integral equation

(16) $\displaystyle
u_{\theta_{j}}+T_{0}^{-1}(au_{\theta_{j}})=T_{0}^{-1}f_{\theta_{j}}-\widetilde{T}_{0}^{-1}(au_{x_{j}})-\widetilde{T}_{0}^{-1}(a_{x_{j}}u)-T_{0}^{-1}(a_{\theta_{j}}u)+\widetilde{T}_{0}^{-1}K_{j}u-T_{0}^{-1}\widehat{K}_{j}u,$

which is of the type (13) with $K=0$. Moreover, since $f\in
W^{2,p}(\Omega\times{{\mathbf{S}}^{1}})$, and, according to (15),
$u_{x_{j}}\in W^{1,p}(\Omega\times{{\mathbf{S}}^{1}}),\;j=1,2$, the right-
hand-side of (16) lies in $W^{1,p}(\Omega\times{{\mathbf{S}}^{1}})$. Again, by
applying part (i), we get

$u_{\theta_{j}}\in W^{1,p}(\Omega\times{{\mathbf{S}}^{1}}),\;j=1,2.$

Thus, $u\in W^{2,p}(\Omega\times{{\mathbf{S}}^{1}})$.

∎

##  3\. Ingredients from $A$-analytic theory

For $0<\mu<1$, and $\varGamma$ (some part of) the boundary $\partial\Omega$ we
consider the Banach spaces:

$\displaystyle l^{1,1}_{{\infty}}(\varGamma)$ $\displaystyle:=\left\\{{\bf
v}\;:\left\|{\bf
v}\right\|_{l^{1,1}_{{\infty}}(\varGamma)}:=\sup_{\xi\in\varGamma}\sum_{j=1}^{{\infty}}j\lvert
v_{-j}(\xi)\rvert<{\infty}\right\\},$ $\displaystyle C^{\mu}(\varGamma;l_{1})$
$\displaystyle:=\left\\{{\bf v}:\sup_{\xi\in\varGamma}\lVert{\bf
v}(\xi)\rVert_{\displaystyle
l_{1}}+\underset{{\begin{subarray}{c}\xi,\eta\in\varGamma\\\
\xi\neq\eta\end{subarray}}}{\sup}\frac{\lVert{\bf v}(\xi)-{\bf
v}(\eta)\rVert_{\displaystyle l_{1}}}{|\xi-\eta|^{\mu}}<{\infty}\right\\}.$

We similarly consider $l^{1,1}_{{\infty}}(\overline{\Omega})$,
$C^{\mu}(\overline{\Omega};l_{1})$, and
$C^{\mu}(\overline{\Omega};l_{{\infty}})$.

A sequence valued map $\Omega\ni z\mapsto{\bf v}(z):=\langle
v_{0}(z),v_{-1}(z),v_{-2}(z),...\rangle$ in
$C(\overline{\Omega};l_{\infty})\cap C^{1}(\Omega;l_{\infty})$ is called
$\mathcal{L}^{2}$-analytic (in the sense of Bukhgeim), if

(17) $\overline{\partial}{\bf v}(z)+\mathcal{L}^{2}\partial{\bf v}(z)=0,\quad
z\in\Omega,$

where $\mathcal{L}$ is the left shift operator,
$\displaystyle\mathcal{L}\langle v_{0},v_{-1},v_{-2},\cdots\rangle=\langle
v_{-1},v_{-2},\cdots\rangle,$ and
$\mathcal{L}^{2}=\mathcal{L}\circ\mathcal{L}$. Note that we use the sequences
of non-positive indexes to conform with the original notation in Bukhgeim’s
work [7].

Analogous to the analytic maps, the $\mathcal{L}^{2}$-analytic maps are
determined by their boundary values via a Cauchy-like integral formula [7].
Following [13], the Bukhgeim-Cauchy operator $\mathcal{B}$ acting on ${\bf
v}=\langle v_{0},v_{-1},v_{-2},...\rangle$ is defined component-wise for
$n\leq 0$ by

(18) $\displaystyle(\mathcal{B}{\bf v})_{n}(z)$
$\displaystyle:=\frac{1}{2\pi\mathrm{i}}\int_{\varGamma}\frac{v_{n}(\zeta)}{\zeta-z}d\zeta+\frac{1}{2\pi\mathrm{i}}\int_{\varGamma}\left\\{\frac{d\zeta}{\zeta-z}-\frac{d\overline{\zeta}}{\overline{\zeta-z}}\right\\}\sum_{j=1}^{\infty}v_{n-2j}(\zeta)\left(\frac{\overline{\zeta-z}}{\zeta-z}\right)^{j},\;z\in\Omega.$

As shown in [32, Theorem 2.2], if ${\bf v}=\langle
v_{0},v_{-1},v_{-2},...\rangle\in l^{1,1}_{{\infty}}(\varGamma)\cap
C^{\mu}(\varGamma;l_{1})$, then $\mathcal{B}{\bf v}\in
C^{1,\mu}(\Omega;l_{\infty})\cap C(\overline{\Omega};l_{\infty})$ is
$\mathcal{L}^{2}$-analytic in $\Omega$. If ${\bf v}\in
C^{1,\mu}(\Omega;l_{\infty})\cap C(\overline{\Omega};l_{\infty})$ is
$\mathcal{L}^{2}$-analytic in $\overline{\Omega}$, then ${\bf
v}(z)=\mathcal{B}{\bf v}(z)$, for $z\in\overline{\Omega}$.

Also similar to the analytic maps, the traces on the boundary of
$\mathcal{L}^{2}$-analytic maps satisfy some constraints, which can be
expressed in terms of a corresponding Hilbert transform introduced in [31].
More precisely, the Bukhgeim-Hilbert transform $\mathcal{H}$ is defined
component-wise for $n\leq 0$ by

(19) $\displaystyle(\mathcal{H}{\bf v})_{n}(\xi)$
$\displaystyle=\frac{1}{\pi}\int_{\varGamma}\frac{v_{n}(\zeta)}{\zeta-\xi}d\zeta+\frac{1}{\pi}\int_{\varGamma}\left\\{\frac{d\zeta}{\zeta-\xi}-\frac{d\overline{\zeta}}{\overline{\zeta-\xi}}\right\\}\sum_{j=1}^{\infty}v_{n-2j}(\zeta)\left(\frac{\overline{\zeta-\xi}}{\zeta-\xi}\right)^{j},\;\xi\in\varGamma,$

and we refer to [31] for its mapping properties. For the proof of the theorem
below we refer to [31, Theorem 3.2].

######  Theorem 3.1.

Let ${\bf v}=\langle v_{0},v_{-1},v_{-2},...\rangle\in
l^{1,1}_{{\infty}}(\varGamma)\cap C^{\mu}(\varGamma;l_{1})$ be defined on the
boundary $\varGamma$. Then ${\bf v}$ is the boundary value of an
$\mathcal{L}^{2}$-analytic function if and only if

(20) $\displaystyle(I+\mathrm{i}\mathcal{H}){\bf v}={\bf{0}}.$

In addition to $\mathcal{L}^{2}$-analytic maps, another ingredient consists in
the one-to-one relation between solutions ${\bf u}:=\langle
u_{0},u_{-1},u_{-2},...\rangle$ to

(21) $\displaystyle\overline{\partial}u_{-n}(z)+\partial
u_{-n-2}(z)+a(z)u_{-n-1}(z)$ $\displaystyle=0,\quad z\in\Omega,\;n\geq 0.$

and the $\mathcal{L}^{2}$-analytic map ${\bf v}$ satisfying

(22) $\displaystyle\overline{\partial}v_{-n}(z)+\partial v_{-n-2}(z)$
$\displaystyle=0,\quad z\in\Omega,\;n\geq 0;$

see [33, Lemma 4.2] for details. The relation can be expressed via the
convolutions

(23) $\displaystyle v_{-n}(z)$
$\displaystyle=\sum_{j=0}^{{\infty}}\alpha_{j}(z)u_{-n-j}(z),\quad
z\in\Omega,\;n\geq 0,$ $\displaystyle u_{-n}(z)$
$\displaystyle=\sum_{j=0}^{{\infty}}\beta_{j}(z)v_{-n-j}(z),\quad
z\in\Omega,\;n\geq 0,$

where $\alpha_{j}$’s and $\beta_{j}$’s are the Fourier modes of $e^{\mp h}$,

(24) $\displaystyle
e^{-h(z,\boldsymbol{\theta})}:=\sum_{m=0}^{{\infty}}\alpha_{m}(z)e^{\mathrm{i}m\theta},\quad
e^{h(z,\boldsymbol{\theta})}:=\sum_{m=0}^{{\infty}}\beta_{m}(z)e^{\mathrm{i}m\theta},\quad(z,\theta)\in\overline{\Omega}\times{{\mathbf{S}}^{1}},$

with $h$ defined by

(25) $\displaystyle
h(z,\boldsymbol{\theta}):=Da(z,\boldsymbol{\theta})-\frac{1}{2}\left(I-\mathrm{i}H\right)Ra(z\cdot\boldsymbol{\theta}^{\perp},\boldsymbol{\theta}^{\perp}).$

In the above formula, $\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 $H$ is the (infinite)
Hilbert transform

(26) $\displaystyle
Hf(x)=\frac{1}{\pi}\int_{-{\infty}}^{{\infty}}\frac{f(s)}{x-s}ds$

taken in the first variable and evaluated at
$s=z\cdotp\boldsymbol{\theta}^{\perp}$. The function $h$ appeared first in
[25] and enjoys the crucial property of having vanishing negative Fourier
modes. We refer to [33, Lemma 4.1] for the properties of $h$ used in here.

##  4\. Reconstruction in $\Omega^{+}$ of a sufficiently smooth isotropic
source $f$

The method of reconstruction proposed here considers the operator
$[I-\mathrm{i}H_{t}]$, where $H_{t}$ is the finite Hilbert transform

(27) $\displaystyle
H_{t}g(x)=\frac{1}{\pi}\int_{-l}^{l}\frac{g(s)}{x-s}ds,\quad x\in(-l,l),$

with the integral understood in the sense of principal value. It is well-known
([39]) that $\mathrm{i}H_{t}$ is a bounded operator in $L^{2}(-l,l)$ with
spectrum $[-1,1]$, see [19, 29]. However, $1$ is not in the point spectrum,
see [41] for a proof based on a Riemann-Hilbert problem. The arguments below
(from [18, Proposition 2.1]) use the unitary property in $L^{2}(\mathbb{R})$
of the (infinite) Hilbert transform $H$ in (26). We repeat them here for
reader’s convenience.

######  Proposition 4.1.

In $L^{2}(-l,l)$, $\displaystyle Ker[I-\mathrm{i}H_{t}]=\\{0\\}$.

###### Proof.

Let $f\in L^{2}(-l,l)$ be extended by zero to the entire real line. If $f\in
Ker[I-\mathrm{i}H_{t}]$, then $|f|=|H_{t}f|=|Hf|$ a.e. in $(-l,l)$, and

$\displaystyle\left\|f\right\|_{L^{2}(-l,l)}^{2}$
$\displaystyle=\left\|f\right\|_{L^{2}(\mathbb{R})}^{2}=\left\|Hf\right\|_{L^{2}(\mathbb{R})}^{2}=\int_{-l}^{l}|H_{t}f|^{2}dx+\int_{|x|>l}|Hf(x)|^{2}dx$
$\displaystyle=\int_{-l}^{l}|f|^{2}dx+\int_{|x|>l}|Hf(x)|^{2}dx.$

The latter identity shows that $|Hf(x)|=0$, for $|x|>l$. Thus, for $|x|>l$,

$\displaystyle 0$
$\displaystyle=\int_{-l}^{l}\frac{f(s)}{x-s}ds=\frac{1}{x}\int_{-l}^{l}f(s)\sum_{j=0}^{{\infty}}\left(\frac{s}{x}\right)^{j}ds=\sum_{j=0}^{{\infty}}\frac{1}{x^{j+1}}\int_{-l}^{l}f(s)s^{j}ds,$

yielding $f$ orthogonal on any polynomial. By density of polynomials in
$L^{2}(-l,l)$, $f=0$. ∎

Recall the boundary value problem (1):

(28) $\displaystyle\boldsymbol{\theta}\cdot\nabla
u(z,\boldsymbol{\theta})+a(z)u(z,\boldsymbol{\theta})-\int_{{{\mathbf{S}}^{1}}}k(z,\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime})u(z,\boldsymbol{\theta}^{\prime})d\boldsymbol{\theta}^{\prime}=f(z),\quad(z,\boldsymbol{\theta})\in\Omega\times{{\mathbf{S}}^{1}},$
$\displaystyle u\lvert_{\Gamma_{-}}=0,$

for an isotropic source $f$ and attenuation coefficient $a$, and with a
scattering kernel $k$ of the type (3),

(29) $\displaystyle
k(z,\cos\theta)=k_{0}(z)+2\sum_{n=1}^{M}k_{-n}(z)\cos(n\theta),$

for some fixed integer $M\geq 1$.

We assume that $a,k_{0},k_{-1},...,k_{-M}\in C^{2}(\overline{\Omega})$ are
such that the forward problem (28) has a unique solution $u\in
L^{p}(\Omega\times{{\mathbf{S}}^{1}})$ for any $f\in L^{p}(\Omega)$, see
Theorem 2.1. We also assume an _unknown_ source of a priori regularity $f\in
W^{2,p}(\overline{\Omega})$, $p>4$. According to Theorem 2.2 part (ii), $u\in
C^{1,\mu}(\Omega\times{{\mathbf{S}}^{1}})$ with $\mu>1/2$. In agreement with
the physics model, the functions $a,f,k$ are further assumed real valued, so
that the solution $u$ is also real valued. Note that, since $k(z,\cos\theta)$
in (29) is both real valued and even in $\theta$, the coefficient $k_{-n}$ in
(29) is the $(-n)^{th}$ Fourier coefficient of $k(z,\cos(\cdot))$. Moreover
$k_{-n}$ is real valued, and
$k_{n}(z)=k_{-n}(z)=\frac{1}{2\pi}\int_{-\pi}^{\pi}k(z,\cos\theta)e^{in\theta}d\theta$.

Let $u(z,\boldsymbol{\theta})=\sum_{-\infty}^{\infty}u_{n}(z)e^{in\theta}$ be
the formal Fourier series representation of the solution of (28) 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

(30) $\displaystyle\Omega\ni z\mapsto{\bf u}(z)$ $\displaystyle:=\langle
u_{0}(z),u_{-1}(z),u_{-2}(z),...\rangle.$

Consider the decomposition of the advection operator
$\boldsymbol{\theta}\cdot\nabla=e^{-\mathrm{i}\theta}\overline{\partial}+e^{\mathrm{i}\theta}\partial$,
where $\overline{\partial}=(\partial_{x}+\mathrm{i}\partial_{y})/2$ and
$\partial=(\partial_{x}-\mathrm{i}\partial_{y})/2$ are derivatives in the
spatial domain. By identifying the Fourier coefficients of the same order, the
equation (28) reduces to the system:

(31) $\displaystyle\overline{\partial}u_{1}(z)+\partial
u_{-1}(z)+a(z)u_{0}(z)$ $\displaystyle=k_{0}(z)u_{0}(z)+f(z),$ (32)
$\displaystyle\overline{\partial}u_{-n}(z)+\partial
u_{-n-2}(z)+a(z)u_{-n-1}(z)$ $\displaystyle=k_{-n-1}(z)u_{-n-1}(z),\quad 0\leq
n\leq M-1,$ (33) $\displaystyle\overline{\partial}u_{-n}(z)+\partial
u_{-n-2}(z)+a(z)u_{-n-1}(z)$ $\displaystyle=0,\qquad n\geq M.$

Without loss of generality we consider cartesian coordinates such that
$\Lambda$ lies in the upper half plane with endpoints on the real axis, and
let $L=(-l,l)$ be the segment joining the endpoints of the arc. Let
$\Omega^{+}=\\{z\in\Omega:\;\operatorname{\mathbb{I}m}{z}>0\\}$ denote the
convex hull of $\Lambda$, and note that $\partial\Omega^{+}=\Lambda\cup L$.

To simplify the statement of the next result, for each $n\geq 0$, let us
introduce the functions $F_{-n}(z)$ defined for $z\neq\pm l$ in $\Lambda\cup
L$ by

(34) $\displaystyle
F_{-n}(z):=\frac{1}{\mathrm{i}\pi}\int_{\Lambda}\frac{u_{-n}(\zeta)}{\zeta-z}d\zeta+\frac{1}{\mathrm{i}\pi}\int_{\Lambda}\left\\{\frac{d\zeta}{\zeta-z}-\frac{d\overline{\zeta}}{\overline{\zeta}-\overline{z}}\right\\}\sum_{j=1}^{\infty}u_{-n-2j}(\zeta)\left(\frac{\overline{\zeta}-\overline{z}}{\zeta-z}\right)^{j}.$

For $z\in\Lambda$, the first integral is in the sense of principal value. Note
that $F_{-n}$ is directly determined by the data $u_{-n}\lvert_{\Lambda}$,
$n\geq 0$.

The proof of the following result is constructive and provides the basis of
the reconstruction method implemented in Section 5.

######  Theorem 4.1.

Let $\Omega\subset\mathbb{R}^{2}$ be a strictly convex bounded domain,
$\Lambda$ be an arc of its boundary, and $\Omega^{+}$ be the convex hull of
$\Lambda$. Consider the boundary value problem (28) for some known real valued
$a,k_{0},k_{-1},...,k_{-M}\in C^{2}(\overline{\Omega})$ such that (28) is
well-posed. If the unknown source $f$ is real valued and
$W^{2,p}(\overline{\Omega})$-regular, with $p>4$, then $u\lvert_{\Lambda_{+}}$
uniquely determines $f$ in $\Omega^{+}$.

###### Proof.

Let $u$ be the solution of the boundary value problem (28) and ${\bf
u}=\langle u_{0},u_{-1},u_{-2},...\rangle$ be the sequence valued map of its
non-positive Fourier modes, Since $f\in W^{2,p}(\overline{\Omega})$, $p>4$,
then by Proposition 2.2 (ii), $u\in W^{2,p}(\Omega\times{{\mathbf{S}}^{1}})$.
By the Sobolev embedding, $u\in
C^{1,\mu}(\overline{\Omega}\times{{\mathbf{S}}^{1}})$ with
$\mu=1-\frac{2}{p}>\frac{1}{2}$, and thus, by [31, Proposition 4.1 (i)], ${\bf
u}\in l^{1,1}_{{\infty}}(\partial\Omega^{+})\cap
C^{\mu}(\partial\Omega^{+};l_{1})$.

We note that the shifted sequence valued map $\mathcal{L}^{M}{\bf u}$ solves

(35) $\displaystyle\overline{\partial}\mathcal{L}^{M}{\bf
u}(z)+\mathcal{L}^{2}\partial\mathcal{L}^{M}{\bf
u}(z)+a(z)\mathcal{L}^{M+1}{\bf u}(z)=\mathbf{0},\quad z\in\Omega,$

and then the associated sequence valued map $\mathcal{L}^{M}{\bf
v}=(v_{-M},v_{-M-1},v_{-M-2}...)$ defined by the convolutions (23) solves

(36) $\displaystyle\overline{\partial}v_{-n}(z)+\partial v_{-n-2}(z)$
$\displaystyle=0,\quad z\in\Omega,\;n\geq M.$

In particular, $\mathcal{L}^{M}{\bf v}$ is $\mathcal{L}^{2}$-analytic.

By (2), the data $u|_{\Lambda_{+}}=g$ on $\Lambda_{+}$ determines
$\mathcal{L}^{M}{\bf u}$ on $\Lambda$. By the convolution formula (23) for
$n\geq M$, $\mathcal{L}^{M}{\bf u}\lvert_{\Lambda}$ determines the traces
$\mathcal{L}^{M}{\bf v}\in l^{1,1}_{{\infty}}(\Lambda)\cap
C^{\mu}(\Lambda;l_{1})$ on $\Lambda$.

Since $\mathcal{L}^{M}{\bf v}\in l^{1,1}_{{\infty}}(\partial\Omega^{+})\cap
C^{\mu}(\partial\Omega^{+};l_{1})$ is the boundary value of an
$\mathcal{L}^{2}$-analytic function in $\Omega^{+}$, then the necessity part
of Theorem 3.1 yields

(37) $\displaystyle[I+\mathrm{i}\mathcal{H}]\mathcal{L}^{M}{\bf v}={\bf{0}},$

where $\mathcal{H}$ is the Bukhgeim-Hilbert transform in (19).

We consider (37) on $L=(-l,l)$, where for each $x\in(-l,l)$ and $n\geq M$, the
$n$-th component yields

$\displaystyle
v_{-n}(x)-\frac{\mathrm{i}}{\pi}\int_{-l}^{l}\frac{v_{-n}(s)}{x-s}ds=-\frac{\mathrm{i}}{\pi}\int_{\Lambda}\frac{v_{-n}(\zeta)}{\zeta-x}d\zeta$
$\displaystyle-\frac{\mathrm{i}}{\pi}\int_{\Lambda}\left\\{\frac{d\zeta}{\zeta-x}-\frac{d\overline{\zeta}}{\overline{\zeta}-x}\right\\}\sum_{j=1}^{\infty}v_{-n-2j}(\zeta)\left(\frac{\overline{\zeta}-x}{\zeta-x}\right)^{j}$
(38)
$\displaystyle-\frac{\mathrm{i}}{\pi}\int_{-l}^{l}\left\\{\frac{d\zeta}{\zeta-x}-\frac{d\overline{\zeta}}{\overline{\zeta}-x}\right\\}\sum_{j=1}^{\infty}v_{-n-2j}(\zeta)\left(\frac{\overline{\zeta}-x}{\zeta-x}\right)^{j}.$

Since the last integral in (38) ranges over the reals, it vanishes. The
remaining two integrals in the right hand side give $F_{-n}(x)$ in (34), and
(38) becomes

(39) $\displaystyle[I-\mathrm{i}H_{t}](v_{-n})(x)=$ $\displaystyle
F_{-n}(x),\quad x\in L,\quad n\geq M,$

where $H_{t}$ is the finite Hilbert transform in (27). For each $n\geq M$, by
Proposition 4.1, $v_{-n}\big{\lvert}_{L}$ is determined as the unique solution
in $L^{2}(-l,l)$ of (39).

Note that the equation (39) may not have any solution for an arbitrary right
hand side in $L^{2}(-l,l)$. However, in our inverse problem, the function
$F_{-n}$ already belongs to the range of $I-\mathrm{i}H_{t}$, so that the
solution exists. Moreover, since the range is open, (39) is uniquely solvable
in a sufficiently small $L^{2}$-neighborhood of $F_{-n}$.

With $v_{-n}$ now known on $\Lambda\cup L$ for $n\geq M$, we apply the
Bukhgeim-Cauchy integral formula (18) to find $v_{-n}$ for $n\geq M$ in
$\Omega^{+}$.

Using again the convolution formula (23), now in $\Omega^{+}$, we determine
$u_{-n}$ for $n\geq M$ inside $\Omega^{+}$. In particular we recovered
$u_{-M-1},u_{-M}$.

Recall that $u_{0},u_{-1},u_{-2},\cdots,u_{-M},u_{-M-1}$ satisfy

(40a) $\displaystyle\overline{\partial}u_{-M+j}$ $\displaystyle=-\partial
u_{-M+j-2}-\left[(a-k_{-M+j-1})u_{-M+j-1}\right],\quad 1\leq j\leq M,$ (40b)
$\displaystyle u_{-M+j}\lvert_{\Lambda}$ $\displaystyle=g_{-M+j}.$

We solve (40) iteratively for $j=1,2,...,M$, as a Cauchy problem for the
$\overline{\partial}$-equation with partial boundary data on $\Lambda$,

(41a) $\displaystyle\overline{\partial}w$
$\displaystyle=\Psi,\quad\text{in}\;\Omega^{+},$ (41b) $\displaystyle w$
$\displaystyle=g\quad\text{on}\;\Lambda,$

via the Cauchy-Pompeiu formula [40]:

(42) $\displaystyle
w(z)=\frac{1}{2\pi\mathrm{i}}\int_{\partial\Omega^{+}}\frac{w(\zeta)}{\zeta-z}d\zeta-\frac{1}{\pi}\iint_{\Omega^{+}}\frac{\Psi(\zeta)}{\zeta-z}d\xi
d\eta,\quad\zeta=\xi+\mathrm{i}\eta,\quad z\in\Omega^{+}.$

For $\Psi\in L^{p}(\Omega)$, $p>2$, and $g\in L^{p}(\Lambda)$, any $w$ defined
by (42) solves (41a). However, for the boundary condition (41b) to be
satisfied the following compatibility condition needs to hold: By taking the
limit $\displaystyle\Omega^{+}\ni z\rightarrow z_{0}\in\partial\Omega^{+}$ in
(42) and using the Sokhotski-Plemelj formula [24] in the first integral, and
the continuous dependence on $z$ of the area integral [40, Theorem 1.19], the
trace $w|_{\partial\Omega^{+}}$ and $\Psi$ must satisfy

$\displaystyle
w(z_{0})=\frac{1}{2\pi\mathrm{i}}\int_{\partial\Omega^{+}}\frac{w(\zeta)}{\zeta-
z_{0}}d\zeta+\frac{1}{2}w(z_{0})-\frac{1}{\pi}\iint_{\Omega^{+}}\frac{\Psi(\xi,\eta)}{(\xi-
z_{0})+\mathrm{i}\eta}d\xi d\eta,\quad z_{0}\in\partial\Omega^{+}.$

In our inverse problem this compatibility condition is already satisfied for
$z_{0}\in\Lambda$. We use this compatibility condition for $z_{0}\in L$, to
recover the missing boundary data $w|_{L}$. More precisely, by Proposition
4.1, $w|_{L}$ is the unique solution of

(43)
$\displaystyle[I-\mathrm{i}H_{t}]w(z_{0})=\frac{1}{\pi\mathrm{i}}\int_{\Lambda}\frac{g(\zeta)}{\zeta-
z_{0}}d\zeta-\frac{2}{\pi}\iint_{\Omega^{+}}\frac{\Psi(\xi,\eta)}{(\xi-
z_{0})+\mathrm{i}\eta}d\xi d\eta,\quad z_{0}\in L.$

If $g\in L^{p}(\Lambda)$, $p\geq 2$, then (43) provides a unique solution
$w|_{L}\in L^{p}(-l,l)$. Moreover, for $\Psi\in L^{p}(\Omega)$, $p>2$, the
solution $w$ of (41a) is provided by (42) and lies in $W^{1,p}(\Omega)$. In
the iteration, the right hand side of (40a) is again in $L^{p}(\Omega)$, and
the iteration can proceed.

We solve repeatedly (40) for $j=1,....M$ to recover $u_{-1}$ and $u_{0}$ in
$\Omega^{+}$. A priori, from the regularity of the solution of the forward
problem, we known that $u_{-1},u_{0}\in C^{1,\mu}(\overline{\Omega})$ so that
the source $f$ is recovered pointwise by

(44) $\displaystyle
f|_{\Omega^{+}}(z)=2\operatorname{\mathbb{R}e}\left(\partial
u_{-1}(z)\right)+\left(a(z)-k_{0}(z)\right)u_{0}(z),\quad z\in\Omega^{+}$

as a $C^{\mu}(\overline{\Omega^{+}})$-map. ∎

We summarize below in a stepwise fashion the reconstruction of $f$ in the
convex hull ${\Omega^{+}}$ of the boundary arc $\Lambda$. Recall that $L$ is
the segment joining the endpoints of $\Lambda$.

### Reconstruction procedure

Consider the data ${\bf u}\lvert_{\Lambda}$.

  1. (1)

Using formula (23), and data $\mathcal{L}^{M}{\bf u}$ on $\Lambda$, determine
the traces $\mathcal{L}^{M}{\bf v}\lvert_{\Lambda}$ on $\Lambda$.

  2. (2)

Recover the traces $\mathcal{L}^{M}{\bf v}\lvert_{L}$ pointwise on $L$ as
follows:

    1. (a)

Using $\mathcal{L}^{M}{\bf v}\lvert_{\Lambda}$, compute by formula (34), the
function $F_{-n}$, for each $n\geq M$.

    2. (b)

Recover for each $n\geq M$, the trace $v_{-n}\big{\lvert}_{L}$ by solving
(39).

  3. (3)

By Bukhgeim-Cauchy formula (18), extend $v_{-n}$ for $n\geq M$ from the
boundary $\Lambda\cup L$ to $\Omega^{+}$.

  4. (4)

Using again formula (23), now in $\Omega^{+}$, recover $u_{-n}$ for $n\geq M$
inside $\Omega^{+}$.

  5. (5)

Using $u_{-M-1},u_{-M}$, recover the modes
$u_{-M+1},u_{-M+2},\cdots,u_{-1},u_{0}$ recuresively as follows:

    1. (a)

Using data $u_{-M+1}\lvert_{\Lambda}$, recover the trace
$u_{-M+1}\big{\lvert}_{L}$ by solving (43).

    2. (b)

Using $u_{-M+1}\lvert_{\Lambda\cup L}$, recover $u_{-M+1}$ inside $\Omega^{+}$
by the Cauchy-Pompeiu formula (42).

    3. (c)

Now iterate the steps (5 a) and (5 b) to find the modes
$u_{-M+2},\cdots,u_{-2},u_{-1},u_{0}$ in $\Omega^{+}$.

  6. (6)

Recover $f|_{\Omega^{+}}$ by formula (44).

##  5\. Numerical results

To illustrate the numerical feasibility of the proposed method and its
extensibility to general settings, in this section we present the results of
two numerical experiments. The rigorous analysis on the numerical methods
employed below requires further study and is left for a separate discussion.

The domain $\Omega$ is the unit disk, the measurement boundary $\Lambda$ is
the upper semicircle, and $\Omega^{+}$ is the upper semidisk. The numerical
experiments consider the boundary value problem (28) with the attenuation

$a(z)=\mu_{\text{s}}(z)+\mu_{\text{a}}(z),$

where $\mu_{\text{s}}$, and $\mu_{\text{a}}$ are the scattering and absorption
coefficients, respectively.

In the first numerical experiment $a\in C^{2}(\overline{\Omega})$ and the
scattering kernel is homogeneous with

(45)
$k(\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime})=\dfrac{\mu_{\text{s}}}{2\pi}\left[1+2g(\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime})+2g^{2}\cos(2\arccos(\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime}))\right].$

In contrast, the second numerical experiment uses a discontinuous absorption
coefficient $\mu_{a}$, and a homogeneous scattering kernel $k$ of infinite
Fourier content. More precisely, we work with the two dimensional Henyey-
Greenstein (Poisson) kernel

(46)
$k(\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime})=\mu_{\text{s}}\dfrac{1}{2\pi}\dfrac{1-g^{2}}{1-2g\boldsymbol{\theta}\cdot\boldsymbol{\theta}^{\prime}+g^{2}}.$

Note that (45) comes from the quadratic truncation of (46).

Throughout this section, $\mu_{\text{s}}\equiv 5$ and $g=1/2$ are used, and
computations are processed in the standard double precision arithmetic. The
parameter $g=1/2$ yields an anisotropic scattering half way between the
ballistic ($g=0$) regime and an isotropic scattering regime ($g=1$). The value
of the parameter $\mu_{\text{s}}$ yields that, on average, a particle scatters
after running straight every $1/5$ of the unit pathlength. These parameter
choices are meaningful in certain optical regimes, where the diffusion
approximation would not hold.

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

$\displaystyle B_{1}$
$\displaystyle=\\{(x,y)\>:\>(x-0.5)^{2}+y^{2}<0.3^{2}\\},$ $\displaystyle
B_{2}$
$\displaystyle=\left\\{(x,y)\>:\>(x+0.25)^{2}+\left(y-\dfrac{\sqrt{3}}{4}\right)<0.2^{2}\right\\},\quad\text{and}$
$\displaystyle B_{3}$
$\displaystyle=\\{(x,y)\>:\>x^{2}+(y+0.6)^{2}<0.3^{2}\\}$

be circular regions inside $\Omega$ as illustrated in Figure 1.

$\Omega$$B_{3}$$R$$B_{1}$$B_{2}$ Figure 1. Domain $\Omega$ and inclusions; the
dotted circles ($B_{1}$ and $B_{2}$) indicates highly absorbing regions, while
gray regions ($B_{2}$, $R$ and $B_{3}$) are support of internal sources.

The source term

$f(z)=\begin{cases}2,\quad&\text{in $R$};\\\ 1,\quad&\text{in $B_{2}\cup
B_{3}$};\\\ 0,\quad&\text{otherwise},\end{cases}$

used to generate the boundary data on $\Lambda_{+}$ is to be reconstructed in
the upper semidisk $\Omega^{+}$.

In the first example the absorption coefficient is a $C^{2}$-smoothen version
of the discontinuous case in (47). Namely, for $\epsilon=0.05$,

$\mu_{\text{a}}(z)=\begin{cases}2,\quad&\text{in
$\\{z\>:\>\mbox{dist}(z,\partial B_{1})\geq\epsilon\\}\cap B_{1}$};\\\
1,\quad&\text{in $\\{z\>:\>\mbox{dist}(z,\partial B_{2})\geq\epsilon\\}\cap
B_{2}$};\\\ 0.1,\quad&\text{in
$\\{z\>:\>\mbox{dist}(z,B_{1})\geq\epsilon\\}\cap\\{z\>:\>\mbox{dist}(z,B_{2})\geq\epsilon\\}$},\end{cases}$

while in the $\epsilon$-neighborhoods of $\partial B_{1}$ and $\partial B_{2}$
we use some translations and scalings of the quintic polynomial
$-(|z|-1)^{3}(6|z|^{2}+3|z|+1)$, $|z|\leq 1$, to define $\mu_{\text{a}}\in
C^{2}(\overline{\Omega})$.

For our inverse problem, the “measured” data (2) is generated by the numerical
computation of the corresponding forward problem (1) in
$\Omega\times{{\mathbf{S}}^{1}}$, where we retain the trace of the solution
$u|_{\Lambda_{+}}$ and disregard the rest. Note the contribution to the data
of radiation coming from the source supported in the lower half of the
rectangle $R$ and from the ball $B_{3}$.

The numerical solution is processed by the piecewise constant approximation
method in [15], where the spatial domain $\Omega$ and ${{\mathbf{S}}^{1}}$ are
divided into $4,823,822$ triangles (the maximum diameter is approximately
$0.0025$) and $360$ velocities at equi-intervals respectively. The boundary
$\partial\Omega$ is approximated by $5,234$ equi-length segments and, thus,
$2,617$ spatial measurement nodes are assigned on $\Lambda$.

The computed data is depicted in Figure 2. In there, for $\zeta\in\Lambda$
(indicated by $\times$), the red closed curves represent
$\zeta+2u(\zeta,\boldsymbol{\theta})$. The zero incoming flux boundary
condition can be observed on the computed radiation, where the curves do not
enter $\Omega$.

Figure 2. Boundary measurement $u|_{\Lambda_{+}}$; (left) the measurement arc
$\Lambda$ lies in the upper half plane. Each red closed curve stands for
$\zeta+2u(\zeta,\boldsymbol{\theta})$ for
$\boldsymbol{\theta}\in{{\mathbf{S}}^{1}}$ at $\zeta\in\Lambda$ indicated by
the cross symbol ($\times$); (right) magnification of
$u(\zeta,\boldsymbol{\theta})$ at $\zeta=(0,1)$.

In our numerical reconstruction, the domain of interest $\Omega^{+}$ is
partitioned into a triangular mesh without any prior information on $R$ and
$B_{i}$, $i=1,2,3$. The number of triangles ($8,631$) in this mesh is much
less than that the number of triangles ($4,823,822$) used in the computation
in the forward problem, thus avoiding an inverse crime. The attenuation
coefficient $a$ is assumed known in $\overline{\Omega}$. The triangular mesh
induces $157$ nodes on $\Lambda$ and $100$ nodes on $L$.

The Hilbert transform in (25) is computed by a method proposed in [17]. All
the integrations are approximated by the composite mid-point rule with equi-
spaced intervals. The integrating factors in (24) are computed with $100$
subintervals and $360$ velocities, and the integral equations (39) and (43)
are computed with $1,666$ subintervals (so that $\Delta x=2/1666\approx
0.0012$ on $L$ is about the same with the length $\pi/2617$ of the partition
on the arc $\Lambda$). Of particular interest, and key to our procedure, is
the numerical computation of the integral equations (39) and (43), which is
done via the collocation method with the numerical integration rule

$[I-\mathrm{i}H_{t}]u(x_{i})\approx
u(x_{i})-\dfrac{\mathrm{i}}{\pi}\sum_{j\neq
i}\dfrac{u(x_{j})}{x_{i}-x_{j}}\Delta x.$

Except for the implicit regularization due to the discretization, no other
regularization method is explicitly employed in numerical implementation. The
numerical reconstruction takes approximately $115$ seconds by OpenMP parallel
computation on two Xeon E5-2650 v4 (2.20GHz, 12 cores) processors.

Figure 3 depicts numerical reconstruction results: on the left is the profile
of reconstructed $f(z)$, and on the right is its section on the dotted line
($y=-\sqrt{3}x$). In the figure the support of $f(z)$ is successfully
reproduced. The reconstructed source $f$ away from the segment $L$ is
quantitatively in fine agreement with the exact one. However the accuracy of
the reconstruction decreases at points on the support of $f$ which are close
to the line segment $L$. One of the reasons is the ill-posedness in equation
(39), which we currently mitigate solely via the discretization. For a more
accurate numerical reconstruction at such points, the choice of a better
regularizer is needed.

Figure 3. Numerically reconstructed $f(z)$ for polynomial type scattering
kernel (45) and $\mu_{\text{a}}\in C^{2}(\Omega)$. (left) the profile of
reconstructed $f$ in the domain of interest $\Omega^{+}$; (right) the section
on the dotted line.

In the second example, the proposed algorithm is applied to the non-polynomial
type scattering kernel (46) and the discontinuous absorption coefficient

(47) $\displaystyle\mu_{\text{a}}(z)=\begin{cases}2,\quad&\text{in
$B_{1}$};\\\ 1,\quad&\text{in $B_{2}$};\\\
0.1,\quad&\text{otherwise}.\end{cases}$

Figure 4. Numerically reconstructed $f(z)$ for the case with Poisson kernel
(non-polynomial type scattering kernel) and discontinuous $\mu_{\text{a}}$.

Figure 4 shows numerical results, where (40) is recursively solved 10 times
($M=10$). Even though the smoothness hypotheses in Theorem 4.1 are violated,
and despite of a scattering kernel of infinite Fourier content, the numerical
results are in agreement with the results in the first example (where all the
hypotheses in Theorem 4.1 were satisfied). This illustrates the robustness of
the proposed reconstruction method, and indicates that the reconstruction
result in Theorem 4.1 may hold under more relaxed hypotheses.

## Acknowledgment

The work of H. Fujiwara was supported by JSPS KAKENHI Grant Numbers JP19H00641
and JP20H01821. The work of K. Sadiq was supported by the Austrian Science
Fund (FWF), Project P31053-N32. The work of A. Tamasan was supported in part
by the NSF grant DMS-1907097.

## References

  * [1] D. S. Anikonov, A. E. Kovtanyuk and I. V. Prokhorov, Transport Equation and Tomography, Inverse and Ill-Posed Problems Series, 30, Brill Academics, (2002).
  * [2] E. V. Arbuzov, A. L. Bukhgeim and S. G. Kazantsev, Two-dimensional tomography problems and the theory of A-analytic functions, Siberian Adv. Math., 8 (1998), 1–20.
  * [3] S. R. Arridge, Optical tomography in medical imaging, Inverse Problems, 15 (1999), R41–R93.
  * [4] G. Bal, On the attenuated Radon transform with full and partial measurements, Inverse Problems 20 (2004), 399–418.
  * [5] G. Bal and A. Tamasan, Inverse source problems in transport equations, SIAM J. Math. Anal., 39(2007), 57–76.
  * [6] J. Boman and J.-O. Strömberg, Novikov’s inversion formula for the attenuated Radon transform–a new approach, J. Geom. Anal., 14 (2004), 185–198.
  * [7] A. L. Bukhgeim, Inversion Formulas in Inverse Problems, chapter in Linear Operators and Ill-Posed Problems by M. M. Lavrentiev and L. Ya. Savalev, Plenum, New York, 1995.
  * [8] S. Chandrasekhar, Radiative Transfer, Dover Publ., New York, 1960.
  * [9] M. Choulli and P. Stefanov, Inverse scattering and inverse boundary value problems for the linear Boltzmann equation, Comm. Partial Differential Equations, 21(5–6): (1996), 763–785.
  * [10] M. Choulli and P. Stefanov, An inverse boundary value problem for the stationary transport equation, Osaka J. Math., 36 (1999), 87–104.
  * [11] R. Dautray and J-L. Lions, Mathematical analysis and numerical methods for science and technology, Vol. 4., Springer-Verlag, Berlin, 1990.
  * [12] N. Dunford and J. Schwartz, Linear Operators, Part I: General Theory, Wiley, Interscience Publ., New York, 1957.
  * [13] D. V. Finch, The attenuated x-ray transform: recent developments, in Inside out: inverse problems and applications, Math. Sci. Res. Inst. Publ., 47, Cambridge Univ. Press, Cambridge, 2003, 47–66.
  * [14] L. Florescu, V. Markel, and J. Schotand, Nonreciprocal Broken-Ray Transforms with Applications to Fluorescence Imaging, Inverse Probl. 34(9), (2018), 12 pp.
  * [15] H. Fujiwara, Piecewise constant upwind approximations to the stationary radiative transport equation, in: Mathematical Analysis of Continuum Mechanics and Industrial Applications III, CoMFoS 2018, Mathematics for Industry, 34, Springer, 2020, 35–45.
  * [16] H. Fujiwara, K. Sadiq and A. Tamasan, A Fourier approach to the inverse source problem in an absorbing and anisotropic scattering medium, Inverse Problems 36(1):015005 (2019).
  * [17] H. Fujiwara, K. Sadiq and A. Tamasan, Numerical reconstruction of radiative sources in an absorbing and non-diffusing scattering medium in two dimensions, SIAM J. Imaging Sci., 13(1) (2020), 535–555.
  * [18] H. Fujiwara, K. Sadiq and A. Tamasan, Partial inversion of the 2D attenuated $X$-ray transform with data on an arc, Inverse Probl. Imaging (2021), to appear.
  * [19] W. Koppelman and J. D. Pincus, Spectral representations for finite Hilbert transformations, Math. Z., 71 (1959), 399–407.
  * [20] E. W. Larsen, The inverse source problem in radiative transfer, J. Quant. Spect. Radiat. Transfer, 15 (1975), 1–5.
  * [21] S. G. Mikhlin and S. Prössdorf, Singular integral operators, Springer-Verlag, Berlin, 1986.
  * [22] M. Mokhtar-Kharroubi, Mathematical topics in neutron transport theory, World Scientific, Singapore, 1997.
  * [23] F. Monard and G. Bal, Inverse source problems in transport via attenuated tensor tomography, arXiv:1908.06508v1.
  * [24] N. I. Muskhelishvili, Singular Integral Equations, Dover, New York, 2008.
  * [25] F. Natterer, The mathematics of computerized tomography, Wiley, New York, 1986.
  * [26] F. Natterer, Inversion of the attenuated Radon transform, Inverse Problems 17 (2001), 113–119.
  * [27] R. G. Novikov, Une formule d’inversion pour la transformation d’un rayonnement X atténué, C. R. Acad. Sci. Paris Sér. I Math., 332 (2001), 1059–1063.
  * [28] R. G. Novikov, On the range characterization for the two-dimensional attenuated x-ray transformation, Inverse Problems 18 (2002), no. 3, 677–700.
  * [29] S. Okada, and D. Elliott, The finite Hilbert transform in $L^{2}$, Math. Nachr., 153 (1991), 43–56.
  * [30] M. Reed and B. Simon, Methods of modern mathematical physics, I: Functional analysis, 2nd ed., Academic Press, New York, 1980.
  * [31] K. Sadiq and A. Tamasan, On the range of the attenuated Radon transform in strictly convex sets, Trans. Amer. Math. Soc., 367(8) (2015), 5375–5398.
  * [32] K. Sadiq and A. Tamasan, On the range characterization of the two dimensional attenuated Doppler transform, SIAM J. Math. Anal., 47(3) (2015), 2001–2021.
  * [33] K. Sadiq, O. Scherzer, and A. Tamasan, On the X-ray transform of planar symmetric 2-tensors, J. Math. Anal. Appl., 442(1) (2016), 31–49.
  * [34] V. A. Sharafutdinov, Inverse problem of determining a source in the stationary transport equation on a Riemannian manifold, Mat. Vopr. Teor. Rasprostr. Voln., 26 (1997), 236–242; translation in J. Math. Sci. (New York) 96 (1999), no. 4, 3430–3433.
  * [35] A. V. Smirnov, M. V. Klibanov, and L. H. Nguyen, On an inverse source problem for the full radiative transfer equation with incomplete data, SIAM J. Sci. Comput., 41(5) (2019), B929–B952.
  * [36] P. Stefanov and G. Uhlmann,  An inverse source problem in optical molecular imaging, Anal. PDE, 1(1) (2008), 115–126.
  * [37] A. Tamasan, An inverse boundary value problem in two-dimensional transport, Inverse Problems 18 (2002), 209–219.
  * [38] A. Tamasan, Optical tomography in weakly anisotropic scattering media, Contemporary Mathematics 333 (2003), 199–207.
  * [39] F. G. Tricomi, Integral Equations, Interscience, New York, 1957.
  * [40] I. N. Vekua, Generalized Analytic Functions, Pergamon Press Ltd. 1962.
  * [41] H. Widom, Singular integral equations in $L_{p}$, Trans. Amer. Math. Soc., 97 (1960), 131–160.