# Uniqueness and numerical inversion in the time-domain fluorescence diffuse
optical tomography

Chunlong Sun<EMAIL_ADDRESS>College of Science, Nanjing University of
Aeronautics and Astronautics, Nanjing 211106, Jiangsu, China  Nanjing Center
for Applied Mathematics, Nanjing 211135, Jiangsu, China  Zhidong Zhang
<EMAIL_ADDRESS>School of Mathematics (Zhuhai), Sun Yat-sen
University, Zhuhai 519082, Guangdong, China

###### Abstract

This work considers the time-domain fluorescence diffuse optical tomography
(FDOT). We recover the distribution of fluorophores in biological tissue by
the boundary measurements. With the Laplace transform and the knowledge of
complex analysis, we build the uniqueness theorem of this inverse problem.
After that, the numerical reconstructions are considered. We introduce a non-
iterative inversion strategy by peak detection and an iterative inversion
algorithm under the framework of regularizing scheme, then give several
numerical examples in three-dimensional space illustrating the performance of
the proposed inversion schemes.

Keywords: inverse source problem, uniqueness, diffusion equation, numerical
inversion.

AMS Subject Classifications: 35R30, 65M32.

##  1 Introduction.

###  1.1 Background and literature.

Recently, fluorescence diffuse optical tomography (FDOT) is rapidly gaining
acceptance as an important diagnostic and monitoring tool of symptoms in
medical applications [28, 30, 32]. It aims to recover the fluorophores
quantitatively from some measurements specified on the medium surface.
Compared to X-ray computed tomography (CT), magnetic resonance imaging (MRI)
and positron emission tomography (PET), FDOT possesses some advantages like
low cost and portable, and therefore it has drawn more and more attention.

For the FDOT, two physical processes are coupled, namely, excitation and
fluorescence (emission), which are described by the photon transport model.
The system of Maxwell’s equations is the rigorous photon transport model to
describe the light propagation. However, since the wave property of the photon
is lost by multiple scattering in strongly turbid medium such as biological
tissue, the photon propagation through turbid media for both excitation and
emission can be governed by the Boltzmann radiative transfer equation or
shortly the radiative transfer equation (RTE) [14, 21, 26]. This equation
describes the propagation of photon in scattering media and is given as:

$\displaystyle c^{-1}\partial_{t}\Phi(x,\theta,t)=$
$\displaystyle-\theta\cdot\nabla\Phi(x,\theta,t)-\big{(}\mu_{a}+\mu_{f}+\mu_{s}\big{)}\Phi(x,\theta,t)$
(1.1)
$\displaystyle+\mu_{s}\int_{\mathbb{S}^{2}}\eta(\theta,\theta^{\prime})\Phi(x,\theta^{\prime},t)d\theta^{\prime}+S(x,\theta,t),$

where $\Phi(x,\theta,t)$ is the radiant intensity at $x$ in time $t$, and
$\theta$ is a unit vector pointing in the direction of interest. The
absorption and scattering coefficients of medium, denoted by $\mu_{a}(x)$ and
$\mu_{s}(x)$, are the inverses of the mean free paths for absorption and
scattering, respectively. Here $\mu_{f}(x)$ denotes the absorption coefficient
of the fluorophores, and $c$ is the speed of photon inside the medium. The
scattering phase function $\eta(\theta,\theta^{\prime})$, which describes the
probability that a photon traveling in direction $\theta$ is scattered within
the unit solid angle around the direction $\theta^{\prime}$, is given by
3-dimensional Henyey-Greenstein function [26].

However, the cost of solving RTE (1.1) is extremely high, and it is shown that
the diffusion equation (DE) can be a sufficiently accurate approximation to
the RTE but has much lower computational cost. Therefore people usually
replace the RTE by the diffusion equation. The basic idea of the diffusion
approximation is: when scattering is much stronger than absorption, the
radiant intensity $\Phi$ can be expressed as an isotropic photon density plus
a small photon flux, and sequentially the transport equation (1.1) can be
reduced to a diffusion equation (see Section (2.1)). We consider the diffusion
approximation for both excitation and emission, and give the propagation of
excitation light and emission (fluorescence) light as

$\begin{cases}\begin{aligned}
\left(c^{-1}\partial_{t}-\mu_{D}\Delta+\mu_{a}\right)u_{e}&=0,&&(x,t)\in\Omega\times(0,\infty),\\\
u_{e}(x,0)&=0,&&x\in\Omega,\\\ \frac{\partial u_{e}}{\partial\nu}+\beta
u_{e}&=\delta(x-x_{s})\delta(t),&&(x,x_{s},t)\in\partial\Omega\times\partial\Omega\times(0,\infty),\end{aligned}\end{cases}$
(1.2)

and

$\begin{cases}\begin{aligned}
\left(c^{-1}\partial_{t}-\mu_{D}\Delta+\mu_{a}\right)u_{m}&=S_{0}[\mu_{f},u_{e}](x,t;x_{s}),&&(x,t)\in\Omega\times(0,\infty),\\\
u_{m}(x,0)&=0,&&x\in\Omega,\\\ \frac{\partial u_{m}}{\partial\nu}+\beta
u_{m}&=0,&&(x,t)\in\partial\Omega\times(0,\infty),\end{aligned}\end{cases}$
(1.3)

where $u_{e}(x,t)$ and $u_{m}(x,t)$ are the photon densities of excitation
light and emission light, respectively. Here $\Omega\subset\mathbb{R}^{d}$ is
an open bounded subset with sufficiently smooth boundary; the vanishing
initial condition and Robin impedance boundary condition with $\beta>0$ are
used; $\mu_{D}$ is the photon diffusion coefficient defined by
$\mu_{D}:=1/(3\mu_{s}(1-g))$ where $g$ is the anisotropy parameter. We set
$\mu_{D}$ and $\mu_{a}$ as positive constants in this work. The source term
$S_{0}$ for $u_{m}$ on the right-hand side of (1.3) contains the excitation
$u_{e}$ and is specified by

$S_{0}[\mu_{f},u_{e}](x,t;x_{s})=\frac{\mu_{f}(x)}{\tau}\int_{0}^{t}e^{-\frac{t-s}{\tau}}u_{e}(x,s;x_{s})\,ds,$
(1.4)

where $\tau>0$ is the fluorescence lifetime.

The FDOT based on RTE model or DE model has drawn extensive attention of
researchers in recent years. For such inverse problems, three kinds of
measurements, like time domain [18, 22], continuous wave [13, 25, 29] and
frequency domain [10, 27] are used. For these optical imaging problems, the
necessary regularization techniques such as Tikhonov regularization, sparse
regularization methods and hybrid regularization methods have been introduced
to overcome the ill-posedness of the inverse problems [11, 15, 25].
Especially, the sparse regularization methods have additional advantages in
promoting sparsity and higher spatial resolution for the cases that the target
is relatively small compared to the background. Further, aiming to the best
reconstruction result, one must consider how to place the source-detector
pairs onto the object surface. To answer this problem, the definition of
admissible set of possible sensor patterns on the boundary and a proper
optimal condition are usually needed [7, 12, 16, 20].

Above all, most existing works focus on improving the quality of the
reconstruction, especially the spatial resolution and the reconstruction
speed. However, in the aspect of theorical analysis, the uniqueness of FDOT
has not been rigorously investigated as far as we know. In this work, we will
study the uniqueness of time-domain FDOT. The details are contained in the
main theorem, which is stated in the next subsection.

###  1.2 The main theorem.

In this work, we consider the time-dependent fluorophores which means that the
absorption coefficient $\mu_{f}$ depends on both $x$ and $t$. We suppose that
$S_{0}[\mu_{f},u_{e}](x,t;x_{s})$ in (1.4) possesses the following semi-
discrete formulation:

$S_{0}[\mu_{f},u_{e}]=\sum_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}},$
(1.5)

where the time mesh $\\{t_{k}\\}_{k=0}^{K}$ is given and the spatial
components $\\{p_{k}(x)\\}_{k=1}^{K}$ are undetermined. For the mesh
$\\{t_{k}\\}_{k=0}^{K}$, $K$ can be infinity and
$\inf\\{|t_{k}-t_{k+1}|:k=0,\cdots,K-1\\}>0$; for the unknown
$\\{p_{k}(x)\\}_{k=1}^{K}$, we consider to recover them in $L^{2}(\Omega)$,
i.e. set $\\{p_{k}(x)\\}_{k=1}^{K}\subset L^{2}(\Omega)$. The representation
(1.5) contains the information of the time evolution process of the unknown
fluorophores. Since the excitation $u_{e}$ is known, we can extract the
information of $\mu_{f}$ from the components $\\{p_{k}(x)\\}_{k=1}^{K}$. The
measurements we used are the boundary data given by

$u_{m}(x_{d},t;x_{s}),\quad x_{d}\in\Gamma_{d}\subset\partial\Omega,\
x_{s}\in\Gamma_{s}\subset\partial\Omega,\ t\in(0,\infty),$ (1.6)

where $\Gamma_{s}$ and $\Gamma_{d}$ are the excitation area and observed area,
respectively. From the Robin boundary condition, we can get the flux
$\frac{\partial
u_{m}}{\partial\nu}\big{|}_{\Gamma_{d}\times(0,\infty)\times\Gamma_{s}}$ from
the measurements (1.6). Hence, in this work, the so-called time-domain FDOT
based on the initial-boundary value problems (1.2)-(1.3) is to solve the
inverse problem:

$\text{recovering}\ \\{p_{k}(x)\\}_{k=1}^{K}\ \text{in (\ref{ua-time}) from
the measurements (\ref{Um-measure}) or the flux}\ \frac{\partial
u_{m}}{\partial\nu}\big{|}_{\Gamma_{d}\times(0,\infty)\times\Gamma_{s}}.$

Moreover, if either of $\mu_{D}$ and $\mu_{a}$ is also unknown and we want to
recover the unknown coefficient and the source $\mu_{f}$ simultaneously, an
extra condition needs to be set. Letting $(\lambda_{1},\varphi(x))$ be one
principal eigenpair of the operator $-\Delta$ on $\Omega$ with Robin boundary
condition, we give the following condition:

$\langle p_{1}(\cdot),\varphi(\cdot)\rangle_{L^{2}(\Omega)}\neq 0.$ (1.7)

Here the notation $\langle\cdot,\cdot\rangle_{L^{2}(\Omega)}$ means the inner
product in space $L^{2}(\Omega)$. For the details of the eigensystem of
$-\Delta$ on $\Omega$, see Section (2.2.1). Also we give Assumption (2.4) in
Section (2.2.1), which is referred in the statement of the uniqueness theorem.
Furthermore, to make our analysis more convenient, we set $c=1$ in models
(1.2)-(1.3). Now it is time to state the uniqueness theorem.

######  Theorem 1.

Set Assumption (2.4) be valid, and $\Gamma_{s}$ and $\Gamma_{d}$ be a single
point and a nonempty open subset of boundary $\partial\Omega$, respectively.
For the two sets of unknowns
$\big{\\{}\mu_{D},\mu_{a},\\{p_{k}(x)\\}_{k=1}^{K}\big{\\}}$ and
$\big{\\{}\tilde{\mu}_{D},\tilde{\mu}_{a},\\{\tilde{p}_{k}(x)\\}_{k=1}^{K}\big{\\}}$,
we denote the corresponding solutions of equation (1.3) by $u$ and $\tilde{u}$
respectively, and are provided

$\frac{\partial
u}{\partial\nu}(x_{d},t;x_{s})=\frac{\partial\tilde{u}}{\partial\nu}(x_{d},t;x_{s}),\quad(x_{d},x_{s},t)\in\Gamma_{d}\times\Gamma_{s}\times(0,\infty).$

Then we conclude that:

  * (i)

given $\mu_{a}=\tilde{\mu}_{a}$ and $\mu_{D}=\tilde{\mu}_{D}$, we have

$p_{k}(x)=\tilde{p}_{k}(x)\ \text{in}\ L^{2}(\Omega)\ \text{for}\
k=1,\cdots,K;$

  * (ii)

with condition (1.7) and $\mu_{D}=\tilde{\mu}_{D}$, we have

$\mu_{a}=\tilde{\mu}_{a},\ p_{k}(x)=\tilde{p}_{k}(x)\ \text{in}\
L^{2}(\Omega)\ \text{for}\ k=1,\cdots,K;$

  * (iii)

with condition (1.7) and $\mu_{a}=\tilde{\mu}_{a}$, we have

$\mu_{D}=\tilde{\mu}_{D},\ p_{k}(x)=\tilde{p}_{k}(x)\ \text{in}\
L^{2}(\Omega)\ \text{for}\ k=1,\cdots,K.$

From Theorem (1), we can derive the following corollary straightforwardly.

######  Corollary 1.1.

Under the conditions of Theorem (1), the boundary data (1.6) can identify the
stationary fluorophores uniquely.

The above corollary considers the case of stationary fluorescence target, i.e.
$\mu_{f}:=\mu_{f}(x)$. This is the common problem in the community of
fluorescence tomography, which is covered by Theorem (1).

###  1.3 Contribution and outline.

The literature on the uniqueness of time-domain FDOT is relatively rare and
the existing conclusions seem not to be convenient in practical applications.
Here we discuss [24] as an example. In this article, the authors aim to
identify the absorption coefficient $\mu_{f}(x)$ from time-resolved boundary
measurements. But this work requires strong prior assumptions on $\mu_{f}(x)$.
The authors prove that for $\Omega=\mathbb{R}_{+}^{3}$,
$x=(\tilde{x},x_{3})\in\mathbb{R}^{2}\times\mathbb{R}_{+}^{1}$ and small
lifetime $\tau$, by supposing $\mu_{f}(x)$ has the variable separable form
$\mu_{f}(x)=p(\tilde{x})q(x_{3})$ with known vertical information $q(x_{3})$,
the horizontal information $p(\tilde{x})$ can be uniquely determined from the
boundary measurements (1.6) with
$(x_{d},x_{s})\in\\{x_{d}^{0}\\}\times\partial\Omega$ or
$(x_{d},x_{s})\in\partial\Omega\times\\{x_{s}^{0}\\}$. Here $x_{d}^{0}$ and
$x_{s}^{0}$ are the single points on $\partial\Omega$. However, we see that
the used measurements require either the source point $x_{s}$ or the detector
point $x_{d}$ to cover the whole boundary surface, which implies the
impractical cost in applications.

In this work, the uniqueness of time-domain FDOT in a more general framework
is established. The concerned unknown $\mu_{f}$ has a more general
formulation, and we use the sparse boundary measurements (we call the used
measurements as sparse boundary measurements since the observed area
$\Gamma_{d}$ can be an arbitrarily open subset of boundary $\partial\Omega$
and the excitation area $\Gamma_{s}$ only needs to be a single point). Theorem
(1) confirms that the semi-discrete unknown $\mu_{f}$ can be uniquely
determined by the sparse boundary measurements. Therefore, this conclusion is
of great significance in practical applications. Moreover, Theorem (1) also
concerns the cases of recovering $\mu_{f}$ and one of $\mu_{a}$ or $\mu_{D}$
simultaneously. The reconstructions of $\mu_{a}$ and $\mu_{D}$ are known as
the diffuse optical tomography (DOT). Some aspects of uniqueness and ill-
posedness of DOT are considered [1, 2, 4, 6, 9, 19, 31]. Theorem (1) gives
partial uniqueness result for the DOT with fluorescence but it requires
constant absorption and scattering. Hence the main theorem in this work is
also useful in DOT to some extent.

The rest of this article is organized as follows. In Section (2), we first
introduce the semi-discretized model for our time-domain FDOT based on
diffusion equations. Then we collect several preliminary works. We prove the
uniqueness theorem in Section (3). The Laplace transform and some knowledge of
complex analysis will be used. The numerical inversions will be considered in
Section (4). The inversion scheme minimizing the regularized cost functional
is implemented by an iterative process. We show the validity of the proposed
scheme by several numerical examples in three-dimensional space.

##  2 The diffusion approximation and some preliminary results.

###  2.1 The diffusion approximation.

We refer to [3, 5] for the details of diffusion approximation (or
$P_{1}$-approximation) and [24] for the error estimates of model
approximation. However, we give a sketch for the derivation from RTE to DE for
the completeness of this article.

First, one can expand $\Phi(x,\theta,t)$ and $S(x,\theta,t)$ in RTE (1.1) into
spherical harmonics as

$\displaystyle\Phi(x,\theta,t)$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left(\frac{2l+1}{4\pi}\right)^{1/2}\Phi_{l,m}(x,t)Y_{l,m}(\theta),$
$\displaystyle S(x,\theta,t)$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left(\frac{2l+1}{4\pi}\right)^{1/2}S_{l,m}(x,t)Y_{l,m}(\theta).$

Keeping only the terms $l\leq 1$ in above expansions, we have the
approximation

$\Phi(x,\theta,t)\approx\frac{1}{4\pi}u(x,t)+\frac{3}{4\pi}\theta\cdot{\mathbf{J}}(x,t),\quad\theta\in\mathbb{S}^{2},$
(2.1)

where the quantity

$u(x,t):=\int_{\mathbb{S}^{2}}\Phi(x,\theta,t)\,d\theta$

is called photon density, and the quantity

$\mathbf{J}(x,t):=\int_{\mathbb{S}^{2}}\theta\cdot\Phi(x,\theta,t)\,d\theta$

is called photon flux. Similarly, we have

$S(x,\theta,t)\approx\frac{1}{4\pi}S_{0}(x,t)+\frac{3}{4\pi}\theta\cdot
S_{1}(x,t),\quad\theta\in\mathbb{S}^{2},$ (2.2)

where $S_{0}(x,t):=\int_{\mathbb{S}^{2}}S(x,\theta,t)\,d\theta$ and
$S_{1}(x,t):=\int_{\mathbb{S}^{2}}\theta\cdot S(x,\theta,t)\,d\theta$ are the
isotropic component and the first angular moment of source term, respectively.

Next, by owing (2.1) and (2.2) into the RTE (1.1), and then integrating the
RTE and the RTE multiplied by $\theta$ over ${\mathbb{S}^{2}}$, we have

$\displaystyle\left(c^{-1}\partial_{t}+\mu_{a}+\mu_{f}\right)u(x,t)+\nabla\cdot\mathbf{J}(x,t)$
$\displaystyle=S_{0},$ (2.3)
$\displaystyle\left(c^{-1}\partial_{t}+\mu_{a}+\mu_{f}+\mu_{s}^{\prime}\right)\mathbf{J}(x,t)+\frac{1}{3}\nabla
u(x,t)$ $\displaystyle=S_{1},$

where $\mu_{s}^{\prime}$ is the reduced scattering coefficient defined by
$\mu_{s}^{\prime}:=\mu_{s}(1-g)$ with $g$ as the anisotropy parameter. By
$\partial|\mathbf{J}|/\partial t\ll c\mu_{s}^{\prime}|\mathbf{J}|$ and the
assumption that sources are isotropic, we arrive at the so-called Fick’s law

$\mathbf{J}=-\frac{1}{3(\mu_{a}+\mu_{f}+\mu_{s}^{\prime})}\nabla u.$

Since $\mu_{a},\,\mu_{f}\ll\mu_{s}^{\prime}$ in the strong scattering medium,
we can approximate

$\mathbf{J}=-\frac{1}{3\mu_{s}^{\prime}}\nabla u.$ (2.4)

Finally, supposing the absorption of the fluorophores is much smaller than the
absorption of the medium, from substituting (2.4) into (2.3) we immediately
have the diffusion equation given by

$\left(c^{-1}\partial_{t}-\nabla\cdot\frac{1}{3\mu_{s}^{\prime}(x)}\nabla+\mu_{a}(x)\right)u=S_{0}(x,t).$

We deduce the above diffusion approximation for both excitation and emission,
and suppose that the optical parameters for emission are constants and same to
those for the excitation. Then the system (1.2)-(1.3) can be derived.

###  2.2 Some preliminary results.

####  2.2.1 The eigensystem of $-\Delta$ on $\Omega$.

For the operator $-\Delta$ on $H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$ with Robin
boundary condition, we denote the eigensystem by
$\\{\lambda_{n},\varphi_{n}(x)\\}_{n=1}^{\infty}$. Then the following
properties will be valid:

  * •

$0<\lambda_{1}\leq\lambda_{2}\leq\cdots$ and $\lambda_{n}\to\infty$ as
$n\to\infty$;

  * •

$\\{\varphi_{n}(x)\\}_{n=1}^{\infty}$ is an orthonormal basis of
$L^{2}(\Omega)$.

Furthermore, if $\varphi_{n}$ is an eigenfunction of $-\Delta$ corresponding
to $\lambda_{n}$, so is $\overline{\varphi_{n}}$, where
$\overline{\varphi_{n}}$ is the complex conjugate of $\varphi_{n}$. Hence we
have that the set $\\{\varphi_{n}(x)\\}_{n=1}^{\infty}$ coincides with
$\\{\overline{\varphi_{n}(x)}\\}_{n=1}^{\infty}$. The trace theorem yields
that
$\\{\frac{\partial\varphi_{n}}{\partial\nu}|_{\partial\Omega}\\}_{n=1}^{\infty}\subset
H^{1/2}(\partial\Omega)$. Also, we denote
$\langle\cdot,\cdot\rangle_{L^{2}(\partial\Omega)}$ as the inner product in
$L^{2}(\partial\Omega)$.

The next lemmas concern the vanishing property and the density of
$\frac{\partial\varphi_{n}}{\partial\nu}$ on $\partial\Omega$.

######  Lemma 2.1.

If $\Gamma$ is a nonempty open subset of $\partial\Omega$, then for each
$n\in\mathbb{N}^{+}$, $\frac{\partial\varphi_{n}}{\partial\nu}$ can not vanish
on $\Gamma$.

###### Proof.

See [23, Lemma 2.1]. ∎

######  Lemma 2.2.

The set
$\text{Span}\\{\frac{\partial\varphi_{n}}{\partial\nu}|_{\partial\Omega}\\}_{n=1}^{\infty}$
is dense in $L^{2}(\partial\Omega)$.

###### Proof.

Not hard to see that $H^{3/2}(\partial\Omega)$ is dense in
$L^{2}(\partial\Omega)$ under the norm $\|\cdot\|_{L^{2}(\partial\Omega)}$. So
it is sufficient to show $\tilde{\psi}\in H^{3/2}(\partial\Omega)$ vanishes
almost everywhere on $\partial\Omega$ if
$\langle\tilde{\psi},\frac{\partial\varphi_{n}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}=0$
for $n\in\mathbb{N}^{+}$.

We set $\psi$ be the weak solution of the system below:

$\begin{cases}\begin{aligned} -\Delta\psi(x)&=0,&&x\in\Omega,\\\
\frac{\partial\psi}{\partial\nu}+\beta\psi&=\tilde{\psi},&&x\in\partial\Omega.\end{aligned}\end{cases}$

We have $\psi\in H^{2}(\Omega)$ from the regularity $\tilde{\psi}\in
H^{3/2}(\partial\Omega)$, sequentially the Green’s identity can be used. For
$n\in\mathbb{N}^{+}$, we have

$\displaystyle\langle-\Delta\psi,\varphi_{n}\rangle_{L^{2}(\Omega)}-\langle\psi,-\Delta\varphi_{n}\rangle_{L^{2}(\Omega)}$
$\displaystyle=\langle\psi,\frac{\partial\varphi_{n}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}-\langle\frac{\partial\psi}{\partial\nu},\varphi_{n}\rangle_{L^{2}(\partial\Omega)}$
$\displaystyle=\langle\psi+\beta^{-1}\frac{\partial\psi}{\partial\nu},\frac{\partial\varphi_{n}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}$
$\displaystyle=\beta^{-1}\langle\tilde{\psi},\frac{\partial\varphi_{n}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}.$

From $\Delta\psi=0$ on $\Omega$ and the fact
$\langle\tilde{\psi},\frac{\partial\varphi_{n}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}=0$,
we have

$\langle\psi,-\Delta\varphi_{n}\rangle_{L^{2}(\Omega)}=\lambda_{n}\langle\psi,\varphi_{n}\rangle_{L^{2}(\Omega)}=0.$

So we have proved that for each $n\in\mathbb{N}^{+}$,
$\langle\psi,\varphi_{n}\rangle_{L^{2}(\Omega)}=0$. Recalling the completeness
of $\\{\varphi_{n}\\}_{n=1}^{\infty}$ in $L^{2}(\Omega)$, it holds that
$\|\psi\|_{L^{2}(\Omega)}=0$. From the definition of weak derivative and
Sobolev space, we have $\|\psi\|_{H^{2}(\Omega)}=0$. By the continuity of the
trace operator, it gives that

$\|\psi\|_{L^{2}(\partial\Omega)}\leq
C\|\psi\|_{H^{2}(\Omega)}=0,\quad\Big{\|}\frac{\partial\psi}{\partial\nu}\Big{\|}_{L^{2}(\partial\Omega)}\leq
C\|\psi\|_{H^{2}(\Omega)}=0.$

This means that $\tilde{\psi}=0$ almost everywhere on $\partial\Omega$ and the
proof is complete. ∎

####  2.2.2 The set $\\{\xi_{l}\\}_{l=1}^{\infty}$ and the coefficients
$\\{c_{z,n}\\}$.

From the above lemma, we are allowed to construct the orthonormal basis
$\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$ in $L^{2}(\partial\Omega)$. Firstly we
set
$\tilde{\xi}_{1}=\frac{\partial\varphi_{1}}{\partial\nu}|_{\partial\Omega}/\|\frac{\partial\varphi_{1}}{\partial\nu}\|_{L^{2}(\partial\Omega)},$
and assume that the orthonormal set $\\{\tilde{\xi}_{j}\\}_{j=1}^{l-1}$ has
been built for $l=2,3,\cdots$. Then we set $n_{l}\in\mathbb{N}^{+}$ be the
smallest number such that
$\frac{\partial\varphi_{n_{l}}}{\partial\nu}|_{\partial\Omega}\notin\text{Span}\\{\tilde{\xi}_{j}\\}_{j=1}^{l-1}$,
and pick
$\tilde{\xi}_{l}\in\text{Span}\\{\frac{\partial\varphi_{n_{l}}}{\partial\nu}|_{\partial\Omega},\
\tilde{\xi}_{1},\cdots,\tilde{\xi}_{l-1}\\}$ satisfying

$\langle\tilde{\xi}_{l},\tilde{\xi}_{j}\rangle_{L^{2}(\partial\Omega)}=0\
\text{for}\ j=1,\cdots,l-1,\ \text{and}\
\|\tilde{\xi}_{l}\|_{L^{2}(\partial\Omega)}=1.$

The density of
$\text{Span}\\{\frac{\partial\varphi_{n}}{\partial\nu}|_{\partial\Omega}\\}_{n=1}^{\infty}$
in $L^{2}(\partial\Omega)$ yields that $\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$
is an orthonormal basis in $L^{2}(\partial\Omega)$. Also, we have
$\tilde{\xi}_{l}\in H^{1/2}(\partial\Omega)$ for each $l\in\mathbb{N}^{+}$.

Next, for $l\in\mathbb{N}^{+}$, we define $\xi_{l}\in H^{1}(\Omega)$ be the
weak solution of the system:

$\begin{cases}\begin{aligned}
(-\mu_{D}\Delta+\mu_{a})\xi_{l}(x)&=0,&&x\in\Omega,\\\
\frac{\partial\xi_{l}}{\partial\nu}+\beta\xi_{l}&=\tilde{\xi}_{l},&&x\in\partial\Omega.\end{aligned}\end{cases}$
(2.5)

Fixing $z\in\partial\Omega$, we define the series $\psi_{z}^{N}\in
H^{1}(\Omega)$ as

$\psi_{z}^{N}(x)=\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\overline{\xi_{l}(x)},\
x\in\Omega.$ (2.6)

The definition of coefficients $\\{c_{z,n}\\}$ is given in the next lemma.

######  Lemma 2.3.

For each $z\in\partial\Omega$ and $n\in\mathbb{N}^{+}$,
$\lim_{N\to\infty}\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{L^{2}(\Omega)}$
exists and we denote the limit by $c_{z,n}$.

###### Proof.

From Green’s identities we have

$\displaystyle\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{L^{2}(\Omega)}$
$\displaystyle=(\mu_{D}\lambda_{n}+\mu_{a})^{-1}\langle\psi_{z}^{N},(-\mu_{D}\Delta+\mu_{a})\overline{\varphi_{n}}\rangle_{L^{2}(\Omega)}$
$\displaystyle=(\mu_{D}\lambda_{n}+\mu_{a})^{-1}\Big{(}\langle(-\mu_{D}\Delta+\mu_{a})\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{L^{2}(\Omega)}-\mu_{D}\langle\psi_{z}^{N},\frac{\partial\overline{\varphi_{n}}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}$
$\displaystyle\qquad\qquad\qquad\qquad\
+\mu_{D}\langle\frac{\partial\psi_{z}^{N}}{\partial\nu},\overline{\varphi_{n}}\rangle_{L^{2}(\partial\Omega)}\Big{)}$
$\displaystyle=-(\mu_{D}\lambda_{n}+\mu_{a})^{-1}\mu_{D}\beta^{-1}\langle\frac{\partial\psi_{z}^{N}}{\partial\nu}+\beta\psi_{z}^{N},\frac{\partial\overline{\varphi_{n}}}{\partial\nu}\rangle_{L^{2}(\partial\Omega)}$
$\displaystyle=-(\mu_{D}\lambda_{n}+\mu_{a})^{-1}\mu_{D}\beta^{-1}\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\langle\frac{\partial\varphi_{n}}{\partial\nu},\tilde{\xi}_{l}\rangle_{L^{2}(\partial\Omega)}=:c_{z,n}^{N},$

where the system (2.5) and the boundary condition of $\varphi_{n}$ are used.
From the definition of $\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$, we have
$\langle\frac{\partial\varphi_{n}}{\partial\nu},\tilde{\xi}_{l}\rangle_{L^{2}(\partial\Omega)}=0$
for large $l$. Hence the value of $c_{z,n}^{N}$ will not change if $N$ is
sufficiently large. This gives that
$\lim_{N\to\infty}\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{L^{2}(\Omega)}$
exists and the proof is complete. ∎

From the above lemma we can see

$c_{z,n}=\lim_{N\to\infty}c_{z,n}^{N}=\lim_{N\to\infty}\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{L^{2}(\Omega)},$

also we denote $p_{k,n}:=\langle
p_{k}(\cdot),\varphi_{n}(\cdot)\rangle_{L^{2}(\Omega)}$ for $k=1,\cdots,K$ and
$n\in\mathbb{N}^{+}$. For the coefficients
$\\{c_{z,n},c_{z,n}^{N},p_{k,n}\\}$, we give the following conditions, which
will be used in the future proof.

######  Assumption 2.4.

  * (a)

For $k\in\\{1,\cdots,K\\}$ and a.e. $z\in\partial\Omega$, we can find $C>0$
which is independent of $N$ such that
$\sum_{n=1}^{\infty}|c^{N}_{z,n}p_{k,n}|<C<\infty$ for $N\in\mathbb{N}^{+}$.

  * (b)

$\sum_{k=1}^{K}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|<\infty$ for a.e.
$z\in\partial\Omega$.

####  2.2.3 Auxiliary lemmas from [23].

At the end of this section, we collect some auxiliary results from the
reference [23].

######  Lemma 2.5.

Assuming that $\\{p_{k}(x)\\}_{k=1}^{K}$ possess appropriate regularities,
then for a.e. $z\in\partial\Omega$ and a.e. $t\in(0,\infty)$, it holds that

$-\int_{0}^{t}\frac{\partial u}{\partial\nu}(z,\tau)\
d\tau=\int_{0}^{t}\sum_{k=1}^{K}\chi_{{}_{t-\tau\in[t_{k-1},t_{k})}}\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})\tau})\Big{]}\
d\tau.$

###### Proof.

This is [23, Corollary 3.2]. ∎

The next two lemmas are included in [23, Section 3.3].

######  Lemma 2.6.

We denote the set of distinct eigenvalues with increasing order by
$\\{\lambda_{j}\\}_{j=1}^{\infty}$. For any nonempty open subset
$\Gamma\subset\partial\Omega$, if

$\sum_{\lambda_{n}=\lambda_{j}}c_{z,n}\eta_{n}=0\ \text{for}\
j\in\mathbb{N}^{+}\ \text{and\ a.e.}\ z\in\Gamma,$

then $\\{\eta_{n}\\}_{n=1}^{\infty}=\\{0\\}$.

######  Lemma 2.7.

Let $\Gamma$ be a nonempty open subset of $\partial\Omega$ and
$\sum_{n=1}^{\infty}c_{z,n}\eta_{n}$ be absolutely convergent for a.e.
$z\in\Gamma$. Given $\epsilon>0$, then

$\lim_{\operatorname{Re}s\to\infty}e^{\epsilon
s}\sum_{n=1}^{\infty}c_{z,n}\eta_{n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}=0\
\text{for a.e.}\ z\in\Gamma$

leads to $\\{\eta_{n}\\}_{n=1}^{\infty}=\\{0\\}$.

##  3 The proof of the uniqueness theorem.

###  3.1 The Laplace transform analysis.

The convolution structure in the result of Lemma (2.5) encourages us to apply
the Laplace transform, which is defined as

$\mathcal{L}\\{\psi(t)\\}(s)=\int_{0}^{\infty}e^{-st}\psi(t)\ dt,\quad
s\in\mathbb{C}.$

Not hard to see that for $\operatorname{Re}s>0$,

$\mathcal{L}\\{1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})t}\\}(s)=(\mu_{D}\lambda_{n}+\mu_{a})s^{-1}(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}.$

From Assumption (2.4), it holds that

$\displaystyle\Big{|}\int_{0}^{t}\sum_{k=1}^{K}\chi_{{}_{t-\tau\in[t_{k-1},t_{k})}}\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})\tau})\Big{]}\
d\tau\Big{|}$ $\displaystyle\leq
2\int_{0}^{t}\sum_{k=1}^{K}\big{|}\chi_{{}_{t-\tau\in[t_{k-1},t_{k})}}\big{|}\
\sum_{n=1}^{\infty}\big{|}c_{z,n}p_{k,n}\big{|}\ d\tau\leq Ct.$

Also we can see $|e^{-st}t|$ is integrable on $(0,\infty)$ if
$\operatorname{Re}s>0$. Then by the Dominated Convergence Theorem, taking
Laplace transform on the result in Lemma (2.5) yields that for each
$s\in\mathbb{C}^{+}:=\\{s\in\mathbb{C}:\operatorname{Re}s>0\\},$

$\displaystyle\mathcal{L}\left\\{-\int_{0}^{t}\frac{\partial
u}{\partial\nu}(z,\tau)\ d\tau\right\\}(s)$
$\displaystyle=\sum_{k=1}^{K}\sum_{n=1}^{\infty}\int_{0}^{\infty}e^{-st}\int_{0}^{t}\chi_{{}_{t-\tau\in[t_{k-1},t_{k})}}c_{z,n}p_{k,n}(1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})\tau})\
d\tau\ dt$
$\displaystyle=s^{-2}\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}\Big{]},$

which with $\mathcal{L}(\int\psi)=s^{-1}\mathcal{L}(\psi)$ implies that for
$s\in\mathbb{C}^{+}$,

$s\mathcal{L}\left\\{-\frac{\partial
u}{\partial\nu}(z,t)\right\\}(s)=\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}\Big{]}.$
(3.1)

After deducing (3.1), we need to show the well-definedness and analyticity of
the complex series in it. In the next lemma, we recall that the distinct
eigenvalues with increasing order are denoted by
$\\{\lambda_{j}\\}_{j=1}^{\infty}$.

######  Lemma 3.1.

Under Assumption (2.4), the series

$\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}\Big{]}$

is analytic on $\mathbb{C}^{+}$.

###### Proof.

Firstly, let us show the analyticity of the series
$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}$
on $\mathbb{C}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$ for
$k=1,\cdots,K$. Obviously we see that
$c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}$
is holomorphic on
$\mathbb{C}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$. So it
is sufficient to show the uniform convergence of the above series.

For $R>0$, we define $\mathbb{C}_{R}:=\\{s\in\mathbb{C}:|s|<R\\}$. Recalling
that $\lambda_{n}\to\infty$, then there exists a large $N_{1}>0$ such that
$\mu_{D}\lambda_{n}+\mu_{a}>2R$ for $n\geq N_{1}$. Sequentially, for
$s\in\mathbb{C}_{R}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$
and $n\geq N_{1}$, we have

$|s+\mu_{D}\lambda_{n}+\mu_{a}|\geq|\operatorname{Re}s+\mu_{D}\lambda_{n}+\mu_{a}|=\operatorname{Re}s+\mu_{D}\lambda_{n}+\mu_{a}\geq\mu_{D}\lambda_{n}+\mu_{a}-R,$

which gives

$|(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}|\leq(\mu_{D}\lambda_{n}+\mu_{a})(\mu_{D}\lambda_{n}+\mu_{a}-R)^{-1}<2.$

Given $\epsilon>0,$ Assumption (2.4) yields that there exists $N_{2}>0$ such
that for $l\geq N_{2}$, $\sum_{n=l}^{\infty}|c_{z,n}p_{k,n}|<\epsilon.$ So,
for $l\geq\max\\{N_{1},N_{2}\\}$ and
$s\in\mathbb{C}_{R}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$,

$\Big{|}\sum_{n=l}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}\Big{|}\leq
2\sum_{n=l}^{\infty}|c_{z,n}p_{k,n}|<2\epsilon,$

which implies the uniform convergence. With this uniform convergence result,
we see that the series
$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}$
is holomorphic on
$\mathbb{C}_{R}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$ for
each $R>0$. Given
$s_{0}\in\mathbb{C}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$,
we can find $R>0$ such that
$s_{0}\in\mathbb{C}_{R}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$,
which means
$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}$
is analytic on
$\mathbb{C}\setminus\\{-\mu_{D}\lambda_{j}-\mu_{a}\\}_{j=1}^{\infty}$.

Now, let us show the analyticity of

$\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}\Big{]}$

on $\mathbb{C}^{+}$. For the case of $K$ is finite, from the result that
$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}$
is holomorphic on $\mathbb{C}^{+}$, we can deduce the desired result
straightforwardly. If $K$ is infinity, for $s\in C^{+}$ we see that
$|(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}|\leq 1$ and
$|e^{-t_{k}s}|\leq 1$, so that

$\sum_{k=1}^{K}|e^{-t_{k-1}s}-e^{-t_{k}s}|\Big{[}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|\
|(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}|\Big{]}\leq
2\sum_{k=1}^{K}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|.$

Then the above proof and Assumption (2.4) give the uniform convergence of the
above series on $\mathbb{C}^{+}$. Sequentially, we can deduce the desired
result and the proof is complete. ∎

###  3.2 Proof of Theorem (1).

Here we will prove the main theorem. To shorten our proof, we define

$\displaystyle P_{z,k}(s)$
$\displaystyle:=\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1},$
$\displaystyle\tilde{P}_{z,k}(s)$
$\displaystyle:=\sum_{n=1}^{\infty}c_{z,n}\tilde{p}_{k,n}(\tilde{\mu}_{D}\lambda_{n}+\tilde{\mu}_{a})(s+\tilde{\mu}_{D}\lambda_{n}+\tilde{\mu}_{a})^{-1}.$

###### Proof of Theorem (1) (i) .

Given $\mu_{a}=\tilde{\mu}_{a}$ and $\mu_{D}=\tilde{\mu}_{D}$, from (3.1) and
Lemma (3.1), we have that for $s\in\mathbb{C}^{+}$ and $z\in\Gamma_{d}$,

$\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})P_{z,k}(s)=\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})\tilde{P}_{z,k}(s).$
(3.2)

We first prove that $p_{1,n}=\tilde{p}_{1,n}$ for $n\in\mathbb{N}^{+}$.
Multiplying $e^{(t_{0}+\epsilon)s}$ with sufficiently small $\epsilon>0$ such
that $\epsilon<t_{1}-t_{0}$ on (3.2) gives

$\displaystyle e^{\epsilon s}[P_{z,1}(s)-\tilde{P}_{z,1}(s)]=$ $\displaystyle
e^{(\epsilon+t_{0}-t_{1})s}P_{z,1}(s)-\sum_{k=2}^{K}(e^{(\epsilon+t_{0}-t_{k-1})s}-e^{(\epsilon+t_{0}-t_{k})s})P_{z,k}(s)$
(3.3)
$\displaystyle-e^{(\epsilon+t_{0}-t_{1})s}\tilde{P}_{z,1}(s)+\sum_{k=2}^{K}(e^{(\epsilon+t_{0}-t_{k-1})s}-e^{(\epsilon+t_{0}-t_{k})s})\tilde{P}_{z,k}(s).$

For $s\in\mathbb{C}^{+}$, not hard to see that

$|(\mu_{D}\lambda_{n}+\mu_{a})(s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}|\leq(\mu_{D}\lambda_{n}+\mu_{a})(\operatorname{Re}s+\mu_{D}\lambda_{n}+\mu_{a})^{-1}\leq
1.$

This with Assumption (2.4) yields that

$\displaystyle\lim_{\operatorname{Re}s\to\infty}\Big{|}\sum_{k=2}^{K}(e^{(\epsilon+t_{0}-t_{k-1})s}$
$\displaystyle-e^{(\epsilon+t_{0}-t_{k})s})P_{z,k}(s)\Big{|}$
$\displaystyle\leq\lim_{\operatorname{Re}s\to\infty}2e^{(\epsilon+t_{0}-t_{1})\operatorname{Re}s}\sum_{k=2}^{K}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|=0.$

Analogously, we can show that other terms in the right side of (3.3) tend to
zero as $\operatorname{Re}s\to\infty$. Now we have

$\lim_{\operatorname{Re}s\to\infty}e^{\epsilon
s}[P_{z,1}(s)-\tilde{P}_{z,1}(s)]=0\ \text{for}\ z\in\Gamma_{d}.$

Then with Lemma (2.7) we have $p_{1,n}=\tilde{p}_{1,n}\ \text{for}\
n\in\mathbb{N}^{+},$ namely $\|p_{1}-\tilde{p}_{1}\|_{L^{2}(\Omega)}=0$.

Next we need to show $p_{k,n}=\tilde{p}_{k,n}$ for $n\in\mathbb{N}^{+}$ and
$k>1$. From the result $p_{1,n}=\tilde{p}_{1,n},\ n\in\mathbb{N}^{+},$ we have
$P_{z,1}(s)=\tilde{P}_{z,1}(s)$. Inserting it into (3.2) yields that for
$s\in\mathbb{C}^{+}$ and $z\in\Gamma_{d}$,

$\sum_{k=2}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})P_{z,k}(s)=\sum_{k=2}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})\tilde{P}_{z,k}(s).$

Following the above proof gives that
$\|p_{2}-\tilde{p}_{2}\|_{L^{2}(\Omega)}=0$. Continuing this argument, we
conclude that $\|p_{k}-\tilde{p}_{k}\|_{L^{2}(\Omega)}=0$ for $k=1,\cdots,K$.
The proof is complete. ∎

Before to show Theorem (1) (ii), we state [33, Lemma 3.5] below.

######  Lemma 3.2.

Let $\\{\tau_{n}\\}_{n=1}^{\infty}$ be an absolutely convergent complex
sequence and $\\{\gamma_{n}\\}_{n=1}^{\infty}$ be a real sequence satisfying
$0\leq\gamma_{1}<\gamma_{2}<\cdots$ and $\gamma_{n}\to\infty.$ For the complex
series $\sum_{n=1}^{\infty}\tau_{n}e^{-\gamma_{n}t}$ defined on
$\mathbb{C}^{+}$, if the set of its zeros on $\mathbb{C}^{+}$ has an
accumulation point, then $\\{\tau_{n}\\}_{n=1}^{\infty}=\\{0\\}$.

###### Proof of Theorem (1) (ii).

Given the condition (1.7) and $\mu_{D}=\tilde{\mu}_{D}$, firstly let us show
$\mu_{a}=\tilde{\mu}_{a}$. Assume not, without loss of generality, we can set
$\mu_{a}<\tilde{\mu}_{a}$. From Lemma (2.5), it holds that for
$t\in(t_{0},t_{1}],$

$\displaystyle-\frac{\partial u}{\partial\nu}(z,t)$
$\displaystyle=\sum_{n=1}^{\infty}c_{z,n}p_{1,n}(1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})(t-t_{0})}),$
$\displaystyle-\frac{\partial\tilde{u}}{\partial\nu}(z,t)$
$\displaystyle=\sum_{n=1}^{\infty}c_{z,n}\tilde{p}_{1,n}(1-e^{-(\mu_{D}\lambda_{n}+\tilde{\mu}_{a})(t-t_{0})}).$

Then we have

$\sum_{n=1}^{\infty}c_{z,n}p_{1,n}(1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})(t-t_{0})})=\sum_{n=1}^{\infty}c_{z,n}\tilde{p}_{1,n}(1-e^{-(\mu_{D}\lambda_{n}+\tilde{\mu}_{a})(t-t_{0})})$
(3.4)

for $t\in(t_{0},t_{1}].$ The condition $\mu_{a}<\tilde{\mu}_{a}$ gives that
$\mu_{D}\lambda_{1}+\mu_{a}<\mu_{D}\lambda_{j}+\mu_{a}$ for $j>1$ and
$\mu_{D}\lambda_{1}+\mu_{a}<\mu_{D}\lambda_{j}+\tilde{\mu}_{a}$ for $j\geq 1$.
Using Lemma (3.2), we have that
$\sum_{\lambda_{n}=\lambda_{1}}c_{z,n}p_{1,n}=0$ for $z\in\Gamma_{d}$. This
together with Lemma (2.6) leads to $p_{1,n}=0$ for $\lambda_{n}=\lambda_{1}$,
which contradicts with condition (1.7). So we have $\mu_{a}=\tilde{\mu}_{a}$.

With $\mu_{a}=\tilde{\mu}_{a}$ and the proof for Theorem (1) (i), we deduce
that $\|p_{k}-\tilde{p}_{k}\|_{L^{2}(\Omega)}=0$ for $k=1,\cdots,K$. The proof
is complete. ∎

###### Proof of Theorem (1) (iii).

With condition (1.7) and $\mu_{a}=\tilde{\mu}_{a}$, let us show
$\mu_{D}=\tilde{\mu}_{D}$. Assuming that $\mu_{D}<\tilde{\mu}_{D}$, from (3.4)
we have

$\sum_{n=1}^{\infty}c_{z,n}p_{1,n}(1-e^{-(\mu_{D}\lambda_{n}+\mu_{a})(t-t_{0})})=\sum_{n=1}^{\infty}c_{z,n}\tilde{p}_{1,n}(1-e^{-(\tilde{\mu}_{D}\lambda_{n}+\mu_{a})(t-t_{0})}).$

With Lemmas (2.6) and (3.2), following the proof of Theorem (1) (ii) we can
get $p_{1,n}=0$ if $\lambda_{n}=\lambda_{1}$, which is a contradiction. Hence,
$\mu_{D}=\tilde{\mu}_{D}$.

With the result $\mu_{D}=\tilde{\mu}_{D}$ and the proof of Theorem (1) (i), we
have $\|p_{k}-\tilde{p}_{k}\|_{L^{2}(\Omega)}=0$ for $k=1,\cdots,K$. The proof
is complete. ∎

##  4 Numerical inversions.

Noting that $\mu_{D}:=1/(3\mu_{s}(1-g))$, where $g$ is the anisotropy
parameter and known, we may consider the numerical inversion of
$\textbf{a}:=(\mu_{a},\mu_{s},\mu_{f})$. In practice, for the FDOT people
inject the light from a laser source to the biological tissue via one of the
source fibers at boundary, then measure the amount of transmitted light at all
the boundary detector locations using the detector fibers. This process is
repeated for all the source locations. Now we show the mathematical
description of our FDOT as follows.

Let $\Gamma_{s}$ and $\Gamma_{d}$ be the finite set of source locations and
detector locations, respectively. Let $\mathcal{T}\subset(0,\infty)$ be a time
interval in which we take the measurements. Then, we can denote the exact
boundary measurements corresponding to any given input
$\textbf{a}\in\mathcal{A}$ with
$\big{(}x_{s}^{(n)},x_{d}^{(k)}\big{)}\in\Gamma_{s}\times\Gamma_{d}$ and
$t_{i}\in\mathcal{T}$ by
$u_{m}[\textbf{a}]\big{(}x_{d}^{(k)},t_{i};x_{s}^{(n)}\big{)}$. Here
$\mathcal{A}$ is the admissible set. Taking the measured data at $N$ different
excitation sources, $K$ different detectors and $I$ different times, the
computational inverse problem is to determine the unknown a from the following
set of measurements

$\left\\{u_{m}\big{(}x_{d}^{(k)},t_{i};x_{s}^{(n)}\big{)}:n=1,\cdots,N,\
k=1,\cdots,K,\ i=1,\cdots,I\right\\},$ (4.1)

where the set $\\{x_{d}^{(k)}\\}_{k=1}^{K}$ and the set
$\\{t_{i}\\}_{i=1}^{I}$ may be varied for different source locations. More
precisely, for each given source $x_{s}^{(n)}$, the process of obtaining the
finite set of measurements ((4.1)) can be described as

$\mathbb{K}_{n}:\mathcal{A}\to\mathbb{R}^{K\times I},\
\mathbb{K}_{n}(\textbf{a})=u_{m}[\textbf{a}]\big{(}x_{d}^{(k)},t_{i};x_{s}^{(n)}\big{)},\
k=1,\cdots,K,\ i=1,\cdots,I.$ (4.2)

For convenience, we set $M=N\times K$ and denote
$\\{\omega_{m}:=\big{(}x_{s}^{(m)},x_{d}^{(m)}\big{)}\\}_{m=1}^{M}$ be the $M$
different S-D pairs we used in (4.2), where the source location or detector
location for different S-D pairs may be same. Then we can rewrite ((4.2)) as

$\mathbb{K}:\mathcal{A}\to\mathbb{R}^{M\times
I},\quad\mathbb{K}(\textbf{a})=h,$ (4.3)

where $h:=u_{m}(\hat{\omega},\hat{t})$ is the measured data corresponding to
S-D pairs $\hat{\omega}:=(\omega_{1},\cdots,\omega_{M})$ and time points
$\hat{t}\in\mathcal{T}^{I}$. Then our FDOT is to solve the equation (4.3).

In the follows, unless in particular cases, we always take the physical
parameters as

$\displaystyle
c=0.219\;\rm{mm}/{ps},\quad\mu_{s}^{\prime}=1.0\;\rm{mm}^{-1},\quad\mu_{a}=0.1\;\rm{mm}^{-1},\quad\beta=0.01\;\rm{mm}^{-1}.$

We use mm as the unit of length. For simplicity, we consider the zero-lifetime
case of (1.3), i.e. the lifetime in the source term $S_{0}[\mu_{f},u_{e}]$ is
$\tau=0$. Further, we are here focusing on thick ($>1$ cm) or large volume
tissue ($>10$ cm3) like chest and thus we may set

$\Omega:=\mathbb{R}^{3}_{+}=\left\\{(x_{1},x_{2},x_{3}):(x_{1},x_{2})\in\mathbb{R}^{2},x_{3}>0\right\\}$

with the boundary
$\partial\Omega:=\left\\{(x_{1},x_{2},0):(x_{1},x_{2})\in\mathbb{R}^{2}\right\\}$.

###  4.1 Identify a stationary target by peak detection.

Supposing $\mu_{s}$ and $\mu_{a}$ are known, we consider to identify the
stationary $\mu_{f}(x)$. For the ideal case that the size of fluorescence
target is very small, we may concern the point target, i.e.

$\mu_{f}(x)=P\delta(x-x_{c}),$

where $x_{c}=(x_{c1},x_{c2},x_{c3})\in\Omega$ is the location of point target
and $P>0$ is the concentration of fluorescent target. Then, our FDOT is to
identify $\textbf{a}=(x_{c1},x_{c2},x_{c3},P)$. For each given
$x_{d}\in\partial\Omega$ and $x_{s}\in\partial\Omega$, we have that

$u_{m}(t)=C(t)P\int_{0}^{t}\frac{1}{\big{[}(t-s)s\big{]}^{3/2}}e^{-\frac{|x_{d}-x_{c}|^{2}}{4c\mu_{D}(t-s)}}e^{-\frac{|x_{s}-x_{c}|^{2}}{4c\mu_{D}s}}K_{3}(0,x_{c3};t-s)K_{3}(x_{c3},0;s)\,ds,$
(4.4)

where $C(t):=\frac{e^{-c\mu_{a}t}}{16\pi^{3}c{\mu_{D}}^{2}}$ and

$K_{3}(x_{3},y_{3};t):=1-\beta\sqrt{\pi
c\mu_{D}t}\,\exp{\left(\left(\frac{x_{3}+y_{3}+2\beta
c\mu_{D}t}{\sqrt{4c\mu_{D}t}}\right)^{2}\right)}\mathop{\mathrm{erfc}}\left(\frac{x_{3}+y_{3}+2\beta
c\mu_{D}t}{\sqrt{4c\mu_{D}t}}\right).$

The expression (4.4) is a temporal point-spread function (TPSF). Before
introducing the inversion strategy by peak detection, we define two geometric
concepts for TPFS as follows (see Figure (1)):

  * •

Peak intensity $u_{m}^{\rm peak}$: the TPSF maximum;

  * •

Peak time $t_{\rm peak}$: the temporal position of TPSF maximum.

The existence and uniqueness of the peak time for (4.4) are mathematically
investigated in [17]. There holds that

$u_{m}^{\rm peak}:=u_{m}^{\rm peak}[x_{s},x_{d}]\approx C(t_{\rm
peak})P\left(\sqrt{\frac{\pi}{A}}+\sqrt{\frac{\pi}{B}}\right)\exp\left(-\frac{2A+2B}{t_{\rm
peak}}\right){t_{\rm peak}}^{-\frac{3}{2}}$ (4.5)

and

$t_{\rm peak}:=t_{\rm
peak}[x_{s},x_{d}]\approx\frac{1}{4c\mu_{a}}\left(-3+\sqrt{9+32c\mu_{a}(A+B)}\right),$
(4.6)

where

$A:=\frac{|x_{d}-x_{c}|^{2}}{4c\mu_{D}},\quad
B:=\frac{|x_{s}-x_{c}|^{2}}{4c\mu_{D}}.$

Figure 1: The profile of TPSF for a point target located at $x_{c}=(0,0,5)$.

From above expressions we can see that the shape of the profile of $u_{m}$ is
actually dominated by the distances $|x_{d}-x_{c}|$ and $|x_{s}-x_{c}|$.
Fixing the distance between detector and source, we have stronger peak
intensity as the S-D pair closing to the target and we can obtain the
strongest one when the S-D pair satisfies $\frac{x_{s1}+x_{d1}}{2}=x_{c1}$ and
$\frac{x_{s2}+x_{d2}}{2}=x_{c2}$. This implies that we can identify
$(x_{c1},x_{c2})$ by comparing the peak intensity corresponding to different
S-D pairs. Since the peak time is not dependent on the concentration $P$, we
can further identify the depth of target ($x_{c3}$) from the expression of
peak time as in (4.6). Finally, the concentration $P$ can be identified by
(4.5). Hence the horizontal location, vertical location and concentration $P$
of the target can be separately determined by scanning the S-D pair at the
boundary surface. This inversion strategy by peak detection is summarized as
follows.

  * •

Step 1 (identify $(x_{c1},x_{c2})$): Fix the distance $|x_{s}-x_{d}|$ and
search an S-D pair $(x_{s}^{*},x_{d}^{*})$ such that

${u_{m,*}^{\rm peak}}:=u_{m}^{\rm
peak}[x_{s}^{*},x_{d}^{*}]=\max_{(x_{s},x_{d})\in\partial\Omega\times\partial\Omega}\big{\\{}u_{m}^{\rm
peak}[x_{s},x_{d}]\big{\\}}.$

Then the horizontal location of the target is identified by

$(x_{c1},x_{c2})=\left(\frac{x_{s1}^{*}+x_{d1}^{*}}{2},\frac{x_{s2}^{*}+x_{d2}^{*}}{2}\right).$

  * •

Step 2 (identify $x_{c3}$): By (4.6), the depth of the target is identified by

$x_{c3}=\sqrt{c^{2}\mu_{D}\mu_{a}\big{(}t_{\rm
peak}^{*}\big{)}^{2}+\frac{3}{2}c\mu_{D}t_{\rm
peak}^{*}-\frac{1}{4}|x_{s}^{*}-x_{d}^{*}|^{2}},$

where $t_{\rm peak}^{*}$ is the measured peak time for the S-D pair
$(x_{s}^{*},x_{d}^{*})$ obtained from Step 1.

  * •

Step 3 (identify $P$): By (4.5), the concentration of the target is identified
by

$P=\left.{u_{m,*}^{\rm peak}}\middle/\left\\{C(t_{\rm
peak}^{*})\left(\sqrt{\frac{\pi}{A^{*}}}+\sqrt{\frac{\pi}{B^{*}}}\right)\exp\left(-\frac{2A^{*}+2B^{*}}{t_{\rm
peak}^{*}}\right){t_{\rm peak}^{*}}^{-\frac{3}{2}}\right\\}\right..$

Now we test the accuracy of (4.5) and (4.6) to ensure the effectiveness of
above inversion strategy. We start with the point target given in the
following example.

Example 1. Set $x_{s}=(-3,0,0),x_{d}=(3,0,0)$ and $P=10^{6}\;{\rm mm}^{-1}$.
Suppose the point target is located at

$x_{c}=(0,0,x_{c3}).$

We take $x_{c3}\in[0.05,6]$ and show the comparisons of (4.5) in Figure (2)
(a) and (4.6) in Figure (2) (b), respectively.

Figure 2: Comparisons of (4.5) and (4.6) with the exact ones. (point target)

From Figure (2), it can be observed that (4.5) and (4.6) provide a good
approximation to the exact peak intensity and peak time, respectively. The
approximate errors decrease as the depth of point target becomes deeper.

By (4.5) and (4.6), we can further give the expressions of approximate peak
intensity and peak time for a small target. Suppose $\Omega_{0}$ is the
distribution domain of this small target and $x_{c}$ is its center point. By
the middle value theorem we have
$u_{m}[\Omega_{0}](t)\approx|\Omega_{0}|u_{m}[x_{c}](t)$, where $|\Omega_{0}|$
is the volume of small target. Then the peak time of small target will be
given by (4.6), while its peak intensity can be approximated by $u_{m}^{\rm
peak}[\Omega_{0}]\approx|\Omega_{0}|u_{m}^{\rm peak}[x_{c}]$, where
$u_{m}^{\rm peak}[x_{c}]$ is the peak intensity of the point target located at
$x_{c}$. Now we have tested the accuracy of (4.5) and (4.6) to approximate the
peak intensity and peak time of a small target. The case of a cubic target
will be given in the following example.

Example 2. Set $x_{s}=(-3,0,0)$, $x_{d}=(3,0,0)$ and $P=10^{6}\;{\rm
mm}^{-1}$. Suppose the cubic target with side length $L$ is located at
$x_{c}=(0,0,5)$. We take $L\in[0.01,1]$ and show the comparisons of the above
approximate expressions with the exact ones in Figure (3) (a) and (b),
respectively.

Figure 3: Comparisons of the approximate peak intensity and peak time with the
exact ones. (cubic target)

Figure (3) shows that the peak intensity and peak time of a small target can
be approximated well by the ones of a point target located at its center. Even
if $L=1$ mm, the relative error for peak time is $1.59e-2$, implying that
(4.6) is still useful for the case of small target.

Summarizing above numerical results, we realize that the explicit expressions
of peak intensity and peak time are accurate, and hence our proposed inversion
strategy is effective to identify the location of a small target. However, we
point out that the validity of this strategy depends on the relative distance
$|A-B|$, which requires $|A-B|$ is small enough [17].

###  4.2 Identify a target by iterative algorithm.

The FDOT is to solve the operator equation (4.3). For exact observation data,
if one of $\mu_{s}$ and $\mu_{a}$ is known, this equation has a unique
solution from Theorem (1). For given noisy data $h^{\delta}$ satisfying

$\|h^{\delta}-h\|_{2}\leq\delta,$

we consider the approximate solution of (4.3) by the minimizer of the cost
functional

$\mathcal{J}_{\delta,\alpha}(\mathbf{a}):=\|\mathbb{K}(\mathbf{a})-h^{\delta}\|^{2}_{2}+\alpha\|\mathbf{a}\|^{2}_{2},$
(4.7)

where $\alpha>0$ is the regularization parameter. Supposing the fluorescence
target is often of regular shape with a constant absorbing coefficient, we can
describe the distribution of $\mu_{f}$ in $\mathbb{R}^{3}_{+}$ as well as its
interface by a finite dimensional vector $\textbf{a}\in\mathbb{R}^{S}$. Here
$S$ is the dimension number of vector a. Then we can transform our FDOT into a
finite dimensional inverse problem. For instance, if we suppose the stationary
target is a cuboid, we have

$\displaystyle\mu_{f}(x)=\begin{cases}P,&x\in\Omega_{0}:=\left\\{x:\,x_{1}\in(a_{1},b_{1}),\,x_{2}\in(a_{2},b_{2}),\,x_{3}\in(a_{3},b_{3})\right\\},\\\
0,&x\not\in\Omega_{0},\end{cases}$

where $P>0$ is the concentration of fluorescent target and $b_{3}>a_{3}>0$.

Now we give a sketch for solving this minimization problem. In fact, solving
(4.7) can be transformed to solve a normal equation combining with an
iteration process:

$\begin{cases}\begin{aligned} (\alpha
I+\mathbb{G}^{T}\mathbb{G})\delta{\mathbf{a}_{j}}&=\mathbb{G}^{T}(h-\mathbb{K}(\mathbf{a}_{j})),\\\
{\mathbf{a}_{j}}+\delta{\mathbf{a}_{j}}&\to{\mathbf{a}_{j+1}},\quad
j=0,1,\cdots,\end{aligned}\end{cases}$

where $\delta{\mathbf{a}_{j}}$ is a perturbation vector for any given
${\mathbf{a}_{j}}\in\mathcal{A}$, $j$ denotes the iterative number and
${\mathbf{a}_{0}}$ is the initial guess. Here $\alpha>0$ is selected by the
discrepancy principle. The sensitivity matrix $\mathbb{G}$ is given by

$\displaystyle\mathbb{G}={(g_{qs})}_{Q\times S},\quad Q=N\times K\times I,$

where $g_{qs}$ is the derivative of $u_{m}$ on the $s$-th component of a,
corresponding to the $q$-th S-D pair of total $N\times K$ S-D pairs. We
terminate the iteration process by $\|\delta\mathbf{a}_{j}\|_{2}\leq\eta$ for
some specified $\eta>0$, which is taken as $\eta=10^{-8}$ for our numerics.

Figure 4: Horizontal projection of the configuration: the middle red cubic is
the projection of a target to the boundary $\partial\Omega$, the red points
are the excitation points, and the blue circles are the detector points.

Example 3. Set $\mu_{s}=10\;\rm{mm}^{-1}$ and $\mu_{a}=2\;\rm{mm}^{-1}$. We
suppose the distribution of cuboid target is

$\displaystyle\Omega_{0}:=\left\\{(x_{1},x_{2},x_{3}):x_{1}\in(-1,1),\
x_{2}\in(-1,1),\ x_{3}\in(9,11)\right\\}$

with the concentration $P=0.5\;\rm{mm}^{-1}$, i.e.

$\mu_{f}(x)=\begin{cases}P,&x\in\Omega_{0},\\\
0,&x\not\in\Omega_{0}.\end{cases}$

For this example, we observe the measurements from the setup shown in Figure
(4) and test our inversion algorithm for the following two cases:

  * (3a)

given $(\mu_{s},\mu_{a})$, identify stationary $\mu_{f}(x)$;

  * (3b)

simultaneously identify $(\mu_{s},\mu_{a},\mu_{f}(x))$.

As in Figure (4), we have 4 excitation sources and observe the measurements at
8 detector points for each source point. Hence we have 32 S-D pairs in total
for this setup. The sources and detectors are located at

$S_{1}=(-13,-13),\;S_{2}=(9,-13),\;S_{3}=(9,9),\;S_{4}=(-13,9)$

and

$D_{1}=(-13,-6),\;D_{2}=(-6,-18),\;D_{3}=(2,-18),\;D_{4}=(9,-6),$
$D_{5}=(9,2),\;D_{6}=(2,14),\;D_{7}=(-6,14),\;D_{8}=(-13,2).$

For each S-D pair, we select the peak time $t_{\rm peak}$ and choose 20 time
points $\hat{t}:=[t_{\rm peak}-10\Delta t:\Delta t:t_{\rm peak}+9\Delta t]$
with time step $\Delta t=2$ ps, so that the measurement $h$ is a
$640$-dimensional vector. The noisy data $h^{\delta}$ of $h$ is described by

$h^{\delta}=h(1+\zeta\epsilon),$ (4.8)

where $\epsilon>0$ is the noise level and $\zeta$ is a random standard
Gaussian noise.

For the case $(3a)$, we set the exact solution as

$\displaystyle\textbf{a}_{\rm exa}=(-1,1,-1,1,9,11,0.5).$

Noting the randomness of the noisy data (4.8), we carry out 10-time recoveries
from different noisy data set for each given noise level and take an average
of the recoveries. The recovered results from the initial guess
$\mathbf{a}_{0}=(-5.1,-4.9,-2.1,-1.9,5.9,6.1,0.1)$ are listed in Table (1),
where $\textbf{a}_{\rm rec}^{\rm avr}$ denotes the average of 10-time
recoveries and $Err$ is the $L^{2}$ relative error in recoveries given by

$Err:=\frac{\left\|\textbf{a}_{\rm exa}-\textbf{a}_{\rm rec}^{\rm
avr}\right\|_{2}}{\left\|\textbf{a}_{\rm exa}\right\|_{2}}.$

The results of the 10-time recoveries are plotted in Figure (5).

Table 1: The average of 10-time recoveries for the case $(3a)$. $\epsilon$ | ${{\textbf{a}}_{\rm rec}^{\rm avr}}=(a_{1}^{*},b_{1}^{*},a_{2}^{*},b_{2}^{*},a_{3}^{*},b_{3}^{*},P^{*})$ | $Err$
---|---|---
$5\%$ | $(-0.85,0.92,-0.85,0.92,9.01,10.78,0.81)$ | $3.93e-2$
$1\%$ | $(-0.98,0.99,-0.98,0.99,9.00,10.97,0.52)$ | $3.50e-3$

Figure 5: The 10-time recoveries with different $\epsilon$ for the case
$(3a)$.

For the case $(3b)$, we need to consider the simultaneous inversion of
$(\mu_{s},\mu_{a},\mu_{f}(x))$. This implies that we have 9 unknown parameters
and the exact solution is

$\displaystyle\textbf{a}_{\rm exa}=(10,2,-1,1,-1,1,9,11,0.5).$

Setting the initial guess
$\mathbf{a}_{0}=(5,0.5,-5.1,-4.9,-2.1,-1.9,5.9,6.1,0.1)$, the recovered
results are shown in Table (2) and Figure (6).

Table 2: The average of 10-time recoveries for the case $(3b)$. $\epsilon$ | ${{\textbf{a}}_{\rm rec}^{\rm avr}}=(\mu_{s}^{*},\mu_{a}^{*},a_{1}^{*},b_{1}^{*},a_{2}^{*},b_{2}^{*},a_{3}^{*},b_{3}^{*},P^{*})$ | $Err$
---|---|---
$5\%$ | $(10.04,2.02,-1.01,0.98,-0.99,0.99,9.04,11.03,0.80)$ | $2.12e-2$
$1\%$ | $(10.01,2.00,-1.00,1.00,-1.00,1.00,9.01,11.01,0.53)$ | $2.20e-3$

Figure 6: The 10-time recoveries with different $\epsilon$ for the case
$(3b)$.

From Figures (5) and (6), we see that the recoveries are satisfactory in view
of the high ill-posedness of our inverse problem, which is caused by the
nonlinearity of this inverse problem and the sparsity of the boundary
measurements. Further, the results for case $(3b)$ indicate that we may
reconstruct $(\mu_{s},\mu_{a},\mu_{f}(x))$ simultaneously, which is not
covered by Theorem (1).

In what follows, we further test the recovery performances of our inversion
strategy for a time-dependent fluorescence target.

Example 4. Set $\mu_{s}=10\;\rm{mm}^{-1},\,\mu_{a}=2\;\rm{mm}^{-1}$ and
$P=0.5\;\rm{mm}^{-1}$. We recover a time-dependent cubic target for the
following two cases:

  * $(4a)$

fixed center $x_{c}=(0,0,10)$, time-dependent side length
$L(t)=\frac{1}{2}+t$;

  * $(4b)$

fixed side length $L=2$ mm, time-dependent center $x_{c}(t)$ with

$\begin{cases}x_{c1}=-4+6t,\\\ x_{c2}=-5+8t,\\\ x_{c3}=8.\end{cases}$

In this example, we suppose either the side length or the center of the target
depends on the time linearly and we are going to recover the time-dependent
target by identifying its expansion coefficients in polynomial space. We also
obtain the measurements from the setup in Figure (4). For the case $(4a)$,
setting the initial guess $(-2,-2,5,0.2)$ for $(x_{c},P)$ and the initial
guess $\frac{1}{10}+\frac{1}{5}t$ for $L(t)$, the recovered results are
plotted in Figures (7) and (8). For the case $(4b)$, setting the initial guess
$(-1+2t,-2+3t,4)$ for $x_{c}(t)$ and the initial guess $(0.2,0.1)$ for
$(L,P)$, the numerical performances of our inversion scheme are shown in
Figures (9) and (10).

Figures (7)–(10) show that the reconstructions are satisfactory. The recovered
targets are in good agreement with the exact shape and its information of the
time evolution process can be reconstructed well. Hence, summarizing the
numerical results obtained from Examples 3–4, our proposed iterative algorithm
is effective to recover the fluorescence target. Further, although we do the
inversion only for a cubic or cuboid target in above examples, the position of
the target with a general shape can be identified by cuboid approximation [8].

Figure 7: The exact and recovered results for the case $(4a)$
($\epsilon=1\%$).
Figure 8: The exact solution (blue) and the recovered solution (red) for the
case $(4a)$ ($\epsilon=1\%$).
Figure 9: The exact and recovered $x_{c}(t)$ for the case $(4b)$
($\epsilon=1\%$).
Figure 10: The exact solution (blue) and recovered solution (red) for the case
$(4b)$ ($\epsilon=1\%$).

##  Acknowledgement

Chunlong Sun is supported by National Natural Science Foundation of China
(Grant No. 11971104) and by Natural Science Foundation of Jiangsu Province,
China (Grant No. BK20210268). Zhidong Zhang is supported by National Natural
Science Foundation of China (Grant No. 12101627).

## References

  * [1] D. S. Anikonov.  Uniqueness of the simultaneous determination of two coefficients of the transport equation.  Dokl. Akad. Nauk SSSR, 277(4):777–779, 1984.
  * [2] D. S. Anikonov.  Uniqueness of the determination of the coefficient of the transport equation with a special type of source.  Dokl. Akad. Nauk SSSR, 284(5):1033–1037, 1985.
  * [3] S. R. Arridge.  Optical tomography in medical imaging.  Inverse Problems, 15(2):R41–R93, 1999.  URL: https://doi.org/10.1088/0266-5611/15/2/022, doi:10.1088/0266-5611/15/2/022.
  * [4] S. R. Arridge and Wrb Lionheart.  Nonuniqueness in diffusion-based optical tomography.  Optics Letters, 23(11):882–884, 1998.
  * [5] Simon R. Arridge and John C. Schotland.  Optical tomography: forward and inverse problems.  Inverse Problems, 25(12):123010, 59, 2009.  URL: https://doi.org/10.1088/0266-5611/25/12/123010, doi:10.1088/0266-5611/25/12/123010.
  * [6] Guillaume Bal.  Inverse transport theory and applications.  Inverse Problems, 25(5):053001, 48, 2009.  URL: https://doi.org/10.1088/0266-5611/25/5/053001, doi:10.1088/0266-5611/25/5/053001.
  * [7] Maïtine Bergounioux, Élie Bretin, and Yannick Privat.  How to position sensors in thermo-acoustic tomography.  Inverse Problems, 35(7):074003, 25, 2019.  URL: https://doi.org/10.1088/1361-6420/ab0e4d, doi:10.1088/1361-6420/ab0e4d.
  * [8] Sun C, Nakamura G, Nishimura G, Jiang Y, Liu J, and Machida M.  Fast and robust reconstruction algorithm for fluorescence diffuse optical tomography assuming a cuboid target.  J Opt Soc Am A Opt Image Sci Vis, 37(2):231–239, 2020.  doi:10.1364/JOSAA.37.000231.
  * [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):763–785, 1996.  URL: https://doi.org/10.1080/03605309608821207, doi:10.1080/03605309608821207.
  * [10] A. Corlu, R. Choe, T. Durduran, M. A. Rosen, and A. G. Yodh.  Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans.  Optics Express, 15(11):6696–716, 2007.
  * [11] Teresa Correia, Maximilian Koch, Angelique Ale, Vasilis Ntziachristos, and Simon Arridge.  Patch-based anisotropic diffusion scheme for fluorescence diffuse optical tomography - part 2: Image reconstruction.  Physics in Medicine and Biology, 61(4):1452–1475, 2016.
  * [12] E. Demidenko, A. Hartov, N. Soni, and K. D. Paulsen.  On optimal current patterns for electrical impedance tomography.  IEEE Trans Biomed Eng, 52(2):238–248, 2005.
  * [13] N. Ducros, C D’Andrea, A. Bassi, and F. Peyrin.  Fluorescence diffuse optical tomography: Time-resolved versus continuous-wave in the reflectance configuration.  IRBM, 32(4):243–250, 2011.
  * [14] T Durduran, R Choe, W B Baker, and A G Yodh.  Diffuse optics for tissue monitoring and tomography.  Reports on Progress in Physics, 73(7):076701, jun 2010.  URL: https://doi.org/10.1088/0034-4885/73/7/076701, doi:10.1088/0034-4885/73/7/076701.
  * [15] J. Dutta, S. Ahn, C. Li, S. R. Cherry, and R. M. Leahy.  Joint l1 and total variation regularization for fluorescence molecular tomography.  Physics in Medicine and Biology, 57(6):1459–1476, 2015.
  * [16] Joyita Dutta, Sangtae Ahn, Anand A Joshi, and Richard M Leahy.  Illumination pattern optimization for fluorescence tomography: theory and simulation studies.  Physics in Medicine and Biology, 55(10):2961–2982, apr 2010.  URL: https://doi.org/10.1088/0031-9155/55/10/011, doi:10.1088/0031-9155/55/10/011.
  * [17] J.Y Eom, M. Machida, G. Nakamura, G. Nishimura, and C.L. Sun.  Expression of the peak time for time-domain boundary measurements in diffuse light.  doi:arXiv:1907.00719v2.
  * [18] Feng Gao, Huijuan Zhao, Limin Zhang, Yukari Tanikawa, Andhi Marjono, and Yukio Yamada.  A self-normalized, full time-resolved method for fluorescence diffuse optical tomography.  Opt. Express, 16(17):13104–13121, Aug 2008.  doi:10.1364/OE.16.013104.
  * [19] Bastian Harrach.  On uniqueness in diffuse optical tomography.  Inverse Problems, 25(5):055010, 14, 2009.  URL: https://doi.org/10.1088/0266-5611/25/5/055010, doi:10.1088/0266-5611/25/5/055010.
  * [20] Nuutti Hyvönen, Aku Seppänen, and Stratos Staboulis.  Optimizing electrode positions in electrical impedance tomography.  SIAM J. Appl. Math., 74(6):1831–1851, 2014.  URL: https://doi.org/10.1137/140966174, doi:10.1137/140966174.
  * [21] H. B. Jiang.  Diffuse optical tomography: principles and applications.  CRC Press, Taylor $\&$ Francis Group, Boca Raton, 2011.
  * [22] S. Lam, F. Lesage, and X. Intes.  Time domain fluorescent diffuse optical tomography: analytical expressions.  Opt. Express, 13(7), 2263-2275, 2005.
  * [23] G. Lin, Z. Zhang, and Z. Zhang.  Theoretical and numerical studies of inverse source problem for the linear parabolic equation with sparse boundary measurements.  2021\.  doi:arXiv:2111.02285.
  * [24] J.J. Liu, M. Machida, G. Nakamura, G. Nishimura, and C.L. Sun.  On fluorescence imaging: The diffusion equation model and recovery of the absorption coefficient of fluorophores.  Sci. China Math, 2020.  doi:10.1007/s11425-020-1731-y.
  * [25] Yan Liu, Wuwei Ren, and Habib Ammari.  Robust reconstruction of fluorescence molecular tomography with an optimized illumination pattern.  Inverse Probl. Imaging, 14(3):535–568, 2020.  URL: https://doi.org/10.3934/ipi.2020025, doi:10.3934/ipi.2020025.
  * [26] F. Marttelli, S. D. Bianco, A. Ismaelli, and G. Zaccanti.  Light propagation through biological tissue and other diffusive media: theory, solutions, software.  SPIE Press, Bellingham, Washington, 2010.
  * [27] V. Nitziachristos and R. Weissleder.  Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized born approximation.  Optics letters, 26(12):893–5, 2001.
  * [28] V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder.  Fluorescence molecular tomography resolves protease activity in vivo.  Nature Medicine, 2002.
  * [29] S.V. Patwardhan, S.R. Bloch, S. Achilefu, and J.P. Culver.  Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice.  Opt. Express, 13(7), 2564-2577, 2005.
  * [30] B. W. Pogue and M. A. Mycek.  Handbook of Biomedical Fluorescence.  Marcel Dekker, New York, 2003.
  * [31] Vladimir G. Romanov and Sailing He.  Some uniqueness theorems for mammography-related time-domain inverse problems for the diffusion equation.  Inverse Problems, 16(2):447–459, 2000.  URL: https://doi.org/10.1088/0266-5611/16/2/312, doi:10.1088/0266-5611/16/2/312.
  * [32] M. Rudin.  Molecular imaging: basic principles and applications in Biomedical research, 2nd ed. Imperial College Press, London, 2013.
  * [33] William Rundell and Zhidong Zhang.  On the identification of source term in the heat equation from sparse data.  SIAM J. Math. Anal., 52(2):1526–1548, 2020.  URL: https://doi.org/10.1137/19M1279915, doi:10.1137/19M1279915.