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

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

$$\begin{equation} \begin{array}{r@{\ }c@{\ }ll} \Delta u+ k^2 \left(1+\eta\right) u +ik\sigma u& =& - 2 k^2 \gamma u^* v, & \mbox{in}\ \ \Omega\\ \Delta v+ \left(2k\right)^2 \left(1+\eta\right) v +i2k\sigma v & = & -\left(2k\right)^2 \gamma u^2, &\mbox{in}\ \ \Omega\\ u = g, & & \quad v = h, & \mbox{on}\ \partial\Omega \end{array} \end{equation} \tag{ 1 }$$

where Δ 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 [14, 15, 17, 23, 35, 37, 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 η and γ, respectively, and an absorption coefficient σ.

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:

$$\begin{align} H\left({\mathbf{x}}\right) = \Gamma\left({\mathbf{x}}\right) \sigma\left({\mathbf{x}}\right) \left(|u|^2+|v|^2\right),\qquad {\mathbf{x}}\in\Omega\,. \end{align} \tag{ 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 [3, 12, 13, 20] 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 Γ, which is mostly ignored in the previous studies of quantitative thermoacoustic imaging [2, 6, 11, 12].

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 ω) and generated (at frequency 2ω) electric wave of the form:

$$\begin{align*} u\left({\mathbf{x}},t\right) & = 2\Re\left(u\left({\mathbf{x}}\right)e^{-i\omega t}\right) = 2|u\left({\mathbf{x}}\right)|\text{cos}\left(\varphi_u\left({\mathbf{x}}\right) - \omega t\right)\,,\\ v\left({\mathbf{x}},t\right) & = 2\Re\left(v\left({\mathbf{x}}\right)e^{-i2\omega t}\right) = 2|v\left({\mathbf{x}}\right)|\text{cos}\left(\varphi_v\left({\mathbf{x}}\right) - 2\omega t\right)\,, \end{align*}$$

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

$$\begin{align*} I\left({\mathbf{x}}\right) & = \frac{1}{T}\int_0^T \frac{1}{2} |u\left({\mathbf{x}},t\right)+v\left({\mathbf{x}},t\right)|^2 \,\textrm{d}t = \frac{1}{2T}\int_0^T u\left({\mathbf{x}},t\right)^2 + v\left({\mathbf{x}},t\right)^2 + 2u\left({\mathbf{x}},t\right)v\left({\mathbf{x}},t\right)\, \textrm{d}t\\ & = 4|u\left({\mathbf{x}}\right)|^2\frac{1}{2T}\int_0^T\text{cos}^2\left(\varphi_u\left({\mathbf{x}}\right) - \omega t\right)\,\textrm{d}t + 4|v\left({\mathbf{x}}\right)|^2\frac{1}{2T}\int_0^T\text{cos}^2\left(\varphi_v\left({\mathbf{x}}\right) - 2\omega t\right)\,\textrm{d}t\\ & \quad + 8|u\left({\mathbf{x}}\right)||v\left({\mathbf{x}}\right)|\frac{1}{2T}\int_0^T \text{cos}\left(\varphi_u\left({\mathbf{x}}\right) - \omega t\right)\text{cos}\left(\varphi_v\left({\mathbf{x}}\right) - 2\omega t\right)\,\textrm{d}t\\ & = \left(|u\left({\mathbf{x}}\right)|^2+|v\left({\mathbf{x}}\right)|^2\right)\,, \end{align*}$$

where we have used the standard trigonometic identity $\text{cos}(x)\text{cos}(y) = \frac{1}{2}(\text{cos}(x+y)+\text{cos}(x-y))$ to simplify the integrals. Therefore, the cross-term vanishes due to orthogonality. The absorbed radiation at location 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:

$$\begin{align} \Lambda_{\Gamma, \eta, \gamma, \sigma}: \left(g,h\right) \mapsto H. \end{align} \tag{ 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 γ (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 γ 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.

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

• ($\mathcal{A}$-i)  
  The domain Ω is bounded with smooth boundary $\partial\Omega$.
• ($\mathcal{A}$-ii)  
  The coefficients $\Gamma,\eta,\sigma,\gamma$ all lie in the set

  $$\begin{align*} {\mathcal{Y}} : = \left\{f\in {\mathcal{C}}^2\left(\overline{\Omega}; {\mathbb{R}}\right): c_1\unicode{x2A7D} f\unicode{x2A7D} c_2\ \mbox{in}\ \overline{\Omega}\right\} \end{align*}$$

  for some $c_1\gt0$ and $c_2\gt0$.

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 ε > 0 and δ > 0 such that for all $g,h\in {\mathcal{C}}^{2,\alpha}(\partial\Omega;\mathbb C)$ with $\|g\|_{{\mathcal{C}}^{2,\alpha}(\partial\Omega)}\lt\varepsilon$ and $\|h\|_{{\mathcal{C}}^{2,\alpha}(\partial\Omega)}\lt\varepsilon$, the boundary value problem (1) has a unique solution

$$\begin{align*} \left(u,v\right)\in \left\{f\in {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega};\mathbb C\right): \|{\,}f{\,}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\unicode{x2A7D} \delta\right\}^2. \end{align*}$$

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

$$\begin{align} \begin{aligned} \|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} &\unicode{x2A7D} C\left( \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right),\\ \|v\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} &\unicode{x2A7D} C\left( \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right). \end{aligned} \end{align} \tag{ 4 }$$

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, the 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, 24–29, 32, 33]. We now document the linearization process.

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

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

with boundary conditions

$$\begin{align} 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\,. \end{align} \tag{ 6 }$$

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

$$\begin{align*} H_\varepsilon = \Gamma \sigma \left(|u_\varepsilon|^2 + |v_\varepsilon|^2\right)\,.\end{align*}$$

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

$$\begin{align} \begin{aligned} u_\varepsilon\left({\mathbf{x}}\right) & = \varepsilon u^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12 \varepsilon^2 u^{\left(2\right)}\left({\mathbf{x}}\right) + o\left(\varepsilon^2\right),\\ v_\varepsilon\left({\mathbf{x}}\right) & = \varepsilon v^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12 \varepsilon^2 v^{\left(2\right)}\left({\mathbf{x}}\right) + o\left(\varepsilon^2\right),\\ H_\varepsilon\left({\mathbf{x}}\right) & = \varepsilon H^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12 \varepsilon^2 H^{\left(2\right)}\left({\mathbf{x}}\right) + \frac16 \varepsilon^3 H^{\left(3\right)}\left({\mathbf{x}}\right) + o\left(\varepsilon^3\right), \end{aligned} \end{align} \tag{ 7 }$$

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

$$\begin{align} \begin{split} u^{\left(1\right)}\left({\mathbf{x}}\right) &: = \lim_{\varepsilon\rightarrow 0} \frac{u_\varepsilon\left({\mathbf{x}}\right)}{\varepsilon},\qquad u^{\left(2\right)}\left({\mathbf{x}}\right) : = \lim_{\varepsilon\rightarrow 0} \frac{u_\varepsilon\left({\mathbf{x}}\right) - \varepsilon u^{\left(1\right)}\left({\mathbf{x}}\right)}{\frac12 \varepsilon^2},\\ v^{\left(1\right)}\left({\mathbf{x}}\right) &: = \lim_{\varepsilon\rightarrow 0} \frac{v_\varepsilon\left({\mathbf{x}}\right)}{\varepsilon},\qquad v^{\left(2\right)}\left({\mathbf{x}}\right) : = \lim_{\varepsilon\rightarrow 0} \frac{v_\varepsilon\left({\mathbf{x}}\right) - \varepsilon v^{\left(1\right)}\left({\mathbf{x}}\right)}{\frac12 \varepsilon^2},\\ H^{\left(1\right)}\left({\mathbf{x}}\right) &: = \lim_{\varepsilon\rightarrow 0} \frac{H_\varepsilon\left({\mathbf{x}}\right)}{\varepsilon},\qquad H^{\left(2\right)}\left({\mathbf{x}}\right) : = \lim_{\varepsilon\rightarrow 0} \frac{H_\varepsilon\left({\mathbf{x}}\right) - \varepsilon H^{\left(1\right)}\left({\mathbf{x}}\right)}{\frac12 \varepsilon^2},\\ H^{(3)}({\mathbf{x}}) &: = \lim_{\varepsilon\rightarrow 0} \frac{H_\varepsilon({\mathbf{x}}) - \varepsilon H^{(1)}({\mathbf{x}}) - \frac12\varepsilon^2 H^{(2)}({\mathbf{x}})}{\frac16 \varepsilon^3}. \end{split}\end{align} \tag{ 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{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u^{\left(1\right)}+ k^2 \left(1+\eta\right) u^{\left(1\right)} +ik\sigma u^{\left(1\right)}& =& 0, &\mbox{in}\ \ \Omega\\ \Delta v^{\left(1\right)}+ \left(2k\right)^2 \left(1+\eta\right) v^{\left(1\right)} +i2k\sigma v^{\left(1\right)}& = & 0, &\mbox{in}\ \ \Omega\\ u^{\left(1\right)} = g_1, \qquad v^{\left(1\right)} = h_1, & & & \mbox{on}\ \partial\Omega \end{array} \end{align} \tag{ 9 }$$

while $H^{(1)}$ satisfies

$$\begin{align} H^{\left(1\right)} = 0\,.\end{align} \tag{ 10 }$$

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

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u^{\left(2\right)}+ k^2 \left(1+\eta\right) u^{\left(2\right)}+ik\sigma u^{(2)}& =& -4 k^2 \gamma u^{\left(1\right)*}v^{\left(1\right)}, &\mbox{in}\ \ \Omega\\ \Delta v^{\left(2\right)}+ \left(2k\right)^2 \left(1+\eta\right) v^{\left(2\right)}+i2k\sigma v^{(2)}& =& -2 \left(2k\right)^2\gamma \left(u^{\left(1\right)}\right)^2, &\mbox{in}\ \ \Omega\\ u^{\left(2\right)} = g_2, && \quad v^{(2)} = h_2, & \mbox{on}\ \partial\Omega\,. \end{array} \end{align} \tag{ 11 }$$

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

$$\begin{align} H^{\left(2\right)} = 2\Gamma \sigma \left( u^{\left(1\right)*} u^{\left(1\right)} + v^{\left(1\right)*} v^{\left(1\right)}\right)\,.\end{align} \tag{ 12 }$$

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

$$\begin{align} H^{\left(3\right)} = 3\Gamma\sigma\left(u^{\left(1\right)*}u^{\left(2\right)}+u^{\left(1\right)}u^{\left(2\right)*}+v^{\left(1\right)*}v^{\left(2\right)}+v^{\left(1\right)}v^{\left(2\right)*}\right)\,.\end{align} \tag{ 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,h_1,h_2\in {\mathcal{C}}^{2,\alpha}(\partial\Omega;\mathbb C)$. For sufficiently small ε, 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) and (6). Then the derivatives (8) are all well-defined. Moreover, (9)–(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.

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₁ = 0, the problem reduces to reconstructing $(\Gamma,\sigma,\eta)$ from the data

$$\begin{align} H^{\left(2\right)} = 2\Gamma \sigma u^{\left(1\right)*} u^{\left(1\right)}, \end{align} \tag{ 14 }$$

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

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

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u + k^2\left(1+\eta\left({\mathbf{x}}\right)\right) u + ik\sigma\left({\mathbf{x}}\right) u & =& 0, & \mbox{in}\ \Omega\\ u & = & g, & \mbox{on}\ \partial\Omega \end{array} \end{align} \tag{ 15 }$$

with internal data

$$\begin{align} H\left({\mathbf{x}}\right) = \Gamma \sigma\left({\mathbf{x}}\right)|u|^2\,, \end{align} \tag{ 16 }$$

it is shown in [11] that when Γ and η are known, one can reconstruct σ 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 σ is given (again, assuming that η and Γ are known). Here, we modify the method in order to deal with the case of unknown refractive index η and unknown Grüneisen coefficient Γ. We use the procedure to develop a uniqueness and stability result from two well-chosen datasets.

Let g₁ and g₂ 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₁ and g₂. 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:

$$\begin{align*} u_1 u_2^* = \frac12\left(|u_1+u_2|^2+i|u_1+iu_2|^2-\left(1+i\right)|u_1|^2-\left(1+i\right)|u_2|^2\right) \end{align*}$$

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

$$\begin{align*} \Gamma\sigma u_1 u_2^* = \frac12\left(\Gamma\sigma|u_1+u_2|^2+i\Gamma\sigma|u_1+iu_2|^2-\left(1+i\right)\Gamma\sigma|u_1|^2-\left(1+i\right)\Gamma\sigma|u_2|^2\right) \end{align*}$$

is known.

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

$$\begin{align} E_j = \Gamma\sigma u_j u_1^*, \end{align} \tag{ 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:

• ($\mathcal{B}-i$)  
  $E_1({\mathbf{x}})\unicode{x2A7E} \alpha_0\gt0$ for some α₀, a.e. ${\mathbf{x}}\in \Omega$.
• ($\mathcal{B}-ii$)  
  The vector field

  $$\begin{align*} \beta\left({\mathbf{x}}\right): = \nabla \frac{E_2}{E_1}\, \end{align*}$$

  is at least $W^{1,\infty}$, and $|\beta|\unicode{x2A7E} \beta_0\gt0$ for some β₀, a.e. ${\mathbf{x}}\in\Omega$.

Then Γ, η, and σ 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₁ by u₂ and multiply the equation for u₂ by u₁. We subtract the results to have

$$\begin{align*} u_1\Delta u_2-u_2\Delta u_1 = 0\,. \end{align*}$$

We can then rewrite this into

$$\begin{align*} \nabla \cdot u_1^2 \nabla \frac{u_2}{u_1} = \nabla \cdot u_1^2 \nabla \frac{E_2}{E_1} = \nabla\cdot \left(u_1^2 \beta\right) = 0\,. \end{align*}$$

The vector field β is known from the data. Therefore the above equation is a transport equation for u₁, that is,

$$\begin{align} \nabla\cdot \left(u_1^2 \beta\right) = 0,\ \ \mbox{in}\ \ \Omega, \qquad u_1 = g_1,\ \ \mbox{on}\ \ \partial\Omega\,. \end{align} \tag{ 18 }$$

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

Now that we have reconstructed u₁, we can use the equation (15) to reconstruct the potential q:

$$\begin{align} q\left({\mathbf{x}}\right): = k^2\left(1+\eta\right) + ik\sigma = -\frac{\Delta u_1}{u_1}\,. \end{align} \tag{ 19 }$$

This gives us η and σ (which are obtained by taking real and imaginary parts of q). The last step is to reconstruct Γ as

$$\begin{align} \Gamma = \frac{H_1}{\sigma|u_1|^2}\,. \end{align} \tag{ 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 shown that the multiplicative coefficient Γ 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₁, computing q using (19) and then computing Γ 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₁ and g₂ are selected such that $\frac{E_2}{E_1}$ is sufficiently small. Then we have that, for some constants c > 0,

$$\begin{align} \|\Gamma\sigma-\widetilde\Gamma\widetilde\sigma\|_{L^2\left(\Omega\right)} \unicode{x2A7D} c\left(\|H_1-\widetilde H_1\|_{L^2\left(\Omega\right)}+ \|E-\widetilde E\|_{\left(\mathcal H^2\left(\Omega\right)\right)^2} \right)\,. \end{align} \tag{ 21 }$$

Proof. We first observe that

$$\begin{align*} \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}\,. \end{align*}$$

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

$$\begin{align} \|\Gamma \sigma-\widetilde\Gamma\widetilde\sigma\|_{L^2\left(\Omega\right)}\unicode{x2A7D} \widetilde c \left(\|H_1-\widetilde H_1\|_{L^2\left(\Omega\right)} + \|\widetilde u_1^2-u_1^2\|_{L^2\left(\Omega\right)}\right)\,, \end{align} \tag{ 22 }$$

for some $\widetilde c\gt0$.

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

$$\begin{align} \nabla \cdot \left(\left(\xi-\widetilde\xi\right) \beta\right) + \nabla\cdot \left(\widetilde \xi \left(\beta-\widetilde\beta\right)\right) = 0\,. \end{align} \tag{ 23 }$$

This can be further rewritten into

$$\begin{align*} \beta\cdot\nabla\left(\xi-\widetilde\xi\right) = -\left(\xi-\widetilde\xi\right)\nabla\cdot\beta - \nabla\cdot \left(\widetilde \xi \left(\beta-\widetilde\beta\right)\right)\,, \end{align*}$$

which immediately leads to the bound

$$\begin{align} \left\|\beta\cdot\nabla\left(\xi-\widetilde\xi\right)\right\|_{L^2\left(\Omega\right)} \unicode{x2A7D} \left\|\left(\xi-\widetilde\xi\right)\nabla\cdot\beta\right\|_{L^2\left(\Omega\right)}^2+\left\|\nabla\cdot \left(\xi \left(\beta-\widetilde\beta\right)\right)\right\|_{L^2\left(\Omega\right)}\,. \end{align} \tag{ 24 }$$

With the same algebra, we can derive the bound

$$\begin{align} \left\|\widetilde \beta\cdot\nabla\left(\xi-\widetilde\xi\right)\right\|_{L^2\left(\Omega\right)}\unicode{x2A7D} \left\|\left(\xi-\widetilde\xi\right)\nabla\cdot\widetilde \beta\right\|_{L^2\left(\Omega\right)}+\left\|\nabla\cdot \left(\widetilde \xi \left(\beta-\widetilde\beta\right)\right)\right\|_{L^2\left(\Omega\right)}\,. \end{align} \tag{ 25 }$$

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

$$\begin{align*} \nabla\cdot \left(\bigg|\xi-\widetilde\xi\bigg|^2 \beta\right) - \left(\xi-\widetilde\xi\right)\beta\cdot \nabla \left(\xi-\widetilde\xi\right)^* + \left(\xi-\widetilde\xi\right)^*\nabla\cdot \widetilde\xi \left(\beta-\widetilde\beta\right) = 0\,. \end{align*}$$

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

$$\begin{align*} & \int_\Omega \bigg|\xi-\widetilde\xi\bigg|^2 \beta\cdot\nabla \phi \textrm{d}{\mathbf{x}} + \int_\Omega \left(\xi-\widetilde\xi\right)\beta \cdot \phi\nabla \left(\xi-\widetilde\xi\right)^* \textrm{d}{\mathbf{x}} \\ & \quad +\int_\Omega \left(\xi-\widetilde\xi\right)^* \widetilde\xi \left(\beta-\widetilde\beta\right) \cdot \nabla \phi \textrm{d}{\mathbf{x}} +\int_\Omega \phi \widetilde\xi \left(\beta-\widetilde\beta\right)\cdot\nabla \left(\xi-\widetilde\xi\right)^* \textrm{d}{\mathbf{x}} = 0\,. \end{align*}$$

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

$$\begin{align*} \int_\Omega \bigg|\xi-\widetilde\xi\bigg|^2 \beta\cdot\nabla \phi \textrm{d}{\mathbf{x}} + \int_\Omega \left(\beta \xi-\widetilde\beta\widetilde\xi\right)\cdot \phi\nabla \left(\xi-\widetilde\xi\right)^* \textrm{d}{\mathbf{x}} +\int_\Omega \left(\xi-\widetilde\xi\right)^* \widetilde\xi \left(\beta-\widetilde\beta\right) \cdot \nabla \phi \textrm{d}{\mathbf{x}} = 0\,. \end{align*}$$

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

$$\begin{align*} \int_\Omega \bigg|\xi-\widetilde\xi\bigg|^2 |\beta|^2 \textrm{d}{\mathbf{x}} = - \int_\Omega \left(\beta \xi-\widetilde\beta\widetilde\xi\right)\cdot \frac{E_2^*}{E_1^*}\nabla \left(\xi-\widetilde\xi\right)^* \textrm{d}{\mathbf{x}} -\int_\Omega \left(\xi-\widetilde\xi\right)^*\beta^* \cdot \widetilde\xi \left(\beta-\widetilde\beta\right) \textrm{d}{\mathbf{x}} \,. \end{align*}$$

This gives us the bound

$$\begin{align} \left\|\left(\xi-\widetilde\xi\right)\beta\right\|_{L^2\left(\Omega\right)}^2\unicode{x2A7D} \int_\Omega \bigg|\left(\beta \xi-\widetilde\beta\widetilde\xi\right)\cdot \frac{E_2^*}{E_1^*}\nabla \left(\xi-\widetilde\xi\right)^*\bigg| \textrm{d}{\mathbf{x}} + \int_\Omega \bigg|\left(\xi-\widetilde\xi\right)^*\beta^* \cdot \widetilde\xi \left(\beta-\widetilde\beta\right)\bigg| \textrm{d}{\mathbf{x}} \,. \end{align} \tag{ 26 }$$

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

$$\begin{align*} & \int_\Omega \bigg|\left(\beta \xi-\widetilde\beta\widetilde\xi\right)\cdot \frac{E_2^*}{E_1^*}\nabla \left(\xi-\widetilde\xi\right)^*\bigg| \textrm{d}{\mathbf{x}} \\ & \quad \unicode{x2A7D} \left\|\frac{E_2^*}{E_1^*}\right\|_{L^\infty\left(\Omega\right)} \int_\Omega \bigg|\left(\beta \xi-\widetilde\beta\widetilde\xi\right)\cdot \nabla \left(\xi-\widetilde\xi\right)^*\bigg| \textrm{d}{\mathbf{x}} \\ & \quad \unicode{x2A7D} \frac{1}{2}\left\|\frac{E_2^*}{E_1^*}\right\|_{L^\infty\left(\Omega\right)}\left[ \int_\Omega \bigg|\beta \xi\cdot \nabla \left(\xi-\widetilde\xi\right)^*\bigg|^2 \textrm{d}{\mathbf{x}} +\int_\Omega \bigg|\widetilde\beta \widetilde \xi\cdot \nabla \left(\xi-\widetilde\xi\right)^*\bigg|^2 \textrm{d}{\mathbf{x}}\right] \\ & \quad \unicode{x2A7D} 2\left\|\frac{E_2^*}{E_1^*}\right\|_{L^\infty\left(\Omega\right)}\left[\left\|\left(\xi-\widetilde\xi\right)\nabla\cdot\beta\right\|_{L^2\left(\Omega\right)}^2+\left\|\nabla\cdot \left(\xi \left(\beta-\widetilde\beta\right)\right)\right\|_{L^2\left(\Omega\right)}^2\right]\,, \end{align*}$$

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:

$$\begin{align*} \int_\Omega \bigg|\left(\xi-\widetilde\xi\right)^*\beta^* \cdot \widetilde\xi \left(\beta-\widetilde\beta\right)\bigg| \textrm{d}{\mathbf{x}} \unicode{x2A7D} \frac{1}{2}\left[\frac{1}{\kappa^2}\left\|\left(\xi-\widetilde\xi\right)\beta\right\|_{L^2\left(\Omega\right)}^2 +\kappa^2 \left\|\left(\beta-\widetilde\beta\right)\widetilde\xi\right\|_{L^2\left(\Omega\right)}^2\right]\,, \end{align*}$$

for any κ > 0.

Under the assumption that $|\frac{E_2}{E_1}|$ is sufficiently small, we can take κ to be sufficiently large so that (26) now implies that

$$\begin{align} \|\xi-\widetilde\xi \|_{L^2\left(\Omega\right)}^2\lesssim \left\|\left(\beta-\widetilde\beta\right)\widetilde\xi\right\|_{L^2\left(\Omega\right)}^2 +\left\|\nabla\cdot \left(\xi \left(\beta-\widetilde\beta\right)\right)\right\|_{L^2\left(\Omega\right)}^2\,.\end{align} \tag{ 27 }$$

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

$$\begin{align*} \beta-\widetilde\beta & = \left(\widetilde E_1-E_1\right)\nabla \frac{\widetilde E_2}{E_1\widetilde E_1}+\frac{\widetilde E_2}{E_1\widetilde E_1}\nabla\left(\widetilde E_1-E_1\right)\\ & \quad +\left(E_2-\widetilde E_2\right)\nabla\frac{1}{E_1}+\frac{1}{E_1^2}\nabla\left(E_2-\widetilde E_2\right)+\frac{1}{E_1\widetilde E_1}\left(\widetilde E_1-E_1\right)\nabla\left(E_2-\widetilde E_2\right)\,, \end{align*}$$

to derive the bound

$$\begin{align*} \|\beta-\widetilde \beta\|_{L^2\left(\Omega\right)}\lesssim \|E-\widetilde E\|_{\left(\mathcal H^1\left(\Omega\right)\right)^2}\,. \end{align*}$$

In a similar manner, we find that

$$\begin{align*} \left\|\nabla\cdot \xi \left(\beta-\widetilde \beta\right)\right\|_{L^2\left(\Omega\right)}\lesssim \|E-\widetilde E\|_{\left(\mathcal H^2\left(\Omega\right)\right)^2}\,. \end{align*}$$

Plugging these results into (27) will give us

$$\begin{align*} \|\xi-\widetilde\xi\|_{L^2\left(\Omega\right)} \lesssim \|E-\widetilde E\|_{\left(\mathcal H^2\left(\Omega\right)\right)^2}\,. \end{align*}$$

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\gt1+\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 Γσ.

The remaining problem is to reconstruct γ 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

$$\begin{align} \Pi_{\gamma}: \left(g, h\right) \mapsto \left(\dfrac{\partial{u}}{\partial{\nu}}\bigg|_{\partial\Omega}, \dfrac{\partial{v}}{\partial{\nu}}\bigg|_{\partial\Omega}\right)\equiv \left(J_u,J_v\right)\,. \end{align} \tag{ 28 }$$

Note that we omit the dependence of Π on Γ, η, and σ 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

$$\begin{align*} u^{\left(1\right)*}u^{\left(2\right)}+u^{\left(1\right)}u^{\left(2\right)*}+v^{\left(1\right)*}v^{\left(2\right)}+v^{\left(1\right)}v^{\left(2\right)*} = \frac{H^{\left(3\right)}}{3\Gamma\sigma}, \end{align*}$$

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

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

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u^{\left(2\right)}+ k^2 \left(1+\eta\right) u^{\left(2\right)}+ik\sigma u^{\left(2\right)} & = & -4 k^2 \gamma u^{\left(1\right)*}v^{\left(1\right)}, &\mbox{in}\ \ \Omega\\ \Delta v^{\left(2\right)}+ \left(2k\right)^2 \left(1+\eta\right) v^{\left(2\right)}+i2k\sigma v^{\left(2\right)} & = & -2 \left(2k\right)^2\gamma \left(u^{\left(1\right)}\right)^2, &\mbox{in}\ \ \Omega\\ u^{\left(2\right)} = 0,& & \quad v^{\left(2\right)} = 0,& \mbox{on}\ \ \partial\Omega \end{array} \end{align} \tag{ 29 }$$

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

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \left(\Delta + q_1\right)u^{\left(2\right)} + 4k^2 u^{\left(1\right)*}v^{\left(1\right)}\gamma & = & 0,&\mbox{in}\ \Omega\\ \left(\Delta + q_1^*\right)u^{\left(2\right)*} + 4k^2 u^{\left(1\right)}v^{\left(1\right)*}\gamma & = & 0,&\mbox{in}\ \Omega\\ \left(\Delta + q_2\right)v^{\left(2\right)} + 2\left(2k\right)^2\left(u^{\left(1\right)}\right)^2\gamma & = & 0, &\mbox{in}\ \Omega\\ \left(\Delta + q_2^*\right)v^{\left(2\right)*} + 2\left(2k\right)^2\left(u^{\left(1\right)*}\right)^2\gamma & = & 0, &\mbox{in}\ \Omega\\ u^{\left(1\right)*}u^{\left(2\right)}+u^{\left(1\right)}u^{\left(2\right)*}+v^{\left(1\right)*}v^{\left(2\right)}+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} \end{align} \tag{ 30 }$$

where we have used the notation

$$\begin{align*} q_1 : = k^2\left(1+\eta\right) + ik\sigma,\qquad q_2 : = \left(2k\right)^2\left(1+\eta\right) + i2k\sigma\,.\end{align*}$$

If we can solve (30), we can reconstruct γ (and the associated $(u^{(2)}, v^{(2)})$). This is a non-iterative reconstruction method. In the rest of this section, we show that γ 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 Γ, η, σ 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)\gt0$ such that for every $\Omega\subset B({\mathbf{x}}_0, \epsilon)$, there exist g₁ and h₁ in (6) such that

$$\begin{align*} \left(H^{\left(3\right)}, J_u^{\left(2\right)}, J_v^{\left(2\right)}\right) = \left(\widetilde H^{\left(3\right)}, \widetilde J_u^{\left(2\right)}, \widetilde J_v^{\left(2\right)}\right) \implies \gamma = \widetilde \gamma\,, \end{align*}$$

and for all p > 1,

$$\begin{align} \|\gamma-\widetilde{\gamma}\|_{L^p\left(\Omega\right)} & \unicode{x2A7D} C\left(\bigg\|H^{\left(3\right)} - \widetilde{H}^{\left(3\right)}\bigg\|_{W^{{\,}2,p}\left(\Omega\right)} + \bigg\|J_u^{\left(2\right)} - \widetilde J_u^{\left(2\right)}\bigg\|_{W^{1-1/p, p}\left(\partial\Omega\right)} + \bigg\|J_v^{\left(2\right)} - \widetilde J_v^{\left(2\right)}\bigg\|_{W^{1-1/p, p}\left(\partial\Omega\right)} \right),\nonumber \end{align} \tag{ 31 }$$

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

Proof. Let us define

$$\begin{align*} \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. \end{align*}$$

The linear system (30) then implies that

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \left(\Delta + q_1\right)\widehat{u}^{\left(2\right)} + 4k^2 u^{\left(1\right)*}v^{\left(1\right)}\widehat{\gamma} & = & 0,&\mbox{in}\ \ \Omega\\ \left(\Delta + q_1^*\right)\widehat{u}^{\left(2\right)*} + 4k^2 u^{\left(1\right)}v^{\left(1\right)*}\widehat{\gamma} & = & 0,&\mbox{in}\ \ \Omega\\ \left(\Delta + q_2\right)\widehat{v}^{\left(2\right)} + 2\left(2k\right)^2\left(u^{\left(1\right)}\right)^2\widehat{\gamma} & = & 0, &\mbox{in}\ \ \Omega\\ \left(\Delta + q_2^*\right)\widehat{v}^{\left(2\right)*} + 2\left(2k\right)^2\left(u^{\left(1\right)*}\right)^2\widehat{\gamma} & = & 0, &\mbox{in}\ \Omega\\ u^{\left(1\right)*}\widehat{u}^{\left(2\right)}+u^{\left(1\right)}\widehat{u}^{\left(2\right)*}+v^{\left(1\right)*}\widehat{v}^{\left(2\right)}+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} \end{align} \tag{ 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

$$\begin{align*} \left(\Delta + q_1\right)\widehat{u}^{\left(2\right)} -\frac{v^{\left(1\right)}}{2u^{\left(1\right)*}}\left(\Delta + q_2^*\right)\widehat{v}^{\left(2\right)*} & = 0,\\ \left(\Delta + q_1^*\right)\widehat{u}^{\left(2\right)*} - \frac{u^{\left(1\right)}v^{\left(1\right)*}}{2\left(u^{\left(1\right)*}\right)^2}\left(\Delta + q_2^*\right)\widehat{v}^{\left(2\right)*} & = 0,\\ \left(\Delta + q_2\right)\widehat{v}^{\left(2\right)} - \frac{\left(u^{\left(1\right)}\right)^2}{\left(u^{\left(1\right)*}\right)^2}\left(\Delta + q_2^*\right)\widehat{v}^{\left(2\right)*} & = 0, \\ \Delta\left(u^{\left(1\right)*}\widehat{u}^{(2)}\right)+\Delta(u^{(1)}\widehat{u}^{(2)*})+\Delta(v^{(1)*}\widehat{v}^{(2)})+\Delta(v^{(1)}\widehat{v}^{(2)*}) & = \Delta\left(\frac{\widehat{H}^{(3)}}{3\Gamma\sigma}\right), \end{align*}$$

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{align} \begin{array}{r@{\ }c@{\ }ll} \mathcal A w & = & \mathcal S, & \mbox{in}\ \ \Omega\\ w& = & \left(0,0,0,0\right),&\mbox{on}\ \partial\Omega\\ \dfrac{\partial{w}}{\partial{\nu}} & = & \left(\widehat J_u^{\left(2\right)},\widehat J_u^{\left(2\right)*},\widehat J_v^{\left(2\right)},\widehat J_v^{\left(2\right)*}\right),&\mbox{on}\ \partial\Omega\,, \end{array} \end{align} \tag{ 33 }$$

where

$$\begin{align*} \mathcal A\left({\mathbf{x}},D\right) : = \begin{pmatrix} \Delta+q_1 & & & -\frac{v^{\left(1\right)}}{2u^{\left(1\right)*}}\left(\Delta + q_2^*\right)\\[4pt] & \Delta+q_1^*& & -\frac{u^{\left(1\right)}v^{\left(1\right)*}}{2\left(u^{\left(1\right)*}\right)^2}\left(\Delta + q_2^*\right)\\[4pt] & & \Delta+q_2 & -\frac{\left(u^{\left(1\right)}\right)^2}{\left(u^{\left(1\right)*}\right)^2}\left(\Delta + q_2^*\right)\\[4pt] \Delta\left(u^{\left(1\right)*}\cdot\right) & \Delta\left(u^{\left(1\right)}\cdot\right) & \Delta\left(v^{\left(1\right)*}\cdot\right) & \Delta\left(v^{\left(1\right)}\cdot\right) \end{pmatrix}\,, \end{align*}$$

and

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

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

$$\begin{align} \left(s_1,s_2,s_3,s_4\right) = \left(0,0,0,0\right),\qquad \left(t_1,t_2,t_3,t_4\right) = \left(2,2,2,2\right). \end{align} \tag{ 34 }$$

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

$$\begin{align*} \mathcal A_0\left({\mathbf{x}},{\boldsymbol{\xi}}\right) : = \begin{pmatrix} |{\boldsymbol{\xi}}|^2 & & & -\frac{v^{\left(1\right)}}{2u^{\left(1\right)*}}|{\boldsymbol{\xi}}|^2\\ & |{\boldsymbol{\xi}}|^2& & -\frac{u^{\left(1\right)}v^{\left(1\right)*}}{2\left(u^{\left(1\right)*}\right)^2}|{\boldsymbol{\xi}}|^2\\ & & |{\boldsymbol{\xi}}|^2 & -\frac{\left(u^{\left(1\right)}\right)^2}{\left(u^{\left(1\right)*}\right)^2}|{\boldsymbol{\xi}}|^2\\ u^{\left(1\right)*}|{\boldsymbol{\xi}}|^2 & u^{\left(1\right)}|{\boldsymbol{\xi}}|^2 & v^{\left(1\right)*}|{\boldsymbol{\xi}}|^2 & v^{\left(1\right)}|{\boldsymbol{\xi}}|^2 \end{pmatrix}\,.\end{align*}$$

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$:

$$\begin{align*} -\frac{v^{\left(1\right)}}{2u^{\left(1\right)*}}\cdot u^{\left(1\right)*} -\frac{u^{\left(1\right)}v^{\left(1\right)*}}{2\left(u^{\left(1\right)*}\right)^2}\cdot u^{\left(1\right)} -\frac{\left(u^{\left(1\right)}\right)^2}{\left(u^{\left(1\right)*}\right)^2}\cdot v^{\left(1\right)*} \neq v^{\left(1\right)}, \end{align*}$$

or equivalently

$$\begin{align*} -\frac{\left(u^{\left(1\right)}\right)^2}{\left(u^{\left(1\right)*}\right)^2} \neq \frac{v^{\left(1\right)}}{v^{\left(1\right)*}} \end{align*}$$

at ${\mathbf{x}}_0$. This condition on $u^{(1)}({\mathbf{x}}_0)$ and $v^{(1)}({\mathbf{x}}_0)$ is easily achieved by selecting g₁ and h₁ appropriately. To be precise, let us consider some ball $B({\mathbf{x}}_0,\epsilon_0) \subset X$ and let u₀ and v₀ be any ${\mathcal{C}}^2$ functions on $B({\mathbf{x}}_0,\epsilon_0)$ satisfying

$$\begin{align*} \Delta u_0+ q_1 u_0 & = 0,\quad \mbox{in}\ B\left({\mathbf{x}}_0,\epsilon_0\right), \\ \Delta v_0 + q_2 v_0 & = 0,\quad \mbox{in}\ B\left({\mathbf{x}}_0,\epsilon_0\right), \\ -\frac{u_0^2}{\left(u_0^*\right)^2}\left({\mathbf{x}}_0\right) &\neq \frac{v_0}{v_0^*}\left({\mathbf{x}}_0\right). \end{align*}$$

The existence of such u₀ and v₀ is obvious as we can take u₀ and v₀ to be any solutions to the first two equations which are nonzero at ${\mathbf{x}}_0$ 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₀ and v₀ depend only on η and σ, not γ.

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

$$\begin{align*} u^{\left(1\right)} = u_0|_\Omega,\qquad v^{\left(1\right)} = v_0|_\Omega\,. \end{align*}$$

This means we can write $\mathcal A$ as

$$\begin{align*} \mathcal A\left({\mathbf{x}},D\right) : = \begin{pmatrix} \Delta+q_1 & & & -\frac{v_0}{2u_0^*}\left(\Delta + q_2^*\right)\\[4pt] & \Delta+q_1^*& & -\frac{u_0 v_0^*}{2\left(u_0^*\right)^2}\left(\Delta + q_2^*\right)\\[4pt] & & \Delta+q_2 & -\frac{u_0^2}{\left(u_0^*\right)^2}\left(\Delta + q_2^*\right)\\[4pt] \Delta\left(u_0^*\cdot\right) & \Delta\left(u_0\cdot\right) & \Delta\left(v_0^*\cdot\right) & \Delta\left(v_0\cdot\right) \end{pmatrix}\,. \end{align*}$$

By construction, the constant-coefficient operator $\mathcal A({\mathbf{x}}_0,D)$ with coefficients frozen at ${\mathbf{x}}_0$ is elliptic. Additionally, from the continuity of u₀ and v₀ we see that there exists $\epsilon_1 = \epsilon_1(\eta,\sigma,{\mathbf{x}}_0)\gt0$ 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 Ω.

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)\gt0$ such that for every $\Omega\subset B({\mathbf{x}}_0,\epsilon_2)$, the boundary value problem (33) has a unique solution.

Now set $\epsilon = \text{min}\{\epsilon_1,\epsilon_2\}\gt0$ and let $\Omega\subset B({\mathbf{x}}_0,\epsilon)$, so that $\mathcal A$ is an elliptic operator on Ω 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

$$\begin{align*} \|\widehat{u}^{\left(2\right)}\|_{W^{{\,}2,p}\left(\Omega\right)} + \|\widehat{v}^{\left(2\right)}\|_{W^{{\,}2,p}\left(\Omega\right)}&\unicode{x2A7D} C\left(\left\|\Delta\left(\frac{\widehat{H}^{\left(3\right)}}{3\Gamma\sigma}\right)\right\|_{L^p\left(\Omega\right)} + \left\|\widehat J_u^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)}\right.\\ & \left.\quad + \left\|\widehat J_v^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)}\right) + C_2\left(\|\widehat{u}^{\left(2\right)}\|_{L^p\left(\Omega\right)} + \|\widehat{v}^{\left(2\right)}\|_{L^p\left(\Omega\right)}\right). \end{align*}$$

In fact, since the solution is unique we may drop the last term on the right-hand side, i.e. set C₂ = 0. This gives

$$\begin{align*} &\|\widehat{u}^{\left(2\right)}\|_{W^{{\,}2,p}\left(\Omega\right)} + \|\widehat{v}^{\left(2\right)}\|_{W^{{\,}2,p}\left(\Omega\right)}\\ & \quad \unicode{x2A7D} C\left(\left\|\Delta\left(\frac{\widehat{H}^{\left(3\right)}}{3\Gamma\sigma}\right)\right\|_{L^p\left(\Omega\right)} + \left\|\widehat J_u^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)} + \left\|\widehat J_v^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)}\right)\\ & \quad \unicode{x2A7D} C\left(\|\widehat{H}^{\left(3\right)}\|_{W^{{\,}2,p}\left(\Omega\right)} + \left\|\widehat J_u^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)} + \left\|\widehat J_v^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)}\right), \end{align*}$$

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

Finally, we substitute this estimate into (32) to obtain

$$\begin{align*} \|\widehat{\gamma}\|_{L^p\left(\Omega\right)} & = \left\|-\frac{\left(\Delta + q_1\right)\widehat{u}^{\left(2\right)}}{4k^2 u_0^* v_0}\right\|_{L^p\left(\Omega\right)} \unicode{x2A7D} C\|\widehat{u}^{\left(2\right)}\|_{W^{{\,}2,p}\left(\Omega\right)}\\ &\unicode{x2A7D} C\left(\|\widehat{H}^{\left(3\right)}\|_{W^{{\,}2,p}\left(\Omega\right)} + \left\|\widehat J_u^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)} + \left\|\widehat J_v^{\left(2\right)}\right\|_{W^{1-1/p,p}\left(\partial\Omega\right)}\right), \end{align*}$$

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 γ requires both the availability of the additional boundary data (28) and the assumption that Ω 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.

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{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u+ k^2 \left(1+\eta\right) u +ik\sigma u& = & 0, & \mbox{in}\ \ \Omega\\ \Delta v+ \left(2k\right)^2 \left(1+\eta\right) v +i2k\sigma v & = & -\left(2k\right)^2 \gamma u^2, &\mbox{in}\ \ \Omega\\ u = g, \qquad v+i2k{\boldsymbol \nu}\cdot \nabla v & = & 0, & \mbox{on}\ \partial\Omega \end{array} \end{align} \tag{ 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₁ = 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\unicode{x2A7E} 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).

5.1. Numerical experiment I: reconstructing γ

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

$$\begin{align*} \Phi\left(\gamma\right) : = \dfrac{1}{2}\sum_{j = 1}^{N_s}\left\|\Gamma\sigma\left(|u_j|^2 + |v_j|^2\right) -H_j\right\|_{L^2\left(\Omega\right)}^2+\frac12\beta \left\|\nabla \gamma\right\|_{L^2\left(\Omega\right)}^2\,, \end{align*}$$

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

$$\begin{align*} \begin{array}{rll} \Delta w_j+ \left(2k\right)^2 \left(1+\eta\right) w_j +i2k\sigma w_j & = -\left[\Gamma\sigma\left(|u_j|^2+|v_j|^2\right) - H_j\right]\Gamma\sigma v_j^*, &\mbox{in}\ \ \Omega\\ w_j+i2k{\boldsymbol \nu}\cdot \nabla w_j & = 0, & \mbox{on}\ \partial\Omega \end{array} \end{align*}$$

It is then straightforward to verify that the Fréchet derivative of Φ in direction δγ can be written as:

$$\begin{align*} \Phi^{^{\prime}}\left(\gamma\right)\left[\delta\gamma\right] = \int_{\Omega} \left(2\left(2k\right)^2 \Re\sum_{j = 1}^{N_s} w_j u_j^2\right)\delta\gamma\,\textrm{d}{\mathbf{x}} -\beta\int_{\Omega} \left(\Delta\gamma\right)\delta\gamma\,\textrm{d}{\mathbf{x}} + \beta\int_{\partial\Omega} \dfrac{\partial{\gamma}}{\partial{\nu}}\delta\gamma\,dS \end{align*}$$

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

Figure 1 shows the reconstruction of a simple profile of γ 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 γ profiles we tested. The regularization parameter is selected in a trial-and-error manner. The value of β 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.

Zoom In Zoom Out Reset image size

5.2. Numerical experiment II: reconstructing $(\eta, \sigma,\gamma)$

In the second numerical example, we consider the case where Γ is known but η, σ, and γ 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 η and γ are simple.

Zoom In Zoom Out Reset image size

5.3. Numerical experiment III: reconstructing $(\eta,\gamma,\Gamma)$

In this example, we assume that σ is known and we are interested in reconstructing η, γ and Γ. Due to the fact that Γ 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 Γ 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 Γ. That is, we minimize the functional

$$\begin{align*} \Psi\left(\eta,\gamma\right) : = \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\left(\Omega\right)}^2 + \frac12\beta_1\|\nabla \eta\|_{L^2\left(\Omega\right)}^2 + \frac12\beta_2\|\nabla \gamma\|_{L^2\left(\Omega\right)}^2\,, \end{align*}$$

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

$$\begin{align*} \begin{array}{r@{\ }c@{\ }ll} \Delta w_j+ \left(2k\right)^2 \left(1+\eta\right) 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} \end{align*}$$

and

$$\begin{align*} \begin{array}{r@{\ }c@{\ }ll} \Delta z_j+ \left(2k\right)^2 \left(1+\eta\right) 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}{\left(|u_1|^2+|v_1|^2\right)^2}v_1^*, &\mbox{in}\ \ \Omega\\ z_j+i2k{\boldsymbol \nu}\cdot \nabla z_j & = & 0, & \mbox{on}\ \partial\Omega. \end{array} \end{align*}$$

It is then straightforward to verify that the Fréchet derivative of Ψ with respect to γ in direction δγ can be written as:

$$\begin{align*} \Psi^{^{\prime}}_{\gamma}\left(\eta,\gamma\right)\left[\delta\gamma\right] = \int_{\Omega} \left(2\left(2k\right)^2\Re\sum_{j = 1}^{N_s} \left[w_j u_j^2 + z_j u_1^2\right]\right)\delta\gamma\,\textrm{d}{\mathbf{x}} -\beta_2\int_{\Omega} \left(\Delta\gamma\right)\delta\gamma\,\textrm{d}{\mathbf{x}} + \beta_2\int_{\partial\Omega} \dfrac{\partial{\gamma}}{\partial{\nu}}\delta\gamma\,dS\,. \end{align*}$$

The Fréchet derivative with respect to η can be computed in a similar fashion. Once η and γ are reconstructed, we can reconstruct Γ as

$$\begin{align} \Gamma = \frac{1}{N_s}\sum_{j = 1}^{N_s} \frac{H_j}{\sigma \left(|u_j|^2+|v_j|^2\right)}\,. \end{align} \tag{ 36 }$$

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

Zoom In Zoom Out Reset image size

5.4. 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² error. Note from the reconstruction formula (36) that any inaccuracies in the reconstruction of σ will directly translate into artifacts in the reconstruction of Γ. This can be observed in figure 4(see columns 2 and 4), most notably near the edges of the square anomaly in σ.

Zoom In Zoom Out Reset image size

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 Ω 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 γ.

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 ε-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.

We would like to thank the anonymous referees for their useful comments that helped us improve the quality of this work. This work is partially supported by the National Science Foundation through Grants DMS-1913309 and DMS-1937254.

No new data were created or analysed in this study.

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 \gt0$, then the boundary value problem

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u + qu & = & f, & \mbox{in}\ \Omega\\ u & = & g, & \mbox{on}\ \partial\Omega \end{array} \end{align} \tag{ 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:

$$\begin{align} \|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\left(\|{\,}f{\,}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right)\,. \end{align} \tag{ 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{align} \begin{array}{r@{\ }c@{\ }ll} \Delta u + qu & = & 0, & \mbox{in}\ \Omega\\ u & = & 0, & \mbox{on}\ \partial\Omega\,. \end{array} \end{align} \tag{ 39 }$$

Multiplying both sides of the PDE by $u^*$ and integrating over Ω results in

$$\begin{align*} \int_\Omega -|\nabla u|^2 + q|u|^2\, \textrm{d}{\mathbf{x}} = 0\, , \end{align*}$$

whereupon taking imaginary parts yields $\int_\Omega \Im q |u|^2\, \textrm{d}{\mathbf{x}} = 0$. The assumption $\Im q \gt 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

$$\begin{align*} \|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\left(\|{\,}f{\,}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}+\|u\|_{{\mathcal{C}}^0\left(\overline{\Omega}\right)}\right)\,. \end{align*}$$

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 × 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 × Y. Finally, for δ > 0 define the complete metric space

$$\begin{align*} \mathcal X_\delta : = \left\{f\in {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right): \|{\,}f{\,}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\unicode{x2A7D} \delta\right\}\,. \end{align*}$$

To start, let ɛ > 0, δ > 0 and let $g,h\in {\mathcal{C}}^{2,\alpha}(\partial\Omega;\mathbb C)$ with $\|g\|_{{\mathcal{C}}^{2,\alpha}(\partial\Omega)}\lt\varepsilon$ and $\|h\|_{{\mathcal{C}}^{2,\alpha}(\partial\Omega)}\lt\varepsilon$. We shall determine how small ɛ and δ 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{align*} \begin{array}{r@{\ }c@{\ }ll} \Delta U+ q_1U & = & f_1, &\mbox{in}\ \ \Omega\\ \Delta V+ q_2V & = & f_2, &\mbox{in}\ \ \Omega\\ U = g,& &\quad V = h, &\mbox{on}\ \ \partial\Omega\,. \end{array} \end{align*}$$

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

$$\begin{align*} T\left(\phi_1,\phi_2\right) : = L^{-1}\left(-2k^2\gamma \phi_1^* \phi_2,\ -\left(2k\right)^2\gamma \phi_1^2\right)\,. \end{align*}$$

It remains to show that for sufficiently small ε > 0 and δ > 0,

• (i)  
  T is a well-defined operator from $\mathcal X_\delta\times \mathcal X_\delta$ to itself, and
• (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

$$\begin{align*} \|fg\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\unicode{x2A7D} \|{\,}f{\,}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|g\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} \end{align*}$$

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

$$\begin{align*} \|-2k^2\gamma \phi_1^* \phi_2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\|\phi_1\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|\phi_2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\|\phi_1\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\|\phi_2\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\delta^2 \end{align*}$$

and similarly

$$\begin{align*} \|-\left(2k\right)^2\gamma \phi_1^2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\|\phi_1\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} ^2 \unicode{x2A7D} C\|\phi_1\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} ^2 \unicode{x2A7D} C\delta^2 \end{align*}$$

Let $(U,V) : = T(\phi_1,\phi_2)$. Combining the above estimates with the Schauder estimate (38) for the Helmholtz equation then leads to

$$\begin{align*} \|U\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\left(\|-2k^2\gamma \phi_1^* \phi_2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \unicode{x2A7D} C\left(\delta^2 + \varepsilon\right) \end{align*}$$

and

$$\begin{align*} \|V\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\left(\|-\left(2k\right)^2\gamma \phi_1^2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \unicode{x2A7D} C\left(\delta^2 + \varepsilon\right). \end{align*}$$

We can force these quantities to be less than δ by choosing ɛ and δ 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

$$\begin{align*} \|-2k^2\gamma \phi_1^*\phi_2 - \left(-2k^2\gamma \left(\phi_1^{^{\prime}}\right)^*\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} &\unicode{x2A7D} C\|\phi_1^*\phi_2 - \left(\phi_1^{^{\prime}}\right)^*\phi_2^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\\ & = C\|\phi_1^*\left(\phi_2-\phi_2^{^{\prime}}\right) + \left(\phi_1^*-\left(\phi_1^{^{\prime}}\right)^*\right)\phi_2^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\\ &\unicode{x2A7D} C\left(\|\phi_1\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|\phi_2-\phi_2^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\right.\nonumber\\ & \quad \left.+ \|\phi_1-\phi_1^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|\phi_2^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\right)\\ &\unicode{x2A7D} C\left(\|\phi_1\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\|\phi_2-\phi_2^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \right.\\ & \quad \left.+ \|\phi_1-\phi_1^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\|\phi_2^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\right)\\ &\unicode{x2A7D} C\delta \|\left(\phi_1,\phi_2\right) - \left(\phi_1^{^{\prime}},\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)\times {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \end{align*}$$

and similarly

$$\begin{align*} \left\|-\left(2k\right)^2\gamma\phi_1^2 - \left(-\left(2k\right)^2\gamma\left(\phi_1^{^{\prime}}\right)^2\right)\right\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} &\unicode{x2A7D} C\|\phi_1^2-\left(\phi_1^{^{\prime}}\right)^2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\\ &\unicode{x2A7D} C\|\phi_1-\phi_1^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|\phi_1+\phi_1^{^{\prime}}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\\ &\unicode{x2A7D} C\delta\|\phi_1-\phi_1^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\\ &\unicode{x2A7D} C\delta \|\left(\phi_1,\phi_2\right) - \left(\phi_1^{^{\prime}},\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)\times {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}. \end{align*}$$

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{align*} \begin{array}{r@{\ }c@{\ }ll} \Delta \left(U-U^{^{\prime}}\right)+ q_1\left(U-U^{^{\prime}}\right) & = & -2k^2\gamma \phi_1^*\phi_2 - \left(-2k^2\gamma \left(\phi_1^{^{\prime}}\right)^*\phi_2^{^{\prime}}\right), &\mbox{in}\ \ \Omega\\ \Delta \left(V-V^{^{\prime}}\right)+ q_2\left(V-V^{^{\prime}}\right) & = & -\left(2k\right)^2\gamma\phi_1^2 - \left(-\left(2k\right)^2\gamma\left(\phi_1^{^{\prime}}\right)^2\right), &\mbox{in}\ \ \Omega\\ U-U^{^{\prime}} = 0, \qquad V-V^{^{\prime}} & = & 0, &\mbox{on}\ \ \partial\Omega\,. \end{array} \end{align*}$$

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

$$\begin{align*} \|U-U^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} & \unicode{x2A7D} C\|-2k^2\gamma \phi_1^*\phi_2 - \left(-2k^2\gamma \left(\phi_1^{^{\prime}}\right)^*\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\\ & \unicode{x2A7D} C\delta \|\left(\phi_1,\phi_2\right) - \left(\phi_1^{^{\prime}},\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)\times {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \end{align*}$$

and

$$\begin{align*} \|V-V^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} & \unicode{x2A7D} C\left\|-\left(2k\right)^2\gamma\phi_1^2 - \left(-\left(2k\right)^2\gamma\left(\phi_1^{^{\prime}}\right)^2\right)\right\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\\ & \unicode{x2A7D} C\delta \|\left(\phi_1,\phi_2\right) - \left(\phi_1^{^{\prime}},\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)\times {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}. \end{align*}$$

We conclude that

$$\begin{align*} \|T\left(\phi_1,\phi_2\right) - T\left(\phi_1^{^{\prime}},\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)\times {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} & = \|U-U^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} + \|V-V^{^{\prime}}\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\\ &\unicode{x2A7D} C\delta \|\left(\phi_1,\phi_2\right) - \left(\phi_1^{^{\prime}},\phi_2^{^{\prime}}\right)\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)\times {\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}. \end{align*}$$

The factor $C\delta$ can be made strictly less than 1 when δ 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

$$\begin{align*} \|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} &\unicode{x2A7D} C\left(\|-2k^2\gamma u^* v\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \\ & \unicode{x2A7D} C\left(\|u\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|v\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right)\\ &\unicode{x2A7D} C\left(\|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\delta + \|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \end{align*}$$

When δ is sufficiently small, this implies that

$$\begin{align*} \|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\,, \end{align*}$$

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

$$\begin{align*} \|v\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} &\unicode{x2A7D} C\left(\|-\left(2k\right)^2 \gamma u^2\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \unicode{x2A7D} C\left(\|u\|_{{\mathcal{C}}^{0,\alpha}\left(\partial\Omega\right)}^2 + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right)\\ &\unicode{x2A7D} C\left(\delta \|u\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \unicode{x2A7D} C\left(\delta\|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right)\\ &\unicode{x2A7D} C\left(\|g\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)} + \|h\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right)\,, \end{align*}$$

as desired. The proof is complete. □

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_ɛ are differentiable with respect to ɛ.

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{align*} \begin{array}{r@{\ }c@{\ }ll} \Delta \widetilde{u}^{\left(1\right)}+ q_1\widetilde{u}^{\left(1\right)} & = & 0, &\mbox{in}\ \ \Omega\\ \Delta \widetilde{v}^{\left(1\right)}+ q_2\widetilde{v}^{\left(1\right)} & = & 0, &\mbox{in}\ \ \Omega\\ \widetilde{u}^{\left(1\right)} = g_1,& &\quad \widetilde{v}^{\left(1\right)} = h_1, &\mbox{on}\ \ \partial\Omega \end{array} \end{align*}$$

and

$$\begin{align*} \begin{array}{r@{\ }c@{\ }ll} \Delta \widetilde{u}^{\left(2\right)}+ q_1\widetilde{u}^{\left(2\right)} & = & -4k^2\gamma \widetilde{u}^{\left(1\right)*}\widetilde{v}^{\left(1\right)}, &\mbox{in}\ \ \Omega\\ \Delta \widetilde{v}^{\left(2\right)}+ q_2\widetilde{v}^{\left(2\right)} & = & -2\left(2k\right)^2\gamma \left(\widetilde{u}^{\left(1\right)}\right)^2, &\mbox{in}\ \ \Omega\\ \widetilde{u}^{\left(2\right)} = g_2,& &\quad \widetilde{v}^{\left(2\right)} = h_2,&\mbox{on}\ \ \partial\Omega\,. \end{array} \end{align*}$$

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

$$\begin{align} \begin{aligned} \mu_\varepsilon &: = u_\varepsilon - \varepsilon\widetilde{u}^{\left(1\right)} - \frac12\varepsilon^2\widetilde{u}^{\left(2\right)},\\ \nu_\varepsilon &: = v_\varepsilon - \varepsilon\widetilde{v}^{\left(1\right)} - \frac12\varepsilon^2\widetilde{v}^{\left(2\right)}. \end{aligned} \end{align} \tag{ 40 }$$

We wish to show that µ_ɛ and ν_ɛ are in a certain sense '$o(\varepsilon^2)$' as $\varepsilon\rightarrow 0$. This will be accomplished in two rounds of estimates on µ_ɛ , ν_ɛ , u_ɛ and v_ɛ .

Round 1 estimates. We begin by using the linearity of the operators $\Delta + q_1$ and $\Delta + q_2$ to find that

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \Delta \mu_\varepsilon + q_1\mu_\varepsilon & = & -2k^2\gamma \left[u_\varepsilon^*v_\varepsilon - \varepsilon^2 \widetilde{u}^{\left(1\right)*}\widetilde{v}^{\left(1\right)}\right], &\mbox{in}\ \ \Omega\\ \Delta \nu_\varepsilon + q_2\nu_\varepsilon & = & -\left(2k\right)^2\gamma \left[u_\varepsilon^2 - \varepsilon^2\left(\widetilde{u}^{\left(1\right)}\right)^2\right], &\mbox{in}\ \ \Omega\\ \mu_\varepsilon = 0,& &\quad \nu_\varepsilon = 0,&\mbox{on}\ \ \partial\Omega\,. \end{array} \end{align} \tag{ 41 }$$

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

$$\begin{align*} \|u_\varepsilon\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\unicode{x2A7D} \|u_\varepsilon\|_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)} \unicode{x2A7D} C\left(\left\|\varepsilon g_1 + \frac12\varepsilon^2 g_2\right\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)} + \left\|\varepsilon h_1 + \frac12\varepsilon^2 h_2\right\|_{{\mathcal{C}}^{2,\alpha}\left(\partial\Omega\right)}\right) \unicode{x2A7D} C \varepsilon, \end{align*}$$

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

$$\begin{align} u_\varepsilon = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon\right),\qquad v_\varepsilon = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon\right). \end{align} \tag{ 42 }$$

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

$$\begin{align*} \|fg\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\unicode{x2A7D} \|{\,}f{\,}\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\|g\|_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)} \end{align*}$$

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

$$\begin{align} -2k^2\gamma \left[u_\varepsilon^*v_\varepsilon - \varepsilon^2 \widetilde{u}^{\left(1\right)*}\widetilde{v}^{\left(1\right)}\right] = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right),\qquad -\left(2k\right)^2\gamma \left[u_\varepsilon^2 - \varepsilon^2\left(\widetilde{u}^{\left(1\right)}\right)^2\right] = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right). \end{align} \tag{ 43 }$$

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

$$\begin{align} \mu_\varepsilon = \mathcal O_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right),\qquad \nu_\varepsilon = \mathcal O_{{\mathcal{C}}^{2,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right), \end{align} \tag{ 44 }$$

and in particular

$$\begin{align} \mu_\varepsilon = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right),\qquad \nu_\varepsilon = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right). \end{align} \tag{ 45 }$$

Round 2 estimates. Using the estimates from Round 1 and recalling the definition (40) of the remainder terms µ_ɛ and ν_ɛ , we can now refine the bounds (43) on the right hand sides of (41):

$$\begin{align*} & -2k^2\gamma \left[u_\varepsilon^*v_\varepsilon - \varepsilon^2 \widetilde{u}^{\left(1\right)*}\widetilde{v}^{\left(1\right)}\right]\\ & \quad = -2k^2\gamma\left[\left(\varepsilon\widetilde{u}^{\left(1\right)*} + \frac12\varepsilon^2\widetilde{u}^{\left(2\right)*}+\mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right)\right)\right.\\ & \qquad \times \left.\left(\varepsilon \widetilde{v}^{\left(1\right)} + \frac12\varepsilon^2\widetilde{v}^{\left(2\right)}+\mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right)\right) - \varepsilon^2 \widetilde{u}^{\left(1\right)*}\widetilde{v}^{\left(1\right)}\right]\\ & \quad = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^3\right) \end{align*}$$

and

$$\begin{align*} -\left(2k\right)^2\gamma \left[u_\varepsilon^2 - \varepsilon^2\left(\widetilde{u}^{\left(1\right)}\right)^2\right] & = -\left(2k\right)^2\gamma \left[\left(\varepsilon\widetilde{u}^{\left(1\right)} + \frac12\varepsilon^2\widetilde{u}^{\left(2\right)} + \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^2\right)\right)^2 - \varepsilon^2\left(\widetilde{u}^{\left(1\right)}\right)^2\right] \\ & = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^3\right). \end{align*}$$

With these improved bounds in hand, we again apply the Schauder estimate (38)–(41) to obtain the following bounds for the remainder terms µ_ɛ and ν_ɛ :

$$\begin{align*} \mu_\varepsilon = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^3\right),\qquad \nu_\varepsilon = \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^3\right). \end{align*}$$

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

Therefore from the definition (40) of µ_ɛ and ν_ɛ we conclude that

$$\begin{align*} u_\varepsilon & = \varepsilon\widetilde{u}^{\left(1\right)} + \frac12\varepsilon^2\widetilde{u}^{\left(2\right)} + \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^3\right),\\ v_\varepsilon & = \varepsilon\widetilde{v}^{\left(1\right)} + \frac12\varepsilon^2\widetilde{v}^{\left(2\right)} + \mathcal O_{{\mathcal{C}}^{0,\alpha}\left(\overline{\Omega}\right)}\left(\varepsilon^3\right), \end{align*}$$

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

$$\begin{align*} u_\varepsilon\left({\mathbf{x}}\right) & = \varepsilon\widetilde{u}^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12\varepsilon^2\widetilde{u}^{\left(2\right)}\left({\mathbf{x}}\right) + \mathcal O\left(\varepsilon^3\right),\\ v_\varepsilon\left({\mathbf{x}}\right) & = \varepsilon\widetilde{v}^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12\varepsilon^2\widetilde{v}^{\left(2\right)}\left({\mathbf{x}}\right) + \mathcal O\left(\varepsilon^3\right), \end{align*}$$

as $\varepsilon\rightarrow 0$.

Asymptotic expansion of H_ɛ . To finish the proof, for each ${\mathbf{x}}\in\overline{\Omega}$ we compute:

$$\begin{align*} H_\varepsilon\left({\mathbf{x}}\right) & = \Gamma\left({\mathbf{x}}\right)\sigma\left({\mathbf{x}}\right)\left[u_\varepsilon\left({\mathbf{x}}\right)u^*_\varepsilon\left({\mathbf{x}}\right) + v_\varepsilon\left({\mathbf{x}}\right)v^*_\varepsilon\left({\mathbf{x}}\right)\right]\\ & = \Gamma\left({\mathbf{x}}\right)\sigma\left({\mathbf{x}}\right)\left[\left(\varepsilon\widetilde{u}^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12\varepsilon^2\widetilde{u}^{\left(2\right)}\left({\mathbf{x}}\right) + \mathcal O\left(\varepsilon^3\right)\right)\left(\varepsilon\widetilde{u}^{\left(1\right)*}\left({\mathbf{x}}\right) + \frac12\varepsilon^2\widetilde{u}^{\left(2\right)*}\left({\mathbf{x}}\right) + \mathcal O\left(\varepsilon^3\right)\right)\right.\\ & \quad \left. + \left(\varepsilon\widetilde{v}^{\left(1\right)}\left({\mathbf{x}}\right) + \frac12\varepsilon^2\widetilde{v}^{\left.(2\right)}({\mathbf{x}}) + \mathcal O(\varepsilon^3)\right)\left(\varepsilon\widetilde{v}^{(1)*}({\mathbf{x}}) + \frac12\varepsilon^2\widetilde{v}^{(2)*}({\mathbf{x}}) + \mathcal O(\varepsilon^3)\right) \right]\\ & = \varepsilon H^{(1)}({\mathbf{x}}) + \frac12\varepsilon^2 H^{(2)}({\mathbf{x}}) + \frac16\varepsilon^3 H^{(3)}({\mathbf{x}}) + \mathcal O(\varepsilon^4), \end{align*}$$

where

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

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

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\unicode{x2A7E} 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 Ω:

$$\begin{align*} \mathcal A\left({\mathbf{x}},D\right) v = \mathcal S\quad\mbox{in}\ \Omega. \end{align*}$$

Here $v = (v_1,\dots, v_N)$, $D = (\partial_{x_1},\dots,\partial_{x_n})$, and $\mathcal A({\mathbf{x}},D)$ is an M × N matrix linear partial differential operator. That is, for each $1\unicode{x2A7D} i\unicode{x2A7D} M$ and $1\unicode{x2A7D} j\unicode{x2A7D} N$, the entry $\mathcal A_{ij}({\mathbf{x}},D)$ is a linear partial differential operator, and the above matrix equation means that

$$\begin{align*} \sum_{j = 1}^N \mathcal A_{ij}\left({\mathbf{x}},D\right)v_j = \mathcal S_i,\qquad 1\unicode{x2A7D} i\unicode{x2A7D} M. \end{align*}$$

To each row $1\unicode{x2A7D} i\unicode{x2A7D} M$ let us now associate an integer s_i , and to each column $1\unicode{x2A7D} j\unicode{x2A7D} 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\lt0$, 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 × 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\unicode{x2A7D} i\unicode{x2A7D} M}$ and $(t_j)_{1\unicode{x2A7D} j\unicode{x2A7D} N}$ exist, and the matrix $\mathcal A_0({\mathbf{x}},{\mathbf{x}}i)$ has full rank N for each ${\mathbf{x}}\in\Omega$ and ${\mathbf{x}}i\in \mathbb S^{n-1}$. (The matrix $\mathcal A_0({\mathbf{x}},{\mathbf{x}}i)$ 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\unicode{x2A7D} i\unicode{x2A7D} M}$ and $(t_j)_{1\unicode{x2A7D} j\unicode{x2A7D} N}$ such that

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} s_i & = & 0, &1\unicode{x2A7D} i\unicode{x2A7D} M,\\ t_j & = & \tau, &1\unicode{x2A7D} j\unicode{x2A7D} N, \end{array} \end{align} \tag{ 46 }$$

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

$$\begin{align} \begin{array}{r@{\ }c@{\ }ll} \mathcal A\left({\mathbf{x}},D\right)v & = & \mathcal S &\mbox{in}\ \Omega,\\ \left(\dfrac{\partial{}}{\partial{\nu}}\right)^q v_j & = & \phi_{qj} &\mbox{on}\ \partial\Omega,\quad 0\unicode{x2A7D} q\unicode{x2A7D} \tau-1, \quad 1\unicode{x2A7D} j\unicode{x2A7D} N. \end{array} \end{align} \tag{ 47 }$$

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

Theorem C.2 ([8] p 13). Let p > 1 and $\ell \unicode{x2A7E} 0$. Assume that $\mathcal S_i\in W^{\ell,p}(\Omega)$ for all $1\unicode{x2A7D} i\unicode{x2A7D} M$ and $\phi_{qj}\in W^{\ell+\tau-q-\frac{1}{p},p}(\partial\Omega)$ for all $0\unicode{x2A7D} q\unicode{x2A7D} \tau-1$ and $1\unicode{x2A7D} j\unicode{x2A7D} 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):

$$\begin{align} \sum_{j = 1}^{N} \|v_j\|_{W^{\ell+\tau,p}\left(\Omega\right)} &\unicode{x2A7D} C\left(\sum_{i = 1}^M \|\mathcal S_i\|_{W^{\ell,p}\left(\Omega\right)} + \sum_{q,j}\|\phi_{qj}\|_{W^{\ell+\tau-q-\frac{1}{p},p}\left(\partial\Omega\right)}\right) + C_2\sum_{j = 1}^N \|v_j\|_{L^p\left(\Omega\right)}. \end{align} \tag{ 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 ε > 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{align} \begin{array}{r@{\ }c@{\ }ll} \mathcal A\left({\mathbf{x}},D\right)v & = & 0 &\mbox{in}\ \Omega^{^{\prime}},\\ \left(\dfrac{\partial{}}{\partial{\nu}}\right)^q v_j & = & 0 &\mbox{on}\ \partial\Omega^{^{\prime}},\quad 0\unicode{x2A7D} q\unicode{x2A7D} \tau-1, \quad 1\unicode{x2A7D} j\unicode{x2A7D} N. \end{array} \end{align} \tag{ 49 }$$

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