# Recovering coefficients in a system of semilinear Helmholtz equations from
internal data

Kui Ren  Department of Applied Physics and Applied Mathematics, Columbia
University, New York, NY 10027<EMAIL_ADDRESS>Nathan Soedjak  Department
of Applied Physics and Applied Mathematics, Columbia University, New York, NY
10027<EMAIL_ADDRESS>
###### Abstract

We study an inverse problem for a coupled system of semilinear Helmholtz
equations where we are interested in reconstructing multiple coefficients in
the system from internal data measured in applications such as thermoacoustic
imaging. We derive results on the uniqueness and stability of the inverse
problem in the case of small boundary data based on the technique of first-
and higher-order linearization. Numerical simulations are provided to
illustrate the quality of reconstructions that can be expected from noisy
data.

Key words. Inverse problems, semilinear Helmholtz equation, uniqueness,
stability, thermoacoustic imaging, nonlinear physics, second harmonic
generation

AMS subject classifications 2020. 35R30, 49M41, 65N21, 78A46

##  1 Introduction

Let $\Omega\subset\mathbb{R}^{d}$ ($d\geq 2$) be a bounded domain with smooth
boundary $\partial\Omega$. We consider the following system of coupled
semilinear Helmholtz equations

$\begin{array}[]{rcll}\Delta u+k^{2}(1+\eta)u+ik\sigma u&=&-k^{2}\gamma
u^{*}v,&\mbox{in}\ \ \Omega\\\ \Delta v+(2k)^{2}(1+\eta)v+i2k\sigma
v&=&-(2k)^{2}\gamma u^{2},&\mbox{in}\ \ \Omega\\\ u=g,&&v=h,&\mbox{on}\
\partial\Omega\end{array}$ (1)

where $\Delta$ denotes the standard Laplacian operator, and $u^{*}$ denotes
the complex conjugate of $u$. This system serves as a simplified model of the
second harmonic generation process in a heterogeneous medium excited by an
incident wave source $g$ [15, 17, 14, 37, 35, 23, 38]. The fields $u$ and $v$
are, respectively, the incident field (with wave number $k$) and the generated
second-harmonics (with wave number $2k$). The medium has first- and second-
order susceptibility $\eta$ and $\gamma$, respectively, and an absorption
coefficient $\sigma$.

We are interested in inverse problems to system (1) where the objective is to
reconstruct the coefficients in the system from data of the form:

$H(\mathbf{x})=\Gamma(\mathbf{x})\sigma(\mathbf{x})(|u|^{2}+|v|^{2}),\qquad\mathbf{x}\in\Omega\,.$
(2)

where $\Gamma(\mathbf{x})$ is an additional physical coefficient that appears
in the data generation process. This inverse problem is motivated by
applications in thermoacoustic imaging, a hybrid imaging modality where
thermoacoustic effect is used to couple high-resolution ultrasound imaging to
microwave imaging to achieve high-resolution and high-contrast imaging of
physical properties of heterogeneous media in the microwave regime. In
thermoacoustic imaging, $H$ is the initial pressure field of the ultrasound
generated by the thermoacoustic effect. It is proportional to the local energy
absorbed by the medium from microwave illumination, that is,
$\sigma(|u|^{2}+|v|^{2})$. The proportional constant $\Gamma(\mathbf{x})$ is
called the Grüneisen coefficient [10]. We refer interested readers to [12, 13,
20, 3] and references therein for the recent development in the modeling and
computational aspects of thermoacoustic imaging.

There are two main differences between the inverse problem we study here and
those that exist in the literature. First, our model (1) takes into account
second-harmonic generation, a nonlinear mechanism that is often used for the
imaging of molecular signatures of particular proteins in biomedical
applications. Second, the objective of our inverse problem includes the
Grüneisen coefficient $\Gamma$, which is mostly ignored in the previous
studies of quantitative thermoacoustic imaging [2, 6, 12, 11].

The fact that the absorbed energy is in the form of
$\sigma(\mathbf{x})(|u|^{2}+|v|^{2})$ has to be understood from the physics of
the thermoacoustic process. In a nutshell, consider the full time-dependent
forms of the incident (at frequency $\omega$) and generated (at frequency
$2\omega$) electric wave of the form:

$\displaystyle u(\mathbf{x},t)$ $\displaystyle=2\Re(u(\mathbf{x})e^{-i\omega
t})=2|u(\mathbf{x})|\cos(\varphi_{u}(\mathbf{x})-\omega t)\,,$ $\displaystyle
v(\mathbf{x},t)$ $\displaystyle=2\Re(v(\mathbf{x})e^{-i2\omega
t})=2|v(\mathbf{x})|\cos(\varphi_{v}(\mathbf{x})-2\omega t)\,,$

where $\varphi_{u}$ (resp. $\varphi_{v}$) is the phase of $u$ (resp. $v$). Let
$I(\mathbf{x})$ denote the energy density of the total electric field at the
location $\mathbf{x}$, averaged over a period of length $T:=2\pi/\omega$. It
is then clear that

$\displaystyle I(\mathbf{x})$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\frac{1}{2}|u(\mathbf{x},t)+v(\mathbf{x},t)|^{2}\,dt=\frac{1}{2T}\int_{0}^{T}u(\mathbf{x},t)^{2}+v(\mathbf{x},t)^{2}+2u(\mathbf{x},t)v(\mathbf{x},t)\,dt$
$\displaystyle=4|u(\mathbf{x})|^{2}\frac{1}{2T}\int_{0}^{T}\cos^{2}(\varphi_{u}(\mathbf{x})-\omega
t)\,dt+4|v(\mathbf{x})|^{2}\frac{1}{2T}\int_{0}^{T}\cos^{2}(\varphi_{v}(\mathbf{x})-2\omega
t)\,dt$
$\displaystyle\qquad+8|u(\mathbf{x})||v(\mathbf{x})|\frac{1}{2T}\int_{0}^{T}\cos(\varphi_{u}(\mathbf{x})-\omega
t)\cos(\varphi_{v}(\mathbf{x})-2\omega t)\,dt$
$\displaystyle=(|u(\mathbf{x})|^{2}+|v(\mathbf{x})|^{2})\,,$

where we have used the standard trigonometic identity
$\cos(x)\cos(y)=\frac{1}{2}(\cos(x+y)+\cos(x-y))$ to simplify the integrals.
Therefore, the cross-term vanishes due to orthogonality. The absorbed
radiation at location $\mathbf{x}$ is thus
$\sigma(\mathbf{x})I(\mathbf{x})=\sigma(\mathbf{x})(|u|^{2}+|v|^{2})$. This
simple calculation provides a (maybe overly-simplified) justification of the
data (2) as the internal data in thermoacoustic imaging with second-harmonic
generation.

The main objective of this paper is to study the problem of determining
information on $(\Gamma,\eta,\gamma,\sigma)$ from information encoded in the
map:

$\Lambda_{\Gamma,\eta,\gamma,\sigma}:(g,h)\mapsto H.$ (3)

We will show that under appropriate conditions, the data (3) allow unique (and
stable, in an appropriate sense) reconstruction of the coefficients
$(\Gamma,\eta,\gamma,\sigma)$. Moreover, there is an explicit reconstruction
method to recover $(\Gamma,\eta,\sigma)$ (see the proof of Theorem 3.1), and
another explicit method to reconstruct $\gamma$ (see the remarks below (30)).

The paper is organized as follows. We first review in Section 2 some of the
elementary properties of the model (1) that we will use in our analysis. We
also introduce the multilinearization method as the basis of the study of the
inverse problems. We then derive the uniqueness and stability of
reconstructing $(\Gamma,\eta,\sigma)$ in Section 3 and study the problem of
reconstructing $\gamma$ in Section 4. Numerical simulations based on synthetic
data will be provided in Section 5 to demonstrate the quality of the
reconstructions that can be achieved in such an inverse problem before we
conclude the paper with additional remarks in Section 6.

##  2 The forward model and its linearization

Throughout the paper, we make the following assumptions on the domain $\Omega$
and the physical coefficients involved in the inverse problem:

  1. ($\mathcal{A}$-i)

The domain $\Omega$ is bounded with smooth boundary $\partial\Omega$.

  2. ($\mathcal{A}$-ii)

The coefficients $\Gamma,\eta,\sigma,\gamma$ all lie in the set

$\displaystyle\mathcal{Y}:=\\{f\in\mathcal{C}^{2}(\overline{\Omega};\mathbb{R}):c_{1}\leq
f\leq c_{2}\ \mbox{in}\ \overline{\Omega}\\}$

for some $c_{1}>0$ and $c_{2}>0$.

While it is clear that such assumptions can be slightly relaxed for the
technical results in the rest of the paper to still hold, we choose the
current form to make the presentation of the paper easy to follow.

###  2.1 Well-posedness of the forward model

We start with the well-posedness of the semilinear system (1) for small
boundary data.

######  Theorem 2.1.

Let $\alpha\in(0,1)$. Under the assumptions in ($\mathcal{A}$-i) and
($\mathcal{A}$-ii), there exist $\varepsilon>0$ and $\delta>0$ such that for
all $g,h\in\mathcal{C}^{2,\alpha}(\partial\Omega;\mathbb{C})$ with
$\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}<\varepsilon$ and
$\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}<\varepsilon$, the boundary
value problem (1) has a unique solution

$\displaystyle(u,v)\in\\{f\in\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C}):\|f\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq\delta\\}^{2}.$

Moreover, there exists a constant $C=C(\alpha,\Omega,\eta,\sigma,\gamma)$ such
that this unique solution satisfies the estimates

$\displaystyle\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C(\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}),$
(4) $\displaystyle\|v\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C(\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}).$

This result comes as a more-or-less straightforward application of the Banach
fixed point theorem in a standard setting. For the convenience of the readers,
we provide the proof in Appendix A.

The above well-posedness result is not satisfactory as it requires that the
boundary data to be small. Currently, we do not have a stronger result. This
result, however, is sufficient for the inverse problem we want to study as our
method in the next sections will be mainly based on the linearization of the
forward model with small boundary data.

In the engineering literature, it is often the case that one drops the $\gamma
u^{*}v$ term in the first equation of system (1). In this case, the system is
only one-way coupled. The solution to the first equation only appears in the
second equation as the source term. In such a case, well-posedness of the
system can be easily established for general boundary conditions. The
corresponding inverse problems are also simplified. We will comment more on
this issue in the next sections.

###  2.2 First- and higher-order linearizations

To deal with the challenge caused by the nonlinearity of the forward model
(1), we use the technique of linearization [7, 18, 25, 26, 32, 33, 28, 29, 24,
27]. We now document the linearization process.

For a given small number $\varepsilon>0$, let
$(u_{\varepsilon},v_{\varepsilon})$ be the solution to the system

$\begin{array}[]{rcll}\Delta
u_{\varepsilon}+k^{2}(1+\eta)u_{\varepsilon}+ik\sigma
u_{\varepsilon}&=&-k^{2}\gamma u_{\varepsilon}^{*}v_{\varepsilon},&\mbox{in}\
\ \Omega\\\ \Delta v_{\varepsilon}+(2k)^{2}(1+\eta)v_{\varepsilon}+i2k\sigma
v_{\varepsilon}&=&-(2k)^{2}\gamma u_{\varepsilon}^{2},&\mbox{in}\ \
\Omega\end{array}$ (5)

with boundary conditions

$u_{\varepsilon}=\varepsilon g_{1}+\dfrac{1}{2}\varepsilon^{2}g_{2},\qquad
v_{\varepsilon}=\varepsilon
h_{1}+\dfrac{1}{2}\varepsilon^{2}h_{2},\qquad\mbox{on}\ \ \partial\Omega\,.$
(6)

We denote by $(u_{0},v_{0})=(0,0)$ the solution for the case of
$\varepsilon=0$, and by $H_{\varepsilon}$ be the data of the form (2)
corresponding to $(u_{\varepsilon},v_{\varepsilon})$, that is

$\displaystyle
H_{\varepsilon}=\Gamma\sigma(|u_{\varepsilon}|^{2}+|v_{\varepsilon}|^{2})\,.$

We expect that the solution $(u_{\varepsilon},v_{\varepsilon})$ varies
sufficiently smoothly with respect to $\varepsilon$ when $\varepsilon$ is
adequately small. Therefore, formally we have expansions of the solution and
the data in the form of:

$\displaystyle u_{\varepsilon}(\mathbf{x})$ $\displaystyle=\varepsilon
u^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}u^{(2)}(\mathbf{x})+o(\varepsilon^{2}),$
(7) $\displaystyle v_{\varepsilon}(\mathbf{x})$ $\displaystyle=\varepsilon
v^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}v^{(2)}(\mathbf{x})+o(\varepsilon^{2}),$
$\displaystyle H_{\varepsilon}(\mathbf{x})$ $\displaystyle=\varepsilon
H^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}H^{(2)}(\mathbf{x})+\frac{1}{6}\varepsilon^{3}H^{(3)}(\mathbf{x})+o(\varepsilon^{3}),$

as $\varepsilon\rightarrow 0$. When this expansion is well-defined, we have
that

$\displaystyle\begin{split}u^{(1)}(\mathbf{x})&:=\lim_{\varepsilon\rightarrow
0}\frac{u_{\varepsilon}(\mathbf{x})}{\varepsilon},\qquad
u^{(2)}(\mathbf{x}):=\lim_{\varepsilon\rightarrow
0}\frac{u_{\varepsilon}(\mathbf{x})-\varepsilon
u^{(1)}(\mathbf{x})}{\frac{1}{2}\varepsilon^{2}},\\\
v^{(1)}(\mathbf{x})&:=\lim_{\varepsilon\rightarrow
0}\frac{v_{\varepsilon}(\mathbf{x})}{\varepsilon},\qquad
v^{(2)}(\mathbf{x}):=\lim_{\varepsilon\rightarrow
0}\frac{v_{\varepsilon}(\mathbf{x})-\varepsilon
v^{(1)}(\mathbf{x})}{\frac{1}{2}\varepsilon^{2}},\\\
H^{(1)}(\mathbf{x})&:=\lim_{\varepsilon\rightarrow
0}\frac{H_{\varepsilon}(\mathbf{x})}{\varepsilon},\qquad
H^{(2)}(\mathbf{x}):=\lim_{\varepsilon\rightarrow
0}\frac{H_{\varepsilon}(\mathbf{x})-\varepsilon
H^{(1)}(\mathbf{x})}{\frac{1}{2}\varepsilon^{2}},\\\
H^{(3)}(\mathbf{x})&:=\lim_{\varepsilon\rightarrow
0}\frac{H_{\varepsilon}(\mathbf{x})-\varepsilon
H^{(1)}(\mathbf{x})-\frac{1}{2}\varepsilon^{2}H^{(2)}(\mathbf{x})}{\frac{1}{6}\varepsilon^{3}}.\end{split}$
(8)

Assuming for the moment that all the derivatives are well-defined,
straightforward formal calculations then show that on the first order, we have
that $(u^{(1)},v^{(1)})$ solves the boundary value problem:

$\begin{array}[]{rcll}\Delta u^{(1)}+k^{2}(1+\eta)u^{(1)}+ik\sigma
u^{(1)}&=&0,&\mbox{in}\ \ \Omega\\\ \Delta
v^{(1)}+(2k)^{2}(1+\eta)v^{(1)}+i2k\sigma v^{(1)}&=&0,&\mbox{in}\ \ \Omega\\\
u^{(1)}=g_{1},\qquad v^{(1)}=h_{1},&&&\mbox{on}\ \partial\Omega\end{array}$
(9)

while $H^{(1)}$ satisfies

$H^{(1)}=0\,.$ (10)

On the second order, we can formally verify that $(u^{(2)},v^{(2)})$ solves
the boundary value problem:

$\begin{array}[]{rcll}\Delta u^{(2)}+k^{2}(1+\eta)u^{(2)}+ik\sigma
u^{(2)}&=&-2k^{2}\gamma u^{(1)*}v^{(1)},&\mbox{in}\ \ \Omega\\\ \Delta
v^{(2)}+(2k)^{2}(1+\eta)v^{(2)}+i2k\sigma
v^{(2)}&=&-2(2k)^{2}\gamma(u^{(1)})^{2},&\mbox{in}\ \ \Omega\\\
u^{(2)}=g_{2},&&v^{(2)}=h_{2},&\mbox{on}\ \partial\Omega\,.\end{array}$ (11)

The corresponding perturbative data $H^{(2)}$ can be expressed as

$H^{(2)}=2\Gamma\sigma(u^{(1)*}u^{(1)}+v^{(1)*}v^{(1)})\,.$ (12)

A little more algebra shows that the third-order data perturbation is in the
form:

$\displaystyle
H^{(3)}=3\Gamma\sigma\left(u^{(1)*}u^{(2)}+u^{(1)}u^{(2)*}+v^{(1)*}v^{(2)}+v^{(1)}v^{(2)*}\right)\,.$
(13)

The whole linearization process can be justified mathematically. We summarize
the result here.

######  Theorem 2.2.

Let $\alpha\in(0,1)$ and
$g_{1},g_{2}\in\mathcal{C}^{2,\alpha}(\partial\Omega;\mathbb{C})$. For
sufficiently small $\varepsilon$, let $(u_{\varepsilon},v_{\varepsilon})$
denote the unique small solution in
$\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})\times\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$
to the system (5). Then the derivatives (8) are all well-defined. Moreover,
(9), (10), (11), (12), and (13) hold.

The proof of this differentiability result is provided in Appendix B.

The multilinearization procedure we outlined here is quite standard. It has
been improved by many authors and utilized to solve various inverse problems
to nonlinear models; see, for instance, [22, 28, 30, 31, 36] and references
therein for some examples of such results.

##  3 The reconstruction of ($\Gamma$, $\sigma$, $\eta$)

The first inverse problem is therefore to reconstruct $(\Gamma,\sigma,\eta)$
from the data $H^{(2)}$ in (12) with the model for $u^{(1)}$ and $v^{(1)}$
given in (9). By taking $h_{1}=0$, the problem reduces to reconstructing
$(\Gamma,\sigma,\eta)$ from the data

$H^{(2)}=2\Gamma\sigma u^{(1)*}u^{(1)},$ (14)

with the model for $u^{(1)}$ given in (9).

When $\Gamma$ and $\eta$ are known, this problem was analyzed in [6, 11]. It
was shown that $\sigma$ can be uniquely recovered with a fixed point
iteration. More precisely, for the model

$\begin{array}[]{rcll}\Delta
u+k^{2}(1+\eta(\mathbf{x}))u+ik\sigma(\mathbf{x})u&=&0,&\mbox{in}\ \Omega\\\
u&=&g,&\mbox{on}\ \partial\Omega\end{array}$ (15)

with internal data

$H(\mathbf{x})=\Gamma\sigma(\mathbf{x})|u|^{2}\,,$ (16)

it is shown in [11] that when $\Gamma$ and $\eta$ are known, one can
reconstruct $\sigma$ uniquely and stably (in appropriate metrics) from one
dataset, provided that the boundary illumination $g$ is appropriately chosen.
(More specifically, the proof requires that $g$ is sufficiently close to a a
function of the form $e^{\rho\cdot x}|_{\partial\Omega}$, for some
$\rho\in\mathbb{C}^{n}$ with $\rho\cdot\rho=0$ and $|\rho|$ sufficiently
large.)

In [6], an explicit procedure for reconstructing $\sigma$ is given (again,
assuming that $\eta$ and $\Gamma$ are known). Here, we modify the method in
order to deal with the case of unknown refractive index $\eta$ and unknown
Grüneisen coefficient $\Gamma$. We use the procedure to develop a uniqueness
and stability result from two well-chosen datasets.

Let $g_{1}$ and $g_{2}$ be two incident sources. We measure data corresponding
to the illuminations $g_{1}+g_{2}$ and $g_{1}+ig_{2}$ in addition to those
corresponding to $g_{1}$ and $g_{2}$. The linearity of (15) means that
solutions corresponding to $g_{1}+g_{2}$ and $g_{1}+ig_{2}$ are $u_{1}+u_{2}$
and $u_{1}+iu_{2}$ respectively. The corresponding data are
$\Gamma\sigma|u_{1}+u_{2}|^{2}$ and $\Gamma\sigma|u_{1}+iu_{2}|^{2}$
respectively. We may now apply the _polarization identity_ to get:

$\displaystyle
u_{1}u_{2}^{*}=\frac{1}{2}\left(|u_{1}+u_{2}|^{2}+i|u_{1}+iu_{2}|^{2}-(1+i)|u_{1}|^{2}-(1+i)|u_{2}|^{2}\right)$

on the inner product space $\mathbb{C}$. This gives us that the quantity

$\displaystyle\Gamma\sigma
u_{1}u_{2}^{*}=\frac{1}{2}\left(\Gamma\sigma|u_{1}+u_{2}|^{2}+i\Gamma\sigma|u_{1}+iu_{2}|^{2}-(1+i)\Gamma\sigma|u_{1}|^{2}-(1+i)\Gamma\sigma|u_{2}|^{2}\right)$

is known.

Henceforth, given illuminations $\\{g_{j}\\}_{j=1}^{2}$, we can reconstruct
from the measured internal data the new data:

$E_{j}=\Gamma\sigma u_{j}u_{1}^{*},$ (17)

for $j=1,2$.

The above construction can be used to develop a uniqueness result
straightforwardly.

######  Theorem 3.1.

Let $\\{g_{j}\\}_{j=1}^{2}$ be a set of incident source functions, and suppose
that the measured data $\\{E_{j}\\}_{j=1}^{2}$ satisfy the following two
conditions:

  1. ($\mathcal{B}$-i)

$E_{1}(\mathbf{x})\geq\alpha_{0}>0$ for some $\alpha_{0}$, a.e.
$\mathbf{x}\in\Omega$.

  2. ($\mathcal{B}$-ii)

The vector field

$\beta(\mathbf{x}):=\nabla\frac{E_{2}}{E_{1}}\,$

is at least $W^{1,\infty}$, and $|\beta|\geq\beta_{0}>0$ for some $\beta_{0}$,
a.e. $\mathbf{x}\in\Omega$.

Then $\Gamma$, $\eta$, and $\sigma$ are uniquely determined from the data
$\\{E_{j}\\}_{j=1}^{2}$.

###### Proof.

We follow the procedures developed in [6, 9]. We multiply the equation for
$u_{1}$ by $u_{2}$ and multiply the equation for $u_{2}$ by $u_{1}$. We
subtract the results to have

$u_{1}\Delta u_{2}-u_{2}\Delta u_{1}=0\,.$

We can then rewrite this into

$\nabla\cdot u_{1}^{2}\nabla\frac{u_{2}}{u_{1}}=\nabla\cdot
u_{1}^{2}\nabla\frac{E_{2}}{E_{1}}=\nabla\cdot(u_{1}^{2}\beta)=0\,.$

The vector field $\beta$ is known from the data. Therefore the above equation
is a transport equation for $u_{1}$, that is,

$\nabla\cdot(u_{1}^{2}\beta)=0,\ \ \mbox{in}\ \ \Omega,\qquad u_{1}=g_{1},\ \
\mbox{on}\ \ \partial\Omega\,.$ (18)

With the assumption in ($\mathcal{B}$-$ii$), classical results in [5, 16, 19,
21] show that there exists a unique weak solution $u_{1}$ to (18). This gives
us the unique reconstruction of $u_{1}$.

Now that we have reconstructed $u_{1}$, we can use the equation (15) to
reconstruct the potential $q$:

$q(\mathbf{x}):=k^{2}(1+\eta)+ik\sigma=-\frac{\Delta u_{1}}{u_{1}}\,.$ (19)

This gives us $\eta$ and $\sigma$ (which are obtained by taking real and
imaginary parts of $q$). The last step is to reconstruct $\Gamma$ as

$\Gamma=\frac{H_{1}}{\sigma|u_{1}|^{2}}\,.$ (20)

The proof is complete. ∎

This uniqueness result shows a dramatic difference between the inverse problem
defined by (15) and (16) and a similar inverse problem in quantitative
photoacoustic tomography in [9] where it is show that the multiplicative
coefficient $\Gamma$ causes non-uniqueness in the reconstructions, independent
of the amount of data available.

The proof of the above uniqueness result is constructive in the sense that it
provides an explicit way to solve the inverse problem: solving (18) for
$u_{1}$, computing $q$ using (19) and then computing $\Gamma$ as in (20).

In fact, the above explicit reconstruction procedure also leads to partial
(weighted) stability results for the inverse problem.

######  Theorem 3.2.

Let $E=(E_{1},E_{2})$ and
$\widetilde{E}=(\widetilde{E}_{1},\widetilde{E}_{2})$ be the data
corresponding to the coefficients $(\Gamma,\eta,\sigma)\in\mathcal{Y}^{3}$ and
$(\widetilde{\Gamma},\widetilde{\eta},\widetilde{\sigma})\in\mathcal{Y}^{3}$
respectively, generated from illumination source pair $(g_{1},g_{2})$. Under
the assumption that $E$ and $\widetilde{E}$ satisfy
($\mathcal{B}$-i)-($\mathcal{B}$-ii), we assume further that $g_{1}$ and
$g_{2}$ are selected such that $\frac{E_{2}}{E_{1}}$ is sufficiently small.
Then we have that, for some constants $c>0$,

$\|\Gamma\sigma-\widetilde{\Gamma}\widetilde{\sigma}\|_{L^{2}(\Omega)}\leq
c\Big{(}\|H_{1}-\widetilde{H}_{1}\|_{L^{2}(\Omega)}+\|E-\widetilde{E}\|_{(\mathcal{H}^{2}(\Omega))^{2}}\Big{)}\,.$
(21)

###### Proof.

We first observe that

$\Gamma\sigma-\widetilde{\Gamma}\widetilde{\sigma}=\frac{H_{1}}{|u_{1}|^{2}}-\frac{\widetilde{H}_{1}}{|\widetilde{u}_{1}|^{2}}=\frac{H_{1}-\widetilde{H}_{1}}{|u_{1}|^{2}}+\frac{\widetilde{H}_{1}}{|\widetilde{u}_{1}|^{2}}\frac{|\widetilde{u}_{1}|^{2}-|u_{1}|^{2}}{|u_{1}|^{2}}\,.$

This, together with the Triangle Inequality and the fact that $|u_{k}|^{2}$ is
bounded from below, gives us

$\|\Gamma\sigma-\widetilde{\Gamma}\widetilde{\sigma}\|_{L^{2}(\Omega)}\leq\widetilde{c}\Big{(}\|H_{1}-\widetilde{H}_{1}\|_{L^{2}(\Omega)}+\|\widetilde{u}_{1}^{2}-u_{1}^{2}\|_{L^{2}(\Omega)}\Big{)}\,,$
(22)

for some $\widetilde{c}>0$.

To bound the second term in (22) by the data, let $\xi=u_{1}^{2}$ and
$\widetilde{\xi}=\widetilde{u}_{1}^{2}$. Then we have from the equations
$\nabla\cdot\xi\beta=0$ and $\nabla\cdot\widetilde{\xi}\widetilde{\beta}=0$
that

$\nabla\cdot\Big{(}(\xi-\widetilde{\xi})\beta\Big{)}+\nabla\cdot\Big{(}\widetilde{\xi}(\beta-\widetilde{\beta})\Big{)}=0\,.$
(23)

This can be further rewritten into

$\beta\cdot\nabla(\xi-\widetilde{\xi})=-(\xi-\widetilde{\xi})\nabla\cdot\beta-\nabla\cdot\Big{(}\widetilde{\xi}(\beta-\widetilde{\beta})\Big{)}\,,$

which immediately leads to the bound

$\|\beta\cdot\nabla(\xi-\widetilde{\xi})\|_{L^{2}(\Omega)}\leq\|(\xi-\widetilde{\xi})\nabla\cdot\beta\|_{L^{2}(\Omega)}^{2}+\|\nabla\cdot\Big{(}\xi(\beta-\widetilde{\beta})\Big{)}\|_{L^{2}(\Omega)}\,.$
(24)

With the same algebra, we can derive the bound

$\|\widetilde{\beta}\cdot\nabla(\xi-\widetilde{\xi})\|_{L^{2}(\Omega)}\leq\|(\xi-\widetilde{\xi})\nabla\cdot\widetilde{\beta}\|_{L^{2}(\Omega)}+\|\nabla\cdot\Big{(}\widetilde{\xi}(\beta-\widetilde{\beta})\Big{)}\|_{L^{2}(\Omega)}\,.$
(25)

We now multiply (23) by $(\xi-\widetilde{\xi})^{*}$ to have the equation,
after a little algebra,

$\nabla\cdot\Big{(}\big{|}\xi-\widetilde{\xi}\big{|}^{2}\beta\Big{)}-(\xi-\widetilde{\xi})\beta\cdot\nabla(\xi-\widetilde{\xi})^{*}+(\xi-\widetilde{\xi})^{*}\nabla\cdot\widetilde{\xi}(\beta-\widetilde{\beta})=0\,.$

Integrating this equation against a test function
$\phi\in\mathcal{H}^{1}(\Omega)$ and using integration-by-parts on the last
term lead us to the identity

$\int_{\Omega}\big{|}\xi-\widetilde{\xi}\big{|}^{2}\beta\cdot\nabla\phi
d\mathbf{x}+\int_{\Omega}(\xi-\widetilde{\xi})\beta\cdot\phi\nabla(\xi-\widetilde{\xi})^{*}d\mathbf{x}\\\
+\int_{\Omega}(\xi-\widetilde{\xi})^{*}\widetilde{\xi}(\beta-\widetilde{\beta})\cdot\nabla\phi
d\mathbf{x}+\int_{\Omega}\phi\widetilde{\xi}(\beta-\widetilde{\beta})\cdot\nabla(\xi-\widetilde{\xi})^{*}d\mathbf{x}=0\,.$

To simplify the presentation, we combine the second and the fourth terms in
the equation to have

$\int_{\Omega}\big{|}\xi-\widetilde{\xi}\big{|}^{2}\beta\cdot\nabla\phi
d\mathbf{x}+\int_{\Omega}(\beta\xi-\widetilde{\beta}\widetilde{\xi})\cdot\phi\nabla(\xi-\widetilde{\xi})^{*}d\mathbf{x}+\int_{\Omega}(\xi-\widetilde{\xi})^{*}\widetilde{\xi}(\beta-\widetilde{\beta})\cdot\nabla\phi
d\mathbf{x}=0\,.$

Taking the test function $\phi=\frac{E_{2}^{*}}{E_{1}^{*}}$ (hence
$\nabla\phi=\nabla\frac{E_{2}^{*}}{E_{1}^{*}}=\beta^{*}$), we have

$\int_{\Omega}\big{|}\xi-\widetilde{\xi}\big{|}^{2}|\beta|^{2}d\mathbf{x}=-\int_{\Omega}(\beta\xi-\widetilde{\beta}\widetilde{\xi})\cdot\frac{E_{2}^{*}}{E_{1}^{*}}\nabla(\xi-\widetilde{\xi})^{*}d\mathbf{x}-\int_{\Omega}(\xi-\widetilde{\xi})^{*}\beta^{*}\cdot\widetilde{\xi}(\beta-\widetilde{\beta})d\mathbf{x}\,.$

This gives us the bound

$\|(\xi-\widetilde{\xi})\beta\|_{L^{2}(\Omega)}^{2}\leq\int_{\Omega}\big{|}(\beta\xi-\widetilde{\beta}\widetilde{\xi})\cdot\frac{E_{2}^{*}}{E_{1}^{*}}\nabla(\xi-\widetilde{\xi})^{*}\big{|}d\mathbf{x}+\int_{\Omega}\big{|}(\xi-\widetilde{\xi})^{*}\beta^{*}\cdot\widetilde{\xi}(\beta-\widetilde{\beta})\big{|}d\mathbf{x}\,.$
(26)

The first term on the right-hand side of (26) can be bounded as follows:

$\int_{\Omega}\big{|}(\beta\xi-\widetilde{\beta}\widetilde{\xi})\cdot\frac{E_{2}^{*}}{E_{1}^{*}}\nabla(\xi-\widetilde{\xi})^{*}\big{|}d\mathbf{x}\leq\|\frac{E_{2}^{*}}{E_{1}^{*}}\|_{L^{\infty}(\Omega)}\int_{\Omega}\big{|}(\beta\xi-\widetilde{\beta}\widetilde{\xi})\cdot\nabla(\xi-\widetilde{\xi})^{*}\big{|}d\mathbf{x}\\\
\leq\frac{1}{2}\|\frac{E_{2}^{*}}{E_{1}^{*}}\|_{L^{\infty}(\Omega)}\Big{[}\int_{\Omega}\big{|}\beta\xi\cdot\nabla(\xi-\widetilde{\xi})^{*}\big{|}^{2}d\mathbf{x}+\int_{\Omega}\big{|}\widetilde{\beta}\widetilde{\xi}\cdot\nabla(\xi-\widetilde{\xi})^{*}\big{|}^{2}d\mathbf{x}\Big{]}\\\
\leq
2\|\frac{E_{2}^{*}}{E_{1}^{*}}\|_{L^{\infty}(\Omega)}\Big{[}\|(\xi-\widetilde{\xi})\nabla\cdot\beta\|_{L^{2}(\Omega)}^{2}+\|\nabla\cdot\Big{(}\xi(\beta-\widetilde{\beta})\Big{)}\|_{L^{2}(\Omega)}^{2}\Big{]}\,,$

where we have used (24) and (25) to get the last inequality. The second term
on the right-hand side of (26) can be bounded as:

$\int_{\Omega}\big{|}(\xi-\widetilde{\xi})^{*}\beta^{*}\cdot\widetilde{\xi}(\beta-\widetilde{\beta})\big{|}d\mathbf{x}\leq\frac{1}{2}\Big{[}\frac{1}{\kappa^{2}}\|(\xi-\widetilde{\xi})\beta\|_{L^{2}(\Omega)}^{2}+\kappa^{2}\|(\beta-\widetilde{\beta})\widetilde{\xi}\|_{L^{2}(\Omega)}^{2}\Big{]}\,,$

for any $\kappa>0$.

Under the assumption that $|\frac{E_{2}}{E_{1}}|$ is sufficiently small, we
can take $\kappa$ to be sufficiently large so that (26) now implies that

$\|\xi-\widetilde{\xi}\|_{L^{2}(\Omega)}^{2}\lesssim\|(\beta-\widetilde{\beta})\widetilde{\xi}\|_{L^{2}(\Omega)}^{2}+\|\nabla\cdot\Big{(}\xi(\beta-\widetilde{\beta})\Big{)}\|_{L^{2}(\Omega)}^{2}\,.$
(27)

The next step is to bound $\|\beta-\widetilde{\beta}\|_{L^{2}(\Omega)}$ and
$\|\nabla\cdot\xi\big{(}\beta-\widetilde{\beta}\big{)}\|_{L^{2}(\Omega)}$. To
this end, we use the expansion

$\beta-\widetilde{\beta}=(\widetilde{E}_{1}-E_{1})\nabla\frac{\widetilde{E}_{2}}{E_{1}\widetilde{E}_{1}}+\frac{\widetilde{E}_{2}}{E_{1}\widetilde{E}_{1}}\nabla(\widetilde{E}_{1}-E_{1})\\\
+(E_{2}-\widetilde{E}_{2})\nabla\frac{1}{E_{1}}+\frac{1}{E_{1}^{2}}\nabla(E_{2}-\widetilde{E}_{2})+\frac{1}{E_{1}\widetilde{E}_{1}}(\widetilde{E}_{1}-E_{1})\nabla(E_{2}-\widetilde{E}_{2})\,,$

to derive the bound

$\|\beta-\widetilde{\beta}\|_{L^{2}(\Omega)}\lesssim\|E-\widetilde{E}\|_{(\mathcal{H}^{1}(\Omega))^{2}}\,.$

In a similar manner, we find that

$\|\nabla\cdot\xi\big{(}\beta-\widetilde{\beta}\big{)}\|_{L^{2}(\Omega)}\lesssim\|E-\widetilde{E}\|_{(\mathcal{H}^{2}(\Omega))^{2}}\,.$

Plugging these results into (27) will give us

$\|\xi-\widetilde{\xi}\|_{L^{2}(\Omega)}\lesssim\|E-\widetilde{E}\|_{(\mathcal{H}^{2}(\Omega))^{2}}\,.$

This, together with the bound (22), will lead us to the stability results of
(21). ∎

######  Remark 3.3.

With the standard techniques complex geometrical solutions, one can show that
for every value of the true coefficients
$\eta,\sigma\in\mathcal{H}^{m}(\Omega)$, where $m>1+\frac{d}{2}$, there exists
a set of illuminations $(g_{j})_{j=1}^{d+1}$ such that the corresponding
measured data $(E_{j})_{j=1}^{d+1}$ satisfies both conditions
($\mathcal{B}$-$i$) and ($\mathcal{B}$-$ii$) [6]. In fact, following [4], it
may be possible to ensure that ($\mathcal{B}$-$i$) and ($\mathcal{B}$-$ii$)
hold with high probability by drawing the boundary illuminations
$(g_{j})_{j=1}^{d+1}$ independently at random from a sub-Gaussian distribution
on $\mathcal{H}^{1/2}(\partial\Omega)$.

######  Remark 3.4.

We observe that the above reconstruction procedure also works in the case when
the internal datum is of the form $H=|u|$ (in which case $|u|^{2}$ is known),
that is, the datum is independent of $\Gamma\sigma$.

##  4 The reconstruction of $\gamma$

The remaining problem is to reconstruct $\gamma$ using third-order
perturbation of the data. In the rest of this section, we assume that in
addition to the internal data (3), we also have access to the Dirichlet-to-
Neumann map

$\Pi_{\gamma}:(g,h)\mapsto\left(\dfrac{\partial{u}}{\partial{\nu}}\Big{|}_{\partial\Omega},\dfrac{\partial{v}}{\partial{\nu}}\Big{|}_{\partial\Omega}\right)\equiv(J_{u},J_{v})\,.$
(28)

Note that we omit the dependence of $\Pi$ on $\Gamma$, $\eta$, and $\sigma$
intentionally here since those coefficients are already known.

The multilinearization of $(J_{u},J_{v})$ can be established with the
calculations in Appendix (B). We will directly use the derivatives
$(J_{u}^{(1)},J_{v}^{(1)})$ and $(J_{u}^{(2)},J_{v}^{(2)})$.

Let us recall that the third-order derivative of the data $H^{(3)}$ is given
in (13). This implies that

$\displaystyle
u^{(1)*}u^{(2)}+u^{(1)}u^{(2)*}+v^{(1)*}v^{(2)}+v^{(1)}v^{(2)*}=\frac{H^{(3)}}{3\Gamma\sigma},$

where $(u^{(1)},v^{(1)})$ and $(u^{(2)},v^{(2)})$ are respectively the
solutions to (9) and (11), is known in $\Omega$.

From now on, we set $g_{2}=h_{2}=0$ in (6). Consequently, the system (11) for
$(u^{(2)},v^{(2)})$ reduces to

$\begin{array}[]{rcll}\Delta u^{(2)}+k^{2}(1+\eta)u^{(2)}+ik\sigma
u^{(2)}&=&-2k^{2}\gamma u^{(1)*}v^{(1)},&\mbox{in}\ \ \Omega\\\ \Delta
v^{(2)}+(2k)^{2}(1+\eta)v^{(2)}+i2k\sigma
v^{(2)}&=&-2(2k)^{2}\gamma(u^{(1)})^{2},&\mbox{in}\ \ \Omega\\\
u^{(2)}=0,&&v^{(2)}=0,&\mbox{on}\ \ \partial\Omega\end{array}$ (29)

We can now take the complex conjugate of (29) and leverage the fact that
$\gamma$ is real-valued to write down the following system of _linear_
equations for $(u^{(2)},v^{(2)})$, $(u^{(2)*},v^{(2)*})$, and $\gamma$:

$\begin{array}[]{rcll}(\Delta+q_{1})u^{(2)}+2k^{2}u^{(1)*}v^{(1)}\gamma&=&0,&\mbox{in}\
\Omega\\\
(\Delta+q_{1}^{*})u^{(2)*}+2k^{2}u^{(1)}v^{(1)*}\gamma&=&0,&\mbox{in}\
\Omega\\\ (\Delta+q_{2})v^{(2)}+2(2k)^{2}(u^{(1)})^{2}\gamma&=&0,&\mbox{in}\
\Omega\\\
(\Delta+q_{2}^{*})v^{(2)*}+2(2k)^{2}(u^{(1)*})^{2}\gamma&=&0,&\mbox{in}\
\Omega\\\
u^{(1)*}u^{(2)}+u^{(1)}u^{(2)*}+v^{(1)*}v^{(2)}+v^{(1)}v^{(2)*}&=&\frac{H^{(3)}}{3\Gamma\sigma},&\mbox{in}\
\Omega\\\ (u^{(2)},u^{(2)*},v^{(2)},v^{(2)*})&=&(0,0,0,0),&\mbox{on}\
\partial\Omega\\\
\left(\dfrac{\partial{u^{(2)}}}{\partial{\nu}},\dfrac{\partial{u^{(2)*}}}{\partial{\nu}},\dfrac{\partial{v^{(2)}}}{\partial{\nu}},\dfrac{\partial{v^{(2)*}}}{\partial{\nu}}\right)&=&(J_{u}^{(2)},J_{u}^{(2)*},J_{v}^{(2)},J_{v}^{(2)*}),&\mbox{on}\
\partial\Omega\,,\end{array}$ (30)

where we have used the notation

$\displaystyle q_{1}:=k^{2}(1+\eta)+ik\sigma,\qquad
q_{2}:=(2k)^{2}(1+\eta)+i2k\sigma\,.$

If we can solve (30), we can reconstruct $\gamma$ (and the associated
$(u^{(2)},v^{(2)})$). This is a non-iterative reconstruction method. In the
rest of this section, we show that $\gamma$ can be uniquely reconstructed from
available data by analyzing the uniqueness of the solution to the linear
system (30). The analysis is based on the uniqueness theory for redundant
elliptic systems reviewed in [8], which we summarize briefly in Appendix C for
the convenience of the readers.

######  Theorem 4.1.

Let $X\subset\mathbb{R}^{n}$ be an open set and $\Gamma$, $\eta$, $\sigma$ be
given positive $\mathcal{C}^{2}$ functions on $\overline{X}$. For every
bounded open subset $\Omega\subset X$ with smooth boundary $\partial\Omega$,
we denote by $(\Lambda_{\gamma}^{(\Omega)},\Pi_{\gamma}^{(\Omega)})$ and
$(\Lambda_{\widetilde{\gamma}}^{(\Omega)},\Pi_{\widetilde{\gamma}}^{(\Omega)})$
the data corresponding to $\gamma\in\mathcal{Y}$ and
$\widetilde{\gamma}\in\mathcal{Y}$ respectively. Let $\mathbf{x}_{0}\in X$ be
arbitrary. Then there exists $\epsilon=\epsilon(\eta,\sigma,\mathbf{x}_{0})>0$
such that for every $\Omega\subset B(\mathbf{x}_{0},\epsilon)$, there exist
$g_{1}$ and $h_{1}$ in (6) such that

$(H^{(3)},J_{u}^{(2)},J_{v}^{(2)})=(\widetilde{H}^{(3)},\widetilde{J}_{u}^{(2)},\widetilde{J}_{v}^{(2)})\implies\gamma=\widetilde{\gamma}\,,$

and for all $p>1$,

$\|\gamma-\widetilde{\gamma}\|_{L^{p}(\Omega)}\leq
C\Big{(}\Big{\|}H^{(3)}-\widetilde{H}^{(3)}\Big{\|}_{W^{2,p}(\Omega)}\\\
+\Big{\|}J_{u}^{(2)}-\widetilde{J}_{u}^{(2)}\Big{\|}_{W^{1-1/p,p}(\partial\Omega)}+\Big{\|}J_{v}^{(2)}-\widetilde{J}_{v}^{(2)}\Big{\|}_{W^{1-1/p,p}(\partial\Omega)}\Big{)},$
(31)

for some constant $C=C(\Omega,\Gamma,\eta,\sigma)>0$.

###### Proof.

Let us define

$\displaystyle\widehat{u}:=u-\widetilde{u},\quad\widehat{v}:=v-\widetilde{v},\quad\widehat{\gamma}:=\gamma-\widetilde{\gamma},\quad\widehat{H}:=H-\widetilde{H},\quad\widehat{J}_{u}:=J_{u}-\widetilde{J}_{u},\quad\widehat{J}_{v}:=J_{v}-\widetilde{J}_{v}.$

The linear system (30) then implies that

$\begin{array}[]{rcll}(\Delta+q_{1})\widehat{u}^{(2)}+2k^{2}u^{(1)*}v^{(1)}\widehat{\gamma}&=&0,&\mbox{in}\
\ \Omega\\\
(\Delta+q_{1}^{*})\widehat{u}^{(2)*}+2k^{2}u^{(1)}v^{(1)*}\widehat{\gamma}&=&0,&\mbox{in}\
\ \Omega\\\
(\Delta+q_{2})\widehat{v}^{(2)}+2(2k)^{2}(u^{(1)})^{2}\widehat{\gamma}&=&0,&\mbox{in}\
\ \Omega\\\
(\Delta+q_{2}^{*})\widehat{v}^{(2)*}+2(2k)^{2}(u^{(1)*})^{2}\widehat{\gamma}&=&0,&\mbox{in}\
\Omega\\\
u^{(1)*}\widehat{u}^{(2)}+u^{(1)}\widehat{u}^{(2)*}+v^{(1)*}\widehat{v}^{(2)}+v^{(1)}\widehat{v}^{(2)*}&=&\frac{\widehat{H}^{(3)}}{3\Gamma\sigma},&\mbox{in}\
\ \Omega\\\
(\widehat{u}^{(2)},\widehat{u}^{(2)*},\widehat{v}^{(2)},\widehat{v}^{(2)*})&=&(0,0,0,0),&\mbox{on}\
\ \partial\Omega\\\
\left(\dfrac{\partial{\widehat{u}^{(2)}}}{\partial{\nu}},\dfrac{\partial{\widehat{u}^{(2)*}}}{\partial{\nu}},\dfrac{\partial{\widehat{v}^{(2)}}}{\partial{\nu}},\dfrac{\partial{\widehat{v}^{(2)*}}}{\partial{\nu}}\right)&=&(\widehat{J}_{u}^{(2)},\widehat{J}_{u}^{(2)*},\widehat{J}_{v}^{(2)},\widehat{J}_{v}^{(2)*}),&\mbox{on}\
\ \partial\Omega\,.\end{array}$ (32)

We first eliminate $\widehat{\gamma}$ by plugging in
$\widehat{\gamma}=-\frac{(\Delta+q_{2}^{*})\widehat{v}^{(2)*}}{2(2k)^{2}(u^{(1)*})^{2}}$
from the fourth equation into the first three equations. We then take the
Laplacian of the last equation. These procedures lead us to the linear system

$\displaystyle(\Delta+q_{1})\widehat{u}^{(2)}-\frac{v^{(1)}}{4u^{(1)*}}(\Delta+q_{2}^{*})\widehat{v}^{(2)*}$
$\displaystyle=0,$
$\displaystyle(\Delta+q_{1}^{*})\widehat{u}^{(2)*}-\frac{u^{(1)}v^{(1)*}}{4(u^{(1)*})^{2}}(\Delta+q_{2}^{*})\widehat{v}^{(2)*}$
$\displaystyle=0,$
$\displaystyle(\Delta+q_{2})\widehat{v}^{(2)}-\frac{(u^{(1)})^{2}}{(u^{(1)*})^{2}}(\Delta+q_{2}^{*})\widehat{v}^{(2)*}$
$\displaystyle=0,$
$\displaystyle\Delta(u^{(1)*}\widehat{u}^{(2)})+\Delta(u^{(1)}\widehat{u}^{(2)*})+\Delta(v^{(1)*}\widehat{v}^{(2)})+\Delta(v^{(1)}\widehat{v}^{(2)*})$
$\displaystyle=\Delta\left(\frac{\widehat{H}^{(3)}}{3\Gamma\sigma}\right),$

in the unknowns $\widehat{u}^{(2)}$, $\widehat{u}^{(2)*}$,
$\widehat{v}^{(2)}$, and $\widehat{v}^{(2)*}$. The system may be written in
the following matrix form, for the quantity
$w:=(\widehat{u}^{(2)},\widehat{u}^{(2)*},\widehat{v}^{(2)},\widehat{v}^{(2)*})$:

$\begin{array}[]{rcll}\mathcal{A}w&=&\mathcal{S},&\mbox{in}\ \ \Omega\\\
w&=&(0,0,0,0),&\mbox{on}\ \partial\Omega\\\
\dfrac{\partial{w}}{\partial{\nu}}&=&(\widehat{J}_{u}^{(2)},\widehat{J}_{u}^{(2)*},\widehat{J}_{v}^{(2)},\widehat{J}_{v}^{(2)*}),&\mbox{on}\
\partial\Omega\,,\end{array}$ (33)

where

$\displaystyle\mathcal{A}(\mathbf{x},D):=\begin{pmatrix}\Delta+q_{1}&&&-\frac{v^{(1)}}{4u^{(1)*}}(\Delta+q_{2}^{*})\\\
&\Delta+q_{1}^{*}&&-\frac{u^{(1)}v^{(1)*}}{4(u^{(1)*})^{2}}(\Delta+q_{2}^{*})\\\
&&\Delta+q_{2}&-\frac{(u^{(1)})^{2}}{(u^{(1)*})^{2}}(\Delta+q_{2}^{*})\\\
\Delta(u^{(1)*}\cdot)&\Delta(u^{(1)}\cdot)&\Delta(v^{(1)*}\cdot)&\Delta(v^{(1)}\cdot)\end{pmatrix}\,,$

and

$\displaystyle
S:=\left(0,0,0,\Delta\left(\frac{\widehat{H}^{(3)}}{3\Gamma\sigma}\right)\right)\,.$

In the rest of the proof, we show that $\mathcal{A}$ is an elliptic operator
in the sense of Douglis and Nirenberg [8]. For the convenience of the reader,
we provide a brief review of elliptic system theory in Appendix C. We choose
the Douglis-Nirenberg numbers

$\displaystyle(s_{1},s_{2},s_{3},s_{4})=(0,0,0,0),\qquad(t_{1},t_{2},t_{3},t_{4})=(2,2,2,2).$
(34)

The principal part $\mathcal{A}_{0}(\mathbf{x},D)$ of $\mathcal{A}$ has symbol

$\displaystyle\mathcal{A}_{0}(\mathbf{x},{\boldsymbol{\xi}}):=\begin{pmatrix}|{\boldsymbol{\xi}}|^{2}&&&-\frac{v^{(1)}}{4u^{(1)*}}|{\boldsymbol{\xi}}|^{2}\\\
&|{\boldsymbol{\xi}}|^{2}&&-\frac{u^{(1)}v^{(1)*}}{4(u^{(1)*})^{2}}|{\boldsymbol{\xi}}|^{2}\\\
&&|{\boldsymbol{\xi}}|^{2}&-\frac{(u^{(1)})^{2}}{(u^{(1)*})^{2}}|{\boldsymbol{\xi}}|^{2}\\\
u^{(1)*}|{\boldsymbol{\xi}}|^{2}&u^{(1)}|{\boldsymbol{\xi}}|^{2}&v^{(1)*}|{\boldsymbol{\xi}}|^{2}&v^{(1)}|{\boldsymbol{\xi}}|^{2}\end{pmatrix}\,.$

One readily sees that $\mathcal{A}_{0}(\mathbf{x}_{0},{\boldsymbol{\xi}})$ has
full rank $4$ for all ${\boldsymbol{\xi}}\neq 0$ if and only if the following
condition holds at $\mathbf{x}_{0}$:

$\displaystyle-\frac{v^{(1)}}{4u^{(1)*}}\cdot
u^{(1)*}-\frac{u^{(1)}v^{(1)*}}{4(u^{(1)*})^{2}}\cdot
u^{(1)}-\frac{(u^{(1)})^{2}}{(u^{(1)*})^{2}}\cdot v^{(1)*}\neq v^{(1)},$

or equivalently

$\displaystyle-\frac{(u^{(1)})^{2}}{(u^{(1)*})^{2}}\neq\frac{v^{(1)}}{v^{(1)*}}$

at $\mathbf{x}_{0}$. This condition on $u^{(1)}(\mathbf{x}_{0})$ and
$v^{(1)}(\mathbf{x}_{0})$ is easily achieved by selecting $g_{1}$ and $h_{1}$
appropriately. To be precise, let us consider some ball
$B(\mathbf{x}_{0},\epsilon_{0})\subset X$ and let $u_{0}$ and $v_{0}$ be any
$\mathcal{C}^{2}$ functions on $B(\mathbf{x}_{0},\epsilon_{0})$ satisfying

$\displaystyle\Delta u_{0}+q_{1}u_{0}$ $\displaystyle=0,\quad\mbox{in}\
B(\mathbf{x}_{0},\epsilon_{0}),$ $\displaystyle\Delta v_{0}+q_{2}v_{0}$
$\displaystyle=0,\quad\mbox{in}\ B(\mathbf{x}_{0},\epsilon_{0}),$
$\displaystyle-\frac{u_{0}^{2}}{(u_{0}^{*})^{2}}(\mathbf{x}_{0})$
$\displaystyle\neq\frac{v_{0}}{v_{0}^{*}}(\mathbf{x}_{0}).$

The existence of such $u_{0}$ and $v_{0}$ is obvious as we can take $u_{0}$
and $v_{0}$ to be any solutions to the first two equations and rescale them by
a suitable complex constant to satisfy the condition
$-\frac{u_{0}^{2}}{(u_{0}^{*})^{2}}(\mathbf{x}_{0})\neq\frac{v_{0}}{v_{0}^{*}}(\mathbf{x}_{0})$.
It is also useful to observe that $u_{0}$ and $v_{0}$ depend only on $\eta$
and $\sigma$, not $\gamma$.

Now suppose $\Omega\subset B(\mathbf{x}_{0},\epsilon_{0})$, and select
$g_{1}=u_{0}|_{\partial\Omega}$ and $h_{1}=v_{0}|_{\partial\Omega}$ in (9), so
that

$\displaystyle u^{(1)}=u_{0}|_{\Omega},\qquad v^{(1)}=v_{0}|_{\Omega}\,.$

This means we can write $\mathcal{A}$ as

$\displaystyle\mathcal{A}(\mathbf{x},D):=\begin{pmatrix}\Delta+q_{1}&&&-\frac{v_{0}}{4u_{0}^{*}}(\Delta+q_{2}^{*})\\\
&\Delta+q_{1}^{*}&&-\frac{u_{0}v_{0}^{*}}{4(u_{0}^{*})^{2}}(\Delta+q_{2}^{*})\\\
&&\Delta+q_{2}&-\frac{u_{0}^{2}}{(u_{0}^{*})^{2}}(\Delta+q_{2}^{*})\\\
\Delta(u_{0}^{*}\cdot)&\Delta(u_{0}\cdot)&\Delta(v_{0}^{*}\cdot)&\Delta(v_{0}\cdot)\end{pmatrix}\,.$

By construction, the constant-coefficient operator
$\mathcal{A}(\mathbf{x}_{0},D)$ with coefficients frozen at $x_{0}$ is
elliptic. Additionally, from the continuity of $u_{0}$ and $v_{0}$ we see that
there exists $\epsilon_{1}=\epsilon_{1}(\eta,\sigma,\mathbf{x}_{0})>0$ such
that $-\frac{u_{0}^{2}}{(u_{0}^{*})^{2}}\neq\frac{v_{0}}{v_{0}^{*}}$ on
$B(\mathbf{x}_{0},\epsilon_{1})$. That is, for every $\Omega\subset
B(\mathbf{x}_{0},\epsilon_{1})$, the operator $\mathcal{A}(\mathbf{x},D)$ is
elliptic on $\Omega$.

Moreover, observe that the Douglis-Nirenberg numbers (34) of $\mathcal{A}$
satisfy $s_{i}=0$ for all $i$ and $t_{j}=2$ is independent of $j$. This means
that the uniqueness theory for elliptic systems presented in [8, Section 3]
applies. More specifically, by Theorem C.3, we conclude that there exists
$\epsilon_{2}=\epsilon_{2}(\eta,\sigma,\mathbf{x}_{0})>0$ such that for every
$\Omega\subset B(\mathbf{x}_{0},\epsilon_{2})$, the boundary value problem
(33) has a unique solution.

Now set $\epsilon=\min\\{\epsilon_{1},\epsilon_{2}\\}>0$ and let
$\Omega\subset B(\mathbf{x}_{0},\epsilon)$, so that $\mathcal{A}$ is an
elliptic operator on $\Omega$ and the problem (33) has a unique solution.

By elliptic regularity estimate Theorem C.2 applied to (33) with $\ell=0$,
there exist constants $C=C(\Omega,\eta,\sigma,u_{0})=C(\Omega,\eta,\sigma)$
and $C_{2}=C_{2}(\Omega,\eta,\sigma)$ such that

$\displaystyle\|\widehat{u}^{(2)}\|_{W^{2,p}(\Omega)}+\|\widehat{v}^{(2)}\|_{W^{2,p}(\Omega)}$
$\displaystyle\leq
C\Bigg{(}\left\|\Delta\left(\frac{\widehat{H}^{(3)}}{3\Gamma\sigma}\right)\right\|_{L^{p}(\Omega)}+\left\|\widehat{J}_{u}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}$
$\displaystyle\qquad\qquad+\left\|\widehat{J}_{v}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}\Bigg{)}+C_{2}(\|\widehat{u}^{(2)}\|_{L^{p}(\Omega)}+\|\widehat{v}^{(2)}\|_{L^{p}(\Omega)}).$

In fact, since the solution is unique we may drop the last term on the right-
hand side, i.e. set $C_{2}=0$. This gives

$\displaystyle\|\widehat{u}^{(2)}\|_{W^{2,p}(\Omega)}+\|\widehat{v}^{(2)}\|_{W^{2,p}(\Omega)}$
$\displaystyle\leq
C\left(\left\|\Delta\left(\frac{\widehat{H}^{(3)}}{3\Gamma\sigma}\right)\right\|_{L^{p}(\Omega)}+\left\|\widehat{J}_{u}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}+\left\|\widehat{J}_{v}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}\right)$
$\displaystyle\leq
C\left(\|\widehat{H}^{(3)}\|_{W^{2,p}(\Omega)}+\left\|\widehat{J}_{u}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}+\left\|\widehat{J}_{v}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}\right),$

where $C=C(\Omega,\Gamma,\eta,\sigma)$.

Finally we substitute this estimate into (32) to obtain

$\displaystyle\|\widehat{\gamma}\|_{L^{p}(\Omega)}$
$\displaystyle=\left\|-\frac{(\Delta+q_{1})\widehat{u}^{(2)}}{2k^{2}u_{0}^{*}v_{0}}\right\|_{L^{p}(\Omega)}\leq
C\|\widehat{u}^{(2)}\|_{W^{2,p}(\Omega)}$ $\displaystyle\leq
C\left(\|\widehat{H}^{(3)}\|_{W^{2,p}(\Omega)}+\left\|\widehat{J}_{u}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}+\left\|\widehat{J}_{v}^{(2)}\right\|_{W^{1-1/p,p}(\partial\Omega)}\right),$

where once again $C=C(\Omega,\Gamma,\eta,\sigma)$. This is precisely the
desired stability estimate (31). ∎

The above theory on the reconstruction of the coefficient $\gamma$ requires
both the availability of the additional boundary data (28) and the assumption
that $\Omega$ is sufficiently small. We made these assumptions merely to
simplify the proof. We believe that they can be removed without breaking the
uniqueness and stability results.

##  5 Numerical experiments

We now present some numerical simulations to demonstrate the quality of
reconstructions that can be achieved for the inverse problem.

We will perform numerical reconstructions with a slightly simplified version
of model (1):

$\begin{array}[]{rcll}\Delta u+k^{2}(1+\eta)u+ik\sigma u&=&0,&\mbox{in}\ \
\Omega\\\ \Delta v+(2k)^{2}(1+\eta)v+i2k\sigma v&=&-(2k)^{2}\gamma
u^{2},&\mbox{in}\ \ \Omega\\\ u=g,\qquad v+i2k{\boldsymbol{\nu}}\cdot\nabla
v&=&0,&\mbox{on}\ \partial\Omega\end{array}$ (35)

In other words, we omit the backward coupling term on the right-hand side of
the first equation in (1). This model is connected to the linearized problem
in (11). Indeed, if we take the boundary condition of $v^{(1)}$ to be $0$,
that is, $h_{1}=0$ in (9), then the first equation in (9) and the second
equation in (11) can be combined to get (35). Note that we intentionally
changed the Dirichlet boundary condition for $v$ to the more realistic Robin
boundary condition. Moreover, due to the fact that the equations in the model
(35) are only one-way coupled, we are not limited to the usage of small
boundary data $g$.

The measured interior data still take the form (2). We will use data generated
from $N_{s}\geq 1$ different boundary conditions $\\{g_{j}\\}_{j=1}^{N_{s}}$:
$\\{H_{j}\\}_{j=1}^{N_{s}}$.

The numerical reconstructions are performed using standard least-squares
optimization procedures that we will outline below. The computational
implementation of the numerical simulations in this section can be found at
https://github.com/nsoedjak/Imaging-SHG. All of the following numerical
experiments can be reproduced by simply running the appropriate example file
(e.g., Experiment_I_gamma.m for Numerical Experiment I).

#### Numerical Experiment I: reconstructing $\gamma$.

We start with reconstructing the coefficient $\gamma$, assuming all other
coefficients are known. The reconstruction is achieved with an optimization
algorithm that finds $\gamma$ by minimizing the functional

$\Phi(\gamma):=\dfrac{1}{2}\sum_{j=1}^{N_{s}}\|\Gamma\sigma(|u_{j}|^{2}+|v_{j}|^{2})-H_{j}\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\beta\|\nabla\gamma\|_{L^{2}(\Omega)}^{2}\,,$

where we assume that we have collected data from $N_{s}$ different boundary
conditions $\\{g_{j}\\}_{j=1}^{N_{s}}$. The regularization parameter $\beta$
will be selected with a trial and error approach. Following the standard
adjoint-state method, we introduce the adjoint equations

$\begin{array}[]{rcll}\Delta w_{j}+(2k)^{2}(1+\eta)w_{j}+i2k\sigma
w_{j}&=&-\Big{[}\Gamma\sigma(|u_{j}|^{2}+|v_{j}|^{2})-H_{j}\Big{]}\Gamma\sigma
v_{j}^{*},&\mbox{in}\ \ \Omega\\\ w_{j}+i2k{\boldsymbol{\nu}}\cdot\nabla
w_{j}&=&0,&\mbox{on}\ \partial\Omega\end{array}$

It is then straightforward to verify that the Fréchet derivative of $\Phi$ in
direction $\delta\gamma$ can be written as:

$\displaystyle\Phi^{\prime}(\gamma)[\delta\gamma]=\int_{\Omega}\left(2(2k)^{2}\Re\sum_{j=1}^{N_{s}}w_{j}u_{j}^{2}\right)\delta\gamma\,d\mathbf{x}-\beta\int_{\Omega}(\Delta\gamma)\delta\gamma\,d\mathbf{x}+\beta\int_{\partial\Omega}\dfrac{\partial{\gamma}}{\partial{\nu}}\delta\gamma\,dS$

Once we have the gradient of the objective function with respect to $\gamma$,
we feed it into a quasi-Newton optimization algorithm with the BFGS updating
rule on the Hessian, implemented in MATLAB.

Figure 1: True $\gamma$ (left), $\gamma$ reconstructed from noise-free data
(middle), and $\gamma$ reconstructed from data containing 1% random noise
(right).

Figure 1 shows the reconstruction of a simple profile of $\gamma$ from both
noise-free and noisy data. The regularization parameter is set to be
$\beta=10^{-7}$ for this particular case. The quality of the reconstructions
is reasonable by visual inspection. Similar levels of reconstruction quality
are observed for various $\gamma$ profiles we tested. The regularization
parameter is selected in a trial-and-error manner. The value of $\beta$ we
used in the simulations may not be the ones to give the best reconstructions.
However, we are not interested in tuning the regularization parameter to
improve the reconstruction quality slightly. Therefore, we will not discuss
this issue here.

#### Numerical Experiment II: reconstructing $(\eta,\sigma,\gamma)$.

In the second numerical example, we consider the case where $\Gamma$ is known
but $\eta$, $\sigma$, and $\gamma$ are unknown. The inversions are done with a
least-squares minimization algorithm that is similar to the one used in
Numerical Experiment I. Figure 2 shows that we are still able to obtain good
numerical reconstructions, at least in the case when the profiles of $\eta$
and $\gamma$ are simple.

Figure 2: True (top) and reconstructed (bottom) $\eta$ (left), $\sigma$
(middle), and $\gamma$ (right) in Numerical Experiment II.

#### Numerical Experiment III: reconstructing $(\eta,\gamma,\Gamma)$.

In this example, we assume that $\sigma$ is known and we are interested in
reconstructing $\eta$, $\gamma$ and $\Gamma$. Due to the fact that $\Gamma$
only appears in the measurement, not the PDE model, a naive least-squares
minimization formulation like the ones in the previous examples will lead to
unbalanced sensitivity between $\Gamma$ and the rest of the parameters. Hence
we instead take a two-step reconstruction approach. In the first step, we use
the ratio between measurements to eliminate $\Gamma$. That is, we minimize the
functional

$\Psi(\eta,\gamma):=\dfrac{1}{2}\sum_{j=2}^{N_{s}}\left\|\frac{|u_{j}|^{2}+|v_{j}|^{2}}{|u_{1}|^{2}+|v_{1}|^{2}}-\frac{H_{j}}{H_{1}}\right\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\beta_{1}\|\nabla\eta\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\beta_{2}\|\nabla\gamma\|_{L^{2}(\Omega)}^{2}\,,$

where we assume that we have collected data from $N_{s}$ different boundary
conditions $\\{g_{j}\\}_{j=1}^{N_{s}}$. It is clear that $\Psi$ only depends
on $\eta$ and $\gamma$, not $\Gamma$. The Fréchet derivatives of $\Psi$ can
again be found using the standard adjoint-state method. For example, for the
derivative with respect to $\gamma$, we introduce the adjoint equations

$\begin{array}[]{rcll}\Delta w_{j}+(2k)^{2}(1+\eta)w_{j}+i2k\sigma
w_{j}&=&-\left(\dfrac{|u_{j}|^{2}+|v_{j}|^{2}}{|u_{1}|^{2}+|v_{1}|^{2}}-\dfrac{H_{j}}{H_{1}}\right)\dfrac{1}{|u_{1}|^{2}+|v_{1}|^{2}}v_{j}^{*},&\mbox{in}\
\ \Omega\\\ w_{j}+i2k{\boldsymbol{\nu}}\cdot\nabla w_{j}&=&0,&\mbox{on}\
\partial\Omega\end{array}$

and

$\begin{array}[]{rcll}\Delta z_{j}+(2k)^{2}(1+\eta)z_{j}+i2k\sigma
z_{j}&=&\left(\dfrac{|u_{j}|^{2}+|v_{j}|^{2}}{|u_{1}|^{2}+|v_{1}|^{2}}-\dfrac{H_{j}}{H_{1}}\right)\dfrac{|u_{j}|^{2}+|v_{j}|^{2}}{(|u_{1}|^{2}+|v_{1}|^{2})^{2}}v_{1}^{*},&\mbox{in}\
\ \Omega\\\ z_{j}+i2k{\boldsymbol{\nu}}\cdot\nabla z_{j}&=&0,&\mbox{on}\
\partial\Omega.\end{array}$

It is then straightforward to verify that the Fréchet derivative of $\Psi$
with respect to $\gamma$ in direction $\delta\gamma$ can be written as:

$\Psi^{\prime}_{\gamma}(\eta,\gamma)[\delta\gamma]=\int_{\Omega}\left(2(2k)^{2}\Re\sum_{j=1}^{N_{s}}\left[w_{j}u_{j}^{2}+z_{j}u_{1}^{2}\right]\right)\delta\gamma\,d\mathbf{x}-\beta_{2}\int_{\Omega}(\Delta\gamma)\delta\gamma\,d\mathbf{x}+\beta_{2}\int_{\partial\Omega}\dfrac{\partial{\gamma}}{\partial{\nu}}\delta\gamma\,dS\,.$

The Fréchet derivative with respect to $\eta$ can be computed in a similar
fashion. Once $\eta$ and $\gamma$ are reconstructed, we can reconstruct
$\Gamma$ as

$\Gamma=\frac{1}{N_{s}}\sum_{j=1}^{N_{s}}\frac{H_{j}}{\sigma(|u_{j}|^{2}+|v_{j}|^{2})}\,.$
(36)

Figure 3: True (top) and reconstructed (bottom) $\eta$ (left), $\gamma$
(middle), and $\Gamma$ (right).

A typical reconstruction is shown in Figure 3. The reconstructions are highly
accurate in this case.

#### Numerical Experiment IV: reconstructing $(\eta,\sigma,\gamma,\Gamma)$.

Figure 4 shows a typical reconstruction of all four coefficients
simultaneously. The reconstruction quality is high in the eyeball norm and can
be characterized more precisely with numbers such as the relative $L^{2}$
error. Note from the reconstruction formula (36) that any inaccuracies in the
reconstruction of $\sigma$ will directly translate into artifacts in the
reconstruction of $\Gamma$. This can be observed in Figure 4 (see columns 2
and 4), most notably near the edges of the square anomaly in $\sigma$.

Figure 4: True (top) and reconstructed (bottom) $(\eta,\sigma,\gamma,\Gamma)$.
From left to right are $\eta$, $\sigma$, $\gamma$, and $\Gamma$.

##  6 Concluding remarks

We performed a systematic study on inverse problems to a system of coupled
semilinear Helmholtz equations as the model for second harmonic generation in
thermoacoustic imaging. We developed uniqueness and stability theory for the
inverse problems utilizing the multilinearization technique. We showed, via
both mathematical analysis and numerical simulations, that it is possible to
reconstruct all four coefficients of interest from noisy interior data.

While our results show great promise for the solution of the inverse problems,
several aspects of our study’s technical side still need to be significantly
improved. For instance, we have assumed the Dirichlet boundary condition for
the generated second harmonic wave $v$ in model (1). This should certainly be
replaced with homogeneous Robin-type boundary conditions that are more
physical (as what we did in the computational experiments). Moreover, in
Theorem 4.1, we should be able to relax the requirement that the domain
$\Omega$ is sufficiently small. In the same theorem, we should be able to
remove the requirement on the additional Neumann boundary data to have a
unique reconstruction of $\gamma$.

We have a few future directions in mind to continue the investigation from the
perspective of practical applications. First, our mathematical results are
mainly based on the assumption that the incident wave, that is, the Dirichlet
boundary condition in system (1), is weak since this is the case where we can
establish the well-posedness of the mathematical model. This assumption,
however, severely limits the applicability of the analysis for practical
applications as one needs to have a sufficiently strong boundary source to
generate strong second-harmonic waves in order to see its impact on the data
used for inversion. Second, the linearization method requires access to a
sequence of datasets generated from $\varepsilon$-dependent boundary source.
This is a large amount of data. It would be interesting to see if our
uniqueness and stability results can be reproduced for a finite number of
measurements. Third, it would be of great interest to see if one can perform a
similar analysis on the same inverse problem to the Maxwell model of second-
harmonic generation, such as the model introduced in [7]. In fact, the
linearization machinery for the Maxwell model has already been built in [7].
However, it is not obvious whether or not our results can be generalized to
the Maxwell model with the same type of data in a straightforward way.

## Acknowledgments

This work is partially supported by the National Science Foundation through
grants DMS-1913309 and DMS-1937254.

##  Appendix A Well-posedness of system (1)

In this appendix, we establish the well-posedness of the boundary value
problem (1) for sufficiently small boundary illuminations $g$ and $h$ using a
standard contraction mapping theorem argument.

We begin by recording a result on the well-posedness of the Helmholtz problem
(15).

######  Theorem A.1.

Let $\alpha\in(0,1)$,
$q,f\in\mathcal{C}^{0,\alpha}(\overline{\Omega};\mathbb{C})$, and
$g\in\mathcal{C}^{2,\alpha}(\partial\Omega;\mathbb{C})$. If $\Im q>0$, then
the boundary value problem

$\begin{array}[]{rcll}\Delta u+qu&=&f,&\mbox{in}\ \Omega\\\ u&=&g,&\mbox{on}\
\partial\Omega\end{array}$ (37)

has a unique solution
$u\in\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$. Moreover, there
exists a constant $C=C(\alpha,\Omega,q)$ such that the following Schauder
estimate holds:

$\displaystyle\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C(\|f\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\,.$
(38)

###### Proof.

The first task is to employ an energy method to show that (37) has at most one
solution $u\in\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$. To this
end, suppose that $u\in\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$
solves the homogeneous problem

$\begin{array}[]{rcll}\Delta u+qu&=&0,&\mbox{in}\ \Omega\\\ u&=&0,&\mbox{on}\
\partial\Omega\,.\end{array}$ (39)

Multiplying both sides of the PDE by $u^{*}$ and integrating over $\Omega$
results in

$\displaystyle\int_{\Omega}-|\nabla u|^{2}+q|u|^{2}\,d\mathbf{x}=0\,,$

whereupon taking imaginary parts yields $\int_{\Omega}\Im
q|u|^{2}\,d\mathbf{x}=0$. The assumption $\Im q>0$ then leads to $u\equiv 0$,
as desired. This completes the proof of uniqueness.

Now that we have shown uniqueness, the existence of a solution
$u\in\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$ to the elliptic
problem (37) follows from the Fredholm alternative: see [1, Theorem 12.7]
(which applies to elliptic operators with not only real-valued but also
complex-valued coefficients).

Finally, from [1, Theorem 7.3] we have the Schauder estimate

$\displaystyle\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C(\|f\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}+\|u\|_{\mathcal{C}^{0}(\overline{\Omega})})\,.$

On account of the problem having a _unique_ solution in
$\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$, the last term
$\|u\|_{\mathcal{C}^{0}(\overline{\Omega})}$ may be dropped (see Remark 2
following [1, Theorem 7.3]) to arrive at the desired (38). ∎

###### Proof of Theorem 2.1.

The proof is a standard argument based on the Banach fixed point theorem.

Before proceeding, we establish some notation. Let us define
$q_{1}:=k^{2}(1+\eta)+ik\sigma$ and $q_{2}:=(2k)^{2}(1+\eta)+i2k\sigma$ for
the sake of brevity of notation. If $X$ and $Y$ are two metric spaces, we
shall equip the Cartesian product $X\times Y$ with any of the standard
metrics, say the metric $d_{X\times
Y}((x_{1},y_{1}),(x_{2},y_{2})):=d_{X}(x_{1},x_{2})+d_{Y}(y_{1},y_{2})$. If
$X$ and $Y$ are both complete, then so is $X\times Y$. Finally, for $\delta>0$
define the complete metric space

$\displaystyle\mathcal{X}_{\delta}:=\\{f\in\mathcal{C}^{2,\alpha}(\overline{\Omega}):\|f\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq\delta\\}\,.$

To start, let $\varepsilon>0$, $\delta>0$ and let
$g,h\in\mathcal{C}^{2,\alpha}(\partial\Omega;\mathbb{C})$ with
$\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}<\varepsilon$ and
$\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}<\varepsilon$. We shall
determine how small $\varepsilon$ and $\delta$ need to be later.

In order to formulate the problem in terms of fixed points, define the
operator
$L^{-1}:\mathcal{C}^{0,\alpha}(\overline{\Omega})\times\mathcal{C}^{0,\alpha}(\overline{\Omega})\rightarrow\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})$
by setting $L^{-1}(f_{1},f_{2})$ as the unique solution
$(U,V)\in\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})$
to the problem

$\begin{array}[]{rcll}\Delta U+q_{1}U&=&f_{1},&\mbox{in}\ \ \Omega\\\ \Delta
V+q_{2}V&=&f_{2},&\mbox{in}\ \ \Omega\\\ U=g,&&V=h,&\mbox{on}\ \
\partial\Omega\,.\end{array}$

Then a solution $(u,v)\in\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$ to
(1) is precisely the same as a fixed point in
$\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$ of the operator $T$ defined
by

$\displaystyle T(\phi_{1},\phi_{2}):=L^{-1}(-k^{2}\gamma\phi_{1}^{*}\phi_{2},\
-(2k)^{2}\gamma\phi_{1}^{2})\,.$

It remains to show that for sufficiently small $\epsilon>0$ and $\delta>0$,

  1. (i)

$T$ is a well-defined operator from
$\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$ to itself, and

  2. (ii)

$T$ is a contraction on $\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$.

In order to perform the next several calculations, recall that
$\mathcal{C}^{0,\alpha}(\overline{\Omega})$ is a Banach algebra, meaning that

$\displaystyle\|fg\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq\|f\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|g\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$

for all $f,g\in\mathcal{C}^{0,\alpha}(\overline{\Omega})$.

#### Proof of (i).

For all
$(\phi_{1},\phi_{2})\in\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$, we
compute

$\displaystyle\|-k^{2}\gamma\phi_{1}^{*}\phi_{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq
C\|\phi_{1}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|\phi_{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq
C\|\phi_{1}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\|\phi_{2}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C\delta^{2}$

and similarly

$\displaystyle\|-(2k)^{2}\gamma\phi_{1}^{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq
C\|\phi_{1}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}^{2}\leq
C\|\phi_{1}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}^{2}\leq C\delta^{2}$

Let $(U,V):=T(\phi_{1},\phi_{2})$. Combining the above estimates with the
Schauder estimate (38) for the Helmholtz equation then leads to

$\displaystyle\|U\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C(\|-k^{2}\gamma\phi_{1}^{*}\phi_{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\leq
C(\delta^{2}+\varepsilon)$

and

$\displaystyle\|V\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C(\|-(2k)^{2}\gamma\phi_{1}^{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\leq
C(\delta^{2}+\varepsilon).$

We can force these quantities to be less than $\delta$ by choosing
$\varepsilon$ and $\delta$ sufficiently small. This implies that
$T(\phi_{1},\phi_{2})\in\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$ as
desired.

#### Proof of (ii).

Let
$(\phi_{1},\phi_{2}),(\phi_{1}^{\prime},\phi_{2}^{\prime})\in\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$.
Then we compute

$\displaystyle\|-k^{2}\gamma\phi_{1}^{*}\phi_{2}-(-k^{2}\gamma(\phi_{1}^{\prime})^{*}\phi_{2}^{\prime})\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C\|\phi_{1}^{*}\phi_{2}-(\phi_{1}^{\prime})^{*}\phi_{2}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$
$\displaystyle=C\|\phi_{1}^{*}(\phi_{2}-\phi_{2}^{\prime})+(\phi_{1}^{*}-(\phi_{1}^{\prime})^{*})\phi_{2}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C(\|\phi_{1}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|\phi_{2}-\phi_{2}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|\phi_{1}-\phi_{1}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|\phi_{2}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})})$
$\displaystyle\leq
C(\|\phi_{1}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\|\phi_{2}-\phi_{2}^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}+\|\phi_{1}-\phi_{1}^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\|\phi_{2}^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})})$
$\displaystyle\leq
C\delta\|(\phi_{1},\phi_{2})-(\phi_{1}^{\prime},\phi_{2}^{\prime})\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})}$

and similarly

$\displaystyle\|-(2k)^{2}\gamma\phi_{1}^{2}-(-(2k)^{2}\gamma(\phi_{1}^{\prime})^{2})\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C\|\phi_{1}^{2}-(\phi_{1}^{\prime})^{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C\|\phi_{1}-\phi_{1}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|\phi_{1}+\phi_{1}^{\prime}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C\delta\|\phi_{1}-\phi_{1}^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C\delta\|(\phi_{1},\phi_{2})-(\phi_{1}^{\prime},\phi_{2}^{\prime})\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})}.$

Let $(U,V):=T(\phi_{1},\phi_{2})$ and
$(U^{\prime},V^{\prime}):=T(\phi_{1}^{\prime},\phi_{2}^{\prime})$. Then
$U-U^{\prime}$ and $V-V^{\prime}$ satisfy

$\begin{array}[]{rcll}\Delta(U-U^{\prime})+q_{1}(U-U^{\prime})&=&-k^{2}\gamma\phi_{1}^{*}\phi_{2}-(-k^{2}\gamma(\phi_{1}^{\prime})^{*}\phi_{2}^{\prime}),&\mbox{in}\
\ \Omega\\\
\Delta(V-V^{\prime})+q_{2}(V-V^{\prime})&=&-(2k)^{2}\gamma\phi_{1}^{2}-(-(2k)^{2}\gamma(\phi_{1}^{\prime})^{2}),&\mbox{in}\
\ \Omega\\\ U-U^{\prime}=0,\qquad V-V^{\prime}&=&0,&\mbox{on}\ \
\partial\Omega\,.\end{array}$

Combining the above estimates with the Schauder estimate (38) for the
Helmholtz equation then leads to

$\displaystyle\|U-U^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C\|-k^{2}\gamma\phi_{1}^{*}\phi_{2}-(-k^{2}\gamma(\phi_{1}^{\prime})^{*}\phi_{2}^{\prime})\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq
C\delta\|(\phi_{1},\phi_{2})-(\phi_{1}^{\prime},\phi_{2}^{\prime})\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})}$

and

$\displaystyle\|V-V^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C\|-(2k)^{2}\gamma\phi_{1}^{2}-(-(2k)^{2}\gamma(\phi_{1}^{\prime})^{2})\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq
C\delta\|(\phi_{1},\phi_{2})-(\phi_{1}^{\prime},\phi_{2}^{\prime})\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})}.$

We conclude that

$\displaystyle\|T(\phi_{1},\phi_{2})-T(\phi_{1}^{\prime},\phi_{2}^{\prime})\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle=\|U-U^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}+\|V-V^{\prime}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle\leq
C\delta\|(\phi_{1},\phi_{2})-(\phi_{1}^{\prime},\phi_{2}^{\prime})\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})\times\mathcal{C}^{2,\alpha}(\overline{\Omega})}.$

The factor $C\delta$ can be made strictly less than $1$ when $\delta$ is
sufficiently small. This makes $T$ into a contraction, as desired.

Having proved that $T$ is a contraction on the complete metric space
$\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$, the Banach fixed point
theorem guarantees that there exists a unique
$(u,v)\in\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$ such that
$T(u,v)=(u,v)$. As discussed earlier, this is equivalent to saying that there
exists a unique $(u,v)\in\mathcal{X}_{\delta}\times\mathcal{X}_{\delta}$
satisfying the boundary value problem (1). This completes the proof of the
first part of the theorem.

#### Proof of the estimates (4).

We perform a calculation similar to those in the proof of (i) to obtain

$\displaystyle\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle\leq C(\|-k^{2}\gamma
u^{*}v\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\leq
C(\|u\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|v\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})$
$\displaystyle\leq
C(\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\delta+\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})$

When $\delta$ is sufficiently small, this implies that

$\displaystyle\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}\,,$

as desired. To get the estimate for $v$, we calculate

$\displaystyle\|v\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}$
$\displaystyle\leq C(\|-(2k)^{2}\gamma
u^{2}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\leq
C(\|u\|_{\mathcal{C}^{0,\alpha}(\partial\Omega)}^{2}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})$
$\displaystyle\leq
C(\delta\|u\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\leq
C(\delta\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})$
$\displaystyle\leq
C(\|g\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}+\|h\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)})\,,$

as desired. The proof is complete. ∎

##  Appendix B Differentiability result for linearization

We provide here the mathematical justification of the linearization process we
outlined in Section 2.2. More precisely, we prove Theorem 2.2 by showing that
$(u_{\varepsilon},v_{\varepsilon})$ and therefore $H_{\varepsilon}$ are
differentiable with respect to $\varepsilon$.

###### Proof of Theorem 2.2.

Let us define $q_{1}:=k^{2}(1+\eta)+ik\sigma$ and
$q_{2}:=(2k)^{2}(1+\eta)+i2k\sigma$ to ease notation.

To start the proof, let $\widetilde{u}^{(1)}$, $\widetilde{v}^{(1)}$,
$\widetilde{u}^{(2)}$, $\widetilde{v}^{(2)}$ denote the unique functions in
$\mathcal{C}^{2,\alpha}(\overline{\Omega};\mathbb{C})$ which satisfy equations
(9) and (11). That is,

$\begin{array}[]{rcll}\Delta\widetilde{u}^{(1)}+q_{1}\widetilde{u}^{(1)}&=&0,&\mbox{in}\
\ \Omega\\\ \Delta\widetilde{v}^{(1)}+q_{2}\widetilde{v}^{(1)}&=&0,&\mbox{in}\
\ \Omega\\\ \widetilde{u}^{(1)}=g_{1},&&\widetilde{v}^{(1)}=h_{1},&\mbox{on}\
\ \partial\Omega\end{array}$

and

$\begin{array}[]{rcll}\Delta\widetilde{u}^{(2)}+q_{1}\widetilde{u}^{(2)}&=&-2k^{2}\gamma\widetilde{u}^{(1)*}\widetilde{v}^{(1)},&\mbox{in}\
\ \Omega\\\
\Delta\widetilde{v}^{(2)}+q_{2}\widetilde{v}^{(2)}&=&-2(2k)^{2}\gamma(\widetilde{u}^{(1)})^{2},&\mbox{in}\
\ \Omega\\\ \widetilde{u}^{(2)}=g_{2},&&\widetilde{v}^{(2)}=h_{2},&\mbox{on}\
\ \partial\Omega\,.\end{array}$

(Note that the existence and uniqueness of these functions is guaranteed by
Theorem A.1.) Define now the “remainder” terms

$\displaystyle\mu_{\varepsilon}$
$\displaystyle:=u_{\varepsilon}-\varepsilon\widetilde{u}^{(1)}-\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)},$
(40) $\displaystyle\nu_{\varepsilon}$
$\displaystyle:=v_{\varepsilon}-\varepsilon\widetilde{v}^{(1)}-\frac{1}{2}\varepsilon^{2}\widetilde{v}^{(2)}.$

We wish to show that $\mu_{\varepsilon}$ and $\nu_{\varepsilon}$ are in a
certain sense “$o(\varepsilon^{2})$” as $\varepsilon\rightarrow 0$. This will
be accomplished in two rounds of estimates on $\mu_{\varepsilon}$,
$\nu_{\varepsilon}$, $u_{\varepsilon}$ and $v_{\varepsilon}$.

#### Round 1 estimates.

We begin by using the linearity of the operators $\Delta+q_{1}$ and
$\Delta+q_{2}$ to find that

$\begin{array}[]{rcll}\Delta\mu_{\varepsilon}+q_{1}\mu_{\varepsilon}&=&-k^{2}\gamma[u_{\varepsilon}^{*}v_{\varepsilon}-\varepsilon^{2}\widetilde{u}^{(1)*}\widetilde{v}^{(1)}],&\mbox{in}\
\ \Omega\\\
\Delta\nu_{\varepsilon}+q_{2}\nu_{\varepsilon}&=&-(2k)^{2}\gamma[u_{\varepsilon}^{2}-\varepsilon^{2}(\widetilde{u}^{(1)})^{2}],&\mbox{in}\
\ \Omega\\\ \mu_{\varepsilon}=0,&&\nu_{\varepsilon}=0,&\mbox{on}\ \
\partial\Omega\,.\end{array}$ (41)

To obtain control on the size of the right hand sides, we utilize the well-
posedness result Theorem 2.1 to see that

$\displaystyle\|u_{\varepsilon}\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq\|u_{\varepsilon}\|_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}\leq
C\left(\left\|\varepsilon
g_{1}+\frac{1}{2}\varepsilon^{2}g_{2}\right\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}+\left\|\varepsilon
h_{1}+\frac{1}{2}\varepsilon^{2}h_{2}\right\|_{\mathcal{C}^{2,\alpha}(\partial\Omega)}\right)\leq
C\varepsilon,$

and similarly for $v_{\varepsilon}$. Here,
$C=C(\alpha,\Omega,\eta,\sigma,\gamma,g_{1},g_{2})$ is a constant not
depending on $\varepsilon$. We can write these bounds succinctly as

$\displaystyle
u_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon),\qquad
v_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon).$
(42)

In order to perform the next several calculations, recall that
$\mathcal{C}^{0,\alpha}(\overline{\Omega})$ is a Banach algebra, meaning that

$\displaystyle\|fg\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\leq\|f\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}\|g\|_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}$

for all $f,g\in\mathcal{C}^{0,\alpha}(\overline{\Omega})$. With the help of
this property, we plug the bounds (42) into the right hand sides of (41) to
discover that

$\displaystyle-k^{2}\gamma[u_{\varepsilon}^{*}v_{\varepsilon}-\varepsilon^{2}\widetilde{u}^{(1)*}\widetilde{v}^{(1)}]=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2}),\qquad-(2k)^{2}\gamma[u_{\varepsilon}^{2}-\varepsilon^{2}(\widetilde{u}^{(1)})^{2}]=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2}).$
(43)

The Schauder estimate (38) for the Helmholtz equation applied to (41) then
gives

$\displaystyle\mu_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}(\varepsilon^{2}),\qquad\nu_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{2,\alpha}(\overline{\Omega})}(\varepsilon^{2}),$
(44)

and in particular

$\displaystyle\mu_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2}),\qquad\nu_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2}).$
(45)

#### Round 2 estimates.

Using the estimates from Round 1 and recalling the definition (40) of the
remainder terms $\mu_{\varepsilon}$ and $\nu_{\varepsilon}$, we can now refine
the bounds (43) on the right hand sides of (41):

$\displaystyle-k^{2}\gamma[u_{\varepsilon}^{*}v_{\varepsilon}-\varepsilon^{2}\widetilde{u}^{(1)*}\widetilde{v}^{(1)}]$
$\displaystyle=-k^{2}\gamma\left[\left(\varepsilon\widetilde{u}^{(1)*}+\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)*}+\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2})\right)\left(\varepsilon\widetilde{v}^{(1)}+\frac{1}{2}\varepsilon^{2}\widetilde{v}^{(2)}+\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2})\right)-\varepsilon^{2}\widetilde{u}^{(1)*}\widetilde{v}^{(1)}\right]$
$\displaystyle=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{3})$

and

$\displaystyle-(2k)^{2}\gamma[u_{\varepsilon}^{2}-\varepsilon^{2}(\widetilde{u}^{(1)})^{2}]$
$\displaystyle=-(2k)^{2}\gamma\left[\left(\varepsilon\widetilde{u}^{(1)}+\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)}+\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{2})\right)^{2}-\varepsilon^{2}(\widetilde{u}^{(1)})^{2}\right]=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{3}).$

With these improved bounds in hand, we again apply the Schauder estimate (38)
to (41) to obtain the following bounds for the remainder terms
$\mu_{\varepsilon}$ and $\nu_{\varepsilon}$:

$\displaystyle\mu_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{3}),\qquad\nu_{\varepsilon}=\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{3}).$

Note that this is a refinement over the previous remainder bounds (45). This
concludes the Round 2 estimates.

Therefore from the definition (40) of $\mu_{\varepsilon}$ and
$\nu_{\varepsilon}$ we conclude that

$\displaystyle u_{\varepsilon}$
$\displaystyle=\varepsilon\widetilde{u}^{(1)}+\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)}+\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{3}),$
$\displaystyle v_{\varepsilon}$
$\displaystyle=\varepsilon\widetilde{v}^{(1)}+\frac{1}{2}\varepsilon^{2}\widetilde{v}^{(2)}+\mathcal{O}_{\mathcal{C}^{0,\alpha}(\overline{\Omega})}(\varepsilon^{3}),$

and in particular for each $\mathbf{x}\in\overline{\Omega}$ we have the
pointwise estimates

$\displaystyle u_{\varepsilon}(\mathbf{x})$
$\displaystyle=\varepsilon\widetilde{u}^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)}(\mathbf{x})+\mathcal{O}(\varepsilon^{3}),$
$\displaystyle v_{\varepsilon}(\mathbf{x})$
$\displaystyle=\varepsilon\widetilde{v}^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}\widetilde{v}^{(2)}(\mathbf{x})+\mathcal{O}(\varepsilon^{3}),$

as $\varepsilon\rightarrow 0$.

#### Asymptotic expansion of $H_{\varepsilon}$.

To finish the proof, for each $\mathbf{x}\in\overline{\Omega}$ we compute:

$\displaystyle H_{\varepsilon}(\mathbf{x})$
$\displaystyle=\Gamma(\mathbf{x})\sigma(\mathbf{x})[u_{\varepsilon}(\mathbf{x})u^{*}_{\varepsilon}(\mathbf{x})+v_{\varepsilon}(\mathbf{x})v^{*}_{\varepsilon}(\mathbf{x})]$
$\displaystyle=\Gamma(\mathbf{x})\sigma(\mathbf{x})\Bigg{[}\left(\varepsilon\widetilde{u}^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)}(\mathbf{x})+\mathcal{O}(\varepsilon^{3})\right)\left(\varepsilon\widetilde{u}^{(1)*}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}\widetilde{u}^{(2)*}(\mathbf{x})+\mathcal{O}(\varepsilon^{3})\right)$
$\displaystyle\qquad+\left(\varepsilon\widetilde{v}^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}\widetilde{v}^{(2)}(\mathbf{x})+\mathcal{O}(\varepsilon^{3})\right)\left(\varepsilon\widetilde{v}^{(1)*}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}\widetilde{v}^{(2)*}(\mathbf{x})+\mathcal{O}(\varepsilon^{3})\right)\Bigg{]}$
$\displaystyle=\varepsilon
H^{(1)}(\mathbf{x})+\frac{1}{2}\varepsilon^{2}H^{(2)}(\mathbf{x})+\frac{1}{6}\varepsilon^{3}H^{(3)}(\mathbf{x})+\mathcal{O}(\varepsilon^{4}),$

where

$\displaystyle H^{(1)}$ $\displaystyle:=0,$ $\displaystyle H^{(2)}$
$\displaystyle:=2\Gamma\sigma(\widetilde{u}^{(1)*}\widetilde{u}^{(1)}+\widetilde{v}^{(1)*}\widetilde{v}^{(1)}),$
$\displaystyle H^{(3)}$
$\displaystyle:=3\Gamma\sigma\left(\widetilde{u}^{(1)*}\widetilde{u}^{(2)}+\widetilde{u}^{(1)}\widetilde{u}^{(2)*}+\widetilde{v}^{(1)*}\widetilde{v}^{(2)}+\widetilde{v}^{(1)}\widetilde{v}^{(2)*}\right).$

These are precisely the equations (10), (12), and (13), respectively. This
completes the proof. ∎

##  Appendix C Elliptic system theory

For the sake of completeness, in this appendix we review the elliptic system
theory that appears in the proof of Theorem 4.1. We mostly follow the
presentation of [8].

Let $M\geq N$ be positive integers and consider the following system of $M$
partial differential equations in $N$ unknown functions $v_{1},\dots,v_{N}$
defined on an open set $\Omega$:

$\displaystyle\mathcal{A}(\mathbf{x},D)v=\mathcal{S}\quad\mbox{in}\ \Omega.$

Here $v=(v_{1},\dots,v_{N})$, $D=(\partial_{x_{1}},\dots,\partial_{x_{n}})$,
and $\mathcal{A}(\mathbf{x},D)$ is an $M\times N$ matrix linear partial
differential operator. That is, for each $1\leq i\leq M$ and $1\leq j\leq N$,
the entry $\mathcal{A}_{ij}(\mathbf{x},D)$ is a linear partial differential
operator, and the above matrix equation means that

$\displaystyle\sum_{j=1}^{N}\mathcal{A}_{ij}(\mathbf{x},D)v_{j}=\mathcal{S}_{i},\qquad
1\leq i\leq M.$

To each row $1\leq i\leq M$ let us now associate an integer $s_{i}$, and to
each column $1\leq j\leq N$ let us associate an integer $t_{j}$. Let us choose
the numbers in such a way that the order of the partial differential operator
$\mathcal{A}_{ij}(\mathbf{x},D)$ is no greater than $s_{i}+t_{j}$. (If
$s_{i}+t_{j}<0$, then we require that $\mathcal{A}_{ij}(\mathbf{x},D)=0$.)

The _principal part_ $\mathcal{A}_{0}(\mathbf{x},D)$ of
$\mathcal{A}(\mathbf{x},D)$ is defined as the $M\times N$ matrix linear
partial differential operator such that $(\mathcal{A}_{0})_{ij}(\mathbf{x},D)$
consists of the terms in $\mathcal{A}_{ij}(\mathbf{x},D)$ of order exactly
$s_{i}+t_{j}$.

######  Definition C.1.

The matrix partial differential operator $\mathcal{A}$ is called _elliptic_ if
such Douglis-Nirenberg numbers $(s_{i})_{1\leq i\leq M}$ and $(t_{j})_{1\leq
j\leq N}$ exist, and the matrix
$\mathcal{A}_{0}(\mathbf{x},{\boldsymbol{\xi}})$ has full rank $N$ for each
$\mathbf{x}\in\Omega$ and ${\boldsymbol{\xi}}\in\mathbb{S}^{n-1}$. (The matrix
$\mathcal{A}_{0}(\mathbf{x},{\boldsymbol{\xi}})$ is called the _symbol_ of the
operator $\mathcal{A}_{0}$.)

We now summarize the parts of [8, Section 3] that are relevant for the proof
of Theorem 4.1. Assume from now on that $\mathcal{A}$ is an elliptic operator
with continuous coefficients and Douglis-Nirenberg numbers $(s_{i})_{1\leq
i\leq M}$ and $(t_{j})_{1\leq j\leq N}$ such that

$\begin{array}[]{rcll}s_{i}&=&0,&1\leq i\leq M,\\\ t_{j}&=&\tau,&1\leq j\leq
N,\end{array}$ (46)

for some positive integer $\tau$. In the following, we will consider the
Dirichlet boundary value problem

$\begin{array}[]{rcll}\mathcal{A}(\mathbf{x},D)v&=&\mathcal{S}&\mbox{in}\
\Omega,\\\
\left(\dfrac{\partial{}}{\partial{\nu}}\right)^{q}v_{j}&=&\phi_{qj}&\mbox{on}\
\partial\Omega,\quad 0\leq q\leq\tau-1,\quad 1\leq j\leq N.\end{array}$ (47)

For this boundary value problem, we have the following estimate.

######  Theorem C.2.

[8, p. 13] Let $p>1$ and $\ell\geq 0$. Assume that $\mathcal{S}_{i}\in
W^{\ell,p}(\Omega)$ for all $1\leq i\leq M$ and $\phi_{qj}\in
W^{\ell+\tau-q-\frac{1}{p},p}(\partial\Omega)$ for all $0\leq q\leq\tau-1$ and
$1\leq j\leq N$. Also, assume that all the coefficients in $\mathcal{A}$ are
in $\mathcal{C}^{\ell}(\overline{\Omega})$ [34, Theorem 1.1]. Then the
following elliptic regularity estimate holds for the boundary value problem
(47):

$\displaystyle\sum_{j=1}^{N}\|v_{j}\|_{W^{\ell+\tau,p}(\Omega)}$
$\displaystyle\leq
C\left(\sum_{i=1}^{M}\|\mathcal{S}_{i}\|_{W^{\ell,p}(\Omega)}+\sum_{q,j}\|\phi_{qj}\|_{W^{\ell+\tau-q-\frac{1}{p},p}(\partial\Omega)}\right)+C_{2}\sum_{j=1}^{N}\|v_{j}\|_{L^{p}(\Omega)}$
(48)

The following result on uniqueness of solutions in sufficiently small domains
is used in the proof of Theorem 4.1.

######  Theorem C.3.

[8, Theorem 3.5] Let $\mathbf{x}_{0}\in\Omega$. Then there exists $\epsilon>0$
such that for every small domain $\Omega^{\prime}\subset
B(\mathbf{x}_{0},\epsilon)$, the only solution to the homogeneous Dirichlet
boundary value problem

$\begin{array}[]{rcll}\mathcal{A}(\mathbf{x},D)v&=&0&\mbox{in}\
\Omega^{\prime},\\\
\left(\dfrac{\partial{}}{\partial{\nu}}\right)^{q}v_{j}&=&0&\mbox{on}\
\partial\Omega^{\prime},\quad 0\leq q\leq\tau-1,\quad 1\leq j\leq
N.\end{array}$ (49)

in $\Omega^{\prime}$ is the trivial solution $v=0$. In particular, (48) holds
with $C_{2}=0$.

## References

  * [1] S. Agmon, A. Douglis, and L. Nirenberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. I, Comm. Pure Appl. Math., 12 (1959), pp. 623–727.
  * [2] H. Akhouayri, M. Bergounioux, A. Da Silva, P. Elbau, A. Litman, and L. Mindrinos, Quantitative thermoacoustic tomography with microwaves sources, J. Inverse Ill-Posed Probl., 25 (2017), pp. 703–717.
  * [3] H. Al Jebawy and A. El Badia, Direct algorithm for reconstructing small absorbers in thermoacoustic tomography problem from a single data, Inverse Problems, 36 (2020), p. 065010.
  * [4] G. S. Alberti, Non-zero constraints in elliptic PDE with random boundary values and applications to hybrid inverse problems, arXiv:2205.00994, (2022).
  * [5] L. Ambrosio, Transport equation and Cauchy problem for BV vector fields, Inventiones Mathematicae, 158 (2004), p. 227.
  * [6] H. Ammari, J. Garnier, W. Jing, and L. Nguyen, Quantitative thermo-acoustic imaging: An exact reconstruction formula, Submitted to Archive for Rational Mechanics and Analysis, (2011).
  * [7] Y. M. Assylbekov and T. Zhou, Inverse problems for nonlinear maxwell’s equations with second harmonic generation, Journal of Differential Equations, 296 (2021), pp. 148–169.
  * [8] G. Bal, Hybrid inverse problems and redundant systems of partial differential equations, in Inverse Problems and Applications, P. Stefanov, A. Vasy, and M. Zworski, eds., vol. 615 of Contemporary Mathematics, American Mathematical Society, 2013, pp. 15–48.
  * [9] G. Bal and K. Ren, Multi-source quantitative PAT in diffusive regime, Inverse Problems, 27 (2011).  075003\.
  * [10] , On multi-spectral quantitative photoacoustic tomography in diffusive regime, Inverse Problems, 28 (2012).  025010\.
  * [11] G. Bal, K. Ren, G. Uhlmann, and T. Zhou, Quantitative thermo-acoustics and related problems, Inverse Problems, 27 (2011).  055007\.
  * [12] G. Bal and T. Zhou, Hybrid inverse problems for a system of Maxwell’s equations, Inverse Problems, 30 (2014).  055013\.
  * [13] M. Bergounioux, É. Bretin, and Y. Privat, How to position sensors in thermo-acoustic tomography, Inverse Problems, 35 (2019), p. 074003.
  * [14] L. Borcea, W. Li, A. Mamonov, and J. C. Schotland, Second-harmonic imaging in random media, Inverse Problems, 33 (2017).  065004\.
  * [15] M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, Cambridge University Press, New York, 1999.
  * [16] F. Bouchut and G. Crippa, Uniqueness, renormalization and smooth approximations for linear transport equations, SIAM J. Math. Anal., 38 (2006), pp. 1316–1328.
  * [17] R. W. Boyd, Nonlinear optics, Academic press, 2020.
  * [18] M. Choulli, Stable determination of the nonlinear term in a quasilinear elliptic equation by boundary measurements, arXiv:2205.16000, (2022).
  * [19] F. Colombini and N. Lerner, Uniqueness of continuous solutions for BV vector fields, Duke Math. J., 111 (2002), pp. 357–384.
  * [20] M. Cristofol, S. Li, and Y. Shang, Carleman estimates and some inverse problems for the coupled quantitative thermoacoustic equations by partial boundary layer data. part ii: Some inverse problems, Mathematical Methods in the Applied Sciences, (2023).
  * [21] R. J. DiPerna and P.-L. Lions, On the Cauchy problem for Boltzmann equations: global existence and weak stability, Ann. Math., 130 (1989), pp. 321–366.
  * [22] A. Feizmohammadi, T. Liimatainen, and Y.-H. Lin, An inverse problem for a semilinear elliptic equation on conformally transversally anisotropic manifolds, arXiv:2112.08305, (2021).
  * [23] J. Francés, J. Tervo, S. Gallego, S. Bleda, C. Neipp, and A. Márquez, Split-field finite-difference time-domain method for second-harmonic generation in two-dimensionally periodic structures, J. Opt. Soc. Am. B, 32 (2015), pp. 664–669.
  * [24] B. Harrach and Y.-H. Lin, Simultaneous recovery of piecewise analytic coefficients in a semilinear elliptic equation, Nonlinear Analysis, 228 (2023), p. 113188.
  * [25] V. Isakov, On uniqueness in inverse problems for semilinear parabolic equations, Arch. Rational Mech. Anal., 124 (1993), pp. 1–12.
  * [26] , Inverse Problems for Partial Differential Equations, Springer-Verlag, New York, second ed., 2006.
  * [27] Y. Kian, Lipschitz and hölder stable determination of nonlinear terms for elliptic equations, Nonlinearity, 36 (2023), p. 1302.
  * [28] K. Krupchyk and G. Uhlmann, Partial data inverse problems for semilinear elliptic equations with gradient nonlinearities, arXiv:1909.08122v1, (2019).
  * [29] , A remark on partial data inverse problems for semilinear elliptic equations, Proceedings of the AMS, (2019).
  * [30] R.-Y. Lai and Y.-H. Lin, Inverse problems for fractional semilinear elliptic equations, Nonlinear Analysis, 216 (2022), p. 112699.
  * [31] R.-Y. Lai, K. Ren, and T. Zhou, Inverse transport and diffusion problems in photoacoustic imaging with nonlinear absorption, SIAM J. Appl. Math., 82 (2022), pp. 602–624.  arXiv:2107.08118.
  * [32] M. Lassas, T. Liimatainen, Y.-H. Lin, and M. Salo, Inverse problems for elliptic equations with power type nonlinearities, Journal de mathématiques pures et appliquées, 145 (2021), pp. 44–82.
  * [33] S. Lu and J. Zhai, Increasing stability of a linearized inverse boundary value problem for a nonlinear schrödinger equation on transversally anisotropic manifolds, arXiv:2301.07875, (2023).
  * [34] V. A. Solonnikov, Overdetermined elliptic boundary-value problems, J. Sov. Math., 1 (1973), pp. 477–512.
  * [35] T. Szarvas and Z. Kis, Numerical simulation of nonlinear second harmonic wave generation by the finite difference frequency domain method, J. Opt. Soc. Am. B, 35 (2018), pp. 731–740.
  * [36] G. Uhlmann and J. Zhai, On an inverse boundary value problem for a nonlinear elastic wave equation, Journal de Mathématiques Pures et Appliquées, 153 (2021), pp. 114–136.
  * [37] J. Yuan and J. Yang, Computational design for efficient second-harmonic generation in nonlinear photonic crystals, J. Opt. Soc. Am. B, 30 (2013), pp. 205–210.
  * [38] Y. Zeng, W. Hoyer, J. Liu, S. W. Koch, and J. V. Moloney, Classical theory for second-harmonic generation from metallic nanoparticles, Phys. Rev. B, 79 (2009).  235109\.