# Path integral molecular dynamics for thermodynamics and Green’s function of
ultracold spinor bosons

Yongle Yu  Innovation Academy for Precision Measurement Science and
Technology, Chinese Academy of Sciences, Wuhan 430071, P. R. China  Shujuan
Liu  College of Science, Zhejiang University of Technology, Hangzhou 31023, P.
R. China  Hongwei Xiong  College of Science, Zhejiang University of
Technology, Hangzhou 31023, P. R. China  Yunuo Xiong<EMAIL_ADDRESS>College of Science, Zhejiang University of Technology, Hangzhou 31023, P. R.
China

###### Abstract

Most recently, the path integral molecular dynamics has been successfully used
to consider the thermodynamics of single-component identical bosons and
fermions. In this work, the path integral molecular dynamics is developed to
simulate the thermodynamics, Green’s function and momentum distribution of
two-component bosons in three dimensions. As an example of our general method,
we consider the thermodynamics of up to sixteen bosons in a three-dimensional
harmonic trap. For noninteracting spinor bosons, our simulation shows a bump
in the heat capacity. As the repulsive interaction strength increases,
however, we find the gradual disappearance of the bump in the heat capacity.
We believe this simulation result can be tested by ultracold spinor bosons
with optical lattices and magnetic-field Feshbach resonance to tune the inter-
particle interaction. We also calculate Green’s function and momentum
distribution of spinor bosons. Our work facilitates the exact numerical
simulation of spinor bosons, whose property is one of the major problems in
ultracold Bose gases.

##  I Introduction

The exact numerical simulation of identical particles at finite temperature or
zero temperature is one of the major problems in quantum physics. The so-
called exact numerical simulation means that we should not assume in advance
the trial wave function Rossi  or assume the usual contact interaction Parisi
; Parisi1  for the convenience of calculations. The exact numerical simulation
of the thermodynamics of identical particles requests that we carry out the
theoretical studies directly from the partition function of the quantum
system. Fortunately, without predetermining any properties of the quantum
system, the partition function in the path integral formalism feynman ;
kleinert ; Tuckerman  can be mapped to an equivalent classical ensemble of
interacting ring polymers chandler ; Parrinello ; Miura ; Cao ; Cao2 ; Jang2 ;
Ram ; Kinugawa ; Roy ; Roy1 ; Roy2 ; Roy3 ; Kinugawa1 ; Kinugawa2 ; Roy4 ;
Roy5 ; Roy6 ; Poly ; Craig ; Braa ; Haber ; Thomas . However, during this
mapping, one of the key obstacles we face is how to give exact consideration
of quantum exchange effects for identical particles. In a pioneering work
Hirshberg  on single-component identical bosons, a recursion formula is found
to express the partition function in a highly efficient way. In Ref. Xiong ,
path integral molecular dynamics (PIMD) and the recursion formula are
generalized to calculate Green’s function and momentum distribution of single-
component bosons. Other recent advances include the application of PIMD to
supersolid phase in high-pressure deuterium Deuterium  and bosons in singular
magnetic field Xiong3 , and the generalization of PIMD to spin-polarized
fermions HirshbergFermi ; Xiong2 ; Xiong4 . In particular, PIMD was used most
recently to consider the fermion sign problem for large fermion system Xiong4
.

Of course, we need to generalize PIMD to two-component (or spinor) bosons,
because this is very common in quantum gases, in particular with intensive
experimental and theoretical studies on ultracold spinor bosons Yuki ; Mie ;
Sad ; Chen ; Wang . The purpose of the present work is to develop PIMD for
two-component bosons with the same intracomponent and intercomponent
interaction, and it is straightforward to generalize our method to different
intracomponent and intercomponent interactions. Different from previous PIMD
simulation Hirshberg ; Xiong ; HirshbergFermi ; Xiong2 ; Xiong3 ; Xiong4  for
identical particles in two dimensions, here we consider spinor bosons in three
dimensions and study the number of beads (in PIMD simulation) needed in
numerical simulation to achieve numerical convergence. As an example, we
consider the thermodynamics of two-component bosons in a three-dimensional
harmonic trap. The energy and heat capacity are calculated to show the role of
temperature and repulsive interaction in the thermodynamics of ultracold
spinor bosons. We also generalize previous method Xiong ; Xiong2  to calculate
Green’s function and momentum distribution of spinor bosons to pave the way of
the general application of PIMD for ultracold spinor Bose gases.

##  II The partition function for two-component spinor bosons

We consider in this work spinor bosons with two internal states denoted by
$\left|\uparrow\right\rangle$ and $\left|\downarrow\right\rangle$. It is
straightforward to generalize to more general M-component spinor bosons. We
consider the following Hamiltonian operator which will not change the internal
state of a boson:

$\hat{H}=\sum_{j=1}^{N_{\uparrow}}\frac{\hat{\mathbf{p}}_{j}^{2}}{2m}+\sum_{j^{\prime}=1}^{N_{\downarrow}}\frac{\hat{\mathbf{p}}_{j^{\prime}}^{2}}{2m}+V\left(\hat{\mathbf{r}}_{1},\cdots,\hat{\mathbf{r}}_{N}\right).$
(1)

Here $N_{\uparrow}$ and $N_{\downarrow}$ denote the particle numbers in the
internal states $|\uparrow>$ and $|\downarrow>$, respectively.
$N=N_{\uparrow}+N_{\downarrow}$ is the total particle number of the quantum
system. Here $V$ includes both the external potential and inter-particle
interaction.

At finite temperature $\beta=1/k_{B}T$, the partition function is

$Z_{B}=Tr\left(e^{-\beta\hat{H}}\right).$ (2)

The purpose of the present work is to develop the PIMD to calculate the
thermodynamics of spinor bosons based on this partition function.

For two-component spinor bosons, we define the following basis which is
invariant about the exchange of any two bosons:

$\left|N_{B}\right\rangle\sim\sum_{p\in
S_{N}}\left|p\left\\{\textbf{r}_{1}^{\uparrow},\cdots,\textbf{r}_{N_{\uparrow}}^{\uparrow};\textbf{r}_{N_{\uparrow}+1}^{\downarrow},\cdots,\textbf{r}_{N}^{\downarrow}\right\\}\right\rangle.$
(3)

Here $p$ represents $N!$ permutation operations on the index
$\left\\{1,2,\cdots,N_{\uparrow},N_{\uparrow}+1,\cdots,N\right\\}$. We do not
change the order of the index for the internal state in the permutation
summation.

In this case, the partition function $Z_{B}$ for bosons is

$Z_{B}\sim\int\prod_{j=1}^{N}d\textbf{r}_{j}\left\langle
N_{B}\right|e^{-\Delta\beta\hat{H}}\cdots
e^{-\Delta\beta\hat{H}}\left|N_{B}\right\rangle.$ (4)

Here $\Delta\beta=\beta/P$ and $P$ is the number of beads in PIMD. One of the
difficulties to calculate the partition function is the $N!\times N!$
summation in the above expression.

Assume

$\left|N\right\rangle=\left|\textbf{r}_{1}^{\uparrow},\cdots,\textbf{r}_{N_{\uparrow}}^{\uparrow};\textbf{r}_{N_{\uparrow}+1}^{\downarrow},\cdots,\textbf{r}_{N}^{\downarrow}\right\rangle$
(5)

without permutation summation, we have

$Z_{B}\sim\int\prod_{j=1}^{N}d\textbf{r}_{j}\left\langle
N\right|e^{-\Delta\beta\hat{H}}\cdots
e^{-\Delta\beta\hat{H}}\left|N_{B}\right\rangle.$ (6)

Now we consider the following identity operator

$\hat{I}_{B}\sim\int\prod_{j=1}^{N}d\textbf{r}_{j}\left|N_{B}\right\rangle\left\langle
N_{B}\right|.$ (7)

By inserting a series of this identity operator in Eq. (6), the partition
function becomes:

$Z_{B}\sim\int\prod_{j=1}^{N}d\textbf{r}_{j}\left\langle
N\right|e^{-\Delta\beta\hat{H}}\hat{I}_{B}e^{-\Delta\beta\hat{H}}\cdots
e^{-\Delta\beta\hat{H}}\hat{I}_{B}e^{-\Delta\beta\hat{H}}\hat{I}_{B}\left|N_{B}\right\rangle.$
(8)

Without considering the exchange symmetry, we also have another form of the
identity operator

$\hat{I}=\int\prod_{j=1}^{N}d\textbf{r}_{j}\left|N\right\rangle\left\langle
N\right|.$ (9)

In this case, it is easy to prove that

$e^{-\beta\hat{H}}\hat{I}_{B}=e^{-\frac{\beta}{P}\hat{H}}\hat{I}e^{-\frac{\beta}{P}\hat{H}}\hat{I}\cdots\hat{I}e^{-\frac{\beta}{P}\hat{H}}\hat{I}_{B}.$
(10)

We get

$Z_{B}\sim\int\prod_{j=1}^{N}d\textbf{r}_{j}\left\langle
N\right|e^{-\Delta\beta\hat{H}}\hat{I}e^{-\Delta\beta\hat{H}}\hat{I}\cdots\hat{I}e^{-\Delta\beta\hat{H}}\hat{I}e^{-\Delta\beta\hat{H}}\left|N_{B}\right\rangle.$
(11)

We see that for Hamiltonian operator which does not change the internal state
of the bosons, in calculating the partition function, we should consider all
the permutations due to the exchange symmetry for bosons in the same internal
state. In Fig. 1, we illustrate an example of 3 bosons in internal state
$\left|\uparrow\right\rangle$ and 3 bosons in internal state
$\left|\downarrow\right\rangle$. The dashed line in this figure illustrates
the interaction between particles, while the solid line shows the connection
of the beads for the same particle. It is worthy to point out that in the
integral $\int\prod_{j=1}^{N}d\textbf{r}_{j}$, the coordinates are the same
for $\left|N\right>$ and $\left|N_{B}\right\rangle$. Hence, the permutation
summation in $\left|N_{B}\right\rangle$ and closed ring polymer should be
considered simultaneously.

Figure 1:  Illustration of the ring polymers for the partition function of
spinor bosons with two internal states $|\uparrow>$ and $|\downarrow>$. The
solid line shows the connection of the beads for the same particle, while the
dashed line shows the inter-particle interaction between particles at the same
imaginary time.

There are $N_{\uparrow}!\times N_{\downarrow}!$ permutation summations for the
partition function shown in Fig. 1, which increases exponentially with $N$.
Fortunately, in a recent pioneering work Hirshberg , Hirshberg, Rizzi and
Parrinello found that the partition function corresponding to a system of
identical bosons without internal freedom can be calculated by a recursion
formula which can be evaluated in $O(PN^{3})$ time, providing an exponential
speedup over simply summing all configurations for different permutations.

Generalizing the recursion formula for single-component bosons, we get the
following expression of the partition function for two-component bosons

$Z_{B}=\left(\frac{mP}{2\pi\hbar^{2}\beta}\right)^{PdN/2}\int
d\textbf{R}_{1}...d\textbf{R}_{N}e^{-\beta U_{B}^{(N)}}.$ (12)

The system under consideration has $d$ spatial dimensions. Here
$\mathbf{R}_{i}$ represents the collection of ring polymer coordinates
$(\mathbf{r}_{i}^{1},...,\mathbf{r}_{i}^{P})$ corresponding to the $i$th
particle. In addition,

$U_{B}^{(N)}=-\frac{1}{\beta}\ln W_{B}^{(N_{\uparrow})}-\frac{1}{\beta}\ln
W_{B}^{(N_{\downarrow})}+\frac{1}{P}\sum_{j=1}^{P}V\left(\textbf{r}_{1}^{j},...,\textbf{r}_{N}^{j}\right).$
(13)
$W_{B}^{(N_{\uparrow})}=\frac{1}{N_{\uparrow}}\sum_{k=1}^{N_{\uparrow}}e^{-\beta
E_{N_{\uparrow}}^{(k)}}W_{B}^{(N_{\uparrow}-k)}.$ (14)

$E_{N_{\uparrow}}^{(k)}=\frac{1}{2}m\omega_{P}^{2}\sum_{l=N_{\uparrow}-k+1}^{N_{\uparrow}}\sum_{j=1}^{P}\left(\textbf{r}_{l}^{j+1}-\textbf{r}_{l}^{j}\right)^{2}.$
(15)

Here $\textbf{r}_{l}^{P+1}=\textbf{r}_{l+1}^{1}$, except for $l=N_{\uparrow}$
for which $\textbf{r}_{N_{\uparrow}}^{P+1}=\textbf{r}_{N_{\uparrow}-k+1}^{1}$.
In addition, $\omega_{P}=\sqrt{P}/\beta\hbar$.

$W_{B}^{(N_{\downarrow})}=\frac{1}{N_{\downarrow}}\sum_{k=1}^{N_{\downarrow}}e^{-\beta
E_{N_{\downarrow}}^{(k)}}W_{B}^{(N_{\downarrow}-k)}.$ (16)

$E_{N_{\downarrow}}^{(k)}=\frac{1}{2}m\omega_{P}^{2}\sum_{l=N_{\downarrow}-k+1}^{N_{\downarrow}}\sum_{j=1}^{P}\left(\textbf{r}_{N_{\uparrow}+l}^{j+1}-\textbf{r}_{N_{\uparrow}+l}^{j}\right)^{2}.$
(17)

Here $\textbf{r}_{N_{\uparrow}+l}^{P+1}=\textbf{r}_{N_{\uparrow}+l+1}^{1}$,
except for $l=N_{\downarrow}$ for which
$\textbf{r}_{N}^{P+1}=\textbf{r}_{N-k+1}^{1}$.

The expression of $W_{B}^{(N_{\uparrow})}$ ($W_{B}^{(N_{\downarrow})}$)
includes the exchange effect by describing all the possible ring polymer
configurations for $N$ particles in the internal state
$\left|\uparrow\right\rangle$ ($\left|\downarrow\right\rangle$). The recursion
formula for $W_{B}^{(N_{\uparrow})}$ and $W_{B}^{(N_{\downarrow})}$ is the
same as that given in Ref. Hirshberg  for bosons without considering internal
state. We see that the evaluation of the partition function still takes
$O(PN^{3})$ time, rather than $O(PN_{\uparrow}^{3}N_{\downarrow}^{3})$. For
two-component spinor bosons there are in total $N_{\uparrow}!N_{\downarrow}!$
terms in the partition function summation for different configurations.
However, because spin degree of freedom doesn’t change, the spin up
contribution and spin down contribution to the partition function are
unrelated, and the partition function can be divided into two parts and
evaluated separately, so the evaluation of the partition function presented in
this work takes $O(PN_{\uparrow}^{3}+PN_{\downarrow}^{3})$ times, which is
still $O(PN^{3})$. However, if the spin degree of freedom has its own dynamics
due to spin-changing interaction, then there is no such shortcut and the
evaluation of the full partition function now takes
$O(PN_{\uparrow}^{3}N_{\downarrow}^{3})$ times.

After we get the final expression of $Z_{B}$, we may use massive Nosé-Hoover
chain Nose1 ; Nose2 ; Hoover ; Martyna ; Jang  to establish constant
temperature for the system, where each degree of freedom of the system has
been coupled to a separate Nosé-Hoover thermostat. From the distribution
generated by this method, we can calculate the density distribution and energy
of the system. For example, the average density is simply given by

$\rho_{\uparrow}(\textbf{x})=\left<\frac{1}{P}\sum_{j=1}^{P}\sum_{k=1}^{N_{\uparrow}}\delta(\mathbf{r}_{k}^{j}-\textbf{x})\right>,\rho_{\downarrow}(\textbf{x})=\left<\frac{1}{P}\sum_{j=1}^{P}\sum_{k=N_{\uparrow}+1}^{N}\delta(\mathbf{r}_{k}^{j}-\textbf{x})\right>.$
(18)

The energy estimator can be obtained from the following expression:

$E(\beta)=-\frac{1}{Z_{B}}\frac{\partial Z_{B}}{\partial\beta}.$ (19)

The energy estimator is

$E(\beta)=\frac{PdN}{2\beta}+\frac{\left<U\right>}{P}+\left<V_{B}^{(N_{\uparrow})}+\beta\frac{\partial
V_{B}^{(N_{\uparrow})}}{\partial\beta}\right>+\left<V_{B}^{(N_{\downarrow})}+\beta\frac{\partial
V_{B}^{(N_{\downarrow})}}{\partial\beta}\right>.$ (20)

where $V_{B}^{(N_{\uparrow})}=-\frac{1}{\beta}\log W_{B}^{(N_{\uparrow})}$ and
$V_{B}^{(N_{\downarrow})}=-\frac{1}{\beta}\log W_{B}^{(N_{\downarrow})}$.
$V_{B}^{(N_{\uparrow})}+\beta\frac{\partial
V_{B}^{(N_{\uparrow})}}{\partial\beta}$ may be evaluated as

$V_{B}^{(N_{\uparrow})}+\beta\frac{\partial
V_{B}^{(N_{\uparrow})}}{\partial\beta}=\frac{\sum_{k=1}^{N_{\uparrow}}\left[V_{B}^{(N_{\uparrow}-k)}+\beta\frac{\partial
V_{B}^{(N_{\uparrow}-k)}}{\partial\beta}-E_{N_{\uparrow}}^{(k)}\right]e^{-\beta(E_{N_{\uparrow}}^{(k)}+V_{B}^{(N_{\uparrow}-k)})}}{\sum_{k=1}^{N_{\uparrow}}e^{-\beta(E_{N_{\uparrow}}^{(k)}+V_{B}^{(N_{\uparrow}-k)})}},$
(21)

with $V_{B}^{(0)}+\beta\frac{\partial V_{B}^{(0)}}{\partial\beta}=0$. Similar
formulas hold for $V_{B}^{(N_{\downarrow})}+\beta\frac{\partial
V_{B}^{(N_{\downarrow})}}{\partial\beta}$. The interaction potential $U$ is
given by

$U=\sum_{j=1}^{P}V\left(\textbf{r}_{1}^{j},...,\textbf{r}_{N}^{j}\right).$
(22)

One may refer to our open-source code with C++ to see the details of the
implementation for the simulation of the partition function, energy estimator
and density estimator, etc..

##  III Results for two-component spinor bosons

In our PIMD simulation, the number of beads used decreases as temperature
increases to ensure numerical stability and assure convergence. Compared with
two dimensions with $P=12/T$ beads in the usual case Xiong ; Xiong2 ; Xiong3 ,
it is worthy to point out that many more beads are needed to assure accuracy
and convergence, because the size of phase space has been greatly enlarged
compared with two dimensional case. Based on the relation
$\Delta\beta=\beta/P$ and the fact that the evolution operator has been
expanded to an order of $O(\Delta\beta^{2})$, we would like to keep
$\Delta\beta$ a constant with different $\beta$, which leads to the choice
P=70/T. In general, of course, the convergence with respect to the number of
beads has to be explicitly verified by increasing the number of beads and
check if the results remain the same. It’s true that in this way one finds
that non-interacting case requires less beads than the interacting case, but
we choose the same number of beads for all cases while ensuring that the
results converge with respect to the number of beads. Moreover, we performed
12 independent simulations each with at least $10^{6}$ MD steps to establish
thermal equilibrium and calculate expectation values. In all of the following
we checked convergence with respect to the number of beads and MD steps
performed, by increasing P and MD steps and check if we obtain the same
results. In our simulation, the statistical fluctuations are negligible.
Hence, there is no error bar in all figures for numerical results.

###  III.1 Noninteracting spinor Bose gas in 3D harmonic trap

To test our algorithm, we first consider noninteracting identical bosons for
$N_{\uparrow}=N_{\downarrow}=7$ in a three-dimensional harmonic trap with
angular frequency $\omega$. We use the natural units of $\hbar=1,m=1$ and
consider $\omega=1$. The dimensionless Hamiltonian operator is

$\hat{H}=\frac{1}{2}\sum_{j=1}^{N}\hat{\mathbf{p}}_{j}^{2}+\frac{1}{2}\sum_{j=1}^{N}\hat{\textbf{r}}_{j}^{2}.$
(23)

In Fig. 2(a), the orange circle gives the average energy in our numerical
simulation for different temperatures, while the blue circle gives the average
energy in grand canonical ensemble for ideal bosons. In this work, the heat
capacity is calculated by the finite-difference derivative for the energy with
respect to the temperature. In Fig. 2(b), we give the heat capacity $C(T)$
based on the average energy calculated by PIMD. Our simulation shows the
presence of a bump in the heat capacity. In the thermodynamic limit, the
critical temperature Bagnato ; Pitaevskii  is determined by
$k_{B}T_{c}=0.94\hbar(\omega_{x}\omega_{y}\omega_{z})^{1/3}N_{\uparrow}^{1/3}$.
We have $T_{c}\approx 1.8$ with our parameters. The peak in our simulation for
the heat capacity is at $T_{c}\approx 1.53$, which is smaller than the result
with the formula of the thermodynamic limit. This comparison agrees
qualitatively with the result that the finite-size effect has the effect of
decreasing the critical temperature Shujuan ; Xiong-finite .

Figure 2:  For $N_{\uparrow}=N_{\downarrow}=7$ non-interacting bosons in
three-dimensional harmonic trap, (a) shows the average energy by PIMD for
spinor bosons (orange circle) and the average energy calculated by grand
canonical ensemble (blue circle). The good agreement proves the validity and
accuracy of our PIMD simulation. (b) shows the heat capacity calculated from
the average energy of PIMD, which shows a clear bump as the temperature
decreases. The statistical fluctuations are negligible.

In Fig. 3, for $k_{B}T=0.3$ (in unit of $\hbar\omega$) we give the numerical
result of the normalized density distribution $\rho_{\uparrow}(x)$ along $x$
direction by solid circle, while the solid line gives the analytical result of
the ground state. For $k_{B}T<<\hbar\omega$, the density distribution has the
following normalized analytical approximation:

$\rho_{\uparrow}(x)\simeq\left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-\frac{m\omega}{\hbar}x^{2}}.$
(24)

Good agreement is found with our simulation results.

Figure 3:  For $N_{\uparrow}=N_{\downarrow}=7$ non-interacting bosons in a
three-dimensional harmonic trap at the temperature $k_{B}T=0.3$, the circle
shows the density distribution $\rho_{\uparrow}(x)$ by PIMD, while the solid
line is the analytical result of the ground state density distribution. The
statistical fluctuations are negligible.

###  III.2 Spinor Bose gas in three-dimensional harmonic trap with repulsive
Gaussian interaction

The great merit of PIMD for identical bosons lies in that the inclusion of
inter-particle interaction will not add too much additional difficulty and
calculation error.

We consider the following Gaussian interaction Mujal  between bosons

$V(\mathbf{r}_{1},...,\mathbf{r}_{N})=\frac{1}{2}\sum_{j\neq
j^{\prime}=1}^{N}\frac{g}{\gamma^{2}\pi}e^{-|{\textbf{r}}_{j}-{\textbf{r}}_{j^{\prime}}|^{2}/\gamma^{2}}.$
(25)

We consider here $\gamma=1$. With length unit $\sqrt{\hbar/m\omega}$ and
energy unit $\hbar\omega$, we consider the average energy for
$N_{\uparrow}=N_{\downarrow}=7$ with different coupling constant $g$. In Fig.
4(a), we show the average energy for different $g$, while in Fig. 4(b) the
heat capacity for different $g$ is given. As expected, at the same
temperature, the repulsive interaction has the effect of increasing the total
energy. Significant role of the inter-particle repulsive interaction is shown
by the heat capacity, which is still quite challenging to reveal in mean-field
theory or the theory beyond the mean-field theory. Our accurate simulations
show that the repulsive inter-particle interaction will gradually erase the
bump in the heat capacity of ideal bosons. This unique result will provide the
benchmark for other methods in future to calculate accurately the
thermodynamics of interacting bosons. The disappearance of the bump in the
heat capacity is due to the strong interaction between particles.

Figure 4:  For $N_{\uparrow}=N_{\downarrow}=7$ bosons in a three-dimensional
harmonic trap with $g=0$ (blue), $g=1$ (orange) and $g=3$ (green), (a) shows
the average energy for different $g$, while (b) shows the heat capacity for
different $g$. As the inter-particle coupling constant increases, we see a
gradual disappearance of the bump in the heat capacity. The statistical
fluctuations are negligible.

We consider here Gaussian interaction, because for single-component bosons
with Gaussian interaction in two dimensions, an analytical expression is given
by Mujal et al. Mujal  for $N\leq 4$ bosons. Our simulation may give the
benchmark for other method to calculate accurately the energy of spinor bosons
in three dimensions in future. The s-wave scattering length $a$ may estimated
with $\int
d^{3}\textbf{r}\frac{g}{\gamma^{2}\pi}e^{-|\textbf{r}|^{2}/\gamma^{2}}=4\pi\hbar^{2}a/m$
Pitaevskii . In this case, we have the s-wave scattering length
$a=\frac{g\gamma}{4\sqrt{\pi}}$. It is worth pointing out that, to compare
with experiments of cold bosonic atoms with s-wave scattering length $a$, the
Gaussian interaction should be chosen to satisfy $a/\bar{l}<<1$, with
$\bar{l}$ the average interparticle distance.

In Fig. 5, we give the simulation result of the average energy for
$N_{\uparrow}=N_{\downarrow}=3,4,5$ with the coupling constant $g=3$.

Figure 5:  The simulated energy for $N_{\uparrow}=N_{\downarrow}=$ 3 (blue), 4
(orange), 5 (green) with the coupling constant $g=3$. The statistical
fluctuations are negligible.

###  III.3 Spinor Bose gas in three-dimensional harmonic trap with dipole
interaction

Now we consider dipole interaction which is the real situation of spinor Bose
gases Yuki . We consider the following simplest dipole interaction

$V(\mathbf{r}_{1},...,\mathbf{r}_{N})=\frac{1}{2}\sum_{j\neq
j^{\prime}=1}^{N}\frac{\lambda}{|{\textbf{r}}_{j}-{\textbf{r}}_{j^{\prime}}|^{3}}.$
(26)

Usually, the realistic dipole interaction comprises of an anisotropic term
Stuhler . The above dipole interaction is also used in path integral Monte
Carlo simulations of fermions Dornheim . We use this dipole interaction so
that it has the chance to provide the benchmark for other relevant studies.

Figure 6:  For $N_{\uparrow}=N_{\downarrow}=7$ bosons in a three-dimensional
harmonic trap with dipole interaction, we show the result of different
coupling constant $\lambda=0$ (blue), $\lambda=0.02$ (orange) and
$\lambda=0.1$ (green). (a) shows the average energy for different $\lambda$,
while (b) shows the heat capacity for different $\lambda$. As the inter-
particle coupling constant increases, we see a gradual disappearance of the
bump in the heat capacity, which is similar to the case of Gaussian
interaction. The statistical fluctuations are negligible.

By varying the interaction strength $\lambda$ we observed similar phenomenon
as we did with Gaussian interaction, as shown in Fig. 6.

##  IV Green’s function and momentum distribution

Now we turn to consider the momentum distribution of spinor bosons. In our
previous work, we calculated the momentum distribution for single-component
bosons in PIMD by first considering Green’s function of the system. Our method
can be generalized directly to the case of spinor bosons.

For spinor bosons, the thermal Green’s function is defined as

$G_{\uparrow\uparrow}(\textbf{x},\tau_{1};\textbf{y},\tau_{2})=\left<\mathcal{T}\left\\{\hat{\psi}_{\uparrow}(\mathbf{y},\tau_{2})\hat{\psi}_{\uparrow}^{\dagger}(\mathbf{x},\tau_{1})\right\\}\right>,$
(27)

where $\left<\cdots\right>$ denotes thermal average, $\mathcal{T}$ is the
imaginary time-ordering operator. In addition,

$\hat{\psi}_{\uparrow}(\mathbf{x},\tau)=e^{\hat{H}\tau}\hat{\psi}_{\uparrow}(\mathbf{x})e^{-\hat{H}\tau},\hat{\psi}^{\dagger}_{\uparrow}(\mathbf{x},\tau)=e^{\hat{H}\tau}\hat{\psi}^{\dagger}_{\uparrow}(\mathbf{x})e^{-\hat{H}\tau}.$
(28)

Of course, we have $G_{\uparrow\downarrow}=0$ and $G_{\downarrow\uparrow}=0$.

To get momentum distribution, we should calculate

$G_{\uparrow\uparrow}(\textbf{x},\tau_{2}+\Delta\beta;\textbf{y},\tau_{2})=\left<\hat{\psi}_{\uparrow}^{\dagger}(\mathbf{x},\tau_{2}+\Delta\beta)\hat{\psi}_{\uparrow}(\mathbf{y},\tau_{2})\right>.$
(29)

In Fig. 7, we show the idea to calculate this Green’s function. One may find
the details to implement the algorithm in our previous work. Green’s function
can be estimated as

$G_{\uparrow\uparrow}(\textbf{x}^{\prime},\tau_{1};\textbf{y}^{\prime},\tau_{2})=\left<\delta(\textbf{x}-\textbf{x}^{\prime})\delta(\textbf{y}-\textbf{y}^{\prime})\right>,$
(30)

where x and y denote the positions of two beads at the end of the gap in Fig.
7.

Figure 7:  Modified ring polymers for Green’s function $G_{\uparrow\uparrow}$.
There is a gap between two beads, from $\tau_{2}$ to $\tau_{2}+\Delta\beta$.
In particular, the coordinates of the beads at the ends of the gap are
$\mathbf{y}$ and $\mathbf{x}$. There is no harmonic spring between those two
beads.

As illustrated in Fig. 7, there would be a gap for the beads in the top row,
due to the field operator
$\hat{\psi}_{\uparrow}^{\dagger}(\mathbf{x},\tau_{2}+\Delta\beta)\hat{\psi}_{\uparrow}(\mathbf{y},\tau_{2})$.
Following the method given in Ref. Xiong , it is straightforward to write the
recursion formula for this Green’s function.

Based on thermal Green’s function, the momentum distribution is

$\rho_{\uparrow}(\textbf{p})=\frac{1}{(2\pi\hbar)^{d}}\int
d\textbf{x}d\textbf{y}G_{\uparrow\uparrow}(\textbf{x},\tau_{2}+\Delta\beta;\textbf{y},\tau_{2})e^{\frac{i}{\hbar}\textbf{p}\cdot(\textbf{x}-\textbf{y})}.$
(31)

It is similar to consider
$G_{\downarrow\downarrow}(\textbf{x},\tau_{2}+\Delta\beta;\textbf{y},\tau_{2})$
and $\rho_{\downarrow}(\textbf{p})$.

As an application of our method, we consider $N_{\uparrow}=N_{\downarrow}=8$
bosons in 3D harmonic trap with different dipole interaction at $k_{B}T=0.3$.
In Fig. 8, we show the momentum distribution along $x$ direction for different
$\lambda$. We see that as $\lambda$ increases, the momentum distribution
becomes narrower, as expected from the physical picture that repulsive
interactions have the effect of increasing the width of the quantum system. At
the same temperature, this effect of inter-particle interaction for momentum
distribution may be also interpreted by considering Green’s function in real
space. The ground state wave function is a Gaussian function for non-
interacting case, while in the Thoms-Fermi approximation the quadratic
behavior is dominating as interaction increases so that it becomes flatter
Pitaevskii . In experiments, one may use optical lattices lattice  so that in
a lattice site there exists a few or a dozen spinor bosons. By using magnetic-
field Feshbach resonance Feshbach  to turn the inter-particle interaction, we
have the chance to show experimentally the behavior of the heat capacity
simulated in this work.

Figure 8:  Shown is the momentum distribution $\rho(p_{x})$ of
$N_{\uparrow}=N_{\downarrow}=8$ bosons in 3D harmonic trap with different
dipole interaction at $k_{B}T=0.3$. The blue circle is the momentum
distribution of $\lambda=0$, while the green circle and orange circle are for
$\lambda=0.5$ and $\lambda=1$, respectively. The solid line is the momentum
distribution of the ground state wave function, which agrees well with our
simulation at low temperature. The statistical fluctuations are negligible.

##  V Conclusion and discussion

As a summary, in this work we showed how to generalize the original PIMD
methodology to consider spinor bosons and extract thermal properties from our
simulations, such as density, heat capacity and momentum distribution. By
using our exact numerical method, we are able to study the behavior of heat
capacity and we observed that as inter-particle interaction strength
increases, the bump in heat capacity curve is gradually erased. Issues worth
considering in the future are the compressibility in a harmonic trap, which
has been calculated by path integral Monte Carlo for a homogeneous dilute Bose
gas and a discontinuous behavior is found at the critical point Spada . We
also suggest possible experiment scheme to observe this behavior
experimentally. Of course, the method we developed here can be generalized
straightforwardly to consider bosons with multiple components and other
systems, providing opportunities for future research.

###### Acknowledgements.

This work is partly supported by the National Natural Science Foundation of
China under grant numbers 11175246, and 11334001.

DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request. The code of this study is openly
available in GitHub (https://github.com/xiongyunuo/PIMD-Pro-SB-MT).

## References

  * (1) M. Rossi, L. Salasnich, F. Ancilotto, and F. Toigo, Monte Carlo simulations of the unitary Bose gas, Phys. Rev. A 89, 041602(R) (2014).
  * (2) L. Parisi, G. E. Astrakharchik, and S. Giorgini, The liquid state of one-dimensional Bose mixtures: a quantum Monte-Carlo study, Phys. Rev. Lett. 122, 105302 (2019).
  * (3) L. Parisi and S. Giorgini, Quantum droplets in one-dimensional Bose mixtures: a quantum Monte-Carlo study, Phys. Rev. A 102, 023318 (2020).
  * (4) R. P. Feynman and A. R. Hibbs, Quantum mechanics and path integrals, Dover Publications, New York (2010).
  * (5) H. Kleinert, Path integrals in quantum mechanics, statistics, polymer physics, and financial markets, World Scientific, Singapore (2009).
  * (6) M. E. Tuckerman, Statistical mechanics: theory and molecular simulation, Oxford University, New York (2010).
  * (7) D. Chandler and P. G. Wolynes, Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids, J. Chem. Phys. 74, 4078 (1981).
  * (8) M. Parrinello and A. Rahman, Study of an F center in molten KCl, J. Chem. Phys. 80, 860 (1984).
  * (9) S. Miura and S. Okazaki, Path integral molecular dynamics for Bose-Einstein and Fermi-Dirac statistics. J. Chem. Phys. 112, 10116 (2000).
  * (10) J. Cao and G. A. Voth, The formulation of quantum statistical mechanics based on the Feynman path centroid density. I. Equilibrium properties, J. Chem. Phys. 100, 5093 (1994).
  * (11) J. Cao and G. A. Voth, The formulation of quantum statistical mechanics based on the Feynman path centroid density. II. Dynamical properties, J. Chem. Phys. 100, 5106 (1994).
  * (12) S. Jang and G. A. Voth, A derivation of centroid molecular dynamics and other approximate time evolution methods for path integral centroid variables, J. Chem. Phys. 111, 2371 (1999).
  * (13) R. RamíRez and T. LóPez-Ciudad, The Schrödinger formulation of the Feynman path centroid density, J. Chem. Phys. 111, 3339 (1999).
  * (14) K. Kinugawa, H. Nagao, and K. Ohta, Path integral centroid molecular dynamics method for Bose and Fermi statistics: formalism and simulation, Chem. Phys. Lett. 307, 187 (1999).
  * (15) Pierre-Nicholas Roy, Seogjoo Jang, and Gregory A. Voth, Feynman path centroid dynamics for Fermi–Dirac statistics, J. Chem. Phys. 111, 5303 (1999).
  * (16) P.-N. Roy and G.A. Voth, On the Feynman path centroid density for Bose-Einstein and Fermi-Dirac statistics, J. Chem. Phys. 110, 3647 (1999).
  * (17) N. Blinov, P.-N. Roy, and G.A. Voth, Path integral formulation of centroid dynamics for systems obeying Bose-Einstein statistics, J. Chem. Phys. 115, 4484 (2001).
  * (18) N. Blinov and P.-N. Roy, Operator formulation of centroid dynamics for Bose-Einstein and Fermi-Dirac statistics, J. Chem. Phys. 115, 7822 (2001).
  * (19) K. Kinugawa, A semiclassical approach to the dynamics of many-body Bose/Fermi systems by the path integral centroid molecular dynamics, J. Chem. Phys. 114, 1454 (2001).
  * (20) K. Kinugawa, H. Nagao, and K. Ohta, A path integral centroid molecular dynamics method for Bose and Fermi statistics, J. Mol. Liq. 90, 11 (2001).
  * (21) Nicholas Blinov and Pierre-Nicholas Roy, An effective centroid Hamiltonian and its associated centroid dynamics for indistinguishable particles in a harmonic trap, J. Chem. Phys. 116, 4808 (2002).
  * (22) Pierre-Nicholas Roy and Nicholas Blinov, Centroid dynamics with quantum statistics, Isr J. Chem. 42, 183 (2002).
  * (23) Paul Moffatt, Nicholas Blinov, and Pierre-Nicholas Roy, On the calculation of single-particle time correlation functions from Bose–Einstein centroid dynamics, J. Chem. Phys. 120, 4614 (2004).
  * (24) E. A. Polyakov, A. P. Lyubartsev, and P. N. Vorontsov-Velyaminov, Centroid molecular dynamics: Comparison with exact results for model systems, J. Chem. Phys. 133, 194103 (2010).
  * (25) I. R. Craig and D. E. Manolopoulos, Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics, J. Chem. Phys. 121, 3368 (2004).
  * (26) B. J. Braams and D. E. Manolopoulos, On the short-time limit of ring polymer molecular dynamics, J. Chem. Phys. 125, 124105 (2006).
  * (27) S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller 3rd, Ring-polymer molecular dynamics: quantum effects in chemical dynamics from classical trajectories in an extended phase space, Annu. Rev. Phys. Chem. 64, 387 (2013).
  * (28) T. E. Markland and M. Ceriotti, Nuclear quantum effects enter the mainstream, Nat. Rev. Chem. 2, 0109 (2018).
  * (29) B. Hirshberg, V. Rizzi, and M. Parrinello, Path integral molecular dynamics for bosons, Proc. Natl. Acad. Sci. U. S. A. 116, 21445 (2019).
  * (30) Y. N. Xiong and H. W. Xiong, Path integral molecular dynamics simulations for Green’s function in a system of identical bosons, J. Chem. Phys. 156, 134112 (2022).
  * (31) C. W. Myung, B. Hirshberg, and M. Parrinello, Prediction of a supersolid phase in high-pressure deuterium, Phys. Rev. Lett. 128, 045301 (2022).
  * (32) Y. N. Xiong and H. W. Xiong, Path integral and winding number in singular magnetic field, Eur. Phys. J. Plus 137 550 (2022).
  * (33) B. Hirshberg, M. Invernizzi, and M. Parrinello, Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality, J. Chem. Phys. 152, 171102 (2020).
  * (34) Y. N. Xiong and H. W. Xiong, Numerical calculation of Green’s function and momentum distribution for spin-polarized fermions by path integral molecular dynamics, J. Chem. Phys. 156, 204117 (2022).
  * (35) Y. N. Xiong and H. W. Xiong, A solution of fermion sign problem for large fermion systems, arXiv:2206.08341 (2022).
  * (36) Yuki Kawaguchia and Masahito Ueda, Spinor Bose–Einstein condensates, Physics Reports 520, 253 (2012).
  * (37) H. -J. Miesner, D. M. Stamper-Kurn, J. Stenger, S. Inouye, A. P. Chikkatur, and W. Ketterle, Observation of Metastable States in Spinor Bose-Einstein Condensates, Phys. Rev. Lett. 82, 2228 (1999).
  * (38) L. Sadler, J. Higbie, S. Leslie, M. Vengalattore, and D. M. Stamper-Kurn, Spontaneous symmetry breaking in a quenched ferromagnetic spinor Bose–Einstein condensate, Nature 443, 312 (2006).
  * (39) Li Chen, Han Pu, and Yunbo Zhang, Spin-orbit angular momentum coupling in a spin-1 Bose-Einstein condensate, Phys. Rev. A 93, 013629 (2016).
  * (40) Chunji Wang, Chao Gao, Chao-Ming Jian, and Hui Zhai, Spin-Orbit Coupled Spinor Bose-Einstein Condensates, Phys. Rev. Lett. 105, 160403 (2010).
  * (41) S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys. 52, 255 (1984).
  * (42) S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81, 511 (1984).
  * (43) W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985).
  * (44) G. J. Martyna, M. L. Klein, and M. Tuckerman, Nosé-Hoover chains: The canonical ensemble via continuous dynamics, J. Chem. Phys. 97, 2635 (1992).
  * (45) S. Jang and G. A. Voth, Simple reversible molecular dynamics algorithms for Nosé-Hoover chain dynamics, J. Chem. Phys. 107, 9514 (1997).
  * (46) V. Bagnato, David E. Pritchard, and D. Kleppner, Bose-Einstein condensation in an external potential, Phys. Rev. A 35, 4354 (1987).
  * (47) Lev. P. Pitaevskii and S. Stringari, Bose Einstein condensation and superfluidity, Oxford University Press, New York (2016).
  * (48) S. Liu, G. Huang, L. Ma, S. Zhu, and H. Xiong, Thermodynamic properties of the trapped interacting Bose gas with finite size, J. Phys. B: At. Mol. Opt. Phys. 33, 3911 (2000).
  * (49) H. Xiong, S. Liu, G. Huang, Z. Xu, and C. Zhang, Critical temperature and condensate fraction of the trapped interacting Bose gas with finite-size effects, J. Phys. B: At. Mol. Opt. Phys. 34, 3013 (2001).
  * (50) P. Mujal, E. Sarlé, A. Polls, B. Juliá-Díaz, Quantum correlations and degeneracy of identical bosons in a two-dimensional harmonic trap, Phys. Rev. A 96, 043614 (2017).
  * (51) T. Dornheim, Fermion sign problem in path integral Monte Carlo simulations: Quantum dots, ultracold atoms, and warm dense matter, Phys. Rev. E 100, 023307 (2019).
  * (52) J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Observation of Dipole-Dipole Interaction in a Degenerate Quantum Gas, Phys. Rev. Lett. 95, 150406 (2005).
  * (53) M. Greiner, O. Mandel, T. Esslinger, Theodor W. Hänsch, and I. Bloch, Nature 415, 39 (2002).
  * (54) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * (55) G. Spada, S. Pilati, and S. Giorgini, Thermodynamics of a dilute Bose gas: A path-integral Monte Carlo study, Phys. Rev. A 105, 013325 (2022).


# Simultaneous recovery of a locally rough interface and the embedded obstacle
with the reverse time migration

Jianliang Li School of Mathematics and Statistics, Hunan Normal University,
Changsha 410081, China<EMAIL_ADDRESS>Jiaqing Yang School of Mathematics
and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
<EMAIL_ADDRESS>
###### Abstract

Consider the inverse acoustic scattering of time-harmonic point sources by an
unbounded locally rough interface with bounded obstacles embedded in the lower
half-space. A novel version of reverse time migration is proposed to
reconstruct both the locally rough interface and the embedded obstacle. By a
modified Helmholtz-Kirchhoff identity associated with a planar interface, we
obtain a modified imaging functional which has been shown that it always peaks
on the local perturbation of the interface and on the embedded obstacle.
Numerical examples are presented to demonstrate the effectiveness of the
method.

Keywords: inverse acoustic scattering, locally rough interfaces, embedded
obstacles, reverse time migration.

##  1 Introduction

This paper concerns the two-dimensional acoustic scattering of time-harmonic
point sources by a locally rough interface with an embedded obstacle in the
lower half-space. Given the incident wave, the direct scattering problem is to
determine the distribution of the scattered wave; while the inverse scattering
problem aims to recover the locally rough interface as well as the embedded
obstacle from the measured scattered wave in a certain domain. These problems
was motivated by significant applications in diverse scientific areas such as
medical imaging [1] and exploration geophysics [3].

The main difficulty of the rough surface scattering problem is the
unboundedness of the rough surface, which makes that the related integral
operators are non-compact and so the classical Fredholm alternative is not
available. Based on a generalized Fredholm theory [23, 24], the well-posedness
of the direct scattering problem by rough surfaces has been established in
[19, 20, 21, 41, 46] by the integral equation method. In addition, we refer to
[18, 34, 42, 43, 44, 48] for the well-posedness of the direct rough surface
scattering problem. Different from all previous works, if the rough surface is
a local perturbation of a planar surface, by introducing a special locally
rough surface, the scattering problem can be transformed into an equivalent
integral equation defined in a bounded domain, for which the well-posedness
follows from the classical Fredholm alternative, for details we refer to [25,
36] for the scattering by the locally rough surface and [37, 45] for the
scattering by the locally rough interface with an embedded obstacle. For
inverse problems, the reference [45] has established a global uniqueness which
shows that the locally rough interface, the embedded obstacle and the wave
number in the lower half-space can be uniquely determined by means of near-
field measurements above the interface. Based on this uniqueness, a modified
linear sampling method has been developed in [37] to solve the inverse problem
of simultaneously reconstructing the locally rough interface and the embedded
obstacle. However, it is worth pointing out that the linear sampling method in
[37] is sensitive to the noise. For the inverse problem, there exists a
variety of numerical algorithms. For the inverse scattering by a planar
surface with buried objects, we refer to the MUSIC-type scheme [2], the
asymptotic factorization method [26], the sampling method [27], and the direct
imaging algorithm [32, 38]. For the inverse scattering by rough surfaces, we
refer to iterative algorithms [7, 8, 9, 22], the algorithm based on the
transformed field expansion [5, 6], the factorization method [28], the
singular source method [31], the direct imaging method [39, 40], and linear
sampling methods [25, 36, 47].

The reverse time migration (RTM) method is a popular sample-type method which
has been extensively applied in seismic imaging [4] and exploration geophysics
[3, 4, 10]. This method first back-propagates the complex conjugated data into
the background medium and then computes the cross-correlation between the
back-propagated field and the incident field to output the imaging indicator,
which has different behaviors when the sampling point is near the scatterer
and far away from the scatterer. Based on this property, it is efficient,
stable, and robust to noise. The key point of the RTM method is to establish
the related Helmholtz-Kirchhoff identity, which plays a crucial role in the
analysis of the indicator. For the inverse obstacle scattering problem, the
mathematical justification of the RTM method has been proved rigorously in
[11], which is based on an usual Helmholtz-Kirchhoff identity associated with
the fundamental solution of the Helmholtz operator in the free space. The
results in [11] has been extended to [12, 13, 14, 15, 16, 17, 29] to solve
some other inverse obstacle scattering problems. However, there are few
results for the RTM method to recover unbounded rough surfaces since the usual
Helmholtz-Kirchhoff identity in [11] is not valid. To overcome this
difficulty, we has established a modified Helmholtz-Kirchhoff identity
associated with a background Green function in a two-layered medium separated
by a special locally rough surface and then extended the RTM method to
reconstruct the locally rough surface without embedded obstacles, see [35] for
details.

Unfortunately, if there is a obstacle embedded in the lower half-space
separated by a locally rough interface, both the usual Helmholtz-Kirchhoff
identity in [11] and the modified Helmholtz-Kirchhoff identity in [35] are not
valid to simultaneously reconstruct the locally rough interface and the
embedded obstacle. Notice that the rough interface is locally perturbed, we
are able to introduce a planar interface and then transform the scattering
problem by the locally rough interface and the embedded obstacle into the
scattering problem by several inhomogeneous mediums and the embedded obstacle.
Based on this observation, we can establish the corresponding Helmholtz-
Kirchhoff identity associated with the background Green function in a two-
layered medium separated by the planar interface. Thus, the mathematical
justification of the RTM method has been proved rigorously, where we
illustrate that the corresponding imaging indicator always peaks on the local
perturbation of the interface and on the embedded obstacle, which is confirmed
in the numerical examples.

The outline of this paper is as follows. In section 2, we introduce the
mathematical model for the forward scattering problem, and present a novel
version of the RTM method based on a modified Helmholtz-Kirchhoff identity. In
section 3, numerical examples are reported to illustrate the effectiveness of
the proposed method. The paper is concluded with some general remarks and
discussions on the future research in section 4.

##  2 The RTM method

In this section, we shall introduce the mathematical model of the scattering
problem by an unbounded, locally rough interface and an obstacle in the lower
half-space, and investigate the RTM method for this model.

Let the scattering interface be described by a curve

$\displaystyle\Gamma:=\\{(x_{1},x_{2})\in\mathbb{R}^{2}:x_{2}=f(x_{1})\\}$
(2.1)

where $f$ is assumed to be a Lipschitz continuous function with compact
support. It means that the interface $\Gamma$ is a local perturbation of the
planar interface

$\displaystyle\Gamma_{0}:=\\{(x_{1},x_{2})\in\mathbb{R}^{2}:x_{2}=0\\}.$

The whole space is separated by $\Gamma$ into two half-spaces denoted by

$\displaystyle\Omega_{1}:=\\{(x_{1},x_{2})\in\mathbb{R}^{2}:x_{2}>f(x_{1})\\}\quad{\rm
and}\quad\Omega_{2}:=\\{(x_{1},x_{2})\in\mathbb{R}^{2}:x_{2}<f(x_{1})\\}.$

In the lower half-space $\Omega_{2}$, we assume that there is an bounded
obstacle $D$ with a smooth boundary $\partial D\in C^{2,\alpha}$ for some
Hölder exponent $0<\alpha\leq 1$. In this paper, for simplicity, $D$ is
assumed to be sound-soft which means that a Dirichlet boundary condition is
imposed on $\partial D$.

Consider the incident field $u^{i}(x,x_{s})$ to be generated by a point source

$\displaystyle u^{i}(x,x_{s})=\Phi_{\kappa}(x,x_{s}):=\frac{\rm
i}{4}H_{0}^{(1)}(\kappa|x-x_{s}|)\quad{\rm
for}\;\;x_{s}\in{\mathbb{R}}^{2}\setminus(\overline{D}\cup\Gamma).$ (2.2)

Here, $H_{0}^{(1)}$ is the Hankel function of the first kind of order zero,
$\kappa$ is the wave number which is a piecewise constant given by
$\kappa=\kappa_{1}$ for $x_{s}\in\Omega_{1}$ and $\kappa=\kappa_{2}$ for
$x_{s}\in\Omega_{2}\setminus\overline{D}$, and for some constant $\lambda$,
$\Phi_{\lambda}$ is the fundamental solution of the Helmholtz equation
satisfying
$\Delta\Phi_{\lambda}(x,x_{s})+\lambda^{2}\Phi_{\lambda}(x,x_{s})=-\delta_{x_{s}}(x)$
in $\mathbb{R}^{2}$ in the distributional sense, where
$\delta_{x_{s}}(x):=\delta(x-x_{s})$ is the Kronecker delta distribution. Then
the scattering of $u^{i}$ by the scatterer $(\Gamma,D)$ can be reduced to the
problem of seeking $u=u(x,x_{s})$ satisfying that

$\left\\{\begin{array}[]{lll}\Delta
u+\kappa^{2}u=-\delta_{x_{s}}&\textrm{in}\;\;{\mathbb{R}}^{2}\setminus\overline{D},\\\\[5.69054pt]
u=0&\textrm{on}\;\;\partial D,\\\\[5.69054pt]
\lim\limits_{|x|\rightarrow\infty}|x|^{\frac{1}{2}}\left(\partial_{|x|}u^{s}-{\rm
i}\kappa u^{s}\right)=0,\end{array}\right.$ (2.3)

in the distributional sense. Here, $u(x,x_{s})$ and $u^{s}(x,x_{s})$ are the
total field and the scattered field, respectively, which are related by

$u(x,x_{s})=\left\\{\begin{array}[]{ll}\Phi_{\kappa_{1}}(x,x_{s})+u^{s}(x,x_{s})&\textrm{in}\;\;\Omega_{1},\\\\[5.69054pt]
u^{s}(x,x_{s})&\textrm{in}\;\;\Omega_{2}\setminus\overline{D},\end{array}\right.$

for $x_{s}\in\Omega_{1}$, and

$u(x,x_{s})=\left\\{\begin{array}[]{ll}u^{s}(x,x_{s})&\textrm{in}\;\;\Omega_{1},\\\\[5.69054pt]
\Phi_{\kappa_{2}}(x,x_{s})+u^{s}(x,x_{s})&\textrm{in}\;\;\Omega_{2}\setminus\overline{D},\end{array}\right.$

for $x_{s}\in\Omega_{2}\setminus\overline{D}$, the wave number
$\kappa=\kappa(x)$ is a piecewise constant defined by

$\kappa(x)=\left\\{\begin{array}[]{ll}\kappa_{1}&\textrm{for}\;\;x\in\Omega_{1},\\\\[5.69054pt]
\kappa_{2}&\textrm{for}\;\;x\in\Omega_{2}\setminus\overline{D},\end{array}\right.$

and the last condition in Problem (2.3) is the well-known Sommerfeld radiation
condition which holds uniformly for all directions
$\hat{x}:=x/|x|\in{\mathbb{S}}:=\\{x\in{\mathbb{R}}^{2}:|x|=1\\}$. Moreover,
the Sommerfeld radiation condition allows that the scattered field $u^{s}$ has
the following asymptotic behavior

$\displaystyle u^{s}(x,x_{s})=\frac{e^{{\rm
i}\kappa|x|}}{|x|^{\frac{1}{2}}}\left\\{u^{\infty}(\hat{x},x_{s})+O\left(\frac{1}{|x|}\right)\right\\}\quad{\rm
for}\;|x|\to\infty$

uniformly in all direction $\hat{x}\in{\mathbb{S}}$, where
$u^{\infty}(\hat{x},x_{s})$ is known as the far-field pattern of
$u^{s}(x,x_{s})$. It is shown by Theorem 3.4 in [45] that Problem (2.3) is
well-posed.

Figure 1: The setting of RTM method for reconstruction of $(\Gamma,D)$.

As shown in Figure 1, we first introduce some notations used in the RTM
approach. Let $\Gamma$ be the locally rough interface given by (2.1), whose
local perturbations are denoted by $B_{1}:={\mathbb{R}}^{2}_{+}\cap\Omega_{2}$
and $B_{2}:={\mathbb{R}}^{2}_{-}\cap\Omega_{1}$ with
${\mathbb{R}}^{2}_{\pm}:=\\{x\in{\mathbb{R}}^{2}:x_{2}\gtrless 0\\}$. For
simplicity, in this paper we consider a simple case that $\Gamma$ has only two
local perturbations as seen in Figure 1. The results obtained in this paper is
also valid for the case of multiple local perturbations. For $j=1,2$, let
$\chi_{j}(x)$ denote the characterization function of the domain $B_{j}$,
defined by $\chi_{j}(x)=1$ in $B_{j}$ and $\chi_{j}(x)=0$ outside of $B_{j}$.
Throughout, denoted by $B:=B_{1}\cup B_{2}$ and
$\chi(x):=\chi_{1}(x)-\chi_{2}(x)$. Let $D$ be the obstacle embedded in
$\Omega_{2}$ and we assume a priori that $D\cap B=\emptyset$. Since the RTM
method is a sample-type method, we choose a rectangle sampling domain $S$ such
that $B\subset S$ and $D\subset S$. And we choose a sufficiently large $R$
such that $S\subset B_{R}$, where $B_{R}$ stands for the disc with the origin
as the center and $R$ as the radius. Assume that there are $N_{s}$ point
sources $x_{s}$ uniformly distributed on $\Gamma_{s}:=\partial B_{R}$ and
$N_{r}$ receivers $x_{r}$ uniformly distributed on $\Gamma_{r}:=\partial
B_{R}$.

To introduce the RTM method, we first consider the scattering of the incident
point source $u^{i}(x,x_{s})$ given by (2.2) by the planar interface
$\Gamma_{0}$, which reads

$\left\\{\begin{array}[]{ll}\Delta
G(x,x_{s})+\kappa_{0}^{2}(x)G(x,x_{s})=-\delta_{x_{s}}(x)&\textrm{in}\;\;{\mathbb{R}}^{2},\\\\[5.69054pt]
\lim\limits_{|x|\rightarrow\infty}|x|^{\frac{1}{2}}\left(\partial_{|x|}G(x,x_{s})-{\rm
i}\kappa_{0}(x)G(x,x_{s})\right)=0,\end{array}\right.$ (2.4)

in the distributional sense with the Sommerfeld radiation condition uniformly
for all $\hat{x}\in{\mathbb{S}}$. Here, the wave number $\kappa_{0}(x)$ is
defined by $\kappa_{0}(x)=\kappa_{1}$ in ${\mathbb{R}}^{2}_{+}$ and
$\kappa_{0}(x)=\kappa_{2}$ in ${\mathbb{R}}^{2}_{-}$. We refer to [30, 45] for
the explicit expression of $G(x,x_{s})$.

For $x_{s}\in\Gamma_{s}$, we define

$V(x,x_{s}):=u(x,x_{s})-G(x,x_{s}),$ (2.5)

then it follows from (2.3) and (2.4) that $V(x,x_{s})$ solves

$\left\\{\begin{array}[]{lll}\Delta
V(x,x_{s})+\kappa^{2}(x)V(x,x_{s})=\beta\chi(x)G(x,x_{s})&\textrm{in}\;\;{\mathbb{R}}^{2}\setminus\overline{D},\\\\[5.69054pt]
V(x,x_{s})=-G(x,x_{s})&\textrm{on}\;\;\partial D,\\\\[5.69054pt]
\lim\limits_{|x|\rightarrow\infty}|x|^{\frac{1}{2}}\left(\partial_{|x|}V(x,x_{s})-{\rm
i}\kappa(x)V(x,x_{s})\right)=0,\end{array}\right.$ (2.6)

where $\beta:=\kappa_{1}^{2}-\kappa_{2}^{2}$. We refer to Theorem 3.1 in [37]
for the unique solvability of Problem (2.6).

Since we can obtain $G(x_{r},x_{s})$ by solving Problem (2.4) through the
Nyström method or the finite element method, we can obtain the data
$V(x_{r},x_{s})$ from the measurement $u(x_{r},x_{s})$ and (2.5). As mentioned
in the introduction, we first back-propagate the complex conjugated data
$\overline{V(x_{r},x_{s})}$ into the domain ${\mathbb{R}}^{2}$, and then
define the indicator as the imaginary part of the cross-correlation of
$G(\cdot,x_{s})$ and the back-propagation field. More precisely, we summarize
it in the following algorithm.

Algorithm 1 (RTM for locally rough interfaces and embedded obstacles): Given
the data $V(x_{r},x_{s})$ for $r=1,2,...,N_{r}$ and $s=1,2,...,N_{s}$.

  * •

Back-propagation: for $s=1,2,...,N_{s}$, solve the problem

$\displaystyle\left\\{\begin{aligned} &\Delta
W(x,x_{s})+\kappa_{0}^{2}(x)W(x,x_{s})=\frac{|\Gamma_{r}|}{N_{r}}\sum_{r=1}^{N_{r}}\overline{V(x_{r},x_{s})}\delta_{x_{r}}(x)\qquad\textrm{in}\;\;{\mathbb{R}}^{2},\\\
&\;\displaystyle\lim_{|x|\rightarrow\infty}|x|^{\frac{1}{2}}\left(\partial_{|x|}W(x,x_{s})-{\rm
i}\kappa_{0}(x)W(x,x_{s})\right)=0,\end{aligned}\right.$

to obtain the solution $W(x,x_{s})$.

  * •

Cross-correlation: for each sampling point $z\in S$, compute the indicator
function

$\displaystyle{\rm Ind}(z)=\kappa(x_{r}){\rm
Im}\left\\{\frac{|\Gamma_{s}|}{N_{s}}\sum_{s=1}^{N_{s}}\kappa(x_{s})G(z,x_{s})W(z,x_{s})\right\\}$

and then plot the mapping ${\rm Ind}(z)$ against $z$.

From Problem (2.4) and the linearity, we immediately see that

$\displaystyle
W(x,x_{s})=-\frac{|\Gamma_{r}|}{N_{r}}\sum_{r=1}^{N_{r}}G(x,x_{r})\overline{V(x_{r},x_{s})},$

which implies that

$\displaystyle{\rm Ind}(z)=-{\rm
Im}\left\\{\frac{|\Gamma_{s}|}{N_{s}}\frac{|\Gamma_{r}|}{N_{r}}\sum_{s=1}^{N_{s}}\sum_{r=1}^{N_{r}}\kappa(x_{r})\kappa(x_{s})G(z,x_{s})G(z,x_{r})\overline{V(x_{r},x_{s})}\right\\}\quad
z\in S.$ (2.7)

Observing that $G(z,x_{s})$, $G(z,x_{r})$ and $V(x_{r},x_{s})$ are continuous
for $z\in S$, $x_{r}\in\Gamma_{r}$, $x_{s}\in\Gamma_{s}$, it follows from the
trapezoid quadrature formula that ${\rm Ind}(z)$ given by (2.7) is a discrete
formula of the following continuous function

$\displaystyle\widetilde{{\rm Ind}}(z)=-{\rm
Im}\int_{\Gamma_{r}}\int_{\Gamma_{s}}\kappa(x_{r})\kappa(x_{s})G(z,x_{s})G(z,x_{r})\overline{V(x_{r},x_{s})}{\rm
d}s(x_{s}){\rm d}s(x_{r}),\quad z\in S.$ (2.8)

In the remaining part of this section, we restrict ourselves to demonstrate
that the indicator $\widetilde{{\rm Ind}}(z)$ enjoys the nice feature that it
always peaks on the local perturbation $B$ and on the embedded obstacle $D$.
To this end, we first introduce the following modified Helmholtz-Kirchhoff
identity associated with the Green function $G$.

######  Lemma 2.1.

Let $G$ be the background Green’s function defined by (2.4). Then for any
$x,z\in B_{R}\setminus\Gamma_{0}$, we have

$\displaystyle\int_{\partial B_{R}}\left(\overline{G(\xi,x)}\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-\frac{\partial\overline{G(\xi,x)}}{\partial\nu(\xi)}G(\xi,z)\right){\rm
d}s(\xi)=2{\rm i}{\rm Im}G(x,z).$ (2.9)

###### Proof.

For any $x,z\in B_{R}\setminus\Gamma_{0}$, we choose a sufficiently small
$\varepsilon>0$ such that the circles $B_{\varepsilon}(x)$,
$B_{\varepsilon}(z)$ with $x,z$ as the center and $\varepsilon$ as the radius
contain in the domain $B_{R}$ and
$B_{\varepsilon}(x)\cap\Gamma_{0}=\emptyset$,
$B_{\varepsilon}(z)\cap\Gamma_{0}=\emptyset$. Note that $G$ and its normal
derivative are continuous across the interface $\Gamma_{0}$, thus a direct
application of the Green theorem to $\overline{G(\cdot,x)}$ and $G(\cdot,z)$
in the domain $(B_{R}^{+}\cup B_{R}^{-})\setminus(B_{\varepsilon}(x)\cup
B_{\varepsilon}(z))$ yields

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\int_{(B_{R}^{+}\cup
B_{R}^{-})\setminus(B_{\varepsilon}(x)\cup
B_{\varepsilon}(z))}\left(\overline{G(\xi,x)}\Delta
G(\xi,z)-\Delta\overline{G(\xi,x)}G(\xi,z)\right){\rm d}\xi$ $\displaystyle=$
$\displaystyle\int_{\partial B_{R}\cup\partial B_{\varepsilon}(x)\cup\partial
B_{\varepsilon}(z)}\left(\overline{G(\xi,x)}\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-\frac{\partial\overline{G(\xi,x)}}{\partial\nu(\xi)}G(\xi,z)\right){\rm
d}s(\xi)$ $\displaystyle:=$ $\displaystyle I_{1}+I_{2}+I_{3}$

where $B_{R}^{+}:=B_{R}\cap{\mathbb{R}}^{2}_{+}$,
$B_{R}^{-}:=B_{R}\cap{\mathbb{R}}^{2}_{-}$, $\nu(\xi)$ denotes the unit
outward normal to $\partial B_{R}$ when $\xi\in\partial B_{R}$, $\nu(\xi)$
denotes the unit normal to $\partial B_{\varepsilon}(y)$ into the interior of
$B_{\varepsilon}(y)$ for $y\in\\{x,z\\}$. For the case
$x\in{\mathbb{R}}^{2}_{+}$, notice that
$G(\xi,x):=\Phi_{\kappa_{1}}(\xi,x)+G^{s}(\xi,x)$, we have

$\displaystyle I_{2}$ $\displaystyle=$ $\displaystyle\int_{\partial
B_{\varepsilon}(x)}\left(\overline{\Phi_{\kappa_{1}}(\xi,x)}\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-\frac{\partial\overline{\Phi_{\kappa_{1}}(\xi,x)}}{\partial\nu(\xi)}G(\xi,z)\right){\rm
d}s(\xi)$ $\displaystyle+\int_{\partial
B_{\varepsilon}(x)}\left(\overline{G^{s}(\xi,x)}\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-\frac{\partial\overline{G^{s}(\xi,x)}}{\partial\nu(\xi)}G(\xi,z)\right){\rm
d}s(\xi)$ $\displaystyle\to$ $\displaystyle-G(x,z),\quad{\rm
as}\;\;\varepsilon\to 0.$

For the case $x\in{\mathbb{R}}^{2}_{-}$, by a similar argument, we obtain
$I_{2}\to-G(x,z)$ as $\varepsilon\to 0$. Similarly, we have

$\displaystyle\lim_{\varepsilon\to 0}I_{3}=\overline{G(z,x)}.$

Thus, we obtain

$\displaystyle\int_{\partial B_{R}}\left(\overline{G(\xi,x)}\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-\frac{\partial\overline{G(\xi,x)}}{\partial\nu(\xi)}G(\xi,z)\right){\rm
d}s(\xi)=G(x,z)-\overline{G(z,x)}=2{\rm i}{\rm Im}G(x,z).$

Here, we use the reciprocity $G(x,z)=G(z,x)$ for $x,z\in
B_{R}\setminus\Gamma_{0}$, which can be proven by a direct application of the
Green theorem and the Sommerfeld radiation condition. The proof is completed.
∎

With the aid of the modified Helmholtz-Kirchhoff identity (2.9), we can obtain
the following lemma which plays an important role in the analysis of
$\widetilde{{\rm Ind}}(z)$.

######  Lemma 2.2.

For any $x,z\in S$, we have

$\displaystyle\int_{\partial B_{R}}\kappa(\xi)\overline{G(x,\xi)}G(\xi,z){\rm
d}s(\xi)={\rm Im}G(x,z)+\zeta(x,z)$

where $\zeta(x,z)$ satisfies

$\displaystyle|\zeta(x,z)|+|\nabla_{x}\zeta(x,z)|\leq CR^{-1}$

uniformly for any $x,z\in S$.

###### Proof.

For $x,z\in S$, using Lemma 2.1, we have

$\displaystyle 2{\rm i}{\rm Im}G(x,z)$ $\displaystyle=$
$\displaystyle\int_{\partial B_{R}}\left(\overline{G(\xi,x)}\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-\frac{\partial\overline{G(\xi,x)}}{\partial\nu(\xi)}G(\xi,z)\right){\rm
d}s(\xi)$ $\displaystyle=$ $\displaystyle\int_{\partial
B_{R}}\Bigg{\\{}\overline{G(\xi,x)}\left[\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-{\rm i}\kappa(\xi)G(\xi,z)\right]{\rm d}s(\xi)$
$\displaystyle-G(\xi,z)\left[\frac{\overline{\partial
G(\xi,x)}}{\partial\nu(\xi)}+{\rm
i}\kappa(\xi)\overline{G(\xi,x)}\right]\Bigg{\\}}{\rm d}s(\xi)$
$\displaystyle+2{\rm i}\int_{\partial
B_{R}}\kappa(\xi)\overline{G(\xi,x)}G(\xi,z){\rm d}s(\xi),$

which implies

$\displaystyle\int_{\partial B_{R}}\kappa(\xi)\overline{G(x,\xi)}G(\xi,z){\rm
d}s(\xi)={\rm Im}G(x,z)+\zeta(x,z)\quad{\rm for}\;\;{\forall x,z\in S},$

with

$\displaystyle\zeta(x,z)$ $\displaystyle=$ $\displaystyle\frac{\rm
i}{2}\int_{\partial B_{R}}\Bigg{\\{}\overline{G(\xi,x)}\left[\frac{\partial
G(\xi,z)}{\partial\nu(\xi)}-{\rm i}\kappa(\xi)G(\xi,z)\right]{\rm d}s(\xi)$
$\displaystyle\qquad\;\;\;\;\;-G(\xi,z)\left[\frac{\overline{\partial
G(\xi,x)}}{\partial\nu(\xi)}+{\rm
i}\kappa(\xi)\overline{G(\xi,x)}\right]\Bigg{\\}}{\rm d}s(\xi).$

For $y\in\\{x,z\\}$, since

$\displaystyle G(\xi,y)=O(|\xi|^{-\frac{1}{2}}),\qquad\frac{\partial
G(\xi,y)}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)G(\xi,y)=O(|\xi|^{-\frac{3}{2}})$

we conclude that

$\displaystyle|\zeta(x,z)|\leq CR^{-1}$

uniformly for $x,z\in S$. Moreover, due to

$\displaystyle\frac{\partial G(\xi,x)}{\partial
x_{j}}=O(|\xi|^{-\frac{1}{2}}),\qquad\frac{\partial}{\partial
x_{j}}\left[\frac{\partial G(\xi,x)}{\partial\nu(\xi)}-{\rm i}\kappa
G(\xi,x)\right]=O(|\xi|^{-\frac{3}{2}})$

for $j=1,2$, we obtain

$\displaystyle|\nabla_{x}\zeta(x,z)|\leq CR^{-1}$

uniformly for $x,z\in S$. The proof is finished. ∎

######  Lemma 2.3.

Let $V(x,x_{s})$ be the solution of Problem (2.6). Then we have

$\displaystyle V(x,x_{s})=\int_{\partial D}\left[\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}G(\xi,x)-\frac{\partial
G(\xi,x)}{\partial\nu(\xi)}V(\xi,x_{s})\right]{\rm
d}s(\xi)-\beta\int_{B}\chi(\xi)G(\xi,x)u(\xi,x_{s}){\rm d}\xi$

for $x\in{\mathbb{R}}^{2}\setminus(\Gamma\cup\overline{D})$, where $\nu(\xi)$
denotes the unit normal at $\xi\in\partial D$ which directs into the interior
of $D$.

###### Proof.

For any $x\in{\mathbb{R}}^{2}\setminus(\Gamma\cup\overline{D})$, we choose a
sufficiently small $\varepsilon>0$ and a sufficiently large $\rho>0$ such that
the disc $B_{\varepsilon}(x)$ with $x$ as the center and $\varepsilon$ as the
radius contains in the domain
$\Omega_{\rho}:=B_{\rho}\setminus(\Gamma\cup\overline{D})$, where $B_{\rho}$
stands for the disc with the origin as the center and $\rho$ as the radius.
Using the Green’s theorem to the functions $V(\xi,x_{s})$ and $G(\xi,x)$ in
the domain $\Omega_{\rho}\setminus B_{\varepsilon}(x)$, we have

$\displaystyle\int_{\Omega_{\rho}\setminus B_{\varepsilon}(x)}\left[\Delta
V(\xi,x_{s})G(\xi,x)-\Delta G(\xi,x)V(\xi,x_{s})\right]{\rm d}\xi$
$\displaystyle=\int_{\partial B_{\rho}\cup\partial D\cup\partial
B_{\varepsilon}(x)}\left[\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}G(\xi,x)-\frac{\partial
G(\xi,x)}{\partial\nu(\xi)}V(\xi,x_{s})\right]{\rm d}s(\xi)$
$\displaystyle:=J_{1}+J_{2}+J_{3},$ (2.10)

where $\nu(\xi)$ denotes the unit normal which directs into the interior of
$B_{\varepsilon}(x)$ when $\xi\in\partial B_{\varepsilon}(x)$, and directs
into the exterior of $B_{\rho}$ when $\xi\in\partial B_{\rho}$. By (2.4) and
(2.6), it is easily seen that

$\displaystyle\lim_{\varepsilon\to 0}\int_{\Omega_{\rho}\setminus
B_{\varepsilon}(x)}\left[\Delta V(\xi,x_{s})G(\xi,x)-\Delta
G(\xi,x)V(\xi,x_{s})\right]{\rm
d}\xi=\beta\int_{B}\chi(\xi)G(\xi,x)u(\xi,x_{s}){\rm d}\xi.$ (2.11)

For the item $J_{1}$, we have

$\displaystyle J_{1}$ $\displaystyle=$ $\displaystyle\int_{\partial
B_{\rho}}\left[\frac{\partial V(\xi,x_{s})}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)V(\xi,x_{s})\right]G(\xi,x){\rm d}s(\xi)$
$\displaystyle-\int_{\partial B_{\rho}}\left[\frac{\partial
G(\xi,x)}{\partial\nu(\xi)}-{\rm i}\kappa(\xi)G(\xi,x)\right]V(\xi,x_{s}){\rm
d}s(\xi)$ $\displaystyle:=$ $\displaystyle J_{11}-J_{12}.$

Since

$\displaystyle G(\xi,x)=O(|\xi|^{-\frac{1}{2}}),\qquad\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)V(\xi,x_{s})=O(|\xi|^{-\frac{3}{2}})$

it follows from the Cauchy-Schwarz inequality that

$\displaystyle\lim_{\rho\to\infty}J_{11}=0.$

For the item $J_{12}$, we claim that

$\displaystyle\int_{\partial B_{\rho}}|V(\xi,x_{s})|^{2}{\rm
d}s(\xi)=O(1)\quad{\rm as}\;\;\rho\to\infty.$ (2.12)

To show this, it is found by the Sommerfeld radiation condition that

$\displaystyle\int_{\partial B_{\rho}}\left[\left|\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}\right|^{2}+\kappa^{2}(\xi)|V(\xi,x_{s})|^{2}+2\kappa(\xi){\rm
Im}\left(V(\xi,x_{s})\frac{\partial\overline{V(\xi,x_{s})}}{\partial\nu(\xi)}\right)\right]{\rm
d}s(\xi)\to 0$ (2.13)

as $\rho\to\infty$. Applying the Green theorem for $V(\xi,x_{s})$ and
$\overline{V(\xi,x_{s})}$ in the domain $B_{\rho}\setminus\overline{D}$
implies that

$\displaystyle\int_{\partial
B_{\rho}}V(\xi,x_{s})\frac{\partial\overline{V(\xi,x_{s})}}{\partial\nu(\xi)}{\rm
d}s(\xi)$ (2.14) $\displaystyle=$
$\displaystyle\int_{B_{\rho}\setminus\overline{D}}\left(|\nabla
V(\xi,x_{s})|^{2}-\kappa^{2}(\xi)|V(\xi,x_{s})|^{2}\right){\rm d}\xi$
$\displaystyle-$ $\displaystyle\int_{\partial
D}V(\xi,x_{s})\frac{\partial\overline{V(\xi,x_{s})}}{\partial\nu(\xi)}{\rm
d}s(\xi)+\beta\int_{B}\chi(\xi)V(\xi,x_{s})\overline{G(\xi,x_{s})}{\rm d}\xi.$

Taking the imaginary part of (2.14) and substituting it to (2.13) gives that

$\displaystyle\lim_{\rho\to\infty}\frac{1}{2}\int_{\partial
B_{\rho}}\left(\frac{1}{\kappa(\xi)}\left|\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}\right|^{2}+\kappa(\xi)|V(\xi,x_{s})|^{2}\right){\rm
d}s(\xi)$ $\displaystyle={\rm Im}\left\\{\int_{\partial
D}V(\xi,x_{s})\frac{\partial\overline{V(\xi,x_{s})}}{\partial\nu(\xi)}{\rm
d}s(\xi)-\beta\int_{B}\chi(\xi)V(\xi,x_{s})\overline{G(\xi,x_{s})}{\rm
d}\xi\right\\}$

which implies that the claim (2.12) holds true. Since

$\displaystyle\frac{\partial G(\xi,x)}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)G(\xi,x)=O(|\xi|^{-\frac{3}{2}}),$ (2.15)

it follows from (2.12), (2.15), and the Cauchy-Schwarz inequality that
$\lim\limits_{\rho\to\infty}J_{12}=0$. Hence, we obtain

$\displaystyle\lim_{\rho\to\infty}J_{1}=0.$ (2.16)

For the item $J_{3}$, if $x\in{\mathbb{R}}^{2}_{+}$, using
$G(\xi,x)=\Phi_{\kappa_{1}}(\xi,x)+G^{s}(\xi,x)$, we have

$\displaystyle J_{3}$ $\displaystyle=$ $\displaystyle\int_{\partial
B_{\varepsilon}(x)}\left(\Phi_{\kappa_{1}}(\xi,x)\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}-\frac{\partial\Phi_{\kappa_{1}}(\xi,x)}{\partial\nu(\xi)}V(\xi,x_{s})\right){\rm
d}s(\xi)$ $\displaystyle+\int_{\partial
B_{\varepsilon}(x)}\left(G^{s}(\xi,x)\frac{\partial
V(\xi,x_{s})}{\partial\nu(\xi)}-\frac{\partial
G^{s}(\xi,x)}{\partial\nu(\xi)}V(\xi,x_{s})\right){\rm d}s(\xi)$
$\displaystyle:=$ $\displaystyle J_{31}+J_{32}.$

A direct calculation, using the mean value theorem, shows that

$\displaystyle\lim_{\varepsilon\to 0}J_{31}=-V(x,x_{s}).$

For the item $J_{32}$, if $x\in{\mathbb{R}}^{2}_{+}\cap\Omega_{1}$, by the
Green theorem, (2.4) and (2.6), we have

$\displaystyle J_{32}=\int_{B_{\varepsilon}(x)}\left(\Delta
V(\xi,x_{s})G^{s}(\xi,x)-\Delta G^{s}(\xi,x)V(\xi,x)\right){\rm d}\xi=0.$

If $x\in B_{1}$, combining the Green theorem, (2.4) and (2.6) shows that

$\displaystyle J_{32}$ $\displaystyle=$
$\displaystyle\int_{B_{\varepsilon}(x)}\left(\Delta
V(\xi,x_{s})G^{s}(\xi,x)-\Delta G^{s}(\xi,x)V(\xi,x)\right){\rm d}\xi$
$\displaystyle=$
$\displaystyle\beta\int_{B_{\varepsilon}(x)}G^{s}(\xi,x)u(\xi,x_{s}){\rm
d}\xi\to 0\quad{\rm as}\;\;\varepsilon\to 0.$

If $x\in{\mathbb{R}}^{2}_{-}$, by a similar argument, we can obtain the same
result. Thus, we conclude that

$\displaystyle\lim_{\varepsilon\to 0}J_{3}=-V(x,x_{s}).$ (2.17)

It is found by (2.10), (2.11), (2.16), and (2.17) that the Lemma holds true.
The proof is completed. ∎

With the help of Lemma 2.1, Lemma 2.2, and Lemma 2.3, we are in position to
present the main result of this paper, which focuses on the resolution
analysis of the RTM approach for simultaneously reconstruction of the locally
rough interface and the embedded obstacle.

######  Theorem 2.4.

For any $z\in S$, let $\psi(\xi,z)$ be the solution of

$\left\\{\begin{array}[]{lll}\Delta\psi(\xi,z)+\kappa^{2}(\xi)\psi(\xi,z)=\beta\chi(\xi){\rm
Im}G(\xi,z)&{\rm in}\;\;{\mathbb{R}}^{2}\setminus\overline{D},\\\\[5.69054pt]
\psi(\xi,z)=-{\rm Im}G(\xi,z)&{\rm on}\;\;\partial D,\\\\[5.69054pt]
\lim\limits_{|\xi|\rightarrow\infty}|\xi|^{\frac{1}{2}}\left(\partial_{|\xi|}\psi(\xi,z)-{\rm
i}\kappa(\xi)\psi(\xi,z)\right)=0,\end{array}\right.$ (2.18)

and $\psi^{\infty}(\hat{\xi},z)$ be the corresponding far-field pattern. Then
we have

$\displaystyle\widetilde{{\rm
Ind}}(z)=\int_{{\mathbb{S}}}\kappa(\hat{\xi})|\psi^{\infty}(\hat{\xi},z)|^{2}{\rm
d}s(\hat{\xi})+\eta(z),\qquad\forall z\in S,$

where $\|\eta(z)\|_{L^{\infty}(S)}\leq CR^{-1}$ with some constant $C$
depending on $B$ and $D$.

###### Proof.

Define

$\displaystyle\widetilde{W}(z,x_{s})=\int_{\Gamma_{r}}\kappa(x_{r})G(z,x_{r})\overline{V(x_{r},x_{s})}{\rm
d}s(x_{r}),$ (2.19)

substituting the Green formula presented in Lemma 2.3 into
$\widetilde{W}(z,x_{s})$ gives that

$\displaystyle\widetilde{W}(z,x_{s})=\int_{\partial
D}\left[\frac{\partial\overline{V(\xi,x_{s})}}{\partial\nu(\xi)}\gamma(\xi,z)-\frac{\partial\gamma(\xi,z)}{\partial\nu(\xi)}\overline{V(\xi,x_{s})}\right]{\rm
d}s(\xi)-\beta\int_{B}\chi(\xi)\overline{u(\xi,x_{s})}\gamma(\xi,z){\rm d}\xi$

with

$\displaystyle\gamma(\xi,z):=\int_{\Gamma_{r}}\kappa(x_{r})G(z,x_{r})\overline{G(\xi,x_{r})}{\rm
d}s(x_{r})={\rm Im}G(\xi,z)+\zeta(\xi,z)$ (2.20)

where we use Lemma 2.2. According to (2.8) and (2.19), we obtain

$\displaystyle\widetilde{\rm Ind}(z)$ $\displaystyle=$ $\displaystyle-{\rm
Im}\int_{\Gamma_{s}}\kappa(x_{s})G(z,x_{s})\widetilde{W}(z,x_{s}){\rm
d}s(x_{s})$ (2.21) $\displaystyle=$ $\displaystyle-{\rm Im}\int_{\partial
D}\left[\frac{\partial\phi_{1}(\xi,z)}{\partial\nu(\xi)}\gamma(\xi,z)-\frac{\partial\gamma(\xi,z)}{\partial\nu(\xi)}\phi_{1}(\xi,z)\right]{\rm
d}s(\xi)$ $\displaystyle+{\rm
Im}\beta\int_{B}\chi(\xi)\phi_{2}(\xi,z)\gamma(\xi,z){\rm d}\xi$

where $\phi_{1}(\xi,z)$ and $\phi_{2}(\xi,z)$ are defined by

$\displaystyle\phi_{1}(\xi,z):=\int_{\Gamma_{s}}\kappa(x_{s})G(z,x_{s})\overline{V(\xi,x_{s})}{\rm
d}s(x_{s})$
$\displaystyle\phi_{2}(\xi,z):=\int_{\Gamma_{s}}\kappa(x_{s})G(z,x_{s})\overline{u(\xi,x_{s})}{\rm
d}s(x_{s}).$

By the identity $u(\xi,x_{s})=V(\xi,x_{s})+G(\xi,x_{s})$ and Lemma 2.2, it is
easily found that

$\displaystyle\phi_{2}(\xi,z)=\phi_{1}(\xi,z)+{\rm
Im}G(\xi,z)+\overline{\zeta(\xi,z)}.$ (2.22)

Since $V(\xi,x_{s})$ satisfies Problem (2.6), it is easily checked that
$\overline{\phi_{1}(\xi,z)}$ solves

$\left\\{\begin{array}[]{lll}\Delta\overline{\phi_{1}(\xi,z)}+\kappa^{2}(\xi)\overline{\phi_{1}(\xi,z)}=\beta\chi(\xi)\left[{\rm
Im}G(\xi,z)+\zeta(\xi,z)\right]&{\rm
in}\;\;{\mathbb{R}}^{2}\setminus\overline{D},\\\\[8.53581pt]
\overline{\phi_{1}(\xi,z)}=-\left[{\rm Im}G(\xi,z)+\zeta(\xi,z)\right]&{\rm
on}\;\;\partial D,\\\\[8.53581pt]
\lim\limits_{|\xi|\rightarrow\infty}|\xi|^{\frac{1}{2}}\left(\partial_{|\xi|}\overline{\phi_{1}(\xi,z)}-{\rm
i}\kappa(\xi)\overline{\phi_{1}(\xi,z)}\right)=0,\end{array}\right.$ (2.23)

where we use Lemma 2.2 to obtain the right hand term. Let $\psi(\xi,z)$ and
$\varphi(\xi,z)$ solve Problem (2.23) expect the right hand term are replaced
by $(\beta\chi(\xi){\rm Im}G(\xi,z),-{\rm Im}G(\xi,z))$ and
$(\beta\chi(\xi)\zeta(\xi,z),-\zeta(\xi,z))$, respectively. Then by the
linearity we have

$\displaystyle\overline{\phi_{1}(\xi,z)}=\psi(\xi,z)+\varphi(\xi,z).$ (2.24)

Substituting (2.20), (2.22) and (2.24) into (2.21) implies that

$\displaystyle\widetilde{\rm Ind}(z)$ $\displaystyle=-{\rm Im}\int_{\partial
D}\bigg{\\{}\left[{\rm
Im}G(\xi,z)+\zeta(\xi,z)\right]\frac{\partial}{\partial\nu(\xi)}\left[\overline{\psi(\xi,z)}+\overline{\varphi(\xi,z)}\right]$
$\displaystyle\quad-\left[\overline{\psi(\xi,z)}+\overline{\varphi(\xi,z)}\right]\frac{\partial}{\partial\nu(\xi)}\left[{\rm
Im}G(\xi,z)+\zeta(\xi,z)\right]\bigg{\\}}{\rm d}s(\xi)$
$\displaystyle\quad+{\rm
Im}\beta\int_{B}\chi(\xi)\left[\overline{\psi(\xi,z)}+\overline{\varphi(\xi,z)}+{\rm
Im}G(\xi,z)+\overline{\zeta(\xi,z)}\right]\left[{\rm
Im}G(\xi,z)+\zeta(\xi,z)\right]{\rm d}\xi$ $\displaystyle=-{\rm
Im}\int_{\partial
D}\left[\frac{\partial\overline{\psi(\xi,z)}}{\partial\nu(\xi)}{\rm
Im}G(\xi,z)-\frac{\partial{\rm
Im}G(\xi,z)}{\partial\nu(\xi)}\overline{\psi(\xi,z)}\right]{\rm d}s(\xi)$
$\displaystyle\quad+{\rm Im}\beta\int_{B}\chi(\xi)\overline{\psi(\xi,z)}{\rm
Im}G(\xi,z){\rm d}\xi+\eta(z)$ $\displaystyle={\rm Im}\int_{\partial
D}\frac{\partial\overline{\psi(\xi,z)}}{\partial\nu(\xi)}\psi(\xi,z){\rm
d}s(\xi)+{\rm Im}\beta\int_{B}\chi(\xi)\overline{\psi(\xi,z)}{\rm
Im}G(\xi,z){\rm d}\xi+\eta(z),$

where $\eta(z)$ is defined by

$\displaystyle\eta(z)$ $\displaystyle:=$ $\displaystyle-{\rm Im}\int_{\partial
D}\left[\frac{\partial\overline{\psi(\xi,z)}}{\partial\nu(\xi)}\zeta(\xi,z)-\frac{\partial\zeta(\xi,z)}{\partial\nu(\xi)}\overline{\psi(\xi,z)}\right]{\rm
d}s(\xi)$ $\displaystyle-{\rm Im}\int_{\partial
D}\left[\frac{\partial\overline{\varphi(\xi,z)}}{\partial\nu(\xi)}{\rm
Im}G(\xi,z)-\frac{\partial{\rm
Im}G(\xi,z)}{\partial\nu(\xi)}\overline{\varphi(\xi,z)}\right]{\rm d}s(\xi)$
$\displaystyle-{\rm Im}\int_{\partial
D}\left[\frac{\partial\overline{\varphi(\xi,z)}}{\partial\nu(\xi)}\zeta(\xi,z)-\frac{\partial\zeta(\xi,z)}{\partial\nu(\xi)}\overline{\varphi(\xi,z)}\right]{\rm
d}s(\xi)$ $\displaystyle+{\rm
Im}\beta\int_{B}\chi(\xi)\left[\overline{\psi(\xi,z)}\zeta(\xi,z)+\overline{\varphi(\xi,z)}{\rm
Im}G(\xi,z)+\overline{\varphi(\xi,z)}\zeta(\xi,z)\right]{\rm d}s(\xi).$

We choose a sufficiently large $\rho>0$ such that $B\cup D\subset B_{\rho}$
and apply the Green theorem for $\psi(\xi,z)$ and $\overline{\psi(\xi,z)}$ in
the domain $B_{\rho}\setminus\overline{D}$ to get

$\displaystyle\int_{B_{\rho}\setminus\overline{D}}\Delta\psi(\xi,z)\overline{\psi(\xi,z)}{\rm
d}\xi$
$\displaystyle=\int_{B_{\rho}\setminus\overline{D}}\left[\beta\chi(\xi){\rm
Im}G(\xi,z)-\kappa^{2}(\xi)\psi(\xi,z)\right]\overline{\psi(\xi,z)}{\rm d}\xi$
$\displaystyle=\int_{\partial
B_{\rho}}\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}\overline{\psi(\xi,z)}{\rm
d}s(\xi)+\int_{\partial
D}\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}\overline{\psi(\xi,z)}{\rm
d}s(\xi)-\int_{B_{\rho}\setminus\overline{D}}|\nabla\psi(\xi,z)|^{2}{\rm
d}\xi.$

Inserting the imaginary part on both sides of the above equation gives that

$\displaystyle{\rm Im}\left\\{\int_{\partial
D}\frac{\partial\overline{\psi(\xi,z)}}{\partial\nu(\xi)}\psi(\xi,z){\rm
d}s(\xi)+\beta\int_{B}\chi(\xi)\overline{\psi(\xi,z)}{\rm Im}G(\xi,z){\rm
d}\xi\right\\}={\rm Im}\int_{\partial
B_{\rho}}\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}\overline{\psi(\xi,z)}{\rm
d}s(\xi).$

Thus, we obtain

$\displaystyle\widetilde{\rm Ind}(z)={\rm Im}\int_{\partial
B_{\rho}}\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}\overline{\psi(\xi,z)}{\rm
d}s(\xi)+\eta(z)$ $\displaystyle={\rm Im}\int_{\partial
B_{\rho}}\left[\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)\psi(\xi,z)\right]\overline{\psi(\xi,z)}{\rm
d}s(\xi)+\int_{\partial B_{\rho}}\kappa(\xi)|\psi(\xi,z)|^{2}{\rm
d}s(\xi)+\eta(z).$

Since

$\displaystyle\psi(\xi,z)=O(|\xi|^{-\frac{1}{2}}),\qquad\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)\psi(\xi,z)=O(|\xi|^{-\frac{3}{2}})$

we have

$\displaystyle\lim_{\rho\to\infty}\int_{\partial
B_{\rho}}\left[\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}-{\rm
i}\kappa(\xi)\psi(\xi,z)\right]\overline{\psi(\xi,z)}{\rm d}s(\xi)=0.$

It follows from the Sommerfeld radiation condition that

$\displaystyle\psi(\xi,z)=\frac{e^{{\rm
i}\kappa|\xi|}}{|\xi|^{\frac{1}{2}}}\left\\{\psi^{\infty}(\hat{\xi},z)+O\left(\frac{1}{|\xi|}\right)\right\\}\quad{\rm
for}\;|\xi|\to\infty$

which leads to

$\displaystyle\lim_{\rho\to\infty}\int_{\partial
B_{\rho}}\kappa(\xi)|\psi(\xi,z)|^{2}{\rm
d}s(\xi)=\int_{\mathbb{S}}\kappa(\hat{\xi})|\psi^{\infty}(\hat{\xi},z)|^{2}{\rm
d}s(\hat{\xi}).$

Thus, we conclude that

$\displaystyle\widetilde{{\rm
Ind}}(z)=\int_{{\mathbb{S}}}\kappa(\hat{\xi})|\psi^{\infty}(\hat{\xi},z)|^{2}{\rm
d}s(\hat{\xi})+\eta(z),\qquad\forall z\in S.$

In the remaining part of the proof, we restrict ourselves to the estimate of
$\eta(z)$. For $\xi\in S$, $z\in S$, due to Lemma 2.2, we have

$\displaystyle|\zeta(\xi,z)|+|\nabla_{\xi}\zeta(\xi,z)|\leq CR^{-1}.$ (2.25)

Notice that

$\displaystyle{\rm
Im}G(\xi,z)=\left\\{\begin{array}[]{ll}\frac{1}{4}J_{0}(\sigma|\xi-z|)+{\rm
Im}G^{s}(\xi,z)&{\rm for}\;\;\xi,z\in{\mathbb{R}}^{2}_{+}\;{\rm
or}\;\xi,z\in{\mathbb{R}}^{2}_{-},\\\\[8.53581pt] {\rm Im}G^{s}(\xi,z)&{\rm
for}\;\;\xi\in{\mathbb{R}}^{2}_{+},z\in{\mathbb{R}}^{2}_{-}\;{\rm
or}\;\xi\in{\mathbb{R}}^{2}_{-},z\in{\mathbb{R}}^{2}_{+}\end{array}\right.$

where $J_{0}$ denotes the Bessel function of order zero, $\sigma=\kappa_{1}$
for $\xi,z\in{\mathbb{R}}^{2}_{+}$ and $\sigma=\kappa_{2}$ for
$\xi,z\in{\mathbb{R}}^{2}_{-}$. By the smoothness of $G^{s}(\xi,z)$ and
$J_{0}$, we obtain

$\displaystyle\left|{\rm Im}G(\xi,z)\right|+\left|\frac{\partial{\rm
Im}G(\xi,z)}{\partial\nu(\xi)}\right|\leq C$ (2.27)

for $\xi\in B\cup D$ and $z\in S$. Since $\psi(\xi,z)$ and $\varphi(\xi,z)$
solve Problem (2.23) with the boundary data $(\beta\chi(\xi){\rm
Im}G(\xi,z),-{\rm Im}G(\xi,z))$ and
$(\beta\chi(\xi)\zeta(\xi,z),-\zeta(\xi,z))$, respectively, a straightforward
application of Theorem 3.1 in [37] leads to

$\displaystyle\|\psi(\xi,z)\|_{H^{\frac{1}{2}}(\partial
D)}+\left\|\frac{\partial\psi(\xi,z)}{\partial\nu(\xi)}\right\|_{H^{-\frac{1}{2}}(\partial
D)}+\|\psi(\xi,z)\|_{L^{2}(B)}$ $\displaystyle\lesssim\|{\rm
Im}G(\xi,z)\|_{L^{2}(B)}+\|{\rm Im}G(\xi,z)\|_{H^{\frac{1}{2}}(\partial D)}$
$\displaystyle\leq C$ (2.28)

and

$\displaystyle\|\varphi(\xi,z)\|_{H^{\frac{1}{2}}(\partial
D)}+\left\|\frac{\partial\varphi(\xi,z)}{\partial\nu(\xi)}\right\|_{H^{-\frac{1}{2}}(\partial
D)}+\|\varphi(\xi,z)\|_{L^{2}(B)}$
$\displaystyle\lesssim\|\zeta(\xi,z)\|_{L^{2}(B)}+\|\zeta(\xi,z)\|_{H^{\frac{1}{2}}(\partial
D)}$ $\displaystyle\leq CR^{-1},$ (2.29)

where we use the trace theorem in the first step of (2.28) and (2.29). Thus,
with the help of (2.25)-(2.29), we obtain

$\displaystyle\|\eta(z)\|_{L^{\infty}(S)}\leq CR^{-1}$

with $C$ depending on $B$ and $D$. The proof is finished. ∎

##  3 Numerical examples

In this section, we first analyze the behavior of the indicator
$\widetilde{\rm Ind}(z)$ when $z$ is close to $B\cup D$, and then present
several numerical examples to show the effectiveness of the RTM method.

Due to Theorem 2.4, it is found that the behavior of $\widetilde{\rm Ind}(z)$
depends on the function $\psi(\xi,z)$ when the source and measurement radius
$R$ is large enough. Here, the function $\psi(\xi,z)$ solves Problem (2.18)
with boundary data $(\beta\chi(\xi){\rm Im}G(\xi,z),-{\rm Im}G(\xi,z))$.
Notice that

$\displaystyle{\rm
Im}G(\xi,z)=\left\\{\begin{array}[]{ll}\frac{1}{4}J_{0}(\sigma|\xi-z|)+{\rm
Im}G^{s}(\xi,z)&{\rm for}\;\;\xi,z\in{\mathbb{R}}^{2}_{+}\;{\rm
or}\;\xi,z\in{\mathbb{R}}^{2}_{-},\\\\[8.53581pt] {\rm Im}G^{s}(\xi,z)&{\rm
for}\;\;\xi\in{\mathbb{R}}^{2}_{+},z\in{\mathbb{R}}^{2}_{-}\;{\rm
or}\;\xi\in{\mathbb{R}}^{2}_{-},z\in{\mathbb{R}}^{2}_{+}\end{array}\right.$

where $\sigma=\kappa_{1}$ for $\xi,z\in{\mathbb{R}}^{2}_{+}$ and
$\sigma=\kappa_{2}$ for $\xi,z\in{\mathbb{R}}^{2}_{-}$. Observing from (a) of
Figure 2, the function $J_{0}(\lambda|\xi-z|)$ achieves a maximum at $\xi=z$,
thus, we guess that ${\rm Im}G(\xi,z)$ will also achieves a maximum at
$\xi=z$, which is conformed numerically in (b) and (c) of Figure 2. Based on
this observation, we can expect that $\widetilde{\rm Ind}(z)$ will peak on
$B\cup D$.

(a) $J_{0}$

(b) ${\rm Im}G(x,z)$ with $z=(0,0.5)$

(c) ${\rm Im}G(x,z)$ with $z=(0,-0.5)$

Figure 2: The image of functions $J_{0}$, ${\rm Im}G(\xi,z)$ with
$\xi\in[-5,5]\times[-8.95,1.05]$, $\kappa_{1}=10$, and $\kappa_{2}=5$.

In all examples, if not stated otherwise, we assume that the wavenumber
$\kappa_{1}=10$, $\kappa_{2}=5$, the locally rough interface function $f$ is
supported in the interval $[-5,5]$, the sample domain
$S=[-5,5]\times[-8.95,1.05]$, $N_{s}=N_{r}=1024$, and $R=100$. The synthetic
data $u^{s}(x_{r},x_{s})$ and $G^{s}(x_{r},x_{s})$ are generated by solving
Problems (2.3) and (2.4) through the Nyström method [33]. For some relative
error $\tau>0$, we inject some noise into the data by defining

$\displaystyle
u_{\tau}^{s}(x_{r},x_{s})=u^{s}(x_{r},x_{s})+\tau\frac{\lambda}{\|\lambda\|_{2}}\|u^{s}(x_{r},x_{s})\|_{2}$

where $\lambda=\lambda_{1}+{\rm i}\lambda_{2}$ is complex-valued with
$\lambda_{1}$ and $\lambda_{2}$ consisting of random numbers obeying standard
normal distribution $N(0,1)$.

Example 1. In this example, we consider two simple cases and exam the RTM
method at different noise level. The first one (see (a) in Figure 3) is
related to a planar interface $\Gamma=\Gamma_{0}$ and a circle given by

$\displaystyle
x(\theta)=(0.5\cos(\theta),-4+0.5\sin(\theta)),\qquad\theta\in[0,2\pi).$

The second case (see (d) in Figure 3) is just related to a locally rough
interface $\Gamma$ which is described by

$\displaystyle f(t)=\Omega_{3}(2t+4)-0.6\Omega_{3}(2t-5)$

with $\Omega_{3}(\cdot)\in C_{0}^{2}(\mathbb{R})$ being a cubic $B$-spline
function given by

$\displaystyle\Omega_{3}(t)=\left\\{\begin{aligned}
\frac{1}{2}|t|^{3}-t^{2}+\frac{2}{3}&\qquad\qquad\textrm{for}\;|t|\leq 1,\\\
-\frac{1}{6}|t|^{3}+t^{2}-2|t|+\frac{4}{3}&\qquad\qquad\textrm{for}\;1<|t|<2\,,\\\
0&\qquad\qquad\textrm{for}\;|t|\geq 2.\end{aligned}\right.$

The reconstruction results are presented in Figure 3, which shows that the RTM
method can provide a satisfied reconstruction for these two simple cases, even
for $10\%$ noisy data.

Example 2. In this example, we test the dependence of the RTM algorithm on the
relative position and distance between the local perturbation $B$ and the
embedded obstacle $D$. The locally rough interface $\Gamma$ and the rounded
square shaped obstacle $D$ are parameterized by

$\displaystyle
f(t)=\left(0.6e^{-6(t+3)^{2}}+0.5e^{-7t^{2}}+0.5e^{-8(t-3)^{2}}\right)\cdot
f_{0}(t)\quad t\in{\mathbb{R}},$ $\displaystyle
x(\theta)=(3+0.3(\cos^{3}(\theta)+\cos(\theta)),-6+0.3(\sin^{3}(\theta)+\sin(\theta))),\;\theta\in[0,2\pi),$

where $f_{0}(t)\in C_{0}^{\infty}({\mathbb{R}})$ is a cut-off function defined
by

$\displaystyle f_{0}(t)=\left\\{\begin{aligned}
1&\qquad\qquad\textrm{for}\;|t|\leq 4,\\\
\;\left(1+e^{\frac{1}{5-|t|}+\frac{1}{4-|t|}}\right)^{-1}&\qquad\qquad\textrm{for}\;4<|t|<5\,,\\\
0&\qquad\qquad\textrm{for}\;|t|\geq 5;\end{aligned}\right.$

see (a) in Figure 4. Next, we fix the locally rough interface $\Gamma$ and
move the embedded obstacle $D$ up by four units and five units, see (b) and
(c) in Figure 4. The numerical results from $10\%$ noisy data are illustrated
in (d), (e), (f) of Figure 4. It is readily seen from Figure 4 that the
quality of the reconstruction will depend on the relative position of $B$ and
$D$. We guess that a possible reason is due to the strong multiple scattering
between $B$ and $D$ when $D$ is close to $\Gamma$.

Example 3. In the last example, as an attempt, we consider a piecewise
continuous locally rough interface and exam the RTM algorithm at different
radiuses $R$. The locally rough interface is described by a piecewise constant
defined by

$f(t)=\left\\{\begin{array}[]{lll}0.2\qquad\;\;|t|\leq 1,\\\\[2.84526pt]
0.3\qquad\;\;3\leq|t|\leq 4,\\\\[2.84526pt] 0\qquad\;\;{\rm
others}.\end{array}\right.$

The embedded obstacle is rounded triangle shaped given by

$\displaystyle
x(\theta)=(-3+\left(0.5+0.1\cos(3\theta)\right)\cos(\theta),-6+\left(0.5+0.1\cos(3\theta)\right)\sin(\theta)),\;\theta\in[0,2\pi),$

see (a) in Figure 5. The radius is set to be $R=20$ and $R=100$, and the noise
level is $10\%$. We present the reconstructions in (b), (c) of Figure 5, which
shows that the RTM algorithm is able to provide high reconstruction quality at
these measurement radiuses.

It is observed from Figure 3, Figure 4, and Figure 5 that the RTM approach
proposed in Theorem 2.4 can provide accurate and stable reconstructions of the
locally rough interface as well as the embedded obstacles for a variety of
interfaces and obstacles. The image quality depends on the interaction of the
locally rough interface and the embedded obstacle. Moreover, as shown in
Figure 3, Figure 4, and Figure 5, the RTM method is robust to noise.

(a) Physical configuration

(b) No noise

(c) $10\%$ noise

(d) Physical configuration

(e) No noise

(f) $10\%$ noise

Figure 3: Reconstructions of the two simple cases in Example 1 from data at
different noise level.

(a) Physical configuration

(b) Physical configuration

(c) Physical configuration

(d) $10\%$ noise

(e) $10\%$ noise

(f) $10\%$ noise

Figure 4: Reconstructions results for different relative position between the
locally rough interface and the embedded obstacle.

(a) Physical configuration

(b) $R=20$, $10\%$ noise

(c) $R=100$, $10\%$ noise

Figure 5: Reconstructions of the locally rough interface and the object given
in Example 3 for different measurement radiuses.

##  4 Conclusion

In this paper, we proposed an extended RTM method to simultaneously recover a
locally rough interface and an embedded object in the lower half-space from
near-field measurements. The main idea is based on constructing a modified
Helmholtz-Kirchhoff identity associated with the planar interface. Numerical
experiments showed that the novel RTM algorithm can provide an accurate and
stable reconstruction for a large number of locally rough interfaces and
embedded obstacles, even for piecewise continuous interfaces. Furthermore, the
reconstructions can be regarded as a good initial guess for an iterative
algorithm to obtain a more accurate result. As far as we know, this is the
first RTM method to simultaneously recover an unbounded scatterer and a
bounded scatterer. It is easily seen that the RTM method proposed in Theorem
2.4 depends crucially on the priori information that the interface is locally
perturbed. However, if the rough surface is non-local, it is unclear how to
establish the related Helmholtz-Kirchhoff identity and then extend this method
to reconstruct the non-local rough surface. We hope to report the progress on
this topic in the future.

Acknowledgments. This work was supported by the NNSF of China grants No.
12171057, 12122114, and Education Department of Hunan Province No. 21B0299.

## References

  * [1] S. Arridge, Optical tomography in medical imaging, Inverse Probl., 15(1999), R41-R93.
  * [2] H. Ammari, E. Iakovleva, and D. Lesselier, A MUSIC algorithm for locating small inclusions buried in a half-space from the scattering amplitude at a fixed frequency, SIAM Multiscale Model. Simul., 3(2005), 597-628.
  * [3] A. Berkhout, Seismic migration: imaging of acoustic energy by wave field extrapolation, Elsevier, 1984.
  * [4] N. Bleistein, J. Cohen, and J. Stockwell, Mathematics of multidimensional seismic imaging, migration, and inversion, Springer, 2001.
  * [5] G. Bao and P. Li, Near-field imaging of infinite rough surfaces, SIAM J. Appl. Math., 73 (2013), 2162-2187.
  * [6] G. Bao and P. Li, Near-field imaging of infinite rough surfaces in dielectric media, SIAM J. Imaging Sci., 7 (2014), 867-899.
  * [7] G. Bao and J. Lin, Imaging of local surface displacement on an infinite ground plane: the multiple frequency case, SIAM J. Appl. Math., 71 (2011), 1733-1752.
  * [8] G. Bao and J. Lin, Imaging of local surface displacement on an infinite ground plane, Inverse Probl. Imag., 7(2013), 377-396.
  * [9] C. Burkard and R. Potthast, A multi-section approach for rough surface reconstruction via the Kirsch-Kress scheme, Inverse Probl., 26 (2010), 045007.
  * [10] J F. Claerbout, Imaging the Earth’s interior, Oxford: Blackwell Scientific Publication, 1985.
  * [11] J. Chen, Z. Chen, and G. Huang, Reverse time migration for extended obstacles: acoustic waves, Inverse Probl., 29 (2013), 085005.
  * [12] J. Chen, Z. Chen, and G. Huang, Reverse time migration for extended obstacles: electromagnetic waves, Inverse Probl., 29 (2013), 085006.
  * [13] Z. Chen and G. Huang, Reverse time migration for extended obstacles: elastic waves, Sci. Sin. Math., 45(2015), 1103-1114.
  * [14] Z. Chen and G. Huang, Reverse time migration for reconstructing extended obstacles in planar acoustic waveguides, Sci. China Math., 58, 1811-1834.
  * [15] Z. Chen and G. Huang, Reverse time migration for reconstructing extended obstacles in the half space, Inverse Probl., 31(2015), 055007.
  * [16] Z. Chen and G. Huang, A direct imaging method for electromagnetic scattering data without phase information, SIAM J. Imaging Sci., 9(2016), 1273-1297.
  * [17] Z. Chen and G. Huang, Phaseless imaging by reverse time migration: acoustic waves. Numer. Math. Theor. Meth. Appl., 10(2017), 1-21.
  * [18] S.N. Chandler-Wilde and J. Elschner, Variational approach in weighted Sobolev spaces to scattering by unbounded rough surfaces, SIAM J. Math. Anal., 42 (2010), 2554-2580.
  * [19] S.N. Chandler-Wilde and A.T. Peplow, A boundary integral equation formulation for the helmholtz equation in a locally perturbed half-plane, ZAMM-J. Appl. Math. Mech., 85(2005), 79-88.
  * [20] S.N. Chandler-Wilde, C.R. Ross, and B. Zhang, Scattering by infinite one-dimensional rough surfaces, Proc. R. Soc. Lond., A455 (1999), 3767-3787.
  * [21] S.N. Chandler-Wilde and B. Zhang, Scattering of electromagnetic waves by rough interfaces and inhomogeneous layers, SIAM J. Math. Anal., 30 (1999), 559-583.
  * [22] L. Chorfi and P. Gaitan, Reconstruction of the interface between two-layered media using far-field measurements, Inverse Probl., 27 (2011), 075001.
  * [23] S.N. Chandler-Wilde and B. Zhang, On the solvability of a class of second kind integral equations on unbounded domains, J. Math. Anal. Appl., 214 (1997), 482-502.
  * [24] S.N. Chandler-Wilde, B. Zhang and C.R. Ross, On the solvability of second kind integral equations on the real line, J. Math. Anal. Appl., 245 (2000), 28-51.
  * [25] M. Ding, J. Li, K. Liu, and J. Yang, Imaging of locally rough surfaces by the linear sampling method with the near-field data, SIAM J. Imaging Sci., 10(2017), 1579-1602.
  * [26] R. Griesmaier, An asymptotic factorization method for inverse electromagnetic scattering in layered media, SIAM J. Appl. Math., 68(2008), 1378-1403.
  * [27] B. Gebauer, M. Hanke, A. Kirsch, W. Muniz, and C. Schneider, A sampling method for detecting buried objects using electromagnetic scattering, Inverse Probl., 21(2005), 2035-2050.
  * [28] A. Lechleiter, Factorization Methods for Photonics and Rough Surfaces, PhD thesis, KIT, Germany, 2008.
  * [29] J. Li, Reverse time migration for inverse obstacle scattering with a generalized impedance boundary condition, Appl. Anal., 101(2022), 48-62.
  * [30] P. Li, Coupling of finite element and boundary integral methods for electromagnetic scattering in a two-layered medium, J. Comput. Phys., 229.2 (2010), 481-497.
  * [31] C. Lines, Inverse scattering by unbounded rough surfaces, PhD thesis, Department of Mathematics, Brunel University, U.K., 2003.
  * [32] J. Li, P. Li, H. Liu, and X. Liu, Recovering multiscale buried anomalies in a two-layered medium, Inverse Probl., 31 (2015), 105006.
  * [33] J. Li, G. Sun, and R. Zhang, The numerical solution of scattering by infinite rough interfaces based on the integral equation method, Comput. Math. Appl., 71 (2016), 1491-1502.
  * [34] P. Li, J. Wang, and L. Zhang, Inverse obstacle scattering for Maxwell’s equations in an unbounded structure, Inverse Probl., 35(2019), 095002.
  * [35] J. Li and J. Yang, Reverse time migration for inverse acoustic scattering by locally rough surfaces, preprint.
  * [36] J. Li, J. Yang, and B. Zhang, A linear sampling method for inverse acoustic scattering by a locally rough interface, Inverse Probl. Imag., 15(2021), 1247-1267.
  * [37] J. Li, J. Yang, and B. Zhang, Near-field imaging of a locally rough interface and buried obstacles with the linear sampling method, arXiv: 2009.00379.
  * [38] L. Li, J. Yang, B. Zhang, and H. Zhang, Imaging of buried obstacles in a two-layered medium with phaseless far-field data, Inverse Probl., 37(2021), 055004.
  * [39] X. Liu, B. Zhang, and H. Zhang, A direct imaging method for inverse scattering by unbounded rough surfaces, SIAM J. Imaging Sci., 11(2018), 1629-1650.
  * [40] X. Liu, B. Zhang, and H. Zhang, Near-field imaging of an unbounded elastic rough surface with a direct imaging method, SIAM J. Appl. Math., 79(2019), 153-176.
  * [41] D. Natroshvili, T. Arens, and S.N. Chandler-Wilde, Uniqueness, existence, and integral equation formulations for interface scattering problems, Mem. Differential Equations Math. Phys., 30(2003), 105-146.
  * [42] F. Qu, B. Zhang, and H. Zhang, A novel integral equation for scattering by locally rough surfaces and application to the inverse problem: The Neumann case, SIAM J. Sci. Comput., 41(2019), A3763-A3702.
  * [43] A. Willers and P. Werner, The helmholtz equation in disturbed half-spaces, Math. Method Appl. Sci., 9(1987), 312-323.
  * [44] M. Thomas, Analysis of Rough Surface Scattering Problems, PhD thesis, Department of Mathematics, The University of Reading, UK, 2006.
  * [45] J. Yang, J. Li, and B. Zhang, Simultaneous recovery of a locally rough interface and the embedded obstacle with its surrounding medium, Inverse Probl., 38(2022), 045011.
  * [46] B. Zhang and S.N. Chandler-Wilde, Integral equation methods for scattering by infinite rough surfaces, Math. Methods Appl. Sci., 26 (2003), 463-488.
  * [47] T. Zhu and J. Yang, A non-iterative sampling method for inverse elastic wave scattering by rough surfaces, Inverse Probl. Imag., accepted for publication.
  * [48] H. Zhang and B. Zhang, A novel integral equation for scattering by locally rough surfaces and application to the inverse problem, SIAM J. Appl. Math., 73 (2013), 1811-1829.