# A robust approach with numerical demonstrations for the inverse scattering
problem using a Carleman contraction map

Phuong M. Nguyen  Department of Mathematics and Statistics, University of
North Carolina at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>Loc H. Nguyen Department of Mathematics and Statistics, University of North
Carolina at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>
###### Abstract

This paper addresses the inverse scattering problem in a domain $\Omega$. The
input data, measured outside $\Omega$, involve the waves generated by the
interaction of plane waves with various directions and unknown scatterers
fully occluded inside $\Omega$. The output of this problem is the spatially
dielectric constant of these scatterers. Our approach to solving this problem
consists of two primary stages. Initially, we eliminate the unknown dielectric
constant from the governing equation, resulting in a system of partial
differential equations. Subsequently, we develop the Carleman contraction
mapping method to effectively tackle this system. It is noteworthy to
highlight this method’s robustness. It does not request a precise initial
guess of the true solution, and its computational cost is not expensive. Some
numerical examples are presented.

Keywords: inverse scattering problem; Carleman contraction mapping method;
globally convergence; low computational cost.

AMS subject classification: 35R25, 35N10, 35R30, 78A46

##  1 Introduction

Let $d\geq 2$ be the spatial dimension. Let $c:\mathbb{R}^{d}\to\mathbb{R}$ be
the spatially distributed dielectric constant of the medium. Let $\Omega$ be a
domain of $\mathbb{R}^{d}$ with a smooth boundary. Assume that $c$ satisfies

$c({\bf x})=\left\\{\begin{array}[]{ll}\geq 1&{\bf x}\in\Omega,\\\ 1&{\bf
x}\in\mathbb{R}^{d}\setminus\Omega.\end{array}\right.$ (1.1)

Setting $c({\bf x})=1$ outside $\Omega$ means that
$\mathbb{R}^{d}\setminus\Omega$ is either vacant or air-filled, whereas
$c({\bf x})\not=1$ inside $\Omega$ indicates the presence of entities such as
obstacles, acting as scatterers. To identify these scatterers, we use optical
waves (known as the incident wave) to illuminate $\Omega$. These incident
waves propagate through space and scatter in all directions. By recording the
scattered waves and using these data to calculate the spatially distributed
dielectric constant $c({\bf x}),$ ${\bf x}\in\Omega$, we can map the internal
geometry of $\Omega$. The knowledge of this function $c$ helps identify the
shape, position, and certain optical characteristics of the obstacles we aim
to detect. Typical examples of those scatterers are anti-personnel explosive
devices buried under the ground, cancerous tumors inside living tissues, and
nanostructures. Therefore, the success in finding $c$ strongly impacts many
other fields such as bio-medical imaging, nondestructive testing, radar,
detection and identification of explosive devices buried under the ground,
optical physics, seismic exploration, and nanoscience. See Figure 1 for a 2D
diagram of the scattering phenomenon.

A unknown scattererScattering waveScattering wave

Incident wave

$\Omega$

Figure 1:  A visualization of the inverse scattering process: An incident wave
hits a scatterer in an inaccessible domain $\Omega$. The interaction between
the wave and the scatterer produces a scattered wave, which contains important
information about the scatterer (shape, location, and some physical
properties). The primary aim of the inverse scattering problem is to use the
measurements of the scattering wave outside $\Omega$ to find this information.

In this paper, we use plane waves with a fixed wave number $k$ and various
angles $\theta\in[0,2\pi]$ as the incident waves. The incident plane wave is
given by

$u_{\rm inc}({\bf x},\theta)=e^{{\rm i}k\hat{\theta}\cdot{\bf
x}},\quad\mbox{for all }{\bf x}\in\mathbb{R}^{d},$ (1.2)

where

$\hat{\theta}=(\cos\theta,0,\dots,0,\sin\theta),\quad\theta\in[0,2\pi].$ (1.3)

The resulting total wave $u=u({\bf x},\theta)$, $({\bf
x},\theta)\in\mathbb{R}^{d}\times[0,2\pi]$, satisfies the Helmholtz equation
and the Sommerfeld radiation condition

$\left\\{\begin{array}[]{ll}\Delta u+k^{2}c({\bf x})u=0&{\bf
x}\in\mathbb{R}^{d},\\\ \partial_{|{\bf x}|}(u-u_{\rm inc})-{\rm i}k(u-u_{\rm
inc})=o(|{\bf x}|^{\frac{1-d}{2}})&|{\bf x}|\to\infty.\end{array}\right.$
(1.4)

We are interested in the following inverse scattering problem:

######  Problem 1.1 (Inverse scattering problem).

Given the boundary data

$f({\bf x},\theta)=u({\bf x},\theta)\quad\mbox{and}\quad g({\bf
x},\theta)=\partial_{\nu}u({\bf x},\theta)$ (1.5)

for all $({\bf x},\theta)\in\partial\Omega\times[0,2\pi],$ determine the
function $c({\bf x}),$ for all ${\bf x}\in\Omega.$

It is worth mentioning that the data in (1.5) are suitably requested. The
target function to be determined, $c({\bf x})$ has $d$ variables
$x_{1},\dots,x_{d}$. The domain $\partial\Omega\times[0,2\pi]$ of the data has
a combined dimension of $(d-1)+1$. From a mathematical standpoint, there is a
precise alignment between the dimensions of the required data and that of the
target function. Consequently, Problem 1.1 is efficiently designed, requiring
no redundant data, and is not over-determined. Our procedure to solve Problem
1.1 consists of two steps. In the first step, we eliminate the unknown
coefficient $c$ from (1.4). The result of this step is a system of partial
differential equations (PDEs), namely $\mathcal{L}{\bf v}=0$ where
$\mathcal{L}$ is a nonlinear differential operator. In the second step, we
solve the system $\mathcal{L}{\bf v}=0$ for the vector-valued function ${\bf
v}$. Once the solution ${\bf v}$ to the system above is found, we can compute
the solution to the governing equation. The solution to the inverse scattering
problem follows. The first step is based on a change of variable. The second
step is challenging since $\mathcal{L}$ is highly nonlinear. Rather than using
the optimization-based method, which requires a good initial guess of the true
solution, we will employ the recently developed method in [1, 37, 38, 49] to
find ${\bf v}$. More precisely, we construct a contraction mapping whose
unique fixed point is an approximation of ${\bf v}$. The solution $c$ to
Problem 1.1 follows.

There are many effective approaches to solve various versions of the inverse
scattering problem:

  * •

the imaging techniques based on sampling and the factorization methods [2, 3,
10, 13, 17, 22, 26, 42, 41, 45, 52];

  * •

the methods based on optimization [4, 11, 14, 20, 21, 53];

  * •

the methods based on Born series [9, 16, 19, 27, 35, 43];

  * •

the method based on linearization [5, 6, 7, 15];

  * •

the convexification method [23, 24, 25];

  * •

see [18] for a relatively complete list.

Although effective, each approach above has drawbacks. For example, the
imaging techniques only provide locations and shapes of the unknown objects;
the methods based on optimization are local in the sense that they require a
good initial guess of the true solution. There are many existing methods to
globalize the local algorithms, one of which is the convexification method,
first introduced in [29]. The key of the convexification method is to include
a Carleman weight function into the mismatch functional. Using Carleman
estimates, we can prove that the mismatch functional with the presence of
suitable Carleman weight functions are strictly convex. Hence, the mismatch
functional has a unique minimizer. We can rigorously prove that this unique
minimizer is an approximation of the true solution of the inverse problem as
the noise contained in the data tends to $0$ with the Lipchitz rate. Our
research team has developed new versions of the convexification method to
solve inverse scattering problem in the frequency domain [23, 24, 25] and in
the time domain [30]. These methods based on convexification were verified by
both simulated and experimental data. In this paper, as mentioned in the
previous paragraph, we will employ the contraction principle rather than
convexifying the mismatch functional whose unique minimizer yields the desired
solution. Given that the recursive iteration procedure to compute that fixed
point converges at a linear rate, by proposing this method, we will improve
the reconstruction quality and to reduce the computational time.

The paper is organized as follows. In Section 2, we derive an approximation
model, which is a system of PDEs. In Section 3, we define a contraction
mapping and propose an iterative procedure to solve that system. In Section 4,
we present some numerical examples and some details in implementation. Section
5 is for some discussions and concluding remarks.

##  2 An approximate model

In this section, we outline the derivation of a system of PDEs, the solution
of which provides the knowledge of $c$. While this step was introduced and
used previously in [23, 24, 25], we recall them here for completeness and for
the reader’s convenience. Introduce the change of variable

$v({\bf x},\theta)=\frac{1}{k^{2}}\log\frac{u({\bf x},\theta)}{u_{\rm
inc}({\bf x},\theta)}\quad\mbox{for all }({\bf
x},\theta)\in\Omega\times[0,2\pi].$ (2.1)

Defining the logarithm of the complex function $\frac{u({\bf
x},\theta)}{u_{\rm inc}({\bf x},\theta)}$, $({\bf
x},\theta)\in\Omega\times[0,2\pi]$, in (2.1) presents a challenge. However,
this logarithm is well-defined under some certain circumstances. For instance,
according to [32, Theorem 3], the function $\frac{u({\bf x},\theta)}{u_{\rm
inc}({\bf x},\theta)}$ can be approximated by $A({\bf x})e^{{\rm i}k(\tau({\bf
x})-{\bf x}\cdot\hat{\theta})}+O(1/k)$ when $k$ is sufficiently large, where
$A$ and $\tau$ are some positive functions. It is important to note that [32,
Theorem 3] was proved for a 3D plane wave in the form of $e^{{\rm i}kx_{3}}$.
However, with a variable change (if necessary), this theorem can be applied to
the incident wave in the form $e^{{\rm i}k{\bf x}\cdot\hat{\theta}}$ with any
direction as used in this paper. More significantly, according to [34, Theorem
3.3], the function $A$ exceeds a positive constant $a_{0}$. Therefore, the
function $v$ in (2.1) is well-defined and can be approximated by
$\frac{1}{k^{2}}[\ln A({\bf x})+{\rm i}k(\tau({\bf x})-{\bf
x}\cdot\hat{\theta})]$. In the numerical implementation, we defined the
logarithm of the above function straightforwardly, encountering no
difficulties.

By a straightforward algebra, we have

$\displaystyle\nabla v({\bf x},\theta)$ $\displaystyle=\frac{u_{\rm inc}({\bf
x},\theta)}{k^{2}u({\bf x},\theta)}\nabla\frac{u({\bf x},\theta)}{u_{\rm
inc}({\bf x},\theta)}=\frac{u_{\rm inc}({\bf x},\theta)}{k^{2}u({\bf
x},\theta)}\frac{u_{\rm inc}({\bf x},\theta)\nabla u({\bf x},\theta)-u({\bf
x},\theta)\nabla u_{\rm inc}({\bf x},\theta)}{u_{\rm inc}^{2}({\bf
x},\theta)}$ $\displaystyle=\frac{1}{k^{2}}\Big{[}\frac{\nabla u({\bf
x},\theta)}{u({\bf x},\theta)}-\frac{\nabla u_{\rm inc}({\bf
x},\theta)}{u_{\rm inc}({\bf
x},\theta)}\Big{]}=\frac{1}{k^{2}}\Big{[}\frac{\nabla u({\bf
x},\theta)}{u({\bf x},\theta)}-{\rm i}k\hat{\theta}\Big{]}$ (2.2)

and

$\Delta v({\bf x},\theta)=\frac{1}{k^{2}}\Big{[}\frac{\Delta u({\bf
x},\theta)}{u({\bf x},\theta)}-\Big{(}\frac{\nabla u({\bf x},\theta)}{u({\bf
x},\theta)}\Big{)}^{2}\Big{]}$ (2.3)

for all $({\bf x},\theta)\in\Omega\times[0,2\pi]$. It follows by the Helmholtz
equation in (1.4) that $\frac{\Delta u({\bf x},\theta)}{u({\bf
x},\theta)}=-k^{2}c({\bf x})$. On the other hand, by (2.2), $\frac{\nabla
u({\bf x},\theta)}{u({\bf x},\theta)}=k^{2}\nabla v({\bf x},\theta)+{\rm
i}k\hat{\theta}$. Combining these identities and (2.3) gives

$\displaystyle\Delta v({\bf x},\theta)$
$\displaystyle=\frac{1}{k^{2}}\Big{[}-k^{2}c({\bf x})-(k^{2}\nabla v({\bf
x},\theta)+{\rm i}k\hat{\theta})^{2}\Big{]}$ $\displaystyle=-c({\bf
x})-k^{2}(\nabla v({\bf x},\theta))^{2}-2{\rm i}k\nabla v({\bf
x},\theta)\cdot\hat{\theta}+1$

for all $({\bf x},\theta)\in\Omega\times[0,2\pi]$. Therefore, the function $v$
is governed by the equation

$\Delta v({\bf x},\theta)+2{\rm i}k\nabla v({\bf
x},\theta)\cdot\hat{\theta}+k^{2}(\nabla v({\bf x},\theta))^{2}=-c({\bf x})+1$
(2.4)

for all $({\bf x},\theta)\in\Omega\times[0,2\pi]$. Differentiate equation
(2.4) with respect to $\theta$ and note (from (1.3)) that

$\frac{d\hat{\theta}}{d\theta}=\hat{\theta}^{\prime}=(-\sin\theta,0,\dots,0,\cos\theta)\quad\mbox{for
all }\theta\in[0,2\pi].$ (2.5)

We have

$\Delta v_{\theta}({\bf x},\theta)+2{\rm i}k\nabla v_{\theta}({\bf
x},\theta)\cdot\hat{\theta}+2{\rm i}k\nabla v({\bf
x},\theta)\cdot\hat{\theta}^{\prime}+2k^{2}\nabla v({\bf x},\theta)\cdot\nabla
v_{\theta}({\bf x},\theta)=0$ (2.6)

for $({\bf x},\theta)\in\Omega\times[0,2\pi]$.

Let $\\{\Psi\\}_{n\geq 1}$ be the polynomial-exponential orthonormal basis of
$L^{2}(0,2\pi)$. This basis was constructed in [28] by using the Gram-Schmidt
process on the complete set
$\\{\phi_{n}(\theta)=\theta^{n-1}e^{\theta}\\}_{n\geq 1}$ of $L^{2}(0,2\pi)$.
See Remark 2.1 for the reason why we choose this basis among many other
orthonormal bases of $L^{2}(0,2\pi)$. For each ${\bf x}\in\Omega$, we write
the Fourier expansion of $v({\bf x},\theta)$ with respect to the basis
$\\{\Psi\\}_{n\geq 1}$

$v({\bf x},\theta)=\sum_{n=1}^{\infty}v_{n}({\bf
x})\Psi_{n}(\theta),\quad\mbox{where }v_{n}({\bf x})=\int_{0}^{2\pi}v({\bf
x},\theta)\Psi_{n}(\theta)d\theta.$ (2.7)

Motivated by the Galerkin approximation and since we pay attention on
computation, we approximate $v({\bf x},\theta)$ by truncating the series in
(2.7)

$v({\bf x},\theta)=\sum_{n=1}^{N}v_{n}({\bf x})\Psi_{n}(\theta),$ (2.8)

for $({\bf x},\theta)\in\Omega\times[0,2\pi]$, for some cut-off number $N$. We
will discuss a data-driven approach to choose $N$ later in Subsection 4.2.
Plugging (2.8) into (2.6), we have

$\sum_{n=1}^{N}\Delta v_{n}({\bf x})\Psi_{n}^{\prime}(\theta)+2{\rm
i}k\sum_{n=1}^{N}\nabla v_{n}({\bf
x})\Psi_{n}^{\prime}(\theta)\cdot\hat{\theta}+2{\rm i}k\sum_{n=1}^{N}\nabla
v_{n}({\bf x})\Psi_{n}(\theta)\cdot\hat{\theta}^{\prime}\\\
+2k^{2}\sum_{n=1}^{N}\sum_{l=1}^{N}\nabla v_{n}({\bf x})\cdot\nabla v_{l}({\bf
x})\Psi_{n}(\theta)\Psi_{l}^{\prime}(\theta)=0$ (2.9)

for $({\bf x},\theta)\in\Omega\times[0,2\pi]$. For each $m\in\\{1,\dots,N\\}$,
we multiply $\Psi_{m}(\theta)$ to both sides of (2.9). Then, integrating the
resulting equation, we obtain

$\sum_{n=1}^{N}s_{mn}\Delta v_{n}({\bf x})+\sum_{n=1}^{N}\nabla v_{n}({\bf
x})\cdot B_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}a_{mnl}\nabla v_{n}({\bf
x})\cdot\nabla v_{l}({\bf x})=0$ (2.10)

for all ${\bf x}\in\Omega$ where

$\displaystyle s_{mn}$
$\displaystyle=\int_{0}^{2\pi}\Psi_{n}^{\prime}(\theta)\Psi_{m}(\theta)d\theta,$
$\displaystyle B_{mn}$ $\displaystyle=2{\rm
i}k\Big{[}\int_{0}^{2\pi}\hat{\theta}\Psi_{n}^{\prime}(\theta)\Psi_{m}(\theta)d\theta+\int_{0}^{2\pi}\hat{\theta}^{\prime}\Psi_{n}(\theta)\Psi_{m}(\theta)d\theta\Big{]},$
$\displaystyle a_{mnl}$
$\displaystyle=2k^{2}\int_{0}^{2\pi}\Psi_{n}(\theta)\Psi_{l}^{\prime}(\theta)\Psi_{m}(\theta)d\theta.$

Due to the presence of $s_{mn}$, the system (2.10) might not be elliptic.
However, it is equivalent to an elliptic system. Recall that the matrix $S$ is
invertible and denote by $(\widetilde{s}_{mn})_{m,n=1}^{N}$ the matrix
$S^{-1}$. Define $\widetilde{B}_{mn}=\sum_{i=1}^{N}\widetilde{s}_{mi}B_{in}$
and $\widetilde{a}_{mnl}=\sum_{i=1}^{N}\widetilde{s}_{mi}a_{inl}$ for
$m,n,l\in\\{1,\dots,N\\}.$ It follows from (2.10) that

$\Delta v_{m}({\bf x})+\sum_{n=1}^{N}\nabla v_{n}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}({\bf x})\cdot\nabla v_{l}({\bf x})=0$ (2.11)

for ${\bf x}\in\Omega$. Coupling (2.11) for $m\in\\{1,\dots,N\\}$, we obtain a
nonlinear system for the vector ${\bf
v}=[\begin{array}[]{lll}v_{1}&\dots&v_{N}\end{array}]^{\rm T}$ where the
superscript T indicates the transpose.

######  Remark 2.1.

The system formed by (2.10) with $m\in\\{1,2,\dots,N\\}$ plays the key role in
our algorithm. It is our approximation model for Problem 1.1. The choice of
$\\{\Psi\\}_{n\geq 1}$ is important in establishing (2.10). For any
$N\in\mathbb{N}$, let $S$ be the $N\times N$ matrix $S$ whose $mn-$entry is
$s_{mn}=\int_{0}^{2\pi}\Psi_{n}^{\prime}(\theta)\Psi_{m}(\theta)d\theta$ as
above. It was proved in [28] that:

  1. 1.

the matrix $S$ is invertible;

  2. 2.

$\Psi^{\prime}_{n}$ is not identically zero in $[0,2\pi],$ for all $n\geq 1$.

It is obvious that the invertibility of $S$ is important. Without the
knowledge of the matrix $S^{-1}$, defining coefficients for (2.10) is
impossible. On the other hand, the second property mentioned above is also
holds significant. Specifically, if $\Psi^{\prime}_{n}=0$ for some $n$,
crucial information regarding $\Delta v_{n}$ in the first term of (2.9) would
be lost. To illustrate, consider replacing the polynomial-exponential basis
with widely used alternatives such as the Legendre polynomial or the
trigonometric basis, denoted by $\\{\phi_{n}\\}_{n\geq 1}$. In this scenario,
$\phi_{1}$ is a constant, resulting in $\phi_{1}^{\prime}\equiv 0$.
Consequently, the principal part $\Delta v_{1}$ would not be present in (2.9)
and, subsequently, in (2.10). As a consequence, the error in computing $v_{1}$
might be substantial, potentially leading to inaccuracies in the recovered
$v({\bf x},\theta)$ using (2.8). This is primarily because the first term
holds significant importance in this truncated series.

We next compute the boundary information for the vector $V$. For each
$m\in\\{1,\dots,N\\},$ due to (2.7),

$v_{m}({\bf x})=\int_{0}^{2\pi}v({\bf
x},\theta)\Psi_{m}(\theta)d\theta=\frac{1}{k^{2}}\int_{0}^{2\pi}\log\frac{u({\bf
x},\theta)}{u_{\rm inc}({\bf
x},\theta)}\Psi_{m}(\theta)d\theta=\frac{1}{k^{2}}\int_{0}^{2\pi}\log\frac{f({\bf
x},\theta)}{u_{\rm inc}({\bf x},\theta)}\Psi_{m}(\theta)d\theta$ (2.12)

for all ${\bf x}\in\partial\Omega.$ Here, we have used (2.1) and the knowledge
of the measurement $f({\bf x},\theta)=u({\bf x},\theta)$ on
$\partial\Omega\times[0,2\pi]$ as in (1.5). Similarly, using (2.2), we can
compute the normal derivative of $v_{m}$, $m\in\\{1,\dots,N\\}$, as

$\partial_{\nu}v_{m}({\bf
x})=\frac{1}{k^{2}}\int_{0}^{2\pi}\Psi_{m}(\theta)\Big{[}\frac{\nabla u({\bf
x},\theta)\cdot\nu}{u({\bf x},\theta)}-{\rm
i}k\hat{\theta}\cdot\nu\Big{]}d\theta=\frac{1}{k^{2}}\int_{0}^{2\pi}\Psi_{m}(\theta)\Big{[}\frac{g({\bf
x},\theta)}{f({\bf x},\theta)}-{\rm i}k\hat{\theta}\cdot\nu\Big{]}d\theta$
(2.13)

for all ${\bf x}\in\partial\Omega.$ Defining

$\displaystyle F_{m}({\bf x})$
$\displaystyle=\frac{1}{k^{2}}\int_{0}^{2\pi}\log\frac{f({\bf
x},\theta)}{u_{\rm inc}({\bf x},\theta)}\Psi_{m}(\theta)d\theta,$ (2.14)
$\displaystyle G_{m}({\bf x})$
$\displaystyle=\frac{1}{k^{2}}\int_{0}^{2\pi}\Psi_{m}(\theta)\Big{[}\frac{g({\bf
x},\theta)}{f({\bf x},\theta)}-{\rm i}k\hat{\theta}\cdot\nu\Big{]}d\theta$
(2.15)

for ${\bf x}\in\partial\Omega$, $m\in\\{1,\dots,N\\}$, we obtain a system for
${\bf v}=[\begin{array}[]{lll}v_{1}&\dots&v_{N}\end{array}]^{\rm T}$

$\left\\{\begin{array}[]{ll}\displaystyle\Delta v_{m}({\bf
x})+\sum_{n=1}^{N}\nabla v_{n}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}({\bf x})\cdot\nabla v_{l}({\bf x})=0&{\bf x}\in\Omega,\\\ v_{m}({\bf
x})=F_{m}({\bf x})&{\bf x}\in\partial\Omega,\\\ \partial_{\nu}v_{m}({\bf
x})=G_{m}({\bf x})&{\bf x}\in\partial\Omega\end{array}\right.$ (2.16)

for $m\in\\{1,\dots,N\\}.$

The derivation of (2.16) suggests a numerical method to solve Problem 1.1.
Initially, we compute the solution to the nonlinear problem (2.16). Next, we
use (2.8) to evaluate $v({\bf x},\theta)$ for $({\bf
x},\theta)\in\Omega\times[0,2\pi].$ Lastly, the desired solution $c({\bf x})$
can be computed using (2.4). Tackling (2.16) poses considerable difficulty due
to its nonlinear nature, prompting the application of the Picard iteration
method. Let ${\bf v}^{(0)}$ be an “initial guess” and assume by induction that
the vector ${\bf
v}^{(p)}=[\begin{array}[]{lll}v_{1}^{(p)}&\dots&v_{N}^{(p)}\end{array}]^{\rm
T}$, $p\geq 0$, is given. We compute ${\bf
v}^{(p+1)}=[\begin{array}[]{lll}v_{1}^{(p+1)}&\dots&v_{N}^{(p+1)}\end{array}]^{\rm
T}$ by solving

$\left\\{\begin{array}[]{ll}\displaystyle\Delta v_{m}^{(p+1)}({\bf
x})+\sum_{n=1}^{N}\nabla v_{n}^{(p+1)}({\bf x})\cdot\widetilde{B}_{mn}\\\
\hskip
142.26378pt\displaystyle+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p+1)}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x})=0&{\bf x}\in\Omega,\\\
v_{m}^{(p+1)}({\bf x})=F_{m}({\bf x})&{\bf x}\in\partial\Omega,\\\
\partial_{\nu}v_{m}^{(p+1)}({\bf x})=G_{m}({\bf x})&{\bf
x}\in\partial\Omega\end{array}\right.$ (2.17)

Given that (2.17) is linear with respect to ${\bf v}^{(p+1)}$, it can be
solved through the quasi-reversibility method or the least-squares
optimization. A critical consideration is the convergence of the sequence
$\\{{\bf v}^{(p)}\\}_{p\geq 0}$ to the actual solution ${\bf v}^{*}$ of
(2.16), particularly when the initial guess ${\bf v}^{(0)}$ far away from
${\bf v}^{*}$. To guarantee the convergence, we will include a Carleman weight
function in the least-squares optimization above. Details about combining a
Carleman weight function and the least squares optimization will be provided
in the next section.

######  Remark 2.2.

The Picard iterative scheme above is different from the ones suggested in [37,
38, 49]. In fact, if we applied the scheme in those papers, problem (2.17)
would be read as

$\left\\{\begin{array}[]{ll}\displaystyle\Delta v_{m}^{(p+1)}({\bf
x})+\sum_{n=1}^{N}\nabla v_{n}^{(p)}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p)}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x})=0&{\bf x}\in\Omega,\\\
v_{m}^{(p+1)}({\bf x})=F_{m}({\bf x})&{\bf x}\in\partial\Omega,\\\
\partial_{\nu}v_{m}^{(p+1)}({\bf x})=G_{m}({\bf x})&{\bf
x}\in\partial\Omega.\end{array}\right.$ (2.18)

In contrast, we solve ${\bf v}^{(p+1)}$ by solving (2.17) rather than (2.18)
because the nonlinear operator

$F_{\rm nonlinear}({\bf
v})=\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla v_{n}({\bf
x})\cdot\nabla v_{l}({\bf x}),\quad{\bf v}\in H,$

is not Lipschitz, as requested in [37, 38, 49]. Since the Lipschitz condition
is not satisfied, we numerically experience the divergence in computation when
applying the Picard scheme in those papers. Therefore, the change of replacing
(2.18) by (2.17) is significant to the success to our numerical scheme.

##  3 The Carleman contraction mapping method for (2.16)

The first version of the Carleman contraction mapping method was introduced in
[38]. We refer the reader to the following works [1, 37, 40, 44, 49] for
following up works and their application in several different inverse problems
for PDEs. The main idea of this method is to combine the well-known Picard
iteration, defined in the proof of the contraction mapping principle, and a
suitable Carleman estimate. For convenience, we recall here a Carleman
estimate, which was established in [49], which is important for the
convergence of our iterative method.

######  Lemma 3.1 (Carleman estimate).

Fix a point ${\bf x}_{0}\in\mathbb{R}^{d}\setminus\Omega$. Define $r({\bf
x})=|{\bf x}-{\bf x}_{0}|$ for all ${\bf x}\in\Omega.$ There exist positive
constants $\beta_{0}$ depending only on ${\bf x}_{0}$, $\Omega$, and $d$ such
that for all function $v\in C^{2}(\overline{\Omega})$ satisfying

$v({\bf x})=\partial_{\nu}v({\bf x})=0\quad\mbox{for all }{\bf
x}\in\partial\Omega,$ (3.1)

the following estimate holds true

$\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}|\Delta v|^{2}d{\bf x}\geq
C\lambda\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}|\nabla v({\bf
x})|^{2}\,d{\bf x}+C\lambda^{3}\int_{\Omega}e^{2\lambda r^{-\beta}({\bf
x})}|v({\bf x})|^{2}\,d{\bf x}$ (3.2)

for all $\beta>\beta_{0}$ and $\lambda\geq\lambda_{0}$. Here,
$\lambda_{0}=\lambda_{0}({\bf x}_{0},\Omega,d,\beta)$ and $C=C({\bf
x}_{0},\Omega,d,\beta)>0$ depend only on the listed parameters.

The derivation of Lemma 3.1 is directly based on [46, Lemma 5], see [39, Lemma
2.1] for a comprehensive explanation. Alternatively, another method to obtain
(3.2) involves using a distinct Carleman weight function through the Carleman
estimate presented in [36, Chapter 4, Section 1, Lemma 3], which applies to
general parabolic operators. It is important to point out to the reader the
diverse spectrum of Carleman estimates for different types of differential
operators (elliptic, parabolic, and hyperbolic) and their application in the
fields of inverse problems and computational mathematics, as discussed in
sources such as [8, 12, 31, 47]. Moreover, it is notable that certain Carleman
estimates are effective for any functions $v$ meeting the conditions
$v|_{\partial\Omega}=0$ and $\partial_{\nu}v|_{\Gamma}=0$, where $\Gamma$ is a
subset of $\partial\Omega$. This is exemplified in works like [33, 50]. These
Carleman estimates were especially useful in tackling quasilinear elliptic
PDEs with data limited to the portion $\Gamma$ of $\partial\Omega$.

Let $s$ be an integer with $s>\lceil\frac{d}{2}\rceil+2$ such that
$H^{s}(\Omega)$ can be continuously embedded into $C^{2}(\overline{\Omega})$.
For the theoretical purpose, we assume that the target function $c$ is
sufficiently smooth such that the true solution to the forward problem $u$ is
smooth with $\|u(\cdot,\theta)\|_{H^{s}(\overline{\Omega})}$ is finite
uniformly in $\theta\in[0,2\pi].$ We also assume that $|u({\bf x},\theta)|$ is
uniformly bounded from below by a positive constant $u_{0}$. Due to the change
of variable (2.1), the function $v({\bf x},\theta)$ is in the class
$H^{s}(\overline{\Omega})$ with $H^{s}$ norm being uniformly bounded for all
$\theta\in[0,2\pi].$ Hence, there is a number $M$, depending on an upper
bounded of $\|u(\cdot,\theta)\|_{H^{s}(\Omega)}$, $\theta\in[0,2\pi]$ and the
set $\\{\Psi_{n}\\}_{n=1}^{N}$, such that $\|{\bf
v}^{*}\|_{H^{s}(\overline{\Omega})^{N}}<M$ where ${\bf v}^{*}$ is the true
solution to (2.16). We find the solution to (2.16) in the set of admissible
solutions

$H=\big{\\{}{\bf v}=[\begin{array}[]{lll}v_{1}&\dots&v_{N}\end{array}]^{\rm
T}\in H^{s}(\Omega)^{N}:\|{\bf v}\|_{C^{2}(\overline{\Omega})^{N}}\leq
M,v_{m}|_{\partial\Omega}=F_{m}\mbox{ and
}\partial_{\nu}v_{m}|_{\partial\Omega}=G_{m}\big{\\}}.$ (3.3)

Throughout the paper, we assume that $H$ is nonempty and the true solution
$v^{*}$ to (2.16) belongs to $H$.

As outlined in the paragraph before Remark 2.2 in Section 2, we construct a
sequence $\\{{\bf v}^{(p)}\\}_{p\geq 0}$ with the hope that it converges to
the true solution to (2.16). Starting with ${\bf v}^{(0)}$ as a vector-valued
function in $H$, we proceed under the induction assumption that ${\bf
v}^{(p)}$, $p\geq 0$, is already determined. The subsequent vector ${\bf
v}^{(p+1)}$ is then derived by solving (2.17). Given that (2.17) is linear
with respect to ${\bf v}^{(p+1)}$, the least squares method (so-called the
quasi-reversibility method in solving PDEs) is applicable to compute ${\bf
v}^{(p+1)}$. To guarantee global convergence, it is important to include the
Carleman weight function $e^{2\lambda r^{-\beta}({\bf x})}$ in the
corresponding least-squares cost functional. More precisely, we set

${\bf v}^{(p+1)}=\underset{\varphi\in H}{\rm argmin}J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)$ (3.4)

where

$J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)=\sum_{m=1}^{N}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}\Big{|}\Delta\varphi_{m}({\bf
x})+\sum_{n=1}^{N}\nabla\varphi_{n}({\bf x})\cdot\widetilde{B}_{mn}\\\
+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla\varphi_{n}({\bf
x})\cdot\nabla v_{l}^{(p)}({\bf x})\Big{|}^{2}d{\bf
x}+\epsilon\|\varphi\|^{2}_{H^{s}(\Omega)^{N}}$ (3.5)

for all
$\varphi=[\begin{array}[]{lll}\varphi_{1}&\dots&\varphi_{N}\end{array}]^{\rm
T}\in H.$ Here, $\epsilon$ is a regularization parameter and
$\epsilon\|\varphi\|^{2}_{H^{s}(\Omega)^{N}}$ serves as a regularization term.
The “Carleman parameters” $\lambda$ and $\beta$ and the regularization
parameter $\epsilon$ will be chosen by a trial-error procedure, see Subsection
4.2 for more details. The following theorem holds.

######  Theorem 3.1.

Let $H$ be defined as in (3.3). Recursively define a sequence $\\{{\bf
v}^{(p)}\\}_{p\geq 0}$ as in (3.4) where ${\bf v}^{(0)}$ is a vector-valued
function in $H$. Assume the true solution ${\bf v}^{*}$ to (2.17) is in $H$
with $\|{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}<M$. Fix $\beta>\beta_{0}$ where
$\beta_{0}$ was defined in Lemma 3.1. Let $\lambda_{0}=\lambda_{0}({\bf
x}_{0},\Omega,d,\beta)$ be the number in Lemma 3.1. Then, for all
$\lambda\geq\lambda_{0}$, we have

$\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}[|\nabla({\bf v}^{(p+1)}-{\bf
v}^{*})|^{2}+|{\bf v}^{(p+1)}-{\bf v}^{*}|^{2}]\,d{\bf x}+\epsilon\|{\bf
v}^{(p+1)}-{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\\\
\leq\eta^{p+1}\Big{(}\int_{\Omega}e^{2\lambda r^{-\beta}({\bf
x})}\big{[}|\nabla\big{(}{\bf v}^{(0)}-{\bf v}^{*}\big{)}|^{2}+|{\bf
v}^{(0)}-{\bf v}^{*}|^{2}\big{]}d{\bf x}+\epsilon\|{\bf v}^{(0)}-{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}\\\ +\frac{\epsilon\eta}{1-\eta}\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}$ (3.6)

where $\eta\in(0,1)$ is a positive constant depending only on ${\bf x}_{0},$
$\Omega,$ $\beta$, $\\{\Psi_{n}\\}_{n=1}^{N}$, $M$, and $d$.

The conclusion of Theorem 3.1 is similar to those found in the convergence
theorems outlined in [37, 38, 49]. However, direct application of the
convergence theorems in those publications to derive Theorem 3.1 is not
possible. This limitation arises primarily because the nonlinearities present
in the PDEs of [37, 38, 49] must be Lipschitz, whereas the nonlinear function
${\bf v}\mapsto\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}({\bf x})\cdot\nabla v_{l}({\bf x})]$, $m=1,\dots,N$, in (2.16) is not.
To overcome this obstacle, we substitute the Lipschitz condition with the
boundedness assumption on $H$. The reader can find the convergence theorem
under the boundedness condition of $H$ in [40], in which we tackle the problem
of thermo/photo-acoustic tomography.

###### Proof of Theorem 3.1.

Fix $p\geq 0$. Recall $J^{{\bf v}^{(p)}}_{\lambda,\beta,\epsilon}$ as in
(3.5). Since ${\bf v}^{(p+1)}$ is the minimizer of $J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}$ in $H$ and $H$ is a convex subset of
$H^{s}(\Omega)^{N}$, for $0<t<1$, we have

${\bf v}^{(p+1)}-t({\bf v}^{(p+1)}-{\bf v}^{*})=(1-t){\bf v}^{(p+1)}+t{\bf
v}^{*}\in H.$ (3.7)

Define ${\bf h}=[\begin{array}[]{ccc}h_{1}&\dots&h_{N}\end{array}]^{\rm
T}={\bf v}^{(p+1)}-{\bf v}^{*}.$ We have

$\frac{J^{{\bf v}^{(p)}}_{\lambda,\beta,\epsilon}\big{(}{\bf v}^{(p+1)}-t{\bf
h}\big{)}-J^{{\bf v}^{(p)}}_{\lambda,\beta,\epsilon}\big{(}{\bf
v}^{(p+1)}\big{)}}{t}\geq 0.$ (3.8)

Using (3.5) and (3.8), we obtain

$\frac{1}{t}\sum_{m=1}^{N}\Big{[}\Big{\langle}e^{2\lambda r^{-\beta}({\bf
x})}\Big{(}2\Delta v_{m}^{(p+1)}({\bf x})+2\sum_{n=1}^{N}\nabla
v_{n}^{(p+1)}({\bf
x})\cdot\widetilde{B}_{mn}+2\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p+1)}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x})\Big{)}\\\ -t\Delta
h_{m}({\bf x})-t\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}-t\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
h_{n}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x}),\\\ -t\Delta h_{m}({\bf
x})-t\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}-t\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
h_{n}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf
x})\Big{\rangle}_{L^{2}(\Omega)^{N}}\\\ +\epsilon\langle 2{\bf
v}^{(p+1)}-t{\bf h},-t{\bf h}\rangle_{H^{s}(\Omega)^{N}}\Big{]}\geq 0.$ (3.9)

Let $t$ in (3.9) tend to $0$. We have

$\sum_{m=1}^{N}\Big{\langle}e^{2\lambda r^{-\beta}({\bf x})}\Big{(}\Delta
v_{m}^{(p+1)}({\bf x})+\sum_{n=1}^{N}\nabla v_{n}^{(p+1)}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p+1)}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x})\Big{)},\\\ -\Delta
h_{m}({\bf x})-\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}-\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
h_{n}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf
x})\Big{\rangle}_{L^{2}(\Omega)^{N}}\\\ -\epsilon\langle{\bf v}^{(p+1)},{\bf
h}\rangle_{H^{s}(\Omega)^{N}}\geq 0.$ (3.10)

or, equivalently,

$\sum_{m=1}^{N}\Big{\langle}e^{2\lambda r^{-\beta}({\bf x})}\Big{(}\Delta
v_{m}^{(p+1)}({\bf x})+\sum_{n=1}^{N}\nabla v_{n}^{(p+1)}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p+1)}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x})\Big{)},\\\ +\Delta
h_{m}({\bf x})+\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
h_{n}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf
x})\Big{\rangle}_{L^{2}(\Omega)^{N}}\\\ +\epsilon\langle{\bf v}^{(p+1)},{\bf
h}\rangle_{H^{s}(\Omega)^{N}}\leq 0.$ (3.11)

On the other hand, since ${\bf v}^{*}$ solves (2.16),

$\sum_{m=1}^{N}\Big{\langle}e^{2\lambda r^{-\beta}({\bf x})}\Big{(}\Delta
v_{m}^{*}({\bf x})+\sum_{n=1}^{N}\nabla v_{n}^{*}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{*}({\bf x})\cdot\nabla v_{l}^{*}({\bf x})\Big{)},\\\ \Delta h_{m}({\bf
x})+\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
h_{n}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf
x})\Big{\rangle}_{L^{2}(\Omega)^{N}}\\\ +\epsilon\langle{\bf v}^{*},{\bf
h}\rangle_{H^{s}(\Omega)^{N}}=\epsilon\langle{\bf v}^{*},{\bf
h}\rangle_{H^{s}(\Omega)^{N}}.$ (3.12)

Combining (3.11) and (3.17) gives

$\sum_{m=1}^{N}\Big{\langle}e^{2\lambda r^{-\beta}({\bf
x})}\Big{(}\Delta(v_{m}^{(p+1)}({\bf x})-v^{*}_{m}({\bf
x}))+\sum_{n=1}^{N}\nabla(v_{n}^{(p+1)}({\bf x})-v_{n}^{*}({\bf
x}))\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p+1)}({\bf x})\cdot\nabla v_{l}^{(p)}({\bf x})\\\
-\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla v_{n}^{*}({\bf
x})\cdot\nabla v_{l}^{*}({\bf x})\Big{)},\Delta h_{m}({\bf
x})+\sum_{n=1}^{N}\nabla h_{n}({\bf x})\cdot\widetilde{B}_{mn}\\\
+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla h_{n}({\bf
x})\cdot\nabla v_{l}^{(p)}({\bf
x})\Big{\rangle}_{L^{2}(\Omega)^{N}}+\epsilon\langle{\bf v}^{(p+1)}-{\bf
v}^{*},{\bf h}\rangle_{H^{s}(\Omega)^{N}}\leq-\epsilon\langle{\bf v}^{*},{\bf
h}\rangle_{H^{s}(\Omega)^{N}}.$ (3.13)

Recall that ${\bf h}={\bf v}^{(p+1)}-{\bf v}^{*}$. We obtain from (3.19) that

$\sum_{m=1}^{N}\Big{\langle}e^{2\lambda r^{-\beta}({\bf x})}\Big{(}\Delta
h_{m}({\bf x})+\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla
v_{n}^{(p+1)}({\bf x})\cdot\nabla[v_{l}^{(p)}({\bf x})-v^{*}_{l}({\bf x})]\\\
+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla h_{n}({\bf
x})\cdot\nabla v_{l}^{*}({\bf x})\Big{)},\Delta h_{m}({\bf
x})+\sum_{n=1}^{N}\nabla h_{n}({\bf x})\cdot\widetilde{B}_{mn}\\\
+\sum_{n=1}^{N}\sum_{l=1}^{N}\widetilde{a}_{mnl}\nabla h_{n}({\bf
x})\cdot\nabla v_{l}^{(p)}({\bf
x})\Big{\rangle}_{L^{2}(\Omega)^{N}}+\epsilon\|{\bf
h}\|_{H^{s}(\Omega)^{N}}^{2}\leq-\epsilon\langle{\bf v}^{*},{\bf
h}\rangle_{H^{s}(\Omega)^{N}}.$

Thus, by using the inequality $2ab\leq\frac{1}{4}a^{2}+4b^{2}$, the Hölder
inequality, and noting that $\|{\bf v}^{(p+1)}\|_{H^{s}(\Omega)^{N}}\leq M$
and $\|{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}\leq M$ we get

$\sum_{m=1}^{N}\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}\Big{|}\Delta
h_{m}({\bf x})+\sum_{n=1}^{N}\nabla h_{n}({\bf
x})\cdot\widetilde{B}_{mn}\Big{|}^{2}d{\bf x}+\epsilon\|{\bf
h}\|_{H^{s}(\Omega)^{N}}^{2}\\\ \leq C\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}|\nabla\big{(}{\bf v}^{(p)}({\bf x})-{\bf v}^{*}({\bf
x})\big{)}|^{2}d{\bf x}+\int_{\Omega}e^{2\lambda r^{-\beta}({\bf
x})}|\nabla{\bf h}|^{2}d{\bf x}\Big{)}+\epsilon C\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}.$ (3.14)

Here, $C$ is a positive constant depending only on $\\{\Psi_{n}\\}_{n=1}^{N}$
and $M$. Using the inequality $(a+b)^{2}\geq\frac{a^{2}}{2}-b^{2}$, we obtain
from (3.14) that

$\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}\Big{|}\Delta{\bf h}({\bf
x})\Big{|}^{2}d{\bf x}+\epsilon\|{\bf h}\|_{H^{s}(\Omega)^{N}}^{2}\\\ \leq
C\Big{(}\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}|\nabla\big{(}{\bf
v}^{(p)}({\bf x})-{\bf v}^{*}({\bf x})\big{)}|^{2}d{\bf
x}+\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}|\nabla{\bf h}|^{2}d{\bf
x}\Big{)}+\epsilon C\|{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}.$ (3.15)

From now on, we allow the universal constant $C$ to depend not only on
$\\{\Psi_{n}\\}_{n=1}^{N}$ and $M$ but also on ${\bf x}_{0}$, $\Omega,$
$\beta$, and $d$. Recall the number $\lambda_{0}$ as in Lemma 3.1. Applying
the Carleman estimate (3.2) for each entry of ${\bf h}$, we obtain

$\displaystyle\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}|\Delta{\bf
h}|^{2}d{\bf x}$ $\displaystyle\geq C\lambda\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}|\nabla{\bf h}({\bf x})|^{2}\,d{\bf
x}+C\lambda^{3}\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}|{\bf h}({\bf
x})|^{2}\,d{\bf x}$ $\displaystyle\geq C\lambda\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}[|\nabla{\bf h}({\bf x})|^{2}+|{\bf h}({\bf
x})|^{2}]\,d{\bf x}.$ (3.16)

It follows from (3.15) and (3.16) that

$\lambda\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}[|\nabla{\bf h}({\bf
x})|^{2}+|{\bf h}({\bf x})|^{2}]\,d{\bf x}+\epsilon\|{\bf
h}\|_{H^{s}(\Omega)^{N}}^{2}\\\ \leq C\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}|\nabla\big{(}{\bf v}^{(p)}({\bf x})-{\bf v}^{*}({\bf
x})\big{)}|^{2}d{\bf x}+\int_{\Omega}e^{2\lambda r^{-\beta}({\bf
x})}|\nabla{\bf h}|^{2}d{\bf x}\Big{)}+\epsilon C\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}.$ (3.17)

Since (3.17) holds true for all $\lambda>\lambda_{0}$, we can choose $\lambda$
large such that the term $C\int_{\Omega}e^{2\lambda r^{-\beta}({\bf
x})}|\nabla{\bf h}({\bf x})|^{2}d{\bf x}$ in the right-hand side of (3.17) is
absorbed into the left-hand side. Recall that ${\bf h}={\bf v}^{(p+1)}-{\bf
v}^{*}$. We obtain

$\displaystyle\int_{\Omega}$ $\displaystyle e^{2\lambda r^{-\beta}({\bf
x})}[|\nabla({\bf v}^{(p+1)}-{\bf v}^{*})|^{2}+|{\bf v}^{(p+1)}-{\bf
v}^{*}|^{2}]\,d{\bf x}+\epsilon\|{\bf v}^{(p+1)}-{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}$
$\displaystyle\leq\frac{C}{\lambda}\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}|\nabla\big{(}{\bf v}^{(p)}-{\bf v}^{*}\big{)}|^{2}d{\bf
x}+\epsilon\|{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}$
$\displaystyle\leq\frac{C}{\lambda}\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}\big{[}|\nabla\big{(}{\bf v}^{(p)}-{\bf
v}^{*}\big{)}|^{2}+|{\bf v}^{(p)}-{\bf v}^{*}|^{2}\big{]}d{\bf
x}+\epsilon\|{\bf v}^{(p)}-{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}.$ (3.18)

Replacing $p$ by $p-1$ in (3.18), we obtain

$\displaystyle\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}$
$\displaystyle[|\nabla({\bf v}^{(p+1)}-{\bf v}^{*})|^{2}+|{\bf v}^{(p+1)}-{\bf
v}^{*}|^{2}]\,d{\bf x}+\epsilon\|{\bf v}^{(p+1)}-{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}$
$\displaystyle\leq\frac{C}{\lambda}\Big{[}\frac{C}{\lambda}\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}\big{[}|\nabla\big{(}{\bf v}^{(p-1)}-{\bf
v}^{*}\big{)}|^{2}+|{\bf v}^{(p-1)}-{\bf v}^{*}|^{2}\big{]}d{\bf
x}+\epsilon\|{\bf v}^{(p-1)}-{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}$
$\displaystyle\hskip 227.62204pt+\epsilon\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}\Big{)}+\epsilon\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}\Big{]}$
$\displaystyle=\Big{(}\frac{C}{\lambda}\Big{)}^{2}\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}\big{[}|\nabla\big{(}{\bf v}^{(p-1)}-{\bf
v}^{*}\big{)}|^{2}+|{\bf v}^{(p-1)}-{\bf v}^{*}|^{2}\big{]}d{\bf
x}+\epsilon\|{\bf v}^{(p-1)}-{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}$
$\displaystyle\hskip
227.62204pt+\Big{(}\Big{(}\frac{C}{\lambda}\Big{)}^{2}+\frac{C}{\lambda}\Big{)}\epsilon\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}.$

Continuing this process, we obtain

$\int_{\Omega}e^{2\lambda r^{-\beta}({\bf x})}[|\nabla({\bf v}^{(p+1)}-{\bf
v}^{*})|^{2}+|{\bf v}^{(p+1)}-{\bf v}^{*}|^{2}]\,d{\bf x}+\epsilon\|{\bf
v}^{(p+1)}-{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\\\
\leq\Big{(}\frac{C}{\lambda}\Big{)}^{p+1}\Big{(}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}\big{[}|\nabla\big{(}{\bf v}^{(0)}-{\bf
v}^{*}\big{)}|^{2}+|{\bf v}^{(0)}-{\bf v}^{*}|^{2}\big{]}d{\bf
x}+\epsilon\|{\bf v}^{(0)}-{\bf v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}\Big{)}\\\
+\sum_{i=1}^{p+1}\Big{(}\frac{C}{\lambda}\Big{)}^{i}\epsilon\|{\bf
v}^{*}\|_{H^{s}(\Omega)^{N}}^{2}.$ (3.19)

Define $\eta=\frac{C}{\lambda}\in(0,1)$. The assertion (3.6) is directly
deduced from (3.19). ∎

In theory, the estimate in (3.6) guarantee that when $p$ is sufficiently
large, ${\bf v}^{(p)}$ well-approximates the solution ${\bf v}^{*}$ to (2.16).
However, in our numerical study, the approximation already meets our
expectations when $p=5$. As mentioned in Section 2, once (2.16) is solved, the
computed solution to the inverse scattering problem follows. Our algorithm to
solve this inverse scattering problem is presented in Algorithm 1. Some
details about implementation of Algorithm 1 will be presented in Subsection
4.2.

Algorithm 1 The Carleman contraction mapping method to compute the numerical
solution to Problem 1.1

1: Choose a cut-off number $N$. Choose Carleman parameters ${\bf x}_{0},$
$\beta$, and $\lambda$ and a regularization parameter $\epsilon$.

2: Choose an initial solution ${\bf v}^{(0)}({\bf x})\in H$ and a number
$P\in\mathbb{N}$.

3:for $p=0$ to $P-1$ do

4: Set ${\bf v}^{(p+1)}=\underset{\varphi\in H}{\rm argmin}J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)$ where $J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}$ was defined in (3.4).

5:end for

6:Set ${\bf v}^{\rm comp}={\bf v}^{(P)}=\left[\begin{array}[]{cccc}v_{1}^{\rm
comp}&v_{2}^{\rm comp}&\dots&v_{N}^{\rm comp}\end{array}\right]^{\rm T}$. By
(2.8), we compute $v^{\rm comp}({\bf x},\theta)$ using  $v^{\rm comp}({\bf
x},\theta)=\sum_{n=1}^{N}v_{n}^{\rm comp}({\bf x})\Psi_{n}(\theta),$ (3.20)
for ${\bf x}\in\Omega,$ $\theta\in[0,2\pi]$.

7:Due to (2.4), a numerical solution to Problem 1.1 can be computed via
$c^{\rm comp}({\bf x})=\frac{1}{2\pi}\int_{0}^{2\pi}\Big{|}\Delta v({\bf
x},\theta)+2{\rm i}k\nabla v({\bf x},\theta)\cdot\hat{\theta}+k^{2}(\nabla
v({\bf x},\theta))^{2}\Big{|}d{\bf x}d\theta+1$ (3.21) for all ${\bf
x}\in\Omega.$

######  Remark 3.1.

For any given ${\bf v}^{(p)}\in H$, it follows from (3.18) that ${\bf
v}^{(p+1)}=\underset{\varphi\in H}{\rm argmin}J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)$ is closer to ${\bf v}^{*}$, in
comparison to ${\bf v}^{(0)}$, with respect to the norm

$\|\varphi\|_{\lambda,\beta,\epsilon}=\Big{[}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}[|\nabla\varphi|^{2}+|\varphi|^{2}]\,d{\bf
x}+\epsilon\|\varphi\|_{H^{s}(\Omega)^{N}}^{2}\Big{]}^{\frac{1}{2}}.$

So, the map $H\ni{\bf v}\mapsto\underset{\varphi\in H}{\rm argmin}J^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)$ is a contractor toward ${\bf
v}^{*}$ with respect to the norm above. This contraction behavior of this map,
together with the involvement of the Carleman weight function $e^{2\lambda
r^{-\beta}({\bf x})}$, suggests the name of our method: the Carleman
contraction mapping method.

##  4 Numerical study

We display some numerical results based on Algorithm 1 in this section. For
simplicity, we only implement this algorithm in 2D.

###  4.1 Data generation

In this section, we discuss how to solve the forward problem to generate the
data $f$ and $g$ for Problem 1.1. Rather than solving (1.4) in the whole
$\mathbb{R}^{2}$, which is complicated, we approximate this model by a Robin-
type boundary problem for the Helmholtz equation

$\left\\{\begin{array}[]{ll}\Delta u({\bf x},\theta)+k^{2}c({\bf x})u({\bf
x},\theta)=0&{\bf x}\in\Omega,\\\ \partial_{\nu}u_{\rm sc}({\bf
x},\theta)-{\rm i}ku_{\rm sc}({\bf x},\theta)=0&{\bf
x}\in\partial\Omega,\end{array}\right.$ (4.1)

for $\theta\in[0,2\pi].$ Here, $u_{\rm sc}({\bf x},\theta)=u({\bf
x},\theta)-u_{\rm inc}({\bf x},\theta)$ is the scattering wave. This
adjustment does not impact our methodology, as Algorithm 1 is constructed
based on the form of the Helmholtz equation within $\Omega$ while the
equation’s form outside $\Omega$ is irrelevant. In this paper, we set
$\Omega=(-1,1)^{2}$. On $\overline{\Omega}$, we arrange a uniform $64\times
64$ grid of points

$\mathcal{G}=\Big{\\{}(x_{i}=-1+(i-1)d_{\bf x},y_{j}=-1+(j-1)d_{\bf x}):1\leq
i,j\leq 64\Big{\\}}$

where $d_{\bf x}=\frac{2}{63}$ is the step size in space. We also discretize
the interval of angle $[0,2\pi]$ by the partition

$\Theta=\big{\\{}\theta_{1},\theta_{2},\dots,\theta_{N_{\theta}}\big{\\}}$

where $N_{\theta}=150$, $d_{\theta}=2\pi/149$, and
$\theta_{i}=(i-1)d_{\theta}.$ It is not hard to verify that for all
$\theta\in[0,2\pi],$ the scattering wave $u_{\rm sc}$ satisfies the following
problem

$\left\\{\begin{array}[]{ll}\Delta u_{\rm sc}({\bf x},\theta)+k^{2}c({\bf
x})u_{\rm sc}({\bf x},\theta)=-k^{2}(c({\bf x})-1)u_{\rm inc}({\bf
x},\theta)&{\bf x}\in\Omega,\\\ \partial_{\nu}u_{\rm sc}({\bf x},\theta)-{\rm
i}ku_{\rm sc}({\bf x},\theta)=0&{\bf x}\in\partial\Omega.\end{array}\right.$
(4.2)

In computation, we solve (4.2) by the finite difference method for $u_{\rm
sc}$. Then, we compute $u$ by $u_{\rm sc}+u_{\rm inc}.$ Then, we extract the
“noiseless” boundary data $f^{*}({\bf x},\theta)$ and $g^{*}({\bf x},\theta)$
for all $({\bf x},\theta)\in\mathcal{G}\times\Theta.$ The noisy data are given
by

$f=f^{*}(1+\delta{\rm rand}),\quad g=g^{*}(1+\delta{\rm rand})$

where rand is a function that generates random complex numbers with the range
of both real part and imaginary part being $[-1,1].$ The number $\delta$ is
the noise level. In this paper, we set $\delta=10\%.$

###  4.2 The implementations

In this subsection, we discuss how we generate the simulated data and how we
implement steps 1, 2, and 4 of Algorithm 1. The implementation of other steps
is straightforward.

Step 1: Choosing artificial parameters. We use the information from the data
$f$ to determine the appropriate cut-off parameter, $N$. Given that $\Omega$
is a square, it is convenient in numerical implementation to examine the data
on one side, denoted as $\Gamma$, of its boundary $\partial\Omega$. For each
$N$, we define the discrepancy between $f|_{\Gamma\times[0,2\pi]}$ and its
approximation obtained by truncating the Fourier series
$\displaystyle\sum_{n=1}^{N}f_{n}(\mathbf{x})\Psi_{n}(\theta)$, where
$f_{n}({\bf x})=\displaystyle\int_{0}^{2\pi}f({\bf
x},\theta)\Psi_{n}(\theta)d\theta$. This difference is:

$e(N)=\frac{\Big{\|}f-\displaystyle\sum_{n=1}^{N}f_{n}({\bf
x})\Psi_{n}(\theta)\Big{\|}_{L^{2}(\Gamma\times[0,2\pi])}}{\|f\|_{L^{2}(\Gamma\times[0,2\pi])}}.$

Figure 2 illustrates the plots of the function $e(N)$, where the data is
sampled from the set

$\Gamma=\\{(x,y=-1):|x|\leq 1\\}\subset\partial\Omega.$

These data are generated for all numerical examples presented in Subsection
4.3. In theory, the function $e(N)$ should decrease with respect to $N$.
However, due to the high oscillation of $\Psi_{n}$ for large $n$, leading to
errors in computing the integral in $\displaystyle f_{n}({\bf
x})=\int_{0}^{2\pi}f({\bf x},\theta)\Psi_{n}(\theta)d\theta$ using the finite
difference method, this function experiences an increase at a certain critical
value of $N$. That means the approximation is no longer true when $N$ is
larger than this critical value. We select this critical value of $N$ as the
cut-off parameter in Step 1. Based on Figure 2, we choose $N=42$.

(a) Test 1

(b) Test 2

(c) Test 3

(d) Test 4

Figure 2:  The graphs of the functions $e(N)$ are computed by the data, each
containing 10% noise, generated for all four numerical tests in Subsection
4.3. These graphs suggest us to choose $N=42$.

A trial and error procedure chooses the other artificial parameters ${\bf
x}_{0},\beta,\lambda$, and $\epsilon$. We use Test 1 in Subsection 4.3 as a
reference test, for which the true solution is known. These parameters are
manually adjusted until the computed dielectric constant achieves an
acceptable level of accuracy. Once determined, these parameters are then
consistently applied across all subsequent tests. In our numerical study, we
set $\mathbf{x}_{0}=(0,-10)$, $\beta=20$, $\lambda=6$, and
$\epsilon=10^{-5.5}$.

Step 2: Computing the initial solution. Choosing the vector-valued function
${\bf v}^{(0)}\in H$ is natural. To so so, we simply solve the over-determined
boundary value problem

$\left\\{\begin{array}[]{ll}\Delta v_{m}^{(0)}=0&\mbox{in }\Omega,\\\
v_{m}^{(0)}({\bf x})=F_{m}({\bf x})&{\bf x}\in\partial\Omega,\\\
\partial_{\nu}v_{m}^{(0)}({\bf x})=G_{m}({\bf x})&{\bf
x}\in\partial\Omega\end{array}\right.$ (4.3)

by the regularized least square method, or so-called the quasi-reversibility
method. That means we minimize the cost functional and set the minimizer as
${\bf v}^{(0)}$

${\bf v}^{(0)}=\underset{\varphi\in H}{{\rm
argmin}}\int_{\Omega}|\Delta\varphi|^{2}d{\bf x}.$

We find that this initial estimate for $\mathbf{v}^{*}$ may not be
particularly close to the true solution. However, it is a robust starting
point for the iterative process in Algorithm 1.

Alternatively, a simpler choice for $\mathbf{v}^{(0)}$ is to set
$\mathbf{v}^{(0)}$ as the zero vector everywhere. Despite the fact that
$\mathbf{v}^{(0)}\equiv 0$ does not belong to $H$, we have observed that
Algorithm 1 still yields surprisingly satisfactory numerical results.

Step 4: Computing ${\bf v}^{(p+1)}$ given the knowledge of ${\bf v}^{(p)}$. We
observe that although $S$ is invertible, the matrix $S^{-1}$ contains
“singular” entries with large number when $N=42$ as in our choice. This might
cause some unnecessary errors in computation. We, therefore, rather than solve
(2.16), we solve its equivalent system in which the partial differential
equations are replaced by (2.10). The corresponding definition of ${\bf
v}^{(p+1)}$ is

${\bf v}^{(p+1)}=\underset{\varphi\in H}{\rm argmin}\widetilde{J}^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)$ (4.4)

where

$\widetilde{J}^{{\bf
v}^{(p)}}_{\lambda,\beta,\epsilon}(\varphi)=\sum_{m=1}^{N}\int_{\Omega}e^{2\lambda
r^{-\beta}({\bf x})}\Big{|}\sum_{n=1}^{N}s_{mn}\Delta\varphi_{n}({\bf x})\\\
+\sum_{n=1}^{N}\nabla\varphi_{n}({\bf x})\cdot
B_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}a_{mnl}\nabla\varphi_{n}({\bf
x})\cdot\nabla v_{l}^{(p)}({\bf x})\Big{|}^{2}d{\bf
x}+\epsilon\|\varphi\|^{2}_{H^{s}(\Omega)^{N}}$ (4.5)

for all
$\varphi=[\begin{array}[]{lll}\varphi_{1}&\dots&\varphi_{N}\end{array}]^{\rm
T}\in H.$

Note that the maps

$H\ni\varphi\mapsto e^{\lambda r^{-\beta}({\bf
x})}\sum_{n=1}^{N}s_{mn}\Delta\varphi_{n}({\bf x})\\\
+\sum_{n=1}^{N}\nabla\varphi_{n}({\bf x})\cdot
B_{mn}+\sum_{n=1}^{N}\sum_{l=1}^{N}a_{mnl}\nabla\varphi_{n}({\bf
x})\cdot\nabla v_{l}^{(p)}({\bf x})$

and

$H\ni\varphi\mapsto\sum_{|\alpha|\leq s}D^{\alpha}\varphi$

are linear. Consequently, we can discretize all functions involved in (4.5)
using a finite difference scheme. For a comprehensive understanding of the
discretization process, we refer interested readers to [48, 51]. The
minimization procedure after discretion can be executed using the optimization
packages available in Matlab. In our implementation, we employ the lsqlin
command, which is a built-in tool in Matlab, for this purpose.

###  4.3 Numerical examples

We display four numerical tests.

Test 1. We test the case when “two scatterers” with disk shape and with
different values are presented inside $\Omega$. The function $c^{\rm true}$ is
given by

$c^{\rm true}(x,y)=\left\\{\begin{array}[]{ll}2&(x+0.5)^{2}+(y+0.5)^{2}\leq
0.04,\\\ 1.5&(x-0.5)^{2}+(y-0.5)^{2}\leq 0.04,\\\
1&\mbox{otherwise}.\end{array}\right.$ (4.6)

In this test, the wave number $k$ is $3\pi.$ The true and constructed
spatially distributed dielectric constants are displayed in Figure 3.

(a) The function $c^{\rm true}$

(b) The function $c^{\rm comp}$

(c)  The consecutive difference $\frac{\|{\bf v}^{(p)}-{\bf
v}^{(p-1)}\|_{L^{2}(\Omega)}}{\|{\bf v}^{(p)}\|_{L^{2}(\Omega)}}$

Figure 3:  The actual and computed spatial dielectric constants $c$ along with
the $L^{2}$ relative difference of the computed $\mathbf{v}$ are presented. It
is seen that the computed function $c$ meets the acceptable criteria.
Additionally, Figure (c) provides numerical evidence supporting the rapid
convergence of Algorithm 1. Notably, it can be observed that only 3 iterations
are adequate to achieve an accurate reconstruction. The data in this test is
corrupted with 10% of noise.

In this test, the presence of two scatterers within $\Omega$ is clearly seen.
The errors in the computed results remain within acceptable bounds. The true
value of $c$ within the lower scatterer is known to be 2, while the maximum
value of $c$ calculated within the computed scatterer reaches 1.78, resulting
in a relative error of 11.0%. For the upper scatterer, the true $c$ value is
1.5, with the maximum computed $c$ value inside the computed scatterer being
1.42, yielding a relative error of 5.3%. Both computed errors are compatible
with the noise level $10\%$ present in the data. We highlight that Figure 3c
serves as numerical verification of the convergence, proved in Theorem 3.1.

Test 2. We test the case when “one rectangle scatterer” is occluded inside
$\Omega$. The function $c^{\rm true}$ is given by

$c^{\rm
true}(x,y)=\left\\{\begin{array}[]{ll}2&\max\Big{\\{}\frac{|x-0.5|}{0.5},\frac{|y|}{1.5}\Big{\\}}\leq
0.3,\\\ 1&\mbox{otherwise}.\end{array}\right.$ (4.7)

In this test, the wave number $k$ is $2\pi.$ The true and constructed
spatially distributed dielectric constants are displayed in Figure 4.

(a) The function $c^{\rm true}$

(b) The function $c^{\rm comp}$

(c)  The consecutive difference $\frac{\|{\bf v}^{(p)}-{\bf
v}^{(p-1)}\|_{L^{2}(\Omega)}}{\|{\bf v}^{(p)}\|_{L^{2}(\Omega)}}$

Figure 4:  The actual and computed spatial dielectric constants $c$ along with
the $L^{2}$ relative difference of the computed $\mathbf{v}$ are presented. It
is seen that the computed function $c$ meets the acceptable criteria.
Additionally, Figure (c) provides numerical evidence supporting the rapid
convergence of Algorithm 1. Notably, it can be observed that only 3 iterations
are adequate to achieve an accurate reconstruction. The data in this test is
corrupted with 10% of noise.

In this test, the presence of the rectangular scatterer $\Omega$ is clearly
seeable, with an acceptable reconstructed shape and location. Regarding the
reconstructed value, the true value of $c$ in the scatterer is 2 while the
maximum value of $c$ calculated within the reconstructed scatterer reaches
1.79, resulting in a relative error of 10.5%. Similarly to Test 1, Figure 4c
not only visually confirms the convergence, as proven in Theorem 3.1, but also
serves as numerical validation of the convergence behavior detailed in the
theorem.

Test 3. We test the case when “one square scatterer” is placed at the center
of $\Omega$, which is away from the measurement sites. The function $c^{\rm
true}$ is given by

$c^{\rm true}(x,y)=\left\\{\begin{array}[]{ll}2.5&\max\\{|x|,|y|\\}\leq
0.09,\\\ 1&\mbox{otherwise}.\end{array}\right.$ (4.8)

In this test, the wave number $k$ is $2\pi.$ The true and constructed
spatially distributed dielectric constants are displayed in Figure 5.

(a) The function $c^{\rm true}$

(b) The function $c^{\rm comp}$

(c)  The consecutive difference $\frac{\|{\bf v}^{(p)}-{\bf
v}^{(p-1)}\|_{L^{2}(\Omega)}}{\|{\bf v}^{(p)}\|_{L^{2}(\Omega)}}$

Figure 5:  The actual and computed spatial dielectric constants $c$ along with
the $L^{2}$ relative difference of the computed $\mathbf{v}$ are presented. It
is seen that the computed function $c$ meets the acceptable criteria.
Additionally, Figure (c) provides numerical evidence supporting the rapid
convergence of Algorithm 1. Notably, it can be observed that only 3 iterations
are adequate to achieve an accurate reconstruction. The data in this test is
corrupted with 10% of noise.

In this test, the presence of the square scatterer $\Omega$ is clearly
seeable, with an acceptable reconstructed shape and location at the center of
$\Omega$. Regarding the reconstructed value, the true value of $c$ in the
square scatterer is 2.5, while the maximum value of $c$ calculated within the
reconstructed scatterer reaches 2.34. The relative error is 6.4%. Similarly to
Test 1 and Test 2, Figure 5c plays the role of the evidence for the
convergence behavior detailed in the theorem.

Test 4. We examine a test involving ”fully occluded” scatterers. We consider
two scatterers: one looking like a ring and the other resembling a disk placed
within the first. The true function $c$ is defined as follows:

$c^{\rm
true}(x,y)=\left\\{\begin{array}[]{ll}1.5&0.5^{2}<x^{2}+y^{2}<0.7^{2}\mbox{ or
}x^{2}+y^{2}<0.2^{2}\\\ 1&\mbox{otherwise}.\end{array}\right.$ (4.9)

In this test, the wave number $k$ is $4\pi.$ The true and constructed
spatially distributed dielectric constants are displayed in Figure 6.

(a) The function $c^{\rm true}$

(b) The function $c^{\rm comp}$

(c)  The consecutive difference $\frac{\|{\bf v}^{(p)}-{\bf
v}^{(p-1)}\|_{L^{2}(\Omega)}}{\|{\bf v}^{(p)}\|_{L^{2}(\Omega)}}$

Figure 6:  The actual and computed spatial dielectric constants $c$ along with
the $L^{2}$ relative difference of the computed $\mathbf{v}$ are presented. It
is seen that the computed function $c$ meets the acceptable criteria.
Additionally, Figure (c) provides numerical evidence supporting the rapid
convergence of Algorithm 1. Notably, it can be observed that only 3 iterations
are adequate to achieve an accurate reconstruction. The data in this test is
corrupted with 10% of noise.

Observing both the outer and inner scatterers proves intriguing. Within the
ring-shaped region, the true value of $c$ is known to be 1.5, while the
maximum reconstructed value reaches 1.31, resulting in a relative error of
12.67%. Similarly, within the inner scatterer, the true value of $c$ remains
at 1.5, with the maximum reconstructed value reaching 1.30, yielding a
relative error of 13.33%. Notably, while this test presents a significant
challenge due to the need to detect one inclusion within another, the relative
errors remain consistent with the 10% noise level. Once again, the convergence
of our algorithm is verified both by Theorem 3.1 and by the observations
presented in Figure 6c.

##  5 Further discussions and concluding remarks

In the first part of this section, we delve into the computational cost of the
Carleman contraction mapping method in Algorithm 1. Remarkably, the
computational expense of the Carleman contraction mapping method employed in
this study is notably modest. In particular, we only need to compute solutions
for several linear systems to obtain the desired solutions to nonlinear
problems in the tests above. As evident from Figures 3c–6c, although solving 3
linear systems suffices, we execute this process 10 times to demonstrate
convergence. Moreover, it is worth highlighting that the Carleman contraction
mapping method notably outperforms the Carleman convexification method
introduced in our prior works such as [23, 24, 25]. The computational Carleman
convexification method incurs significant expense due to its involvement in
the intensive procedure of minimizing a cost functional using the gradient
descent method. On the other hand, the Carleman contraction mapping method
does not require a precise initial guess of the true solution to the inverse
scattering problem. Consequently, we assert that the Carleman contraction
mapping method is not only fast but also globally convergent.

The inverse scattering problem is widely recognized as severely ill-posed.
Given the high noise levels typically encountered, one might question how we
achieve satisfactory numerical results. The primary reason for our success
lies in the truncation technique employed in (2.8), which renders the original
inverse scattering problem distinct from the approximation model (2.16). An
intriguing aspect emerges from this disparity: we effectively trade off
equivalence for the challenge of solving a system of elliptic equations with
Cauchy boundary data, a problem known for its stability. Consequently, our
method demonstrates reduced sensitivity to noise, thereby facilitating robust
numerical solutions.

In summary, this paper focuses on addressing the inverse scattering problem,
which involves using waves generated through the interaction of plane waves
with various directions and scatterers to detect these scatterers. These
scatterers are characterized by their dielectric constant $c$. Our approach to
solving the inverse scattering problem comprises two key steps. Initially, we
eliminate the unknown coefficient $c$ from the governing equation, resulting
in a system of partial differential equations. Subsequently, we employ the
Carleman contraction mapping method to solve this system effectively. Once the
solution to the aforementioned system is obtained, we can then compute the
solution to the governing equation, thereby yielding the solution to the
inverse scattering problem.

## Declarations

This manuscript is an original work. It has not been published previously and
is not under consideration for publication elsewhere. Both authors have
contributed significantly and equally to the paper.

## Acknowledgement

This work was partially supported by National Science Foundation grant
DMS-2208159.

## References

  * [1] R. Abney, T. T. Le, L. H. Nguyen, and C. Peters.  A Carleman-Picard approach for reconstructing zero-order coefficients in parabolic equations with limited data.  preprint ArXiv:2309:14599, 2023.
  * [2] H. Ammari, Y. Chow, and J. Zou.  Phased and phaseless domain reconstruction in inverse scattering problem via scattering coefficients.  SIAM J. Appl. Math., 76:1000–1030, 2016.
  * [3] H. Ammari and H. Kang.  Reconstruction of Small Inhomogeneities from Boundary Measurements, volume 1846.  Lecture Notes in Mathematics, Springer, 2004.
  * [4] A. B. Bakushinskii and M. Y. Kokurin.  Iterative Methods for Approximate Solutions of Inverse Problems.  Springer, New York, 2004.
  * [5] G. Bao and P. Li.  Inverse medium scattering for the Helmholtz equation at fixed frequency.  Inverse Problems, 21:1621–1641, 2005.
  * [6] G. Bao and P. Li.  Inverse medium scattering problems for electromagnetic waves.  SIAM J. Appl. Math., 65:2049–2066, 2005.
  * [7] Gang Bao, P. Li, J. Lin, and F. Triki.  Inverse scattering problems with multi-frequencies.  Inverse Problems, 31:093001, 2015.
  * [8] L. Beilina and M. V. Klibanov.  Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems.  Springer, New York, 2012.
  * [9] N. Bleistein.  Mathematical Methods for Wave Phenomena.  Academic Press, Orlando, 1984.
  * [10] T. Bui-Thanh and O. Ghattas.  Analysis of the Hessian for inverse scattering problems: I. Inverse shape scattering of acoustic waves.  Inverse Problems, 28:055001, 2012.
  * [11] T. Bui-Thanh and O. Ghattas.  Analysis of the Hessian for inverse scattering problems: II. Inverse medium scattering of acoustic waves.  Inverse Problems, 28:055002, 2012.
  * [12] A. L. Bukhgeim and M. V. Klibanov.  Uniqueness in the large of a class of multidimensional inverse problems.  Soviet Math. Doklady, 17:244–247, 1981.
  * [13] M. Burger and S. Osher.  A survey on level set methods for inverse problems and optimal design.  European J. of Appl. Math., 16:263–301, 2005.
  * [14] G. Chavent.  Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications, Scientic Computation.  Springer, New York, 2009.
  * [15] Y. Chen.  Inverse scattering via Heisenberg’s uncertainty principle.  Inverse Problems, 13:253–282, 1997.
  * [16] W. Chew.  Waves and Fields in Inhomogeneous Media.  Van Nostrand Reinhold, New York, 1990.
  * [17] D. Colton and A. Kirsch.  A simple method for solving inverse scattering problems in the resonance region.  Inverse Problems, 12:383–393, 1996.
  * [18] David Colton and Rainer Kress.  Inverse acoustic and electromagnetic scattering theory. Applied Mathematical Sciences.  Springer, New York, 3rd edition, 2013.
  * [19] A. J. Devaney.  Mathematical Foundations of Imaging, Tomography and Wavefield Inversion.  Cambridge University Press, Cambridge, 2012.
  * [20] H. W. Engl, M. Hanke, and A. Neubauer.  Regularization of Inverse Problems, Mathematics and its Applications.  Kluwer Academic Publishers Group, Dordrecht, 1996.
  * [21] A. V. Goncharsky and S. Y. Romanov.  Supercomputer technologies in inverse problems of ultrasound tomography.  Inverse Problems, 29:075004, 2013.
  * [22] I. Harris and D-L. Nguyen.  Orthogonality sampling method for the electromagnetic inverse scattering problem.  SIAM Journal on Scientific Computing, 42:B722–B737, 2020.
  * [23] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  Convexification and experimental data for a 3D inverse scattering problem with the moving point source.  Inverse Problems, 36:085007, 2020.
  * [24] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  An inverse problem of a simultaneous reconstruction of the dielectric constant and conductivity from experimental backscattering data.  Inverse Problems in Science and Engineering, 29(5):712–735, 2021\.
  * [25] V. A. Khoa, M. V. Klibanov, and L. H. Nguyen.  Convexification for a 3D inverse scattering problem with the moving point source.  SIAM J. Imaging Sci., 13(2):871–904, 2020.
  * [26] A. Kirsch.  Characterization of the shape of a scattering obstacle using the spectral data of the far field operator.  Inverse Problems, 14:1489–1512, 1998.
  * [27] A. Kirsch.  Remarks on the Born approximation and the factorization method.  Appl. Anal., 96:70–84, 2017.
  * [28] M. V. Klibanov.  Convexification of restricted Dirichlet to Neumann map.  J. Inverse and Ill-Posed Problems, 25(5):669–685, 2017.
  * [29] M. V. Klibanov and O. V. Ioussoupova.  Uniform strict convexity of a cost functional for three-dimensional inverse scattering problem.  SIAM J. Math. Anal., 26:147–179, 1995.
  * [30] M. V. Klibanov, T. T. Le, L. H. Nguyen, A. Sullivan, and L. Nguyen.  Convexification-based globally convergent numerical method for a 1D coefficient inverse problem with experimental data.  Inverse Problems and Imaging, 16:1579–1618, 2022.
  * [31] M. V. Klibanov and J. Li.  Inverse Problems and Carleman Estimates: Global Uniqueness, Global Convergence and Experimental Data.  De Gruyter, 2021.
  * [32] M. V. Klibanov, D-L Nguyen, and L. H. Nguyen.  A coefficient inverse problem with a single measurement of phaseless scattering data.  SIAM Journal of Applied Mathematics, 79:1–27, 2019.
  * [33] M. V. Klibanov, L. H. Nguyen, and H. V. Tran.  Numerical viscosity solutions to Hamilton-Jacobi equations via a Carleman estimate and the convexification method.  Journal of Computational Physics, 451:110828, 2022.
  * [34] M. V. Klibanov and V. G. Romanov.  A hölder stability estimate for a 3d coefficient inverse problem for a hyperbolic equation with a plane wave.  Journal of Inverse and Ill-posed Problems, 31:223–242, 2023.
  * [35] K. J. Langenberg.  Applied inverse problems for acoustic, electromagnetic and elastic wave scattering.  In Sabatier, editor, Basic Methods of Tomography and Inverse Problems, pages 127–467. Adam Hilger, 1987.
  * [36] M. M. Lavrent’ev, V. G. Romanov, and S. P. Shishat$\cdot$skiĭ.  Ill-Posed Problems of Mathematical Physics and Analysis.  Translations of Mathematical Monographs. AMS, Providence: RI, 1986.
  * [37] T. T. Le.  Global reconstruction of initial conditions of nonlinear parabolic equations via the Carleman-contraction method.  In D-L. Nguyen, L. H. Nguyen, and T-P. Nguyen, editors, Advances in Inverse problems for Partial Differential Equations, volume 784 of Contemporary Mathematics, pages 23–42. American Mathematical Society, 2023.
  * [38] T. T. Le and L. H. Nguyen.  A convergent numerical method to recover the initial condition of nonlinear parabolic equations from lateral Cauchy data.  Journal of Inverse and Ill-posed Problems, 30(2):265–286, 2022\.
  * [39] T. T. Le, L. H. Nguyen, and H. V. Tran.  A Carleman-based numerical method for quasilinear elliptic equations with over-determined boundary data and applications.  Computers and Mathematics with Applications, 125:13–24, 2022.
  * [40] T. T. Le, L. V. Nguyen, L. H. Nguyen, and H. Park.  The time dimensional reduction method to determine the initial conditions without the knowledge of damping coefficients.  preprint arXiv:2308.13152, 2023.
  * [41] J. Li, H. Liu, and Q. Wang.  Enhanced multilevel linear sampling methods for inverse scattering problems.  J. Comput. Phys., 257:554–571, 2014.
  * [42] Z. Li, H. Liu, and J. Zou.  Locating multiple multiscale acoustic scatterers.  SIAM Multiscale Model. Simul., 12:927–952, 2014.
  * [43] S. Moskow and J. Schotland.  Convergence and stability of the inverse Born series for diffuse waves.  Inverse Problems, 24:065004, 2008.
  * [44] D-L Nguyen, L. H. Nguyen, and T. Truong.  The Carleman-based contraction principle to reconstruct the potential of nonlinear hyperbolic equations.  Computers and Mathematics with Applications, 128:239–248, 2022\.
  * [45] D-L. Nguyen and T. Truong.  Imaging of bi-anisotropic periodic structures from electromagnetic near field data.  Journal of Inverse and Ill-posed Problems, 30:205–219, 2022.
  * [46] H. M. Nguyen and L. H. Nguyen.  Cloaking using complementary media for the Helmholtz equation and a three spheres inequality for second order elliptic equations.  Transaction of the American Mathematical Society, 2:93–112, 2015\.
  * [47] L. H. Nguyen.  An inverse space-dependent source problem for hyperbolic equations and the Lipschitz-like convergence of the quasi-reversibility method.  Inverse Problems, 35:035007, 2019.
  * [48] L. H. Nguyen.  A new algorithm to determine the creation or depletion term of parabolic equations from boundary measurements.  Computers and Mathematics with Applications, 80:2135–2149, 2020\.
  * [49] L. H. Nguyen.  The Carleman contraction mapping method for quasilinear elliptic equations with over-determined boundary data.  Acta Mathematica Vietnamica, 48:401–422, 2023.
  * [50] L. H. Nguyen, Q. Li, and M. V. Klibanov.  A convergent numerical method for a multi-frequency inverse source problem in inhomogenous media.  Inverse Problems and Imaging, 13:1067–1094, 2019.
  * [51] P. M. Nguyen and L. H. Nguyen.  A numerical method for an inverse source problem for parabolic equations and its application to a coefficient inverse problem.  Journal of Inverse and Ill-posed Problems, 38:232–339, 2020.
  * [52] M. Soumekh.  Synthetic Aperture Radar Signal Processing with MATLAB Algorithms (https://www.mathworks.com/matlabcentral/fileexchange/2188-synthetic-aperture-radar-signal-processing-with-matlab-algorithms), MATLAB Central File Exchange. John Wiley & Sons, 1999.
  * [53] A. N. Tikhonov, A. Goncharsky, V. V. Stepanov, and A. G. Yagola.  Numerical Methods for the Solution of Ill-Posed Problems.  Kluwer Academic Publishers Group, Dordrecht, 1995.