# Theoretical and numerical studies of inverse source problem for the linear
parabolic equation with sparse boundary measurements

Guang Lin111Department of Mathematics and Mechanical Engineering, Purdue
University, West Lafayette, IN 47906, USA. (Email: guanglin@purdue.edu),
Zecheng Zhang 222Department of Mathematics, Purdue University, West Lafayette,
IN 47906, USA. (Email: zecheng.zhang.math@gmail.com), Zhidong Zhang 333School
of Mathematics (Zhuhai), Sun Yat-sen University, Zhuhai 519082, Guangdong,
China. (Email<EMAIL_ADDRESS>
###### Abstract

We consider the inverse source problem in the parabolic equation, where the
unknown source possesses the semi-discrete formulation. Theoretically, we
prove that the flux data from any nonempty open subset of the boundary can
uniquely determine the semi-discrete source. This means the observed area can
be extremely small, and that is why we call the data as sparse boundary data.
For the numerical reconstruction, we formulate the problem from the Bayesian
sequential prediction perspective and conduct the numerical examples which
estimate the space-time-dependent source state by state. To better demonstrate
the performance of the method, we solve two common multiscale problems from
two models with a long sequence of the source. The numerical results
illustrate that the inversion is accurate and efficient.

Keywords: inverse source problem, parabolic equation, sparse boundary
measurements, uniqueness, numerical reconstruction.

AMS Subject Classifications: 35R30, 65M32, 93E11.

##  1 Introduction.

###  1.1 Background and literature.

As a classical type of PDEs, the parabolic equation is widely applied in
physics, engineering, finance, and so on. The inverse source problems in the
parabolic equation have various applications in the real world, and the
corresponding research has a long history. We list several representative
works here [2, 13, 15, 16, 17, 26]. Denoting the unknown source by $F(x,t)$,
to recover it, we need the observations of the solution $u$ on the whole
domain $\mathbb{R}^{d}\times(0,\infty)$, which is impractical in almost every
situation. Therefore, in the research of the inverse problem, people usually
consider some special cases of the unknown source. For instance, a popular
case is the variable separable source, i.e., $F(x,t):=p(x)q(t)$. Given the
spatial component $p(x)$ or the temporal component $q(t)$, recovering the
other unknown part is a classical field in the inverse problems. See [3, 13,
14, 22] and the references therein. Furthermore, the work [23] considered the
case when $p(x)$ and $q(t)$ are both unknown; in [17] the authors recovered
the unknown source $p(x^{\prime},t)$ where $x^{\prime}\in\mathbb{R}^{d-1}$.
Comparing with the variable separable form, the semi-discrete formulation
below simulates the general source $F(x,t)$ better:

$F(x,t):=\sum_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}}.$

We can see the above formulation can approximate the general form $F(x,t)$
accurately if the time mesh $\\{0\leq t_{0}<t_{1}<\cdots\\}$ is fine enough.
The parabolic equation with a space-time-dependent source has many
applications. For instance, in the area of medical research, one needs to
trace the blood distribution in some tissues of the human body; in the
reservoir simulation area, one example is to trace the amount of liquid
injected into oil field consisting of the impermeable rocks; in the ocean,
people may need to determine the location of a leaking oil tanker and so on.
In [19], the inverse source problem with the semi-discrete source is
investigated.

In this work, the unknown source still has semi-discrete formulation. The
measurements we used are the boundary flux data, which means the observed area
will be the subset of the boundary of the domain. To save cost, absolutely we
want the observed area to be as small as possible. In [19, 23], the authors
considered the heat equation on the two-dimensional unit disc, and proved the
uniqueness theorem under the boundary flux data from two chosen points on the
boundary. The proof depended on the explicit representation of the eigensystem
of the Laplacian $\Delta$ on the two-dimensional unit disc. The conclusions in
[19, 23] confirm that in the heat equation, if the domain has a smooth shape,
then it is possible to recover the source from sparse boundary data.
Sequentially, there is a natural question that can we solve the similar
inverse source problem in the parabolic equation with a general domain, in
which the explicit representation of the eigensystem is not applicable.

###  1.2 Mathematical statement and main theorem.

We give the mathematical model as follows:

$\begin{cases}\begin{aligned}
(\partial_{t}+\mathcal{A})u(x,t)&=\sum_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}},&&(x,t)\in\Omega\times(0,\infty),\\\
u(x,t)&=0,&&(x,t)\in\partial\Omega\times(0,\infty)\cup\Omega\times\\{0\\},\end{aligned}\end{cases}$
(1.1)

where $\Omega\subset\mathbb{R}^{d}$ is open and bounded. The operator
$\mathcal{A}$ is a symmetric elliptic operator, defined as

$\mathcal{A}\psi=-\nabla\cdot(\kappa(x)\nabla\psi)+q(x)\psi,\quad\psi\in
H^{2}(\Omega)\cap H_{0}^{1}(\Omega).$ (1.2)

Here $\kappa$ and $q$ possess appropriate regularities to support the future
proof and they satisfy

$0<C_{1}\leq\kappa(x)\leq C_{2}<\infty\ \text{and}\ q(x)\geq 0\ \text{for
a.e.}\ x\in\Omega.$

In the source term, the spatial components $\\{p_{k}(x)\\}_{k=1}^{K}$ and the
time mesh $\\{t_{k}\\}_{k=0}^{K}$ are all undetermined. The measurements we
used are the boundary flux data:

$\frac{\partial u}{\partial{\overrightarrow{\bf n}}}(x,t),\
(x,t)\in\Gamma\times(0,\infty)\subset\partial\Omega\times(0,\infty),$

where $\overrightarrow{\bf n}$ means the outward normal unit vector of the
boundary and $\Gamma\subset\partial\Omega$ is the observed area. Hence, the
interested inverse problem in this work is given as follows:

$\text{recovering}\ \\{t_{0},t_{k},p_{k}(x)\\}_{k=1}^{K}\ \text{in the source
from the data}\ \frac{\partial u}{\partial{\overrightarrow{\bf
n}}}\big{|}_{\Gamma\times(0,\infty)}.$

For this inverse problem, we have two goals: first, proving that the boundary
flux data generated from a small observed area $\Gamma$ can uniquely determine
the source; second, recovering the unknown source from the sparse boundary
data numerically. In the aspect of theoretical analysis, we attempt to build
the uniqueness theorem. Firstly, we give some prior assumptions for the semi-
discrete unknown source
$\sum_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}}$.

######  Assumption 1.1.

  * (a)

For the time mesh $\\{t_{k}\\}_{k=0}^{K}$, we have
$K\in\mathbb{N}^{+}\cup\\{\infty\\},$ and there exists $\eta>0$ such that
$\inf\\{|t_{k}-t_{k+1}|:k=0,\cdots,K-1\\}\geq\eta.$

  * (b)

$\\{p_{k}(x)\\}_{k=1}^{K}\subset L^{2}(\Omega)$,
$\|p_{k}\|_{L^{2}(\Omega)}\neq 0$ for $k=1,\cdots,K$, and
$\|p_{k}-p_{k+1}\|_{L^{2}(\Omega)}\neq 0$ for $k=1,\cdots,K-1$.

In Assumption 1.1 (a), we do not require the time mesh $\\{t_{k}\\}_{k=0}^{K}$
be finite, which means $K$ can be infinity. Also, Assumption 1.1 (b) is set to
make sure the source
$\sum\nolimits_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}}$ can not be
simplified further. Otherwise, if
$\|p_{k_{0}}\|_{L^{2}(\Omega)}=\|p_{k_{1}-1}-p_{k_{1}}\|_{L^{2}(\Omega)}=0$,
then we can write the source as

$\sum\nolimits_{k\notin\\{k_{0},k_{1}-1,k_{1}\\}}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}}+p_{k_{1}}(x)\chi_{{}_{t\in[t_{k_{1}-2},t_{k_{1}})}}.$

In the next section, we will give Assumption 2.4. With Assumptions 1.1 and
2.4, we can state the main theorem of this work.

######  Theorem 1.

Under Assumptions 1.1 and 2.4, the flux data from any nonempty open subset of
$\partial\Omega$ can uniquely determine the semi-discrete unknown source
$\sum\nolimits_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}}$.

More precisely, given two sets of unknowns
$\\{t_{0},t_{k},p_{k}(x)\\}_{k=1}^{K}$ and
$\\{\tilde{t}_{0},\tilde{t}_{k},\tilde{p}_{k}(x)\\}_{k=1}^{\tilde{K}}$, we
denote the corresponding solutions by $u$ and $\tilde{u}$ respectively, assume
Assumptions 1.1 and 2.4 be valid, and let $\Gamma\subset\partial\Omega$ be an
arbitrary nonempty open subset. Provided

$\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}\big{|}_{\Gamma\times(0,\infty)}=\frac{\partial\tilde{u}}{\partial{\overrightarrow{\bf
n}}}\big{|}_{\Gamma\times(0,\infty)},$

we have
$\\{t_{0},t_{k},p_{k}(x)\\}_{k=1}^{K}=\\{\tilde{t}_{0},\tilde{t}_{k},\tilde{p}_{k}(x)\\}_{k=1}^{\tilde{K}}$.

Theorem 1 confirms that the data from any nonempty open subset of the boundary
is sufficient to support the uniqueness of this inverse source problem.

###  1.3 Bayesian formulation and outline.

For this inverse source problem, after proving Theorem 1, it is time to
consider reconstructing the unknown source numerically. The design of the
algorithm will be challenging since there are too many unknowns in the source.
The conventional methods are hard to handle problems with high dimensional
unknowns and are very sensitive to the observation locations. However, from
the semi-discrete formulation, one needs to estimate a sequence of unknown
states $\\{p_{k}(x)\\}_{k=1}^{K}$ recursively in time. One choice is to
estimate each state $p_{k}(x)$ given a sequence of observations up to $k$ in
time by the posterior distribution. This task of sequential prediction based
on the online observations can then be categorized as the standard Bayesian
filtering problems which have been thoroughly studied in the control theory
[10, 5, 24, 20]. This motivates us to formulate the inversion problem in the
Bayesian framework.

The numerical reconstruction is a time series tracking problem and consists of
constructing the posterior probability distribution of a sequence of states
based on observations. The resulting model typically is highly nonlinear with
non-Gaussian distribution, and hence obtaining an explicit exact solution is a
hard work. This reminds us of estimating the state recursively via a filtering
approach [18]. In this work, we are going to use the particle filter (PF)
algorithm to compute the posterior of the state in a sequence. The benefit of
the particle filter algorithm is: it can handle any nonlinear properties and
is capable to solve any distribution (no Gaussian assumption) [25, 9]. The PF
algorithms usually involve sampling from a probabilistic model which depends
on the previous states and observations. The samples are called particles and
constitute an empirical approximation of the model. The trajectories are then
extended by sampling the next state particles based on the existing particles
and new observations. Degeneration issue is a concern of the PF algorithm; but
the resampling techniques can reduce the variance of the particles. In this
work, we use the systematical resampling [7, 9] and please check [7] for a
comprehensive review of the resampling schemes. Another concern of the PF
algorithms is the design of the proposal. Many techniques like the
linearization [8], extended Kalman filter, and unscented PF [25] are among the
most important options to this issue. In this work, we use the transition
prior, and the results are satisfactory. The details can be seen in Section 4.

The contribution of this work can be summarized as follow:

  1. 1.

Prove the uniqueness theorem rigorously, and confirm that the semi-discrete
source can be recovered by the flux data from any nonempty open subset of the
boundary. This conclusion can decrease the cost of the real projects and is
meaningful in practical applications.

  2. 2.

To solve the inverse problem numerically, we formulate the process as a
sequential prediction and hence use the filter algorithm to predict the semi-
discrete source. The method enables us to solve the models with a long
sequence of states and is very accurate.

Finally, we give the outline of this paper. In Section 2, we collect several
useful results and show some lemmas which will be used in future proofs. In
Section 3 we will prove the main result–Theorem 1. The proof depends on the
Laplace transform and the knowledge of complex analysis. In Section 4, we
consider the numerical reconstruction of this inverse source problem. The
algorithm is given, and some satisfactory numerical results are provided. In
Section 5, we list some future works for this inverse source problem.

##  2 Preliminaries.

###  2.1 Eigensystem of $\mathcal{A}$ on $\Omega$.

For the operator $\mathcal{A}$ defined in (1.2), we denote the eigensystem by
$\\{\lambda_{n},\varphi_{n}(x)\\}_{n=1}^{\infty}$. Since $\mathcal{A}$ is
positive definite and self-adjoint, we have that
$0<\lambda_{1}\leq\lambda_{2}\leq\cdots$ and $\lambda_{n}\to\infty$ as
$n\to\infty$, and the set of eigenfunctions
$\\{\varphi_{n}(x)\\}_{n=1}^{\infty}$ constitutes an orthonormal basis of the
space $L^{2}(\Omega)$. Furthermore, if $\varphi_{n}$ is an eigenfunction of
$\mathcal{A}$ corresponding to $\lambda_{n}$, so is $\overline{\varphi_{n}}$.
Here $\overline{\varphi_{n}}$ means the complex conjugate of $\varphi_{n}$.
This gives that the set $\\{\varphi_{n}(x)\\}_{n=1}^{\infty}$ coincides with
$\\{\overline{\varphi_{n}(x)}\\}_{n=1}^{\infty}$. Also, from the trace
theorem, we have $\\{\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\\}_{n=1}^{\infty}\subset H^{1/2}(\partial\Omega)$, and
the next lemma concerns the non-vanishing property of
$\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\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{\overrightarrow{\bf
n}}}$ can not vanish almost everywhere on $\Gamma$.

###### Proof.

Assume that there exists $k\in\mathbb{N}^{+}$ and a nonempty open subset
$\Gamma\subset\partial\Omega$ such that
$\frac{\partial\varphi_{k}}{\partial{\overrightarrow{\bf n}}}=0$ a.e. on
$\Gamma$. This implies that we can find $r>0$ and $x_{0}\in\partial\Omega$
such that $\frac{\partial\varphi_{k}}{\partial{\overrightarrow{\bf n}}}\equiv
0$ a.e. on $B(x_{0},r)\cap\partial\Omega$, where $B(x_{0},r)$ is the ball
centered at $x_{0}$ with radius $r$ in $\mathbb{R}^{d}.$ Now define the
extension $\Omega_{e}$ of $\Omega$ as $\Omega_{e}=\Omega\cup B(x_{0},r)$ and
the extended function $\varphi_{k,e}$ of $\varphi_{k}$ on $\Omega_{e}$ as

$\varphi_{k,e}(x)=\begin{cases}\begin{aligned}
&\varphi_{k}(x),&&x\in\Omega,\\\ &0,&&x\in
B(x_{0},r)\setminus\Omega.\end{aligned}\end{cases}$

Obviously, $\varphi_{k,e}\in H^{1}(\Omega_{e})$. Now for any $\psi\in
C_{c}^{\infty}(\Omega_{e}),$ we have

$\displaystyle\int_{\Omega_{e}}(\kappa\nabla\varphi_{k,e}\cdot\nabla\psi-\lambda_{k}\varphi_{k,e}\psi)\
dx$
$\displaystyle=\int_{\Omega}(\kappa\nabla\varphi_{k,e}\cdot\nabla\psi-\lambda_{k}\varphi_{k,e}\psi)\
dx+\int_{B(x_{0},r)\setminus\Omega}0\ dx$
$\displaystyle=\int_{\Omega}[(\mathcal{A}-\lambda_{k})\varphi_{k}]\psi\ dx=0.$

Hence, $\varphi_{k,e}\in H^{1}(\Omega_{e})$ is a weak solution of the equation

$(\mathcal{A}-\lambda_{k})v(x)=0,\ x\in\Omega_{e},$

and we can find $B(x_{0}^{\prime},r^{\prime})\subset
B(x_{0},r)\subset\Omega_{e}$ s.t. $\varphi_{k,e}\equiv 0$ on
$B(x_{0}^{\prime},r^{\prime})$. By the unique continuation principle, we have
$\varphi_{k,e}\equiv 0$ on $\Omega_{e},$ i.e. $\varphi_{k}$ vanishes on
$\Omega,$ which contradicts with the fact that
$\|\varphi_{k}\|_{L^{2}(\Omega)}=1.$ Therefore we conclude that
$\frac{\partial\varphi_{k}}{\partial{\overrightarrow{\bf n}}}$ does not vanish
almost everywhere on any nonempty open subset of $\partial\Omega$, and the
proof is complete. ∎

###  2.2 The set $\\{\xi_{l}\\}_{l=1}^{\infty}$.

In this part, we build the auxiliary set of functions
$\\{\xi_{l}\\}_{l=1}^{\infty}$ which will be used later. For
$\psi_{1},\psi_{2}\in L^{2}(\Omega),\ \psi_{3},\psi_{4}\in
L^{2}(\partial\Omega),$ we define

$\langle\psi_{1},\psi_{2}\rangle_{\Omega}=\int_{\Omega}\psi_{1}(\cdot)\overline{\psi_{2}(\cdot)}\
dx,\ \
\langle\psi_{3},\psi_{4}\rangle_{\kappa,\partial\Omega}=\int_{\partial\Omega}\kappa(\cdot)\psi_{3}(\cdot)\overline{\psi_{4}(\cdot)}\
dx,$

and accordingly we give the induced norm
$\|\cdot\|_{L^{2}(\Omega)},\|\cdot\|_{L^{2}(\kappa,\partial\Omega)}$ as

$\|\cdot\|^{2}_{L^{2}(\Omega)}=\langle\cdot,\cdot\rangle_{\Omega},\
\|\cdot\|^{2}_{L^{2}(\kappa,\partial\Omega)}=\langle\cdot,\cdot\rangle_{\kappa,\partial\Omega}.$

Not hard to see $\|\cdot\|_{L^{2}(\kappa,\partial\Omega)}$ is equivalent to
the usual $L^{2}$ norm $\|\cdot\|_{L^{2}(\partial\Omega)}$. The next lemma
considers a density property.

######  Lemma 2.2.

For the eigenfunctions $\\{\varphi_{n}\\}_{n=1}^{\infty}$ defined in Section
2.1, we have
$\text{Span}\\{\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\\}_{n=1}^{\infty}$ is dense in
$L^{2}(\kappa,\partial\Omega)$.

###### Proof.

With the density of $H^{3/2}(\partial\Omega)$ in $L^{2}(\partial\Omega)$ under
the norm $\|\cdot\|_{L^{2}(\partial\Omega)}$, it is sufficient to prove that
$\tilde{\psi}\in H^{3/2}(\partial\Omega)$ vanishes almost everywhere on
$\partial\Omega$ if
$\langle\tilde{\psi},\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\rangle_{\kappa,\partial\Omega}=0$ for each
$n\in\mathbb{N}^{+}$.

We let $\psi$ satisfy the following system:

$\begin{cases}\begin{aligned} \mathcal{A}\psi(x)&=0,&&x\in\Omega,\\\
\psi&=\tilde{\psi},&&x\in\partial\Omega.\end{aligned}\end{cases}$

The regularity $\tilde{\psi}\in H^{3/2}(\partial\Omega)$ gives that $\psi\in
H^{2}(\Omega)$, then the Green’s identity can be used. Fixing one
$n\in\mathbb{N}^{+}$, we have

$\langle\mathcal{A}\psi,\varphi_{n}\rangle_{\Omega}-\langle\psi,\mathcal{A}\varphi_{n}\rangle_{\Omega}=\langle\tilde{\psi},\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}\rangle_{\kappa,\partial\Omega}-\langle\frac{\partial\psi}{\partial{\overrightarrow{\bf
n}}},\varphi_{n}\rangle_{\kappa,\partial\Omega}.$

With the results $\mathcal{A}\psi=0$ on $\Omega$, $\varphi_{n}=0$ on
$\partial\Omega$ and
$\langle\tilde{\psi},\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}\rangle_{\kappa,\partial\Omega}=0$, the above equality gives that

$\langle\psi,\mathcal{A}\varphi_{n}\rangle_{\Omega}=\lambda_{n}\langle\psi,\varphi_{n}\rangle_{\Omega}=0.$

So we see that for each $n\in\mathbb{N}^{+}$,
$\langle\psi,\varphi_{n}\rangle_{\Omega}=0$. This together with the
completeness of $\\{\varphi_{n}\\}_{n=1}^{\infty}$ in $L^{2}(\Omega)$ yields
that $\|\psi\|_{L^{2}(\Omega)}=0$. From the definition of weak derivative, it
is not hard to see that the first-order weak derivative of $\psi$ is also
vanishing on $\Omega$. This gives that $\|\psi\|_{H^{1}(\Omega)}=0$. By the
continuity of the trace operator, it holds that

$\|\psi|_{\partial\Omega}\|_{L^{2}(\partial\Omega)}=\|\tilde{\psi}\|_{L^{2}(\partial\Omega)}\leq
C\|\psi\|_{H^{1}(\Omega)}=0.$

So we have $\tilde{\psi}=0$ almost everywhere on $\partial\Omega$ and the
proof is complete. ∎

Now we can build the auxiliary set $\\{\xi_{l}\\}_{l=1}^{\infty}$. With the
density of
$\text{Span}\\{\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\\}_{n=1}^{\infty}$ in $L^{2}(\kappa,\partial\Omega)$,
we can construct the orthonormal basis $\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$
in $L^{2}(\kappa,\partial\Omega)$ as follows. Firstly we set
$\tilde{\xi}_{1}=\frac{\partial\varphi_{1}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}/\|\frac{\partial\varphi_{1}}{\partial{\overrightarrow{\bf
n}}}\|_{L^{2}(\kappa,\partial\Omega)},$ and for $l=2,3,\cdots$, we assume that
the orthonormal set $\\{\tilde{\xi}_{j}\\}_{j=1}^{l-1}$ has been constructed.
Then we let $n_{l}\in\mathbb{N}^{+}$ be the smallest number such that
$\frac{\partial\varphi_{n_{l}}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\notin\text{Span}\\{\tilde{\xi}_{j}\\}_{j=1}^{l-1}$, and
after that we can pick
$\tilde{\xi}_{l}\in\text{Span}\\{\frac{\partial\varphi_{n_{l}}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega},\ \tilde{\xi}_{1},\cdots,\tilde{\xi}_{l-1}\\}$ such
that

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

The density of
$\text{Span}\\{\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\\}_{n=1}^{\infty}$ in $L^{2}(\kappa,\partial\Omega)$
ensures that $\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$ is an orthonormal basis in
$L^{2}(\kappa,\partial\Omega)$. Also, from the smoothness property
$\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}|_{\partial\Omega}\in H^{1/2}(\partial\Omega)$ we have
$\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}\subset H^{1/2}(\partial\Omega)$.

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

$\begin{cases}\begin{aligned} \mathcal{A}\xi_{l}(x)&=0,&&x\in\Omega,\\\
\xi_{l}&=\tilde{\xi}_{l},&&x\in\partial\Omega.\end{aligned}\end{cases}$ (2.1)

Then we get the auxiliary set $\\{\xi_{l}\\}_{l=1}^{\infty}$ on $\Omega$.

###  2.3 The coefficients $\\{c_{z,n}\\}$ and Assumption 2.4.

The proof of Theorem 1 relies on Assumption 2.4. Before to state Assumption
2.4, we need to build the coefficients $\\{c_{z,n}\\}$. For a fixed point
$z\in\partial\Omega$, we define the series $\psi_{z}^{N}\in
H^{1}(\overline{\Omega})$ as

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

where $\\{\tilde{\xi}_{l},\xi_{l}\\}_{l=1}^{\infty}$ are defined in the above
subsection. Then we consider the following system and denote the solutions by
$\\{u_{z}^{N}\\}$:

$\begin{cases}\begin{aligned}
(\partial_{t}+\mathcal{A})u_{z}^{N}(x,t)&=0,&&(x,t)\in\Omega\times(0,\infty),\\\
u_{z}^{N}(x,t)&=0,&&(x,t)\in\partial\Omega\times(0,\infty),\\\
u_{z}^{N}(x,0)&=-\psi_{z}^{N},&&x\in\Omega.\end{aligned}\end{cases}$ (2.3)

The finite summation $\psi_{z}^{N}$ is constructed following the role of Dirac
delta function, which reflects the information of the targeted function on a
specific point via integration. With the following lemma, we can give the
coefficients $\\{c_{z,n}\\}$.

######  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_{\Omega}$
exists.

###### Proof.

Firstly from Green’s identities we have

$\displaystyle\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{\Omega}$
$\displaystyle=\lambda_{n}^{-1}\langle\psi_{z}^{N},\mathcal{A}\overline{\varphi_{n}}\rangle_{\Omega}$
$\displaystyle=\lambda_{n}^{-1}\Big{(}\langle\mathcal{A}\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{\Omega}-\langle\psi_{z}^{N},\frac{\partial\overline{\varphi_{n}}}{\partial{\overrightarrow{\bf
n}}}\rangle_{\kappa,\partial\Omega}+\langle\frac{\partial\psi_{z}^{N}}{\partial{\overrightarrow{\bf
n}}},\overline{\varphi_{n}}\rangle_{\kappa,\partial\Omega}\Big{)}$
$\displaystyle=-\lambda_{n}^{-1}\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\langle\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}},\tilde{\xi}_{l}\rangle_{\kappa,\partial\Omega}=:c_{z,n}^{N},$

where the system (2.1) and the boundary condition of $\varphi_{n}$ are used.
From the definition of $\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$, we see that
$\langle\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}},\tilde{\xi}_{l}\rangle_{\kappa,\partial\Omega}=0$ for large $l$. Hence if
$N$ is sufficiently large, the value
$\lambda_{n}^{-1}\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\langle\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}},\tilde{\xi}_{l}\rangle_{\kappa,\partial\Omega}$ will not depend on the
value of $N$. This means that
$\lim_{N\to\infty}\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{\Omega}$
exists and the proof is complete. ∎

Now we define

$c_{z,n}:=\lim_{N\to\infty}\langle\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{\Omega}=\lim_{N\to\infty}c_{z,n}^{N},\quad
p_{k,n}:=\langle p_{k}(\cdot),\varphi_{n}(\cdot)\rangle_{\Omega},$ (2.4)

it is time to state Assumption 2.4.

######  Assumption 2.4.

  * (a)

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

  * (b)

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

  * (c)

For a.e. $t\in[0,\infty)$, the spatial components $\\{p_{k}(x)\\}_{k=1}^{K}$
are smooth enough such that the series
$\sum_{l=1}^{\infty}\tilde{\xi}_{l}(x)\langle\frac{\partial
u}{\partial{\overrightarrow{\bf
n}}}(\cdot,t),\tilde{\xi}_{l}(\cdot)\rangle_{\kappa,\partial\Omega}$ converges
to $\frac{\partial u}{\partial{\overrightarrow{\bf n}}}(x,t)$ pointwisely for
a.e. $x\in\partial\Omega$.

##  3 Uniqueness theorem.

###  3.1 Representation for the boundary flux data.

In this part, we will consider how to present the boundary flux
$\frac{\partial u}{\partial{\overrightarrow{\bf n}}}|_{\partial\Omega}$ by the
unknown source. Firstly, we give the next lemma.

######  Lemma 3.1.

For a.e. $z\in\partial\Omega$, we define $w_{z}^{N}=u_{z}^{N}+\psi_{z}^{N},$
where $u_{z}^{N}$ and $\psi_{z}^{N}$ are given by (2.3) and (2.2),
respectively. Then we have:

$\displaystyle\text{for a.e.}\ t\in[0,t_{0}],$
$\displaystyle\qquad\qquad-\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}(z,t)=0;$ $\displaystyle\text{for a.e.}\ t\in(t_{m_{0}-1},t_{m_{0}}]\
\text{with}\ 1\leq m_{0}\leq K,$ $\displaystyle\qquad\qquad-\frac{\partial
u}{\partial{\overrightarrow{\bf
n}}}(z,t)=\lim_{N\to\infty}\int_{\Omega}p_{m_{0}}(x)w_{z}^{N}(x,t-t_{m_{0}-1})\
dx$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\sum_{m=1}^{m_{0}-1}\lim_{N\to\infty}\int_{\Omega}p_{m}(x)[w_{z}^{N}(x,t-t_{m-1})-w_{z}^{N}(x,t-t_{m})]\
dx;$ $\displaystyle\text{for a.e.}\ t\in(t_{K},\infty)\ \text{(when K is
finite)},$ $\displaystyle\qquad\qquad-\frac{\partial
u}{\partial{\overrightarrow{\bf
n}}}(z,t)=\sum_{m=1}^{K}\lim_{N\to\infty}\int_{\Omega}p_{m}(x)[w_{z}^{N}(x,t-t_{m-1})-w_{z}^{N}(x,t-t_{m})]\
dx.$

###### Proof.

From the definition of $\psi_{z}^{N}$, we see that
$(\partial_{t}+\mathcal{A})\psi_{z}^{N}=0.$ This result and (2.3) show that
$w_{z}^{N}$ satisfies the equation

$(\partial_{t}+\mathcal{A})w_{z}^{N}=0,\ (x,t)\in\Omega\times(0,\infty),$

with the boundary condition
$w_{z}^{N}|_{\partial\Omega}=\psi_{z}^{N}|_{\partial\Omega}$ and the initial
condition $w_{z}^{N}(x,0)=0$. Therefore, Green’s identities give that for each
$v\in H_{0}^{1}(\Omega)$,

$\int_{\Omega}(\partial_{t}+q)w_{z}^{N}(x,t)v(x)+\kappa(x)\nabla
w_{z}^{N}(x,t)\cdot\nabla v(x)\ dx=0,\ \ t\in(0,\infty).$

From (1.1), we obtain

$\displaystyle\int_{0}^{t}\int_{\Omega}F(x,\tau)w_{z}^{N}(x,t-\tau)\ dx\
d\tau=\int_{0}^{t}\int_{\Omega}(\partial_{t}+\mathcal{A})u(x,\tau)\
w_{z}^{N}(x,t-\tau)\ dx\ d\tau.$

Green’s identities and the vanishing initial conditions of $u$ and $w_{z}^{N}$
give that

$\displaystyle\int_{0}^{t}\int_{\Omega}(\partial_{t}+q)u(x,\tau)\
w_{z}^{N}(x,t-\tau)\ dx\ d\tau=$
$\displaystyle\int_{0}^{t}\int_{\Omega}(\partial_{t}+q)w_{z}^{N}(x,t-\tau)\
u(x,\tau)\ dx\ d\tau,$
$\displaystyle\int_{0}^{t}\int_{\Omega}-\nabla\cdot(\kappa\nabla u(x,\tau))\
w_{z}^{N}(x,t-\tau)\ dx\ d\tau=$
$\displaystyle\int_{0}^{t}\int_{\Omega}\kappa(x)\nabla u(x,\tau)\cdot\nabla
w_{z}^{N}(x,t-\tau)\ dx\ d\tau$
$\displaystyle-\int_{0}^{t}\int_{\partial\Omega}\kappa(x)\frac{\partial
u}{\partial{\overrightarrow{\bf n}}}(x,\tau)\psi_{z}^{N}(x)\ dx\ d\tau.$

Hence we have

$\displaystyle\int_{0}^{t}\int_{\Omega}F(x,\tau)w_{z}^{N}(x,t-\tau)\ dx\
d\tau$
$\displaystyle=\int_{0}^{t}\int_{\Omega}(\partial_{t}+q)w_{z}^{N}(x,t-\tau)\
u(x,\tau)+\kappa(x)\nabla w_{z}^{N}(x,t-\tau)\cdot\nabla u(x,\tau)\ dx\ d\tau$
$\displaystyle\quad-\int_{0}^{t}\int_{\partial\Omega}\kappa(x)\frac{\partial
u}{\partial{\overrightarrow{\bf n}}}(x,\tau)\psi_{z}^{N}(x)\ dx\ d\tau$
$\displaystyle=-\int_{0}^{t}\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\
\big{\langle}\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}(\cdot,\tau),\tilde{\xi}_{l}(\cdot)\big{\rangle_{\kappa,\partial\Omega}}\
d\tau.$

For $t\in(t_{m_{0}-1},t_{m_{0}}]$ with $1\leq m_{0}\leq K$, the left side of
the above equality can be written as

$\displaystyle\int_{0}^{t}\int_{\Omega}F(x,t-\tau)w_{z}^{N}(x,\tau)\ dx\
d\tau=$
$\displaystyle\int_{0}^{t-t_{m_{0}-1}}\int_{\Omega}p_{m_{0}}(x)w_{z}^{N}(x,\tau)\
dx\ d\tau$
$\displaystyle+\sum_{m=1}^{m_{0}-1}\int_{t-t_{m}}^{t-t_{m-1}}\int_{\Omega}p_{m}(x)w_{z}^{N}(x,\tau)\
dx\ d\tau.$

Similarly, if $t>t_{K}$, we have

$\displaystyle\int_{0}^{t}\int_{\Omega}F(x,t-\tau)w_{z}^{N}(x,\tau)\ dx\
d\tau=\sum_{m=1}^{K}\int_{t-t_{m}}^{t-t_{m-1}}\int_{\Omega}p_{m}(x)w_{z}^{N}(x,\tau)\
dx\ d\tau.$

Then for the above equality, differentiating at time $t$ gives that for
$t\in(t_{m_{0}-1},t_{m_{0}}]$ with $1\leq m_{0}\leq K$,

$\displaystyle-\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\ \big{\langle}\frac{\partial
u}{\partial{\overrightarrow{\bf
n}}}(\cdot,t),\tilde{\xi}_{l}(\cdot)\big{\rangle_{\kappa,\partial\Omega}}=$
$\displaystyle\int_{\Omega}p_{m_{0}}(x)w_{z}^{N}(x,t-t_{m_{0}-1})\ dx$
$\displaystyle+\sum_{m=1}^{m_{0}-1}\int_{\Omega}p_{m}(x)[w_{z}^{N}(x,t-t_{m-1})-w_{z}^{N}(x,t-t_{m})]\
dx,$

and for $t>t_{K}$ (From Assumption 1.1, we see that $K$ is finite),

$-\sum_{l=1}^{N}\tilde{\xi}_{l}(z)\ \big{\langle}\frac{\partial
u}{\partial{\overrightarrow{\bf
n}}}(\cdot,t),\tilde{\xi}_{l}(\cdot)\big{\rangle_{\kappa,\partial\Omega}}=\sum_{m=1}^{K}\int_{\Omega}p_{m}(x)[w_{z}^{N}(x,t-t_{m-1})-w_{z}^{N}(x,t-t_{m})]\
dx.$

From Assumption 2.4 $(c)$ we have for a.e. $t>t_{K}$,

$-\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}(z,t)=\sum_{m=1}^{K}\lim_{N\to\infty}\int_{\Omega}p_{m}(x)[w_{z}^{N}(x,t-t_{m-1})-w_{z}^{N}(x,t-t_{m})]\
dx,$

and for a.e. $t\in(t_{m_{0}-1},t_{m_{0}}]$ with $1\leq m_{0}\leq K$,

$\displaystyle-\frac{\partial u}{\partial{\overrightarrow{\bf n}}}(z,t)=$
$\displaystyle\lim_{N\to\infty}\int_{\Omega}p_{m_{0}}(x)w_{z}^{N}(x,t-t_{m_{0}-1})\
dx$
$\displaystyle+\sum_{m=1}^{m_{0}-1}\lim_{N\to\infty}\int_{\Omega}p_{m}(x)[w_{z}^{N}(x,t-t_{m-1})-w_{z}^{N}(x,t-t_{m})]\
dx.$

Assumption 1.1 $(a)$ ensures that for a fixed $t$, the summations in the right
sides of the above results are finite. Therefore, the order of summation and
limit can be exchanged. Moreover, the proof for $t\in[0,t_{0}]$ is trivial.
The proof is complete. ∎

From the above lemma, the following corollary can be deduced.

######  Corollary 3.2.

For a.e. $z\in\partial\Omega$ and $t\in(0,\infty)$, it holds that

$-\int_{0}^{t}\frac{\partial u}{\partial{\overrightarrow{\bf n}}}(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^{-\lambda_{n}\tau})\Big{]}\
d\tau.$

Here $p_{k,n}$ and $c_{z,n}$ are defined in (2.4).

###### Proof.

Firstly, let us evaluate the limit $\lim_{N\to\infty}\langle
p_{k}(\cdot),\overline{w_{z}^{N}(\cdot,t)}\rangle_{\Omega}$ for $k=1,\cdots,K$
and $t\in(0,\infty)$. Fixing $N\in\mathbb{N}^{+}$, from the definition of
$\psi_{z}^{N}$ we have $\psi_{z}^{N}\in L^{2}(\Omega)$. So the Fourier
expansion of $\psi_{z}^{N}$ can be given as
$\psi_{z}^{N}=\sum_{n=1}^{\infty}c_{z,n}^{N}\overline{\varphi_{n}}$, where
$c_{z,n}^{N}$ is defined in the proof of Lemma 2.4. Moreover, from
$w_{z}^{N}=u_{z}^{N}+\psi_{z}^{N}$ and equation (2.3), we have

$w_{z}^{N}(x,t)=\sum_{n=1}^{\infty}c_{z,n}^{N}(1-e^{-\lambda_{n}t})\overline{\varphi_{n}(x)}.$

The above result together with the regularity $\psi_{z}^{N}\in L^{2}(\Omega)$
yields that $w_{z}^{N}(\cdot,t)\in L^{2}(\Omega)$ for $t\in[0,\infty)$. Also
recall that $p_{k}\in L^{2}(\Omega)$, then

$\langle
p_{k}(\cdot),\overline{w_{z}^{N}(\cdot,t)}\rangle_{\Omega}=\sum_{n=1}^{\infty}c_{z,n}^{N}p_{k,n}(1-e^{-\lambda_{n}t}).$

The condition in Assumption 2.4 implies that the Dominated Convergence Theorem
can be applied to the above series. So we have

$\lim_{N\to\infty}\sum_{n=1}^{\infty}c_{z,n}^{N}p_{k,n}(1-e^{-\lambda_{n}t})=\sum_{n=1}^{\infty}\lim_{N\to\infty}c_{z,n}^{N}p_{k,n}(1-e^{-\lambda_{n}t}).$

From (2.4), we have

$\lim_{N\to\infty}\langle
p_{k}(\cdot),\overline{w_{z}^{N}(\cdot,t)}\rangle_{\Omega}=\lim_{N\to\infty}\sum_{n=1}^{\infty}c_{z,n}^{N}p_{k,n}(1-e^{-\lambda_{n}t})=\sum_{n=1}^{\infty}c_{z,n}p_{k,n}(1-e^{-\lambda_{n}t}).$

Now we claim that for $t>0$,

$-\int_{0}^{t}\frac{\partial u}{\partial{\overrightarrow{\bf n}}}(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^{-\lambda_{n}\tau})\Big{]}\
d\tau.$ (3.1)

If $t\leq t_{0}$, from Lemma 3.1 we have $\frac{\partial
u}{\partial{\overrightarrow{\bf n}}}(z,\tau)=0$ for $\tau\leq t$. This means
(3.1) is valid.

For $t\in(t_{m_{0}-1},t_{m_{0}}]$ with $1\leq m_{0}\leq K$, Lemma 3.1 yields
that

$\displaystyle-\int_{0}^{t}\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}(z,\tau)\ d\tau=$
$\displaystyle-\sum_{m=1}^{m_{0}-1}\int_{t_{m-1}}^{t_{m}}\frac{\partial
u}{\partial{\overrightarrow{\bf n}}}(z,\tau)\
d\tau-\int_{t_{m_{0}-1}}^{t}\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}(z,\tau)\ d\tau$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{m_{0}-1}\int_{t_{m-1}}^{t_{m}}\sum_{l=1}^{m-1}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(e^{-\lambda_{n}(\tau-
t_{l})}-e^{-\lambda_{n}(\tau-t_{l-1})})\ d\tau$
$\displaystyle+\sum_{m=1}^{m_{0}-1}\int_{t_{m-1}}^{t_{m}}\sum_{n=1}^{\infty}c_{z,n}p_{m,n}(1-e^{-\lambda_{n}(\tau-
t_{m-1})})\ d\tau$
$\displaystyle+\int_{t_{m_{0}-1}}^{t}\sum_{l=1}^{m_{0}-1}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(e^{-\lambda_{n}(\tau-
t_{l})}-e^{-\lambda_{n}(\tau-t_{l-1})})\ d\tau$
$\displaystyle+\int_{t_{m_{0}-1}}^{t}\sum_{n=1}^{\infty}c_{z,n}p_{m_{0},n}(1-e^{-\lambda_{n}(\tau-
t_{m_{0}-1})})\ d\tau$ $\displaystyle=:$ $\displaystyle
S_{1}+S_{2}+S_{3}+S_{4}.$

The straightforward computation gives that

$\displaystyle S_{1}=$
$\displaystyle\sum_{m=1}^{m_{0}-1}\sum_{l=1}^{m-1}\int_{t_{m-1}-t_{l-1}}^{t_{m}-t_{l-1}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau-\sum_{m=1}^{m_{0}-1}\sum_{l=1}^{m-1}\int_{t_{m-1}-t_{l}}^{t_{m}-t_{l}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{m_{0}-2}\sum_{m=l+1}^{m_{0}-1}\int_{t_{m-1}-t_{l-1}}^{t_{m}-t_{l-1}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau-\sum_{l=1}^{m_{0}-2}\sum_{m=l+1}^{m_{0}-1}\int_{t_{m-1}-t_{l}}^{t_{m}-t_{l}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{m_{0}-2}\int_{t_{l}-t_{l-1}}^{t_{m_{0}-1}-t_{l-1}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau-\sum_{l=1}^{m_{0}-2}\int_{0}^{t_{m_{0}-1}-t_{l}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau.$

Similarly, we have

$\displaystyle S_{2}=$
$\displaystyle\sum_{m=1}^{m_{0}-1}\int_{0}^{t_{m}-t_{m-1}}\sum_{n=1}^{\infty}c_{z,n}p_{m,n}(1-e^{-\lambda_{n}\tau})\
d\tau,$ $\displaystyle S_{3}=$
$\displaystyle\sum_{l=1}^{m_{0}-1}\int_{t_{m_{0}-1}-t_{l-1}}^{t-t_{l-1}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau-\sum_{l=1}^{m_{0}-1}\int_{t_{m_{0}-1}-t_{l}}^{t-t_{l}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau,$ $\displaystyle S_{4}=$
$\displaystyle\int_{0}^{t-t_{m_{0}-1}}\sum_{n=1}^{\infty}c_{z,n}p_{m_{0},n}(1-e^{-\lambda_{n}\tau})\
d\tau=\int_{0}^{t}\chi_{{}_{t-\tau\in[t_{m_{0}-1},t_{m_{0}})}}\sum_{n=1}^{\infty}c_{z,n}p_{m_{0},n}(1-e^{-\lambda_{n}\tau})\
d\tau.$

These give that

$\displaystyle S_{1}+S_{2}+S_{3}=$
$\displaystyle\sum_{l=1}^{m_{0}-2}\int_{0}^{t-t_{l-1}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau-\sum_{l=1}^{m_{0}-2}\int_{0}^{t-t_{l}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau$
$\displaystyle+\int_{0}^{t-t_{m_{0}-2}}\sum_{n=1}^{\infty}c_{z,n}p_{m_{0}-1,n}(1-e^{-\lambda_{n}\tau})\
d\tau-\int_{0}^{t-t_{m_{0}-1}}\sum_{n=1}^{\infty}c_{z,n}p_{m_{0}-1,n}(1-e^{-\lambda_{n}\tau})\
d\tau$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{m_{0}-1}\int_{t-t_{l}}^{t-t_{l-1}}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\
d\tau$ $\displaystyle=$
$\displaystyle\int_{0}^{t}\sum_{l=1}^{m_{0}-1}\chi_{{}_{t-\tau\in[t_{l-1},t_{l})}}\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\Big{]}\
d\tau.$

Hence, it holds that

$\displaystyle-\int_{0}^{t}\frac{\partial u}{\partial{\overrightarrow{\bf
n}}}(z,\tau)\ d\tau$ $\displaystyle=S_{1}+S_{2}+S_{3}+S_{4}$
$\displaystyle=\int_{0}^{t}\sum_{l=1}^{m_{0}}\chi_{{}_{t-\tau\in[t_{l-1},t_{l})}}\Big{[}\sum_{n=1}^{\infty}c_{z,n}p_{l,n}(1-e^{-\lambda_{n}\tau})\Big{]}\
d\tau,$

which confirms the claim.

When $K$ is finite, we can pick $t>0$ such that $t>t_{K}$. In this case, the
term $S_{4}$ is eliminated, and $m_{0}-1$ is replaced by $K$ in the terms
$S_{1},S_{2},S_{3}$. So the claim (3.1) is still valid.

Above all, we have proved the claim (3.1) and the proof is complete. ∎

###  3.2 The analysis of Laplace transform.

The proof of Theorem 1 will rely on the Laplace transform $\mathcal{L}$,
defined as

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

Recalling the result in Corollary 3.2, it is not hard to see that

$\mathcal{L}\\{1-e^{-\lambda_{n}t}\\}(s)=\lambda_{n}s^{-1}(s+\lambda_{n})^{-1},\
\operatorname{Re}s>0.$

Also, Assumption 2.4 gives 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^{-\lambda_{n}\tau})\Big{]}\
d\tau\Big{|}\leq 2t\sum_{k=1}^{K}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|\leq Ct,$

and we have $|e^{-st}t|$ is integrable on $(0,\infty)$ for
$\operatorname{Re}s>0$. These mean the Dominated Convergence Theorem can be
used and it gives that for $\operatorname{Re}s>0,$

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

Hence, we see that for $\operatorname{Re}s>0$,

$s^{2}\mathcal{L}\Big{\\{}-\int_{0}^{t}\frac{\partial
u}{\partial{\overrightarrow{\bf n}}}(z,\tau)\
d\tau\Big{\\}}(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}\lambda_{n}(s+\lambda_{n})^{-1}\Big{]}.$
(3.2)

We index the set of distinct eigenvalues as $\\{\lambda_{j}\\}_{j=1}^{\infty}$
with increasing order. Then we give the next lemma, which contains the well-
definedness and the analyticity of the above complex-valued series.

######  Lemma 3.3.

Under Assumption 2.4, we have the following properties:

  * (a)

For $R>0$ we define $\Lambda_{R}:=\\{s\in\mathbb{C}:|s|<R\\}$, then for each
$k\in\\{1,\cdots,K\\}$,
$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}$ is
uniformly convergent for
$s\in\Lambda_{R}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$;

  * (b)

$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}$ is
holomorphic on $\mathbb{C}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$;

  * (c)

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}\lambda_{n}(s+\lambda_{n})^{-1}\big{]}$
is analytic on
$\mathbb{C}^{+}:=\\{s\in\mathbb{C}:\operatorname{Re}s>0\\}$.

###### Proof.

For $(a)$, fixing $R>0$, since $\lim_{n\to\infty}\lambda_{n}=\infty$, then we
can find large enough $N_{1}\in\mathbb{N}^{+}$ such that $\lambda_{n}\geq 2R$
if $n>N_{1}$. For $s\in\Lambda_{R}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$
and $n>N_{1}$,

$|s+\lambda_{n}|\geq|\operatorname{Re}s+\lambda_{n}|=\operatorname{Re}s+\lambda_{n}\geq\lambda_{n}-R.$

This means that $|\lambda_{n}(s+\lambda_{n})^{-1}|\leq 2$. Given $\epsilon>0$,
from Assumption 2.4, we can find $N_{2}>0$ such that
$\sum_{n=n_{0}}^{\infty}|c_{z,n}p_{k,n}|<\epsilon/2$ if $n_{0}>N_{2}$. Hence,
for $n_{0}>\max\\{N_{1},N_{2}\\}$, we have

$\Big{|}\sum_{n=n_{0}}^{\infty}c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}\Big{|}\leq
2\sum_{n=n_{0}}^{\infty}|c_{z,n}p_{k,n}|<\epsilon\ \text{for}\
s\in\Lambda_{R}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}.$

This gives the uniform convergence.

For $(b)$, fixing $R>0$, it is clear that
$c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}$ is holomorphic on
$\Lambda_{R}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$. Then the uniform
convergence gives that
$\sum_{n=1}^{\infty}c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}$ is
holomorphic on $\Lambda_{R}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$. For
each $s_{0}\in\mathbb{C}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$, we can
find a sufficiently large $R>0$ such that
$s_{0}\in\Lambda_{R}\setminus\\{-\lambda_{j}\\}_{j=1}^{\infty}$. This means
that $\sum_{n=1}^{\infty}c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}$ is
holomorphic on $s=s_{0}$, which leads to the desired result.

For $(c)$, for $s\in\mathbb{C}^{+}$, we have $|e^{-t_{k-1}s}-e^{-t_{k}s}|\leq
2$ and $|\lambda_{n}(s+\lambda_{n})^{-1}|\leq 1$. Also Assumption 2.4 gives
that $\sum_{k=1}^{K}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|<\infty$, then
following the proof for $(a)$, we obtain that 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}\lambda_{n}(s+\lambda_{n})^{-1}\Big{]}$

is uniformly convergent on $\mathbb{C}^{+}$. Sequentially, from the proof for
$(b)$, the holomorphic result can be derived and the proof is complete. ∎

###  3.3 Some auxiliary lemmas.

Before to show Theorem 1, we list and prove several auxiliary lemmas.

######  Lemma 3.4.

Recall that $\\{\lambda_{j}\\}_{j=1}^{\infty}$ is the set of distinct
eigenvalues with increasing order. 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\\}$.

###### Proof.

Fixing $j\in\mathbb{N}^{+}$, from the proof of Lemma 2.4, we have

$\displaystyle c_{z,n}$
$\displaystyle=\lambda_{j}^{-1}\lim_{N\to\infty}\langle\psi_{z}^{N},\mathcal{A}\overline{\varphi_{n}}\rangle_{\Omega}$
$\displaystyle=\lambda_{j}^{-1}\lim_{N\to\infty}\Big{(}\langle\mathcal{A}\psi_{z}^{N},\overline{\varphi_{n}}\rangle_{\Omega}-\langle\psi_{z}^{N},\frac{\partial\overline{\varphi_{n}}}{\partial{\overrightarrow{\bf
n}}}\rangle_{\kappa,\partial\Omega}+\langle\frac{\partial\psi_{z}^{N}}{\partial{\overrightarrow{\bf
n}}},\overline{\varphi_{n}}\rangle_{\kappa,\partial\Omega}\Big{)}$
$\displaystyle=-\lambda_{j}^{-1}\sum_{l=1}^{\infty}\tilde{\xi}_{l}(z)\langle\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}},\tilde{\xi}_{l}\rangle_{\kappa,\partial\Omega}.$

With the definition of the orthonormal basis
$\\{\tilde{\xi}_{l}\\}_{l=1}^{\infty}$, the above series is actually finite.
This gives that for a.e. $z\in\Gamma$,

$\sum_{l=1}^{\infty}\tilde{\xi}_{l}(z)\langle\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}},\tilde{\xi}_{l}\rangle_{\kappa,\partial\Omega}=\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}(z),$

which leads to
$\sum_{\lambda_{n}=\lambda_{j}}\eta_{n}\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}(z)=0$ for a.e. $z\in\Gamma$. Assume that
$\\{\eta_{n}:\lambda_{n}=\lambda_{j}\\}\neq\\{0\\}$, then the linear
independence of the set $\\{\varphi_{n}(x):\lambda_{n}=\lambda_{j}\\}$ yields
that $\sum_{\lambda_{n}=\lambda_{j}}\eta_{n}\varphi_{n}(x)$ is not vanishing
on $\Omega$. However, we see that
$\sum_{\lambda_{n}=\lambda_{j}}\eta_{n}\varphi_{n}(x)$ is an eigenfunction
corresponding to the eigenvalue $\lambda_{j}$. These and Lemma 2.1 give that
$\sum_{\lambda_{n}=\lambda_{j}}\eta_{n}\frac{\partial\varphi_{n}}{\partial{\overrightarrow{\bf
n}}}$ can not vanish almost everywhere on $\Gamma$, which is a contradiction.
Hence, we prove that $\\{\eta_{n}:\lambda_{n}=\lambda_{j}\\}=\\{0\\}$ for each
$j\in\mathbb{N}^{+}$, namely, $\eta_{n}=0$ for $n\in\mathbb{N}^{+}$. The proof
is complete. ∎

######  Lemma 3.5.

Set $\Gamma\subset\partial\Omega$ be a nonempty open subset and let
$\sum_{n=1}^{\infty}c_{z,n}\eta_{n}$ be absolute convergent for a.e.
$z\in\Gamma$. Then given $\epsilon>0$, the result

$\sum_{n=1}^{\infty}c_{z,n}\eta_{n}(1-e^{-\lambda_{n}t})=0\ \text{for}\
t\in(0,\epsilon)\ \text{and\ a.e.}\ z\in\Gamma$

implies that $\eta_{n}=0$ for $n\in\mathbb{N}^{+}$.

###### Proof.

From [23, Lemma 3.5] we have that
$\sum_{\lambda_{n}=\lambda_{j}}c_{z,n}\eta_{n}=0$ for each
$j\in\mathbb{N}^{+}$ and a.e. $z\in\Gamma$. Then Lemma 3.4 gives the desired
result and completes the proof. ∎

######  Lemma 3.6.

Let the conditions in Lemma 3.5 be valid. For $s\in\mathbb{C}^{+}$, if
$\exists\epsilon>0$ such that

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

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

###### Proof.

For a.e. $z\in\Gamma,$ we define

$F_{z}(t):=\chi_{{}_{t\geq
0}}\sum_{n=1}^{\infty}c_{z,n}\eta_{n}(1-e^{-\lambda_{n}t}),\
t\in(-\infty,\infty).$

From the absolute convergence of the series
$\sum_{n=1}^{\infty}c_{z,n}\eta_{n}$, we see that

$|F_{z}(t)|\leq C\sum_{n=1}^{\infty}|c_{z,n}\eta_{n}|\leq C<\infty.$

Also, the direct computation yields that

$\displaystyle\int_{-\epsilon}^{\infty}|e^{-st}|\int_{0}^{t+\epsilon}\
|F_{z}(\tau)|\ d\tau\ dt$ $\displaystyle\leq
C\int_{-\epsilon}^{\infty}e^{-t\operatorname{Re}s}(t+\epsilon)\ dt$
$\displaystyle=Ce^{\epsilon\operatorname{Re}s}\int_{0}^{\infty}e^{-t\operatorname{Re}s}t\
dt$
$\displaystyle=Ce^{\epsilon\operatorname{Re}s}(\operatorname{Re}s)^{-2}<\infty.$

This implies the well-definedness of the integration
$\int_{-\epsilon}^{\infty}|e^{-st}|\int_{0}^{t+\epsilon}\ |F_{z}(\tau)|\
d\tau\ dt$ for $\operatorname{Re}s>0$ and a.e. $z\in\Gamma$.

Introducing the Heaviside function $H(t):=\chi_{{}_{t\geq 0}}$, we see that

$\int_{-\epsilon}^{\infty}e^{-st}\int_{0}^{t+\epsilon}\ F_{z}(\tau)\ d\tau\
dt=\int_{-\infty}^{\infty}e^{-st}\int_{-\infty}^{\infty}H(t-\tau+\epsilon)F_{z}(\tau)\
d\tau\ dt.$

The above argument provides the well-definedness of the integral

$\int_{-\infty}^{\infty}e^{-st}\int_{-\infty}^{\infty}H(t-\tau+\epsilon)F_{z}(\tau)\
d\tau\ dt$

when $\operatorname{Re}s>0$. With Fubini Theorem and the Dominated Convergence
Theorem, we have

$\displaystyle\int_{-\infty}^{\infty}e^{-ts}\int_{-\infty}^{\infty}H(t-\tau+\epsilon)F_{z}(\tau)\
d\tau\ dt$ $\displaystyle=\int_{-\epsilon}^{\infty}\\!e^{-ts}\,dt\
\int_{0}^{\infty}e^{-\tau
s}\Big{[}\sum_{n=1}^{\infty}c_{z,n}\eta_{n}(1-e^{-\lambda_{n}\tau})\Big{]}\,d\tau$
$\displaystyle=s^{-2}e^{\epsilon
s}\Big{[}\sum_{n=1}^{\infty}c_{z,n}\lambda_{n}\eta_{n}(s+\lambda_{n})^{-1}\Big{]},\quad\operatorname{Re}s>0,$

which leads to

$\lim_{\operatorname{Re}s\to\infty}\int_{-\infty}^{\infty}e^{-ts}\int_{-\infty}^{\infty}H(t-\tau+\epsilon)F_{z}(\tau)\
d\tau\ dt=0.$

Also, from the direct calculation, we have

$\displaystyle\int_{-\infty}^{\infty}e^{-ts}\int_{-\infty}^{\infty}H(t-\tau+\epsilon)F_{z}(\tau)\
d\tau\ dt$
$\displaystyle=\int_{-\epsilon}^{\infty}e^{-ts}\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t+\epsilon}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}\ dt$ $\displaystyle=\int_{0}^{\epsilon}e^{(\epsilon-t)s}\
\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}\ dt$ $\displaystyle\quad+\int_{0}^{\infty}e^{-ts}\
\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t+\epsilon}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}\ dt$ $\displaystyle=:S_{z,1}(s)+S_{z,2}(s).$

For $S_{z,2}(s)$, the Dominated Convergence Theorem and the absolute
convergence of $\sum_{n=1}^{\infty}c_{z,n}\eta_{n}$ give that

$|S_{z,2}(s)|\leq
C\sum_{n=1}^{\infty}|c_{z,n}\eta_{n}|\int_{0}^{\infty}e^{-t\operatorname{Re}s}|t+\epsilon|\
dt,$

and clearly the right side converges to $0$ as $\operatorname{Re}s\to\infty$.
So $\lim_{\operatorname{Re}s\to\infty}S_{z,2}(s)=0$. With the limit assumption
in this lemma, we have

$\lim_{\operatorname{Re}s\to\infty}S_{z,1}(s)=\lim_{\operatorname{Re}s\to\infty}[S_{z,1}(s)+S_{z,2}(s)]-\lim_{\operatorname{Re}s\to\infty}S_{z,2}(s)=0.$

For each $R_{0}\in\mathbb{R}$, when $\operatorname{Re}s\in(-\infty,R_{0})$, we
can see

$|S_{z,1}(s)|\leq C\int_{0}^{\epsilon}e^{(\epsilon-t)\operatorname{Re}s}t\
\sum_{n=1}^{\infty}|c_{z,n}\eta_{n}|\ dt\leq Ce^{\epsilon
R_{0}}\int_{0}^{\epsilon}t\ dt=C_{\epsilon,R_{0}}<\infty.$

This with the limit result $\lim_{\operatorname{Re}s\to\infty}S_{z,1}(s)=0$
yields that $S_{z,1}(s)$ is well-defined and bounded on the whole complex
plane $\mathbb{C}$. The definition of $S_{z,1}(s)$ gives that $S_{z,1}(s)$ is
holomorphic on $\mathbb{C}$. Hence, we have that $S_{z,1}(s)$ is a bounded
entire function for a.e. $z\in\Gamma$, which with Liouville’s theorem yields
that $S_{z,1}\equiv C$ on $\mathbb{C}$. Considering the limit result
$\lim_{\operatorname{Re}s\to\infty}S_{z,1}(s)=0$, it gives that

$S_{z,1}(s)=0\ \text{for}\ s\in\mathbb{C}\ \text{and a.e.}\ z\in\Gamma.$

Then we can see

$\displaystyle S_{z,1}(s)=$
$\displaystyle\int_{0}^{\epsilon}e^{(\epsilon-t)s}\
\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}\ dt$ $\displaystyle=$ $\displaystyle e^{\epsilon
s}\int_{0}^{\epsilon}e^{-ts}\
\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}\ dt\equiv 0,$

which means for $\operatorname{Re}s>0$,

$\displaystyle 0\equiv$ $\displaystyle\int_{0}^{\infty}e^{-st}H(\epsilon-t)\
\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}\,dt$ $\displaystyle=$
$\displaystyle\mathcal{L}\Big{\\{}H(\epsilon-t)\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{\\}}(s).$

Using [11, Corollary 8.1], it follows that

$\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau=0,\ t\in(0,\epsilon).$

The absolute convergence of $\sum_{n=1}^{\infty}c_{z,n}\eta_{n}$ supports the
uniform convergence of

$\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\ \text{and}\ \sum_{n=1}^{\infty}c_{z,n}\eta_{n}(1-e^{-\lambda_{n}t})$

on $(0,\epsilon)$. Thus, from [21, Theorem 7.17], differentiating the series
$\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau$ on $t$ yields that for $t\in(0,\epsilon)$ and a.e. $z\in\Gamma$,

$0=\frac{d}{dt}\Big{[}\sum_{n=1}^{\infty}c_{z,n}\int_{0}^{t}\eta_{n}(1-e^{-\lambda_{n}\tau})\
d\tau\Big{]}=\sum_{n=1}^{\infty}c_{z,n}\eta_{n}(1-e^{-\lambda_{n}t}).$

This together with Lemma 3.5 leads to $\eta_{n}=0$ for $n\in\mathbb{N}^{+}$,
and completes the proof. ∎

###  3.4 Proof of Theorem 1.

Now, it is time to show Theorem 1. For shorten the proof, the following
notations will be used:

$\displaystyle p_{k,n}=\langle
p_{k}(\cdot),\varphi_{n}(\cdot)\rangle_{\Omega},$ $\displaystyle
Q_{z,k}(s)=\sum_{n=1}^{\infty}c_{z,n}p_{k,n}\lambda_{n}(s+\lambda_{n})^{-1},$
$\displaystyle\tilde{p}_{k,n}=\langle\tilde{p}_{k}(\cdot),\varphi_{n}(\cdot)\rangle_{\Omega},$
$\displaystyle\tilde{Q}_{z,k}(s)=\sum_{n=1}^{\infty}c_{z,n}\tilde{p}_{k,n}\lambda_{n}(s+\lambda_{n})^{-1}.$

###### Proof of Theorem 1.

The result (3.2), Lemma 3.3 and the analytic continuation give that for
$s\in\mathbb{C}^{+}$ and a.e. $z\in\Gamma$,

$\sum_{k=1}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})Q_{z,k}(s)=\sum_{k=1}^{\tilde{K}}(e^{-\tilde{t}_{k-1}s}-e^{-\tilde{t}_{k}s})\tilde{Q}_{z,k}(s).$
(3.3)

Firstly let us show $t_{0}=\tilde{t}_{0}$. Without loss of generality, we can
assume that $t_{0}<\tilde{t}_{0}$. Picking $\epsilon>0$ such that
$\epsilon<\min\\{\tilde{t}_{0}-t_{0},t_{1}-t_{0}\\}$, from (3.3), we have for
$s\in\mathbb{C}^{+}$ and a.e. $z\in\Gamma$,

$\displaystyle e^{\epsilon s}Q_{z,1}(s)=$ $\displaystyle
e^{(\epsilon+t_{0}-t_{1})s}Q_{z,1}(s)-\sum_{k=2}^{K}(e^{(\epsilon+t_{0}-t_{k-1})s}-e^{(\epsilon+t_{0}-t_{k})s})Q_{z,k}(s)$
(3.4)
$\displaystyle+\sum_{k=1}^{\tilde{K}}(e^{(\epsilon+t_{0}-\tilde{t}_{k-1})s}-e^{(\epsilon+t_{0}-\tilde{t}_{k})s})\tilde{Q}_{z,k}(s).$

We see that $\lambda_{n}(s+\lambda_{n})^{-1}$ is uniformly bounded on
$\mathbb{C}^{+}$. Recalling the convergence result
$\sum_{k=1}^{K}\sum_{n=1}^{\infty}|c_{z,n}p_{k,n}|<\infty$ from Assumption
2.4, the uniform convergence of the series
$\sum_{k=2}^{K}(e^{(\epsilon+t_{0}-t_{k-1})s}-e^{(\epsilon+t_{0}-t_{k})s})Q_{z,k}(s)$
on the half plane $\mathbb{C}^{+}$ can be derived. Then we have

$\displaystyle\lim_{\operatorname{Re}s\to\infty}\sum_{k=2}^{K}(e^{(\epsilon+t_{0}-t_{k-1})s}-e^{(\epsilon+t_{0}-t_{k})s})Q_{z,k}(s)$
$\displaystyle=\sum_{k=2}^{K}\lim_{\operatorname{Re}s\to\infty}(e^{(\epsilon+t_{0}-t_{k-1})s}-e^{(\epsilon+t_{0}-t_{k})s})Q_{z,k}(s)$
$\displaystyle=0.$

Similarly, we see that other terms in the right side of (3.4) also tend to
zero as $\operatorname{Re}s\to\infty$. So we have
$\lim_{\operatorname{Re}s\to\infty}e^{\epsilon s}Q_{z,1}(s)=0$ for a.e.
$z\in\Gamma.$ This with Lemma 3.6 gives that $p_{1,n}=0$ for
$n\in\mathbb{N}^{+}$, i.e. $\|p_{1}\|_{L^{2}(\Omega)}=0$, which contradicts
with Assumption 1.1. Hence, we have $t_{0}=\tilde{t}_{0}$.

Inserting $t_{0}=\tilde{t}_{0}$ into (3.3) and picking
$0<\epsilon<\min\\{t_{1}-t_{0},\tilde{t}_{1}-t_{0}\\},$ we have

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

From the previous arguments, we can prove

$\lim_{\operatorname{Re}s\to\infty}e^{\epsilon
s}[Q_{z,1}(s)-\tilde{Q}_{z,1}(s)]=0,\ \text{a.e.}\ z\in\Gamma.$

Lemma 3.6 gives that $p_{1,n}=\tilde{p}_{1,n}$ for $n\in\mathbb{N}^{+}$,
namely $\|p_{1}-\tilde{p}_{1}\|_{L^{2}(\Omega)}=0$.

Next we want to show $t_{1}=\tilde{t}_{1}$. Owing the results
$p_{1,n}=\tilde{p}_{1,n}$ and $t_{0}=\tilde{t}_{0}$ into (3.3), it follows
that for $s\in\mathbb{C}^{+}$ and a.e. $z\in\Gamma$,

$\displaystyle e^{-t_{1}s}[Q_{z,2}(s)-Q_{z,1}(s)]=$ $\displaystyle
e^{-\tilde{t}_{1}s}[\tilde{Q}_{z,2}(s)-\tilde{Q}_{z,1}(s)]+e^{-t_{2}s}Q_{z,2}(s)-e^{-\tilde{t}_{2}s}\tilde{Q}_{z,2}(s)$
$\displaystyle-\sum_{k=3}^{K}(e^{-t_{k-1}s}-e^{-t_{k}s})Q_{z,k}(s)+\sum_{k=3}^{\tilde{K}}(e^{-\tilde{t}_{k-1}s}-e^{-\tilde{t}_{k}s})\tilde{Q}_{z,k}(s).$

Without loss of generality we assume that $t_{1}<\tilde{t}_{1}$ and pick
$0<\epsilon<\min\\{\tilde{t}_{1}-t_{1},t_{2}-t_{1},\tilde{t}_{2}-t_{1}\\}$.
Following the proof for $t_{0}=\tilde{t}_{0}$, we can show that
$\|p_{2}-p_{1}\|_{L^{2}(\Omega)}=0$. This is a contradiction. Hence
$t_{1}=\tilde{t}_{1}$.

Now with the results $t_{0}=\tilde{t}_{0}$, $t_{1}=\tilde{t}_{1}$ and
$\|p_{1}-\tilde{p}_{1}\|_{L^{2}(\Omega)}=0$, from (3.3), we have that for
$\operatorname{Re}s>0$ and a.e. $z\in\Gamma$,

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

From the previous proof, we can derive that $t_{2}=\tilde{t}_{2}$ and
$\|p_{2}-\tilde{p}_{2}\|_{L^{2}(\Omega)}=0$. Using this technique recursively,
we can conclude that

$t_{0}=\tilde{t}_{0},\ t_{k}=\tilde{t}_{k},\
\|p_{k}-\tilde{p}_{k}\|_{L^{2}(\Omega)}=0\ \text{for}\
k=1,\cdots,\min\\{K,\tilde{K}\\}.$ (3.5)

Finally we need to show $K=\tilde{K}$ (For the case of $K=\infty$ and
$\tilde{K}=\infty$, we regard them as $K=\tilde{K}$). Without loss of
generality assuming that $K<\tilde{K}$, from (3.3) and (3.5), we have for
$\operatorname{Re}s>0$ and a.e. $z\in\Gamma$,

$\sum_{k=K+1}^{\tilde{K}}(e^{-\tilde{t}_{k-1}s}-e^{-\tilde{t}_{k}s})\tilde{Q}_{z,k}(s)=0.$

The former proof gives that $\|\tilde{p}_{K+1}\|_{L^{2}(\Omega)}=0$, which is
a contradiction. Hence, we have that $K=\tilde{K}$, which together with (3.5)
leads to the desired result. The proof is complete. ∎

##  4 Numerical reconstructions.

In this section, we are going to conduct several numerical experiments which
realize the inversion process described in Theorem 1. Here we take the example
of locating the leaking oil tanker in the ocean. To be more specific, suppose
an oil tanker is planning to travel from harbor A to harbor B and the planned
route is known. The ship is leaking oil during the trip and the exact
trajectory of the ship is unknown. The target of the engineers is to find the
exact location where the ship leaks the oil during the trip. The physics of
the oil diffusion in the ocean can be modelled as a parabolic equation (1.1)
and the leaking oil is the heat source in the PDE. By Theorem 1, one can find
these positions given the fluxes from a nonempty open subset of the boundary.

The oil tanker problem can be formulated as estimating the state of the system
(position of the ship) when a sequence of observations in time becomes
available. The popular choice to solve the sequential prediction problem is
the sequential Monte Carlo (SMC) method, which is also called the particle
filter (PF). In this work, we are going to use one variant of the SMC as the
approach to solving the inverse problem. Before presenting the results of the
numerical experiments, we will first review some preliminaries of the SMC
algorithm.

###  4.1 Sequential Monte Carlo (SMC).

Consider a discrete-time Markov process $\\{x_{t}\\}_{t=0}^{n}$:

$\displaystyle x_{0}\sim p(x_{0})\text{ and }x_{n}|x_{n-1}\sim
p(x_{n}|x_{n-1}),$ (4.1)

where $\sim$ means distributed as, $p(x_{0})$ is the prior of the initial
state $x_{0}$ and $p(x^{\prime}|x)$ is the given transition probability prior
from current state $x$ to the next state $x^{\prime}$. We are interested in
estimating $\\{x_{t}\\}_{t=0}^{n}$ which is also denoted as $x_{0:n}$, given
the observations $\\{y_{t}\\}_{t=0}^{n}$ at each time step. In the oil tanker
application, the state $x_{t}$ is the position of the tanker and $y_{t}$ is
the flux on the boundary.

The equation (4.1) together with the likelihood $p(y_{t}|x_{t})$ will define a
Bayesian model. We have the prior on the trajectory
$x_{0:n}:=\\{x_{0},...,x_{n}\\}$ as follow:

$\displaystyle p(x_{0:n})=p(x_{0})\Pi_{k=1}^{n}p(x_{k}|x_{k-1}).$ (4.2)

Moreover,

$\displaystyle p(y_{0:n}|x_{0:n})=\Pi_{k=1}^{n}p(y_{k}|x_{k}).$ (4.3)

From a Bayesian perspective, the posterior distribution of $x_{0:n}$ given a
sequence of the observations $y_{0:n}$ is:

$\displaystyle p(x_{0:n}|y_{0:n})=\frac{p(x_{0:n},y_{0:n})}{p(y_{0:n})},$
(4.4)

where $p(x_{0:n},y_{0:n})=p(y_{0:n}|x_{0:n})p(x_{0:n})$ and $p(y_{0:n})=\int
p(x_{0:n},y_{0:n})dx_{0:n}$.

In our application, the most time-consuming step is to evaluate the likelihood
function $p(y_{0:n}|x_{0:n})$, and we will follow the classical framework as
in [24]. Denoting the forward solver which maps the input state $u$ to the
observation $y$ as $\mathcal{G}$, it follows that

$\displaystyle y=\mathcal{G}(u)+\eta,$

where $\eta$ is the error associated with the process. There are various
sources of error. For example, the error in the observation and the error of
the forward solver. The total error $\eta$ is then assumed to follow the
Gaussian distribution $\mathcal{N}(0,B)$ where $B$ is the covariance with a
proper size. The likelihood function then has the form:

$\displaystyle p(y|u):=p(y-\mathcal{G}(u))=p(\eta)\sim\mathcal{N}(0,B).$

Now the area of interesting is the tracking problem: find the current state
given the observations. Theoretically, this means that one needs to find a
group of $p(x_{0:n}|y_{0:n})$. By the prior (4.2) and the likelihood (4.3),
the joint distribution $p(x_{0:n},y_{0:n})$ in (4.4) then satisfies

$\displaystyle
p(x_{0:n},y_{0:n})=p(x_{1:n},y_{1:n})p(y_{n}|x_{n})p(x_{n}|x_{n-1}),$

and it follows that

$\displaystyle
p(x_{0:n}|y_{0:n})=p(x_{0:n-1}|y_{0:n-1})\frac{p(x_{n}|x_{n-1})p(y_{n}|x_{n})}{p(y_{n}|y_{n-1})},$
(4.5)

where $p(y_{n}|y_{n-1})=\int
p(x_{n-1}|y_{0:n-1})p(x_{n}|x_{n-1})p(y_{n}|x_{n})dx_{{n-1:n}}$.

Most particle filtering methods are created by a numerical approximation to
(4.5). A common and powerful algorithm is the sequential Monte Carlo (SMC)
approximation. One can use a set of samples (particles) which are drawn from a
posterior distribution to approximate the target posterior. More precisely,
the posterior is the average calculated by

$\displaystyle\hat{p}(x_{0:n}|y_{0:n})=\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}\delta_{x_{0:t}^{i}}dx_{0:t}.$

Here, $\\{x_{0:t}^{i}\\}_{i=1}^{N_{p}}$ are the particles and $\delta(d\cdot)$
is the standard delta function. As a result, when the number of samples
$N_{p}$ is large enough, the average will approximate the true distribution.
One of the advantages of the method is that one can prove the asymptotic
convergence to the target distribution of interest. Besides, compared to some
other approaches such as the extended Kalman filter and unscented Kalman
filter, the SMC approach does not assume that the state models are Gaussian
[9, 25]. One hence can apply the method in a broader area.

Unfortunately, it is not easy to draw samples from the posterior. An
alternative way of drawing is to sample from a proposal distribution
$q(x_{0:t}|y_{0:t})$ and to approximate the posterior using a weighted sum of
particles. One usually calls this methodology as the Bayesian important
sampling (IS).

Now, two questions remain unsolved. One is how to design a proposal
distribution; another one is how to calculate the weight for each sample. How
to select $q(x_{0:t}|y_{0:t})$ is an interesting research topic. One usually
requires that the $q(x_{0:t}|y_{0:t})$ satisfies the following structure:

$\displaystyle
q(x_{0:t}|y_{0:t})=q(x_{0:t-1}|y_{1:t-1})q(x_{t}|x_{0:t-1},y_{1:t}).$

In this work, we are going to set $q(x_{0}|y_{0})=p(x_{0})$ and
$q(x_{t}|x_{0:t-1},y_{1:t})=p(x_{t}|x_{t-1})$. According to [12, 1], if the
latent variable dimension is not large and the observations are not too
informative, setting the proposal to be equal to transition probability is
good enough. To verify our method, we will conduct experiments with the
transition probabilities $p(x_{t}|x_{t-1})$; the detailed design is shown in
Section 4.4 and Section 4.5. The weights can then be calculated as follows:

$\displaystyle
w_{t}=w_{t-1}\frac{p(y_{t}|x_{t})p(x_{t}|x_{t-1})}{q(x_{t}|x_{0:t-1},y_{1:t})},$

and $w_{1}=\frac{p(x_{1})p(y_{1}|x_{1})}{q(x_{1}|y_{1})}$, where
$p(y_{t}|x_{t})$ is the likelihood. We omit the derivation and one can refer
to [9, 25, 1] for the details.

One issue of the algorithm is the degeneration of the weights, that is, the
variance of the importance weights will become larger and larger with respect
to time. Consequently, there are only a few samples which have a meaningful
weight when $t$ is large. One way to reduce the effect of degeneration is to
use the resampling method. As a result, the algorithm is called the sequential
importance resampling (SIR). We provide an SIR algorithm in the Appendix A.
There are various of resampling methods; in this work, we apply the systematic
resampling as suggested in [9]; please check Appendix A for the details.

###  4.2 Sources.

In the next two subsections, we are going to demonstrate how we set up the
testing problems. For equation (1.1), we set $\Omega=[0,1]^{2}$ and the
terminal simulation time is $T=0.1$. Also $\kappa(x)$ is the given
permeability which will be defined later in Section 4.3, and the source
$\sum_{k=1}^{K}p_{k}(x)\chi_{{}_{t\in[t_{k-1},t_{k})}}$ is partially known.
More specifically, the time mesh grid $\\{t_{k}\\}_{k=0}^{K}$ is given and the
spatial components $\\{p_{k}(x)\\}_{k=1}^{K}$ are set to be characteristic
functions. The k-th source $p_{k}(x)$ is defined as:

$p_{k}(x)=\begin{cases}p_{k},\quad x\in A_{k},\\\ 0,\quad\
\text{otherwise},\end{cases}$ (4.6)

Here $p_{k}$ is a known constant; $A_{k}$ is a square with a known size (is
equal to $0.06\times 0.06$ in all our experiments), but the location is not
given. Our target work is to find one vertex coordinate of the $A_{k}$ and
$k=1,...,K$. One possible $p_{k}(x)$ is shown in Figure 1 (left). In the
following experiments, we will consider 2 sources. In both cases, $K=10$ and
$p_{k}=1250$ for all $k$.

We here list the vertices (left top) coordinates for both sources and we will
trace these two sequences of the coordinates later by our proposed method. The
first list of vertices is: (0.12, 0.12), (0.20, 0.20), (0.28, 0.28), (0.36,
0.36), (0.44, 0.44), (0.52, 0.52), (0.60, 0.60), (0.68, 0.68), (0.76, 0.76),
(0.84, 0.84); the second list of vertices is: (0.12, 0.12), (0.20, 0.24),
(0.28, 0.36), (0.36, 0.48), (0.44, 0.56), (0.52, 0.64), (0.60, 0.72), (0.68,
0.78), (0.76, 0.84), (0.84, 0.90). For better illustration, we plot the
trajectories of the two sources in Figure 1 and we will discuss these two
cases in Section 4.4 and Section 4.5 respectively. In the oil tanker
application, the ship leaks oil in $A_{k}$ for $k=1,2...,K$ and we need to
trace the position of all $A_{k}$ (the left top vertices) given some prior
information about how this ship is moving (velocity of the ship).

Figure 1: Demonstration of the source. Left: example of $p_{k}(x)$ which is
the characteristic function defined in equation (4.6). Middle and right images
show the trajectories of the source $\\{p_{k}(x)\\}_{k=1}^{K}$. Middle: source
one, which is used in Section 4.4. Right: source two, which is used in Section
4.5.

###  4.3 Permeability fields.

In equation (1.1), the operator $\mathcal{A}$ includes the permeability field
$\kappa(x)$. In order to better demonstrate the idea of the theorem and the
algorithm, we will solve several multiscale problems. Two different
permeability fields are used here. The first $\kappa_{1}(x)$ has multiple
frequencies and is defined as follow:

$\displaystyle\kappa_{1}(x,y)=15\sin(2\pi\cdot 0.01x)\cdot\sin(2\pi\cdot
0.05y)+0.05\sin(2\pi\cdot 6.5x)\cdot\sin(2\pi\cdot 6y)+0.1,$

where $(x,y)\in\Omega:=[0,1]^{2}$ and it is demonstrated in Figure 2.

Figure 2: The first permeability field 2D demonstration.
$\max\kappa_{1}(x)=0.3978$ and $\min\kappa_{1}(x)=0.0521$.

The second permeability field $\kappa_{2}(x)$ (see Figure 3 for the
illustration) has multiple high contrast channels and is widely used in the
multiscale finite element method society [6, 4].

Figure 3: The second permeability field 2D demonstration. The permeability in
the yellow channels and dots is equal to $10^{4}$, while the permeability in
the purple background is equal to $1$.

The multiple frequencies in $\kappa_{1}(x)$ and high contrast in
$\kappa_{2}(x)$ bring difficulty in solving the corresponding equations
numerically. In this work, we will use the standard finite element solver with
spatial discretization $\Delta x=10^{-2}$. For the temporal approximation, the
backward Euler scheme is adopted, and the time discretization is $\Delta
t=10^{-3}$.

###  4.4 The first set of experiments.

In this part, we present the first set of experiments whose source is defined
in the middle image of Figure 1. As we have discussed before, the proposal
distribution $q(x^{\prime}|y^{\prime},x)$ is the same with the transition
probability $p(x^{\prime}|x)$ where $x^{\prime},y^{\prime}$ denote the new
sample and the observation of the next time step, respectively. We hence only
need to specify $p(x^{\prime}|x)$, which is defined in Algorithm 1.

1 Input: current $x=(x_{1},x_{2})$;

2 Initialize: Set the size related parameters $r_{1},r_{2}>0$ and speed
related parameter $s_{0}$, and draw a random number
$s\sim\mathcal{U}(0,s_{t})$. We then define $l_{1}=s_{0}+s+r_{1}$ and
$l_{2}=s_{0}+s+r_{2}$;

3 Create a rectangle: the top right vertex coordinate is $(x_{1}+l_{1}\Delta
x,x_{2}+l_{2}\Delta x)$ and the size is equal to $(2r_{1}\Delta x,2r_{2}\Delta
x)$;

Draw a point $x^{\prime}$ from the rectangle above uniformly and this point
will be the next sample.

Algorithm 1 Transition algorithm

Here we set $s_{0}=6,\ s_{t}=5$ and $r_{1}=r_{2}=4$. It should be noted that
if one is given the exact current state $x=(x_{1},x_{2})$, the next exact
state is $x^{\prime}=(x_{1}+8\Delta x,x_{2}+8\Delta x)$ (please check the
trajectory of the first source in Section 4.2); however, the center (mean) of
all proposed regions is deviated from the $x^{\prime}$. This simulates the
real life scenario: one only knows the ship’s route plan; however the ship may
depart from the original route, and one needs to find the real route and
positions when the ship leaks.

The observation is the flux of the boundary. For the multiple frequencies
example, we measure the flux in the interval $[0.45,0.52]$ from each boundary;
while for the high contrast permeability example, we only measure the flux in
the interval $[0.5,0.55]$. In both cases, the number of particles is $320$. To
evaluate our results, we calculate the mean of the samples at each time step
$k$ (this indicates the time step in the source) and compare it with the pre-
set values. The demonstration and the relative error are shown in Figure 4.

Figure 4: Samples mean and the relative error of the multiple frequencies
permeability with the first source. Left: samples mean and true source at each
time step $k$. Right: relative error of the samples mean.

Figure 5: Samples mean and the relative error of the high contrast
permeability with the first source. Left: samples mean and true source at each
time step $k$. Right: relative error of the samples mean.

###  4.5 The second set of experiments.

In this set of experiments, we will use the second source which is
demonstrated as the right image in Figure 1. We use the same transition
algorithm in Algorithm 1 but we use a different set of parameters,
specifically, $s_{0}=4,\ s_{t}=5,\ r_{1}=4,\ r_{2}=5$. Similar to the first
source, the center (mean) of all proposed regions is deviated from
$x^{\prime}$. This is a real situation in which the ship’s actual route is
different from the original plan. The setting makes the predictions more
challenging. For the flux measurement, we choose an interval $[0.40,0.60]$
from each boundary for the multiple frequencies permeability and $[0.43,0.58]$
from each boundary for the high contract permeability. The results are shown
in Figures 6 and 7, respectively.

Figure 6: Samples mean and the relative error of the multiple frequencies
permeability with the second source. Left: samples mean and true source at
each time step $k$. Right: relative error of the samples mean.

Figure 7: Samples mean and the relative error of the high contrast
permeability with the second source. Left: samples mean and true source at
each time step $k$. Right: relative error of the samples mean.

From the four examples presented above, we conclude that our algorithm is able
to trace the exact source positions even if the proposal is far away from the
exact route. The next challenge for us is to reduce the amount of observation
data since by the theory, we can trace the source with flux on any open subset
of the boundary. We will study this topic in the future.

##  5 Concluding remarks and future works.

In this work, we consider the inverse source problem in the parabolic
equation. The unknown source has a semi-discrete formulation, which can be
used to approximate the general form $F(x,t)$. We prove the uniqueness
theorem–Theorem 1, which says the data from any nonempty open subset of the
boundary can support the uniqueness of the source. This conclusion is of
significance in practical applications since it indicates that the source can
be recovered from sparse boundary data and then save the cost. For the
theoretical analysis, there is an interesting and meaningful work in the
future, which is to determine the minimal observed area. Theorem 1 illustrates
that the nonempty open subset of the boundary is sufficient to support the
uniqueness. However, only from this theorem, we do not know whether we can
minimize the observed area further. To find the minimal observed area is one
of our future works.

In the aspect of numerical reconstructions, we propose an oil diffusion
problem which is modelled by the parabolic equation. One needs to trace the
positions of the oil tanker given the flux of the boundary. Since the model
problem is a sequential prediction with a sequence of observations up to the
time $n$, it is natural to formulate the problem to the Bayesian filtering
problem. We apply one particle filter algorithm to solve the inverse problem
with two multiscale permeabilities. Evaluating the likelihood function is the
most time consuming step since it requires solving a forward problem with a
fine resolution solver. We use the finite element solver in this work,
however, the solver can be improved by the multiscale finite element solver or
the deep learning solver. These solvers are more efficient meanwhile
preserving accuracy and we will study these topics in the future. We can also
apply deep learning to train a solver. This is equivalent to solving the
stochastic parametric PDE with the deep neural network. More precisely, we can
train a mapping from the source to the flux on the boundary. One benefit of
the method is: it will increase the computation efficiency in evaluating the
likelihood function. Another benefit is: we can linearize the network and
hence more advanced particle filter algorithms can be used.

## References

  * [1] Christophe Andrieu, Arnaud Doucet, and Roman Holenstein.  Particle markov chain monte carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342, 2010.
  * [2] J. R. Cannon.  Determination of an unknown heat source from overspecified boundary data source in the heat equation.  SIAM Journal on Numerical Analysis, 5(2), 1968.
  * [3] Jin Cheng and Jijun Liu.  An inverse source problem for parabolic equations with local measurements.  Appl. Math. Lett., 103:106213, 7, 2020.  doi:10.1016/j.aml.2020.106213.
  * [4] Boris Chetverushkin, Eric Chung, Yalchin Efendiev, Sai-Mang Pun, and Zecheng Zhang.  Computational multiscale methods for quasi-gas dynamic equations.  Journal of Computational Physics, 440:110352, 2021.
  * [5] Eric Chung, Yalchin Efendiev, Wing Tat Leung, Sai-Mang Pun, and Zecheng Zhang.  Multi-agent reinforcement learning accelerated mcmc on multiscale inversion problem.  arXiv preprint arXiv:2011.08954, 2020.
  * [6] Eric Chung, Yalchin Efendiev, Sai-Mang Pun, and Zecheng Zhang.  Computational multiscale methods for parabolic wave approximations in heterogeneous media.  arXiv preprint arXiv:2104.02283, 2021.
  * [7] Randal Douc and Olivier Cappé.  Comparison of resampling schemes for particle filtering.  In ISPA 2005. Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis, 2005., pages 64–69. IEEE, 2005.
  * [8] Arnaud Doucet, Simon Godsill, and Christophe Andrieu.  On sequential monte carlo sampling methods for bayesian filtering.  Statistics and computing, 10(3):197–208, 2000.
  * [9] Arnaud Doucet and Adam M Johansen.  A tutorial on particle filtering and smoothing: Fifteen years later.  Handbook of nonlinear filtering, 12(656-704):3, 2009.
  * [10] Yalchin Efendiev, Thomas Hou, and Wuan Luo.  Preconditioning markov chain monte carlo simulations using coarse-scale models.  SIAM Journal on Scientific Computing, 28(2):776–803, 2006.
  * [11] Gerald B. Folland.  Fourier analysis and its applications.  The Wadsworth & Brooks/Cole Mathematics Series. Wadsworth & Brooks/Cole Advanced Books & Software, Pacific Grove, CA, 1992.
  * [12] Neil J Gordon, David J Salmond, and Adrian FM Smith.  Novel approach to nonlinear/non-gaussian bayesian state estimation.  In IEE proceedings F (radar and signal processing), volume 140, pages 107–113. IET, 1993.
  * [13] F. Hettlich and W. Rundell.  Identification of a discontinuous source in the heat equation.  Inverse Problems, 17(5):1465–1482, 2001.  URL: http://dx.doi.org/10.1088/0266-5611/17/5/315, doi:10.1088/0266-5611/17/5/315.
  * [14] X. Huang, O. Yu. Imanuvilov, and M. Yamamoto.  Stability for inverse source problems by Carleman estimates.  Inverse Problems, 36(12):125006, 20, 2020.  doi:10.1088/1361-6420/aba892.
  * [15] Masaru Ikehata.  An inverse source problem for the heat equation and the enclosure method.  Inverse Problems, 23(1):183–202, 2007.  doi:10.1088/0266-5611/23/1/009.
  * [16] Victor Isakov.  Inverse source problems, volume 34 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 1990.  doi:10.1090/surv/034.
  * [17] Bangti Jin, Yavar Kian, and Zhi Zhou.  Reconstruction of a space-time-dependent source in subdiffusion models via a perturbation approach.  SIAM J. Math. Anal., 53(4):4445–4473, 2021.  doi:10.1137/21M1397295.
  * [18] Marine Jouin, Rafael Gouriveau, Daniel Hissel, Marie-Cécile Péra, and Noureddine Zerhouni.  Particle filter-based prognostics: Review, discussion and perspectives.  Mechanical Systems and Signal Processing, 72:2–31, 2016.
  * [19] Zhiyuan Li and Zhidong Zhang.  Unique determination of fractional order and source term in a fractional diffusion equation from sparse boundary data.  Inverse Problems, 36(11):115013, 20, 2020.  doi:10.1088/1361-6420/abbc5d.
  * [20] Guang Lin, Yating Wang, and Zecheng Zhang.  Multi-variance replica exchange stochastic gradient mcmc for inverse and forward bayesian physics-informed neural network.  arXiv preprint arXiv:2107.06330, 2021.
  * [21] Walter Rudin.  Principles of mathematical analysis.  McGraw-Hill Book Co., New York-Auckland-Düsseldorf, third edition, 1976.  International Series in Pure and Applied Mathematics.
  * [22] William Rundell and Zhidong Zhang.  Recovering an unknown source in a fractional diffusion problem.  J. Comput. Phys., 368:299–314, 2018.  doi:10.1016/j.jcp.2018.04.046.
  * [23] 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.  doi:10.1137/19M1279915.
  * [24] Andrew M Stuart.  Inverse problems: a bayesian perspective.  Acta numerica, 19:451–559, 2010.
  * [25] Rudolph Van Der Merwe, Arnaud Doucet, Nando De Freitas, and Eric A Wan.  The unscented particle filter.  2000\.
  * [26] Masahiro Yamamoto.  Carleman estimates for parabolic equations and applications.  Inverse Problems, 25(12):123013, 75, 2009.  doi:10.1088/0266-5611/25/12/123013.

##  Appendix A Particle filter algorithm.

Here we give the particle filter algorithm (PF) which we used in the
simulation.

1 Initialize the number of time steps $n$ and the number of particles $N$;

2 Step 0: at $t=0$, draw the states $x_{0}^{i}$ for $i=1,...,N$ from the prior
$p(x_{0})$;

3 for  _$t=1,2,...,n$_ do

4 for  _$i=1,...,N$_ do

5 Sample $\hat{x}^{i}_{t}$ from the prior $q(x_{t}|x_{0:t-1}^{i},y_{1:t})$ and
extend the current trajectory by adding the temporary proposed state
$\hat{x}_{0:t}^{i}=(x_{0:t-1}^{i},\hat{x}_{t}^{i})$;

6 Calculate the importance weights recursively and normalize the resulting
weights as
$w_{t}^{i}=w_{t-1}^{i}\frac{p(y_{t}|\hat{x}_{t}^{i})p(x_{t}^{i}|x_{t-1}^{i})}{q(\hat{x}_{t}^{i}|x_{0:t-1}^{i},y_{1:t})}.\;$

7 end for

8 for  _$i=1,...,N$_ do

9 Normalize the importance weights for the resampling purpose as
$\tilde{w}_{t}^{i}=\frac{w_{t}^{i}}{\sum_{j=1}^{N}w_{t}^{j}}.\;$

10 end for

11 Resampling: multiply samples $\hat{x}_{0:t}^{i}$ with the normalized
resampling weight $\tilde{w}_{t}^{i}$ to obtain $N$ random samples
$x_{0:t}^{i}$ which is roughly distributed following $p(x_{0:t}^{i}|y_{1:t})$;

12 Set $w_{t}^{i}=\tilde{w}_{t}^{i}=\frac{1}{N}$ for all $i=1,...,N$;

13 The algorithm will finally return a set of samples whose average is an
approximation of the target distribution:
$p(x_{0:t}|y_{0:t})\approx\hat{p}(x_{0:t}|y_{0:t})=\frac{1}{N}\sum_{i=1}^{N}\delta_{(x_{0:t}^{i})}(dx_{0:t}).\;$

14 end for

Algorithm 2 Sequential importance resampling (SIR)

It should be noted that we use the systematic resampling method in this work.
To implement this method, one first samples
$U_{1}\sim\mathcal{U}[0,\frac{1}{N}]$ and defines $U_{i}=U_{1}+\frac{i-1}{N}$
for $i=2,...,N$, then we set

$N_{n}^{i}=\big{|}\\{U_{j}:\sum_{k=1}^{i-1}W_{n}^{k}\leq
U_{j}\leq\sum_{k=1}^{i}W_{n}^{k}\\}\big{|}$

with the convention $\sum_{k=1}^{0}:=0$.