# Simultaneous uniqueness and numerical inversion for an inverse problem in
the time-domain diffuse optical tomography with fluorescence

Zhiyuan Li  School of Mathematics and Statistics, Ningbo University, Ningbo
315211, Zhejiang, China  Chunlong Sun Corresponding author: Dr. Chunlong Sun,
email<EMAIL_ADDRESS>School of Mathematics, Nanjing University of
Aeronautics and Astronautics, Nanjing 211106, Jiangsu, China  Key Laboratory
of Mathematical Modelling and High Performance Computing of Air Vehicles
(NUAA), MIIT, Nanjing, 211106, China  Nanjing Center for Applied Mathematics,
Nanjing 211135, Jiangsu, China

###### Abstract

In this work, an inverse problem on the determination of multiple coefficients
arising from the time-domain diffuse optical tomography with fluorescence
(DOT-FDOT) is investigated. We simultaneously recover the distribution of
background absorption coefficient, photon diffusion coefficient as well as the
fluorescence absorption in biological tissue by the time-dependent boundary
measurements. We build the uniqueness theorem of this multiple coefficients
simultaneous inverse problem. After that, the numerical inversions are
considered. We introduce an accelerated Landweber iterative algorithm and give
several numerical examples illustrating the performance of the proposed
inversion schemes.

Keywords: diffusion equation, simultaneous inverse problem, single
measurement, uniqueness,

AMS Subject Classifications: 35R30, 65M32.

##  1 Introduction.

###  1.1 Background and literature.

Optical tomography is rapidly gaining acceptance as an important diagnostic
and monitoring tool of symptoms in medical applications [32, 33, 35]. Compared
to X-ray computed tomography (CT), magnetic resonance imaging (MRI) and
positron emission tomography (PET), optical tomography possesses some
advantages like low cost and portable, and therefore it has drawn more and
more attention. Light propagation in highly scattering medium, such as
biological tissue, is dominated by multiple scattering, and can be governed by
the Boltzmann radiative transfer equation (RTE) [14, 21, 29]. However, the
cost of solving RTE is extremely high, and therefore people approximately
expand the radiant intensity as an isotropic photon density plus a small
photon flux, and sequentially the transport equation can be reduced to a
diffusion equation. The diffusion equation can be a sufficiently accurate
approximation to the RTE, and has been applied for the quantitative analysis
of optical properties of random media [3, 5] due to its much lower computation
than RTE. This kind of analysis includes diffuse optical spectroscopy (DOS)
and diffuse optical tomography (DOT). Suppose the medium has the diffusion
coefficient $D$ and absorption coefficient $\mu_{a}$. Both of DOS and DOT are
to identify the unknown $(D,\mu_{a})$ from the information of light measured
on the boundary.

Fluorescence diffuse optical tomography (FDOT) is one type of diffuse optical
tomography that uses fluorescence light from fluorophores in biological
tissue, for which two processes are coupled, namely, excitation and emission
(fluorescence). The excitation photons injected on the boundary of the medium
propagate to the fluorophores and then some photons are absorbed by them which
excite the fluorophore molecules. After a moment of the absorption, the
fluorophores emit other photons, fluorescence, at more longer wavelength than
the wavelength of the excitation photons. Let $\Omega\subset\mathbb{R}^{d}$ be
the background medium with its boundary $\partial\Omega$. We denote $\mu_{f}$
the distributions of fluorophores inside $\Omega$. FDOT is a method to achieve
imaging of $\mu_{f}$ from the measured fluorescence on boundary. However, most
existing works focus on the linear FDOT problem, which is predicated on
knowing the background absorption $\mu_{a}$ and diffusion coefficient $D$, and
the nonlinear dependence of excitation on $\mu_{f}$ is ignored. In this work,
we study diffuse optical tomography with fluorescence (DOT-FDOT) in a more
general framework where $(D,\mu_{a})$ are unknown and the excitation also
depends on $\mu_{f}$. It is very nature to consider this simultaneous inverse
problem but really challenging.

Let $u_{e}(x,t)$ and $u_{m}(x,t)$ be the photon density of excitation light
and emission light, respectively. As a forward problem for DOT-FDOT, we
consider the following initial-boundary value problems for $u_{e}$ and $u_{m}$

$\begin{cases}\begin{aligned} c^{-1}\partial_{t}u_{e}-\text{div}(D(x)\nabla
u_{e})+(\mu_{a}(x)+\mu_{f}(x))u_{e}&=0,&\quad&(x,t)\in\Omega\times(0,T),\\\
u_{e}&=0,&\quad&(x,t)\in\Omega\times\\{0\\},\\\
u_{e}&=g,&\quad&(x,t)\in\partial\Omega\times(0,T),\end{aligned}\end{cases}$
(1.1)

and

$\begin{cases}\begin{aligned} c^{-1}\partial_{t}u_{m}-\text{div}(D(x)\nabla
u_{m})+\mu_{a}(x)u_{m}&=\mu_{f}(x)u_{e},&&(x,t)\in\Omega\times(0,T),\\\
u_{m}&=0,&&(x,t)\in\Omega\times\\{0\\},\\\
u_{m}&=0,&&(x,t)\in\partial\Omega\times(0,T).\end{aligned}\end{cases}$ (1.2)

Here $c>0$ is the speed of light in the medium; $D$ is the photon diffusion
coefficient defined by $D:=1/(3\mu_{s}^{\prime})$, where $\mu_{s}^{\prime}>0$
is the reduced scattering coefficient. The DOT-FDOT in this work is formulated
as follows.

Inverse problem: Let $\Gamma_{\rm in}$ and $\Gamma$ be open sub-boundaries of
$\partial\Omega$ satisfying

$\Gamma_{\rm in}\cup\Gamma=\partial\Omega,\quad\Gamma_{\rm
in}\cap\Gamma\neq\emptyset.$ (1.3)

Suppose $D,\mu_{a}$ and $\mu_{f}$ are unknown in (1.1)-(1.2). By a suitable
choice of the input $g$ supported on $\Gamma_{\rm in}\times[0,T]$, we consider
to recover the $(D,\mu_{a},\mu_{f})$ simultaneously from the time-dependent
measurements on the sub-boundary

$D\frac{\partial u_{e}}{\partial\nu}|_{\Gamma\times(0,T)},\quad
D\frac{\partial u_{m}}{\partial\nu}|_{\Gamma\times(0,T)},$ (1.4)

where $\nu=\nu(x)$ is the outward unit normal vector at $x\in\Gamma$.

For these optical imaging problems, it is also possible to collect the
observation data in the frequency domain at one or several modulation
frequencies [3]. Further, to overcome the ill-posedness of the inverse
problems, the necessary regularization techniques such as Tikhonov
regularization, sparse regularization methods and hybrid regularization
methods have been introduced [12, 15, 28]. 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. Above all, there already many works focus on
improving the quality of the reconstruction and the reconstruction speed.
However, in the aspect of theoretical analysis, especially the uniqueness of
DOT-FDOT has not been rigorously investigated as far as we know. In this work,
we will study the uniqueness of time-domain DOT-FDOT. The details are
contained in the main theorem, which is stated in the next subsection.

###  1.2 The main theorem.

In this work, the so-called time-domain DOT-FDOT based on the initial-boundary
value problems (1.1)-(1.2) is to solve the inverse problem:

$\text{recovering}\ \ (D,\mu_{a},\mu_{f})\ \ \text{from the measurements
(\ref{measure})}.$ (1.5)

The DOT-FDOT of (1.5) is equivalent to the two successive inverse problems
given by

$\text{recovering}\ D\ \text{and}\ (\mu_{a}+\mu_{f})\ \text{from the
measurement}\ D\frac{\partial u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$ (1.6)

and

$\text{recovering}\ \mu_{a}\ and\ \mu_{f}\ \text{from the measurement}\
D\big{(}\frac{\partial u_{e}}{\partial\nu}+\frac{\partial
u_{m}}{\partial\nu}\big{)}|_{\Gamma\times(0,T)}.$ (1.7)

To make our analysis more convenient, we set $c=1$ in models (1.1)-(1.2). Now
it is time to state the uniqueness theorem.

######  Theorem 1.

Let $\Gamma_{\rm in}$ and $\Gamma$ be nonempty open sub-boundaries of
$\partial\Omega$, respectively, and fulfill (1.3). For the two sets of
unknowns $\big{\\{}D,\mu_{a},\mu_{f}\big{\\}}$ and
$\big{\\{}\tilde{D},\tilde{\mu}_{a},\tilde{\mu}_{f}\big{\\}}$, by
$\\{u_{e},u_{m}\\}$ and $\\{\tilde{u}_{e},\tilde{u}_{m}\\}$ we respectively
denote the corresponding solutions of initial-boundary value problems
(1.1)-(1.2) with boundary input $g$ satisfying (2.3). Then we conclude that:

  * (i)

$D\frac{\partial
u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}=\tilde{D}\frac{\partial\tilde{u}_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$
with a priori $\nabla D=\nabla\tilde{D}$ on $\partial\Omega$ implies

$D=\tilde{D},\quad\mu_{a}+\mu_{f}=\tilde{\mu}_{a}+\tilde{\mu}_{f};$ (1.8)

  * (ii)

with result (1.8) and $D\frac{\partial
u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}=\tilde{D}\frac{\partial\tilde{u}_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$,
$D\frac{\partial
u_{m}}{\partial\nu}|_{\Gamma\times(0,T)}=\tilde{D}\frac{\partial\tilde{u}_{m}}{\partial\nu}|_{\Gamma\times(0,T)}$,
we have

$\mu_{a}=\tilde{\mu}_{a},\quad\mu_{f}=\tilde{\mu}_{f}.$ (1.9)

######  Remark 1.1.

If $\mu_{f}=0$ in the initial boundary value problem (1.1), the
reconstructions of $\mu_{a}$ and $D$ from the sub-boundary measurement
$D\frac{\partial u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$ are known as the
DOT. Given $(D,\mu_{a})$ and suppose $\mu_{f}=0$ in (1.1), then the inverse
problem of recovering $\mu_{f}$ from the data $D\frac{\partial
u_{m}}{\partial\nu}|_{\Gamma\times(0,T)}$ is the linear FDOT, which is the
common problem in the community of fluorescence tomography. Hence, Theorem (1)
covers the uniqueness result of these two common optical imaging problems but
needs the both boundary flux of excitation $u_{e}$ and fluorescence $u_{m}$.

###  1.3 Contribution and outline.

The literature on the uniqueness of time-domain DOT-FDOT is relatively rare
and the difficulty to prove the uniqueness of such inverse problem is
addressed in the overview [3], especially the identifiability of multiple
coefficients from the boundary measurement. Some aspects of uniqueness and
ill-posedness of DOT are considered [1, 2, 4, 6, 9, 11, 13, 18, 34], and the
uniqueness results are related with the applied DOT or FDOT models. It is in
general not possible to simultaneously recover the diffusion and absorption
coefficients in steady-state (dc) DOT [4]. Still for steady-state DOT, if
restricting the unknown parameters to piecewise constant diffusions and
piecewise analytic absorption coefficients, both parameters can simultaneously
be determined from complete measurement data on an arbitrarily small part of
the boundary [18]. Furthermore, although DOT is a nonlinear inverse problem, a
linearized version is often considered. The uniqueness of linearized time-
domain DOT is considered [34]. For the nonliner time-domain DOT, the
uniqueness can be regained but need the Dirichlet-to-Neumann (DtN) mapping
data on the sub-boundary [9, 13]. For linear time-domain FDOT, an
identifiability result of $\mu_{f}$ is established in [27], but it requires
strong prior assumptions on $\mu_{f}$. Recently, the uniqueness of identifying
unknown $\mu_{f}$ from arbitrary sub-boundary data is established [37] by
assuming the semi-discrete form of right-hand source. However, we see that the
known uniqueness results require either the linearized simple model or the
impractical DtN mapping measurement.

In this work, the uniqueness of time-domain DOT-FDOT in a more general
framework is established. The inverse problem under consideration is more
complicated and challenged in the following sense

  * (i)

The initial-boundary value problems (1.1)-(1.2) are coupled and the excitation
$u_{e}$ depends on $\mu_{f}$ nonlinearly. Therefore the concerned inverse
problems (1.6) and (1.7) are nonlinear.

  * (ii)

The problem of recovering the parameters $D,\mu_{a}$ and $\mu_{f}$ is a
"mixed" inverse problem. The used data is one single measurement on the sub-
boundary $\Gamma\times(0,T)$, not the DtN mapping data.

Theorem (1) confirms that the unknowns $(D,\mu_{a},\mu_{f})$ can be uniquely
determined by the arbitrary sub-boundary measurements $D\frac{\partial
u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$ and $D\frac{\partial
u_{m}}{\partial\nu}|_{\Gamma\times(0,T)}$. The observation area $\Gamma$ can
be an arbitrarily open subset of boundary $\partial\Omega$. Moreover, Theorem
(1) covers the uniqueness results of the common DOT problem and FDOT problem
as mentioned in Remark (1.1).

The rest of this article is organized as follows. In Section (2), we collect
several preliminary works. We prove the uniqueness theorem in Section (3). The
numerical inversions will be considered in Section (4). We show the validity
of the proposed scheme by several numerical examples. Finally, Section (5) is
devoted to conclusion and future works.

##  2 Some preliminary results.

In this section, we first fix some general settings and notations, and show
several lemmas which will be used in the forthcoming discussion.

###  2.1 Settings.

Now we give several settings. We first define the operator $A_{e}$ by

$A_{e}u:=-\text{div}\left(D(x)\nabla u\right)+(\mu_{a}(x)+\mu_{f}(x))u,\quad
u\in D(A):=H^{2}(\Omega)\cap H_{0}^{1}(\Omega).$

Since $D(x)>0$ and $\mu_{e}+\mu_{f}\geq 0$, the operator $A_{e}$ is symmetric
uniformly elliptic. Let $\\{\lambda_{n},\varphi_{n}\\}_{n=1}^{\infty}$ be the
Dirichlet eigensystem of the operator $A_{e}$, where
$0<\lambda_{1}\leq\lambda_{2}\leq\cdots$ and the corresponding eigenfunctions
$\\{\varphi_{n}\\}_{n=1}^{\infty}$ forms an orthonormal basis of
$L^{2}(\Omega)$.

If considering the algebraic multiplicity, we need to make slight changes to
the above symbols. Indeed, by $\\{\lambda_{k}\\}_{k\in\mathbb{N}}$ and
$m_{k}\in\mathbb{N}$ we denote the strictly increasing sequence of the
eigenvalues of $A_{e}$ and the algebraic multiplicity of $\lambda_{k}$,
respectively. For each eigenvalue $\lambda_{k}$, we introduce a family
$\\{\varphi_{k,\ell}\\}_{\ell=1}^{m_{k}}$ of eigenfunctions of $A_{e}$, i.e.,

$A_{e}\varphi_{k,\ell}=\lambda_{k}\varphi_{k,\ell},\quad\ell=1,\ldots,m_{k},$

which forms an orthonormal basis in $L^{2}(\Omega)$ of the algebraic
eigenspace of the operator $A_{e}$ associated with $\lambda_{k}$.

Next we specify the choice of the Dirichlet boundary input $g(x,t)$, which
plays an essential role in the consideration of our inverse problem. Following
the settings in the paper from Kian et al. [26], we choose a sequence of real
and nonnegative functions $\\{\psi_{k}\\}_{k\in\mathbb{N}}$ define on
$\overline{\mathbb{R}_{+}}:=[0,\infty)$ such that $0\leq\psi_{k}\leq 1$ and

$C^{\infty}(\overline{\mathbb{R}_{+}})\ni\psi_{k}=\begin{cases}0&\mbox{on
}[0,t_{2k-2}],\\\ 1&\mbox{on }[t_{2k-1},\infty),\end{cases}$ (2.1)

where $\\{t_{k}\\}_{k=0}^{\infty}$ is a strictly increasing sequence
satisfying $t_{0}=0$ and $\lim_{k\to\infty}t_{k}=T$. We fix also
$\\{c_{k}\\}_{k=0}^{\infty}$ a sequence of $\mathbb{R}_{+}:=(0,\infty)$ such
that

$\sum_{k=1}^{\infty}c_{k}\|\psi_{k}\|_{W^{2,\infty}(\mathbb{R}_{+})}<\infty.$
(2.2)

For example, one can choose $c_{k}:=e^{-k^{2}}$. We select a sequence of
functions $\\{\eta_{k}\\}_{k\in\mathbb{N}}$ of
$H^{\frac{3}{2}}(\partial\Omega)$ such that
$\|\eta_{k}\|_{H^{\frac{3}{2}}(\partial\Omega)}\leq 1$ and
span$\\{\eta_{k}\\}$ is dense in $H^{\frac{3}{2}}(\partial\Omega)$. Finally,
we can construct the Dirichlet boundary input $g$ of the system (1.1)

$g(x,t):=\chi(x)\sum_{k=1}^{\infty}c_{k}\,\psi_{k}(t)\eta_{k}(x),\quad(x,t)\in\partial\Omega\times\mathbb{R}_{+},$
(2.3)

where $0\leq\chi\leq 1$ is such that

$C_{0}^{\infty}(\partial\Omega)\ni\chi=\begin{cases}1&\mbox{ on }\Gamma_{\rm
in}^{\prime},\\\ 0&\mbox{ on }\partial\Omega\setminus\Gamma_{\rm
in},\end{cases}$ (2.4)

where $\Gamma_{\rm in}^{\prime}\subset\Gamma_{\rm in}\subset\partial\Omega$.
Based on the above settings of boundary input $g$, it is not difficult to
check that the function $g$ belongs to the function space
$C^{2}(\overline{\mathbb{R}_{+}};H^{\frac{3}{2}}(\partial\Omega))$.

###  2.2 Time-analyticity of solutions.

Now based on the characteristics of boundary condition $g$ in (2.3), we
construct the following subproblems corresponding to the system (1.1) by
treating each term of the series as a new boundary input, that is,

$\left\\{\begin{aligned} \partial_{t}u_{e,k}+A_{e}u_{e,k}&=0&\quad&\mbox{in
}\Omega\times\mathbb{R}_{+},\\\ u_{e,k}&=0&\quad&\mbox{in
}\Omega\times\\{0\\},\\\ u_{e,k}&=\chi
c_{k}\,\psi_{k}\,\eta_{k}&\quad&\mbox{on
}\partial\Omega\times\mathbb{R}_{+}.\end{aligned}\right.$ (2.5)

We investigate the analyticity of the solutions $u_{e,k}$ to (2.5) in time for
each $k\in\mathbb{N}$. To this end, by $\mathcal{H}(r,\infty;X)$ we denote the
set of holomophic functions on $(r,+\infty)$ taking values in the Banach space
$X$. Then we state the following intermediate result.

######  Lemma 2.1.

For any integer $k\in\mathbb{N}$, the solution $u_{e,k}$ to the sub-problem
(2.5) belongs to $L^{\infty}(\mathbb{R}_{+};H^{2}(\Omega))$ and is such that
$u_{e,k}|_{(t_{2k-1},\infty)}\in\mathcal{H}(t_{2k-1},\infty;H^{2}(\Omega))$.

###### Proof.

Fix $G_{k}$ solving the boundary value problem

$\left\\{\begin{array}[]{rcll}A_{e}G_{k}&=&0&\mbox{in }\Omega,\\\
G_{k}&=&\chi\eta_{k}&\mbox{on}\ \partial\Omega.\end{array}\right.$

It is not difficult to check from the theory of elliptic equations that
$G_{k}\in H^{2}(\Omega)$ and $G_{k}$ is such that the following regularity
estimate

$\|G_{k}\|_{H^{2}(\Omega)}\leq
C\|\eta_{k}\|_{H^{\frac{3}{2}}(\partial\Omega)}\leq C.$ (2.6)

Here and henceforth, by $C>0$ we denote a general constant which is
independent of $k$ and may change by lines. Then, we can decompose $u_{e,k}$
into $u_{e,k}=w_{e,k}+v_{e,k}$ with $v_{e,k}(t,\cdot)=c_{k}\psi_{k}(t)G_{k}$
and $w_{e,k}$ solving

$\left\\{\begin{aligned}
(\partial_{t}+A_{e})w_{e,k}&=-c_{k}\psi_{k}^{\prime}(t)G_{k}&\quad&\mbox{in
}\Omega\times\mathbb{R}_{+},\\\ w_{e,k}(x,t)&=0&\quad&\mbox{on
}\partial\Omega\times\mathbb{R}_{+},\\\ w_{e,k}(\cdot,0)&=0&\quad&\mbox{in
}\Omega.\end{aligned}\right.$

It is clear that $v_{e,k}\in\mathcal{H}(t_{2k-1},\infty;H^{2}(\Omega))$ by the
notation of (2.1). Therefore, we only need to prove that
$w_{e,k}\in\mathcal{H}(t_{2k-1},\infty;H^{2}(\Omega))$. For this, by noticing
the argument of eigenfunction expansion, we can arrive at the representation
formula of the solution $w_{e,k}$ as follows.

$w_{e,k}(\cdot,t)=-c_{k}\sum_{\ell=1}^{\infty}\int_{0}^{t}\exp(-\lambda_{\ell}(t-s))\psi_{k}^{\prime}(s)ds\left\langle
G_{k},\varphi_{\ell}\right\rangle_{L^{2}(\Omega)}\varphi_{\ell}.$ (2.7)

Using the fact that $\psi_{k}^{\prime}=0$ on $(t_{2k-1},+\infty)$, it follows
that

$w_{e,k}(\cdot,t)=-c_{k}\sum_{\ell=1}^{\infty}e^{-\lambda_{\ell}t}\int_{0}^{t_{2k-1}}e^{\lambda_{\ell}s}\psi_{k}^{\prime}(s)ds\left\langle
G_{k},\varphi_{\ell}\right\rangle_{L^{2}(\Omega)}\varphi_{\ell},\quad t\geq
t_{2k-1},$

from which by further noting the analyticity of the exponential function we
conclude that
$w_{e,k}|_{(t_{2k-1},\infty)}\in\mathcal{H}(t_{2k-1},+\infty;H^{2}(\Omega))$.
Therefore, $u_{e,k}$ restricted to $(t_{2k-1},\infty)$ belongs to
$\mathcal{H}(t_{2k-1},\infty;H^{2}(\Omega))$.

We next show that the solution $u_{e,k}\in
L^{\infty}(\mathbb{R}_{+};H^{2}(\Omega))$. It is not difficult to check that
$v_{e,k}\in L^{\infty}(\mathbb{R}_{+};H^{2}(\Omega))$, so it remains to deal
with $w_{e,k}$. For this, on the basis of the representation formula (2.7) and
noting $\|u\|_{H^{2}(\Omega)}\leq C\|A_{e}u\|_{L^{2}(\Omega)}$, by a directly
calculation we arrive at the inequalities

$\displaystyle\|w_{e,k}(\cdot,t)\|_{H^{2}(\Omega)}^{2}\leq
C\|A_{e}w_{e,k}(\cdot,t)\|_{L^{2}(\Omega)}^{2}$ $\displaystyle\leq
C|c_{k}|^{2}\sum_{\ell=1}^{\infty}\lambda_{\ell}^{2}\left|\int_{0}^{t}\exp\left\\{-\lambda_{\ell}(t-s)\right\\}\psi_{k}^{\prime}(s)ds\right|^{2}\left|\langle
G_{k},\varphi_{\ell}\rangle_{L^{2}(\Omega)}\right|^{2}$ $\displaystyle\leq
C|c_{k}|^{2}\|\psi_{k}\|_{W^{1,\infty}(\mathbb{R}_{+})}^{2}\sum_{\ell=1}^{\infty}\lambda_{\ell}^{2}\left|\int_{0}^{t}e^{-\lambda_{\ell}s}ds\right|^{2}\left|\langle
G_{k},\varphi_{\ell}\rangle_{L^{2}(\Omega)}\right|^{2}$
$\displaystyle=C|c_{k}|^{2}\|\psi_{k}\|_{W^{1,\infty}(\mathbb{R}_{+})}^{2}\sum_{\ell=1}^{\infty}\left|1-e^{-\lambda_{\ell}t}\right|^{2}\left|\langle
G_{k},\varphi_{\ell}\rangle_{L^{2}(\Omega)}\right|^{2}\leq
C\|G_{k}\|_{L^{2}(\Omega)}^{2}\leq C.$

Here the last inequality is due to the assumption (2.2) and the regularity
estimate (2.6). Therefore we have the function $t\longmapsto w_{e,k}(\cdot,t)$
is in $L^{\infty}(\mathbb{R}_{+};H^{2}(\Omega))$. Collecting all the above
estimates, we can finish the proof of the lemma. ∎

##  3 The proof of the uniqueness theorem.

As above mentioned, the unique solvability of our DOT-FDOT problem (1.5) is
equivalent to the uniqueness of the two sub inverse problems: (1.6) and (1.7).
Consequently, we first treat the sub-problem (1.6), that is, the determination
of the diffusion coefficient $D$ and mixture absorption $\mu_{a}+\mu_{f}$ in
(1.1). Then we proceed to the uniqueness of the reconstruction of background
absorption $\mu_{a}$ and fluorescence absorption $\mu_{f}$ (i.e, the sub-
problem (1.7)) by employing the similar argument used in solving (1.6).

###  3.1 Recovering the diffusion coefficient and mixture absorption.

In this part, we consider the determination of the diffusion coefficient $D$
and the mixture absorption $\mu_{a}+\mu_{f}$ from the boundary measurement

$D(x)\frac{\partial u_{e}}{\partial\nu}(x,t),\quad(x,t)\in\Gamma\times(0,T).$
(3.1)

Firstly, we transfer the observation information

$D\frac{\partial
u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}=\tilde{D}\frac{\partial\tilde{u}_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$
(3.2)

of original problem (1.1) to a sequential identities of Neumann boundary of
the solutions to the corresponding subproblems (2.5), that is, we have the
following useful lemma.

######  Lemma 3.1.

The assumption (3.2) implies that

$D\frac{\partial
u_{e,k}}{\partial\nu}=\tilde{D}\frac{\partial\tilde{u}_{e,k}}{\partial\nu}\quad\mbox{on
}\Gamma\times(0,\infty),\ k\in\mathbb{N},$ (3.3)

where $u_{e,k}$ and $\tilde{u}_{e,k}$ are solutions of (2.5) with respect to
$(D,\mu_{a}+\mu_{f})$ and $(\tilde{D},\tilde{\mu}_{a}+\tilde{\mu}_{f})$.

###### Proof.

We will prove the conclusion (3.3) by applying inductive argument. Firstly,
for $k=1$, from the choice of the sequence $\\{\psi_{k}\\}_{k\in\mathbb{N}}$
in (2.1), it is not difficult to see that

$\psi_{k}=0\quad\mbox{in }(0,t_{2}),\quad\forall\,k\geq 2.$

Therefore, in the case of $t\in(0,t_{2})$, the boundary condition $g$ can be
deduce into

$g(x,t)=c_{1}\psi_{1}(t)\chi(x)\eta_{1}(x),\quad(x,t)\in\partial\Omega\times(0,t_{2}).$

Now we conclude from the wellposeness of the problem (1.1) that $u=u_{e,1}$
and $\tilde{u}=\tilde{u}_{e,1}$ in $\Omega\times(0,t_{2})$, which further
combined with the condition (3.2) implies

$D\frac{\partial
u_{e,1}}{\partial\nu}=\tilde{D}\frac{\partial\tilde{u}_{e,1}}{\partial\nu}\quad\mbox{on
}\Gamma\times(0,t_{2}).$ (3.4)

On the other hand, from Lemma (2.1) we know that
$u_{e,1},\tilde{u}_{e,1}\in\mathcal{H}(t_{1},\infty;H^{2}(\Omega))$. Thus, the
trace theorem and properties of the analytic functions yield that the
following two functions

$t\longmapsto D\frac{\partial u_{e,1}}{\partial\nu}(\cdot,t)|_{\Gamma}\mbox{
and
}t\longmapsto\tilde{D}\frac{\partial\tilde{u}_{e,1}}{\partial\nu}(\cdot,t)|_{\Gamma},\quad
t>t_{1}$

belong to $\mathcal{H}(t_{1},\infty;L^{2}(\Gamma))$, and the condition (3.4)
implies (3.3) for $k=1$. We finish the proof of (3.3) with $k=1$. Next
assuming that (3.3) is fulfilled for all $k=1,2,\ldots,\ell$ with some
$\ell\in\mathbb{N}$, we will show that the conclusion (3.3) is valid also for
$k=\ell+1$. For this, again noting the definition of $\psi_{k}$ in (2.1), it
can be easily seen that

$\psi_{k}=0\quad\mbox{in }(0,t_{2\ell+2}),\quad\forall\,k\geq\ell+2.$

Therefore, for $t\in(0,t_{2\ell+2})$, the Dirichlet boundary input $g$ can be
rewritten as follows

$g(x,t)=\chi(x)\sum_{k=1}^{\ell+1}c_{k}\psi_{k}(t)\eta_{k}(x),\quad(x,t)\in\partial\Omega\times(0,t_{2\ell+2}).$

By an argument similar to the case of $k=1$, we can obtain

$\sum_{k=1}^{\ell+1}u_{e,k}=u_{e},\quad\sum_{k=1}^{\ell+1}\tilde{u}_{e,k}=\tilde{u}_{e}\quad\mbox{in
}\Omega\times(0,t_{2\ell+2}).$

Therefore, (3.2) implies

$D\sum_{k=1}^{\ell+1}\frac{\partial
u_{e,k}}{\partial\nu}=\tilde{D}\sum_{k=1}^{\ell+1}\frac{\partial\tilde{u}_{e,k}}{\partial\nu}\quad\mbox{on
}\Gamma\times(0,t_{2\ell+2}).$

Now by the inductive assumption, that is, (3.3) holds for $k=1,\cdots,\ell$,
we further obtain that

$D\frac{\partial
u_{e,\ell+1}}{\partial\nu}=\tilde{D}\frac{\partial\tilde{u}_{e,\ell+1}}{\partial\nu}\quad\mbox{on
}\Gamma\times(0,t_{2\ell+2}).$

Again by the use of the $t$-analyticity result in Lemma (2.1), we conclude
that the following two functions are such that

$t\longmapsto D\frac{\partial
u_{e,\ell+1}}{\partial\nu}(\cdot,t)|_{\Gamma}\mbox{ and
}t\longmapsto\tilde{D}\frac{\partial\tilde{u}_{e,\ell+1}}{\partial\nu}(\cdot,t)|_{\Gamma}\in\mathcal{H}(t_{2\ell+1},\infty;L^{2}(\Gamma)),$

from which we further use the continuation principle of analytic functions to
imply (3.3) holds true for $k=\ell+1$. This proves that (3.3) is valid for all
$k\in\mathbb{N}$. We complete the proof of the lemma. ∎

######  Remark 3.2.

Here and henceforth, some parts of the proofs may be similar to the argument
used in Kian et al. [26]. But we retain this proof for the completeness and
consistency of the article because our paper deals with a totally different
model.

From the results in Lemma (2.1) and Sobolev trace theorem, we can directly
deduce that

######  Corollary 3.3.

For any $k\in\mathbb{N}$ and any $\xi>0$, we have

$t\longmapsto e^{-\xi t}u_{e,k}(\cdot,t)\in
L^{1}(\mathbb{R}_{+};H^{2}(\Omega))$

and

$t\longmapsto e^{-\xi t}\frac{\partial u_{e,k}}{\partial\nu}(\cdot,t)\in
L^{1}(\mathbb{R}_{+};L^{2}(\Gamma)).$

The above statement ensures the existence of the Laplace transforms of the two
functions $u_{e,k}$ and $\tilde{u}_{e,k}$, which allows us to use Laplace
transform to consider our inverse problem in the frequency domain. Firstly,
for any $\xi>0$ we introduce the following boundary value problem for elliptic
equation

$\left\\{\begin{aligned} (A_{e}+\xi)U_{\xi}(x)&=0&\quad&\mbox{in }\Omega,\\\
U_{\xi}(x)&=\chi(x)h(x)&\quad&\mbox{on }\partial\Omega,\end{aligned}\right.$
(3.5)

where $h\in H^{\frac{3}{2}}(\partial\Omega)$ and the cut-off function $\chi$
is defined in (2.4). We first assert that the function $U_{\xi}$ is analytic
with respect to the parameter $\xi>0$ and can be analytically extended to the
set $\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$. We have the
following lemma.

######  Lemma 3.4.

The function $\mathbb{R}_{+}\ni\xi\mapsto U_{\xi}\in H^{2}(\Omega)$ is
analytic with respect to $\xi>0$ and can be analytically extended to
$\xi\in\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$.

###### Proof.

Since $h\in H^{\frac{3}{2}}(\partial\Omega)$, noting the Sobolev trace
theorem, we can find a function $H\in H^{2}(\Omega)$ such that
$H|_{\partial\Omega}=-\chi h$ in the trace sense and satisfies the estimate

$\|H\|_{H^{2}(\Omega)}\leq C\|h\|_{H^{\frac{3}{2}}(\partial\Omega)}.$ (3.6)

Letting $W_{\xi}:=U_{\xi}+H$, it is not difficult to see that the function
$W_{\xi}$ reads the following homogeneous boundary value problem

$\left\\{\begin{array}[]{rcll}(A_{e}+\xi)W_{\xi}&=&(A_{e}+\xi)H&\mbox{in
}\Omega,\\\ W_{\xi}&=&0&\mbox{on}\ \partial\Omega.\end{array}\right.$

Now from the argument of the Fourier series expansion, we assert that the
function $W_{\xi}$ can be represented by

$W_{\xi}=\sum_{k=1}^{\infty}\sum_{\ell=1}^{m_{k}}\langle
W_{\xi},\varphi_{k,\ell}\rangle_{L^{2}(\Omega)}\varphi_{k,\ell}=\sum_{k=1}^{\infty}\frac{1}{\xi+\lambda_{k}}\sum_{\ell=1}^{m_{k}}\big{\langle}(A_{e}+\xi)H,\varphi_{k,\ell}\big{\rangle}_{L^{2}(\Omega)}\varphi_{k,\ell},$
(3.7)

which combined with the relation $W_{\xi}=U_{\xi}+H$ gives the representation
formula of the solution $U_{\xi}$:

$U_{\xi}=\sum_{k=1}^{\infty}\sum_{\ell=1}^{m_{k}}\langle
U_{\xi},\varphi_{k,\ell}\rangle_{L^{2}(\Omega)}\varphi_{k,\ell}=\sum_{k=1}^{\infty}\frac{1}{\xi+\lambda_{k}}\sum_{\ell=1}^{m_{k}}\big{\langle}(A_{e}+\xi)H,\varphi_{k,\ell}\big{\rangle}_{L^{2}(\Omega)}\varphi_{k,\ell}-H.$
(3.8)

We next show that the function $\mathbb{R}_{+}\ni\xi\mapsto W_{\xi}\in
H^{2}(\Omega)$) can be analytically extended to the complex domain
$\xi\in\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$. For this, fixing
any compact subset $K$ of
$\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$, on the basis of the
representation formula (3.8) and the estimate $\|u\|_{H^{2}(\Omega)}\leq
C\|A_{e}u\|_{L^{2}(\Omega)}$, by a directly calculation we arrive at the
inequalities

$\displaystyle\|W_{\xi}(\cdot\,)\|_{H^{2}(\Omega)}^{2}$ $\displaystyle\leq
C\|A_{e}W_{\xi}(\cdot\,)\|_{L^{2}(\Omega)}^{2}$ $\displaystyle\leq
C\sum_{k=1}^{\infty}\frac{\lambda_{k}^{2}}{(\xi+\lambda_{k})^{2}}\left|\sum_{\ell=1}^{m_{k}}\langle(A_{e}+\xi)H,\varphi_{k,\ell}\rangle_{L^{2}(\Omega)}\right|^{2}.$

Since $K\subset\subset\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$,
we see that $\frac{\lambda_{k}^{2}}{(\xi+\lambda_{k})^{2}}\leq C_{K}$, where
the constant $C_{K}$ denotes a general positive constant which is independent
of $k$ and may change by lines, which further implies that

$\displaystyle\|W_{\xi}(\cdot\,)\|_{H^{2}(\Omega)}^{2}\leq
C_{K}\|(A_{e}+\xi)H\|_{L^{2}(\Omega)}^{2}\leq
C_{K}\|h\|_{H^{\frac{3}{2}}(\partial\Omega)}.$

Here the last inequality is due to the estimate (3.6). Consequently, we assert
that the series in (3.7) is uniformly convergent for any $\xi\in K$. Moreover,
noting that the compact subset $K$ can be arbitrarily chosen in
$\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$, we have the required
conclusion of the function $W_{\xi}$. Finally, since $H$ does not depend on
$\xi$, we obtain that the map $\mathbb{R}_{+}\ni\xi\mapsto U_{\xi}\in
H^{2}(\Omega)$ can be analytical extended to the set
$\mathbb{C}\setminus\\{-\lambda_{k}\\}_{k=1}^{\infty}$. We finish the proof of
the lemma. ∎

We associate the problem (3.5) with the DtN mapping
$\Lambda_{\xi}:H^{\frac{3}{2}}(\partial\Omega)\to L^{2}(\Gamma)$ which is
defined by

$\Lambda_{\xi}:h\longmapsto D\frac{\partial
U_{\xi}}{\partial\nu}|_{\Gamma},\quad\xi>0.$

Similarly, we can define the DtN map

$\tilde{\Lambda}_{\xi}:h\in
H^{\frac{3}{2}}(\partial\Omega)\longmapsto\tilde{D}\frac{\partial\tilde{U}_{\xi}}{\partial\nu}|_{\Gamma}$

corresponding to the elliptic operator $\tilde{A}_{e}$:

$\tilde{A}_{e}:=-\text{div}(\tilde{D}\nabla)+\tilde{\mu}_{a}+\tilde{\mu}_{f}.$

On the basis of the results in Lemma (3.3), we can obtain the following lemma.

######  Lemma 3.5.

The condition (3.3) implies $\Lambda_{\xi}=\tilde{\Lambda}_{\xi}$ for any
$\xi>0$.

###### Proof.

Applying the Laplace transform $\mathcal{L}$ in time on both sides of (3.3)
implies

$\Lambda_{\xi}(c_{k}\eta_{k}\mathcal{L}\psi_{k})=\tilde{\Lambda}_{\xi}(c_{k}\eta_{k}\mathcal{L}\psi_{k}),\quad\xi>0,\
k\in\mathbb{N}.$

Moreover, in view of the assumption $\psi_{k}\geq 0,\not\equiv 0$, we see that
the Laplace transform $(\mathcal{L}\psi_{k})(\xi)>0$ is valid for all $\xi>0$,
which combined with the fact $c_{k}>0$ and the linearity of $\Lambda_{\xi}$
and $\tilde{\Lambda}_{\xi}$ implies

$\Lambda_{\xi}h=\tilde{\Lambda}_{\xi}h,\quad\forall\xi>0,\ \forall\,h\in{\rm
span}\\{\eta_{k}\\}_{k\in\mathbb{N}}.$

Finally, $\Lambda_{\xi}=\tilde{\Lambda}_{\xi}$ is valid for any $\xi>0$ due to
the density of ${\rm span}\\{\eta_{k}\\}$ in
$H^{\frac{3}{2}}(\partial\Omega)$. We finish the proof of the lemma. ∎

Recalling the notations of the Dirichlet eigensystem
$\\{\lambda_{k},\\{\varphi_{k,\ell}\\}_{\ell=1}^{m_{k}}\\}_{k=1}^{\infty}$ of
the elliptic operator $A_{e}$, we define the sequence of functions:

$\Theta_{k}(x,x^{\prime}):=D(x)D(x^{\prime})\sum_{\ell=1}^{m_{k}}\frac{\partial\varphi_{k,\ell}}{\partial\nu}(x)\frac{\partial\varphi_{k,\ell}}{\partial\nu}(x^{\prime}),\quad
x,x^{\prime}\in\partial\Omega,\ k\in\mathbb{N}.$

######  Lemma 3.6.

The observation data (3.2) implies $\lambda_{k}=\tilde{\lambda}_{k}$ and

$\Theta_{k}(x,x^{\prime})=\tilde{\Theta}_{k}(x,x^{\prime}),\quad x\in\Gamma,\
x^{\prime}\in\Gamma_{\rm in}^{\prime},$ (3.9)

where $k\in\mathbb{N}$ and
$\\{\lambda_{k},\\{\varphi_{k,\ell}\\}_{\ell=1}^{m_{k}}\\}_{k=1}^{\infty}$ and
$\\{\tilde{\lambda}_{k},\\{\tilde{\varphi}_{k,\ell}\\}_{\ell=1}^{\tilde{m}_{k}}\\}_{k=1}^{\infty}$
are the Dirichlet eigensystems of the elliptic operators $A_{e}$ and
$\tilde{A}_{e}$.

###### Proof.

In the first step, we will give other type of representation formula in (3.8).
Taking the scalar product on both sides of (3.5) with eigenfunction of
$\varphi_{k,\ell}$, and by integration by parts we can see that the Fourier
coefficient $\langle U_{\xi},\varphi_{k,\ell}\rangle_{L^{2}(\Omega)}$ in (3.8)
also reads

$\langle
U_{\xi},\varphi_{k,\ell}\rangle_{L^{2}(\Omega)}=-\frac{1}{\xi+\lambda_{k}}\big{\langle}\,\chi
h,D\frac{\partial\varphi_{k,\ell}}{\partial\nu}\big{\rangle}_{L^{2}(\partial\Omega)},$

which further implies the representation of $U_{\xi}$ as follows

$U_{\xi}=-\sum_{k=1}^{\infty}\frac{1}{\xi+\lambda_{k}}{\sum_{\ell=1}^{m_{k}}\big{\langle}\,\chi
h,D\frac{\partial\varphi_{k,\ell}}{\partial\nu}\big{\rangle}_{L^{2}(\partial\Omega)}\varphi_{k,\ell}}.$

On the other hand, in view of the result in Lemma (3.5), we have

$D\frac{\partial
U_{\xi}}{\partial\nu}\big{|}_{\Gamma}=\tilde{D}\frac{\partial\tilde{U}_{\xi}}{\partial\nu}\big{|}_{\Gamma},\quad\xi>0,\
h\in H^{\frac{3}{2}}(\partial\Omega).$

Now from the definition of the functions $\Theta_{k}$ and
$\tilde{\Theta}_{k}$, noting the notation of the scalar product
$\langle\cdot,\cdot\rangle_{L^{2}(\Omega)}$, we can rephrase the above
identity as follows

$\displaystyle\sum_{k=1}^{\infty}\frac{1}{\xi+\lambda_{k}}\int_{\partial\Omega}\Theta_{k}(x,x^{\prime})\chi(x^{\prime})h(x^{\prime})dx^{\prime}=\sum_{k=1}^{\infty}\frac{1}{\xi+\tilde{\lambda}_{k}}\int_{\partial\Omega}\tilde{\Theta}_{k}(x,x^{\prime})\chi(x^{\prime})h(x^{\prime})dx^{\prime},\quad
x\in\Gamma.$

In view of Lemma (3.4), the above series are uniformly convergent for $\xi\in
K$, where $K$ is any compact subset of
$\mathbb{C}\setminus\\{-\lambda_{k},-\tilde{\lambda}_{k}|\ k\in\mathbb{N}\\}$,
therefore from the unique continuation of analytic functions, it follows that

$\displaystyle\sum_{k=1}^{\infty}\frac{1}{z+\lambda_{k}}\int_{\partial\Omega}\Theta_{k}(x,x^{\prime})\chi(x^{\prime})h(x^{\prime})dx^{\prime}=\sum_{k=1}^{\infty}\frac{1}{z+\tilde{\lambda}_{k}}\int_{\partial\Omega}\tilde{\Theta}_{k}(x,x^{\prime})\chi(x^{\prime})h(x^{\prime})dx^{\prime}$
(3.10)

for any $x\in\Gamma$ and
$z\in\mathbb{C}\setminus\\{-\lambda_{k},-\tilde{\lambda}_{k}\\}_{k=1}^{\infty}$.
Now one can show that $\lambda_{k}=\tilde{\lambda}_{k}$, $k\in\mathbb{N}$ in a
similar way to [36, Theorem 4.2]. Indeed, if not, we assume there exists
$k_{0}\in\mathbb{N}$ such that $-\lambda_{k_{0}}\neq\tilde{\lambda}_{k_{0}}$,
then we can take a sufficiently small disk which includes $-\lambda_{k_{0}}$
and does not include
$\\{-\tilde{\lambda}_{k}\\}_{k=1}^{\infty}\cup\\{-\lambda_{k}\\}_{k\neq
k_{0}}$. Integrating (3.10) in this disk, we see that

$\displaystyle\int_{\partial\Omega}\Theta_{k_{0}}(x,x^{\prime})\chi(x^{\prime})h(x^{\prime})dx^{\prime}=0,\quad
x\in\Gamma.$

Since here $h\in H^{\frac{3}{2}}(\partial\Omega)$ is arbitrary chosen, we see
from the above identity that $\chi\Theta_{k_{0}}$ must be vanished on
$\Gamma\times\partial\Omega$, which is impossible. This means that
$\lambda_{k}=\tilde{\lambda}_{k}$ for any $k\in\mathbb{N}$. Repeating the
argument used in the proof of [26, (2.8)] implies that
$\Theta_{k}=\tilde{\Theta}_{k}$ on $\Gamma\times\Gamma_{\rm in}^{\prime}$.
This finishes the proof of the lemma. ∎

On the basis of the results in Lemma (3.6), we next show a lemma related to
the algebraic multiplicity of the eigenvalue $\lambda_{k}$ and the
corresponding eigenfunctions.

######  Corollary 3.7.

Under the same assumption in Lemma (3.6), then for any $k\in\mathbb{N}$, there
exist orthonormal sets of eigenfunctions
$\\{\phi_{k,\ell}\\}_{\ell=1}^{m_{k}}$ and
$\\{\tilde{\phi}_{k,\ell}\\}_{\ell=1}^{\tilde{m}_{k}}$ related to
$\lambda_{k}$ of the elliptic operators $A_{e}$ and $\tilde{A}_{e}$
respectively such that

$m_{k}=\tilde{m}_{k},\quad
D\frac{\partial\phi_{k,\ell}}{\partial\nu}=\tilde{D}\frac{\partial\tilde{\phi}_{k,\ell}}{\partial\nu}\mbox{
on }\partial\Omega,\quad k\in\mathbb{N},\ \ell=1,\ldots,m_{k}.$ (3.11)

Before giving a proof of the lemma, we introduce an algebraic result in
[Canuto01] which takes the following form in our context.

######  Lemma 3.8.

Let $M_{1},M_{2}\in\mathbb{N}$ and consider $f_{\ell}\in C(\partial\Omega),$
$\ell=1,\ldots,M_{1}$ and $g_{\ell}\in C(\partial\Omega),$
$\ell=1,\ldots,M_{2}$. Assume that

$\sum_{\ell=1}^{M_{1}}f_{\ell}(x)f_{\ell}(x^{\prime})=\sum_{\ell=1}^{M_{2}}g_{\ell}(x)g_{\ell}(x^{\prime}),\quad
x\in\Gamma,\ x^{\prime}\in\Gamma^{\prime},$ (3.12)

and that the restrictions
$\\{f_{\ell}|_{\Gamma\cap\Gamma^{\prime}}\\}_{\ell=1}^{M_{1}}$ and
$\\{g_{\ell}|_{\Gamma\cap\Gamma^{\prime}}\\}_{\ell=1}^{M_{2}}$ are linearly
independent respectively, where $\Gamma$ and $\Gamma^{\prime}$ are open
subboundaries of $\partial\Omega$ and such that
$\Gamma\cup\Gamma^{\prime}=\partial\Omega$ and
$\Gamma\cap\Gamma^{\prime}\neq\emptyset$. Then there holds $M_{1}=M_{2}$ and
there exists an $M_{1}\times M_{1}$ orthogonal matrix $O$ such that

$(f_{1},\ldots,f_{M_{1}})^{\mathrm{T}}(x)=O(g_{1},\ldots,g_{M_{1}})^{\mathrm{T}}(x),\quad
x\in\partial\Omega,$

where $(\,\cdot\,)^{\mathrm{T}}$ stands for the transpose.

Now on the basis of the above lemma, by an argument similarly to [26], we can
give the proof of Lemma (3.7).

###### Proof of Lemma (3.7).

For any $k\in\mathbb{N}$, according to the settings in Lemma (3.8), we denote
$M_{1}:=m_{k}$, $M_{2}:=\tilde{m}_{k}$ and define the functions $f_{\ell}$ and
$g_{\ell}$ as follows

$\begin{aligned}
f_{\ell}(x)=&D(x)\frac{\partial\varphi_{k,\ell}}{\partial\nu}(x),\
\ell=1,\ldots,M_{1},\\\
g_{\ell}(x)=&\tilde{D}(x)\frac{\partial\tilde{\varphi}_{k,\ell}}{\partial\nu}(x),\
\ell=1,\ldots,M_{2},\end{aligned}\quad x\in\partial\Omega.$

It is not difficult to check that the above settings fulfill the assumptions
in Lemma (3.8). Indeed, by noting that $\varphi_{k,\ell}$ and
$\tilde{\varphi}_{k,\ell}$ are the eigenfunctions to the elliptic operators
$A_{e}$ and $\tilde{A}_{e}$ respectively, we deduce that $\varphi_{k,\ell}$
and $\tilde{\varphi}_{k,\ell}$ are actually in the function space
$C^{1}(\overline{\Omega})$ according to the regularity theories of elliptic
equations, which further gives that the functions $f_{\ell}$,
$\ell=1,\ldots,M_{1}$ and $g_{\ell}$, $\ell=1,\ldots,M_{2}$ lie in
$C(\partial\Omega)$. Moreover, in view of the identities (3.10) in Lemma
(3.6), the condition (3.12) is fulfilled. On the other hand, on the basis of
the unique continuation principle of the solutions for elliptic equations, we
see that $\\{f_{\ell}\\}_{\ell=1}^{M_{1}}$ and
$\\{g_{\ell}\\}_{\ell=1}^{M_{2}}$ restricted to the subboundary $\Gamma_{\rm
in}^{\prime}\cap\Gamma$ are linearly independent respectively. Therefore,
Lemma (3.8) implies that $m_{k}=\tilde{m}_{k}$ and there exists an
$m_{k}\times m_{k}$ orthogonal matrix $O_{k}$ such that

$\left(D\frac{\partial\varphi_{k,1}}{\partial\nu},\ldots,D\frac{\partial\varphi_{k,m_{k}}}{\partial\nu}\right)^{\mathrm{T}}=O_{k}\left(\tilde{D}\frac{\partial\tilde{\varphi}_{k,1}}{\partial\nu},\ldots,\tilde{D}\frac{\partial\tilde{\varphi}_{k,\tilde{m}_{k}}}{\partial\nu}\right)^{\mathrm{T}}\mbox{
on }\partial\Omega.$

Finally, setting $\phi_{k,\ell}:=\varphi_{k,\ell}$, $\ell=1,\ldots,m_{k}$ and

$\left(\tilde{\phi}_{k,1},\ldots,\tilde{\phi}_{k,m_{k}}\right)^{\mathrm{T}}:=O_{k}\left(\tilde{\varphi}_{k,1},\ldots,\tilde{\varphi}_{k,m_{k}}\right)^{\mathrm{T}},$

we see that $\\{\phi_{k,\ell}\\}_{\ell=1}^{m_{k}}$ and
$\\{\tilde{\phi}_{k,\ell}\\}_{\ell=1}^{m_{k}}$ are the required sequences in
Lemma (3.7). This completes the proof of the lemma. ∎

We will complete the proof of the first part of Theorem (1), that is, the
uniqueness determination of the mixture absorption $\mu_{a}+\mu_{f}$. To this
end, we refer to the following inverse spectral result in [10, Corollaries
1.5–1.7].

######  Lemma 3.9.

Under the conditions of Theorem (1), assume that $\nabla D=\nabla\tilde{D}$ on
boundary $\partial\Omega$ holds. Then (3.10) and (3.11) imply
$(D,\mu_{a}+\mu_{f})=(\tilde{D},\tilde{\mu}_{a}+\tilde{\mu}_{f})$.

Applying this lemma, we can show the uniqueness of determining the mixture
absorption $\mu_{a}+\mu_{f}$.

###  3.2 Recovering the background absorption and fluorescence absorption.

In the above subsection, we can uniquely determine the diffusion coefficient
$D$, the mixture absorption $\mu_{a,f}:=\mu_{a}+\mu_{f}$ and $u_{e}$ by the
boundary measurement $D\frac{\partial
u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$. In this part, by substituting
$\mu_{f}=\mu_{a,f}-\mu_{a}$, we rewrite the problem (1.2) into the following
form:

$\begin{cases}\begin{aligned}
\partial_{t}u_{m}+A_{0}u_{m}&+\mu_{a}(x)(u_{m}+u_{e})=\mu_{a,f}u_{e},&&(x,t)\in\Omega\times(0,T),\\\
u_{m}&=0,&&x\in\Omega,\,t=0,\\\
u_{m}&=0,&&(x,t)\in\partial\Omega\times(0,T).\end{aligned}\end{cases}$ (3.13)

Here the operator $A_{0}u:=-\text{div}(D\nabla u)$, $u\in H^{2}(\Omega)$. By
letting $U_{e,m}:=u_{e}+u_{m}$, and noting that $u_{e}$ solves the system
(1.1), that is $\partial_{t}u_{e}+A_{0}u_{e}+\mu_{a,f}u_{e}=0$, we can see
that $U_{e,m}$ reads the system

$\begin{cases}\begin{aligned}
\partial_{t}U_{e,m}+A_{0}U_{e,m}+\mu_{a}(x)U_{e,m}&=0,&&(x,t)\in\Omega\times(0,T),\\\
U_{e,m}&=0,&&x\in\Omega,\,t=0,\\\
U_{e,m}&=g,&&(x,t)\in\partial\Omega\times(0,T).\end{aligned}\end{cases}$
(3.14)

Now the second observation $D\frac{\partial u_{m}}{\partial\nu}$ on
$\Gamma\times(0,T)$ becomes

$D\frac{\partial U_{e,m}}{\partial\nu}=D\left(\frac{\partial
u_{e}}{\partial\nu}+\frac{\partial u_{m}}{\partial\nu}\right)\mbox{ on
}\Gamma\times(0,T).$

We can follow an argument similar to the recovery of the mixture absorption in
the above subsection to show that the boundary data $D\frac{\partial
U_{e,m}}{\partial\nu}|_{\Gamma\times(0,T)}$ can uniquely determine the
background absorption $\mu_{a}$ in the system (3.14). Since the mixture
absorption $\mu_{a,f}=\mu_{a}+\mu_{f}$ is already determined in the above
subsection (3.1), $\mu_{f}$ can be calculated by $\mu_{f}=\mu_{a,f}-\mu_{a}$.

##  4 Numerical inversions.

We will consider the numerical inversions in this section. Since solving the
simultaneous inverse problems (1.6)-(1.7) is computationally demanding, we
have to make a compromise for the DOT-FDOT problem by selecting the parameters
we desire to reconstruct. This work focuses on the determination of both
background absorption coefficient $\mu_{a}$ and fluorescence absorption
coefficient $\mu_{f}$.

###  4.1 The Fréchet derivatives of observation operators.

Recall that $\mu_{a,f}:=\mu_{a}+\mu_{f}$. We consider the following initial
boundary value problems

$\begin{cases}\begin{aligned}
\left(\partial_{t}-\Delta+\mu_{a,f}(x)\right)u_{e}&=0,&&(x,t)\in\Omega\times(0,T),\\\
u_{e}&=0,&&x\in\Omega,\ t=0,\\\
u_{e}&=g,&&(x,t)\in\partial\Omega\times(0,T),\end{aligned}\end{cases}$ (4.1)

and

$\begin{cases}\begin{aligned}
\left(\partial_{t}-\Delta+\mu_{a}(x)\right)u_{m}&=\mu_{f}(x)u_{e},&&(x,t)\in\Omega\times(0,T),\\\
u_{m}&=0,&&x\in\Omega,\ t=0,\\\
u_{m}&=0,&&(x,t)\in\partial\Omega\times(0,T).\end{aligned}\end{cases}$ (4.2)

Define $U_{e,m}:=u_{e}+u_{m}$. Then it satisfies

$\begin{cases}\begin{aligned}
\left(\partial_{t}-\Delta+\mu_{a}(x)\right)U_{e,m}&=0,&&(x,t)\in\Omega\times(0,T),\\\
U_{e,m}&=0,&&x\in\Omega,\ t=0,\\\
U_{e,m}&=g,&&(x,t)\in\partial\Omega\times(0,T).\end{aligned}\end{cases}$ (4.3)

Introduce the admissible set of unknowns by

$\mathcal{A}:=\big{\\{}q\in L^{2}(\Omega):0\leq q\leq C_{+}\ \text{a.e.\ in}\
\Omega\big{\\}},$ (4.4)

where $C_{+}>0$ is a constant. The coming inverse problem is to recover the
$\mu_{a}\in\mathcal{A}$ as well as $\mu_{f}\in\mathcal{A}$ in (4.1)-(4.2) from
the time-dependent sub-boundary data

$\frac{\partial
u_{e}}{\partial\nu}\big{|}_{\Gamma\times(0,T)},\quad\frac{\partial
u_{m}}{\partial\nu}\big{|}_{\Gamma\times(0,T)}.$ (4.5)

The above simultaneous inverse problem can be divided into following two
steps:

  * •

Step 1: determine the $\mu_{a,f}$, i.e., $\mu_{a}+\mu_{f}$, from the data
$\frac{\partial u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$;

  * •

Step 2: separate the $\mu_{a}$ and $\mu_{f}$ from the recovery of Step 1 by
$\big{(}\frac{\partial u_{e}}{\partial\nu}+\frac{\partial
u_{m}}{\partial\nu}\big{)}\big{|}_{\Gamma\times(0,T)}$.

More precisely, we need to solve the following two operator equations

$\mathbb{G}_{e}(\mu_{a,f})=\frac{\partial
u_{e}}{\partial\nu}\big{|}_{\Gamma\times(0,T)}=:\psi_{e},\quad(x,t)\in\Gamma\times(0,T)$
(4.6)

and

$\mathbb{G}_{e,m}(\mu_{a})=\big{(}\frac{\partial
u_{e}}{\partial\nu}+\frac{\partial
u_{m}}{\partial\nu}\big{)}\big{|}_{\Gamma\times(0,T)}=:\psi_{e,m},\quad(x,t)\in\Gamma\times(0,T),$
(4.7)

where $\mathbb{G}_{e}:\mu_{a,f}\to\psi_{e}$ and
$\mathbb{G}_{e,m}:\mu_{a}\to\psi_{e,m}$ are the observation operators, and
defined by (4.1) and (4.3) respectively.

To solve the nonlinear equations (4.6) and (4.7) numerically, the gradient-
type iteration and Newton-type iteration are common choices. To use such
iterations, the differentiability of the operator
${\mathbb{G}}_{e},{\mathbb{G}}_{e,m}$ and the computations of its
Fr$\acute{\text{e}}$chet derivatives are the crucial points. The
differentiability of ${\mathbb{G}}_{e},{\mathbb{G}}_{e,m}$ for
$\mu_{a},\mu_{f}\in\mathcal{A}$ can be proven in the same way in [30]. Here we
only show the expressions of the Fréchet derivatives of
${\mathbb{G}}_{e},{\mathbb{G}}_{e,m}$ in the following. Given the
perturbations of $\mu_{a,f}$ and $\mu_{a}$ as $\tilde{\mu}_{a,f}$ and
$\tilde{\mu}_{a}$ respectively, we can compute the Fréchet derivatives of
${\mathbb{G}}_{e},{\mathbb{G}}_{e,m}$ along the corresponding perturbations as

$\begin{cases}\begin{aligned}
\mathbb{G}^{\prime}_{e}(\mu_{a,f})\tilde{\mu}_{a,f}&:=\lim_{\epsilon\to
0}\frac{\mathbb{G}_{e}(\mu_{a,f}+\epsilon\tilde{\mu}_{a,f})-\mathbb{G}_{e}(\mu_{a,f})}{\epsilon}=\frac{\partial
w_{e}}{\partial\nu}\big{|}_{\Gamma\times(0,T)},\\\
\mathbb{G}^{\prime}_{e,m}(\mu_{a})\tilde{\mu}_{a}&:=\lim_{\epsilon\to
0}\frac{\mathbb{G}_{e,m}(\mu_{a}+\epsilon\tilde{\mu}_{a})-\mathbb{G}_{e,m}(\mu_{a})}{\epsilon}=\frac{\partial
w_{e,m}}{\partial\nu}\big{|}_{\Gamma\times(0,T)},\end{aligned}\end{cases}$
(4.8)

where $w_{e}$ and $w_{e,m}$ satisfy the initial–boundary value problems

$\begin{cases}\begin{aligned}
\left(\partial_{t}-\Delta+\mu_{a,f}(x)\right)w_{e}&=-\tilde{\mu}_{a,f}(x)u_{e}[\mu_{a,f}],&&(x,t)\in\Omega\times(0,T),\\\
w_{e}&=0,&&x\in\Omega,\ t=0,\\\
w_{e}&=0,&&(x,t)\in\partial\Omega\times(0,T),\end{aligned}\end{cases}$ (4.9)

and

$\begin{cases}\begin{aligned}
\left(\partial_{t}-\Delta+\mu_{a}(x)\right)w_{e,m}&=-\tilde{\mu}_{a}(x)U_{e,m}[\mu_{a}],&&(x,t)\in\Omega\times(0,T),\\\
w_{e,m}&=0,&&x\in\Omega,\ t=0,\\\
w_{e,m}&=0,&&(x,t)\in\partial\Omega\times(0,T).\end{aligned}\end{cases}$
(4.10)

For the exact data $\psi_{e}$ and $\psi_{m}$, the equations (4.6) and (4.7)
admit unique solutions from Theorem (1). However, in practice, we can only
have the noisy data $\psi_{e}^{\delta_{e}}$ and $\psi_{e,m}^{\delta_{e,m}}$
satisfying
${\|\psi_{e}^{\delta_{e}}-\psi_{e}\|_{2}}/{\|\psi_{e}\|_{2}}\leq\delta_{e}$
and
${\|\psi_{e,m}^{\delta_{e,m}}-\psi_{m}\|_{2}}/{\|\psi_{m}\|_{2}}\leq\delta_{e,m}$
respectively, where $\delta_{e},\,\delta_{e,m}\geq 0$ are the noise levels.
Then, solving (4.6) and (4.7) are equivalent to solve the operator equation of
the form

$\mathbb{G}(q)=\psi^{\delta},$ (4.11)

where
$\mathbb{G}=\mathbb{G}_{e},q=\mu_{a,f},\psi^{\delta}=\psi_{e}^{\delta_{e}}$
for (4.6) and
$\mathbb{G}=\mathbb{G}_{e,m},q=\mu_{a},\psi^{\delta}=\psi_{e,m}^{\delta_{e,m}}$
for (4.7). Due to the ill-posedness of these inverse problems, we will use the
regularized inversion algorithms.

The gradient-type iteration such as Landweber iteration [17] and the Newton-
type iteration such as Levenberg–Maquardt method [16], the iteratively
regularized Gauss-Newton method [23, 24] are common regularization methods to
solve (4.11). For instance, the nonlinear Landweber iteration in Hilbert space
is defined via

$q_{n+1}^{\delta}=q_{n}^{\delta}+\big{(}\mathbb{G}^{\prime}(q_{n}^{\delta})\big{)}^{*}(\psi^{\delta}-\mathbb{G}(q_{n}^{\delta})),$
(4.12)

where $(\mathbb{G}^{\prime})^{*}$ means the conjugate operator of
$\mathbb{G}^{\prime}$, the Levenberg–Maquardt method is defined via

$q_{n+1}^{\delta}=q_{n}^{\delta}+\big{[}(\mathbb{G}^{\prime}(q_{n}^{\delta}))^{*}\mathbb{G}^{\prime}(q_{n}^{\delta})+\alpha_{n}I\big{]}^{-1}\big{(}\mathbb{G}^{\prime}(q_{n}^{\delta})\big{)}^{*}(\psi^{\delta}-\mathbb{G}(q_{n}^{\delta})),$
(4.13)

and the iteratively regularized Gauss–Newton method is defined via

$q_{n+1}^{\delta}=q_{n}^{\delta}+\big{[}(\mathbb{G}^{\prime}(q_{n}^{\delta}))^{*}\mathbb{G}^{\prime}(q_{n}^{\delta})+\alpha_{n}I\big{]}^{-1}\big{(}\big{(}\mathbb{G}^{\prime}(q_{n}^{\delta})\big{)}^{*}(\psi^{\delta}-\mathbb{G}(q_{n}^{\delta}))+\alpha_{n}(q_{0}-q_{n}^{\delta})\big{)}.$
(4.14)

The Landweber iteration is quite slow. As an alternative to Landweber-type
methods, the Levenberg–Maquardt method or the iteratively regularized
Gauss–Newton method requires less number of iterations to satisfy the stopping
rule than the Landweber iteration, however each update step of those
iterations might take considerably longer than one step of Landweber
iteration. In fact, the operator
$[(\mathbb{G}^{\prime}(q_{n}^{\delta}))^{*}\mathbb{G}^{\prime}(q_{n}^{\delta})+\alpha_{n}I]^{-1}$
in (4.13) or (4.14) has to be solved in practical applications, implying a
huge linear system of equations has to be solved, which often proves to be
costly, if not infeasible. Therefore, it may be more desirable to accelerate
Landweber iteration by preserving its simple implementation feature. We will
introduce an accelerated Landweber iteration scheme in the next subsection.

###  4.2 Nesterov acceleration of Landweber iteration.

Suppose that $\Omega$ contains a finite number of inclusions (fluorophores)
$\\{E_{j}\\}_{j=1}^{m}$, and denote the collection of inclusions by
$E=\cup_{j=1}^{m}E_{j}$, and assume that there exist constants $d_{1},d_{2}>0$
such that

${\rm dist}(E,\partial\Omega)>d_{1},\;{\rm dist}(E_{i},E_{j})>d_{2}\
\text{for}\ 1\leq i\neq j\leq m.$

The unknown $\mu_{f}$ is a non-smooth function supported on $E$. Furthermore,
suppose each inclusion $E_{j}$ has the fluorescence intensity $\mu_{j}$. Then
we have the fluorescence absorption of the form

$\displaystyle\mu_{f}=\begin{cases}\mu_{j},&x\in E_{j},\;j=1,2,\cdots,m,\\\
0,&x\in\Omega\setminus\overline{E}.\end{cases}$

The boundary of $E_{j}$ denotes the shape of $j$-th fluorophore (object)
inside the medium. Therefore, when to identify $\mu_{a,f}$, $\mu_{a}$ and
$\mu_{f}$ by solving the nonlinear ill-posed equations (4.6) and (4.7), the
interfaces of $\mu_{f}$ should be kept.

However, the classical Landweber iteration in Hilbert spaces and the classical
$L^{2}$ penalty term in Tikhonov regularization have the tendency to over-
smooth solutions, which makes it difficult to capture special features of the
sought solution such as sparsity and discontinuity [31]. To overcome this
drawback, various regularization schemes in Banach spaces or in a manner of
incorporating general non-smooth convex penalty terms have been proposed [7,
22, 25]. Further, a so called two-point gradient (TPG) method was proposed in
[19], which can be considered as some generalization of an accelerated version
of Landweber iteration based on Nesterov’s strategy. Some convergence results
of this method were also concerned in [19, 20, 31]. Due to the usefulness and
acceleration of TPG method in dealing with nonlinear ill-posed problems, we
will apply this method to solve our inverse problems.

Now we give a sketch for solving (4.11) by TPG method for the completeness of
this article, which essentially consists of considering

$\min_{q\in\mathcal{A}}\mathcal{J}(q),\quad\min_{q\in\mathcal{A}}\mathcal{R}(q)$

alternatively in the iteration process, where the date-match functional
$\mathcal{J}(q)$ is given by

$\mathcal{J}(q):=\frac{1}{2}\|\mathbb{G}(q)-\psi^{\delta}\|_{2}^{2},$ (4.15)

and $\mathcal{R}$ is the $L^{2}$ combined with TV penalty term given by

$\mathcal{R}(q)=\beta_{1}\|q\|_{2}^{2}+\beta_{2}|q|_{\rm TV}.$ (4.16)

Pick $q_{-1}=q_{0}\in\mathcal{A}$ and
$\xi_{-1}=\xi_{0}\in\partial\mathcal{R}(q_{0})$ as the initial guesses, where
$\partial\mathcal{R}(q):=\big{\\{}\xi\in
L^{2}(\Omega):\mathcal{R}(\tilde{q})-\mathcal{R}(q)-\langle\xi,\tilde{q}-q\rangle\geq
0\ \text{for all}\ \tilde{q}\in L^{2}(\Omega)\big{\\}}$. For $n\geq 0$, we
define

$\begin{cases}\begin{aligned}
\eta_{n}&=\xi_{n}+\lambda_{n}(\xi_{n}-\xi_{n-1}),\\\ z_{n}&=\arg\min_{z\in
L^{2}(\Omega)}\big{\\{}\mathcal{R}(z)-\langle\eta_{n},z\rangle\big{\\}},\\\
\xi_{n+1}&=\xi_{n}+\alpha_{n}\big{(}\mathbb{G}^{\prime}(z_{n})\big{)}^{*}(\psi^{\delta}-\mathbb{G}(z_{n})),\\\
q_{n+1}&=\arg\min_{q\in
L^{2}(\Omega)}\big{\\{}\mathcal{R}(q)-\langle\xi_{n+1},q\rangle\big{\\}},\\\
\end{aligned}\end{cases}$ (4.17)

where $\lambda_{n}\geq 0$ is the combination parameter, $\alpha_{n}$ is the
step sizes defined by

$\alpha_{n}=\begin{cases}\min\left\\{\frac{\bar{\alpha}_{0}\|\mathbb{G}(z_{n})-\psi^{\delta}\|_{2}^{2}}{\|(\mathbb{G}^{\prime}(z_{n}))^{*}(\mathbb{G}(z_{n})-\psi^{\delta})\|_{2}^{2}},\bar{\alpha}_{1}\right\\}&\text{if}\
\|\mathbb{G}(z_{n})-\psi^{\delta}\|_{2}>\tau\delta,\\\ 0&\text{if}\
\|\mathbb{G}(z_{n})-\psi^{\delta}\|_{2}\leq\tau\delta,\end{cases}$

for some positive constants $\bar{\alpha}_{0}$ and $\bar{\alpha}_{1}$. The
combination parameter $\lambda_{n}\geq 0$ satisfies $\lambda_{0}=0,\
\lambda_{n}\in[0,1]$ for all $n\in\mathbb{N}$. Further, the iteration (4.17)
will be stopped by the discrepancy principle with respect to $z_{n}^{\delta}$
i.e., for a given $\tau>1$, we will terminate the iteration after $n_{*}$
steps, where $n_{*}:=n(\delta,\psi^{\delta})$ is the integer such that

$\|\mathbb{G}(z_{n_{*}})-\psi^{\delta}\|_{2}\leq\tau\delta<\|\mathbb{G}(z_{n_{*}-1})-\psi^{\delta}\|_{2},\quad
0\leq n<n_{*}.$ (4.18)

Unless otherwise specified, we let the combination parameter $\lambda_{n}$ be
chosen by $\lambda_{n}=\frac{n}{n+10}$ in this paper, which leads to a refined
version of the Nesterov acceleration of Landweber iteration [31].

######  Remark 4.1.

Let $q^{\dagger}\in\mathcal{A}$ be the unique solution of (4.11) with
$\delta=0$. Let the initial guesses $q_{0}\in\mathcal{A}$ and
$\xi_{0}\in\partial\mathcal{R}(q_{0})$ satisfy
$\mathcal{R}(q^{\dagger})-\mathcal{R}(q_{0})-\langle\xi_{0},q^{\dagger}-q_{0}\rangle\leq
c_{0}\rho^{2}$ for some $\rho>0$. Moreover, let
$n_{*}=n_{*}(\delta,\psi^{\delta})$ be chosen according to the stopping rule
(4.18) with $\tau>1$. If the combination parameter $\lambda_{n}^{\delta}$ in
(4.17) is chosen by the discrete backtracking search algorithm [31] and the
observation operators $\mathbb{G}$ satisfy the tangential condition

$\|\mathbb{G}(q)-\mathbb{G}(\tilde{q})-\mathbb{G}^{\prime}(\tilde{q})(q-\tilde{q})\|_{2}\leq\eta\|\mathbb{G}(q)-\mathbb{G}(\tilde{q})\|_{2},\
\ \text{for all}\ \ q,\tilde{q}\in\mathcal{A},$ (4.19)

there holds $\|q_{n_{*}}^{\delta}-q^{\dagger}\|_{2}\to 0$ as $\delta\to 0$,
where $q_{n_{*}}^{\delta}$ is the result of (4.17) at the $n_{*}$-th
iteration.

###  4.3 Numerical examples.

We will provide several examples for the implementation of the iterative
algorithm. To run the iteration scheme (4.17) for solving our inverse
problems, we need to derive
$\big{(}\mathbb{G}_{e}^{\prime}(\mu_{a,f})\big{)}^{*}$ and
$\big{(}\mathbb{G}_{e,m}^{\prime}(\mu_{a})\big{)}^{*}$ respectively. By the
definitions of $w_{e}$ and $w_{e,m}$ in (4.8), we have

$\begin{cases}\begin{aligned}
\big{\langle}\mathbb{G}_{e}^{\prime}(\mu_{a,f})\tilde{\mu}_{a,f},\gamma\big{\rangle}_{L^{2}(\Gamma\times(0,T))}&=\int_{\Gamma\times(0,T)}\frac{\partial
w_{e}[\tilde{\mu}_{a,f}]}{\partial\nu}(x,t)\gamma(x,t)\,dxdt,\\\
\big{\langle}\mathbb{G}_{e,m}^{\prime}(\mu_{a})\tilde{\mu}_{a},\gamma\big{\rangle}_{L^{2}(\Gamma\times(0,T))}&=\int_{\Gamma\times(0,T)}\frac{\partial
w_{e,m}[\tilde{\mu}_{a}]}{\partial\nu}(x,t)\gamma(x,t)\,dxdt,\end{aligned}\end{cases}$
(4.20)

where $w_{e}$ and $w_{e,m}$ are the solutions to (4.9)-(4.10) respectively.
Then we can have the representations of
$(\mathbb{G}_{e}^{\prime}(\mu_{a,f}))^{*}$ and
$(\mathbb{G}_{e,m}^{\prime}(\mu_{f}))^{*}$ by the identities

$\begin{cases}\begin{aligned}
\big{\langle}\mathbb{G}_{e}^{\prime}(\mu_{a,f})\tilde{\mu}_{a,f},\gamma\big{\rangle}_{L^{2}(\Gamma\times(0,T))}&\equiv\big{\langle}\tilde{\mu}_{a,f},(\mathbb{G}_{e}^{\prime}(\mu_{a,f}))^{*}\gamma\big{\rangle}_{L^{2}(\Omega)},\\\
\big{\langle}\mathbb{G}_{e,m}^{\prime}(\mu_{a})\tilde{\mu}_{a},\gamma\big{\rangle}_{L^{2}(\Gamma\times(0,T))}&\equiv\big{\langle}\tilde{\mu}_{a},(\mathbb{G}_{e,m}^{\prime}(\mu_{f}))^{*}\gamma\big{\rangle}_{L^{2}(\Omega)},\end{aligned}\end{cases}$
(4.21)

for all $\mu_{a},\mu_{f}\in\mathcal{A}$ and $\gamma\in
L^{2}(\Gamma\times(0,T))$. We state the results as the following lemma.

######  Lemma 4.2.

There hold the variational identities given by

$(\mathbb{G}_{e}^{\prime}(\mu_{a,f}))^{*}\gamma=-\int_{0}^{T}u_{e}[\mu_{a,f}]\phi_{e}[\gamma]\,dt,$
(4.22)

and

$(\mathbb{G}_{e,m}^{\prime}(\mu_{a}))^{*}\gamma=-\int_{0}^{T}\big{(}u_{m}[\mu_{a}]+u_{e}[\mu_{a,f}]\big{)}\phi_{e,m}[\gamma]\,dt,$
(4.23)

where $\phi_{e}[\gamma]$ and $\phi_{e,m}[\gamma]$ are the solutions to the
initial boundary value problems

$\begin{cases}\begin{aligned}
\left(-\partial_{t}-\Delta+\mu_{a,f}(x)\right)\phi_{e}&=0,\qquad\qquad\qquad
x\in\Omega,\,t\in(0,T),\\\ \phi_{e}&=0,\qquad\qquad\qquad x\in\Omega,\,t=T,\\\
\phi_{e}&=\begin{cases}\gamma(x,t),&\qquad\;x\in\Gamma,\,t\in(0,T),\\\
0,&\qquad\;x\in\partial\Omega\setminus\Gamma,\,t\in(0,T),\end{cases}\end{aligned}\end{cases}$
(4.24)

and

$\begin{cases}\begin{aligned}
\left(-\partial_{t}-\Delta+\mu_{a}(x)\right)\phi_{e,m}&=0,\qquad\qquad\qquad
x\in\Omega,\,t\in(0,T),\\\ \phi_{e,m}&=0,\qquad\qquad\qquad
x\in\Omega,\,t=T,\\\
\phi_{e,m}&=\begin{cases}\gamma(x,t),&\qquad\;x\in\Gamma,\,t\in(0,T),\\\
0,&\qquad\;x\in\partial\Omega\setminus\Gamma,\,t\in(0,T).\end{cases}\end{aligned}\end{cases}$
(4.25)

###### Proof.

Multiplying two sides of the equation in (4.9) by
$\phi_{e}(x,t):=\phi_{e}[\gamma](x,t)$ and integrating in $\Omega\times(0,T)$,
we get

$\displaystyle\int_{\Omega\times(0,T)}\left(\partial_{t}-\Delta+\mu_{a,f}\right)w_{e}\phi_{e}\,dxdt=-\int_{\Omega\times(0,T)}\tilde{\mu}_{a,f}u_{e}[\mu_{a,f}]\phi_{e}\,dxdt.$
(4.26)

By $\phi_{e}(x,T)=0$, $w_{e}(x,0)=0$ and $\phi_{e},w_{e}\in
L^{2}(0,T;H^{2}(\Omega))$, we have

$\displaystyle\int_{\Omega\times(0,T)}\ \partial_{t}w_{e}\,\phi_{e}\,dxdt$
$\displaystyle={\int_{\Omega}(w_{e}\phi_{e})|_{t=0}^{t=T}\;dx}-\int_{\Omega\times(0,T)}w_{e}{\partial_{t}}\phi_{e}\,dxdt$
(4.27)
$\displaystyle=-\int_{\Omega\times(0,T)}w_{e}{\partial_{t}}\phi_{e}\,dxdt.$

We also have from the Green’s formula that

$\displaystyle\int_{\Omega\times(0,T)}(-\Delta+\mu_{a,f})w_{e}\phi_{e}\,dxdt$
$\displaystyle=\int_{\Omega\times(0,T)}(-\Delta+\mu_{a,f})\phi_{e}w_{e}\,dxdt$
(4.28)
$\displaystyle\qquad+\int_{0}^{T}\int_{\partial\Omega}\left(\phi_{e}\frac{\partial
w_{e}}{\partial\nu}-w_{e}\frac{\partial\phi_{e}}{\partial\nu}\right)\,dxdt.$

Combining (4.26)-(4.28) together and using (4.24), we finally have

$\displaystyle-\int_{\Omega\times(0,T)}\tilde{\mu}_{a,f}\,u_{e}[\mu_{a,f}]\,\phi_{e}[\gamma]\,dxdt$
$\displaystyle=\int_{0}^{T}\int_{\partial\Omega}\left(\phi_{e}\frac{\partial
w_{e}}{\partial\nu}-w_{e}\frac{\partial\phi_{e}}{\partial\nu}\right)\,dxdt$
(4.29) $\displaystyle=\int_{0}^{T}\int_{\Gamma}\frac{\partial
w_{e}[\tilde{\mu}_{a,f}]}{\partial\nu}\gamma dxdt.$

Then, (4.20) together with (4.21) lead to
$(\mathbb{G}_{e}^{\prime}(\mu_{a,f}))^{*}\gamma=\int_{0}^{T}u_{e}[\mu_{a,f}]\phi_{e}[\gamma]\,dxdt$.
In the same way, we can derive the expression of
$(\mathbb{G}_{m}^{\prime}(\mu_{a}))^{*}\gamma$ which is given as in (4.23).
The proof is complete. ∎

For the noisy data $\psi_{e}^{\delta_{e}}$ and $\psi_{e,m}^{\delta_{e,m}}$, we
describe them as

$\psi_{e}^{\delta_{e}}=\psi_{e}(1+\delta_{e}\zeta)$

and

$\psi_{e,m}^{\delta_{e,m}}=\psi_{e,m}(1+\delta_{e,m}\zeta),$ (4.30)

where $\delta_{e},\,\delta_{e,m}\geq 0$ are the noise levels and $\zeta$ is a
random standard Gaussian noise. Noting the randomness of the measured noise,
we use the mean value of the realizations of the measurements (4.30).
Considering the difficulty of obtaining the medical data, the amount of
realizations can not be too large and we set it be 10 in our experiments. For
the numerical examples, we solve the forward problems (4.1)-(4.3) and dual
problems (4.24)-(4.25) by the finite difference method and give an uniform
mesh with the time mesh size $0.01$. The primal dual hybrid gradient method
will be applied to solve (4.17) iteratively. We set $\tau=1.05$ in (4.18) to
stop the iteration. In the sequel, denote by
$\mu_{a}^{\dagger},\,\mu_{f}^{\dagger},\,\mu_{a,f}^{\dagger}$ the exact
solutions and $\mu_{a}^{N},\,\mu_{f}^{N},\,\mu_{a,f}^{N}$ the recoveries,
where $N$ is the final iteration number.

Now we are ready to show the performances of above iterative algorithm for
recovering the unknown absorption coefficients. Set $T=1$,
$\Omega:=\\{x=(x_{1},x_{2})\in\mathbb{R}^{2}|\,0\leq x_{1}\leq 1,\,0\leq
x_{2}\leq 1\\}$ and suppose
$u_{e}(0,x_{2},t)=e^{\frac{x_{2}}{2}-t},\,u_{e}(1,x_{2},t)=e^{\frac{1+x_{2}}{2}-t},\,u_{e}(x_{1},0,t)=e^{\frac{x_{1}}{2}-t}$
and $u_{e}(x_{1},1,t)=e^{\frac{1+x_{1}}{2}-t}$. We consider the numerical
inversions for the following two cases:

  * •

(1a) the homogeneous background absorption and the non-smooth fluoresence
absorption given by

$\displaystyle\mu_{a}^{\dagger}\equiv
1,\quad\mu_{f}^{\dagger}=\begin{cases}2,&x\in E,\\\ 0,&{\rm else},\end{cases}$
(4.31)

where $E$ is a disc with radius $0.2$, centered at $(0.4,0.4)$;

  * •

(1b) the non-smooth background absorption and the non-smooth fluoresence
absorption given by

$\displaystyle\mu_{a}^{\dagger}=\begin{cases}2,&x\in E_{1},\\\ 1,&{\rm
else},\end{cases}\quad\mu_{f}^{\dagger}=\begin{cases}2,&x\in E_{2},\\\ 0,&{\rm
else},\end{cases}$ (4.32)

where $E_{1},\ E_{2}$ are two discs with radius $0.15$, centered at
$(0.3,0.3)$ and $(0.7,0.7)$, respectively.

First, set $\beta_{1}=0.1,\beta_{2}=1$ and pick the initial guess as
$\xi_{0}=0$. The results of Step 1 (i.e., recovering
$\mu_{a,f}:=\mu_{a}+\mu_{f}$) from different noisy level data are listed in
Table (1) and the recovered solutions are plotted in Figure (1) and Figure
(2), respectively.

Next, fix the noise level of boundary measurement $\frac{\partial
u_{e}}{\partial\nu}|_{\Gamma\times(0,T)}$ as $\delta_{e}=0.1\%$ for Case (1a)
and $\delta_{e}=0.01\%$ for Case (1b). Still pick the initial guess as
$\xi_{0}=0$, the results of Step 2 from the measurement data with different
noisy level $\delta_{e,m}$ are listed in Table (2) and the recovered solutions
are plotted in Figure (3) and Figure (4), respectively.

Table 1: Numerical results of Step 1. $\mu_{a,f}^{\dagger}:=\mu_{a}^{\dagger}+\mu_{f}^{\dagger}$ | $\delta_{e}$ | $\|\mathbb{G}_{e}(\mu_{a,f}^{N})-\psi_{e}^{\delta_{e}}\|_{2}$ | $\frac{\|\mu_{a,f}^{N}-\mu_{a,f}^{\dagger}\|_{2}}{\|\mu_{a,f}^{\dagger}\|_{2}}$ | $N$
---|---|---|---|---
Case (1a) | $1\%$ | 1.5485e-5 | 0.1035 | 605
| $0.1\%$ | 1.6742e-7 | 0.0725 | 914
| $0$ | 1.1929e-7 | 0.0697 | 6588
Case (1b) | $0.1\%$ | 1.3534e-6 | 0.1957 | 2373
| $0.01\%$ | 1.9077e-7 | 0.1110 | 4772
| $0$ | 8.4075e-8 | 0.0718 | 14980
Figure 1: The recoveries of $\mu_{a,f}$ in Step 1 for the Case (1a). Figure 2: The recoveries of $\mu_{a,f}$ in Step 1 for the Case (1b). Table 2: Numerical results of Step 2 ($\delta_{e}=0.1\%$). $\mu_{a}^{\dagger},\;\mu_{f}^{\dagger}$ | $\delta_{e}$ | $\delta_{e,m}$ | $\|\mathbb{G}_{e,m}(\mu_{f}^{N})-\psi_{e,m}^{\delta_{e,m}}\|_{2}$ | $\frac{\|\mu_{a}^{N}-\mu_{a}^{\dagger}\|_{2}}{\|\mu_{a}^{\dagger}\|_{2}}$ | $\frac{\|\mu_{f}^{N}-\mu_{f}^{\dagger}\|_{2}}{\|\mu_{f}^{\dagger}\|_{2}}$ | $N$
---|---|---|---|---|---|---
Case (1a) | $0.1\%$ | $5\%$ | 3.6800e-2 | 0.1315 | 0.1931 | 5
|  | $1\%$ | 7.4000e-3 | 0.0304 | 0.1344 | 21
|  | $0.1\%$ | 7.3593e-4 | 0.0094 | 0.1334 | 129
Case (1b) | $0.01\%$ | $1\%$ | 1.2439e-5 | 0.1398 | 0.2690 | 597
|  | $0.1\%$ | 1.2997e-6 | 0.1111 | 0.2740 | 2262
|  | $0.01\%$ | 2.7656e-7 | 0.0794 | 0.2675 | 7950

From Tables (1)-(2) and Figures (1)-(4), we see that the recoveries are
satisfactory in view of the high ill-posedness of our inverse problems, which
is caused by the nonlinearity of this inverse problems and the applied
boundary measurements. The algorithm can recover the overall shape of the
interface of unknown absorptions. These results indicate that we may recover
the piecewise constant $(\mu_{a},\mu_{f})$ simultaneously from the boundary
data excited by general smooth sources, which is not covered by Theorem (1)
and worth our further investigation.

However, numerically, due to the non-smoothness, the interfaces of
$\mu_{a,f}$, $\mu_{a}$ or $\mu_{f}$ are much more challenging to recover and
hence less accurately resolved. Especially, as shown in Figure (2), the
recovery of $\mu_{a,f}$ around the part between disc $E_{1}$ and disc $E_{2}$
has lower accuracy. Furthermore, as shown in Figure (4), the recovery error of
$\mu_{a,f}$ in Step 1 will influence the recovery accuracy of $\mu_{f}$ in
Step 2, since we need capture the interface of $\mu_{f}$ by
$\mu_{a,f}-\mu_{a}$. By comparing the results for Case (1a) and Case (1b) as
shown in Figures (3)-(4), we can see that the non-smooth background absorption
and non-smooth fluoresence absorption is more difficult to be accurately
separated than the case with homogeneous background absorption. Furthermore,
the numerical results indicate that the algorithm does converge, but it may
take many iterations to yield satisfactory recoveries, even if the Nesterov
acceleration is applied in the iteration process, reflecting the severe ill-
posedness of nonlinear multiple coefficients simultaneous inversion.

######  Remark 4.3.

It is well-known that the convergence and convergence speed of iterative
algorithm heavily depend on the initial guess. In above numerical examples, we
fixed the initial guess as $\xi_{0}=0$ and investigated the performance of
proposed algorithm. In fact, one can extract the prior information of unknown
inclusions $\\{E_{j}\\}_{j=1}^{m}$ by applying our shape-constraint strategy
in [8], which can estimate the positions of unknown inclusions by cuboid
approximation. Then we can get a good initial guess for the iterative
algorithm in this work and save the number of iteration.

Figure 3: The recoveries of $\mu_{f}$ and $\mu_{a}$ in Step 2 for the Case
(1a). Figure 4: The recoveries of Step 2 for the Case (1b).

##  5 Conclusion and future works.

This work studies a multiple coefficients inversion problem arising from the
time-domain diffuse optical tomography with fluorescence, which is to recover
the diffusion coefficient, background absorption as well as fluorescence
absorption from one single boundary measurement. We deduce the uniqueness
theorem and do the numerical inversion by the Nesterov accelerated Landweber
iteration. The numerical results are satisfactory in view of the high ill-
posedness of our simultaneous inverse problem. However, to make the theory and
the numerical experiments in this work more applicable, we will focus on the
following two aspects of future works.

Firstly, in this work, we assumed that the absorption and scattering
coefficients for the excitation and emission lights are the same, which leads
to the same $D$ and $\mu_{a}$ in (1.1) and (1.2). However, in many practical
applications, the absorption and scattering coefficients of the background
medium are not same for the lights with different wave lengths. Hence one of
our next works is to prove the uniqueness theorem of recovering
$(D_{e},D_{m},\mu_{a}^{e},\mu_{a}^{m},\mu_{f})$ for the DOT-FDOT with the
governing equations given by

$c^{-1}\partial_{t}u_{e}-\text{div}(D_{e}\nabla
u_{e})+(\mu_{a}^{e}+\mu_{f})u_{e}=0,$

and

$c^{-1}\partial_{t}u_{m}-\text{div}(D_{m}\nabla
u_{e})+\mu_{a}^{m}u_{m}=\mu_{f}u_{e}.$

Secondly, to save the computational costs, we only consider the simultaneous
inversion of background absorption and fluoresence absorption in the numerical
experiments in this work. The fast algorithms of recovering the background
absorption, fluoresence absorption as well as the diffusion coefficient are
important and worth our further investigations.

## Acknowledgment.

Zhiyuan Li thanks National Natural Science Foundation of China (Grant
No.12271277). Chunlong Sun thanks National Natural Science Foundation of China
(Grant No.12201298), Natural Science Foundation of Jiangsu Province, China
(Grant No.BK20210269) and “Double Innovation” Doctor of Jiangsu Province,
China (Grant No.JSSCBS20220227). This work is partially supported by the Open
Research Fund of Key Laboratory of Nonlinear Analysis & Applications (Central
China Normal University), Ministry of Education, P. R. China.

## 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.  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.  doi:10.1088/0266-5611/25/12/123010.
  * [6] Guillaume Bal.  Inverse transport theory and applications.  Inverse Problems, 25(5):053001, 48, 2009.  doi:10.1088/0266-5611/25/5/053001.
  * [7] Radu Ioan Boţ and Torsten Hein.  Iterative regularization with a general penalty term—theory and application to l1 and tv regularization.  Inverse Problems, 28(10):104010, 2012.  doi:10.1088/0266-5611/28/10/104010.
  * [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] Bruno Canuto and Otared Kavian.  Determining coefficients in a class of heat equations via boundary measurements.  SIAM Journal on Mathematical Analysis, 32(5):963–986, 2001.
  * [10] Bruno Canuto and Otared Kavian.  Determining two coefficients in elliptic operators via boundary spectral data: a uniqueness result.  Bollettino Della Unione Matematica Italiana B, 1(1), 2004.
  * [11] 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.  doi:10.1080/03605309608821207.
  * [12] 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.
  * [13] Guillaume Dollé.  Diffuse optical tomography and fluorescence for tumour detection.  PhD thesis, Université de Strasbourg, 2018.
  * [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.  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] Martin Hanke.  A regularizing levenberg-marquardt scheme, with applications to inverse groundwater filtration problems.  Inverse problems, 13(1):79, 1997.
  * [17] Neubauer A. Hanke M. and Scherzer O.  A convergence analysis of the landweber iteration for nonlinear ill-posed problems.  Numerische Mathematik, 72(1):21–37, 1995.  doi:10.1007/s002110050158.
  * [18] Bastian Harrach.  On uniqueness in diffuse optical tomography.  Inverse Problems, 25(5):055010, 14, 2009.  doi:10.1088/0266-5611/25/5/055010.
  * [19] Simon Hubmer and Ronny Ramlau.  Convergence analysis of a two-point gradient method for nonlinear ill-posed problems.  Inverse Problems, 33(9):095004, 2017.  doi:10.1088/1361-6420/aa7ac7.
  * [20] Simon Hubmer and Ronny Ramlau.  Nesterov’s accelerated gradient method for nonlinear ill-posed problems with a locally convex residual functional.  Inverse Problems, 34(9):095003, 2018.  doi:10.1088/1361-6420/aacebe.
  * [21] H. B. Jiang.  Diffuse optical tomography: principles and applications.  CRC Press, Taylor $\&$ Francis Group, Boca Raton, 2011.
  * [22] Qinian Jin and Wei Wang.  Landweber iteration of kaczmarz type with general non-smooth convex penalty functionals.  Inverse Problems, 29(8):085011, 2013.  doi:10.1088/0266-5611/29/8/085011.
  * [23] Zhong M. Jin Q.  On the iteratively regularized gauss-newton method in banach spaces with applications to parameter identification problems.  Numerische Mathematik, 124(2013):647–683, 2013.  doi:10.1007/s00211-013-0529-5.
  * [24] B Kaltenbacher.  A convergence rates result for an iteratively regularized gauss–newton–halley method in banach space.  Inverse Problems, 31(1):015007, 2015.  doi:10.1088/0266-5611/31/1/015007.
  * [25] Barbara Kaltenbacher, Frank Schöpfer, and Thomas Schuster.  Iterative methods for nonlinear ill-posed problems in banach spaces: convergence and applications to parameter identification problems.  Inverse Problems, 25(6):065003, 2009.  doi:10.1088/0266-5611/25/6/065003.
  * [26] Yavar Kian, Zhiyuan Li, Yikan Liu, and Masahiro Yamamoto.  The uniqueness of inverse problems for a fractional equation with a single measurement.  Mathematische Annalen, pages 1–31, 2020.
  * [27] 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.
  * [28] 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.  doi:10.3934/ipi.2020025.
  * [29] 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.
  * [30] M.M. Zhang, J.J. Liu .  Recovery of non-smooth radiative coefficient from nonlocal observation by diffusion system.  J. Inverse Ill-Posed Probl., 28(3):389–410, 2020.
  * [31] Wei Wang Min Zhong and Qinian Jin.  Regularization of inverse problems by two-point gradient methods in banach spaces.  Numerische Mathematik, 143(2019):713–747, 2019.  doi:10.1007/s00211-019-01068-0.
  * [32] V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder.  Fluorescence molecular tomography resolves protease activity in vivo.  Nature Medicine, 2002.
  * [33] B. W. Pogue and M. A. Mycek.  Handbook of Biomedical Fluorescence.  Marcel Dekker, New York, 2003.
  * [34] 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.  doi:10.1088/0266-5611/16/2/312.
  * [35] M. Rudin.  Molecular imaging: basic principles and applications in Biomedical research, 2nd ed. Imperial College Press, London, 2013.
  * [36] Kenichi Sakamoto and Masahiro Yamamoto.  Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems.  J. Math. Anal. Appl., 382(1):426–447, 2011.  doi:10.1016/j.jmaa.2011.04.058.
  * [37] C. Sun and Z. Zhang.  Uniqueness and numerical inversion in the time-domain fluorescence diffuse optical tomography.  Inverse Problems, 38(10):104001, 23, 2022.